diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index dd5a644b0..29796da1e 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -763,6 +763,33 @@ unevaluated are now answered ([#718](https://github.com/asc-community/AngouriMat | `"e^(4*atanh(a*x))*(k-a*k*x)^p".ToEntity().Integrate("x")` | `integral(...)` | 434 characters | | `"e^atanh(x)*x*(1+x)^2".ToEntity().Integrate("x")` | `integral(...)` | 1,110 characters | +### A rational function of a sine over powers of linears in it is split by the residues at their roots + +**Answers where there were none.** `cot(c + d x)^6/(a + b sin(c + d x))^3` was declined, and so were eight more of +Rubi's 4.1.1.3, 4.1.2.1, 4.1.2.2 and 4.3.1.3 ([#233](https://github.com/asc-community/AngouriMath/issues/233)): +- In `s = sin(x)` that integrand is `(1 - s^2)^3/(s^6 (a + b s)^3)`, over powers of linears in `s`: the sine, `1 ± s` + from an even power of the cosine, and `a + b s`. It is split by the residues at their roots, each block a Taylor + series at a root, and every piece is a constant over a power of one linear, which the rules for `csc(x)^n`, + `1/(1 ± sin(x))^n` and `1/(a + b sin(x))^n` answer. Under the half-angle tangent, where the unreleased master took + them, every linear but the sine is a quadratic, and the cube of a symbolic one ran out of time. Only where the + half-angle is that slow: a linear with a symbolic root to the third power or more, and the half-angle's denominator + past the eleventh degree. Not for the cosine, whose half-angle answer is the shorter. +- A repeated quadratic with Gaussian rational roots `r = ρ + σ i` and `r̄` is split by the residues at `r`, one + division of series per order, with those at `r̄` their conjugates. `sin(c + d x)^6/(a + b tan(c + d x))^2` is + `t^6/((1 + t^2)^4 (a + b t)^2)` in `t = tan(c + d x)`, and was declined. + +`csc(c + d x)^2 sec(c + d x)^2/(a + b sin(c + d x))^3` was declined by 2.5.0 too, and the unreleased master answered +it in eleven seconds and 104,032 characters, which are a third of a second and 7,926 now. + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"tan(c + d*x)^4/(a + b*sin(c + d*x))^3".ToEntity().Integrate("x")` | `integral(...)` | 9,096 characters | +| `"cot(c + d*x)^6/(a + b*sin(c + d*x))^3".ToEntity().Integrate("x")` | `integral(...)` | 9,316 characters | +| `"1/((a + b*sin(k + f*x))^3*(c + d*sin(k + f*x))^3)".ToEntity().Integrate("x")` | `integral(...)` | 15,245 characters | +| `"sin(c + d*x)^6/(a + b*tan(c + d*x))^2".ToEntity().Integrate("x")` | `integral(...)` | 1,818 characters | +| `"sin(c + d*x)^4/(a + b*tan(c + d*x))^4".ToEntity().Integrate("x")` | `integral(...)` | 2,586 characters | +| `"csc(c + d*x)^2*sec(c + d*x)^2/(a + b*sin(c + d*x))^3".ToEntity().Integrate("x")` | `integral(...)` | 7,926 characters, where the unreleased master gave 104,032 | + ### A sine's binomial cube below the bar is split over the cube's two factors **Answers where there were none.** `csc(c + d x)/(a + b sin(c + d x)^3)` was declined, and so were the second to fifth diff --git a/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs b/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs index 105e047ad..f9f92cb32 100644 --- a/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs +++ b/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs @@ -1036,6 +1036,24 @@ private static bool TrySplitOverSymbolicLinearFactors( if (!TreeAnalyzer.TryGetPolynomial(quadraticFactor, x, out var quadratic) || !quadratic.ContainsKey(EInteger.FromInt32(2))) return false; var leading = quadratic[EInteger.FromInt32(2)]; + // A repeated quadratic over the rationals whose roots are Gaussian rationals, the + // tangent's `1 + u^2`, by its residues at one root: Newton's iteration below takes + // the symbols of the other factors through every round, and Rubi's + // `sin(x)^6/(a + b tan(x))^2`, `u^6/((1 + u^2)^4 (a + b u)^2)` under the tangent, + // ran past thirty seconds in it. + // https://github.com/asc-community/AngouriMath/issues/233 + if (multiplicity > 1 && OverTheGaussianRoots(numerator, factors, written, quadratic, multiplicity, x) is { } overThePowers) + { + // Over `leading^i Q^(j-i)` as the digits below are, for the same answer. + for (var power = multiplicity; power >= 1; power--) + { + var below = multiplicity - power; + var scale = below > 0 && leading != Integer.One ? (below == 1 ? leading : MathS.Pow(leading, below)) : Integer.One; + if (OverAPowerOfTheQuadratic(overThePowers[power - 1].P, overThePowers[power - 1].Q, quadraticFactor, power, scale, x) is { } term) + sum += term; + } + continue; + } var s = InLowestTermsOverTheSymbols((quadratic.TryGetValue(EInteger.One, out var q1) ? q1 : Integer.Zero) / leading); var t = InLowestTermsOverTheSymbols((quadratic.TryGetValue(EInteger.Zero, out var q0) ? q0 : Integer.Zero) / leading); // The numerators over the powers of the quadratic are the digits of @@ -1072,26 +1090,12 @@ private static bool TrySplitOverSymbolicLinearFactors( var (p, q) = digits[i]; if (p.Nodes.Any(node => node == MathS.NaN) || q.Nodes.Any(node => node == MathS.NaN)) return false; - Entity over = q == Integer.Zero ? p : p == Integer.Zero ? q * x : p + q * x; - if (over == Integer.Zero) - continue; // The digits are of `N/others` in powers of the monic quadratic `Q/leading`, // and `1/(Q/leading)^(j-i)` is `leading^(j-i)/Q^(j-i)`, so over `Q` itself the // `i`th digit stands over `leading^i Q^(j-i)`. - // Written with the coefficients' common denominator in front: left inside, - // `(P/D + Q x/D)/(1 + x^2)` was combined into one quotient over `D (1 + x^2)`, - // a quadratic with symbols in every coefficient, and integrated as a - // piecewise on the sign of its discriminant. - var (pAbove, pBelow) = SingleQuotient.Of(p); - var (qAbove, qBelow) = SingleQuotient.Of(q); - var (above, below) = pBelow == qBelow - ? (pAbove + qAbove * x, pBelow) - : (pAbove * qBelow + qAbove * pBelow * x, pBelow * qBelow); - var remaining = multiplicity - i; - Entity under = remaining == 1 ? quadraticFactor : MathS.Pow(quadraticFactor, remaining); - if (i > 0 && leading != Integer.One) - below *= i == 1 ? leading : MathS.Pow(leading, i); - sum += below == Integer.One ? above / under : MathS.Pow(below, Integer.MinusOne) * (above / under); + var scale = i > 0 && leading != Integer.One ? (i == 1 ? leading : MathS.Pow(leading, i)) : Integer.One; + if (OverAPowerOfTheQuadratic(p, q, quadraticFactor, multiplicity - i, scale, x) is { } term) + sum += term; } } @@ -1101,6 +1105,274 @@ private static bool TrySplitOverSymbolicLinearFactors( return true; } + /// + /// (p + q x)/(scale Q^power) for the written quadratic Q, with the + /// coefficients' common denominator in front; null where the numerator is zero. + /// + private static Entity? OverAPowerOfTheQuadratic(Entity p, Entity q, Entity quadraticFactor, int power, Entity scale, Variable x) + { + Entity over = q == Integer.Zero ? p : p == Integer.Zero ? q * x : p + q * x; + if (over == Integer.Zero) + return null; + // Written with the coefficients' common denominator in front: left inside, + // `(P/D + Q x/D)/(1 + x^2)` was combined into one quotient over `D (1 + x^2)`, + // a quadratic with symbols in every coefficient, and integrated as a + // piecewise on the sign of its discriminant. + var (pAbove, pBelow) = SingleQuotient.Of(p); + var (qAbove, qBelow) = SingleQuotient.Of(q); + var (above, below) = pBelow == qBelow + ? (pAbove + qAbove * x, pBelow) + : (pAbove * qBelow + qAbove * pBelow * x, pBelow * qBelow); + Entity under = power == 1 ? quadraticFactor : MathS.Pow(quadraticFactor, power); + if (scale != Integer.One) + below *= scale; + return below == Integer.One ? above / under : MathS.Pow(below, Integer.MinusOne) * (above / under); + } + + /// + /// The numerators p + q x over l^(k - m) Q^m, l the leading + /// coefficient and k the , for m from one to + /// k at index m - 1, of the block at + /// Q^multiplicity -- the written factor -- of + /// over the product of , where + /// Q is a quadratic over the rationals with the roots r = ρ + σ i and its + /// conjugate, ρ and σ rational. Null where Q is not such a + /// quadratic, or a side does not read as a polynomial over the rationals in the symbols. + /// + /// + /// + /// The residue formula at the root: with Q^k = l^k (x - r)^k (x - r̄)^k, the + /// coefficient c_j of 1/(x - r)^j is the (k - j)th Taylor coefficient + /// at r of N/(R l^k (x - r̄)^k), R the other factors, which is the + /// quotient of two series in x - r -- one division per order. The coefficients + /// at r̄ are the conjugates: conjugation fixes the rationals in the symbols and so + /// fixes N/D, which holds as polynomial identities whatever the symbols are. So + /// the block is Σ 2 Re(c_j (x - r̄)^j)/M^j for the monic M = Q/l, and each + /// numerator, written in powers of M, gives its digits to the powers below it. + /// + /// + /// Over the Gaussian rationals as pairs of polynomials, the unit a variable j + /// with j^2 = -1 read off as the parity of its power: every order of the series + /// is a pair over a power of the norm of the series' leading term, and nothing is put in + /// lowest terms until the last coefficient. The digits of (x - r̄)^j in powers of + /// M are numbers. + /// + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + private static (Entity P, Entity Q)[]? OverTheGaussianRoots( + Entity numerator, List<(Entity Factor, int Degree)> factors, Entity written, + Dictionary quadratic, int multiplicity, Variable x) + { + ERational? Read(int degree) + => !quadratic.TryGetValue(EInteger.FromInt32(degree), out var read) ? ERational.Zero + : read.Evaled is Rational rational ? rational.ERational : null; + if (quadratic.Keys.Any(degree => degree.Sign < 0 || degree.CompareTo(EInteger.FromInt32(2)) > 0) + || Read(2) is not { } leading || leading.IsZero || Read(1) is not { } linear || Read(0) is not { } free) + return null; + // M = x^2 + s x + t = (x - ρ)^2 + σ^2. + var s = linear.Divide(leading).ToLowestTerms(); + var t = free.Divide(leading).ToLowestTerms(); + var rho = s.Negate().Divide(ERational.FromInt32(2)).ToLowestTerms(); + var sigmaSquared = t.Subtract(rho.Multiply(rho)).ToLowestTerms(); + if (sigmaSquared.Sign <= 0) + return null; + var (sigmaAbove, sigmaBelow) = (sigmaSquared.Numerator.Sqrt(), sigmaSquared.Denominator.Sqrt()); + if (!sigmaAbove.Multiply(sigmaAbove).Equals(sigmaSquared.Numerator) || !sigmaBelow.Multiply(sigmaBelow).Equals(sigmaSquared.Denominator)) + return null; + var sigma = ERational.Create(sigmaAbove, sigmaBelow); + + // Both sides at x = r + o, the unit named: x - r̄ is o + 2σ i there. + Entity nameSource = numerator; + foreach (var (factor, _) in factors) + nameSource += factor; + var offset = Variable.CreateUnique(nameSource, "o"); + var unit = Variable.CreateUnique(nameSource + offset, "j"); + var at = Rational.Create(rho) + Rational.Create(sigma) * unit + offset; + Entity rest = MathS.Pow(offset + Rational.Create(sigma.Multiply(ERational.FromInt32(2))) * unit, multiplicity); + if (leading.CompareTo(ERational.One) != 0) + rest *= MathS.Pow(Rational.Create(leading), multiplicity); + // Each side over one bar first, the bars free of x going to the other side: a linear + // made monic is `a/b + x`, and a quotient of symbols is no polynomial to the reader. + Entity shifted = Integer.One; + bool TakeIn(Entity side, bool above) + { + var (sideBase, sidePower) = side is Powf(var raised, Integer { EInteger.Sign: > 0 } times) && times.EInteger.CanFitInInt32() + ? (raised, times.EInteger.ToInt32Unchecked()) : (side, 1); + var (overOneBar, under) = SingleQuotient.Of(SingleQuotient.Combine(sideBase).InnerSimplified); + if (under.ContainsNode(x)) + return false; + var inOffset = MathS.Pow(Bare(overOneBar).Substitute(x, at), sidePower); + if (above) + shifted *= inOffset; + else + rest *= inOffset; + if (under != Integer.One) + { + if (above) + rest *= MathS.Pow(under, sidePower); + else + shifted *= MathS.Pow(under, sidePower); + } + return true; + } + if (!TakeIn(numerator, above: true)) + return null; + foreach (var (other, _) in factors) + if (other != written && !TakeIn(other, above: false)) + return null; + var variables = shifted.Vars.Concat(rest.Vars).Distinct().OrderBy(v => v.Name, System.StringComparer.Ordinal).ToArray(); + if (variables.Length > MultivariatePolynomial.MaxVariables || variables.Contains(x)) + return null; + var indices = new Dictionary(variables.Length); + for (var i = 0; i < variables.Length; i++) + indices[variables[i]] = i; + if (MultivariatePolynomial.TryParse(shifted, indices) is not { } top || MultivariatePolynomial.TryParse(rest, indices) is not { } bottom) + return null; + var zero = MultivariatePolynomial.Zero(variables.Length); + (MultivariatePolynomial Re, MultivariatePolynomial Im)[] TaylorCoefficients(MultivariatePolynomial polynomial) + { + var byPower = polynomial.CoefficientsIn(indices[offset]); + var coefficients = new (MultivariatePolynomial, MultivariatePolynomial)[multiplicity]; + for (var order = 0; order < multiplicity; order++) + coefficients[order] = byPower.TryGetValue(order, out var atOrder) ? RealAndImaginaryParts(atOrder, indices[unit]) : (zero, zero); + return coefficients; + } + const int largest = 4096; + (MultivariatePolynomial, MultivariatePolynomial)? Times((MultivariatePolynomial Re, MultivariatePolynomial Im) left, (MultivariatePolynomial Re, MultivariatePolynomial Im) right) + => left.Re.Multiply(right.Re, largest) is { } rr && left.Im.Multiply(right.Im, largest) is { } ii + && left.Re.Multiply(right.Im, largest) is { } ri && left.Im.Multiply(right.Re, largest) is { } ir + ? (rr.Subtract(ii), ri.Add(ir)) : null; + var above = TaylorCoefficients(top); + var below = TaylorCoefficients(bottom); + // A root the other factors share is not this block's alone. + var (leadRe, leadIm) = below[0]; + if (leadRe.IsZero && leadIm.IsZero) + return null; + var conjugate = (leadRe, leadIm.ScaleBy(ERational.FromInt32(-1))); + if (leadRe.Multiply(leadRe, largest) is not { } leadReSquared || leadIm.Multiply(leadIm, largest) is not { } leadImSquared) + return null; + // The series of N/(R l^k (x - r̄)^k) at r is S_m/norm^(m + 1), order m, with + // S_m = (A_m norm^m - Σ B_n S_(m-n) norm^(n-1)) conj(B_0). + var norm = leadReSquared.Add(leadImSquared); + var normPowers = new MultivariatePolynomial[multiplicity + 1]; + normPowers[0] = MultivariatePolynomial.One(variables.Length); + for (var power = 1; power <= multiplicity; power++) + if (normPowers[power - 1].Multiply(norm, largest) is { } next) + normPowers[power] = next; + else + return null; + var series = new (MultivariatePolynomial Re, MultivariatePolynomial Im)[multiplicity]; + for (var order = 0; order < multiplicity; order++) + { + if (above[order].Re.Multiply(normPowers[order], largest) is not { } valueRe || above[order].Im.Multiply(normPowers[order], largest) is not { } valueIm) + return null; + (MultivariatePolynomial Re, MultivariatePolynomial Im) value = (valueRe, valueIm); + for (var n = 1; n <= order; n++) + { + if (Times(below[n], series[order - n]) is not var (productRe, productIm) + || productRe.Multiply(normPowers[n - 1], largest) is not { } scaledRe || productIm.Multiply(normPowers[n - 1], largest) is not { } scaledIm) + return null; + value = (value.Re.Subtract(scaledRe), value.Im.Subtract(scaledIm)); + } + if (Times(value, conjugate) is not var (seriesRe, seriesIm)) + return null; + series[order] = (seriesRe, seriesIm); + } + + // (x - r̄)^j = ((x - ρ) + σ i)^j, by its real and imaginary parts, lowest power first. + var binomial = new (ERational Re, ERational Im)[] { (ERational.One, ERational.Zero) }; + var real = new ERational[multiplicity + 1][]; + var imaginary = new ERational[multiplicity + 1][]; + for (var j = 1; j <= multiplicity; j++) + { + var next = new (ERational Re, ERational Im)[binomial.Length + 1]; + for (var m = 0; m < next.Length; m++) + next[m] = (ERational.Zero, ERational.Zero); + for (var m = 0; m < binomial.Length; m++) + { + // times x, and times -ρ + σ i + next[m + 1] = (next[m + 1].Re.Add(binomial[m].Re), next[m + 1].Im.Add(binomial[m].Im)); + next[m] = (next[m].Re.Add(binomial[m].Re.Multiply(rho.Negate())).Subtract(binomial[m].Im.Multiply(sigma)), + next[m].Im.Add(binomial[m].Im.Multiply(rho.Negate())).Add(binomial[m].Re.Multiply(sigma))); + } + binomial = next; + real[j] = binomial.Select(c => c.Re.ToLowestTerms()).ToArray(); + imaginary[j] = binomial.Select(c => c.Im.ToLowestTerms()).ToArray(); + } + + // The digits in powers of M: c_j (x - r)^(-j) and its conjugate are 2 Re(c_j (x - r̄)^j)/M^j, + // and the i-th digit of that numerator stands over M^(j - i), which is l^(j - i)/Q^(j - i), + // written l^k/(l^(k - j + i) Q^(j - i)) as the digits of Newton's iteration are. + // Every c_j is over norm^(k - j + 1), and is put over norm^k. + var numerators = new (MultivariatePolynomial Constant, MultivariatePolynomial Slope)[multiplicity + 1]; + for (var m = 1; m <= multiplicity; m++) + numerators[m] = (zero, zero); + for (var j = 1; j <= multiplicity; j++) + { + var (cRe, cIm) = series[multiplicity - j]; + if (cRe.Multiply(normPowers[j - 1], largest) is not { } overRe || cIm.Multiply(normPowers[j - 1], largest) is not { } overIm) + return null; + var realDigits = DigitsOverTheMonicQuadratic(real[j], s, t); + var imaginaryDigits = DigitsOverTheMonicQuadratic(imaginary[j], s, t); + for (var i = 0; 2 * i <= j; i++) + { + var m = j - i; + var scale = ERational.FromInt32(2).Multiply(Power(leading, multiplicity)); + var (re0, re1) = i < realDigits.Count ? realDigits[i] : (ERational.Zero, ERational.Zero); + var (im0, im1) = i < imaginaryDigits.Count ? imaginaryDigits[i] : (ERational.Zero, ERational.Zero); + numerators[m] = ( + numerators[m].Constant.Add(overRe.ScaleBy(scale.Multiply(re0))).Subtract(overIm.ScaleBy(scale.Multiply(im0))), + numerators[m].Slope.Add(overRe.ScaleBy(scale.Multiply(re1))).Subtract(overIm.ScaleBy(scale.Multiply(im1)))); + } + } + var denominator = MathS.Pow(norm.ToEntity(variables), multiplicity); + var result = new (Entity P, Entity Q)[multiplicity]; + for (var m = 1; m <= multiplicity; m++) + { + var p = numerators[m].Constant.IsZero ? Integer.Zero : InLowestTermsOverTheSymbols(numerators[m].Constant.ToEntity(variables) / denominator); + var q = numerators[m].Slope.IsZero ? Integer.Zero : InLowestTermsOverTheSymbols(numerators[m].Slope.ToEntity(variables) / denominator); + if (p.ContainsNode(x) || q.ContainsNode(x) || p.ContainsNode(unit) || q.ContainsNode(unit) || p.ContainsNode(offset) || q.ContainsNode(offset)) + return null; + result[m - 1] = (p, q); + } + return result; + } + + /// to the whole power . + private static ERational Power(ERational value, int exponent) + { + var result = ERational.One; + for (var k = 0; k < exponent; k++) + result = result.Multiply(value).ToLowestTerms(); + return result; + } + + /// + /// The digits (r0, r1), lowest first, of the polynomial over the rationals with + /// , lowest power first, in powers of the monic quadratic + /// x^2 + s x + t. + /// + private static List<(ERational, ERational)> DigitsOverTheMonicQuadratic(ERational[] coefficients, ERational s, ERational t) + { + var digits = new List<(ERational, ERational)>(); + var current = (ERational[])coefficients.Clone(); + while (current.Length > 0) + { + var remainder = (ERational[])current.Clone(); + var quotient = new ERational[System.Math.Max(0, current.Length - 2)]; + for (var k = current.Length - 1; k >= 2; k--) + { + var coefficient = remainder[k]; + quotient[k - 2] = coefficient; + remainder[k - 1] = remainder[k - 1].Subtract(coefficient.Multiply(s)).ToLowestTerms(); + remainder[k - 2] = remainder[k - 2].Subtract(coefficient.Multiply(t)).ToLowestTerms(); + } + digits.Add((remainder[0], remainder.Length > 1 ? remainder[1] : ERational.Zero)); + current = quotient; + } + return digits; + } + /// /// The first Taylor coefficients of the polynomial /// in at , diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 5392434ec..f6787daff 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -32148,6 +32148,223 @@ private static Entity OverTheFunctionLessARoot(Entity root, bool ofTheSine, Enti return -2 * MathS.Arctan(t / sigma) / (rate * (1 + root) * sigma); } + /// + /// A rational function of a sine whose denominator in it is written as powers of + /// linears -- of the sine itself, of 1 ± sin(x) from an even power of the + /// cosine, and of a + b sin(x) with symbols in it: Rubi's + /// cot(x)^6/(a + b sin(x))^3 is (1 - s^2)^3/(s^6 (a + b s)^3) in + /// s = sin(x). Split over those factors by the residues at their roots, every + /// piece is a constant over a power of one linear in the function, which the rules for + /// csc(x)^n, 1/(1 ± sin(x))^n and 1/(a + b sin(x))^n answer. + /// + /// + /// + /// Under the half-angle tangent, which took these before, every linear but the function + /// itself is a quadratic in t -- a + b s is (a t^2 + 2 b t + a)/(1 + t^2) -- + /// and the cube of a symbolic one is split by Newton's iteration modulo its powers, with + /// the symbols in every round. Rubi's cot(x)^6/(a + b sin(x))^3, + /// tan(x)^4/(a + b sin(x))^3 and 1/((a + b sin(x))^3 (c + d sin(x))^3) ran + /// past thirty seconds there. In s each block is a Taylor series at a root. + /// + /// + /// Not of a cosine: a + b cos(x) is ((a + b) + (a - b) t^2)/(1 + t^2), even in + /// t, and the half-angle answers tan(x)^6/(a + b cos(x))^3 in two seconds and + /// ten thousand characters, where the pieces here, each a piecewise of its own, came to + /// three hundred thousand. + /// + /// + /// Only where the half-angle is that slow: a linear with a symbolic root to the third power + /// or more, and the half-angle's denominator past the eleventh degree, which is 2 d - k for a + /// denominator of degree d in the function with s^k among its factors. + /// cot(x)^4/(a + b sin(x))^3, of the tenth, and cos(x)^6 csc(x)^5/(a + b sin(x))^3, + /// of the eleventh, are answered there within a second, and are left to it; past it, + /// csc(x)^2 sec(x)^2/(a + b sin(x))^3 took seven seconds there and a quarter of one here. + /// + /// + /// The decomposition is checked at sampled points by the split, and the sum of the + /// integrated pieces by its derivative before it is returned. A piece is over one linear, + /// which this declines, so nothing it asks comes back here. + /// + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveARationalFunctionOfASineOverPowersOfLinearsInIt(Entity expr, Entity.Variable x, bool integrateByParts) + { + Entity? argument = null; + foreach (var node in expr.Nodes) + { + if (TrigonometricArgument(node) is not { } thisArgument || !thisArgument.ContainsNode(x)) + continue; + if (argument is null) + argument = thisArgument; + else if (argument != thisArgument) + return null; + } + if (argument is null || !TreeAnalyzer.TryGetPolyLinear(argument, x, out var rate, out _) + || rate.ContainsNode(x) || TreeAnalyzer.IsZero(rate)) + return null; + if (!expr.Nodes.All(node => !node.ContainsNode(x) + || node is Variable or Sumf or Minusf or Mulf or Divf or Sinf or Cosf or Secantf or Cosecantf or Tanf or Cotanf + || node is Powf(_, Number.Integer))) + return null; + var sine = MathS.Sin(argument); + var cosine = MathS.Cos(argument); + var inTheTwo = expr.Replace(node => node switch + { + Secantf(var a) when a == argument => 1 / cosine, + Cosecantf(var a) when a == argument => 1 / sine, + Tanf(var a) when a == argument => sine / cosine, + Cotanf(var a) when a == argument => cosine / sine, + _ => node, + }); + var s = Variable.CreateUnique(inTheTwo, "s_trig"); + var w = Variable.CreateUnique(inTheTwo + s, "w_trig"); + var inS = inTheTwo.Substitute(sine, s).Substitute(cosine, w); + if (inS.ContainsNode(x)) + return null; + return OverPowersOfLinearsInTheFunction(inS, s, w, sine, x, integrateByParts) is { } answer + && Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) + ? answer : null; + } + + /// + /// The rule above in , which stands for , + /// for the cosine: null where the integrand is not a rational + /// function of once the even powers of are + /// 1 - s^2, or its denominator's written factors are not powers of linears in + /// and of , or it is not past the half-angle's bound. + /// + private static Entity? OverPowersOfLinearsInTheFunction(Entity inS, Variable s, Variable w, Entity function, Entity.Variable x, bool integrateByParts) + { + var combined = Functions.SingleQuotient.Combine(inS).InnerSimplified; + if (combined is Providedf(var bare, _)) + combined = bare; + var (numerator, denominator) = Functions.SingleQuotient.Of(combined); + // The powers of the function and of its complement the two sides share, cancelled + // first: `sec(x)^5 (A + B sec(x))/(a + b sec(x))^3` comes over one bar as + // `c^3 (A c + B)/(c^6 (a c + b)^3)` in `c = cos(x)`, and the bound below is on the + // degrees of the fraction, not of a spelling of it. + if (WithASharedPowerCancelled(numerator, denominator, s, even: false) is var (fewerAbove, fewerBelow)) + (numerator, denominator) = (fewerAbove, fewerBelow); + if (WithASharedPowerCancelled(numerator, denominator, w, even: true) is var (fewerAboveInW, fewerBelowInW)) + (numerator, denominator) = (fewerAboveInW, fewerBelowInW); + if (WithTheComplementSquared(numerator, s, w) is not { } above) + return null; + Entity constant = Number.Integer.One; + Entity written = Number.Integer.One; + var degree = 0; + var powerOfTheFunction = 0; + var blocks = 0; + var aSymbolicCube = false; + foreach (var factor in Mulf.LinearChildren(denominator)) + { + if (!factor.ContainsNode(s) && !factor.ContainsNode(w)) + { + constant = constant * factor; + continue; + } + var (toRead, multiplicity) = factor is Powf(var repeatedBase, Number.Integer { EInteger.Sign: > 0 } repeated) && repeated.EInteger.CanFitInInt32() + ? (repeatedBase, repeated.EInteger.ToInt32Unchecked()) : (factor, 1); + // An even power of the complement is `(1 - s)^h (1 + s)^h`. + if (toRead == w) + { + if (multiplicity % 2 != 0) + return null; + var half = multiplicity / 2; + written = written * (half == 1 ? 1 - s : MathS.Pow(1 - s, half)) * (half == 1 ? 1 + s : MathS.Pow(1 + s, half)); + degree += multiplicity; + blocks += 2; + continue; + } + if (toRead.ContainsNode(w) || !TreeAnalyzer.TryGetPolynomial(toRead, s, out var read) || read.Count == 0 + || read.Keys.Any(power => power.Sign < 0) || !read.Keys.Max()!.Equals(EInteger.One) + || read.Values.Any(coefficient => coefficient.ContainsNode(s) || coefficient.ContainsNode(w))) + return null; + if (toRead == s) + powerOfTheFunction += multiplicity; + // A linear whose root is a symbol's is a symbolic quadratic under the half-angle; + // `a + a s` is `a (1 + s)`, a square of a linear there. + aSymbolicCube |= multiplicity >= 3 + && Functions.PartialFractions.InLowestTermsOverTheSymbols(-(read.TryGetValue(EInteger.Zero, out var free) ? free : Number.Integer.Zero) / read[EInteger.One]).Vars.Any(); + written = written * factor; + degree += multiplicity; + blocks++; + } + if (!aSymbolicCube || blocks < 2 || 2 * degree - powerOfTheFunction < 12) + return null; + + // Divided down where the numerator's degree reaches the denominator's; the whole + // part is a polynomial in the function, the rules' own. + Entity total = Number.Integer.Zero; + if (TreeAnalyzer.TryGetPolynomial(above, s, out var aboveRead) && aboveRead.Count > 0 && aboveRead.Keys.Max()!.CompareTo(EInteger.FromInt32(degree)) >= 0) + { + var expandedBelow = (constant * written).Expand().InnerSimplified; + if (TreeAnalyzer.PolynomialLongDivision(above, expandedBelow, genericCase: true, inTermsOf: s) is not var (whole, rest)) + return null; + if (Integration.ComputeIndefiniteIntegral(whole.InnerSimplified.Substitute(s, function), x, integrateByParts) is not { } wholePart) + return null; + total += wholePart; + above = (rest is Divf(var top, _) ? top : rest).InnerSimplified; + if (above is Providedf(var bareRest, _)) + above = bareRest; + if (above.Evaled is Number.Complex { IsZero: true }) + return total; + } + if (!Functions.PartialFractions.TrySplitOverWrittenFactors(above, constant == Number.Integer.One ? written : constant * written, s, out var decomposition)) + return null; + var (terms, over) = decomposition is Divf(var sum, var divisor) && !divisor.ContainsNode(s) + ? (sum, divisor) : (decomposition, Number.Integer.One as Entity); + Entity overTheFactors = Number.Integer.Zero; + foreach (var term in Sumf.LinearChildren(terms)) + { + if (term.Evaled is Number.Complex { IsZero: true }) + continue; + if (Integration.ComputeIndefiniteIntegral(term.Substitute(s, function), x, integrateByParts) is not { } integrated) + return null; + overTheFactors += integrated; + } + return total + (over == Number.Integer.One ? overTheFactors : overTheFactors / over); + } + + /// + /// over with the power of + /// written among the denominator's factors and the lowest power of it + /// in the numerator cancelled as far as they share it -- an even number of times where + /// is set; null where nothing is shared. + /// + private static (Entity Numerator, Entity Denominator)? WithASharedPowerCancelled(Entity numerator, Entity denominator, Variable v, bool even) + { + var below = 0; + var rest = new List(); + foreach (var factor in Mulf.LinearChildren(denominator)) + { + if (factor == v) + below++; + else if (factor is Powf(var @base, Number.Integer { EInteger.Sign: > 0 } power) && @base == v && power.EInteger.CanFitInInt32()) + below += power.EInteger.ToInt32Unchecked(); + else + rest.Add(factor); + } + if (below == 0 || !TreeAnalyzer.TryGetPolynomial(numerator, v, out var read) || read.Count == 0 + || read.Keys.Any(power => power.Sign < 0 || !power.CanFitInInt32()) || read.Values.Any(coefficient => coefficient.ContainsNode(v))) + return null; + var taken = System.Math.Min(read.Keys.Min()!.ToInt32Unchecked(), below); + if (even) + taken -= taken % 2; + if (taken == 0) + return null; + Entity? above = null; + foreach (var pair in read) + { + var power = pair.Key.ToInt32Unchecked() - taken; + var term = power == 0 ? pair.Value : power == 1 ? pair.Value * v : pair.Value * MathS.Pow(v, power); + above = above is null ? term : above + term; + } + var left = below - taken; + if (left > 0) + rest.Add(left == 1 ? v : MathS.Pow(v, left)); + return (above!, rest.Count == 0 ? Number.Integer.One : rest.Aggregate((product, factor) => product * factor)); + } + /// /// , a polynomial in and , /// with every even power of written as the power of diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 84a708934..84948298f 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -1105,6 +1105,11 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // among the coefficients, split over the quadratic's two roots: under the // half-angle it is a quartic in t with a symbol in every coefficient. if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionOfASineOverASymbolicQuadratic(expr, x, integrateByParts)) is { }) return answer; + // And of the sine over written powers of linears in it, one with a symbolic root cubed, + // split in the sine itself: under the half-angle that linear is a symbolic quadratic in + // t, and the split of its cube took the budget. Before the substitution search, which + // with a shifted argument asked it again in its own variable and spent the budget. + if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionOfASineOverPowersOfLinearsInIt(expr, x, integrateByParts)) is { }) return answer; // And over a binomial in a higher power of it, over the binomial's roots. if ((answer = IndefiniteIntegralSolver.SolveAReciprocalOfABinomialInAPowerOfASineOverItsRoots(expr, x)) is { }) return answer; // Fractional powers of `a ± a sin(y)` beside anything rational in the sine and diff --git a/Sources/Tests/UnitTests/Calculus/PowersOfLinearsByTheirResiduesIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/PowersOfLinearsByTheirResiduesIntegralTest.cs new file mode 100644 index 000000000..9b2ae4d06 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/PowersOfLinearsByTheirResiduesIntegralTest.cs @@ -0,0 +1,72 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// Rational functions over written powers of linears, split by the residues at their roots. + /// A rational function of the sine over powers of linears in it, one with a symbolic root + /// cubed, is split in the sine itself, where under the half-angle that linear is a symbolic + /// quadratic whose cube took the budget; and a repeated quadratic over the rationals with Gaussian + /// roots, the tangent's 1 + u^2, is split by its residue at i. Rubi's 4.1.1.3, + /// 4.1.2.1, 4.1.2.2 and 4.3.1.3. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class PowersOfLinearsByTheirResiduesIntegralTest + { + [Theory] + [InlineData("tan(c + d*x)^4/(a + b*sin(c + d*x))^3")] + [InlineData("cot(c + d*x)^6/(a + b*sin(c + d*x))^3")] + [InlineData("cos(c + d*x)^6*csc(c + d*x)^8/(a + b*sin(c + d*x))^3")] + [InlineData("csc(c + d*x)^3*sec(c + d*x)^2/(a + b*sin(c + d*x))^3")] + [InlineData("1/((a + b*sin(k + f*x))^3*(c + d*sin(k + f*x))^3)")] + public void OverPowersOfLinearsInTheSine(string integrand) + { + AntiderivativeHolds(integrand, (1.3, 0.7)); + AntiderivativeHolds(integrand, (0.6, 1.7)); + } + + [Theory] + [InlineData("sin(c + d*x)^6/(a + b*tan(c + d*x))^2")] + [InlineData("sin(c + d*x)^6/(a + b*tan(c + d*x))^3")] + [InlineData("sin(c + d*x)^4/(a + b*tan(c + d*x))^4")] + [InlineData("x^6/((1 + x^2)^4*(a + b*x)^2)")] + [InlineData("x^3/((2 + 2*x + x^2)^3*(a + b*x)^2)")] + public void OverARepeatedQuadraticWithGaussianRoots(string integrand) + => AntiderivativeHolds(integrand, (1.3, 0.7)); + + private static void AntiderivativeHolds(string integrand, (double A, double B) pins) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 40000, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("a", pins.A).Substitute("b", pins.B) + .Substitute("c", 0.4).Substitute("d", 0.3).Substitute("k", 0.2).Substitute("f", 0.9); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -0.5, 0.3, 0.8, 1.3, 2.2, 3.1 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + } +}