#include <iostream>
#include <vector>
#include <cmath>
using namespace std;
ostream &operator<<(ostream &o, vector<double> a) {
o << '[';
for (size_t i = 0; i < a.size(); i++) {
if (i > 0) o << ',';
o << a[i];
}
o << ']';
return o;
}
vector<double> dadd(vector<double> a, vector<double> b) {
const int n = a.size();
for (int i = 1; i < n; i++) {
a[i] += b[i];
}
return a;
}
vector<double> dsub(vector<double> a, vector<double> b) {
const int n = a.size();
for (int i = 1; i < n; i++) {
a[i] -= b[i];
}
return a;
}
vector<double> dmul(vector<double> a, vector<double> b) {
const int n = a.size();
vector<double> c(n);
for (int i = 1; i < n; i++) {
for (int j = 1; i * j < n; j++) {
c[i * j] += a[i] * b[j];
}
}
return c;
}
vector<double> ddiff(vector<double> a) {
const int n = a.size();
a[1] = 0;
for (int i = 2; i < n; i++) {
a[i] *= -log(i);
}
return a;
}
vector<double> dint(vector<double> a) {
const int n = a.size();
for (int i = 2; i < n; i++) {
a[i] /= -log(i);
}
return a;
}
/*
a[1]b[1] = 1
a[1]b[2] + a[2]b[1] = 0
a[1]b[3] + a[3]b[1] = 0
a[1]b[4] + a[2]b[2] + a[4]b[1] = 0
a[1]b[5] + a[5]b[1] = 0
a[1]b[6] + a[2]b[3] + a[3]b[2] + a[6]b[1] = 0
*/
vector<double> dinv(vector<double> a) {
const int n = a.size();
vector<double> b(n);
b[1] = 1 / a[1];
for (int i = 2; i < n; i++) {
for (int j = 1; j < i; j++) {
if (i % j == 0) {
b[i] -= b[j] * a[i / j];
}
}
b[i] /= a[1];
}
return b;
}
/*
g(x) = log(f(x))
g'(x) = f'(x)/f(x)
*/
vector<double> dlog(vector<double> a) {
vector<double> b = dint(dmul(ddiff(a), dinv(a)));
return b;
}
/*
g(x) = exp(f(x))
g'(x) = f'(x)g(x)
g[1] = 1
-log(1)g[1] = f'[1]g[1]
-log(2)g[2] = f'[1]g[2] + f'[2]g[1]
-log(3)g[3] = f'[1]g[3] + f'[3]g[1]
-log(4)g[4] = f'[1]g[4] + f'[2]g[2] + f'[4]g[1]
-log(5)g[5] = f'[1]g[5] + f'[5]g[1]
-log(6)g[6] = f'[1]g[6] + f'[2]g[3] + f'[3]g[2] + f'[6]g[1]
-log(1)g[1] = -log(1)f[1]g[1]
-log(2)g[2] = -log(1)f[1]g[2] - log(2)f[2]g[1]
-log(3)g[3] = -log(1)f[1]g[3] - log(3)f[3]g[1]
-log(4)g[4] = -log(1)f[1]g[4] - log(2)f[2]g[2] - log(4)f[4]g[1]
-log(5)g[5] = -log(1)f[1]g[5] - log(5)f[5]g[1]
-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]
*/
vector<double> dexp(vector<double> a) {
const int n = a.size();
vector<double> b(n);
b[1] = 1;
a = ddiff(a);
for (int i = 2; i < n; i++) {
for (int j = 1; j < i; j++) {
if (i % j == 0) {
b[i] += b[j] * a[i / j];
}
}
b[i] /= -log(i);
}
return b;
}
vector<double> dexp2(vector<double> a) {
const int n = a.size();
vector<double> f(30);
f[0] = 1;
for (int i = 1; i < 30; i++) {
f[i] = f[i - 1] / i;
}
vector<double> res(n);
res[1] = 1;
for (int i = 2; i < n; i++) {
for (int j = n / i; j >= 1; j--) {
int e = 0;
double p = 1;
for (long long k = i; j * k < n; k *= i) {
e++;
p *= a[i];
res[j * k] += res[j] * p * f[e];
}
}
}
return res;
}
vector<double> zeta(int n) {
vector<double> res(n);
for (int i = 1; i < n; i++) {
res[i] = 1;
}
return res;
}
vector<double> mebius(int n) {
return dinv(zeta(n));
}
int main() {
vector<double> x = {0,1,2,3,4,5,6,7,8,9};
vector<double> y = {0,1,1,4,1,5,9,2,6,5};
cout << dexp(x) << endl;
cout << dlog(x) << endl;
cout << dmul(x, y) << endl;
cout << dexp(dadd(dlog(x), dlog(y))) << endl;
cout << endl;
cout << dexp(x) << endl;
cout << dexp2(x) << endl;
}