Add poly_complex_exp

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<double> a;
-    a.a = {1, 2, 3};
-    a *= a;
+    a.a = {0,{1,1},0,0};
+
+    auto t = a.poly_exp();
+    auto sol = vector<complex<double>> {{1,0}, {1,1}, {0,1}, {-1.0/3.0, 1.0/3.0}};
 
-    return a.a != vector<complex<double>> {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;
 }
+