Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 34 additions & 0 deletions BREAKING-CHANGES.md
Original file line number Diff line number Diff line change
Expand Up @@ -273,6 +273,40 @@ Each of these was checked by differentiating it back with the parameters pinned
four points. The Rubi sample is unchanged at 231 of 463 with no wrong answers, so nothing that
already had an antiderivative moves.

### A root of an even power is the modulus, and is no longer read as the power

**Wrong answers, silent.** `(u^2)^(3/2)` is `|u|^3`, and two readers took it for `u^3`: the
polynomial parser read `(x^q)^v` as `x^(q v)` for any `v`, and the integrator's `sin^p cos^q`
reader read `(sin^2)^(3/2)` as `sin^3`. Through `Simplify`'s long division the first made
`u^2` of `(u^2)^(3/2)/u` and `-1` of `sqrt(u^2)/(-u)`, wrong for every negative `u`; through
the second, six of Rubi's family-4 integrands were integrated as if the sign of the function
were positive on the whole line.
[#1387](https://github.com/asc-community/AngouriMath/issues/1387)

| | Was | Is |
|---|---|---|
| `"(u^2)^(3/2)/u".Simplify()` | `u^2 provided not u = 0` | `(u^2)^(3/2)/u`, unchanged |
| `"sqrt(u^2)/(-u)".Simplify()` | `-1 provided not u = 0` | unchanged |
| `"(u^3)^2/u".Simplify()`, a whole outer power | `u^5 provided not u = 0` | the same |
| `"x*sqrt(sin(x)^2)".Integrate("x")` | `-x cos(x) + sin(x)`, the integral of `x sin(x)` | unevaluated |
| `"(csc(x)^2)^(3/2)".Integrate("x")`, `"1/(csc(x)^2)^(7/2)".Integrate("x")` | the integral of `csc^3`, of `sin^7` | unevaluated |
| `"(a*sin(x)^2)^(5/2)".Integrate("x")`, `"1/sqrt(a*cot(x)^2)".Integrate("x")` | the integral of `a^(5/2) sin^5`, of `tan/sqrt(a)` | the antiderivative, with `sgn(sin x)` and `sgn(tan x)` in it |
| `"(b*tan(x)^3)^(3/2)".Integrate("x")`, an odd inner power | the integral of `b^(3/2) tan^(9/2)`: right where the tangent is positive, which is where the integrand is real, and off the principal branch where it is negative (`(u^3)^(3/2)` is `i` at `u = -1`, `u^(9/2)` is `-i`) | unevaluated |
| `"sec(c+d*x)^(3/2)/(b*sec(c+d*x))^(5/2)".Integrate("x")` | `sin(c+d*x)/(b^(5/2) d)`, right for `b > 0` | unevaluated |
| `"asin(sqrt(1 - x^2))/sqrt(1 - x^2)".Integrate("x")` | `-arcsin(sqrt(1 - x^2))^2/2`, right for `x > 0` and off by a sign for `x < 0` | unevaluated |
| `"asin(x/a)^(3/2)/sqrt(a^2 - x^2)".Integrate("x")` | `arcsin(x/a)^(5/2)/(5/2)`, right for `a > 0` | unevaluated |

A whole outer power still multiplies, since `(x^q)^3` is `x^q x^q x^q` however `x` is signed; a
fractional one over a first power is still read, `(a x)^(3/2)` as `a^(3/2) x^(3/2)`, the generic
reading every rule uses; over an odd power it is refused too, since `(u^3)^(3/2)` and `u^(9/2)`
differ at every negative `u`, and the tangent row is that reading integrated. The tangent,
`asin(x/a)` and `sec` rows are the cost — each right where the integrand is real, or for a
positive symbol, and only there: `a sqrt(1 - x^2/a^2)` was being collected to `sqrt(a^2 - x^2)` through
`(a^2)^(1/2) = a`, which is the same identity taken on a symbol, and it is not taken on a symbol
either now; and the secant quotient was simplified to `cos/(b^(5/2) d)` through
`(1/cos)^(3/2) = cos^(-3/2)`, which is wrong for every negative cosine and cancelled against the
same wrong step on the denominator. On the 1774-problem independent suites this is 1707 to 1705, the two `asin` rows, with the sign error gone.

---

## 2.5.0 — since 2.4.0
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -3339,6 +3339,13 @@ bool Read(Entity node, ERational multiplicity)
// whether the exponent it carries along is whole — the power rule takes
// <c>u^(-3/2)</c> as readily as <c>u^(-3)</c>. Which of the three cases
// applies still turns on a whole exponent, and that is checked there.
// Not a fractional power of an even power, though: `(sin^2)^(3/2)` is
// `|sin|^3`, and read as `sin^3` it was integrated as one, wrong on every
// other half-turn; `sqrt(a sin^2)^5`, `1/sqrt(a cot^2)`, `(csc^2)^(3/2)`,
// `x sqrt(sin^2)` were four of Rubi's answered so.
// https://github.com/asc-community/AngouriMath/issues/1387
case Powf(var @base, Number.Rational power) when power is not Number.Integer && HasAnEvenPowerOfATrigonometricFunction(@base):
return false;
case Powf(var @base, Number.Rational power):
return Read(@base, multiplicity * power.ERational);
// A rational factor rides along; anything else is declined rather than
Expand All @@ -3358,6 +3365,14 @@ bool Read(Entity node, ERational multiplicity)
}
}

