Integer powers and n-th roots can be hundreds of floats away from the tightest bounds

[Jordan08]

GAOL computes pow(I, n) by repeated rounded products (uipow()), and nth_root(I, n) with mathlib's upow() and a rounded exponent 1/n. The bounds enclose the exact values, but can be several floats away from the tightest ones.

tests/arithmetic.cpp compares them exactly on 5000 random doubles. Largest number of floats between GAOL's bounds and the tightest ones, over all the platforms of the continuous integration:

	exponents from -30 to 30	any double
pow([a], n) for n = 3, 4, 5, 6, 7	2, 3, 5, 6, 8	2, 3, 5, 6, 7
pow([a], -n) for n = 2, 3	2, 4	3, 4 (when a^n is a double)
nth_root([a], n) for n = 3, 5, 6, 7	8, 4, 4, 4	234, 106, 117, 127
nth_root([a], n) for n = 2, 4, and sqrt([a])	1	1

pow([a], -n) is computed as 1/a^n. When a^n is beyond the largest double, the result is [0, 1/max]: pow([0x1.346d964e3694dp+536], -2) = [0, 0x0.4000000000001p-1022], while the exact value is near 2^-1072.

The tests allow about twice these distances. Two possible improvements:

• Roots could be corrected as the square root now is (3257c31): check each bound with exact products, and move it by a float when it is not a bound yet.
• Small integer powers could be computed with exact products (double-double).

[Jordan08]

Negative powers fixed when a^n overflows

pow(I, -n) is now computed as (1/I)^n when I^n is beyond the largest double (17f7662). Elsewhere it is still 1/I^n, which has a single rounded operation when I^n is exact.

	before	now
pow([0x1.346d964e3694dp+536], -2) (the example)	[0, 0x0.4000000000001p-1022]	[0x0.0000000000002p-1022, 0x0.0000000000003p-1022]
pow([10], -400)	[0, 5.6e-309]	[0, 2^-1074]
pow([2], -1050)	[0, 5.6e-309]	[2^-1050]

• Tests. tests/arithmetic.cpp now checks these powers within 8 floats, like the other negative powers, rather than only checking that they enclose the exact value. It also requires exact results on special cases, such as the negative powers of intervals containing 0 or of [1, +oo].
• Roots of 0, 1 and -1 are now exact (7d99807): nth_root([1], 3) was [1 - 2^-53, 1 + 2^-52].

Otherwise the distances in the table are unchanged: nth_root([27], 3) is still [0x1.7ffffffffffffp+1, 0x1.8000000000002p+1]. doc/accuracy.md (0536f71) documents them.

[Jordan08]

Where this stands (master, 17 September)

Unchanged since the last comment for the distances of the table. What is new is about upow(), which nth_root() takes its values from: it was compared with mpmath at 240000 arguments (results near the overflow, subnormal results, bases next to 1 with large exponents, exact powers and their neighbouring doubles, which take its slow path) and was correctly rounded at all of them. pow(I, J) now relies on it too (#8).

So the distances of nth_root() (up to 234 doubles over all doubles) come from the exponent alone: 1/n rounded to a double moves the root by about |log x|·2^-53/n relatively.

Proposals

• Roots corrected with integer powers, as the square root is with a product (3257c31). r being upow(x, 1/n):
  • a lower bound l is proved by pow_up(l, n) <= x (the integer power rounded upward, which GAOL has: ipow_up()), an upper bound u by pow_dn(u, n) >= x;
  • start from the doubles around r, move outward while the proof fails, then inward while it still holds.
  • The integer power of n − 1 rounded products is within (n − 1)·2^-53 relatively, which divided by n is less than one double of the root: the bounds would be within about two doubles of the tightest for every n and every double, instead of 8 to 234. Cost: a few integer powers, log2(n) products each.
• Small integer powers with exact products. x³ to x⁷ computed in double-double (a product and its error, with fma() where the processor has one, Dekker's product elsewhere), rounded once at the end, would be the tightest or one double wider, instead of 2 to 8. About twice the products of now.

Neither is done: both change the time of pow(x, n) and nth_root(), which the benchmark of doc/compare measures (14.3 ns for pow(x, 3)).