Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 25 additions & 0 deletions BREAKING-CHANGES.md
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,31 @@ and the complex arcsine's real part moved two units away. A caller comparing a h
result to a pinned string will see the last digits differ; a caller comparing to a tolerance
will not. [#1338](https://github.com/asc-community/AngouriMath/issues/1338).

### The last digits of an evaluated logarithm, exponential, real power or complex power

**Silent, and in the hundredth digit.** `ln`, `log`, `e^x`, a real base to a real power and the
complex power and inverse trigonometric functions that are built on them are series in fixed
point now rather than PeterO's `Log`, `Exp` and `Pow`, and `e^x` is the exponential of `x` rather
than the constant's hundred digits raised to `x`. Measured on `eb2a5ac2` and on this change,
against mpmath at a hundred and forty digits, in units of the last of a hundred places:

| | Was | Is |
|---|---|---|
| `ln(pi)` | 5.3 off | 0.3 |
| `log(-2, 3 + 2i)` | 7.7 in the real part, 33.5 in the imaginary | 0.7 and 0.5 |
| `arcsin`, `arccos` of `3`, `4`, `i`, `3 ± 2i` | 5 to 28 off in the imaginary part | within 1.4 |
| `arctan(3 + 2i)` | 33 off in the imaginary part | 3.1 |
| `2^i` | 1.3 and 1.2 | 0.3 and 0.2 |
| `(3 ± 2i)^(3 ± 2i)` | 22 to 289 off | 0.5 to 5.6 |
| `e^700` | 97 correct digits | 100 |

Nineteen of the values pinned in `NumericDigits` moved, every one towards the true value. A caller
comparing a hundred-digit result to a pinned string will see the last digits differ; a caller
comparing to a tolerance will not. A value within ten to the minus fifty of an integer -- half
the working digits, or `MathS.Settings.PrecisionErrorZeroRange` where a caller has set it -- is
still that integer, so `e^(-123.456)` is 0 as it was.
[#1338](https://github.com/asc-community/AngouriMath/issues/1338).

### A quotient of polynomials is split into coprime blocks before a root is peeled off

`SolveByPartialFractions` tried `TrySplitOffRationalRoot` before `TrySplitIntoCoprimeParts`. Both
Expand Down
62 changes: 48 additions & 14 deletions Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs
Original file line number Diff line number Diff line change
Expand Up @@ -386,7 +386,7 @@ public static Complex Sqrt(Complex num)
public static Complex Exp(Complex x)
{
var context = MathS.Settings.DecimalPrecisionContext;
var expReal = x.RealPart.EDecimal.Exp(context);
var expReal = x.RealPart.EDecimal.Exponential(context);
var imaginary = x.ImaginaryPart.EDecimal;
if (imaginary.IsZero)
return expReal;
Expand Down Expand Up @@ -418,6 +418,13 @@ static Complex BinaryIntPow(Complex num, EInteger val)
var squared = half * half;
return divRem[1].IsZero ? squared : squared * BinaryIntPow(num, divRem[1]);
}
// e to a real power is the exponential of the power: raising the constant's
// hundred digits loses to the power what it multiplies their error by --
// e^700 came back to ninety-seven digits -- and the exponential is a series
// where the power was a logarithm and a series.
if (@base is Real { EDecimal: var maybeE } && power is Real { EDecimal: var realExponent }
&& maybeE.Equals(InternalAMExtensions.ConstantCache.Lookup(MathS.Settings.DecimalPrecisionContext).E))
return Real.Create(realExponent.Exponential(MathS.Settings.DecimalPrecisionContext));
// TODO: make it more detailed (e. g. +oo ^ +oo = +oo)
if (@base.IsFinite && power is Integer { EInteger: var pow })
{
Expand Down Expand Up @@ -456,7 +463,7 @@ static Complex BinaryIntPow(Complex num, EInteger val)

var context = MathS.Settings.DecimalPrecisionContext;
if (@base is Real { EDecimal: { IsNegative: false } realBase } && power is Real { EDecimal: var realPower })
return realBase.Pow(realPower, context);
return PowerOfReals(realBase, realPower, context);
// From https://source.dot.net/#System.Runtime.Numerics/System/Numerics/Complex.cs,7dc9c2ee4f99814a
// NOTE: System.Numerics.Complex.Pow(0, System.Numerics.Complex(-2, 1)) gives 0 + 0i despite being mathematically undefined
// NOTE: System.Numerics.Complex.Pow(0, 0) gives 1 + 0i despite being mathematically undefined
Expand All @@ -472,8 +479,8 @@ static Complex BinaryIntPow(Complex num, EInteger val)

var rho = @base.Abs().EDecimal;
var theta = baseImaginary.Arctan2(baseReal, context);
var newRho = powerReal.MultiplyAndAdd(theta, powerImaginary.Multiply(rho.Log(context), context), context);
var t = rho.Pow(powerReal, context).Multiply(powerImaginary.Multiply(-theta, context).Exp(context), context);
var newRho = powerReal.MultiplyAndAdd(theta, powerImaginary.Multiply(rho.NaturalLogarithm(context), context), context);
var t = PowerOfReals(rho, powerReal, context).Multiply(powerImaginary.Multiply(-theta, context).Exponential(context), context);

return Complex.Create(t.Multiply(newRho.Cos(context), context), t.Multiply(newRho.Sin(context), context));
}
Expand All @@ -493,11 +500,38 @@ public static Complex Log(Complex @base, Complex x)
if (x is Real real && real.EDecimal.CompareTo(EDecimal.Zero) > 0
&& @base is Real realBase && realBase.EDecimal.CompareTo(EDecimal.Zero) > 0
&& realBase.EDecimal.CompareTo(EDecimal.One) != 0)
return real.EDecimal.LogN(realBase.EDecimal, MathS.Settings.DecimalPrecisionContext);
return LogOfReals(real.EDecimal, realBase.EDecimal, MathS.Settings.DecimalPrecisionContext);
// From https://source.dot.net/#System.Runtime.Numerics/System/Numerics/Complex.cs,cf15f2e5cc49cef1
return Ln(x) / Ln(@base);
}

/// <summary>
/// <c>base^power</c> for a nonnegative real base and a real power, as
/// <c>exp(power ln base)</c> with the logarithm and the product carrying eight
/// guard digits. Zero, an infinity and NaN are PeterO's answers.
/// </summary>
private static EDecimal PowerOfReals(EDecimal @base, EDecimal power, EContext context)
{
if (@base.IsZero || !@base.IsFinite || !power.IsFinite)
return @base.Pow(power, context);
var working = InternalAMExtensions.WithGuardDigits(context, 8);
return power.Multiply(@base.NaturalLogarithm(working), working).Exponential(working).RoundToPrecision(context);
}

/// <summary>
/// <c>log_base(x)</c> for positive reals as <c>ln x / ln base</c>, with the base
/// <c>e</c> -- the value the constant evaluates to in this context -- recognised, so
/// that <c>ln(x)</c>, which arrives here as a logarithm to the base <c>e</c>, is one
/// logarithm and not two and a division.
/// </summary>
private static EDecimal LogOfReals(EDecimal x, EDecimal @base, EContext context)
{
var log = x.NaturalLogarithm(context);
if (@base.Equals(InternalAMExtensions.ConstantCache.Lookup(context).E))
return log;
return log.Divide(@base.NaturalLogarithm(context), context);
}

/// <summary>
/// Whether a number is exactly known as a ratio, yet its
/// <see cref="Real.EDecimal"/> -- which is rounded into
Expand Down Expand Up @@ -528,12 +562,12 @@ private static EDecimal LnOfEInteger(EInteger n, EContext context)
var keep = (context.Precision.IsZero ? EInteger.FromInt32(100) : context.Precision).Add(10);
var digits = n.GetDigitCountAsEInteger();
if (digits.CompareTo(keep) <= 0)
return EDecimal.FromEInteger(n).Log(context);
return EDecimal.FromEInteger(n).NaturalLogarithm(context);
var shift = digits.Subtract(keep);
// Create(n, -shift) is n with its decimal point moved, not a division:
// the digits are untouched, only the exponent comes back into range.
return EDecimal.Create(n, shift.Negate()).Log(context)
.Add(EDecimal.FromInt32(10).Log(context)
return EDecimal.Create(n, shift.Negate()).NaturalLogarithm(context)
.Add(EDecimal.FromInt32(10).NaturalLogarithm(context)
.Multiply(EDecimal.FromEInteger(shift), context), context);
}

Expand All @@ -544,9 +578,9 @@ public static Complex Ln(Complex x)
if (LostToExponentRange(x) && x is Rational { ERational: { Sign: > 0 } ratio })
return LnOfRatio(ratio, context);
if (x is Real { EDecimal: { IsNegative: false } real })
return real.Log(context);
return real.NaturalLogarithm(context);
// From https://source.dot.net/#System.Runtime.Numerics/System/Numerics/Complex.cs,cf15f2e5cc49cef1
return Complex.Create(x.Abs().EDecimal.Log(context), x.ImaginaryPart.EDecimal.Arctan2(x.RealPart.EDecimal, context));
return Complex.Create(x.Abs().EDecimal.NaturalLogarithm(context), x.ImaginaryPart.EDecimal.Arctan2(x.RealPart.EDecimal, context));
}

