diff --git a/polynomials/polycomplex.11.cpp b/polynomials/polycomplex.11.cpp index 7bb7946..ee53459 100644 --- a/polynomials/polycomplex.11.cpp +++ b/polynomials/polycomplex.11.cpp @@ -193,13 +193,78 @@ struct poly_complex { x *= scale; return rn; } + + poly_complex derive(){ + int n = size(); + poly_complex b(n-1); + for (int i = 1; i < n; ++i){ + b[i-1] = a[i] * C(i); + } + + if (n == 1) + return b = {0}; + + return b; + } + + poly_complex integrate(){ + int n = size(); + poly_complex b(n+1); + b[0] = 0; + for (int i = 0; i < n; ++i){ + b[i+1] = a[i] / C(i+1.0); + } + + return b; + } + + void trim(int n){ + if ((int)a.size() > n) + a.resize(n); + } + + // @n=size() - power of two + // @a[0] = 0 + poly_complex poly_exp(){ + poly_complex f; f.a = {1}; + poly_complex g; g.a = {1}; + int m = 1, n = size(); + + while (m <= n/2){ + g = (g * C(2) - f * g * g); + g.trim(m); + + auto q = derive(); + q.trim(m-1); + + + auto w = q + g * (f.derive() - f * q); + w.trim(2*m-1); + + f = f + f * (*this - w.integrate()); + f.trim(2*m); + + m *= 2; + } + return f; + } }; /*snippet-end*/ int main() { poly_complex a; - a.a = {1, 2, 3}; - a *= a; + a.a = {0,{1,1},0,0}; + + auto t = a.poly_exp(); + auto sol = vector> {{1,0}, {1,1}, {0,1}, {-1.0/3.0, 1.0/3.0}}; - return a.a != vector> {1, 4, 10, 12, 9}; + if ((int)t.size() != (int)sol.size()) + return 1; + + for (int i = 0; i < (int)t.size(); ++i) + if (abs(sol[i] - t.a[i]) >= 1e-12) + return 1; + + return 0; } +