1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
|
# Algorithms
This `bc` uses the math algorithms below:
### Addition
This `bc` uses brute force addition, which is linear (`O(n)`) in the number of
digits.
### Subtraction
This `bc` uses brute force subtraction, which is linear (`O(n)`) in the number
of digits.
### Multiplication
This `bc` uses two algorithms: [Karatsuba][1] and brute force.
Karatsuba is used for "large" numbers. ("Large" numbers are defined as any
number with `BC_NUM_KARATSUBA_LEN` digits or larger. `BC_NUM_KARATSUBA_LEN` has
a sane default, but may be configured by the user.) Karatsuba, as implemented in
this `bc`, is superlinear but subpolynomial (bounded by `O(n^log_2(3))`).
Brute force multiplication is used below `BC_NUM_KARATSUBA_LEN` digits. It is
polynomial (`O(n^2)`), but since Karatsuba requires both more intermediate
values (which translate to memory allocations) and a few more additions, there
is a "break even" point in the number of digits where brute force multiplication
is faster than Karatsuba. There is a script (`$ROOT/karatsuba.py`) that will
find the break even point on a particular machine.
***WARNING: The Karatsuba script requires Python 3.***
### Division
This `bc` uses Algorithm D ([long division][2]). Long division is polynomial
(`O(n^2)`), but unlike Karatsuba, any division "divide and conquer" algorithm
reaches its "break even" point with significantly larger numbers. "Fast"
algorithms become less attractive with division as this operation typically
reduces the problem size.
While the implementation of long division may appear to use the subtractive
chunking method, it only uses subtraction to find a quotient digit. It avoids
unnecessary work by aligning digits prior to performing subtraction and finding
a starting guess for the quotient.
Subtraction was used instead of multiplication for two reasons:
1. Division and subtraction can share code (one of the less important goals of
this `bc` is small code).
2. It minimizes algorithmic complexity.
Using multiplication would make division have the even worse algorithmic
complexity of `O(n^(2*log_2(3)))` (best case) and `O(n^3)` (worst case).
### Power
This `bc` implements [Exponentiation by Squaring][3], which (via Karatsuba) has
a complexity of `O((n*log(n))^log_2(3))` which is favorable to the
`O((n*log(n))^2)` without Karatsuba.
### Square Root
This `bc` implements the fast algorithm [Newton's Method][4] (also known as the
Newton-Raphson Method, or the [Babylonian Method][5]) to perform the square root
operation. Its complexity is `O(log(n)*n^2)` as it requires one division per
iteration.
### Sine and Cosine (`bc` Only)
This `bc` uses the series
```
x - x^3/3! + x^5/5! - x^7/7! + ...
```
to calculate `sin(x)` and `cos(x)`. It also uses the relation
```
cos(x) = sin(x + pi/2)
```
to calculate `cos(x)`. It has a complexity of `O(n^3)`.
**Note**: this series has a tendency to *occasionally* produce an error of 1
[ULP][6]. (It is an unfortunate side effect of the algorithm, and there isn't
any way around it; [this article][7] explains why calculating sine and cosine,
and the other transcendental functions below, within less than 1 ULP is nearly
impossible and unnecessary.) Therefore, I recommend that users do their
calculations with the precision (`scale`) set to at least 1 greater than is
needed.
### Exponentiation (`bc` Only)
This `bc` uses the series
```
1 + x + x^2/2! + x^3/3! + ...
```
to calculate `e^x`. Since this only works when `x` is small, it uses
```
e^x = (e^(x/2))^2
```
to reduce `x`. It has a complexity of `O(n^3)`.
**Note**: this series can also produce errors of 1 ULP, so I recommend users do
their calculations with the precision (`scale`) set to at least 1 greater than
is needed.
### Natural Logarithm (`bc` Only)
This `bc` uses the series
```
a + a^3/3 + a^5/5 + ...
```
(where `a` is equal to `(x - 1)/(x + 1)`) to calculate `ln(x)` when `x` is small
and uses the relation
```
ln(x^2) = 2 * ln(x)
```
to sufficiently reduce `x`. It has a complexity of `O(n^3)`.
**Note**: this series can also produce errors of 1 ULP, so I recommend users do
their calculations with the precision (`scale`) set to at least 1 greater than
is needed.
### Arctangent (`bc` Only)
This `bc` uses the series
```
x - x^3/3 + x^5/5 - x^7/7 + ...
```
to calculate `atan(x)` for small `x` and the relation
```
atan(x) = atan(c) + atan((x - c)/(1 + x * c))
```
to reduce `x` to small enough. It has a complexity of `O(n^3)`.
**Note**: this series can also produce errors of 1 ULP, so I recommend users do
their calculations with the precision (`scale`) set to at least 1 greater than
is needed.
### Bessel (`bc` Only)
This `bc` uses the series
```
x^n/(2^n * n!) * (1 - x^2 * 2 * 1! * (n + 1)) + x^4/(2^4 * 2! * (n + 1) * (n + 2)) - ...
```
to calculate the bessel function (integer order only).
It also uses the relation
```
j(-n,x) = (-1)^n * j(n,x)
```
to calculate the bessel when `x < 0`, It has a complexity of `O(n^3)`.
**Note**: this series can also produce errors of 1 ULP, so I recommend users do
their calculations with the precision (`scale`) set to at least 1 greater than
is needed.
### Modular Exponentiation (`dc` Only)
This `dc` uses the [Memory-efficient method][8] to compute modular
exponentiation. The complexity is `O(e*n^2)`, which may initially seem
inefficient, but `n` is kept small by maintaining small numbers. In practice, it
is extremely fast.
[1]: https://en.wikipedia.org/wiki/Karatsuba_algorithm
[2]: https://en.wikipedia.org/wiki/Long_division
[3]: https://en.wikipedia.org/wiki/Exponentiation_by_squaring
[4]: https://en.wikipedia.org/wiki/Newton%27s_method#Square_root_of_a_number
[5]: https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method
[6]: https://en.wikipedia.org/wiki/Unit_in_the_last_place
[7]: https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT
[8]: https://en.wikipedia.org/wiki/Modular_exponentiation#Memory-efficient_method
|