/// <summary>Calculates the exact value of sine of num</summary>
Expand All @@ -560,7 +594,7 @@ public static Complex Sin(Complex num)
// We need both sinh and cosh of imaginary part.
// To avoid multiple calls to Exp with the same value,
// we compute them both here from a single call to Exp.
var p = im.Exp(context);
var p = im.Exponential(context);
var q = EDecimal.One.Divide(p, context);
var sinh = p.Subtract(q, context).Divide(2, context);
var cosh = p.Add(q, context).Divide(2, context);
Expand Down Expand Up @@ -621,7 +655,7 @@ public static Complex Cos(Complex num)
return real.Cos(context);
// From https://source.dot.net/#System.Runtime.Numerics/System/Numerics/Complex.cs,cf15f2e5cc49cef1
var (re, im) = (num.RealPart.EDecimal, num.ImaginaryPart.EDecimal);
var p = im.Exp(context);
var p = im.Exponential(context);
var q = EDecimal.One.Divide(p, context);
var sinh = p.Subtract(q, context).Divide(2, context);
var cosh = p.Add(q, context).Divide(2, context);
Expand Down Expand Up @@ -649,7 +683,7 @@ public static Complex Tan(Complex num)

var x2 = num.RealPart.EDecimal.Multiply(2, context);
var y2 = num.ImaginaryPart.EDecimal.Multiply(2, context);
var p = y2.Exp(context);
var p = y2.Exponential(context);
var q = EDecimal.One.Divide(p, context);
var cosh = p.Add(q, context).Divide(2, context);
if (num.ImaginaryPart.EDecimal.Abs().LessThanOrEquals(4))
Expand Down Expand Up @@ -707,7 +741,7 @@ public static Complex Cotan(Complex num)
var sigma = xm1.MultiplyAndAdd(xm1, y.Multiply(y, context), context).Sqrt(context);
var alpha = rho.Add(sigma, context).Divide(2, context);
return (rho.Subtract(sigma, context).Divide(2, context),
alpha.MultiplyAndSubtract(alpha, EDecimal.One, context).Sqrt(context).Add(alpha, context).Log(context).Multiply((y.IsNegative || y.IsZero) ? -1 : 1, context));
alpha.MultiplyAndSubtract(alpha, EDecimal.One, context).Sqrt(context).Add(alpha, context).NaturalLogarithm(context).Multiply((y.IsNegative || y.IsZero) ? -1 : 1, context));
}

