What happens
PowByLogarithm (SignificantNumber/SignificantNumber.cs:642) starts with
double logValue = Math.Log(Math.Abs(value.To<double>()));
Every non-integer power goes through this line, and so does every integer power above MaxExactIntegerPower (1024, line 744). Rounding the base to a double leaves a relative error of about 1e-16 in x. Near 1, ln x ≈ x − 1, so that becomes a relative error of about 1e-16 / (x − 1) in ln x, and it carries straight into the result. The result is then labelled with the base's significant digits, up to 15.
Repro (reference values from Python decimal at prec=40):
SignificantNumber S(string s) => SignificantNumber.Parse(s, NumberStyles.Float, CultureInfo.InvariantCulture);
S("1.00000000000001").Pow(S("1E14")) // 2.71611003408702 (SignificantDigits 15) true 2.71828182845903
S("1.0000000000001").Pow(S("1E13")) // 2.7161100340869 (14) true 2.71828182845891
S("1.00000000000001").Pow(S("100000000000000.5")) // 2.71611003408704 (15) true 2.71828182845905
S("1.00000000000001").Pow(S("1024")) // 1.00000000001024 exact path, correct
S("1.00000000000001").Pow(S("1025")) // 1.00000000001024 (15) true 1.00000000001025
S("1.0000000000000000001").Pow(S("1E19")) // 1 true 2.71828…
At the 1024/1025 boundary, the result goes from exact to wrong in the 12th digit, and x^1025 == x^1024 for x > 1.
Why it matters
Compounding a small rate many times is an ordinary use of Pow. The library's core promise is that the reported significant digits are trustworthy. Here the result carries 15 claimed digits of which only 2 or 3 are right.
This is distinct from #113. #113 covers Math.Exp amplifying error when |p·ln x| is in the hundreds. Here p·ln x ≈ 1, and the error comes from ln x itself. The cap proposed in #113 would still report 14–15 digits for these inputs.
Suggested fix
- Compute
ln|x| without rounding x to a double first:
- When
|x| is near 1, form d = |x| − 1 exactly in PreciseNumber and use double.LogP1((double)d). Relative precision in d is kept.
- Otherwise use
ln(significand in [1,10)) + exponent·ln 10, in the style of Log10OfMagnitude.
- In addition, bound the reported digits by the conditioning of the input, roughly
min(significantDigits, 15 − log10(1 / |x − 1|)), so the result never claims more digits than a double pipeline can deliver.
Acceptance criteria
- The cases above match the reference values to the number of significant digits reported, or report fewer digits.
x^1025 / x^1024 == x holds for x = 1.00000000000001 to the reported precision.
- Tests cover a base within 1e-10 of 1 on both the integer path (power above 1024) and the non-integer path.
What happens
PowByLogarithm(SignificantNumber/SignificantNumber.cs:642) starts withEvery non-integer power goes through this line, and so does every integer power above
MaxExactIntegerPower(1024, line 744). Rounding the base to adoubleleaves a relative error of about 1e-16 inx. Near 1,ln x ≈ x − 1, so that becomes a relative error of about1e-16 / (x − 1)inln x, and it carries straight into the result. The result is then labelled with the base's significant digits, up to 15.Repro (reference values from Python
decimalat prec=40):At the 1024/1025 boundary, the result goes from exact to wrong in the 12th digit, and
x^1025 == x^1024forx > 1.Why it matters
Compounding a small rate many times is an ordinary use of
Pow. The library's core promise is that the reported significant digits are trustworthy. Here the result carries 15 claimed digits of which only 2 or 3 are right.This is distinct from #113. #113 covers
Math.Expamplifying error when |p·ln x| is in the hundreds. Here p·ln x ≈ 1, and the error comes fromln xitself. The cap proposed in #113 would still report 14–15 digits for these inputs.Suggested fix
ln|x|without roundingxto a double first:|x|is near 1, formd = |x| − 1exactly inPreciseNumberand usedouble.LogP1((double)d). Relative precision indis kept.ln(significand in [1,10)) + exponent·ln 10, in the style ofLog10OfMagnitude.min(significantDigits, 15 − log10(1 / |x − 1|)), so the result never claims more digits than a double pipeline can deliver.Acceptance criteria
x^1025 / x^1024 == xholds forx = 1.00000000000001to the reported precision.