1252884aeSStefan Eßer# Algorithms 2252884aeSStefan Eßer 3252884aeSStefan EßerThis `bc` uses the math algorithms below: 4252884aeSStefan Eßer 5252884aeSStefan Eßer### Addition 6252884aeSStefan Eßer 7252884aeSStefan EßerThis `bc` uses brute force addition, which is linear (`O(n)`) in the number of 8252884aeSStefan Eßerdigits. 9252884aeSStefan Eßer 10252884aeSStefan Eßer### Subtraction 11252884aeSStefan Eßer 12252884aeSStefan EßerThis `bc` uses brute force subtraction, which is linear (`O(n)`) in the number 13252884aeSStefan Eßerof digits. 14252884aeSStefan Eßer 15252884aeSStefan Eßer### Multiplication 16252884aeSStefan Eßer 17252884aeSStefan EßerThis `bc` uses two algorithms: [Karatsuba][1] and brute force. 18252884aeSStefan Eßer 19252884aeSStefan EßerKaratsuba is used for "large" numbers. ("Large" numbers are defined as any 20252884aeSStefan Eßernumber with `BC_NUM_KARATSUBA_LEN` digits or larger. `BC_NUM_KARATSUBA_LEN` has 21252884aeSStefan Eßera sane default, but may be configured by the user.) Karatsuba, as implemented in 22252884aeSStefan Eßerthis `bc`, is superlinear but subpolynomial (bounded by `O(n^log_2(3))`). 23252884aeSStefan Eßer 24252884aeSStefan EßerBrute force multiplication is used below `BC_NUM_KARATSUBA_LEN` digits. It is 25252884aeSStefan Eßerpolynomial (`O(n^2)`), but since Karatsuba requires both more intermediate 26252884aeSStefan Eßervalues (which translate to memory allocations) and a few more additions, there 27252884aeSStefan Eßeris a "break even" point in the number of digits where brute force multiplication 28252884aeSStefan Eßeris faster than Karatsuba. There is a script (`$ROOT/karatsuba.py`) that will 29252884aeSStefan Eßerfind the break even point on a particular machine. 30252884aeSStefan Eßer 31252884aeSStefan Eßer***WARNING: The Karatsuba script requires Python 3.*** 32252884aeSStefan Eßer 33252884aeSStefan Eßer### Division 34252884aeSStefan Eßer 35252884aeSStefan EßerThis `bc` uses Algorithm D ([long division][2]). Long division is polynomial 36252884aeSStefan Eßer(`O(n^2)`), but unlike Karatsuba, any division "divide and conquer" algorithm 37252884aeSStefan Eßerreaches its "break even" point with significantly larger numbers. "Fast" 38252884aeSStefan Eßeralgorithms become less attractive with division as this operation typically 39252884aeSStefan Eßerreduces the problem size. 40252884aeSStefan Eßer 41252884aeSStefan EßerWhile the implementation of long division may appear to use the subtractive 42252884aeSStefan Eßerchunking method, it only uses subtraction to find a quotient digit. It avoids 43252884aeSStefan Eßerunnecessary work by aligning digits prior to performing subtraction and finding 44252884aeSStefan Eßera starting guess for the quotient. 45252884aeSStefan Eßer 46252884aeSStefan EßerSubtraction was used instead of multiplication for two reasons: 47252884aeSStefan Eßer 48252884aeSStefan Eßer1. Division and subtraction can share code (one of the less important goals of 49252884aeSStefan Eßer this `bc` is small code). 50252884aeSStefan Eßer2. It minimizes algorithmic complexity. 51252884aeSStefan Eßer 52252884aeSStefan EßerUsing multiplication would make division have the even worse algorithmic 53252884aeSStefan Eßercomplexity of `O(n^(2*log_2(3)))` (best case) and `O(n^3)` (worst case). 54252884aeSStefan Eßer 55252884aeSStefan Eßer### Power 56252884aeSStefan Eßer 57252884aeSStefan EßerThis `bc` implements [Exponentiation by Squaring][3], which (via Karatsuba) has 58252884aeSStefan Eßera complexity of `O((n*log(n))^log_2(3))` which is favorable to the 59252884aeSStefan Eßer`O((n*log(n))^2)` without Karatsuba. 60252884aeSStefan Eßer 61252884aeSStefan Eßer### Square Root 62252884aeSStefan Eßer 63252884aeSStefan EßerThis `bc` implements the fast algorithm [Newton's Method][4] (also known as the 64252884aeSStefan EßerNewton-Raphson Method, or the [Babylonian Method][5]) to perform the square root 65252884aeSStefan Eßeroperation. Its complexity is `O(log(n)*n^2)` as it requires one division per 66252884aeSStefan Eßeriteration. 67252884aeSStefan Eßer 68252884aeSStefan Eßer### Sine and Cosine (`bc` Only) 69252884aeSStefan Eßer 70252884aeSStefan EßerThis `bc` uses the series 71252884aeSStefan Eßer 72252884aeSStefan Eßer``` 73252884aeSStefan Eßerx - x^3/3! + x^5/5! - x^7/7! + ... 74252884aeSStefan Eßer``` 75252884aeSStefan Eßer 76252884aeSStefan Eßerto calculate `sin(x)` and `cos(x)`. It also uses the relation 77252884aeSStefan Eßer 78252884aeSStefan Eßer``` 79252884aeSStefan Eßercos(x) = sin(x + pi/2) 80252884aeSStefan Eßer``` 81252884aeSStefan Eßer 82252884aeSStefan Eßerto calculate `cos(x)`. It has a complexity of `O(n^3)`. 83252884aeSStefan Eßer 84252884aeSStefan Eßer**Note**: this series has a tendency to *occasionally* produce an error of 1 85252884aeSStefan Eßer[ULP][6]. (It is an unfortunate side effect of the algorithm, and there isn't 86252884aeSStefan Eßerany way around it; [this article][7] explains why calculating sine and cosine, 87252884aeSStefan Eßerand the other transcendental functions below, within less than 1 ULP is nearly 88252884aeSStefan Eßerimpossible and unnecessary.) Therefore, I recommend that users do their 89252884aeSStefan Eßercalculations with the precision (`scale`) set to at least 1 greater than is 90252884aeSStefan Eßerneeded. 91252884aeSStefan Eßer 92252884aeSStefan Eßer### Exponentiation (`bc` Only) 93252884aeSStefan Eßer 94252884aeSStefan EßerThis `bc` uses the series 95252884aeSStefan Eßer 96252884aeSStefan Eßer``` 97252884aeSStefan Eßer1 + x + x^2/2! + x^3/3! + ... 98252884aeSStefan Eßer``` 99252884aeSStefan Eßer 100252884aeSStefan Eßerto calculate `e^x`. Since this only works when `x` is small, it uses 101252884aeSStefan Eßer 102252884aeSStefan Eßer``` 103252884aeSStefan Eßere^x = (e^(x/2))^2 104252884aeSStefan Eßer``` 105252884aeSStefan Eßer 106252884aeSStefan Eßerto reduce `x`. It has a complexity of `O(n^3)`. 107252884aeSStefan Eßer 108252884aeSStefan Eßer**Note**: this series can also produce errors of 1 ULP, so I recommend users do 109252884aeSStefan Eßertheir calculations with the precision (`scale`) set to at least 1 greater than 110252884aeSStefan Eßeris needed. 111252884aeSStefan Eßer 112252884aeSStefan Eßer### Natural Logarithm (`bc` Only) 113252884aeSStefan Eßer 114252884aeSStefan EßerThis `bc` uses the series 115252884aeSStefan Eßer 116252884aeSStefan Eßer``` 117252884aeSStefan Eßera + a^3/3 + a^5/5 + ... 118252884aeSStefan Eßer``` 119252884aeSStefan Eßer 120252884aeSStefan Eßer(where `a` is equal to `(x - 1)/(x + 1)`) to calculate `ln(x)` when `x` is small 121252884aeSStefan Eßerand uses the relation 122252884aeSStefan Eßer 123252884aeSStefan Eßer``` 124252884aeSStefan Eßerln(x^2) = 2 * ln(x) 125252884aeSStefan Eßer``` 126252884aeSStefan Eßer 127252884aeSStefan Eßerto sufficiently reduce `x`. It has a complexity of `O(n^3)`. 128252884aeSStefan Eßer 129252884aeSStefan Eßer**Note**: this series can also produce errors of 1 ULP, so I recommend users do 130252884aeSStefan Eßertheir calculations with the precision (`scale`) set to at least 1 greater than 131252884aeSStefan Eßeris needed. 132252884aeSStefan Eßer 133252884aeSStefan Eßer### Arctangent (`bc` Only) 134252884aeSStefan Eßer 135252884aeSStefan EßerThis `bc` uses the series 136252884aeSStefan Eßer 137252884aeSStefan Eßer``` 138252884aeSStefan Eßerx - x^3/3 + x^5/5 - x^7/7 + ... 139252884aeSStefan Eßer``` 140252884aeSStefan Eßer 141252884aeSStefan Eßerto calculate `atan(x)` for small `x` and the relation 142252884aeSStefan Eßer 143252884aeSStefan Eßer``` 144252884aeSStefan Eßeratan(x) = atan(c) + atan((x - c)/(1 + x * c)) 145252884aeSStefan Eßer``` 146252884aeSStefan Eßer 147252884aeSStefan Eßerto reduce `x` to small enough. It has a complexity of `O(n^3)`. 148252884aeSStefan Eßer 149252884aeSStefan Eßer**Note**: this series can also produce errors of 1 ULP, so I recommend users do 150252884aeSStefan Eßertheir calculations with the precision (`scale`) set to at least 1 greater than 151252884aeSStefan Eßeris needed. 152252884aeSStefan Eßer 153252884aeSStefan Eßer### Bessel (`bc` Only) 154252884aeSStefan Eßer 155252884aeSStefan EßerThis `bc` uses the series 156252884aeSStefan Eßer 157252884aeSStefan Eßer``` 158252884aeSStefan Eßerx^n/(2^n * n!) * (1 - x^2 * 2 * 1! * (n + 1)) + x^4/(2^4 * 2! * (n + 1) * (n + 2)) - ... 159252884aeSStefan Eßer``` 160252884aeSStefan Eßer 161252884aeSStefan Eßerto calculate the bessel function (integer order only). 162252884aeSStefan Eßer 163252884aeSStefan EßerIt also uses the relation 164252884aeSStefan Eßer 165252884aeSStefan Eßer``` 166252884aeSStefan Eßerj(-n,x) = (-1)^n * j(n,x) 167252884aeSStefan Eßer``` 168252884aeSStefan Eßer 169252884aeSStefan Eßerto calculate the bessel when `x < 0`, It has a complexity of `O(n^3)`. 170252884aeSStefan Eßer 171252884aeSStefan Eßer**Note**: this series can also produce errors of 1 ULP, so I recommend users do 172252884aeSStefan Eßertheir calculations with the precision (`scale`) set to at least 1 greater than 173252884aeSStefan Eßeris needed. 174252884aeSStefan Eßer 175252884aeSStefan Eßer### Modular Exponentiation (`dc` Only) 176252884aeSStefan Eßer 177252884aeSStefan EßerThis `dc` uses the [Memory-efficient method][8] to compute modular 178252884aeSStefan Eßerexponentiation. The complexity is `O(e*n^2)`, which may initially seem 179252884aeSStefan Eßerinefficient, but `n` is kept small by maintaining small numbers. In practice, it 180252884aeSStefan Eßeris extremely fast. 181252884aeSStefan Eßer 182252884aeSStefan Eßer[1]: https://en.wikipedia.org/wiki/Karatsuba_algorithm 183252884aeSStefan Eßer[2]: https://en.wikipedia.org/wiki/Long_division 184252884aeSStefan Eßer[3]: https://en.wikipedia.org/wiki/Exponentiation_by_squaring 185252884aeSStefan Eßer[4]: https://en.wikipedia.org/wiki/Newton%27s_method#Square_root_of_a_number 186252884aeSStefan Eßer[5]: https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method 187252884aeSStefan Eßer[6]: https://en.wikipedia.org/wiki/Unit_in_the_last_place 188252884aeSStefan Eßer[7]: https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT 189252884aeSStefan Eßer[8]: https://en.wikipedia.org/wiki/Modular_exponentiation#Memory-efficient_method 190