/// <summary>
/// Whether an even whole power of a sine, cosine, tangent, cotangent, secant or
/// cosecant stands anywhere in <paramref name="expr"/>: a fractional power of it is a
/// power of the function's modulus, not of the function.
/// </summary>
private static bool HasAnEvenPowerOfATrigonometricFunction(Entity expr)
=> expr.Nodes.Any(node => node is Powf(TrigonometricFunction, Number.Integer even) && even.EInteger.IsEven && !even.EInteger.IsZero);

/// <summary>
/// Whether <paramref name="value"/> is a whole number small enough to count with, and
/// what that number is.
Expand Down Expand Up @@ -4475,6 +4490,15 @@ node is Powf(var radical, Number.Rational e) && e is not Number.Integer
return null;
if (!TreeAnalyzer.TryGetPolyLinear(first, x, out var a, out var b) || !TreeAnalyzer.TryGetPolyLinear(second, x, out var c, out var d))
return null;
// Two proportional linears are one radical with a constant in it, and not this
// rule's: their determinant a d - b c is zero, the second linear in t is 0/(a - c t^q)
// and everything it multiplies vanished -- Rubi's
// `sin(a + b (c + d x)^(1/3))/(c e + d e x)^(1/3)` came back as `0 provided ...`.
// https://github.com/asc-community/AngouriMath/issues/1386
var determinant = (a * d - b * c).InnerSimplified;
if (determinant.Evaled is Number.Complex { IsZero: true }
|| determinant.Vars.Any() && Functions.PartialFractions.Bare(determinant.Simplify()).Evaled is Number.Complex { IsZero: true })
return null;

var t = Variable.CreateUnique(expr, "t_rad");
var w = Variable.CreateUnique(expr, "w_rad");
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -294,6 +294,15 @@ internal static (Dictionary<T, Entity> poly, Entity rem) GatherMonomialInformati
}


/// <summary>Whether a monomial's power, integer, rational or decimal, is one.</summary>
private static bool IsOne<T>(T power) => power switch
{
EInteger whole => whole.Equals(EInteger.One),
ERational rational => rational.Equals(ERational.One),
EDecimal decimalPower => decimalPower.CompareTo(EDecimal.One) == 0,
_ => false,
};

