fork download
copy 
  1. #include <iostream>
  2. #include <vector>
  3. #include <cmath>
  4.  
  5. using namespace std;
  6.  
  7. ostream &operator<<(ostream &o, vector<double> a) {
  8.   o << '[';
  9.   for (size_t i = 0; i < a.size(); i++) {
  10.     if (i > 0) o << ',';
  11.     o << a[i];
  12.   }
  13.   o << ']';
  14.   return o;
  15. }
  16.  
  17. vector<double> dadd(vector<double> a, vector<double> b) {
  18.   const int n = a.size();
  19.   for (int i = 1; i < n; i++) {
  20.     a[i] += b[i];
  21.   }
  22.   return a;
  23. }
  24.  
  25. vector<double> dsub(vector<double> a, vector<double> b) {
  26.   const int n = a.size();
  27.   for (int i = 1; i < n; i++) {
  28.     a[i] -= b[i];
  29.   }
  30.   return a;
  31. }
  32.  
  33. vector<double> dmul(vector<double> a, vector<double> b) {
  34.   const int n = a.size();
  35.   vector<double> c(n);
  36.   for (int i = 1; i < n; i++) {
  37.     for (int j = 1; i * j < n; j++) {
  38.       c[i * j] += a[i] * b[j];
  39.     }
  40.   }
  41.   return c;
  42. }
  43.  
  44. vector<double> ddiff(vector<double> a) {
  45.   const int n = a.size();
  46.   a[1] = 0;
  47.   for (int i = 2; i < n; i++) {
  48.     a[i] *= -log(i);
  49.   }
  50.   return a;
  51. }
  52.  
  53. vector<double> dint(vector<double> a) {
  54.   const int n = a.size();
  55.   for (int i = 2; i < n; i++) {
  56.     a[i] /= -log(i);
  57.   }
  58.   return a;
  59. }
  60.  
  61. /*
  62.   a[1]b[1] = 1
  63.   a[1]b[2] + a[2]b[1] = 0
  64.   a[1]b[3] + a[3]b[1] = 0
  65.   a[1]b[4] + a[2]b[2] + a[4]b[1] = 0
  66.   a[1]b[5] + a[5]b[1] = 0
  67.   a[1]b[6] + a[2]b[3] + a[3]b[2] + a[6]b[1] = 0
  68. */
  69. vector<double> dinv(vector<double> a) {
  70.   const int n = a.size();
  71.   vector<double> b(n);
  72.   b[1] = 1 / a[1];
  73.   for (int i = 2; i < n; i++) {
  74.     for (int j = 1; j < i; j++) {
  75.       if (i % j == 0) {
  76.         b[i] -= b[j] * a[i / j];
  77.       }
  78.     }
  79.     b[i] /= a[1];
  80.   }
  81.   return b;
  82. }
  83.  
  84. /*
  85.   g(x) = log(f(x))
  86.   g'(x) = f'(x)/f(x)
  87. */
  88. vector<double> dlog(vector<double> a) {
  89.   vector<double> b = dint(dmul(ddiff(a), dinv(a)));
  90.   return b;
  91. }
  92.  
  93. /*
  94.   g(x) = exp(f(x))
  95.   g'(x) = f'(x)g(x)
  96.  
  97.   g[1] = 1
  98.   -log(1)g[1] = f'[1]g[1]
  99.   -log(2)g[2] = f'[1]g[2] + f'[2]g[1]
  100.   -log(3)g[3] = f'[1]g[3] + f'[3]g[1]
  101.   -log(4)g[4] = f'[1]g[4] + f'[2]g[2] + f'[4]g[1]
  102.   -log(5)g[5] = f'[1]g[5] + f'[5]g[1]
  103.   -log(6)g[6] = f'[1]g[6] + f'[2]g[3] + f'[3]g[2] + f'[6]g[1]
  104.  
  105.   -log(1)g[1] = -log(1)f[1]g[1]
  106.   -log(2)g[2] = -log(1)f[1]g[2] - log(2)f[2]g[1]
  107.   -log(3)g[3] = -log(1)f[1]g[3] - log(3)f[3]g[1]
  108.   -log(4)g[4] = -log(1)f[1]g[4] - log(2)f[2]g[2] - log(4)f[4]g[1]
  109.   -log(5)g[5] = -log(1)f[1]g[5] - log(5)f[5]g[1]
  110.   -log(6)g[6] = -log(1)f[1]g[6] - log(2)f[2]g[3] - log(3)f[3]g[2] - log(6)f[6]g[1]
  111. */
  112. vector<double> dexp(vector<double> a) {
  113.   const int n = a.size();
  114.   vector<double> b(n);
  115.   b[1] = 1;
  116.   a = ddiff(a);
  117.   for (int i = 2; i < n; i++) {
  118.     for (int j = 1; j < i; j++) {
  119.       if (i % j == 0) {
  120.         b[i] += b[j] * a[i / j];
  121.       }
  122.     }
  123.     b[i] /= -log(i);
  124.   }
  125.   return b;
  126. }
  127.  
  128. vector<double> dexp2(vector<double> a) {
  129.   const int n = a.size();
  130.   vector<double> f(30);
  131.   f[0] = 1;
  132.   for (int i = 1; i < 30; i++) {
  133.     f[i] = f[i - 1] / i;
  134.   }
  135.   vector<double> res(n);
  136.   res[1] = 1;
  137.   for (int i = 2; i < n; i++) {
  138.     for (int j = n / i; j >= 1; j--) {
  139.       int e = 0;
  140.       double p = 1;
  141.       for (long long k = i; j * k < n; k *= i) {
  142.         e++;
  143.         p *= a[i];
  144.         res[j * k] += res[j] * p * f[e];
  145.       }
  146.     }
  147.   }
  148.   return res;
  149. }
  150.  
  151. vector<double> zeta(int n) {
  152.   vector<double> res(n);
  153.   for (int i = 1; i < n; i++) {
  154.     res[i] = 1;
  155.   }
  156.   return res;
  157. }
  158.  
  159. vector<double> mebius(int n) {
  160.   return dinv(zeta(n));
  161. }
  162.  
  163. int main() {
  164.   vector<double> x = {0,1,2,3,4,5,6,7,8,9};
  165.   vector<double> y = {0,1,1,4,1,5,9,2,6,5};
  166.   cout << dexp(x) << endl;
  167.   cout << dlog(x) << endl;
  168.   cout << dmul(x, y) << endl;
  169.   cout << dexp(dadd(dlog(x), dlog(y))) << endl;
  170.   cout << endl;
  171.   cout << dexp(x) << endl;
  172.   cout << dexp2(x) << endl;
  173. }
  174.  
  175.  
Success #stdin #stdout 0s 15240KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
[0,1,2,3,6,5,12,7,17.3333,13.5]
[0,0,2,3,2,5,-0,7,2.66667,4.5]
[0,1,3,7,7,10,26,9,20,26]
[0,1,3,7,7,10,26,9,20,26]

[0,1,2,3,6,5,12,7,17.3333,13.5]
[0,1,2,3,6,5,12,7,17.3333,13.5]
https://ideone.com/my3iM2
language:
C++14 (gcc 8.3)
created:
5 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!