A biquadratic denominator is decomposed over the reals (#233)

BREAKING-CHANGES.md

@@ -37,6 +37,7 @@ read first.
 | | `new Entity[0].MultiplyAll()`, and `Mulf.Multiply` on an empty list | `AngouriBugException` | `1` |
 | | `MathS.Vector()` and `new Entity[0].ToVector()` | `IndexOutOfRangeException` — outside the documented hierarchy | `InvalidMatrixOperationException` |
 | | `"x3 - 1".ToEntity().Factorize()`, and every polynomial no rewrite rule has a rule for | `x ^ 3 - 1` — handed back whole | `(x - 1) * (x ^ 2 + x + 1)` |
+| | `"x^2/(x^4 + 1)".Integrate("x")`, and every quotient whose denominator is a biquadratic irreducible over `Q` | `integral(x ^ 2 / (x ^ 4 + 1), x)` — left unevaluated | the antiderivative, over the real quadratic factors |
 | | `Entity.DomainConditionIn(Domain)` | did not exist | the domain of definition for a **stated** reading, through the whole tree |
 | | `MathS.Polynomials.Factor("(y * x3 + 1) * (x4 - y3)", "x")`, and bivariate polynomials whose leading coefficient in the main variable is a polynomial | `null` — a refusal | `(x ^ 4 - y ^ 3) * (x ^ 3 * y + 1)` |
 | | `MathS.Polynomials.Factor("x7 - y7", "x")`, and bivariate polynomials whose substituted image over-factors | `null` — a refusal | `(x - y) * (x ^ 6 + x ^ 5 * y + … + y ^ 6)` |
@@ -2344,13 +2345,66 @@ no integration rule reads: an irreducible of degree three or more, or a quadrati
 `(x^2 + 1)^2` is declined, and the ladder that would decompose it is deliberately not built, because
 every term it produces is over `(x^2 + c)^k` and would come back unevaluated in turn.
 
+*The `x^4 + 1` half of that is no longer true: allowing those real coefficients is exactly what the
+entry below does, and `x^2/(x^4 + 1)` is now answered. The rest of the paragraph stands — an
+irreducible of degree three or more, and a repeated quadratic, are still declined, and `(x^2 + 1)^2`
+still is.*
+
 Deciding that from the factorisation rather than by trying is what keeps the cost of declining to the
 one factorisation. Splitting regardless and recursing made `(1 - x^4)/(1 + x^4 + x^8)`, whose
 factorisation holds the irreducible quartic `x^4 - x^2 + 1`, take 18s to return the same unevaluated
 integral it returns in 203ms — measured, and the reason the guard is there rather than a preference.
 
 [#919](https://github.com/asc-community/AngouriMath/issues/919).
 
+### A biquadratic denominator is decomposed over the reals
+
+The same blindness one level further along, and the last place it reaches. The step above factors
+over `Q` and stops where `Q` does, so `x^4 + 1` — irreducible over the rationals — was left whole
+and `x^2/(x^4 + 1)` came back unevaluated. Over the reals it is
+`(x^2 - sqrt(2)x + 1)(x^2 + sqrt(2)x + 1)`, and both halves are read by the rule for a linear
+numerator over a quadratic. Nothing was missing but a factorisation the rational step is right to
+refuse.
+
+```
+"x^2/(x^4 + 1)".Integrate("x")
+
+was  integral(x ^ 2 / (x ^ 4 + 1), x)
+is   -1/2 * sqrt(2) * 1/2 / 2 * ln(x ^ 2 + sqrt(2) * x + 1)
+     + 1/2 * arctan((2 * x + sqrt(2)) * 1/2 * sqrt(2)) * 1/2 * sqrt(2)
+     + 1/2 * sqrt(2) * 1/2 / 2 * ln(x ^ 2 - sqrt(2) * x + 1)
+     + 1/2 * arctan((2 * x + -sqrt(2)) * 1/2 * sqrt(2)) * 1/2 * sqrt(2) + C
+```
+
+This is the integral [#233](https://github.com/asc-community/AngouriMath/issues/233) names as
+wanting "partial fractioning", and it is the first of that issue's list to need a factorisation
+rather than a rule. `1/(x^4 + 1)`, `1/(x^4 - 2)` and `1/(x^4 + 3x^2 + 1)` come with it.
+
+**Biquadratic only, and that is a boundary rather than a first cut.** A general quartic factors into
+real quadratics through its resolvent cubic, whose roots carry Cardano's nested radicals; a
+biquadratic `x^4 + px^2 + q` is the case where the resolvent is solvable by inspection and the two
+factors stay inside one square root. Two shapes come out of it, by the sign of `p^2 - 4q`: negative
+gives `(x^2 + ax + b)(x^2 - ax + b)` with `b = sqrt(q)` and `a = sqrt(2b - p)`, and positive gives
+the even `(x^2 + u)(x^2 + v)` with `u, v = (p -+ sqrt(p^2 - 4q))/2`. Zero is `(x^2 + p/2)^2`, a
+repeated quadratic, declined for the reason the step above declines one. **A quartic with an odd
+power in it — `x^4 + x^3 + 1`, `x^4 + x + 1` — is still declined**, and so is everything of degree
+five and up that does not factor over `Q`.
+
+**No condition is attached**, on the same argument as the step above: the two factors are distinct,
+so their product is zero exactly where the original denominator is, and nothing is cancelled.
+
+It is tried **after** both rational steps, which is what keeps a denominator that factors over `Q`
+in exact arithmetic: `x^4 + 3x^2 + 2` is decomposed by the step above and never arrives here to be
+given a square root it does not need. Declining stays as cheap as it was —
+`(1 - x^4)/(1 + x^4 + x^8)` returns the same unevaluated integral in the same fraction of a second,
+because every guard here is rational arithmetic on coefficients already read.
+
+`sqrt(tan(x))`, the one remaining entry on #233's list, is **not** answered by this. It reduces
+under `u = sqrt(tan x)` to `2 * integral(u^2/(u^4 + 1), u)`, which is now integrable — but the
+substitution that gets there is a separate capability and is not built here.
+
+[#233](https://github.com/asc-community/AngouriMath/issues/233).
+
 ### A polynomial equation that factors is solved through its factors
 
 A polynomial of degree four or more may factor over the rationals with no rational root anywhere in

Sources/AngouriMath/Functions/Algebra/Polynomials/PartialFractions.cs

@@ -5,6 +5,7 @@
 // Website: https://am.angouri.org.
 //
 
+using PeterO.Numbers;
 using System.Diagnostics.CodeAnalysis;
 using static AngouriMath.Entity;
 using static AngouriMath.Entity.Number;
@@ -133,6 +134,181 @@ internal static bool TrySplitIntoCoprimeParts(
             return true;
         }
 
+        /// <summary>
+        /// <c>N/D</c> written as two fractions over the quadratic factors of a <b>biquadratic</b>
+        /// <paramref name="denominator"/> — one with no odd power in it — or
+        /// <see langword="false"/> where it is not one, or does not split into two distinct
+        /// quadratics.
+        /// </summary>
+        /// <remarks>
+        /// <para>
+        /// <b>Why this exists next to the step above.</b> That one factors over the rationals,
+        /// and stops where the rationals do: <c>x^4 + 1</c> is irreducible over <c>Q</c>, so it
+        /// is left whole and <c>x^2/(x^4 + 1)</c> has no antiderivative — the case
+        /// <a href="https://github.com/asc-community/AngouriMath/issues/233">#233</a> names as
+        /// wanting "partial fractioning". Over the reals it is
+        /// <c>(x^2 - sqrt(2)x + 1)(x^2 + sqrt(2)x + 1)</c>, and both halves are read by the rule
+        /// for a linear numerator over a quadratic. Nothing was missing but a factorisation the
+        /// rational one is right to refuse.
+        /// </para>
+        /// <para>
+        /// <b>Biquadratic only, and that is a real boundary rather than a first cut.</b> A
+        /// general quartic factors into real quadratics through its resolvent cubic, whose roots
+        /// carry Cardano's nested radicals; a biquadratic <c>x^4 + px^2 + q</c> is the case where
+        /// the resolvent is solvable by inspection, and the two factors stay in one square root.
+        /// Two shapes come out of it, by the sign of <c>p^2 - 4q</c>:
+        /// </para>
+        /// <list type="bullet">
+        /// <item>
+        /// <b>Negative</b> — no real root in <c>x^2</c>. Matching
+        /// <c>(x^2 + ax + b)(x^2 - ax + b) = x^4 + (2b - a^2)x^2 + b^2</c> gives
+        /// <c>b = sqrt(q)</c> and <c>a = sqrt(2b - p)</c>, both real because <c>q &gt; 0</c> and
+        /// <c>p^2 &lt; 4q</c> forces <c>p &lt; 2sqrt(q)</c>. This is <c>x^4 + 1</c>, at
+        /// <c>a = sqrt(2)</c>, <c>b = 1</c>.
+        /// </item>
+        /// <item>
+        /// <b>Positive</b> — two distinct real roots in <c>x^2</c>, so
+        /// <c>(x^2 + M)(x^2 + N)</c> with <c>M, N = (p +- sqrt(p^2 - 4q))/2</c>. Both factors are
+        /// even, and the split is two independent pairs of equations rather than four.
+        /// </item>
+        /// <item>
+        /// <b>Zero</b> — <c>(x^2 + p/2)^2</c>, a repeated quadratic, declined for the same reason
+        /// the guard above declines one: there is no rule for a numerator over
+        /// <c>(x^2 + c)^k</c>, so decomposing it ends in the integral it started from.
+        /// </item>
+        /// </list>
+        /// <para>
+        /// <b>No condition is owed</b>, on the same argument as the step above: the two factors
+        /// are distinct and coprime, so their product vanishes exactly where the original
+        /// denominator does. <c>q &gt; 0</c> is required rather than assumed, which is what keeps
+        /// <c>b</c> real; a negative <c>q</c> puts a real root in <c>x^2</c> of either sign and
+        /// is left to the rational step, which reaches it whenever the root is rational.
+        /// </para>
+        /// <para>
+        /// Reached only after the rational split has declined, so a biquadratic that factors over
+        /// <c>Q</c> — <c>x^4 + 3x^2 + 2</c> — is decomposed there, in exact arithmetic, and never
+        /// arrives here to be given a square root it does not need.
+        /// </para>
+        /// </remarks>
+        internal static bool TrySplitBiquadraticOverTheReals(
+            Entity numerator, Entity denominator, Variable x,
+            [NotNullWhen(true)] out Entity? left,
+            [NotNullWhen(true)] out Entity? right)
+        {
+            left = right = null;
+
+            // Every guard below is rational arithmetic on coefficients already in hand, so a
+            // denominator this does not apply to costs one polynomial read to decline.
+            if (!PolynomialFactoring.TryGetRationalCoefficients(
+                    denominator, x, leastTerms: 2, leastDegree: 4, maxDegree: 4, out var d)
+                || d.Length != 5
+                || !d[1].IsZero || !d[3].IsZero)
+                return false;
+
+            // A proper fraction only, as above: an improper one is a polynomial plus a proper
+            // fraction and has to be divided out first, which is not done here. The degree
+            // ceiling of three is what says so, the denominator's being four.
+            if (!PolynomialFactoring.TryGetRationalCoefficients(
+                    numerator, x, leastTerms: 1, leastDegree: 0, maxDegree: 3, out var c))
+                return false;
+
+            var lead = d[4];
+            var p = d[2].Divide(lead);
+            var q = d[0].Divide(lead);
+
+            // At q = 0 the quartic is x^2(x^2 + p), whose rational root zero the step above has
+            // already had. Nothing else about the sign of q is required here: the branch that
+            // needs sqrt(q) real is the one below with a negative discriminant, and p^2 < 4q
+            // makes q positive on its own.
+            if (q.IsZero)
+                return false;
+
+            var discriminant = p.Multiply(p).Subtract(q.Multiply(ERational.FromInt32(4)));
+            if (discriminant.IsZero)
+                return false;
+
+            // The numerator, padded to four coefficients so the two branches can index it
+            // without asking how many terms it happened to have.
+            var n = new Entity[4];
+            for (var i = 0; i < n.Length; i++)
+                n[i] = Rational.Create(i < c.Length ? c[i] : ERational.Zero);
+
+            Entity leftNumerator, leftFactor, rightNumerator, rightFactor;
+            if (discriminant.Sign < 0)
+            {
+                // (x^2 + ax + b)(x^2 - ax + b). Writing the split as (alpha x + beta)/A +
+                // (gamma x + delta)/B and equating the four coefficients of
+                // (alpha x + beta)B + (gamma x + delta)A against the numerator gives, using
+                // B's -a where A has +a:
+                //
+                //   x^3:  alpha + gamma          = n3
+                //   x^2:  a(gamma - alpha) + beta + delta = n2
+                //   x^1:  b(alpha + gamma) + a(delta - beta) = n1
+                //   x^0:  b(beta + delta)        = n0
+                //
+                // which is two sums and two differences rather than a linear solve.
+                var b = MathS.Sqrt(Rational.Create(q)).InnerSimplified;
+                var a = MathS.Sqrt(2 * b - Rational.Create(p)).InnerSimplified;
+
+                var sum = (n[0] / b).InnerSimplified;                       // beta + delta
+                var difference = ((n[1] - b * n[3]) / a).InnerSimplified;   // delta - beta
+                var spread = ((n[2] - sum) / a).InnerSimplified;            // gamma - alpha
+
+                leftFactor = MathS.Sqr(x) + a * x + b;
+                rightFactor = MathS.Sqr(x) - a * x + b;
+                leftNumerator = ((n[3] - spread) / 2 * x + (sum - difference) / 2).InnerSimplified;
+                rightNumerator = ((n[3] + spread) / 2 * x + (sum + difference) / 2).InnerSimplified;
+            }
+            else
+            {
+                // (x^2 + u)(x^2 + v), both even, so the odd and even halves of the numerator
+                // separate and each gives its own pair rather than one system of four. Writing
+                // the split as (alpha x + beta)/(x^2 + u) + (gamma x + delta)/(x^2 + v), the
+                // coefficients of (alpha x + beta)(x^2 + v) + (gamma x + delta)(x^2 + u) are
+                //
+                //   x^3: alpha + gamma = n3     x^1: v*alpha + u*gamma = n1
+                //   x^2: beta  + delta = n2     x^0: v*beta  + u*delta = n0
+                //
+                // so alpha = (n1 - u*n3)/(v - u) and beta = (n0 - u*n2)/(v - u), with v - u the
+                // square root of the discriminant. Note which of the two the numerators divide
+                // by: pairing a numerator with the wrong factor flips the sign of the answer
+                // and still satisfies the x^3 and x^2 rows, so it is not something the identity
+                // check further down would catch on every numerator.
+                var root = MathS.Sqrt(Rational.Create(discriminant)).InnerSimplified;
+                var v = ((Rational.Create(p) + root) / 2).InnerSimplified;
+                var u = ((Rational.Create(p) - root) / 2).InnerSimplified;
+
+                var alpha = ((n[1] - u * n[3]) / root).InnerSimplified;
+                var beta = ((n[0] - u * n[2]) / root).InnerSimplified;
+
+                leftFactor = MathS.Sqr(x) + u;
+                rightFactor = MathS.Sqr(x) + v;
+                leftNumerator = (alpha * x + beta).InnerSimplified;
+                rightNumerator = ((n[3] - alpha) * x + (n[2] - beta)).InnerSimplified;
+            }
+
+            // Two identities, checked rather than assumed -- the numerators against the numerator
+            // they decompose, and the factors against the denominator they came from. The step
+            // above checks one because its factorisation is exact by construction; here the
+            // factors were built by matching coefficients through a square root, so the
+            // factorisation is a claim of its own.
+            //
+            // Neither implies the other, so both are made. A wrong term common to the two factors
+            // -- MathS.Sqr(x) misread as C#'s x ^ 2, which on an Entity is exclusive or and not a
+            // power -- cancels between the two halves of the first identity and passes it, while
+            // the second sees it immediately.
+            if ((leftNumerator * rightFactor + rightNumerator * leftFactor
+                    - numerator).Simplify() != Integer.Create(0))
+                return false;
+            if ((Rational.Create(lead) * leftFactor * rightFactor - denominator).Simplify()
+                    != Integer.Create(0))
+                return false;
+
+            left = (leftNumerator / (Rational.Create(lead) * leftFactor)).InnerSimplified;
+            right = (rightNumerator / (Rational.Create(lead) * rightFactor)).InnerSimplified;
+            return true;
+        }
+
         /// <summary>
         /// A vanishing numerator is answered as zero rather than as a quotient, so that a
         /// numerator sharing a factor with the denominator does not leave the integrator

Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs

@@ -45,6 +45,15 @@ internal static class IndefiniteIntegralSolver
         /// applies and leaves something that will not integrate — it takes the denominator
         /// apart differently, so a failure of the one is no evidence about the other.
         /// </para>
+        /// <para>
+        /// Both of those factor over the rationals and stop where the rationals do, which left
+        /// x^2/(x^4 + 1) unevaluated: x^4 + 1 is irreducible over Q, and over the reals it is
+        /// (x^2 - sqrt(2)x + 1)(x^2 + sqrt(2)x + 1). The third step
+        /// (<see cref="Functions.PartialFractions.TrySplitBiquadraticOverTheReals"/>) reads a
+        /// biquadratic denominator that way. It is last because it is the only one that
+        /// introduces a radical, and a denominator that factors over Q should be taken apart in
+        /// exact arithmetic by one of the two above.
+        /// </para>
         /// </remarks>
         internal static Entity? SolveByPartialFractions(Entity expr, Entity.Variable x, bool integrateByParts)
         {
@@ -63,6 +72,12 @@ internal static class IndefiniteIntegralSolver
                 && Integration.ComputeIndefiniteIntegral(right, x, integrateByParts) is { } overOther)
                 return overOne + overOther;
 
+            if (Functions.PartialFractions.TrySplitBiquadraticOverTheReals(
+                    numerator, denominator, x, out var overOneReal, out var overOtherReal)
+                && Integration.ComputeIndefiniteIntegral(overOneReal, x, integrateByParts) is { } realFirst
+                && Integration.ComputeIndefiniteIntegral(overOtherReal, x, integrateByParts) is { } realRest)
+                return realFirst + realRest;
+
             return null;
         }