From 97169ebfa2bd2405eec5ae087537fdd3fc089efc Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 16 Oct 2026 14:27:23 +0000 Subject: [PATCH 1/3] Integrate a power of a sine or hyperbolic sine over a binomial in a power of it over its roots s^k/(a + b s^n) for a whole k of either sign is divided down, or up, to 0 <= k < n over the binomial, the powers of s beside going to the integrator, and is then (1/(n a)) sum r^(k + 1)/(r - s) over all n roots. Asked first where a power stands above or below, since the rules before it spent the budget: Rubi's 6.1.7 sinh(c + d x)^4/(a + b sinh(c + d x)^3), csch(c + d x)^3/(a + b sinh(c + d x)^3) and the like timed out. DerivativeDiffersWhereReal now decides a difference, or no value, found at the default precision again in forty and eighty digits: the derivative of the right answer to sinh(c + d x)^5/(a + b sinh(c + d x)^3) was NaN there at a sampled point, which counted as a difference. Part of https://github.com/asc-community/AngouriMath/issues/718 Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../Algebra/Polynomials/PartialFractions.cs | 29 +++++-- .../Integration/IndefiniteIntegralSolver.cs | 79 +++++++++++++++---- .../Integration/Integration.Definition.cs | 3 + ...lOfABinomialInAPowerOfASineIntegralTest.cs | 15 ++++ 4 files changed, 106 insertions(+), 20 deletions(-) diff --git a/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs b/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs index 8091a8004..c0877cbc6 100644 --- a/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs +++ b/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs @@ -2147,18 +2147,37 @@ internal static bool DiffersWhereReal(Entity left, Entity right, Variable x) using var _ = MathS.Settings.DowncastingEnabled.Set(false); (left, right) = Pinned(left, right, x); var tolerance = EDecimal.FromString("1e-9"); + var inDecimals = new System.Lazy[] { new(() => PreciseEvaluation.In(40)), new(() => PreciseEvaluation.In(80)) }; foreach (var at in DefaultSamplePoints) { var point = Real.Create(EDecimal.FromString(at)); - if (right.Substitute(x, point).Evaled is not Complex r || r.IsNaN) + var rightAt = right.Substitute(x, point); + if (rightAt.Evaled is not Complex r || r.IsNaN) continue; var scale = EDecimal.Max(EDecimal.One, r.Abs().EDecimal); if (r.ImaginaryPart.EDecimal.Abs().CompareTo(scale.Multiply(tolerance)) > 0) continue; - if (left.Substitute(x, point).Evaled is not Complex l || l.IsNaN) - return true; - if ((l - r).Abs().EDecimal.CompareTo(scale.Multiply(tolerance)) > 0) - return true; + var leftAt = left.Substitute(x, point); + if (leftAt.Evaled is Complex l && !l.IsNaN && (l - r).Abs().EDecimal.CompareTo(scale.Multiply(tolerance)) <= 0) + continue; + // A difference, or no value, at the default precision is decided again in forty and + // eighty digits: `tanh` written as `(e^(2 y) - 1)/(e^(2 y) + 1)` is inf/inf there at + // y near 30, and a right answer to `sinh(c + d x)^5/(a + b sinh(c + d x)^3)` was + // refused for it. Still no value where the integrand has one is a difference; a + // precision that cannot tell says nothing. + foreach (var precise in inDecimals) + { + if (precise.Value.ValueOf(rightAt) is not { } rightPrecise) + continue; + if (precise.Value.ValueOf(leftAt) is not { } leftPrecise) + return true; + if (PreciseEvaluation.Agree(leftPrecise, rightPrecise, 1e-9) is { } decided) + { + if (!decided) + return true; + break; + } + } } return false; } diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 1c6458f74..a9442a720 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -25956,9 +25956,19 @@ static bool IsASumOfPowersOfX(Entity sum, Entity.Variable x) /// 1/(a + b sin(x)^5), 1/(a - b cos(x)^6) and the like were declined. /// https://github.com/asc-community/AngouriMath/issues/718 /// + /// + /// asked first where a power of + /// the function stands above or below the binomial: the rules before it in the dispatcher spent + /// the budget on `sinh(y)^4/(a + b sinh(y)^3)`, and the bare reciprocal they answer is left to them. + /// + internal static Entity? SolveAPowerOverABinomialInAPowerOfASineOverItsRoots(Entity expr, Entity.Variable x) + => Integration.AnsweringTheQuestionAsked && TryReadAsQuotient(expr, out var above, out _) && above.ContainsNode(x) + ? SolveAReciprocalOfABinomialInAPowerOfASineOverItsRoots(expr, x) + : null; + internal static Entity? SolveAReciprocalOfABinomialInAPowerOfASineOverItsRoots(Entity expr, Entity.Variable x) { - if (!TryReadAsQuotient(expr, out var numerator, out var denominator) || numerator.ContainsNode(x)) + if (!TryReadAsQuotient(expr, out var numerator, out var denominator)) return null; // The hyperbolic sine and cosine are parsed as `(e^y -+ e^(-y))/2`, and read back from it. static (Entity Argument, bool Sine, bool Hyperbolic)? Read(Entity node) => node switch @@ -25981,31 +25991,69 @@ when two.Equals(EInteger.FromInt32(2)) && e1 == MathS.e && e2 == MathS.e if (!TreeAnalyzer.TryGetPolyLinear(argument, x, out var slope, out _) || Functions.PartialFractions.IsZeroAsAValue(slope)) return null; var s = Variable.CreateUnique(expr, "s_bin"); - var inS = denominator.Substitute(function, s); - if (inS.ContainsNode(x) || !TreeAnalyzer.TryGetPolynomial(inS, s, out var read) || read.Count != 2 - || !read.ContainsKey(EInteger.Zero) || read.Keys.Max() is not { } top || !top.CanFitInInt32()) + // In s the integrand is C s^k/(a + b s^n) for a whole k of either sign: a power of the + // function above the bar, or of its reciprocal, `csch(y)^3/(a + b sinh(y)^3)`. + var (aboveInS, belowInS) = Functions.SingleQuotient.Of(expr.Substitute(function, s)); + if (aboveInS.ContainsNode(x) || belowInS.ContainsNode(x) + || !TreeAnalyzer.TryGetPolynomial(aboveInS, s, out var readAbove) || readAbove.Count != 1 + || !TreeAnalyzer.TryGetPolynomial(belowInS, s, out var readBelow) || readBelow.Count != 2) + return null; + var lowest = readBelow.Keys.Min()!; + var top = readBelow.Keys.Max()!; + var aboveKey = readAbove.Keys.Single(); + if (!top.CanFitInInt32() || !lowest.CanFitInInt32() || !aboveKey.CanFitInInt32()) + return null; + var n = (top - lowest).ToInt32Unchecked(); + var k = (aboveKey - lowest).ToInt32Unchecked(); + Entity front0 = readAbove[aboveKey]; + var a = readBelow[lowest]; + var b = readBelow[top]; + // The binomial alone, read off the denominator as it is written, as it always was. + if (!numerator.ContainsNode(x) && TreeAnalyzer.TryGetPolynomial(denominator.Substitute(function, s), s, out var asWritten) + && asWritten.Count == 2 && asWritten.ContainsKey(EInteger.Zero)) + (front0, a, b, k) = (numerator, asWritten[EInteger.Zero], asWritten[asWritten.Keys.Max()!], 0); + else if (numerator.ContainsNode(x) == false && k != 0) return null; - var n = top.ToInt32Unchecked(); // Past what the rules before it answer: they write `1/(a + b cosh(x)^3)` through `e^x` in - // 655 characters, where this writes 1491. - if (n < 5 || n > 12) + // 655 characters, where this writes 1491. With a power beside it they timed out at n = 3: + // `sinh(y)^4/(a + b sinh(y)^3)`. + if (n > 12 || n < 3 || n < 5 && k == 0 || System.Math.Abs(k) > 24) return null; - var a = read[EInteger.Zero]; - var b = read[top]; if (!(a + b).Vars.Any() || Functions.PartialFractions.IsZeroAsAValue(a) || Functions.PartialFractions.IsZeroAsAValue(b)) return null; - var q = n % 2 == 0 ? 2 : 1; + // s^k/(a + b s^n) is s^(k - n)/b - (a/b) s^(k - n)/(a + b s^n) down to k < n, and + // s^k/a - (b/a) s^(k + n)/(a + b s^n) up to k >= 0; the powers of s go to the integrator. + Entity powers = Number.Integer.Zero; + while (k >= n) + { + powers += front0 / b * MathS.Pow(function, Number.Integer.Create(k - n)); + front0 = -front0 * a / b; + k -= n; + } + while (k < 0) + { + powers += front0 / a * MathS.Pow(function, Number.Integer.Create(k)); + front0 = -front0 * b / a; + k += n; + } + Entity? powersIntegrated = Number.Integer.Zero; + if (powers != Number.Integer.Zero + && (powersIntegrated = Integration.ComputeAsAQuestionOfItsOwn(powers.InnerSimplified, x, integrateByParts: true)) is null) + return null; + numerator = front0; + // Above the bar s^k with 0 < k < n: over all n roots, (1/(n a)) sum r^(k + 1)/(r - s). + var q = n % 2 == 0 && k == 0 ? 2 : 1; var m = n / q; var half = q == 2 ? argument : argument / 2; var u = hyperbolic ? MathS.Hyperbolic.Tanh(half) : MathS.Tan(half); var radius = MathS.Pow(-a / b, Number.Rational.Create(1, m)); Entity sum = Number.Integer.Zero; - for (var k = 0; k < m; k++) + for (var j = 0; j < m; j++) { // The root of unity as a cosine and a sine: `(-1)^(2/3)` is 1 here, the real cube // root of 1, where it means the principal `e^(2 pi i/3)` to Rubi. - var turn = 2 * MathS.pi * k / m; - var r = (k == 0 ? radius : radius * (MathS.Cos(turn) + MathS.i * MathS.Sin(turn))).InnerSimplified; + var turn = 2 * MathS.pi * j / m; + var r = (j == 0 ? radius : radius * (MathS.Cos(turn) + MathS.i * MathS.Sin(turn))).InnerSimplified; // K^2, the factor in front, and the argument of the arctangent over K. var (square, front, shift) = (q, sine, hyperbolic) switch { @@ -26023,9 +26071,10 @@ when two.Equals(EInteger.FromInt32(2)) && e1 == MathS.e && e2 == MathS.e // An even power under tan(y), where K multiplies u; an odd one under tan(y/2), where it divides. // The hyperbolic ones are the same with the inverse hyperbolic tangent. Entity Inverse(Entity z) => hyperbolic ? MathS.Hyperbolic.Artanh(z) : MathS.Arctan(z); - sum += q == 2 ? front * Inverse(root * u) / root : front * Inverse((u + shift) / root) / root; + // r^0 is not folded into 1 when r is a root written out, so it is left off. + sum += (k == 0 ? Number.Integer.One : MathS.Pow(r, Number.Integer.Create(k))) * (q == 2 ? front * Inverse(root * u) / root : front * Inverse((u + shift) / root) / root); } - var answer = (numerator / (m * a * slope) * sum).InnerSimplified; + var answer = (powersIntegrated + numerator / (m * a * slope) * sum).InnerSimplified; // Refused where it differs, not where it cannot be told: the intervals do not settle // the arctangents of complex arguments at every sampled point, and `1/(a + b sin(x)^5)`, // right at every point checked exactly, was declined for too few comparisons. diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 8b11c7b49..7df763d43 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -664,6 +664,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // A sum free of the variable that comes to a number multiplied out, as that number, // before anything reads a coefficient or an exponent in it: `2 + 2 a - 2 (1 + a)`. if ((answer = IndefiniteIntegralSolver.SolveByFoldingConstantSumsToTheirNumbers(expr, x, integrateByParts)) is { }) return answer; + // A power of a sine, cosine or hyperbolic one over a binomial in a power of it, over the + // binomial's roots, before the rules that spent the budget on `sinh(y)^4/(a + b sinh(y)^3)`. + if ((answer = IndefiniteIntegralSolver.SolveAPowerOverABinomialInAPowerOfASineOverItsRoots(expr, x)) is { }) return answer; // A whole power of a linear in the inverse sine, cosine or hyperbolic cosine of // `c + d x^2` with `c^2 = 1`, by its recurrence, which closes: by parts it went round. if ((answer = IndefiniteIntegralSolver.SolveAPowerOfAnInverseFunctionOfAQuadraticByTheRecurrence(expr, x)) is { }) return answer; diff --git a/Sources/Tests/UnitTests/Calculus/AReciprocalOfABinomialInAPowerOfASineIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AReciprocalOfABinomialInAPowerOfASineIntegralTest.cs index 40e43811b..8d93f2e93 100644 --- a/Sources/Tests/UnitTests/Calculus/AReciprocalOfABinomialInAPowerOfASineIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/AReciprocalOfABinomialInAPowerOfASineIntegralTest.cs @@ -61,5 +61,20 @@ public void IsAnsweredAndDifferentiatesBack(string integrand) DifferentiatesBack(integrand, ("a", 0.6), ("b", 1.25)); DifferentiatesBack(integrand, ("a", -0.75), ("b", 0.4)); } + + // A power of the function above or below the binomial: divided down, or up, to a power between + // 0 and n over the binomial, over all n roots, and the powers of the function beside. + // `sinh(y)^5` is answered only because a value the default precision cannot form is decided + // in decimals before it counts as a difference. + [Theory] + [InlineData("sinh(c+d*x)^5/(a+b*sinh(c+d*x)^3)")] + [InlineData("sinh(c+d*x)^2/(a+b*sinh(c+d*x)^3)")] + [InlineData("csch(c+d*x)^3/(a+b*sinh(c+d*x)^3)")] + [InlineData("sin(x)^2/(a+b*sin(x)^5)")] + public void APowerOverItIsAnsweredAndDifferentiatesBack(string integrand) + { + DifferentiatesBack(integrand, ("a", 0.6), ("b", 1.25), ("c", 0.3), ("d", 0.75)); + DifferentiatesBack(integrand, ("a", -0.75), ("b", 0.4), ("c", -1), ("d", 0.5)); + } } } From 6594fce039843cd9d53d8db06120dd389de1b0e8 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 16 Oct 2026 19:31:33 +0000 Subject: [PATCH 2/3] Leave a power over a binomial in the fourth power to the rules before it They answer sin(y)^3/(a - b sin(y)^4) in a tenth of a second, where asked first the split took half of one. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../Continuous/Integration/IndefiniteIntegralSolver.cs | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index a9442a720..bd921fc79 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -26017,7 +26017,9 @@ when two.Equals(EInteger.FromInt32(2)) && e1 == MathS.e && e2 == MathS.e // Past what the rules before it answer: they write `1/(a + b cosh(x)^3)` through `e^x` in // 655 characters, where this writes 1491. With a power beside it they timed out at n = 3: // `sinh(y)^4/(a + b sinh(y)^3)`. - if (n > 12 || n < 3 || n < 5 && k == 0 || System.Math.Abs(k) > 24) + // And at n = 4 with a power beside it too: they answer `sin(y)^3/(a - b sin(y)^4)` in a tenth + // of a second, where asked first this took half of one. + if (n > 12 || n < 3 || n < 5 && k == 0 || n == 4 || System.Math.Abs(k) > 24) return null; if (!(a + b).Vars.Any() || Functions.PartialFractions.IsZeroAsAValue(a) || Functions.PartialFractions.IsZeroAsAValue(b)) return null; From 15bc1fa18101e7019e91674dc883aa51f76b0892 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 16 Oct 2026 19:42:30 +0000 Subject: [PATCH 3/3] Record the power over a binomial in a power of a sine in the breaking changes Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 8dc7f4401..6f06b1ea1 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -396,6 +396,21 @@ in Rubi's 1.1.3.4 come out shorter as well, and those that took two seconds take | `"1/(x^8*(a + b*x^3)*(c + d*x^3))".ToEntity().Integrate("x")` | `integral(...)` | 1,217 characters | | `"1/(x^6*(a + b*x^4)*(c + d*x^4))".ToEntity().Integrate("x")` | `integral(...)` | 795 characters | +### A power of a sine or hyperbolic sine over a binomial in a power of it is split over the binomial's roots + +**Answers where there were none.** `sinh(c + d x)^4/(a + b sinh(c + d x)^3)` was declined, and so were eight more of +Rubi's 6.1.7 with a power of the hyperbolic sine or cosecant beside `a + b sinh(c + d x)^3`. `s^k/(a + b s^n)` is now +divided down, or up, to a power between 0 and `n` over the binomial, the powers of `s` beside going to the integrator, +and split over all `n` roots of the binomial as the reciprocal alone already was +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sinh(c+d*x)^4/(a+b*sinh(c+d*x)^3)".ToEntity().Integrate("x")` | `integral(...)` | 1,497 characters | +| `"csch(c+d*x)^3/(a+b*sinh(c+d*x)^3)".ToEntity().Integrate("x")` | `integral(...)` | 1,566 characters | +| `"sinh(c+d*x)^6/(a+b*sinh(c+d*x)^3)".ToEntity().Integrate("x")` | `integral(...)` | 1,517 characters | +| `"sin(x)^2/(a+b*sin(x)^5)".ToEntity().Integrate("x")` | `integral(...)` | 2,475 characters | + ### A rational function times a logarithm of nested powers of linears is integrated by parts **Answers where there were none.** `ln(e (f (a + b x)^p (c + d x)^q)^r)` was declined, and so were four more of Rubi's