Add monic polynomial (#453)

* Add monic polynomial

* Adjust monic to be normal function instead of return ref-counted struct

* Remove commented Monic struct

* Fixups for monic polynomial

* Add poly function

* Remove poly/monic slice overloads and reconfigure unit tests

* Remove monic

* Unify poly into a single function

* Unify poly with Polynomial

* Adjust poly signature to scope const F

* Update polynomial.d

* Update polynomial.d

---------

Co-authored-by: Ilia Ki <ki9ilia@gmail.com>

Author: jmh530

source/mir/polynomial.d

@@ -41,27 +41,7 @@ struct Polynomial(F)
         +/
         @optmath typeof(F.init * X.init * 1f + F.init) opCall(X)(in X x) const
         {
-            import mir.internal.utility: Iota;
-            auto ret = cast(typeof(return))0;
-            if (coefficients)
-            {
-                ptrdiff_t i = coefficients.length - 1;
-                assert(i >= 0);
-                auto c = cast()coefficients[i];
-                static foreach (d; Iota!derivative)
-                    c *= i - d;
-                ret = cast(typeof(return)) c;
-                while (--i >= cast(ptrdiff_t)derivative)
-                {
-                    assert(i < coefficients.length);
-                    c = cast()coefficients[i];
-                    static foreach (d; Iota!derivative)
-                        c *= i - d;
-                    ret *= x;
-                    ret += c;
-                }
-            }
-            return ret;
+            return x.poly!derivative(this.coefficients[]);
         }
     }
 }
@@ -81,9 +61,9 @@ version (mir_test) @safe pure nothrow @nogc unittest
     auto a = rcarray!(const double)(3.0, 4.5, 1.9, 2);
     auto p = a.polynomial;
 
-    alias f = (x) => 3.0 + 4.5 * x^^1 + 1.9 * x^^2 + 2 * x^^3;
-    alias df = (x) => 4.5 + 2 * 1.9 * x^^1 + 3 * 2 * x^^2;
-    alias d2f = (x) => 2 * 1.9  + 6 * 2 * x^^1;
+    alias   f = (x) =>     3.0 +     4.5 * x^^1 +   1.9 * x^^2 + 2 * x^^3;
+    alias  df = (x) =>     4.5 + 2 * 1.9 * x^^1 + 3 * 2 * x^^2;
+    alias d2f = (x) => 2 * 1.9 + 6 *   2 * x^^1;
 
     assert(p(3.3).approxEqual(f(3.3)));
     assert(p(7.2).approxEqual(f(7.2)));
@@ -94,3 +74,153 @@ version (mir_test) @safe pure nothrow @nogc unittest
     assert(p.opCall!2(3.3).approxEqual(d2f(3.3)));
     assert(p.opCall!2(7.2).approxEqual(d2f(7.2)));
 }
