1# Algorithms 2 3This `bc` uses the math algorithms below: 4 5### Addition 6 7This `bc` uses brute force addition, which is linear (`O(n)`) in the number of 8digits. 9 10### Subtraction 11 12This `bc` uses brute force subtraction, which is linear (`O(n)`) in the number 13of digits. 14 15### Multiplication 16 17This `bc` uses two algorithms: [Karatsuba][1] and brute force. 18 19Karatsuba is used for "large" numbers. ("Large" numbers are defined as any 20number with `BC_NUM_KARATSUBA_LEN` digits or larger. `BC_NUM_KARATSUBA_LEN` has 21a sane default, but may be configured by the user.) Karatsuba, as implemented in 22this `bc`, is superlinear but subpolynomial (bounded by `O(n^log_2(3))`). 23 24Brute force multiplication is used below `BC_NUM_KARATSUBA_LEN` digits. It is 25polynomial (`O(n^2)`), but since Karatsuba requires both more intermediate 26values (which translate to memory allocations) and a few more additions, there 27is a "break even" point in the number of digits where brute force multiplication 28is faster than Karatsuba. There is a script (`$ROOT/scripts/karatsuba.py`) that 29will find the break even point on a particular machine. 30 31***WARNING: The Karatsuba script requires Python 3.*** 32 33### Division 34 35This `bc` uses Algorithm D ([long division][2]). Long division is polynomial 36(`O(n^2)`), but unlike Karatsuba, any division "divide and conquer" algorithm 37reaches its "break even" point with significantly larger numbers. "Fast" 38algorithms become less attractive with division as this operation typically 39reduces the problem size. 40 41While the implementation of long division may appear to use the subtractive 42chunking method, it only uses subtraction to find a quotient digit. It avoids 43unnecessary work by aligning digits prior to performing subtraction and finding 44a starting guess for the quotient. 45 46Subtraction was used instead of multiplication for two reasons: 47 481. Division and subtraction can share code (one of the less important goals of 49 this `bc` is small code). 502. It minimizes algorithmic complexity. 51 52Using multiplication would make division have the even worse algorithmic 53complexity of `O(n^(2*log_2(3)))` (best case) and `O(n^3)` (worst case). 54 55### Power 56 57This `bc` implements [Exponentiation by Squaring][3], which (via Karatsuba) has 58a complexity of `O((n*log(n))^log_2(3))` which is favorable to the 59`O((n*log(n))^2)` without Karatsuba. 60 61### Square Root 62 63This `bc` implements the fast algorithm [Newton's Method][4] (also known as the 64Newton-Raphson Method, or the [Babylonian Method][5]) to perform the square root 65operation. 66 67Its complexity is `O(log(n)*n^2)` as it requires one division per iteration, and 68it doubles the amount of correct digits per iteration. 69 70### Sine and Cosine (`bc` Math Library Only) 71 72This `bc` uses the series 73 74``` 75x - x^3/3! + x^5/5! - x^7/7! + ... 76``` 77 78to calculate `sin(x)` and `cos(x)`. It also uses the relation 79 80``` 81cos(x) = sin(x + pi/2) 82``` 83 84to calculate `cos(x)`. It has a complexity of `O(n^3)`. 85 86**Note**: this series has a tendency to *occasionally* produce an error of 1 87[ULP][6]. (It is an unfortunate side effect of the algorithm, and there isn't 88any way around it; [this article][7] explains why calculating sine and cosine, 89and the other transcendental functions below, within less than 1 ULP is nearly 90impossible and unnecessary.) Therefore, I recommend that users do their 91calculations with the precision (`scale`) set to at least 1 greater than is 92needed. 93 94### Exponentiation (`bc` Math Library Only) 95 96This `bc` uses the series 97 98``` 991 + x + x^2/2! + x^3/3! + ... 100``` 101 102to calculate `e^x`. Since this only works when `x` is small, it uses 103 104``` 105e^x = (e^(x/2))^2 106``` 107 108to reduce `x`. 109 110It has a complexity of `O(n^3)`. 111 112**Note**: this series can also produce errors of 1 ULP, so I recommend users do 113their calculations with the precision (`scale`) set to at least 1 greater than 114is needed. 115 116### Natural Logarithm (`bc` Math Library Only) 117 118This `bc` uses the series 119 120``` 121a + a^3/3 + a^5/5 + ... 122``` 123 124(where `a` is equal to `(x - 1)/(x + 1)`) to calculate `ln(x)` when `x` is small 125and uses the relation 126 127``` 128ln(x^2) = 2 * ln(x) 129``` 130 131to sufficiently reduce `x`. 132 133It has a complexity of `O(n^3)`. 134 135**Note**: this series can also produce errors of 1 ULP, so I recommend users do 136their calculations with the precision (`scale`) set to at least 1 greater than 137is needed. 138 139### Arctangent (`bc` Math Library Only) 140 141This `bc` uses the series 142 143``` 144x - x^3/3 + x^5/5 - x^7/7 + ... 145``` 146 147to calculate `atan(x)` for small `x` and the relation 148 149``` 150atan(x) = atan(c) + atan((x - c)/(1 + x * c)) 151``` 152 153to reduce `x` to small enough. It has a complexity of `O(n^3)`. 154 155**Note**: this series can also produce errors of 1 ULP, so I recommend users do 156their calculations with the precision (`scale`) set to at least 1 greater than 157is needed. 158 159### Bessel (`bc` Math Library Only) 160 161This `bc` uses the series 162 163``` 164x^n/(2^n * n!) * (1 - x^2 * 2 * 1! * (n + 1)) + x^4/(2^4 * 2! * (n + 1) * (n + 2)) - ... 165``` 166 167to calculate the bessel function (integer order only). 168 169It also uses the relation 170 171``` 172j(-n,x) = (-1)^n * j(n,x) 173``` 174 175to calculate the bessel when `x < 0`, It has a complexity of `O(n^3)`. 176 177**Note**: this series can also produce errors of 1 ULP, so I recommend users do 178their calculations with the precision (`scale`) set to at least 1 greater than 179is needed. 180 181### Modular Exponentiation (`dc` Only) 182 183This `dc` uses the [Memory-efficient method][8] to compute modular 184exponentiation. The complexity is `O(e*n^2)`, which may initially seem 185inefficient, but `n` is kept small by maintaining small numbers. In practice, it 186is extremely fast. 187 188### Non-Integer Exponentiation (`bc` Math Library 2 Only) 189 190This is implemented in the function `p(x,y)`. 191 192The algorithm used is to use the formula `e(y*l(x))`. 193 194It has a complexity of `O(n^3)` because both `e()` and `l()` do. 195 196### Rounding (`bc` Math Library 2 Only) 197 198This is implemented in the function `r(x,p)`. 199 200The algorithm is a simple method to check if rounding away from zero is 201necessary, and if so, adds `1e10^p`. 202 203It has a complexity of `O(n)` because of add. 204 205### Ceiling (`bc` Math Library 2 Only) 206 207This is implemented in the function `ceil(x,p)`. 208 209The algorithm is a simple add of one less decimal place than `p`. 210 211It has a complexity of `O(n)` because of add. 212 213### Factorial (`bc` Math Library 2 Only) 214 215This is implemented in the function `f(n)`. 216 217The algorithm is a simple multiplication loop. 218 219It has a complexity of `O(n^3)` because of linear amount of `O(n^2)` 220multiplications. 221 222### Permutations (`bc` Math Library 2 Only) 223 224This is implemented in the function `perm(n,k)`. 225 226The algorithm is to use the formula `n!/(n-k)!`. 227 228It has a complexity of `O(n^3)` because of the division and factorials. 229 230### Combinations (`bc` Math Library 2 Only) 231 232This is implemented in the function `comb(n,r)`. 233 234The algorithm is to use the formula `n!/r!*(n-r)!`. 235 236It has a complexity of `O(n^3)` because of the division and factorials. 237 238### Logarithm of Any Base (`bc` Math Library 2 Only) 239 240This is implemented in the function `log(x,b)`. 241 242The algorithm is to use the formula `l(x)/l(b)` with double the `scale` because 243there is no good way of knowing how many digits of precision are needed when 244switching bases. 245 246It has a complexity of `O(n^3)` because of the division and `l()`. 247 248### Logarithm of Base 2 (`bc` Math Library 2 Only) 249 250This is implemented in the function `l2(x)`. 251 252This is a convenience wrapper around `log(x,2)`. 253 254### Logarithm of Base 10 (`bc` Math Library 2 Only) 255 256This is implemented in the function `l10(x)`. 257 258This is a convenience wrapper around `log(x,10)`. 259 260### Root (`bc` Math Library 2 Only) 261 262This is implemented in the function `root(x,n)`. 263 264The algorithm is [Newton's method][9]. The initial guess is calculated as 265`10^ceil(length(x)/n)`. 266 267Like square root, its complexity is `O(log(n)*n^2)` as it requires one division 268per iteration, and it doubles the amount of correct digits per iteration. 269 270### Cube Root (`bc` Math Library 2 Only) 271 272This is implemented in the function `cbrt(x)`. 273 274This is a convenience wrapper around `root(x,3)`. 275 276### Greatest Common Divisor (`bc` Math Library 2 Only) 277 278This is implemented in the function `gcd(a,b)`. 279 280The algorithm is an iterative version of the [Euclidean Algorithm][10]. 281 282It has a complexity of `O(n^4)` because it has a linear number of divisions. 283 284This function ensures that `a` is always bigger than `b` before starting the 285algorithm. 286 287### Least Common Multiple (`bc` Math Library 2 Only) 288 289This is implemented in the function `lcm(a,b)`. 290 291The algorithm uses the formula `a*b/gcd(a,b)`. 292 293It has a complexity of `O(n^4)` because of `gcd()`. 294 295### Pi (`bc` Math Library 2 Only) 296 297This is implemented in the function `pi(s)`. 298 299The algorithm uses the formula `4*a(1)`. 300 301It has a complexity of `O(n^3)` because of arctangent. 302 303### Tangent (`bc` Math Library 2 Only) 304 305This is implemented in the function `t(x)`. 306 307The algorithm uses the formula `s(x)/c(x)`. 308 309It has a complexity of `O(n^3)` because of sine, cosine, and division. 310 311### Atan2 (`bc` Math Library 2 Only) 312 313This is implemented in the function `a2(y,x)`. 314 315The algorithm uses the [standard formulas][11]. 316 317It has a complexity of `O(n^3)` because of arctangent. 318 319[1]: https://en.wikipedia.org/wiki/Karatsuba_algorithm 320[2]: https://en.wikipedia.org/wiki/Long_division 321[3]: https://en.wikipedia.org/wiki/Exponentiation_by_squaring 322[4]: https://en.wikipedia.org/wiki/Newton%27s_method#Square_root_of_a_number 323[5]: https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method 324[6]: https://en.wikipedia.org/wiki/Unit_in_the_last_place 325[7]: https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT 326[8]: https://en.wikipedia.org/wiki/Modular_exponentiation#Memory-efficient_method 327[9]: https://en.wikipedia.org/wiki/Root-finding_algorithms#Newton's_method_(and_similar_derivative-based_methods) 328[10]: https://en.wikipedia.org/wiki/Euclidean_algorithm 329[11]: https://en.wikipedia.org/wiki/Atan2#Definition_and_computation 330