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/karatsuba.py`) that will 29find 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. Its complexity is `O(log(n)*n^2)` as it requires one division per 66iteration. 67 68### Sine and Cosine (`bc` Only) 69 70This `bc` uses the series 71 72``` 73x - x^3/3! + x^5/5! - x^7/7! + ... 74``` 75 76to calculate `sin(x)` and `cos(x)`. It also uses the relation 77 78``` 79cos(x) = sin(x + pi/2) 80``` 81 82to calculate `cos(x)`. It has a complexity of `O(n^3)`. 83 84**Note**: this series has a tendency to *occasionally* produce an error of 1 85[ULP][6]. (It is an unfortunate side effect of the algorithm, and there isn't 86any way around it; [this article][7] explains why calculating sine and cosine, 87and the other transcendental functions below, within less than 1 ULP is nearly 88impossible and unnecessary.) Therefore, I recommend that users do their 89calculations with the precision (`scale`) set to at least 1 greater than is 90needed. 91 92### Exponentiation (`bc` Only) 93 94This `bc` uses the series 95 96``` 971 + x + x^2/2! + x^3/3! + ... 98``` 99 100to calculate `e^x`. Since this only works when `x` is small, it uses 101 102``` 103e^x = (e^(x/2))^2 104``` 105 106to reduce `x`. It has a complexity of `O(n^3)`. 107 108**Note**: this series can also produce errors of 1 ULP, so I recommend users do 109their calculations with the precision (`scale`) set to at least 1 greater than 110is needed. 111 112### Natural Logarithm (`bc` Only) 113 114This `bc` uses the series 115 116``` 117a + a^3/3 + a^5/5 + ... 118``` 119 120(where `a` is equal to `(x - 1)/(x + 1)`) to calculate `ln(x)` when `x` is small 121and uses the relation 122 123``` 124ln(x^2) = 2 * ln(x) 125``` 126 127to sufficiently reduce `x`. It has a complexity of `O(n^3)`. 128 129**Note**: this series can also produce errors of 1 ULP, so I recommend users do 130their calculations with the precision (`scale`) set to at least 1 greater than 131is needed. 132 133### Arctangent (`bc` Only) 134 135This `bc` uses the series 136 137``` 138x - x^3/3 + x^5/5 - x^7/7 + ... 139``` 140 141to calculate `atan(x)` for small `x` and the relation 142 143``` 144atan(x) = atan(c) + atan((x - c)/(1 + x * c)) 145``` 146 147to reduce `x` to small enough. It has a complexity of `O(n^3)`. 148 149**Note**: this series can also produce errors of 1 ULP, so I recommend users do 150their calculations with the precision (`scale`) set to at least 1 greater than 151is needed. 152 153### Bessel (`bc` Only) 154 155This `bc` uses the series 156 157``` 158x^n/(2^n * n!) * (1 - x^2 * 2 * 1! * (n + 1)) + x^4/(2^4 * 2! * (n + 1) * (n + 2)) - ... 159``` 160 161to calculate the bessel function (integer order only). 162 163It also uses the relation 164 165``` 166j(-n,x) = (-1)^n * j(n,x) 167``` 168 169to calculate the bessel when `x < 0`, It has a complexity of `O(n^3)`. 170 171**Note**: this series can also produce errors of 1 ULP, so I recommend users do 172their calculations with the precision (`scale`) set to at least 1 greater than 173is needed. 174 175### Modular Exponentiation (`dc` Only) 176 177This `dc` uses the [Memory-efficient method][8] to compute modular 178exponentiation. The complexity is `O(e*n^2)`, which may initially seem 179inefficient, but `n` is kept small by maintaining small numbers. In practice, it 180is extremely fast. 181 182[1]: https://en.wikipedia.org/wiki/Karatsuba_algorithm 183[2]: https://en.wikipedia.org/wiki/Long_division 184[3]: https://en.wikipedia.org/wiki/Exponentiation_by_squaring 185[4]: https://en.wikipedia.org/wiki/Newton%27s_method#Square_root_of_a_number 186[5]: https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method 187[6]: https://en.wikipedia.org/wiki/Unit_in_the_last_place 188[7]: https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT 189[8]: https://en.wikipedia.org/wiki/Modular_exponentiation#Memory-efficient_method 190