71 double x0_2 = x[0] * x[0], x1_2 = x[1] * x[1], x2_2 = x[2] * x[2],
73 double x0_3 = x[0] * x[0] * x[0], x1_3 = x[1] * x[1] * x[1],
74 x2_3 = x[2] * x[2] * x[2], x3_3 = x[3] * x[3] * x[3];
76 a[0] = -(-x0_3 * x2_2 * x[3] *
f[1] + x0_3 * x[2] * x3_2 *
f[1] -
77 x0_3 *
f[3] * x[2] * x1_2 + x0_3 * x[3] *
f[2] * x1_2 +
78 x0_3 *
f[3] * x2_2 * x[1] - x0_3 * x3_2 *
f[2] * x[1] +
79 x0_2 * x1_3 *
f[3] * x[2] - x0_2 * x1_3 *
f[2] * x[3] +
80 x0_2 * x3_3 *
f[2] * x[1] + x0_2 *
f[1] * x2_3 * x[3] -
81 x0_2 *
f[1] * x3_3 * x[2] - x0_2 *
f[3] * x2_3 * x[1] +
82 x[0] * x3_2 *
f[2] * x1_3 - x[0] *
f[3] * x2_2 * x1_3 +
83 x[0] * x1_2 *
f[3] * x2_3 - x[0] * x1_2 *
f[2] * x3_3 -
84 x[0] *
f[1] * x3_2 * x2_3 + x[0] *
f[1] * x2_2 * x3_3 -
85 f[0] * x2_3 * x1_2 * x[3] +
f[0] * x2_2 * x1_3 * x[3] +
86 f[0] * x3_2 * x2_3 * x[1] -
f[0] * x3_3 * x2_2 * x[1] +
87 f[0] * x3_3 * x1_2 * x[2] -
f[0] * x3_2 * x1_3 * x[2]) /
88 (-x2_2 * x[0] * x3_3 + x2_2 * x[0] * x1_3 - x0_2 * x[3] * x2_3 +
89 x0_2 * x3_3 * x[2] - x0_2 * x[1] * x3_3 + x0_2 * x[1] * x2_3 +
90 x0_2 * x1_3 * x[3] - x0_2 * x1_3 * x[2] + x3_2 * x[0] * x2_3 -
91 x3_2 * x[0] * x1_3 + x[3] * x2_3 * x1_2 - x3_2 * x2_3 * x[1] +
92 x3_3 * x2_2 * x[1] - x3_3 * x[2] * x1_2 + x[0] * x3_3 * x1_2 -
93 x[0] * x2_3 * x1_2 - x0_3 * x3_2 * x[2] - x0_3 * x[3] * x1_2 +
94 x0_3 * x3_2 * x[1] + x[2] * x3_2 * x1_3 - x2_2 * x[3] * x1_3 +
95 x0_3 * x2_2 * x[3] - x0_3 * x2_2 * x[1] + x0_3 * x[2] * x1_2);
96 a[1] = (-x2_3 * x1_2 *
f[0] + x3_2 * x2_3 *
f[0] + x2_3 * x0_2 *
f[1] +
97 x1_2 *
f[3] * x2_3 - x2_3 * x0_2 *
f[3] -
f[1] * x3_2 * x2_3 -
98 f[3] * x2_2 * x1_3 - x3_3 * x2_2 *
f[0] +
f[1] * x2_2 * x3_3 +
99 x2_2 * x1_3 *
f[0] -
f[1] * x2_2 * x0_3 +
f[3] * x2_2 * x0_3 -
100 x1_3 * x3_2 *
f[0] - x0_2 * x1_3 *
f[2] -
f[3] * x0_3 * x1_2 +
101 f[1] * x3_2 * x0_3 + x1_2 *
f[2] * x0_3 + x3_3 *
f[0] * x1_2 -
102 x3_2 *
f[2] * x0_3 -
f[1] * x3_3 * x0_2 + x0_2 * x3_3 *
f[2] -
103 x1_2 *
f[2] * x3_3 + x3_2 *
f[2] * x1_3 + x0_2 * x1_3 *
f[3]) /
105 (-x2_2 * x0_2 * x[3] - x2_2 * x[1] * x3_2 + x2_2 * x1_2 * x[3] +
106 x2_2 * x[0] * x3_2 - x2_2 * x[0] * x1_2 + x2_2 * x0_2 * x[1] +
107 x[2] * x[0] * x1_3 + x[2] * x0_3 * x[3] - x[2] * x0_3 * x[1] -
108 x[2] * x1_3 * x[3] - x[2] * x0_2 * x3_2 + x[2] * x1_2 * x3_2 -
109 x[2] * x[3] * x[0] * x1_2 + x[2] * x[3] * x0_2 * x[1] +
110 x0_3 * x1_2 - x0_2 * x1_3 + x[3] * x[0] * x1_3 -
111 x[3] * x0_3 * x[1] - x3_2 * x[0] * x1_2 + x3_2 * x0_2 * x[1]);
112 a[2] = -(-x1_3 *
f[3] * x[2] + x1_3 *
f[2] * x[3] + x1_3 * x[0] *
f[3] +
113 x1_3 *
f[0] * x[2] - x1_3 * x[0] *
f[2] - x1_3 *
f[0] * x[3] +
114 f[3] * x2_3 * x[1] -
f[3] * x0_3 * x[1] - x[1] * x2_3 *
f[0] +
115 x[1] *
f[2] * x0_3 + x3_3 *
f[0] * x[1] - x3_3 *
f[2] * x[1] +
116 f[1] * x[3] * x0_3 -
f[1] * x[0] * x3_3 - x3_3 *
f[0] * x[2] +
117 x3_3 * x[0] *
f[2] -
f[2] * x0_3 * x[3] + x2_3 *
f[0] * x[3] +
118 f[1] * x3_3 * x[2] -
f[1] * x2_3 * x[3] + x2_3 * x[0] *
f[1] -
119 x2_3 * x[0] *
f[3] +
f[3] * x0_3 * x[2] - x[2] *
f[1] * x0_3) /
120 (x[3] * x2_2 - x3_2 * x[2] + x[1] * x3_2 - x[1] * x2_2 -
121 x1_2 * x[3] + x1_2 * x[2]) /
122 (-x[2] * x[1] * x[3] + x[1] * x[2] * x[0] - x0_2 * x[1] +
123 x[1] * x[3] * x[0] + x[2] * x[3] * x[0] + x0_3 - x0_2 * x[2] -
125 a[3] = (x[0] *
f[3] * x1_2 - x0_2 * x[3] *
f[2] + x2_2 * x[0] *
f[1] +
126 x0_2 *
f[3] * x[2] - x3_2 *
f[2] * x[1] -
f[0] * x3_2 * x[2] -
127 f[3] * x[2] * x1_2 - x2_2 * x[0] *
f[3] -
f[0] * x2_2 * x[1] +
128 f[3] * x2_2 * x[1] + x0_2 *
f[1] * x[3] + x[2] * x3_2 *
f[1] -
129 x0_2 *
f[1] * x[2] +
f[0] * x[2] * x1_2 + x[3] *
f[2] * x1_2 +
130 f[0] * x3_2 * x[1] + x3_2 * x[0] *
f[2] - x[0] *
f[2] * x1_2 -
131 f[0] * x[3] * x1_2 - x0_2 * x[1] *
f[3] + x0_2 * x[1] *
f[2] +
132 f[0] * x2_2 * x[3] - x2_2 * x[3] *
f[1] - x3_2 * x[0] *
f[1]) /
133 (-x2_2 * x[0] * x3_3 + x2_2 * x[0] * x1_3 - x0_2 * x[3] * x2_3 +
134 x0_2 * x3_3 * x[2] - x0_2 * x[1] * x3_3 + x0_2 * x[1] * x2_3 +
135 x0_2 * x1_3 * x[3] - x0_2 * x1_3 * x[2] + x3_2 * x[0] * x2_3 -
136 x3_2 * x[0] * x1_3 + x[3] * x2_3 * x1_2 - x3_2 * x2_3 * x[1] +
137 x3_3 * x2_2 * x[1] - x3_3 * x[2] * x1_2 + x[0] * x3_3 * x1_2 -
138 x[0] * x2_3 * x1_2 - x0_3 * x3_2 * x[2] - x0_3 * x[3] * x1_2 +
139 x0_3 * x3_2 * x[1] + x[2] * x3_2 * x1_3 - x2_2 * x[3] * x1_3 +
140 x0_3 * x2_2 * x[3] - x0_3 * x2_2 * x[1] + x0_3 * x[2] * x1_2);
145 const int npts_per_task = std::numeric_limits<int>::max(),
146 World* world_ptr =
nullptr) {
147 const bool use_threads = npts_per_task <
npt && world_ptr;
150 const auto iend =
npt - 2;
151 for (
int i = 1; i < iend; i += npts_per_task) {
152 auto task = [istart = i,
this, &x, &
p, npts_per_task, iend]() {
153 const auto ifence = std::min(istart + npts_per_task, iend);
154 for (
int i = istart; i < ifence; ++i) {
155 double mid = (x[i] + x[i + 1]) * 0.5;
156 double y[4] = {x[i - 1] - mid, x[i] - mid, x[i + 1] - mid,
158 this->a[i * 5] = mid;
163 world_ptr->taskq.add(
task);
168 world_ptr->taskq.fence();
171 for (
int j = 0; j < 5; ++j) {
173 a[5 *
npt - 5 + j] =
a[5 *
npt - 10 + j] =
a[5 *
npt - 15 + j];