1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
|
#include <getopt.h>
#include <vector>
#include <cmath>
#include <iostream>
#include <fftw3.h>
#include <complex>
#include <fstream>
using Real = double;
using Complex = std::complex<Real>;
using namespace std::complex_literals;
inline Real fP(unsigned p, Real q) {
return 0.5 * pow(q, p);
}
inline Real dfP(unsigned p, Real q) {
return 0.5 * p * pow(q, p - 1);
}
inline Real ddfP(unsigned p, Real q) {
return 0.5 * p * (p - 1) * pow(q, p - 2);
}
inline Real f(Real λ, unsigned p, unsigned s, Real q) {
return (1 - λ) * fP(p, q) + λ * fP(s, q);
}
inline Real df(Real λ, unsigned p, unsigned s, Real q) {
return (1 - λ) * dfP(p, q) + λ * dfP(s, q);
}
inline Real ddf(Real λ, unsigned p, unsigned s, Real q) {
return (1 - λ) * ddfP(p, q) + λ * ddfP(s, q);
}
int main(int argc, char* argv[]) {
unsigned p = 2;
unsigned s = 2;
Real λ = 0.5;
Real τ₀ = 0;
Real y₀ = 0;
Real yₘₐₓ = 0.5;
Real Δy = 0.05;
unsigned log2n = 8;
Real τₘₐₓ = 20;
int opt;
while ((opt = getopt(argc, argv, "p:s:2:T:t:0:y:d:")) != -1) {
switch (opt) {
case 'p':
p = atoi(optarg);
break;
case 's':
s = atoi(optarg);
break;
case '2':
log2n = atoi(optarg);
break;
case 'T':
τₘₐₓ = atof(optarg);
break;
case 't':
τ₀ = atof(optarg);
break;
case '0':
y₀ = atof(optarg);
break;
case 'y':
yₘₐₓ = atof(optarg);
break;
case 'd':
Δy = atof(optarg);
break;
default:
exit(1);
}
}
unsigned n = pow(2, log2n);
Real Δτ = τₘₐₓ / M_PI / n;
Real y = y₀;
while (y += Δy, y <= yₘₐₓ) {
std::vector<Real> C(2 * n);
std::vector<Real> R(2 * n);
std::string file_end = std::to_string(p) + "_" + std::to_string(s) + "_" + std::to_string(λ) + "_" + std::to_string(τ₀) + "_" + std::to_string(y) + "_" + std::to_string(log2n) + "_" + std::to_string(τₘₐₓ) + ".dat";
std::ifstream cfile("C_"+file_end, std::ios::binary);
cfile.read((char*)(C.data()), C.size() * sizeof(Real));
cfile.close();
std::ifstream rfile("R_"+file_end, std::ios::binary);
rfile.read((char*)(R.data()), R.size() * sizeof(Real));
rfile.close();
Real energy = 0;
for (unsigned i = 0; i < n; i++) {
energy += y * R[i] * df(λ, p, s, C[i]) * M_PI * Δτ;
}
std::cerr << y << " " << energy << std::endl;
}
return 0;
}
|
