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 2844d4804dSStefan Eßeris faster than Karatsuba. There is a script (`$ROOT/scripts/karatsuba.py`) that 2944d4804dSStefan Eßerwill find 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 6544d4804dSStefan Eßeroperation. 66252884aeSStefan Eßer 6744d4804dSStefan EßerIts complexity is `O(log(n)*n^2)` as it requires one division per iteration, and 6844d4804dSStefan Eßerit doubles the amount of correct digits per iteration. 6944d4804dSStefan Eßer 7044d4804dSStefan Eßer### Sine and Cosine (`bc` Math Library Only) 71252884aeSStefan Eßer 72252884aeSStefan EßerThis `bc` uses the series 73252884aeSStefan Eßer 74252884aeSStefan Eßer``` 75252884aeSStefan Eßerx - x^3/3! + x^5/5! - x^7/7! + ... 76252884aeSStefan Eßer``` 77252884aeSStefan Eßer 78252884aeSStefan Eßerto calculate `sin(x)` and `cos(x)`. It also uses the relation 79252884aeSStefan Eßer 80252884aeSStefan Eßer``` 81252884aeSStefan Eßercos(x) = sin(x + pi/2) 82252884aeSStefan Eßer``` 83252884aeSStefan Eßer 84252884aeSStefan Eßerto calculate `cos(x)`. It has a complexity of `O(n^3)`. 85252884aeSStefan Eßer 86252884aeSStefan Eßer**Note**: this series has a tendency to *occasionally* produce an error of 1 87252884aeSStefan Eßer[ULP][6]. (It is an unfortunate side effect of the algorithm, and there isn't 88252884aeSStefan Eßerany way around it; [this article][7] explains why calculating sine and cosine, 89252884aeSStefan Eßerand the other transcendental functions below, within less than 1 ULP is nearly 90252884aeSStefan Eßerimpossible and unnecessary.) Therefore, I recommend that users do their 91252884aeSStefan Eßercalculations with the precision (`scale`) set to at least 1 greater than is 92252884aeSStefan Eßerneeded. 93252884aeSStefan Eßer 9444d4804dSStefan Eßer### Exponentiation (`bc` Math Library Only) 95252884aeSStefan Eßer 96252884aeSStefan EßerThis `bc` uses the series 97252884aeSStefan Eßer 98252884aeSStefan Eßer``` 99252884aeSStefan Eßer1 + x + x^2/2! + x^3/3! + ... 100252884aeSStefan Eßer``` 101252884aeSStefan Eßer 102252884aeSStefan Eßerto calculate `e^x`. Since this only works when `x` is small, it uses 103252884aeSStefan Eßer 104252884aeSStefan Eßer``` 105252884aeSStefan Eßere^x = (e^(x/2))^2 106252884aeSStefan Eßer``` 107252884aeSStefan Eßer 10844d4804dSStefan Eßerto reduce `x`. 10944d4804dSStefan Eßer 11044d4804dSStefan EßerIt has a complexity of `O(n^3)`. 111252884aeSStefan Eßer 112252884aeSStefan Eßer**Note**: this series can also produce errors of 1 ULP, so I recommend users do 113252884aeSStefan Eßertheir calculations with the precision (`scale`) set to at least 1 greater than 114252884aeSStefan Eßeris needed. 115252884aeSStefan Eßer 11644d4804dSStefan Eßer### Natural Logarithm (`bc` Math Library Only) 117252884aeSStefan Eßer 118252884aeSStefan EßerThis `bc` uses the series 119252884aeSStefan Eßer 120252884aeSStefan Eßer``` 121252884aeSStefan Eßera + a^3/3 + a^5/5 + ... 122252884aeSStefan Eßer``` 123252884aeSStefan Eßer 124252884aeSStefan Eßer(where `a` is equal to `(x - 1)/(x + 1)`) to calculate `ln(x)` when `x` is small 125252884aeSStefan Eßerand uses the relation 126252884aeSStefan Eßer 127252884aeSStefan Eßer``` 128252884aeSStefan Eßerln(x^2) = 2 * ln(x) 129252884aeSStefan Eßer``` 130252884aeSStefan Eßer 13144d4804dSStefan Eßerto sufficiently reduce `x`. 13244d4804dSStefan Eßer 13344d4804dSStefan EßerIt has a complexity of `O(n^3)`. 134252884aeSStefan Eßer 135252884aeSStefan Eßer**Note**: this series can also produce errors of 1 ULP, so I recommend users do 136252884aeSStefan Eßertheir calculations with the precision (`scale`) set to at least 1 greater than 137252884aeSStefan Eßeris needed. 138252884aeSStefan Eßer 13944d4804dSStefan Eßer### Arctangent (`bc` Math Library Only) 140252884aeSStefan Eßer 141252884aeSStefan EßerThis `bc` uses the series 142252884aeSStefan Eßer 143252884aeSStefan Eßer``` 144252884aeSStefan Eßerx - x^3/3 + x^5/5 - x^7/7 + ... 145252884aeSStefan Eßer``` 146252884aeSStefan Eßer 147252884aeSStefan Eßerto calculate `atan(x)` for small `x` and the relation 148252884aeSStefan Eßer 149252884aeSStefan Eßer``` 150252884aeSStefan Eßeratan(x) = atan(c) + atan((x - c)/(1 + x * c)) 151252884aeSStefan Eßer``` 152252884aeSStefan Eßer 153252884aeSStefan Eßerto reduce `x` to small enough. It has a complexity of `O(n^3)`. 154252884aeSStefan Eßer 155252884aeSStefan Eßer**Note**: this series can also produce errors of 1 ULP, so I recommend users do 156252884aeSStefan Eßertheir calculations with the precision (`scale`) set to at least 1 greater than 157252884aeSStefan Eßeris needed. 158252884aeSStefan Eßer 15944d4804dSStefan Eßer### Bessel (`bc` Math Library Only) 160252884aeSStefan Eßer 161252884aeSStefan EßerThis `bc` uses the series 162252884aeSStefan Eßer 163252884aeSStefan Eßer``` 164252884aeSStefan Eßerx^n/(2^n * n!) * (1 - x^2 * 2 * 1! * (n + 1)) + x^4/(2^4 * 2! * (n + 1) * (n + 2)) - ... 165252884aeSStefan Eßer``` 166252884aeSStefan Eßer 167252884aeSStefan Eßerto calculate the bessel function (integer order only). 168252884aeSStefan Eßer 169252884aeSStefan EßerIt also uses the relation 170252884aeSStefan Eßer 171252884aeSStefan Eßer``` 172252884aeSStefan Eßerj(-n,x) = (-1)^n * j(n,x) 173252884aeSStefan Eßer``` 174252884aeSStefan Eßer 175252884aeSStefan Eßerto calculate the bessel when `x < 0`, It has a complexity of `O(n^3)`. 176252884aeSStefan Eßer 177252884aeSStefan Eßer**Note**: this series can also produce errors of 1 ULP, so I recommend users do 178252884aeSStefan Eßertheir calculations with the precision (`scale`) set to at least 1 greater than 179252884aeSStefan Eßeris needed. 180252884aeSStefan Eßer 181252884aeSStefan Eßer### Modular Exponentiation (`dc` Only) 182252884aeSStefan Eßer 183252884aeSStefan EßerThis `dc` uses the [Memory-efficient method][8] to compute modular 184252884aeSStefan Eßerexponentiation. The complexity is `O(e*n^2)`, which may initially seem 185252884aeSStefan Eßerinefficient, but `n` is kept small by maintaining small numbers. In practice, it 186252884aeSStefan Eßeris extremely fast. 187252884aeSStefan Eßer 18844d4804dSStefan Eßer### Non-Integer Exponentiation (`bc` Math Library 2 Only) 18944d4804dSStefan Eßer 19044d4804dSStefan EßerThis is implemented in the function `p(x,y)`. 19144d4804dSStefan Eßer 19244d4804dSStefan EßerThe algorithm used is to use the formula `e(y*l(x))`. 19344d4804dSStefan Eßer 19444d4804dSStefan EßerIt has a complexity of `O(n^3)` because both `e()` and `l()` do. 19544d4804dSStefan Eßer 19644d4804dSStefan Eßer### Rounding (`bc` Math Library 2 Only) 19744d4804dSStefan Eßer 19844d4804dSStefan EßerThis is implemented in the function `r(x,p)`. 19944d4804dSStefan Eßer 20044d4804dSStefan EßerThe algorithm is a simple method to check if rounding away from zero is 20144d4804dSStefan Eßernecessary, and if so, adds `1e10^p`. 20244d4804dSStefan Eßer 20344d4804dSStefan EßerIt has a complexity of `O(n)` because of add. 20444d4804dSStefan Eßer 20544d4804dSStefan Eßer### Ceiling (`bc` Math Library 2 Only) 20644d4804dSStefan Eßer 20744d4804dSStefan EßerThis is implemented in the function `ceil(x,p)`. 20844d4804dSStefan Eßer 20944d4804dSStefan EßerThe algorithm is a simple add of one less decimal place than `p`. 21044d4804dSStefan Eßer 21144d4804dSStefan EßerIt has a complexity of `O(n)` because of add. 21244d4804dSStefan Eßer 21344d4804dSStefan Eßer### Factorial (`bc` Math Library 2 Only) 21444d4804dSStefan Eßer 21544d4804dSStefan EßerThis is implemented in the function `f(n)`. 21644d4804dSStefan Eßer 21744d4804dSStefan EßerThe algorithm is a simple multiplication loop. 21844d4804dSStefan Eßer 21944d4804dSStefan EßerIt has a complexity of `O(n^3)` because of linear amount of `O(n^2)` 22044d4804dSStefan Eßermultiplications. 22144d4804dSStefan Eßer 22244d4804dSStefan Eßer### Permutations (`bc` Math Library 2 Only) 22344d4804dSStefan Eßer 22444d4804dSStefan EßerThis is implemented in the function `perm(n,k)`. 22544d4804dSStefan Eßer 22644d4804dSStefan EßerThe algorithm is to use the formula `n!/(n-k)!`. 22744d4804dSStefan Eßer 22844d4804dSStefan EßerIt has a complexity of `O(n^3)` because of the division and factorials. 22944d4804dSStefan Eßer 23044d4804dSStefan Eßer### Combinations (`bc` Math Library 2 Only) 23144d4804dSStefan Eßer 23244d4804dSStefan EßerThis is implemented in the function `comb(n,r)`. 23344d4804dSStefan Eßer 23444d4804dSStefan EßerThe algorithm is to use the formula `n!/r!*(n-r)!`. 23544d4804dSStefan Eßer 23644d4804dSStefan EßerIt has a complexity of `O(n^3)` because of the division and factorials. 23744d4804dSStefan Eßer 23844d4804dSStefan Eßer### Logarithm of Any Base (`bc` Math Library 2 Only) 23944d4804dSStefan Eßer 24044d4804dSStefan EßerThis is implemented in the function `log(x,b)`. 24144d4804dSStefan Eßer 24244d4804dSStefan EßerThe algorithm is to use the formula `l(x)/l(b)` with double the `scale` because 24344d4804dSStefan Eßerthere is no good way of knowing how many digits of precision are needed when 24444d4804dSStefan Eßerswitching bases. 24544d4804dSStefan Eßer 24644d4804dSStefan EßerIt has a complexity of `O(n^3)` because of the division and `l()`. 24744d4804dSStefan Eßer 24844d4804dSStefan Eßer### Logarithm of Base 2 (`bc` Math Library 2 Only) 24944d4804dSStefan Eßer 25044d4804dSStefan EßerThis is implemented in the function `l2(x)`. 25144d4804dSStefan Eßer 25244d4804dSStefan EßerThis is a convenience wrapper around `log(x,2)`. 25344d4804dSStefan Eßer 25444d4804dSStefan Eßer### Logarithm of Base 10 (`bc` Math Library 2 Only) 25544d4804dSStefan Eßer 25644d4804dSStefan EßerThis is implemented in the function `l10(x)`. 25744d4804dSStefan Eßer 25844d4804dSStefan EßerThis is a convenience wrapper around `log(x,10)`. 25944d4804dSStefan Eßer 26044d4804dSStefan Eßer### Root (`bc` Math Library 2 Only) 26144d4804dSStefan Eßer 26244d4804dSStefan EßerThis is implemented in the function `root(x,n)`. 26344d4804dSStefan Eßer 26444d4804dSStefan EßerThe algorithm is [Newton's method][9]. The initial guess is calculated as 26544d4804dSStefan Eßer`10^ceil(length(x)/n)`. 26644d4804dSStefan Eßer 26744d4804dSStefan EßerLike square root, its complexity is `O(log(n)*n^2)` as it requires one division 26844d4804dSStefan Eßerper iteration, and it doubles the amount of correct digits per iteration. 26944d4804dSStefan Eßer 27044d4804dSStefan Eßer### Cube Root (`bc` Math Library 2 Only) 27144d4804dSStefan Eßer 27244d4804dSStefan EßerThis is implemented in the function `cbrt(x)`. 27344d4804dSStefan Eßer 27444d4804dSStefan EßerThis is a convenience wrapper around `root(x,3)`. 27544d4804dSStefan Eßer 27644d4804dSStefan Eßer### Greatest Common Divisor (`bc` Math Library 2 Only) 27744d4804dSStefan Eßer 27844d4804dSStefan EßerThis is implemented in the function `gcd(a,b)`. 27944d4804dSStefan Eßer 28044d4804dSStefan EßerThe algorithm is an iterative version of the [Euclidean Algorithm][10]. 28144d4804dSStefan Eßer 28244d4804dSStefan EßerIt has a complexity of `O(n^4)` because it has a linear number of divisions. 28344d4804dSStefan Eßer 28444d4804dSStefan EßerThis function ensures that `a` is always bigger than `b` before starting the 28544d4804dSStefan Eßeralgorithm. 28644d4804dSStefan Eßer 28744d4804dSStefan Eßer### Least Common Multiple (`bc` Math Library 2 Only) 28844d4804dSStefan Eßer 28944d4804dSStefan EßerThis is implemented in the function `lcm(a,b)`. 29044d4804dSStefan Eßer 29144d4804dSStefan EßerThe algorithm uses the formula `a*b/gcd(a,b)`. 29244d4804dSStefan Eßer 29344d4804dSStefan EßerIt has a complexity of `O(n^4)` because of `gcd()`. 29444d4804dSStefan Eßer 29544d4804dSStefan Eßer### Pi (`bc` Math Library 2 Only) 29644d4804dSStefan Eßer 29744d4804dSStefan EßerThis is implemented in the function `pi(s)`. 29844d4804dSStefan Eßer 29944d4804dSStefan EßerThe algorithm uses the formula `4*a(1)`. 30044d4804dSStefan Eßer 30144d4804dSStefan EßerIt has a complexity of `O(n^3)` because of arctangent. 30244d4804dSStefan Eßer 30344d4804dSStefan Eßer### Tangent (`bc` Math Library 2 Only) 30444d4804dSStefan Eßer 30544d4804dSStefan EßerThis is implemented in the function `t(x)`. 30644d4804dSStefan Eßer 30744d4804dSStefan EßerThe algorithm uses the formula `s(x)/c(x)`. 30844d4804dSStefan Eßer 30944d4804dSStefan EßerIt has a complexity of `O(n^3)` because of sine, cosine, and division. 31044d4804dSStefan Eßer 31144d4804dSStefan Eßer### Atan2 (`bc` Math Library 2 Only) 31244d4804dSStefan Eßer 31344d4804dSStefan EßerThis is implemented in the function `a2(y,x)`. 31444d4804dSStefan Eßer 31544d4804dSStefan EßerThe algorithm uses the [standard formulas][11]. 31644d4804dSStefan Eßer 31744d4804dSStefan EßerIt has a complexity of `O(n^3)` because of arctangent. 31844d4804dSStefan Eßer 319252884aeSStefan Eßer[1]: https://en.wikipedia.org/wiki/Karatsuba_algorithm 320252884aeSStefan Eßer[2]: https://en.wikipedia.org/wiki/Long_division 321252884aeSStefan Eßer[3]: https://en.wikipedia.org/wiki/Exponentiation_by_squaring 322252884aeSStefan Eßer[4]: https://en.wikipedia.org/wiki/Newton%27s_method#Square_root_of_a_number 323252884aeSStefan Eßer[5]: https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method 324252884aeSStefan Eßer[6]: https://en.wikipedia.org/wiki/Unit_in_the_last_place 325252884aeSStefan Eßer[7]: https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT 326252884aeSStefan Eßer[8]: https://en.wikipedia.org/wiki/Modular_exponentiation#Memory-efficient_method 32744d4804dSStefan Eßer[9]: https://en.wikipedia.org/wiki/Root-finding_algorithms#Newton's_method_(and_similar_derivative-based_methods) 32844d4804dSStefan Eßer[10]: https://en.wikipedia.org/wiki/Euclidean_algorithm 32944d4804dSStefan Eßer[11]: https://en.wikipedia.org/wiki/Atan2#Definition_and_computation 330