|
| 1 | +// |
| 2 | +// Copyright (c) 2019-2026 Angouri. |
| 3 | +// AngouriMath is licensed under MIT. |
| 4 | +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. |
| 5 | +// Website: https://am.angouri.org. |
| 6 | +// |
| 7 | + |
| 8 | +using PeterO.Numbers; |
| 9 | +using static AngouriMath.Entity; |
| 10 | +using static AngouriMath.Entity.Number; |
| 11 | + |
| 12 | +namespace AngouriMath.Functions |
| 13 | +{ |
| 14 | + /// <summary> |
| 15 | + /// A product whose body is a monomial in the index, written in closed form: |
| 16 | + /// <c>product(k, k, 1, n)</c> is <c>factorial(n)</c>. |
| 17 | + /// </summary> |
| 18 | + /// <remarks> |
| 19 | + /// <para> |
| 20 | + /// The sibling of <see cref="PolynomialSummation"/> and a narrower thing than it, because a |
| 21 | + /// product has no linearity to take apart: a sum of two terms is the sum of their sums, and a |
| 22 | + /// product of two terms is not the product of their products in any way that helps. What is |
| 23 | + /// left is the body that <em>is</em> a single term — <c>c * k^p</c> — over which the product |
| 24 | + /// separates into a power of the constant and a power of the factorial. |
| 25 | + /// </para> |
| 26 | + /// <para> |
| 27 | + /// <b>The condition is the same one, for the same reason.</b> An empty range multiplies to |
| 28 | + /// <c>1</c>, so <c>product(k, k, 1, n)</c> is <c>1</c> at every <c>n < 1</c>, while |
| 29 | + /// <c>factorial(n)</c> is not — it is undefined at the negative integers, being the gamma |
| 30 | + /// function's poles. Answering the one with the other unconditionally would turn a value into |
| 31 | + /// an undefinedness, which is the failure the contract's O4 is about. The identity holds |
| 32 | + /// where <c>to >= from - 1</c>, and that is what is attached. |
| 33 | + /// </para> |
| 34 | + /// <para> |
| 35 | + /// <b>A positive lower bound where the index is in the body</b>, and that one cannot go in the |
| 36 | + /// condition. <c>product(k, k, a, b)</c> is <c>b! / (a-1)!</c> only where <c>a >= 1</c>; below |
| 37 | + /// that the range runs through zero and the product is <c>0</c>, while <c>(a-1)!</c> is |
| 38 | + /// undefined. But <c>a < 1</c> does not mean the range is empty, so it cannot share a |
| 39 | + /// branch with the empty-range case — a piecewise saying "identity otherwise" would be wrong |
| 40 | + /// there. So it is decided before anything is built, and a lower bound that is not a concrete |
| 41 | + /// integer of at least one is declined instead. A constant body has no such restriction, |
| 42 | + /// there being no factorial in its answer. |
| 43 | + /// </para> |
| 44 | + /// </remarks> |
| 45 | + internal static class MonomialProduct |
| 46 | + { |
| 47 | + /// <summary> |
| 48 | + /// The largest power of the index this will take. The answer carries the factorial to |
| 49 | + /// that power, so the ceiling is about how big an expression is worth returning. |
| 50 | + /// </summary> |
| 51 | + private const int MaxPower = 16; |
| 52 | + |
| 53 | + /// <summary> |
| 54 | + /// <c>product(expression, index, from, to)</c> written in closed form, or |
| 55 | + /// <see langword="null"/> where the body is not a monomial in the index, or the bounds |
| 56 | + /// are not ones the identity speaks about. |
| 57 | + /// </summary> |
| 58 | + internal static Entity? ClosedForm(Entity expression, Entity var, Entity from, Entity to) |
| 59 | + { |
| 60 | + if (var is not Variable index) |
| 61 | + return null; |
| 62 | + if (!PolynomialSummation.IsWholeOrSymbolic(from) || !PolynomialSummation.IsWholeOrSymbolic(to)) |
| 63 | + return null; |
| 64 | + if (!TreeAnalyzer.TryGetPolynomial(expression, index, out var monomials)) |
| 65 | + return null; |
| 66 | + |
| 67 | + // One term, or this is not a monomial and the product does not separate. |
| 68 | + if (monomials.Count != 1) |
| 69 | + return null; |
| 70 | + var only = System.Linq.Enumerable.Single(monomials); |
| 71 | + if (only.Key.Sign < 0 || !only.Key.CanFitInInt32()) |
| 72 | + return null; |
| 73 | + var power = only.Key.ToInt32Checked(); |
| 74 | + if (power > MaxPower) |
| 75 | + return null; |
| 76 | + var coefficient = only.Value; |
| 77 | + if (coefficient.ContainsNode(index)) |
| 78 | + return null; |
| 79 | + |
| 80 | + // How many times the body is multiplied: to - from + 1, which the condition below |
| 81 | + // guarantees is not negative. |
| 82 | + var count = (to - from + 1).InnerSimplified; |
| 83 | + Entity closed = MathS.Pow(coefficient, count); |
| 84 | + |
| 85 | + if (power > 0) |
| 86 | + { |
| 87 | + // b! / (a-1)! needs a - 1 to be a whole number that factorial is defined at, |
| 88 | + // which is decided here rather than asked in the condition -- see the remarks. |
| 89 | + if (from.Evaled is not Integer lower || lower.EInteger.Sign < 1) |
| 90 | + return null; |
| 91 | + var factorials = (MathS.Factorial(to) / MathS.Factorial(from - 1)).InnerSimplified; |
| 92 | + closed *= power == 1 ? factorials : MathS.Pow(factorials, Integer.Create(power)); |
| 93 | + } |
| 94 | + |
| 95 | + // `to >= from` and not `to >= from - 1`, which is where the sum's condition sits. |
| 96 | + // The two agree everywhere but the empty range itself, and there the closed form is |
| 97 | + // c^0 -- which is 1 for every c except 0, where it is undefined, while the empty |
| 98 | + // product is 1 for every c including 0. Handing the empty range to the identity |
| 99 | + // branch instead keeps a value from becoming an undefinedness, and costs nothing: |
| 100 | + // both branches say 1 there. |
| 101 | + return MathS.Piecewise(new[] |
| 102 | + { |
| 103 | + new Providedf(closed.InnerSimplified, new GreaterOrEqualf(to, from)), |
| 104 | + new Providedf(Integer.Create(1), Entity.Boolean.True), |
| 105 | + }).InnerSimplified; |
| 106 | + } |
| 107 | + } |
| 108 | +} |
0 commit comments