|
| 1 | +/++ |
| 2 | +Polynomial ref-counted structure. |
| 3 | +
|
| 4 | +License: $(HTTP boost.org/LICENSE_1_0.txt, Boost License 1.0). |
| 5 | +Authors: Ilya Yaroshenko |
| 6 | ++/ |
| 7 | +module mir.polynomial; |
| 8 | + |
| 9 | +import mir.math.common: optmath; |
| 10 | +import mir.rc.array; |
| 11 | + |
| 12 | +@optmath: |
| 13 | + |
| 14 | +/++ |
| 15 | +Polynomial callable ref-counted structure. |
| 16 | ++/ |
| 17 | +struct Polynomial(F) |
| 18 | +{ |
| 19 | + /// |
| 20 | + RCArray!(const F) coefficients; |
| 21 | + |
| 22 | + /++ |
| 23 | + Params: |
| 24 | + coefficients = coefficients `c[i]` for polynomial function `f(x)=c[0]+c[1]*x^^1+...+c[n]*x^^n` |
| 25 | + +/ |
| 26 | + this(RCArray!(const F) coefficients) |
| 27 | + { |
| 28 | + import core.lifetime: move; |
| 29 | + this.coefficients = coefficients.move; |
| 30 | + } |
| 31 | + |
| 32 | + /++ |
| 33 | + Params: |
| 34 | + derivative = derivative order |
| 35 | + +/ |
| 36 | + template opCall(uint derivative = 0) |
| 37 | + { |
| 38 | + /++ |
| 39 | + Params: |
| 40 | + x = `x` point |
| 41 | + +/ |
| 42 | + @optmath typeof(F.init * X.init * 1f + F.init) opCall(X)(in X x) |
| 43 | + { |
| 44 | + import mir.internal.utility: Iota; |
| 45 | + auto ret = cast(typeof(return))0; |
| 46 | + if (coefficients) |
| 47 | + { |
| 48 | + ptrdiff_t i = coefficients.length - 1; |
| 49 | + auto c = cast()coefficients[i]; |
| 50 | + static foreach (d; Iota!derivative) |
| 51 | + c *= i - d; |
| 52 | + ret = cast(typeof(return)) c; |
| 53 | + while (--i >= derivative) |
| 54 | + { |
| 55 | + c = cast()coefficients[i]; |
| 56 | + static foreach (d; Iota!derivative) |
| 57 | + c *= i - d; |
| 58 | + ret *= x; |
| 59 | + ret += c; |
| 60 | + } |
| 61 | + } |
| 62 | + return ret; |
| 63 | + } |
| 64 | + } |
| 65 | +} |
| 66 | + |
| 67 | +/// ditto |
| 68 | +Polynomial!F polynomial(F)(RCArray!(const F) coefficients) |
| 69 | +{ |
| 70 | + import core.lifetime: move; |
| 71 | + return typeof(return)(coefficients.move); |
| 72 | +} |
| 73 | + |
| 74 | +/// |
| 75 | +version (mir_test) @safe pure nothrow @nogc unittest |
| 76 | +{ |
| 77 | + import mir.math.common: approxEqual; |
| 78 | + import mir.rc.array; |
| 79 | + auto a = rcarray!(const double)(3.0, 4.5, 1.9, 2); |
| 80 | + auto p = a.polynomial; |
| 81 | + |
| 82 | + alias f = (x) => 3.0 + 4.5 * x^^1 + 1.9 * x^^2 + 2 * x^^3; |
| 83 | + alias df = (x) => 4.5 + 2 * 1.9 * x^^1 + 3 * 2 * x^^2; |
| 84 | + alias d2f = (x) => 2 * 1.9 + 6 * 2 * x^^1; |
| 85 | + |
| 86 | + assert(p(3.3).approxEqual(f(3.3))); |
| 87 | + assert(p(7.2).approxEqual(f(7.2))); |
| 88 | + |
| 89 | + assert(p.opCall!1(3.3).approxEqual(df(3.3))); |
| 90 | + assert(p.opCall!1(7.2).approxEqual(df(7.2))); |
| 91 | + |
| 92 | + assert(p.opCall!2(3.3).approxEqual(d2f(3.3))); |
| 93 | + assert(p.opCall!2(7.2).approxEqual(d2f(7.2))); |
| 94 | +} |
0 commit comments