internal static (Entity FreeMono, TPrimitive Power)? ParseMonomial<T, TPrimitive>(Variable aVar, Entity expr)
where TPrimitive : TreeAnalyzer.IPrimitive<T>, new()
{
Expand All @@ -315,8 +324,18 @@ internal static (Entity FreeMono, TPrimitive Power)? ParseMonomial<T, TPrimitive
return null;
if (mp == aVar)
power.Add(value);
// (x^q)^v is x^(q v) for a whole v, and for a fractional v only where x^q
// is not negative: (x^2)^(3/2) is |x|^3, and read as x^3 it let the long
// division make u^2 of (u^2)^(3/2)/u, whereupon (sin(x)^2)^(3/2) under
// u = cos(x) was integrated as sin(x)^3, wrong on every other half-turn.
// Any q but one is refused, an odd one too: (u^3)^(3/2)/u is not u^(7/2)
// at u = -1, where the first is i and the second -i, and Simplify does not
// change a value off the real line either (#752).
// https://github.com/asc-community/AngouriMath/issues/1387
else if (ParseMonomial<T, TPrimitive>(aVar, @base) is var (tmpFree, q))
{
if (value is not Integer && !IsOne(q.Value))
return null;
freeMono *= MathS.Pow(tmpFree, value);
power.AddMp(q.Value, value);
}
Expand Down
51 changes: 51 additions & 0 deletions Sources/Tests/UnitTests/Calculus/TrigonometricPowerProductTest.cs
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,57 @@ private static void DifferentiatesBack(string integrand)
[InlineData("sin(2*x)^3*cos(2*x)^2")]
public void ALinearArgument(string integrand) => DifferentiatesBack(integrand);

/// <summary>
/// A fractional power of an even power is a power of the modulus, not of the
/// function: <c>(sin^2)^(3/2)</c> is <c>|sin|^3</c>, and read as <c>sin^3</c> it was
/// integrated as one, wrong on every other half-turn; and <c>(2 tan^3)^(3/2)</c>, an
/// odd inner power, was integrated as <c>2^(3/2) tan^(9/2)</c>, which agrees only where
/// the tangent is positive. Each of these is either declined or right at points on
/// both sides of a zero of the function.
/// https://github.com/asc-community/AngouriMath/issues/1387
/// </summary>
[Theory]
[InlineData("(sin(x)^2)^(3/2)")]
[InlineData("(2*tan(x)^3)^(3/2)")]
[InlineData("x*sqrt(sin(x)^2)")]
[InlineData("(csc(x)^2)^(3/2)")]
[InlineData("1/(csc(x)^2)^(7/2)")]
[InlineData("(2*sin(x)^2)^(5/2)")]
[InlineData("1/sqrt(3*cot(x)^2)")]
public void ARootOfAnEvenPowerIsNotThePower(string integrand)
{
var integral = integrand.ToEntity().Integrate("x");
Assert.DoesNotContain("NaN", integral.Stringize());
if (integral.Stringize().Contains("integral("))
return;
var derivative = integral.Substitute("C", 0).Differentiate("x");
var original = integrand.ToEntity();
foreach (var at in new[] { 0.7, 2.1, 3.9, 4.5, 5.6 })
{
var got = derivative.Substitute("x", at).EvalNumerical();
var want = original.Substitute("x", at).EvalNumerical();
if (got.IsNaN || want.IsNaN)
continue;
var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart);
var scale = Math.Max(1.0, Math.Abs((double)want.RealPart));
Assert.True(difference / scale < 1e-9, $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
}
}

/// <summary>
/// Two proportional linears under roots are one radical with a constant in it: the
/// quotient-of-two-radicals rule took them, its determinant was zero, and
/// <c>sin(a + b (c + d x)^(1/3))/(c e + d e x)^(1/3)</c> came back as
/// <c>0 provided ...</c>. https://github.com/asc-community/AngouriMath/issues/1386
/// </summary>
[Fact]
public void TwoProportionalRadicalsAreNotAQuotientOfTwo()
{
var integral = "sin(a + b*(c + d*x)^(1/3))/(c*e2 + d*e2*x)^(1/3)".ToEntity().Integrate("x");
Assert.DoesNotContain("NaN", integral.Stringize());
Assert.NotEqual("0", (integral.Substitute("C", 0) is Entity.Providedf(var inner, _) ? inner : integral.Substitute("C", 0)).Simplify().Stringize());
}

/// <summary>
/// The exponent that is <b>not</b> the one deciding the case need only be rational. It
/// rides through the substitution untouched, and the power rule takes <c>u^(-3/2)</c> as
Expand Down
41 changes: 41 additions & 0 deletions Sources/Tests/UnitTests/PatternsTest/PowerProductBranchTest.cs
Original file line number Diff line number Diff line change
Expand Up @@ -86,5 +86,46 @@ public void TheSoundCasesStillGather(string expr, string expected)
var simplified = expr.ToEntity().Simplify();
Assert.Equal(expected.ToEntity().Simplify(), simplified);
}

/// <summary>
/// A third leak of the same identity: the polynomial parser read <c>(u^2)^(3/2)</c> as
/// <c>u^3</c>, and through the long-division rule <c>(u^2)^(3/2)/u</c> simplified to
/// <c>u^2</c> and <c>sqrt(u^2)/(-u)</c> to <c>-1</c>, wrong for every negative <c>u</c>.
/// An odd inner power is not read either: <c>(u^3)^(3/2)</c> is <c>i</c> at <c>u = -1</c>
/// and <c>u^(9/2)</c> is <c>-i</c> there. Checked at negative and positive values
/// against the expression it came from.
/// https://github.com/asc-community/AngouriMath/issues/1387
/// </summary>
[Theory]
[InlineData("(u^2)^(3/2)/u")]
[InlineData("sqrt(u^2)/(-u)")]
[InlineData("(u^2)^(5/2)/u^3")]
[InlineData("(u^4)^(1/2)*u")]
[InlineData("(u^3)^(3/2)/u")]
[InlineData("(u^3)^(1/2)/u")]
public void ARootOfAnEvenPowerKeepsItsValue(string expr)
{
var simplified = expr.ToEntity().Simplify();
foreach (var u in new[] { -2.0, -0.5, 0.7, 3.0 })
{
var before = expr.ToEntity().Substitute("u", u).EvalNumerical();
var after = simplified.Substitute("u", u).EvalNumerical();
var difference = Math.Abs(before.RealPart.EDecimal.ToDouble() - after.RealPart.EDecimal.ToDouble())
+ Math.Abs(before.ImaginaryPart.EDecimal.ToDouble() - after.ImaginaryPart.EDecimal.ToDouble());
Assert.True(difference < 1e-9,
$"{expr} simplified to {simplified.Stringize()}, which at u = {u} is {after.Stringize()} rather than {before.Stringize()}");
}
}

/// <summary>A whole outer power still multiplies, whatever the sign of the base.</summary>
[Theory]
[InlineData("(u^3)^2/u", "u^5")]
[InlineData("(u^2)^3/u", "u^5")]
public void AWholeOuterPowerStillMultiplies(string expr, string expected)
{
var simplified = expr.ToEntity().Simplify();
var bare = simplified is Entity.Providedf(var inner, _) ? inner : simplified;
Assert.Equal(expected.ToEntity().Simplify(), bare.Simplify());
}
}
}
Loading