+
+/++
+Evaluate polynomial.
+
+Coefficients assumed to be in the order a0 + a1 * x ^^ 1 + ... + aN * x ^^ N
+
+Params:
+    F = controls type of output
+    derivative = order of derivatives (default = 0)
+
+Returns:
+    Value of the polynomial, evaluated at `x`
+
+See_also:
+    $(WEB en.wikipedia.org/wiki/Polynomial, Polynomial).
++/
+template poly(uint derivative = 0)
+{
+    import std.traits: ForeachType;
+    /++
+    Params:
+        x = value to evaluate
+        coefficients = coefficients of polynomial
+    +/
+    @optmath typeof(F.init * X.init * 1f + F.init) poly(X, F)(in X x, scope const F[] coefficients...)
+    {
+        import mir.internal.utility: Iota;
+        auto ret = cast(typeof(return))0;
+        if (coefficients.length > 0)
+        {
+            ptrdiff_t i = coefficients.length - 1;
+            assert(i >= 0);
+            auto c = cast()coefficients[i];
+            static foreach (d; Iota!derivative)
+                c *= i - d;
+            ret = cast(typeof(return)) c;
+            while (--i >= cast(ptrdiff_t)derivative)
+            {
+                assert(i < coefficients.length);
+                c = cast()coefficients[i];
+                static foreach (d; Iota!derivative)
+                    c *= i - d;
+                ret *= x;
+                ret += c;
+            }
+        }
+        return ret;
+    }
+}
+
+///
+version (mir_test) @safe pure nothrow unittest
+{
+    import mir.math.common: approxEqual;
+
+    double[] x = [3.0, 4.5, 1.9, 2];
+
+    alias   f = (x) =>     3.0 +     4.5 * x^^1 +   1.9 * x^^2 + 2 * x^^3;
+    alias  df = (x) =>     4.5 + 2 * 1.9 * x^^1 + 3 * 2 * x^^2;
+    alias d2f = (x) => 2 * 1.9 + 6 *   2 * x^^1;
+
+    assert(poly(3.3, x).approxEqual(f(3.3)));
+    assert(poly(7.2, x).approxEqual(f(7.2)));
+
+    assert(poly!1(3.3, x).approxEqual(df(3.3)));
+    assert(poly!1(7.2, x).approxEqual(df(7.2)));
+
+    assert(poly!2(3.3, x).approxEqual(d2f(3.3)));
+    assert(poly!2(7.2, x).approxEqual(d2f(7.2)));
+}
+
+// static array test
+version (mir_test) @safe pure @nogc nothrow unittest
+{
+    import mir.math.common: approxEqual;
+
+    double[4] x = [3.0, 4.5, 1.9, 2];
+
+    alias   f = (x) =>     3.0 +     4.5 * x^^1 +   1.9 * x^^2 + 2 * x^^3;
+    alias  df = (x) =>     4.5 + 2 * 1.9 * x^^1 + 3 * 2 * x^^2;
+    alias d2f = (x) => 2 * 1.9 + 6 *   2 * x^^1;
+
+    assert(poly(3.3, x).approxEqual(f(3.3)));
+    assert(poly(7.2, x).approxEqual(f(7.2)));
+
+    assert(poly!1(3.3, x).approxEqual(df(3.3)));
+    assert(poly!1(7.2, x).approxEqual(df(7.2)));
+
+    assert(poly!2(3.3, x).approxEqual(d2f(3.3)));
+    assert(poly!2(7.2, x).approxEqual(d2f(7.2)));
+}
+
+// Check coefficient.length = 3
+version (mir_test) @safe pure nothrow unittest
+{
+    import mir.math.common: approxEqual;
+
+    double[] x = [3.0, 4.5, 1.9];
+
+    alias   f = (x) =>     3.0 +     4.5 * x^^1 + 1.9 * x^^2;
+    alias  df = (x) =>     4.5 + 2 * 1.9 * x^^1;
+    alias d2f = (x) => 2 * 1.9;
+
+    assert(poly(3.3, x).approxEqual(f(3.3)));
+    assert(poly(7.2, x).approxEqual(f(7.2)));
+
+    assert(poly!1(3.3, x).approxEqual(df(3.3)));
+    assert(poly!1(7.2, x).approxEqual(df(7.2)));
+
+    assert(poly!2(3.3, x).approxEqual(d2f(3.3)));
+    assert(poly!2(7.2, x).approxEqual(d2f(7.2)));
+}
+
+// Check coefficient.length = 2
+version (mir_test) @safe pure nothrow unittest
+{
+    import mir.math.common: approxEqual;
+
+    double[] x = [3.0, 4.5];
+
+    alias   f = (x) => 3.0 + 4.5 * x^^1;
+    alias  df = (x) => 4.5;
+    alias d2f = (x) => 0.0;
+
+    assert(poly(3.3, x).approxEqual(f(3.3)));
+    assert(poly(7.2, x).approxEqual(f(7.2)));
+
+    assert(poly!1(3.3, x).approxEqual(df(3.3)));
+    assert(poly!1(7.2, x).approxEqual(df(7.2)));
+
+    assert(poly!2(3.3, x).approxEqual(d2f(3.3)));
+    assert(poly!2(7.2, x).approxEqual(d2f(7.2)));
+}
+
+// Check coefficient.length = 1
+version (mir_test) @safe pure nothrow unittest
+{
+    import mir.math.common: approxEqual;
+
+    double[] x = [3.0];
+
+    alias  f = (x) => 3.0;
+    alias df = (x) => 0.0;
+
+    assert(poly(3.3, x).approxEqual(f(3.3)));
+    assert(poly(7.2, x).approxEqual(f(7.2)));
+
+    assert(poly!1(3.3, x).approxEqual(df(3.3)));
+    assert(poly!1(7.2, x).approxEqual(df(7.2)));
+}