diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 2c4e2e5a7..78b746d20 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -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 diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs b/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs index 4b3a3a993..6b9773958 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Number/Operators.cs @@ -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; @@ -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 }) { @@ -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 @@ -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)); } @@ -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); } + /// + /// base^power for a nonnegative real base and a real power, as + /// exp(power ln base) with the logarithm and the product carrying eight + /// guard digits. Zero, an infinity and NaN are PeterO's answers. + /// + 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); + } + + /// + /// log_base(x) for positive reals as ln x / ln base, with the base + /// e -- the value the constant evaluates to in this context -- recognised, so + /// that ln(x), which arrives here as a logarithm to the base e, is one + /// logarithm and not two and a division. + /// + 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); + } + /// /// Whether a number is exactly known as a ratio, yet its /// -- which is rounded into @@ -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); } @@ -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)); } /// Calculates the exact value of sine of num @@ -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); @@ -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); @@ -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)) @@ -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)); } /// Calculates the exact value of arcsine of num diff --git a/Sources/AngouriMath/Docs/WhatsNew/version_performance_control.md b/Sources/AngouriMath/Docs/WhatsNew/version_performance_control.md index f12b3ff85..834cbfa55 100644 --- a/Sources/AngouriMath/Docs/WhatsNew/version_performance_control.md +++ b/Sources/AngouriMath/Docs/WhatsNew/version_performance_control.md @@ -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 diff --git a/Sources/AngouriMath/Functions/InternalAMExtensions.cs b/Sources/AngouriMath/Functions/InternalAMExtensions.cs index a75b8d633..b294118fe 100644 --- a/Sources/AngouriMath/Functions/InternalAMExtensions.cs +++ b/Sources/AngouriMath/Functions/InternalAMExtensions.cs @@ -198,9 +198,31 @@ public static ConstantCache Lookup(EContext context) HalfPi = Pi.Multiply(Half, context); QuarterPi = HalfPi.Multiply(Half, context); E = EDecimal.One.Exp(context); + // The fixed point the logarithm's and the exponential's series run in: the + // context's digits in bits, and forty-eight bits over for the series' own + // truncations and the exponential's ten squarings. + FixedBits = (int)(context.Precision.ToInt32Checked() * 3.32192809488736 + 48); + FixedOne = EInteger.One.ShiftLeft(FixedBits); + FiveToFixedBits = EInteger.FromInt32(5).Pow(FixedBits); + FixedSqrt2 = ToFixed(EDecimal.FromString("1.41421356237309504880168872420969807856967187537694"), FixedBits); + // ln 2 = 2 artanh(1/3), and ln 10 = ln 8 + ln 5/4 = 3 ln 2 + 2 artanh(1/9). + FixedLn2 = TwiceArtanh(FixedOne.Divide(3), FixedBits); + FixedLn10 = FixedLn2.Multiply(3).Add(TwiceArtanh(FixedOne.Divide(9), FixedBits)); } /// Represents public EDecimal Pi { get; } + /// The bits after the point of the fixed-point numbers below + public int FixedBits { get; } + /// One, in fixed point: 2 to the + public EInteger FixedOne { get; } + /// 5 to the , by which a fixed-point number is a decimal exactly + public EInteger FiveToFixedBits { get; } + /// The square root of 2, in fixed point, to fifty digits + public EInteger FixedSqrt2 { get; } + /// The natural logarithm of 2, in fixed point + public EInteger FixedLn2 { get; } + /// The natural logarithm of 10, in fixed point + public EInteger FixedLn10 { get; } /// Represents 2 * public EDecimal TwoPi { get; } /// Represents / 2 @@ -345,13 +367,20 @@ public static EDecimal Cos(this EDecimal x, EContext context) /// context instance, so a working context built afresh on every call would have pi /// recomputed on every call, and would leave an entry behind each time. /// - private static EContext WithGuardDigits(EContext context, int digits) + internal static EContext WithGuardDigits(EContext context, int digits) { - var table = digits == 8 ? eightMoreDigits : fiveMoreDigits; + var table = digits switch + { + 5 => fiveMoreDigits, + 8 => eightMoreDigits, + 12 => twelveMoreDigits, + _ => throw new AngouriBugException("A guard-digit count without its table"), + }; return table.GetValue(context, c => c.WithPrecision(c.Precision.ToInt32Checked() + digits)); } [ConstantField] private static readonly System.Runtime.CompilerServices.ConditionalWeakTable eightMoreDigits = new(); [ConstantField] private static readonly System.Runtime.CompilerServices.ConditionalWeakTable fiveMoreDigits = new(); + [ConstantField] private static readonly System.Runtime.CompilerServices.ConditionalWeakTable twelveMoreDigits = new(); /// /// The sine and the cosine of together, each to the precision of @@ -642,10 +671,168 @@ public static EDecimal Arctan2(this EDecimal y, EDecimal x, EContext context) } + // The logarithm and the exponential are series, and 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 why PeterO's Log and Exp are + // three hundred to seven hundred microseconds. The series here run in fixed point: + // an EInteger holding the value times 2^FixedBits, 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 multiplying with 5^FixedBits and moving the point. + // https://github.com/asc-community/AngouriMath/issues/1338 + + /// times 2^, truncated to an integer. + private static EInteger ToFixed(EDecimal x, int bits) + { + var exponent = x.Exponent.ToInt32Checked(); + var mantissa = x.Mantissa; + return exponent >= 0 + ? mantissa.Multiply(EInteger.FromInt32(10).Pow(exponent)).ShiftLeft(bits) + : mantissa.ShiftLeft(bits).Divide(EInteger.FromInt32(10).Pow(-exponent)); + } + + /// A fixed-point as a decimal, exactly, rounded to . + private static EDecimal FromFixed(EInteger n, ConstantCache consts, EContext working) + => EDecimal.Create(n.Multiply(consts.FiveToFixedBits), -consts.FixedBits).RoundToPrecision(working); + + /// + /// The product of two fixed-point numbers, truncated towards zero -- a shift alone + /// floors, and a negative term of a series floored never reaches zero. + /// + private static EInteger MultiplyFixed(EInteger a, EInteger b, int bits) + { + var product = a.Multiply(b); + return product.Sign < 0 ? product.Negate().ShiftRight(bits).Negate() : product.ShiftRight(bits); + } + + /// + /// 2 artanh(y) as its series 2 (y + y^3/3 + y^5/5 + ...), which is + /// ln((1 + y)/(1 - y)), in fixed point; for |y| below 0.18, + /// seventy terms at a hundred digits. + /// + private static EInteger TwiceArtanh(EInteger y, int bits) + { + var square = MultiplyFixed(y, y, bits); + var term = y; + var sum = y; + for (var k = 1; k < 100000; k++) + { + term = MultiplyFixed(term, square, bits); + if (term.IsZero) + break; + sum = sum.Add(term.Divide(2 * k + 1)); + } + return sum.ShiftLeft(1); + } + + [ConstantField] private static readonly EDecimal sqrt2 = EDecimal.FromString("1.4142135623730950488"); + [ConstantField] private static readonly EDecimal halfSqrt2 = EDecimal.FromString("0.70710678118654752440"); + + /// + /// The natural logarithm of to the precision of + /// : the argument written as a mantissa times a power of ten + /// and of two, the mantissa brought between 1/sqrt(2) and sqrt(2), and + /// 2 artanh((m - 1)/(m + 1)) there, a series in a square below 0.03, in + /// fixed point. PeterO's is seven hundred + /// microseconds at a hundred digits; this is about twenty. An argument within + /// [1/sqrt(2), sqrt(2)] goes to the series as it is, so a value near 1 keeps + /// every digit rather than losing them to ln 10 - 3 ln 2 - ... cancelling. + /// Zero, a negative, an infinity and NaN are PeterO's answers. + /// https://github.com/asc-community/AngouriMath/issues/1338 + /// + public static EDecimal NaturalLogarithm(this EDecimal x, EContext context) + { + if (!x.IsFinite || x.IsNegative || x.IsZero) + return x.Log(context); + var working = WithGuardDigits(context, 8); + var consts = ConstantCache.Lookup(working); + var bits = consts.FixedBits; + EInteger mantissa; + var tens = 0; + var twos = 0; + if (x.CompareTo(sqrt2) <= 0 && x.CompareTo(halfSqrt2) >= 0) + mantissa = ToFixed(x, bits); + else + { + // m * 10^tens with 1 <= m < 10, by moving the point: the digits are untouched; + // then halved, exactly, until it is at most sqrt(2). + tens = x.Exponent.Add(x.Precision()).Subtract(1).ToInt32Checked(); + mantissa = ToFixed(x.MovePointLeft(tens), bits); + while (mantissa.CompareTo(consts.FixedSqrt2) > 0) + { + mantissa = mantissa.ShiftRight(1); + twos++; + } + } + var y = mantissa.Subtract(consts.FixedOne).ShiftLeft(bits).Divide(mantissa.Add(consts.FixedOne)); + var log = TwiceArtanh(y, bits); + if (twos != 0) + log = log.Add(consts.FixedLn2.Multiply(twos)); + if (tens != 0) + log = log.Add(consts.FixedLn10.Multiply(tens)); + return FromFixed(log, consts, working).RoundToPrecision(context); + } + + /// + /// e to the to the precision of : + /// x = k ln 2 + r with |r| at most half of ln 2, r halved + /// ten times, the Taylor series there -- twenty-five terms at a hundred digits -- and + /// ten squarings and the power of two back, in fixed point. PeterO's + /// is three hundred microseconds at a hundred + /// digits; this is about ten. An argument beyond a thousand in magnitude, and an + /// infinity or NaN, are PeterO's answers. + /// https://github.com/asc-community/AngouriMath/issues/1338 + /// + public static EDecimal Exponential(this EDecimal x, EContext context) + { + if (!x.IsFinite || x.Abs().CompareTo(EDecimal.FromInt32(1000)) > 0) + return x.Exp(context); + if (x.IsZero) + return EDecimal.One; + // Twelve guard digits: the ten squarings each double the error, three digits, + // and the reduction by up to fifteen hundred times ln 2 costs four. + var working = WithGuardDigits(context, 12); + var consts = ConstantCache.Lookup(working); + var bits = consts.FixedBits; + var fixedX = ToFixed(x, bits); + var k = fixedX.Divide(consts.FixedLn2); + var r = fixedX.Subtract(consts.FixedLn2.Multiply(k)); + var halfLn2 = consts.FixedLn2.ShiftRight(1); + if (r.CompareTo(halfLn2) > 0) + { + k = k.Add(1); + r = r.Subtract(consts.FixedLn2); + } + else if (r.CompareTo(halfLn2.Negate()) < 0) + { + k = k.Subtract(1); + r = r.Add(consts.FixedLn2); + } + const int halvings = 10; + r = r.Sign < 0 ? r.Negate().ShiftRight(halvings).Negate() : r.ShiftRight(halvings); + var term = consts.FixedOne; + var sum = consts.FixedOne; + for (var n = 1; n < 100000; n++) + { + term = MultiplyFixed(term, r, bits).Divide(n); + if (term.IsZero) + break; + sum = sum.Add(term); + } + for (var i = 0; i < halvings; i++) + sum = MultiplyFixed(sum, sum, bits); + var power = k.ToInt32Checked(); + if (power >= 0) + return FromFixed(sum.ShiftLeft(power), consts, working).RoundToPrecision(context); + return FromFixed(sum, consts, working) + .Divide(EDecimal.FromEInteger(EInteger.One.ShiftLeft(-power)), working) + .RoundToPrecision(context); + } + /// Analogy of public static EDecimal Sinh(this EDecimal x, EContext context) { - var y = x.Exp(context); + var y = x.Exponential(context); var yy = EDecimal.One.Divide(y, context); return y.Subtract(yy, context).Divide(2, context); } @@ -653,7 +840,7 @@ public static EDecimal Sinh(this EDecimal x, EContext context) /// Analogy of public static EDecimal Cosh(this EDecimal x, EContext context) { - var y = x.Exp(context); + var y = x.Exponential(context); var yy = EDecimal.One.Divide(y, context); return y.Add(yy, context).Divide(2, context); } @@ -665,7 +852,7 @@ public static EDecimal Tanh(this EDecimal x, EContext context) return EDecimal.NaN; if (x.IsInfinity()) return x.Sign; - var y = x.Exp(context); + var y = x.Exponential(context); var yy = EDecimal.One.Divide(y, context); return y.Subtract(yy, context).Divide(y.Add(yy, context), context); } diff --git a/Sources/Tests/DotnetBenchmark/performance-baseline.json b/Sources/Tests/DotnetBenchmark/performance-baseline.json index 738d91f6a..f363dda35 100644 --- a/Sources/Tests/DotnetBenchmark/performance-baseline.json +++ b/Sources/Tests/DotnetBenchmark/performance-baseline.json @@ -1,94 +1,94 @@ { "comment": "Allocated bytes and mean nanoseconds per operation for the popular use cases of https://github.com/asc-community/AngouriMath/issues/746. Checked by PerformanceGate.cs; see Sources/AngouriMath/Docs/WhatsNew/version_performance_control.md for when updating it is legitimate.", "benchmark": "CommonFunctionsInterVersion", - "commit": "22395775e8fffb57d30b9ca8b0c338907895471a", - "measuredOn": "2026-09-14", + "commit": "bfb60e9cdcdc855c13b5ac58001db0ec98e3f3b2", + "measuredOn": "2026-09-15", "runtime": ".NET 10.0.10", "machine": "Ubuntu 26.04 LTS, X64, 8 logical cores", "cases": { "CompileEasy": { - "allocatedBytes": 11028, - "meanNanoseconds": 184281.8096842448 + "allocatedBytes": 11023, + "meanNanoseconds": 184521.07338867188 }, "CompileHard": { - "allocatedBytes": 20434, - "meanNanoseconds": 311105.99281529017 + "allocatedBytes": 20737, + "meanNanoseconds": 309526.83076171874 }, "Derivate": { - "allocatedBytes": 53111, - "meanNanoseconds": 14165.4219959804 + "allocatedBytes": 53099, + "meanNanoseconds": 14289.474793570382 }, "EvalEasy": { "allocatedBytes": 0, - "meanNanoseconds": 1.7091758726164699 + "meanNanoseconds": 1.8900332285889558 }, "EvalPolynomialFresh": { "allocatedBytes": 12328, - "meanNanoseconds": 4087.811504618327 + "meanNanoseconds": 4151.610092798869 }, "EvalPolynomialFresh15Digits": { "allocatedBytes": 11656, - "meanNanoseconds": 4086.574955531529 + "meanNanoseconds": 4141.208742959158 }, "EvalTranscendentalFresh": { - "allocatedBytes": 2246889, - "meanNanoseconds": 1328053.0219029018 + "allocatedBytes": 423032, + "meanNanoseconds": 245377.50509207588 }, "EvalTrig": { "allocatedBytes": 1192089, - "meanNanoseconds": 556000.5556640625 + "meanNanoseconds": 551794.5797293527 }, "EvalTrigPrecise": { "allocatedBytes": 7935226, - "meanNanoseconds": 13768291.227163462 + "meanNanoseconds": 13781485.40625 }, "ParseEasy": { "allocatedBytes": 18125, - "meanNanoseconds": 5926.468636067709 + "meanNanoseconds": 5981.867161090558 }, "ParseHard": { "allocatedBytes": 3555057, - "meanNanoseconds": 1378406.7677176339 + "meanNanoseconds": 1367037.3726283482 }, "RunEasy": { "allocatedBytes": 0, - "meanNanoseconds": 20.059395309289297 + "meanNanoseconds": 20.056596242464504 }, "RunHard": { "allocatedBytes": 0, - "meanNanoseconds": 291.68687656947543 + "meanNanoseconds": 295.51566811970304 }, "RunMedium": { "allocatedBytes": 0, - "meanNanoseconds": 163.2946068899972 + "meanNanoseconds": 163.17946696281433 }, "SimplifyEasy": { "allocatedBytes": 78578, - "meanNanoseconds": 45866.643728402945 + "meanNanoseconds": 45092.335445149736 }, "SimplifyHard": { - "allocatedBytes": 302423792, - "meanNanoseconds": 171960155.3846154 + "allocatedBytes": 189726120, + "meanNanoseconds": 109370614.94736843 }, "SolveEasy": { - "allocatedBytes": 6284124, - "meanNanoseconds": 3545737.0346354167 + "allocatedBytes": 4425194, + "meanNanoseconds": 2884867.8782552085 }, "SolveEasyMedium": { - "allocatedBytes": 64974, - "meanNanoseconds": 18541.630039760046 + "allocatedBytes": 65184, + "meanNanoseconds": 18597.25531768799 }, "SolveHard": { - "allocatedBytes": 11714880, - "meanNanoseconds": 70353131.95454545 + "allocatedBytes": 11706656, + "meanNanoseconds": 64677752.72 }, "SolveMedium": { - "allocatedBytes": 444046, - "meanNanoseconds": 266698.74176897324 + "allocatedBytes": 446445, + "meanNanoseconds": 266309.50540597097 }, "SolveMediumHard": { - "allocatedBytes": 1365686, - "meanNanoseconds": 7870781.697916667 + "allocatedBytes": 1367798, + "meanNanoseconds": 8561109.407291668 } }, "ungated": { diff --git a/Sources/Tests/UnitTests/Core/HighPrecisionFunctionsTest.cs b/Sources/Tests/UnitTests/Core/HighPrecisionFunctionsTest.cs index caffe748b..04caf98f3 100644 --- a/Sources/Tests/UnitTests/Core/HighPrecisionFunctionsTest.cs +++ b/Sources/Tests/UnitTests/Core/HighPrecisionFunctionsTest.cs @@ -84,6 +84,40 @@ public void WhereDigitsWereLost(string expression, string reference) public void WholeAndHalfPowersOfADecimal(string expression, string reference) => AgreesTo98Digits(expression.ToEntity(), reference); + /// + /// The logarithm, the exponential, a real power and the hyperbolic functions, which are + /// series in fixed point now rather than PeterO's: near one and far from it, a power of + /// ten and a power of e -- e^700 came back to ninety-seven digits as the + /// constant's hundred raised to the seven hundredth. Not below ten to the minus fifty: + /// a value within half the working digits of an integer is that integer by the + /// downcasting, so e^(-123.456), which is ten to the minus fifty-four, is 0. + /// + [Theory] + [InlineData("ln(0.37)", "-0.994252273343866923667887238337281251302125390089909788421743120742155472364343895072972876507027935902127")] + [InlineData("ln(1.1369)", "0.128305260152923841561498035928794898057419676174790423366932255144985972701456882347699236277855162624175")] + [InlineData("ln(0.9999)", "-0.000100005000333358335333500014286964396835397734571075514089865762716344756611402617022495986337466528754179")] + [InlineData("ln(1.00001)", "0.00000999995000033333083335333316666809522559534920534921544003210755133040854707452208102937253322237289734812")] + [InlineData("ln(750)", "6.62007320653035612461475535805926519129979475498855787159331801755342487831127697636988071608968960872987")] + [InlineData("ln(0.000001234)", "-13.6052496324810780327471192920909104756199783466598383294738650174947053425106213063054082875112674176785")] + [InlineData("ln(123456789.123456789)", "18.6314017671680180326939333482965375427970151745537353083517566119017412766551613015767513407252233304900")] + [InlineData("log(10, 0.37)", "-0.431798275933005003191549310460870552017027309833687453382320089206414574578527530257797656739549285246099")] + [InlineData("log(3, 81.5)", "4.00560148984118758502141372215507044442357868461710473567971265348027971343638223614214033083963425256794")] + [InlineData("e^0.37", "1.44773461466332446158475233551922961456683194184845491456061206892295956926655134553108541144999797619663")] + [InlineData("e^(-0.5)", "0.606530659712633423603799534991180453441918135487186955682892158735056519413748423998647611507989456026424")] + [InlineData("e^25.5", "118716009132.169650965201023040233373526449091282754098342267442638097384414108824078431787894306908281050")] + [InlineData("e^700", "1.01423205473500450945532959523126761520467957224307334878053628124935170250752368304548160316182971369539e+304")] + [InlineData("e^(-100.5)", "2.25634013591703631320281828241534366417520785936464023567989659964288626428305623878961101522065173672292e-44")] + [InlineData("e^0.00001234", "1.00001234007613811318111682981596300497517047048619342510199268785525707213639963880173472004809302592953")] + [InlineData("2^0.37", "1.29235283063749224450556503197070707880867324808527754395903593457909430919166483018351136908382587001355")] + [InlineData("10^(-0.37)", "0.426579518801592658004523896600975319499453929539193328227549626145447937991082374795110805295935878724416")] + [InlineData("1.7^2.9", "4.65909828378606751689142946026153054833060144518731308290895183497457393121344828497028520145573250549740")] + [InlineData("0.37^0.37", "0.692204850015734267986295538874892165757969134156257794701426036412631395756556597414370238390864299729380")] + [InlineData("sinh(0.37)", "0.378500142012984901014906185639889080855936538762531408933184424856263739634050081675804991499874453185980")] + [InlineData("cosh(-2.9)", "9.11458429474973408585310155248112110830188451116796598581181237851004903258039556059848985280824921641633")] + [InlineData("tanh(0.37)", "0.353991712477045994713511975971816782189871283978107285505039728240177889286936341011158953779536534740937")] + public void LogarithmsExponentialsAndPowers(string expression, string reference) + => AgreesTo98Digits(expression.ToEntity(), reference); + /// /// The downcasting to a rational is decided cheaply first and exactly after: what was /// a rational is still one, and what is not is not. diff --git a/Sources/Tests/UnitTests/Core/NumericDigits.cs b/Sources/Tests/UnitTests/Core/NumericDigits.cs index c915c92c9..849f34e28 100644 --- a/Sources/Tests/UnitTests/Core/NumericDigits.cs +++ b/Sources/Tests/UnitTests/Core/NumericDigits.cs @@ -48,22 +48,22 @@ public sealed class NumericDigits [Fact] public void CotanHalfPi() => Test(Cotan(0.5m * pi), "0"); [Fact] public void CotanI() => Test(Cotan(i), "-1.313035285499331303636161246930847832912013941240452655543152967567084270461874382674679241480856303i"); [Fact] public void Cotan3P2I() => Test(Cotan(3+2*i), "-0.0106047834703371017503168962077792923972675909391407246475487930531248300454417108979266267185038253 - 1.035746637764995396112758656897908320248306959923214042664373819742665974520371256149665849444234282i"); - [Fact] public void ArcsinM2() => Test(Arcsin(-2), "-1.570796326794896619231321691639751442098584699687552910487472296153908203143104499314017412671058534 - 1.316957896924816708625046347307968444026981971467516479768472256920460185416443976074219013450101788i"); + [Fact] public void ArcsinM2() => Test(Arcsin(-2), "-1.570796326794896619231321691639751442098584699687552910487472296153908203143104499314017412671058534 - 1.316957896924816708625046347307968444026981971467516479768472256920460185416443976074219013450101784i"); [Fact] public void Arcsin1() => Test(Arcsin(1), "1.570796326794896619231321691639751442098584699687552910487472296153908203143104499314017412671058534"); - [Fact] public void Arcsin3() => Test(Arcsin(3), "1.570796326794896619231321691639751442098584699687552910487472296153908203143104499314017412671058534 - 1.762747174039086050465218649959584618056320656523270821506591217306754368444052175667413783820512091i"); - [Fact] public void Arcsin4() => Test(Arcsin(4), "1.570796326794896619231321691639751442098584699687552910487472296153908203143104499314017412671058534 - 2.063437068895560546727281172620131871456591449883392499836032692765902842847409911780353006403580318i"); - [Fact] public void Arcsin3M2I() => Test(Arcsin(3-2*i), "0.9646585044076027920454110594995323555197773725073316527132580297155508786089335572049608301897631772 - 1.968637925793096291788665095245498189520731012682010573842811017352748255492485345887875752070076217i"); - [Fact] public void Arcsin3P2I() => Test(Arcsin(3+2*i), "0.9646585044076027920454110594995323555197773725073316527132580297155508786089335572049608301897631772 + 1.968637925793096291788665095245498189520731012682010573842811017352748255492485345887875752070076217i"); - [Fact] public void ArcsinI() => Test(Arcsin(i), "0.8813735870195430252326093249797923090281603282616354107532956086533771842220260878337068919102560457i"); - [Fact] public void ArccosM2() => Test(Arccos(-2), "3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068 + 1.316957896924816708625046347307968444026981971467516479768472256920460185416443976074219013450101788i"); + [Fact] public void Arcsin3() => Test(Arcsin(3), "1.570796326794896619231321691639751442098584699687552910487472296153908203143104499314017412671058534 - 1.762747174039086050465218649959584618056320656523270821506591217306754368444052175667413783820512086i"); + [Fact] public void Arcsin4() => Test(Arcsin(4), "1.570796326794896619231321691639751442098584699687552910487472296153908203143104499314017412671058534 - 2.063437068895560546727281172620131871456591449883392499836032692765902842847409911780353006403580326i"); + [Fact] public void Arcsin3M2I() => Test(Arcsin(3-2*i), "0.9646585044076027920454110594995323555197773725073316527132580297155508786089335572049608301897631772 - 1.968637925793096291788665095245498189520731012682010573842811017352748255492485345887875752070076231i"); + [Fact] public void Arcsin3P2I() => Test(Arcsin(3+2*i), "0.9646585044076027920454110594995323555197773725073316527132580297155508786089335572049608301897631772 + 1.968637925793096291788665095245498189520731012682010573842811017352748255492485345887875752070076231i"); + [Fact] public void ArcsinI() => Test(Arcsin(i), "0.8813735870195430252326093249797923090281603282616354107532956086533771842220260878337068919102560430i"); + [Fact] public void ArccosM2() => Test(Arccos(-2), "3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068 + 1.316957896924816708625046347307968444026981971467516479768472256920460185416443976074219013450101784i"); [Fact] public void Arccos1() => Test(Arccos(1), "0"); - [Fact] public void Arccos3() => Test(Arccos(3), "1.762747174039086050465218649959584618056320656523270821506591217306754368444052175667413783820512091i"); - [Fact] public void Arccos4() => Test(Arccos(4), "2.063437068895560546727281172620131871456591449883392499836032692765902842847409911780353006403580318i"); - [Fact] public void ArccosI() => Test(Arccos(i), "1.570796326794896619231321691639751442098584699687552910487472296153908203143104499314017412671058534 - 0.8813735870195430252326093249797923090281603282616354107532956086533771842220260878337068919102560457i"); - [Fact] public void Arccos3P2I() => Test(Arccos(3+2*i), "0.6061378223872938271859106321402190865788073271802212577742142664383573245341709421090565824812953568 - 1.968637925793096291788665095245498189520731012682010573842811017352748255492485345887875752070076217i"); + [Fact] public void Arccos3() => Test(Arccos(3), "1.762747174039086050465218649959584618056320656523270821506591217306754368444052175667413783820512086i"); + [Fact] public void Arccos4() => Test(Arccos(4), "2.063437068895560546727281172620131871456591449883392499836032692765902842847409911780353006403580326i"); + [Fact] public void ArccosI() => Test(Arccos(i), "1.570796326794896619231321691639751442098584699687552910487472296153908203143104499314017412671058534 - 0.8813735870195430252326093249797923090281603282616354107532956086533771842220260878337068919102560430i"); + [Fact] public void Arccos3P2I() => Test(Arccos(3+2*i), "0.6061378223872938271859106321402190865788073271802212577742142664383573245341709421090565824812953568 - 1.968637925793096291788665095245498189520731012682010573842811017352748255492485345887875752070076231i"); [Fact] public void Arctan1() => Test(Arctan(1), "0.7853981633974483096156608458198757210492923498437764552437361480769541015715522496570087063355292670"); [Fact] public void ArctanI() => Test(Arctan(i), "+ooi"); - [Fact] public void Arctan3P2I() => Test(Arctan(3+2*i), "1.338972522294493561124193575909144241084316172544492778582005751793809271060233646663717270678614588 + 0.1469466662255297520474327851547159424423449403442452953891851939502023996823900422792744078835711450i"); + [Fact] public void Arctan3P2I() => Test(Arctan(3+2*i), "1.338972522294493561124193575909144241084316172544492778582005751793809271060233646663717270678614588 + 0.1469466662255297520474327851547159424423449403442452953891851939502023996823900422792744078835711420i"); [Fact] public void Arccotan1() => Test(Arccotan(1), "0.7853981633974483096156608458198757210492923498437764552437361480769541015715522496570087063355292670"); [Fact] public void ArccotanI() => Test(Arccotan(i), "-ooi"); [Fact] public void Arccotan3P2I() => Test(Arccotan(3+2*i), "0.2318238045004030581071281157306072010142685271430601319054665443600989320828708526503001419924439463 - 0.1469466662255297520474327851547159424423449403442452953891851939502023996823900422792744078835711419i"); @@ -76,19 +76,19 @@ public sealed class NumericDigits [Fact] public void LogM2_M2() => Test(Log(-2, -2), "1"); [Fact] public void LogM2_0D5() => Test(Log(-2, 0.5m), "-0.04642032354540812805739794967614531473517486433229878976878405934128995144281958985056393738882100657 + 0.2103936242079302074624500066187461838994455574074656160548588507113367960711277189504016452320136900i"); [Fact] public void LogM2_0D4() => Test(Log(-2, 0.4m), "-0.06136432986843633646517722130553639639897513966073574085079108126276417408725123183789917075670358045 + 0.2781252428256368191899557377147340479778729417292667461789200592366343982454240662260594052903331773i"); - [Fact] public void LogM2_3P2I() => Test(Log(-2, 3+2*i), "0.2643664916612894989375747779599512529839212984780690192522150852389360266747092012500969179019942126 - 0.3498957094724448147885913933302844456684320929581102596014846703961864239635800903604529923905776475i"); + [Fact] public void LogM2_3P2I() => Test(Log(-2, 3+2*i), "0.2643664916612894989375747779599512529839212984780690192522150852389360266747092012500969179019942119 - 0.3498957094724448147885913933302844456684320929581102596014846703961864239635800903604529923905776442i"); [Fact] public void LnE() => Test(Ln(e), "1"); - [Fact] public void LnPi() => Test(Ln(pi), "1.144729885849400174143427351353058711647294812915311571513623071472137769884826079783623270275489713"); + [Fact] public void LnPi() => Test(Ln(pi), "1.144729885849400174143427351353058711647294812915311571513623071472137769884826079783623270275489708"); [Fact] public void Pow2_3() => Test(Pow(2, 3), "8"); - [Fact] public void Pow2_I() => Test(Pow(2, i), "0.7692389013639721265783299936612707014408959949119638531698715074290813468073407890597897424260168074 + 0.6389612763136348011500329114647017842572305378305797294955869566463245224485447499034456609322567384i"); + [Fact] public void Pow2_I() => Test(Pow(2, i), "0.7692389013639721265783299936612707014408959949119638531698715074290813468073407890597897424260168073 + 0.6389612763136348011500329114647017842572305378305797294955869566463245224485447499034456609322567385i"); [Fact] public void PowI_2() => Test(Pow(i, 2), "-1"); [Fact] public void PowI_I() => Test(Pow(i, i), "0.2078795763507619085469556198349787700338778416317696080751358830554198772854821397886002778654260353"); [Fact] public void Pow0D5M0D8I_0D5() => Test(Pow(0.5m-0.8m*i, 0.5m), "0.8495287261787150364907695626370349126367049237188233590072351305588742785045041120988684216395800854 - 0.4708492928770629369416080841473292873459093837847060526082076263776372676635818831502937646218684427i"); [Fact] public void Sqrt0D5M0D8I() => Test(Sqrt(0.5m-0.8m*i), "0.8495287261787150364907695626370349126367049237188233590072351305588742785045041120988684216395800853 - 0.4708492928770629369416080841473292873459093837847060526082076263776372676635818831502937646218684426i"); - [Fact] public void Pow3P2I_3M2I() => Test(Pow(3+2*i,3-2*i), "105.7489754858859776057429151808913270583911891854438642007164406579203854865458413612560669739008029 - 109.0885464095191985068230453514141449073149001480712525014939233398460315014292268679635166428957493i"); - [Fact] public void Pow3P2I_3P2I() => Test(Pow(3+2*i,3+2*i), "-5.409738793917671859317094187110414923431153201422932321801121451217085143160320149482544245235856010 - 13.41044237041274830872806422189982707487514435448700312177401211574208044465073698802176970443887716i"); - [Fact] public void Pow3M2I_3M2I() => Test(Pow(3-2*i,3-2*i), "-5.409738793917671859317094187110414923431153201422932321801121451217085143160320149482544245235856010 + 13.41044237041274830872806422189982707487514435448700312177401211574208044465073698802176970443887716i"); - [Fact] public void Pow3M2I_3P2I() => Test(Pow(3-2*i,3+2*i), "105.7489754858859776057429151808913270583911891854438642007164406579203854865458413612560669739008029 + 109.0885464095191985068230453514141449073149001480712525014939233398460315014292268679635166428957493i"); + [Fact] public void Pow3P2I_3M2I() => Test(Pow(3+2*i,3-2*i), "105.7489754858859776057429151808913270583911891854438642007164406579203854865458413612560669739008053 - 109.0885464095191985068230453514141449073149001480712525014939233398460315014292268679635166428957470i"); + [Fact] public void Pow3P2I_3P2I() => Test(Pow(3+2*i,3+2*i), "-5.409738793917671859317094187110414923431153201422932321801121451217085143160320149482544245235856305 - 13.41044237041274830872806422189982707487514435448700312177401211574208044465073698802176970443887704i"); + [Fact] public void Pow3M2I_3M2I() => Test(Pow(3-2*i,3-2*i), "-5.409738793917671859317094187110414923431153201422932321801121451217085143160320149482544245235856305 + 13.41044237041274830872806422189982707487514435448700312177401211574208044465073698802176970443887704i"); + [Fact] public void Pow3M2I_3P2I() => Test(Pow(3-2*i,3+2*i), "105.7489754858859776057429151808913270583911891854438642007164406579203854865458413612560669739008053 + 109.0885464095191985068230453514141449073149001480712525014939233398460315014292268679635166428957470i"); [Fact] public void FactorialM1() => Test(Factorial(-1), "NaN"); [Fact] public void Factorial1() => Test(Factorial(1), "1"); [Fact] public void FactorialI() => Test(Factorial(i), "0.4980156681183560427136911174621980919529629675876500928926429549984583004359819345078945042826705814 - 0.1549498283018106851249551304838866051958796520793249302658802767988608014911385390129513664794630707i");