Fun math algorithms

[Albert Chan]

(09-06-2020 12:46 AM)Albert Chan Wrote:  Another trick is to scale denominator to 1 - ε, 1/(1-ε) = 1 + ε + ε² + ...

We can reduce summing terms, from O(n), to O(ln(n))

1/(1-ε) = 1 + ε + ε² + ... = (1+ε) + ε²(1+ε) + ε⁴[(1+ε) + ε²(1+ε)] + ...

For binary system, scaling by powers-of-2 can be done with shift in exponent.

If denominator is between 2/3 and 4/3, |ε| ≤ 1/3 → ε^32 < 5.4e-16
Thus, maximum of 5 sums is needed to reach IEEE double precision.

Example: 123/456 = (123/512) / (456/512) = 0.240234375 / 0.890625

>>> x, eps = 0.240234375, 1 - 0.890625
>>> x += x*eps; eps *= eps       # x = 0.266510009765625
>>> x += x*eps; eps *= eps       # x = 0.26969823986291885
>>> x += x*eps; eps *= eps       # x = 0.26973683658086722
>>> x += x*eps; eps *= eps       # x = 0.26973684210526305
>>> x += x*eps; eps *= eps       # x = 0.26973684210526316 = 123/456

We have quadratic convergence, but results is not self-correcting.
In other words, all calculations must be done with full precisions.

For arbitrary precision calculation of reciprocals of n, we can use Newton's method.
Note: we do not let f = 1/n - x, f' = -1, since the whole point is to avoid division.

f = n - 1/x       → f' = 1/x²       → x - f/f' = x + x*(1-n*x)

For 0.5 ≤ n < 1.0, guess x = 2.9-n-n work well, converged within 4 iterations.

Redo 123/456 as 0.240234375 * (1/0.890625)

>>> n = 0.890625
>>> x = 2.9 - n - n         # x = 1.1187499999999999
>>> x += x*(1-n*x)       # x = 1.1227923583984374
>>> x += x*(1-n*x)       # x = 1.1228070173524727
>>> x += x*(1-n*x)       # x = 1.1228070175438596 = 1/n
>>> 0.240234375 * x     # = 123/456
0.26973684210526316