diff --git a/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs b/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs index 85c509f02..054d5f1af 100644 --- a/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs +++ b/Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs @@ -499,7 +499,14 @@ internal static bool TrySplitOverWrittenFactors( // in every term and twenty kilobytes of answer, right and unreadable, and its // fourth power did not evaluate. Numeric coefficients keep the system, which is // exact over the rationals and answers them as before. - if (factors.All(f => TreeAnalyzer.TryGetPolynomial(f.Factor is Powf(var b, _) ? b : f.Factor, x, out var linear) && linear.Count > 0 && linear.Keys.Max()!.Equals(EInteger.One)) + // And beside a quadratic factor, the same way: the linear blocks by their Taylor + // coefficients at the root, and the quadratic's numerator in the ring modulo the + // quadratic. Rubi's `tan^4 (A + B tan)/(a + b tan)^4` is + // `u^4 (A + B u)/((a + b u)^4 (1 + u^2))` under the tangent, and the system + // answered it in `a^63 b^10` and did not evaluate within its budget; this way it + // is a line of arctangents and logarithms. + // https://github.com/asc-community/AngouriMath/issues/718 + if (factors.Any(f => IsLinear(f.Factor, x)) && (numerator + denominator).Vars.Any(v => v != x) && TrySplitOverSymbolicLinearFactors(numerator, denominator, constant, factors, x, out decomposition)) return true; @@ -585,11 +592,64 @@ internal static Entity InLowestTermsOverTheSymbols(Entity constant) if (expandedAbove == Integer.Zero || expandedAbove.Evaled is Complex { IsZero: true } || IsZeroAtPinnedSymbols(expandedAbove)) return Integer.Zero; var expandedBelow = Bare(below.Expand().InnerSimplified); - if (PolynomialGcd.TryCancel(expandedAbove, expandedBelow, out var cancelled)) + // With room for the coefficient a fourth-order block leaves: the series of + // `x^4 (A + B x)/(1 + x^2)` at `-a/b` had `(a^2 + b^2)^3` in common at the third + // order and eleven hundred nodes, and left uncancelled it carried into every + // term of the answer. + if (PolynomialGcd.TryCancel(expandedAbove, expandedBelow, out var cancelled, maxComplexity: 4096)) return Bare(cancelled); + // The gcd declines past its degree and step bounds, and a denominator that comes + // out of a series is a product of written factors -- `b^9 (a^2 + b^2)^7` -- so the + // common factor is one of them to some power: each is divided out of the + // numerator exactly for as long as that goes. + if (CancelledByTheWrittenFactors(expandedAbove, below) is { } byFactors) + return byFactors; return expandedBelow == Integer.One ? expandedAbove : expandedAbove / expandedBelow; } + /// + /// over , a product of + /// powers as written, with every written factor that divides the numerator exactly + /// divided out of both; null where none does, or where either side does not read as + /// a polynomial over the rationals. + /// + private static Entity? CancelledByTheWrittenFactors(Entity numerator, Entity denominator) + { + var variables = numerator.Vars.Concat(denominator.Vars).Distinct().OrderBy(v => v.Name, System.StringComparer.Ordinal).ToArray(); + if (variables.Length == 0 || variables.Length > MultivariatePolynomial.MaxVariables) + return null; + var indices = new Dictionary(variables.Length); + for (var i = 0; i < variables.Length; i++) + indices[variables[i]] = i; + if (MultivariatePolynomial.TryParse(numerator, indices) is not { } top) + return null; + var cancelledAny = false; + Entity left = Integer.One; + foreach (var written in Mulf.LinearChildren(denominator)) + { + var (@base, power) = written is Powf(var repeated, Integer { EInteger.Sign: > 0 } times) && times.EInteger.CanFitInInt32() + ? (repeated, times.EInteger.ToInt32Unchecked()) : (written, 1); + if (@base is Number || !@base.Vars.Any() || MultivariatePolynomial.TryParse(@base, indices) is not { } factor || factor.IsConstant) + { + left *= written; + continue; + } + while (power > 0 && top.DivideExact(factor) is { } quotient) + { + top = quotient; + power--; + cancelledAny = true; + } + if (power > 0) + left *= power == 1 ? @base : MathS.Pow(@base, power); + } + if (!cancelledAny) + return null; + var reduced = top.ToEntity(variables); + var leftExpanded = Bare(left.Expand().InnerSimplified); + return leftExpanded == Integer.One ? reduced : reduced / leftExpanded; + } + /// /// Whether , free of every variable but symbols, is zero at /// two sets of values for them, off the integers and distinct; a value that is not @@ -634,12 +694,19 @@ private static bool TryReadLinear(Entity factor, Variable x, out Entity alpha, o return true; } + /// Whether a written factor, or the base of a written power, is linear in . + private static bool IsLinear(Entity factor, Variable x) + => TreeAnalyzer.TryGetPolynomial(factor is Powf(var b, _) ? b : factor, x, out var linear) && linear.Count > 0 && linear.Keys.Max()!.Equals(EInteger.One); + /// - /// for linear factors only, each with its - /// multiplicity: the coefficient of 1/F^j for F = alpha + beta x of - /// multiplicity k is beta^(j-k) h^(k-j)(r)/(k-j)! with r = -alpha/beta - /// and h the numerator over the other factors. Checked at sampled points like - /// the system's answer. + /// for linear factors with their + /// multiplicities and quadratic factors to the first power, with a symbol somewhere. + /// The coefficient of 1/F^j for F = alpha + beta x of multiplicity + /// k is beta^(j-k) h_(k-j) with h_m the mth Taylor + /// coefficient at r = -alpha/beta of the numerator over the other factors; + /// the numerator over a quadratic is the residue of the numerator times the inverses + /// of the other factors modulo the quadratic. Checked at sampled points like the + /// system's answer. /// private static bool TrySplitOverSymbolicLinearFactors( Entity numerator, Entity denominator, Entity constant, List<(Entity Factor, int Degree)> factors, Variable x, @@ -650,42 +717,178 @@ private static bool TrySplitOverSymbolicLinearFactors( for (var i = 0; i < factors.Count; i++) { var (written, multiplicity) = factors[i]; + if (!IsLinear(written, x)) + continue; var factor = written is Powf(var repeatedBase, _) ? repeatedBase : written; if (!TreeAnalyzer.TryGetPolynomial(factor, x, out var read)) return false; var beta = read[EInteger.One]; var alpha = read.TryGetValue(EInteger.Zero, out var constantTerm) ? constantTerm : Integer.Zero; var root = Bare((-alpha / beta).InnerSimplified); - Entity rest = numerator; + Entity others = Integer.One; for (var j = 0; j < factors.Count; j++) if (j != i) - rest /= factors[j].Factor; - var derivative = rest; - var factorial = EInteger.One; + others *= factors[j].Factor; + // The Taylor coefficients of the numerator over the other factors at the + // root, by the series quotient of their own: each is a polynomial in x, so + // shifted to the root and expanded it is a polynomial in the offset with + // coefficients that are small quotients in the symbols, and the quotient of + // the two series is one division per order. Differentiating the quotient + // symbolically and evaluating at the root gives the same coefficients as + // quotients of derivatives of the other factors' product, which the gcd does + // not cancel: the third derivative of `x^4 (A + B x)/(1 + x^2)` at `-a/b` comes + // out in `b^249`. + if (TaylorCoefficientsAtTheRoot(numerator, root, multiplicity, x) is not { } above + || TaylorCoefficientsAtTheRoot(others, root, multiplicity, x) is not { } below) + return false; + // A root the other factors share is not this block's alone, and the generic + // case has nothing to say about it. + if (below[0] == Integer.Zero || below[0].Evaled is Complex { IsZero: true } || IsZeroAtPinnedSymbols(below[0])) + return false; + var series = new Entity[multiplicity]; for (var order = 0; order < multiplicity; order++) { - if (order > 0) - { - derivative = derivative.Differentiate(x); - factorial = factorial.Multiply(order); - } - var j = multiplicity - order; - var value = Bare(derivative.Substitute(x, root).InnerSimplified); - var coefficient = InLowestTermsOverTheSymbols(MathS.Pow(beta, j - multiplicity) * value / Integer.Create(factorial)); - if (coefficient.Nodes.Any(node => node == MathS.NaN) || coefficient.ContainsNode(x)) + Entity value = above[order]; + for (var k = 1; k <= order; k++) + value -= below[k] * series[order - k]; + series[order] = InLowestTermsOverTheSymbols(value / below[0]); + if (series[order].Nodes.Any(node => node == MathS.NaN) || series[order].ContainsNode(x)) return false; + } + for (var order = 0; order < multiplicity; order++) + { + var j = multiplicity - order; + var coefficient = InLowestTermsOverTheSymbols(MathS.Pow(beta, j - multiplicity) * series[order]); if (coefficient == Integer.Zero || coefficient.Evaled is Complex { IsZero: true }) continue; var term = j == 1 ? coefficient / factor : coefficient / MathS.Pow(factor, j); sum += term; } } + // A quadratic factor takes a numerator of degree one, and it is computed where + // it lives: in the ring of polynomials modulo the quadratic, where every other + // factor is a residue with an inverse, and the numerator is the numerator of the + // whole times those inverses. Two coefficients, each a small quotient in the + // symbols, and no expansion of anything but the quadratic's own reductions -- + // taking the linear blocks off the numerator and dividing by their factors + // needs expansions of nested products of `b (a^2 + b^2)^3` past the term bound. + foreach (var (written, _) in factors) + { + if (IsLinear(written, x)) + continue; + if (!TreeAnalyzer.TryGetPolynomial(written, x, out var quadratic) || !quadratic.ContainsKey(EInteger.FromInt32(2))) + return false; + var leading = quadratic[EInteger.FromInt32(2)]; + 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); + if (ResidueModuloTheQuadratic(numerator, s, t, x) is not { } residue) + return false; + foreach (var (other, _) in factors) + { + if (other == written) + continue; + if (ResidueModuloTheQuadratic(other, s, t, x) is not { } factorResidue + || InverseModuloTheQuadratic(factorResidue, s, t) is not { } inverse) + return false; + residue = ProductModuloTheQuadratic(residue, inverse, s, t); + } + var (p, q) = residue; + 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; + // The residue is of `N/others` modulo the monic quadratic `Q/leading`; over + // `Q` itself the leading coefficient cancels: `(p + q x)/(Q/leading)` is the + // term for `N/(others Q/leading)`, so `(p + q x)/Q` is the term for `N/(others Q)`. + // 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. + if (over != Integer.Zero) + { + 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); + sum += below == Integer.One ? above / written : MathS.Pow(below, Integer.MinusOne) * (above / written); + } + } if (!HoldsAtSampledPoints(numerator / denominator, sum / constant, x)) return false; decomposition = sum / constant; return true; } + /// + /// The first Taylor coefficients of the polynomial + /// in at , + /// each in lowest terms over the symbols; null where it does not read as one. + /// + private static Entity[]? TaylorCoefficientsAtTheRoot(Entity polynomial, Entity root, int count, Variable x) + { + var offset = Variable.CreateUnique(polynomial + root, "t"); + var shifted = Bare(polynomial.Substitute(x, root + offset).Expand().InnerSimplified); + if (!TreeAnalyzer.TryGetPolynomial(shifted, offset, out var read)) + return null; + var coefficients = new Entity[count]; + for (var order = 0; order < count; order++) + { + var coefficient = read.TryGetValue(EInteger.FromInt32(order), out var at) ? InLowestTermsOverTheSymbols(at) : Integer.Zero; + if (coefficient.ContainsNode(x) || coefficient.ContainsNode(offset)) + return null; + coefficients[order] = coefficient; + } + return coefficients; + } + + /// + /// reduced modulo the monic quadratic + /// x^2 + s x + t, as the pair (p, q) of p + q x, each in lowest + /// terms over the symbols; null where it does not read as a polynomial in + /// . + /// + private static (Entity, Entity)? ResidueModuloTheQuadratic(Entity polynomial, Entity s, Entity t, Variable x) + { + if (!TreeAnalyzer.TryGetPolynomial(polynomial, x, out var read) || read.Count == 0) + return null; + var degree = read.Keys.Max()!; + if (!degree.CanFitInInt32() || read.Keys.Any(power => power.Sign < 0)) + return null; + // Horner's scheme, with `x^2` written as `-s x - t` at each step. + (Entity, Entity) residue = (Integer.Zero, Integer.Zero); + for (var power = degree.ToInt32Unchecked(); power >= 0; power--) + { + residue = ProductModuloTheQuadratic(residue, (Integer.Zero, Integer.One), s, t); + if (read.TryGetValue(EInteger.FromInt32(power), out var coefficient)) + residue = (InLowestTermsOverTheSymbols(residue.Item1 + coefficient), residue.Item2); + } + return residue; + } + + /// The product of two residues modulo x^2 + s x + t, in lowest terms. + private static (Entity, Entity) ProductModuloTheQuadratic((Entity, Entity) left, (Entity, Entity) right, Entity s, Entity t) + { + var (a, b) = left; + var (c, d) = right; + // (a + b x)(c + d x) = ac + (ad + bc) x + bd x^2, and x^2 is -s x - t. + return (InLowestTermsOverTheSymbols(a * c - b * d * t), InLowestTermsOverTheSymbols(a * d + b * c - b * d * s)); + } + + /// + /// The inverse of a residue modulo x^2 + s x + t, or null where it has none + /// generically: its norm a^2 - s a b + t b^2 is zero, which is the residue + /// sharing a root with the quadratic. + /// + private static (Entity, Entity)? InverseModuloTheQuadratic((Entity, Entity) residue, Entity s, Entity t) + { + var (a, b) = residue; + // (a + b x)(c + d x) = 1 is ac - t b d = 1 and a d + b c - s b d = 0. + var norm = InLowestTermsOverTheSymbols(a * a - s * a * b + t * b * b); + if (norm == Integer.Zero || norm.Evaled is Complex { IsZero: true } || IsZeroAtPinnedSymbols(norm)) + return null; + return (InLowestTermsOverTheSymbols((a - s * b) / norm), InLowestTermsOverTheSymbols(-b / norm)); + } + /// /// Gaussian elimination on a square system whose entries may be symbols. A pivot is a /// number that is not zero where one is available in its column, and otherwise an entry diff --git a/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialGcd.cs b/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialGcd.cs index ee3cc164a..fd3482a14 100644 --- a/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialGcd.cs +++ b/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialGcd.cs @@ -76,7 +76,7 @@ internal static class PolynomialGcd /// x + 1 provided not x - 1 = 0. /// internal static bool TryCancel(Entity numerator, Entity denominator, - [NotNullWhen(true)] out Entity? cancelled) + [NotNullWhen(true)] out Entity? cancelled, int maxComplexity = MaxComplexity) { cancelled = null; @@ -84,7 +84,11 @@ internal static bool TryCancel(Entity numerator, Entity denominator, // whatever they are, and there is no point parsing them to find that out. if (!numerator.Vars.Any(denominator.Vars.Contains)) return false; - if (numerator.Complexity + denominator.Complexity > MaxComplexity) + // The bound is the simplifier's, for a quotient it meets on every rewrite; a + // caller that has built one quotient it needs in lowest terms -- the + // coefficients of a partial fraction, whose whole point is to be small -- may + // raise it. + if (numerator.Complexity + denominator.Complexity > maxComplexity) return false; var variables = numerator.Vars diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index e220d1dad..297f0ac25 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -605,7 +605,7 @@ is var (multiple, leftover) if (Mulf.LinearChildren(denominator).Any(f => f.ContainsNode(x) && f is Powf(_, Number.Integer { EInteger.Sign: > 0 } e) && e != Number.Integer.One) && IsAProductOfSymbolicLinearFactors(denominator, x) && Functions.PartialFractions.TrySplitOverWrittenFactors(numerator, denominator, x, out var overSymbolicLinears) - && Integration.ComputeIndefiniteIntegral(overSymbolicLinears, x, integrateByParts) is { } overTheLinears) + && IntegratedTermByTerm(overSymbolicLinears, x, integrateByParts) is { } overTheLinears) return overTheLinears; // Every rule below reads the denominator **as written**: the Hermite reduction wants @@ -698,7 +698,7 @@ is var (multiple, leftover) // factors -- `(x + a)(x^2 + b)` -- is decomposed by undetermined coefficients, // checked, and each piece is a shape the rules here read. if (Functions.PartialFractions.TrySplitOverWrittenFactors(numerator, denominator, x, out var overWrittenFactors) - && Integration.ComputeIndefiniteIntegral(overWrittenFactors, x, integrateByParts) is { } termByTerm) + && IntegratedTermByTerm(overWrittenFactors, x, integrateByParts) is { } termByTerm) return termByTerm; // Last, because everything above answers in exact arithmetic where it can: a @@ -9168,25 +9168,69 @@ private static Entity WithLogarithmsOfProductsApart(Entity expr, Entity.Variable return sum; }); + /// + /// A partial-fraction decomposition, a sum over a constant, integrated one term at a + /// time as each is written, with the constant on the result. + /// + /// + /// Handed to the chain whole, the sum is split as written only while every term is + /// small (), and past that by the expanding + /// gather, which combines a term's constant with its quotient: `D^(-1) (P + Q x)/(1 + x^2)` + /// with `P`, `Q` polynomials in the symbols became `(P + Q x)/(D + D x^2)`, a quadratic + /// with symbols in every coefficient, and was integrated as a piecewise on the sign of + /// its discriminant. A decomposition's terms are already the shapes the rules read. + /// + private static Entity? IntegratedTermByTerm(Entity decomposition, Entity.Variable x, bool integrateByParts) + { + var (terms, over) = decomposition is Divf(var sum, var constant) && !constant.ContainsNode(x) + ? (sum, constant) : (decomposition, Number.Integer.One as Entity); + Entity total = Number.Integer.Zero; + foreach (var term in Sumf.LinearChildren(terms)) + { + if (term.Evaled is Number.Complex { IsZero: true }) + continue; + if (Integration.ComputeIndefiniteIntegral(term, x, integrateByParts) is not { } integrated) + return null; + total += integrated; + } + return over == Number.Integer.One ? total : total / over; + } + /// /// Whether is written as a product of two or more - /// distinct linear factors in , to whole powers, with a symbol in - /// a coefficient somewhere. + /// distinct factors in , at least one of them linear and to a + /// whole power, the rest linear or quadratic to the first power, with a symbol in a + /// coefficient somewhere. /// + /// + /// A quadratic beside the linears is allowed since the decomposition takes the linear + /// blocks by their Taylor coefficients and the quadratic's numerator in the ring + /// modulo the quadratic; with the gate asking for linears only, + /// u^4 (A + B u)/((a + b u)^4 (1 + u^2)) -- every rational function of the + /// tangent with a power of a linear in it -- went to the Hermite reduction below, + /// which answered in a^63 b^10. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// private static bool IsAProductOfSymbolicLinearFactors(Entity denominator, Entity.Variable x) { var linears = 0; + var quadratics = 0; foreach (var factor in Mulf.LinearChildren(denominator)) { if (!factor.ContainsNode(x)) continue; var @base = factor is Powf(var b, Number.Integer { EInteger.Sign: > 0 }) ? b : factor; - if (!TreeAnalyzer.TryGetPolynomial(@base, x, out var read) || read.Count == 0 || !read.Keys.Max()!.Equals(EInteger.One) + if (!TreeAnalyzer.TryGetPolynomial(@base, x, out var read) || read.Count == 0 || read.Values.Any(coefficient => coefficient.ContainsNode(x))) return false; - linears++; + if (read.Keys.Max()!.Equals(EInteger.One)) + linears++; + else if (read.Keys.Max()!.Equals(EInteger.FromInt32(2)) && @base == factor) + quadratics++; + else + return false; } - return linears >= 2 && denominator.Vars.Any(v => v != x); + return linears >= 1 && linears + quadratics >= 2 && denominator.Vars.Any(v => v != x); } /// diff --git a/Sources/Tests/UnitTests/Calculus/SymbolicFactorsPartialFractionsTest.cs b/Sources/Tests/UnitTests/Calculus/SymbolicFactorsPartialFractionsTest.cs index 70bff6699..0f0efc095 100644 --- a/Sources/Tests/UnitTests/Calculus/SymbolicFactorsPartialFractionsTest.cs +++ b/Sources/Tests/UnitTests/Calculus/SymbolicFactorsPartialFractionsTest.cs @@ -173,6 +173,25 @@ public void RepeatedSymbolicLinearFactorsAreDecomposedByDerivatives(string integ Assert.True(integrand.ToEntity().Integrate("x").Stringize().Length < 4000, "the answer is a page"); } + /// + /// And beside a quadratic factor the linear blocks are still taken by the derivatives, + /// with what they leave of the numerator over the quadratic alone: Rubi's + /// tan^4 (A + B tan)/(a + b tan)^4 is u^4 (A + B u)/((a + b u)^4 (1 + u^2)) + /// under the tangent, and solved as a system it came out in a^63 b^10 and did + /// not evaluate within the corpus's budget. + /// + [Theory] + [InlineData("x^4*(A + B*x)/((a + b*x)^4*(1 + x^2))")] + [InlineData("1/((a + b*x)^2*(1 + x^2))")] + [InlineData("x/((a + b*x)*(x^2 + 2))")] + [InlineData("x^3/((a + b*x)^2*(1 + x^2)*(4 + x^2))")] + [InlineData("tan(x)^4*(A + B*tan(x))/(a + b*tan(x))^4")] + public void SymbolicLinearBlocksBesideAQuadratic(string integrand) + { + DifferentiatesBack(integrand, "x", ("a", 1.7), ("b", 2.3), ("A", 0.4), ("B", 1.1)); + Assert.True(integrand.ToEntity().Integrate("x").Stringize().Length < 4000, "the answer is a page"); + } + /// /// Two linear factors with one root are one factor: (a + b x)(a + x b)^2 is /// how the rules that make a factor monic and gather its powers write it, and as