/// <summary>Calculates the exact value of arcsine of num</summary>
Expand Down
37 changes: 37 additions & 0 deletions Sources/AngouriMath/Docs/WhatsNew/version_performance_control.md
Original file line number Diff line number Diff line change
Expand Up @@ -256,6 +256,43 @@ cent — its measured run-to-run spread on mean time is up to 51.8%. A move of 8
outside that band by a wide margin and agrees in sign and rough size with the allocation column
beside it. The small rows from the same run are still not worth reading, and are not quoted.

## The 2054th: the logarithm and the exponential as series in fixed point

`EvalTranscendentalFresh` -- `ln(1 + x^2) arctan(x)` at a fresh point -- was 1.3 ms after the
2050th, and a `ln` alone 0.7 ms: PeterO's `Log` is seven hundred microseconds at a hundred
digits and its `Exp` three hundred, and `ln(x)` arrived as `log_e(x)`, which was two of them
and a division. A series in `EDecimal` pays for each term an alignment of exponents, an exact
sum and a rounding, two to three microseconds a term at a hundred digits, which is what those
are made of.

The logarithm and the exponential are their own series now, run in fixed point: an `EInteger`
holding the value times `2^bits`, with the context's digits in bits and forty-eight over, where
a product is a big-integer multiply and a shift, a division by a term's index an integer
division and a sum an addition -- a third of a microsecond a term -- and the value is a decimal
again exactly, by `5^bits` and a move of the point. `ln x` is `2 artanh((m - 1)/(m + 1))` with
`m` the mantissa brought between `1/sqrt(2)` and `sqrt(2)` by powers of ten and two, seventy
terms; `e^x` is `x = k ln 2 + r`, `r` halved ten times, twenty-five terms, ten squarings and the
power of two back. `ln 2` and `ln 10` are in the constant cache, per context. `log_b(x)` is
`ln x / ln b` with the base `e` recognised, a nonnegative real to a real power is
`exp(power ln base)` with eight guard digits, and `e` to any real power is the exponential of
the power -- `e^700` came back to ninety-seven digits as the constant's hundred raised to the
seven hundredth. The hyperbolic functions and the complex logarithm, exponential, power and
inverse trigonometric functions go through the same two.

| benchmark | 2053rd | 2054th | allocation | time |
|---|--:|--:|--:|--:|
| `EvalTranscendentalFresh` | 2,246,889 | **423,032** | **−81.2%** | 1.33 ms → **245 µs** |
| `SimplifyHard` | 302,423,792 | **189,726,120** | **−37.3%** | 172 → **109 ms** |
| `SolveEasy` | 6,284,124 | **4,425,194** | **−29.6%** | 3.55 → **2.88 ms** |
| every other entry | | | within 1.5% | |

Bytes allocated per call, same machine, both columns by the gate in one session; the gate's
baseline was taken from this run. The simplifier and the solver move because they evaluate
exponentials and logarithms numerically underneath, to compare candidates. Alone: `ln(x)` at
a hundred digits is 35 µs where it was 685, `e^x` 30 where it was 740, `sinh(x)` 62 where it
was 1,570, and at five hundred digits `ln(x)` is 2.9 ms where it was 15.4. Nineteen pinned
hundred-digit values moved, every one towards mpmath's; `BREAKING-CHANGES.md` has the table.

## The 2053rd: equality without the caches, and a rational's decimal written out when asked

Two things found under the evaluator's remaining cost, each measured alone and each reaching
Expand Down
Loading
Loading