1 // 2007-02-04 Edward Smith-Rowland <3dw4rd@verizon.net>
2 //
3 // Copyright (C) 2007-2016 Free Software Foundation, Inc.
4 //
5 // This file is part of the GNU ISO C++ Library. This library is free
6 // software; you can redistribute it and/or modify it under the
7 // terms of the GNU General Public License as published by the
8 // Free Software Foundation; either version 3, or (at your option)
9 // any later version.
10 //
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License along
17 // with this library; see the file COPYING3. If not see
18 // <http://www.gnu.org/licenses/>.
19
20 // sph_neumann
21
22
23 // Compare against values generated by the GNU Scientific Library.
24 // The GSL can be found on the web: http://www.gnu.org/software/gsl/
25
26 #include <tr1/cmath>
27 #if defined(__TEST_DEBUG)
28 #include <iostream>
29 #define VERIFY(A) \
30 if (!(A)) \
31 { \
32 std::cout << "line " << __LINE__ \
33 << " max_abs_frac = " << max_abs_frac \
34 << std::endl; \
35 }
36 #else
37 #include <testsuite_hooks.h>
38 #endif
39 #include "../testcase.h"
40
41
42 // Test data for n=0.
43 testcase_sph_neumann<double> data001[] = {
44 { -3.8756496868425789, 0, 0.25000000000000000 },
45 { -1.7551651237807455, 0, 0.50000000000000000 },
46 { -0.97558515849842786, 0, 0.75000000000000000 },
47 { -0.54030230586813977, 0, 1.0000000000000000 },
48 { -0.25225788991621495, 0, 1.2500000000000000 },
49 { -0.047158134445135273, 0, 1.5000000000000000 },
50 { 0.10185488894256690, 0, 1.7500000000000000 },
51 { 0.20807341827357120, 0, 2.0000000000000000 },
52 { 0.27918827676566177, 0, 2.2500000000000000 },
53 { 0.32045744621877348, 0, 2.5000000000000000 },
54 { 0.33610995586635040, 0, 2.7500000000000000 },
55 { 0.32999749886681512, 0, 3.0000000000000000 },
56 { 0.30588605417862963, 0, 3.2500000000000000 },
57 { 0.26755905351165610, 0, 3.5000000000000000 },
58 { 0.21881582862388288, 0, 3.7500000000000000 },
59 { 0.16341090521590299, 0, 4.0000000000000000 },
60 { 0.10496176233265714, 0, 4.2500000000000000 },
61 { 0.046843510984617719, 0, 4.5000000000000000 },
62 { -0.0079162427132582220, 0, 4.7500000000000000 },
63 { -0.056732437092645263, 0, 5.0000000000000000 },
64 };
65
66 // Test function for n=0.
67 template<typename Tp>
68 void
test001()69 test001()
70 {
71 bool test [[gnu::unused]] = true;
72 const Tp eps = std::numeric_limits<Tp>::epsilon();
73 Tp max_abs_diff = -Tp(1);
74 Tp max_abs_frac = -Tp(1);
75 unsigned int num_datum = sizeof(data001)
76 / sizeof(testcase_sph_neumann<double>);
77 for (unsigned int i = 0; i < num_datum; ++i)
78 {
79 const Tp f = std::tr1::sph_neumann(Tp(data001[i].n), Tp(data001[i].x));
80 const Tp f0 = data001[i].f0;
81 const Tp diff = f - f0;
82 if (std::abs(diff) > max_abs_diff)
83 max_abs_diff = std::abs(diff);
84 if (std::abs(f0) > Tp(10) * eps
85 && std::abs(f) > Tp(10) * eps)
86 {
87 const Tp frac = diff / f0;
88 if (std::abs(frac) > max_abs_frac)
89 max_abs_frac = std::abs(frac);
90 }
91 }
92 VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
93 }
94
95 // Test data for n=1.
96 testcase_sph_neumann<double> data002[] = {
97 { -16.492214584388407, 1, 0.25000000000000000 },
98 { -4.4691813247698970, 1, 0.50000000000000000 },
99 { -2.2096318913623492, 1, 0.75000000000000000 },
100 { -1.3817732906760363, 1, 1.0000000000000000 },
101 { -0.96099400741744090, 1, 1.2500000000000000 },
102 { -0.69643541403279308, 1, 1.5000000000000000 },
103 { -0.50407489024649721, 1, 1.7500000000000000 },
104 { -0.35061200427605527, 1, 2.0000000000000000 },
105 { -0.22172663116544869, 1, 2.2500000000000000 },
106 { -0.11120587915407318, 1, 2.5000000000000000 },
107 { -0.016564013158538646, 1, 2.7500000000000000 },
108 { 0.062959163602315973, 1, 3.0000000000000000 },
109 { 0.12740959652576553, 1, 3.2500000000000000 },
110 { 0.17666922320036457, 1, 3.5000000000000000 },
111 { 0.21076723929766045, 1, 3.7500000000000000 },
112 { 0.23005335013095779, 1, 4.0000000000000000 },
113 { 0.23528261660264485, 1, 4.2500000000000000 },
114 { 0.22763858414438104, 1, 4.5000000000000000 },
115 { 0.20871085184465679, 1, 4.7500000000000000 },
116 { 0.18043836751409864, 1, 5.0000000000000000 },
117 };
118
119 // Test function for n=1.
120 template<typename Tp>
121 void
test002()122 test002()
123 {
124 bool test [[gnu::unused]] = true;
125 const Tp eps = std::numeric_limits<Tp>::epsilon();
126 Tp max_abs_diff = -Tp(1);
127 Tp max_abs_frac = -Tp(1);
128 unsigned int num_datum = sizeof(data002)
129 / sizeof(testcase_sph_neumann<double>);
130 for (unsigned int i = 0; i < num_datum; ++i)
131 {
132 const Tp f = std::tr1::sph_neumann(Tp(data002[i].n), Tp(data002[i].x));
133 const Tp f0 = data002[i].f0;
134 const Tp diff = f - f0;
135 if (std::abs(diff) > max_abs_diff)
136 max_abs_diff = std::abs(diff);
137 if (std::abs(f0) > Tp(10) * eps
138 && std::abs(f) > Tp(10) * eps)
139 {
140 const Tp frac = diff / f0;
141 if (std::abs(frac) > max_abs_frac)
142 max_abs_frac = std::abs(frac);
143 }
144 }
145 VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
146 }
147
148 // Test data for n=2.
149 testcase_sph_neumann<double> data003[] = {
150 { -194.03092532581832, 2, 0.25000000000000000 },
151 { -25.059922824838637, 2, 0.50000000000000000 },
152 { -7.8629424069509692, 2, 0.75000000000000000 },
153 { -3.6050175661599688, 2, 1.0000000000000000 },
154 { -2.0541277278856431, 2, 1.2500000000000000 },
155 { -1.3457126936204509, 2, 1.5000000000000000 },
156 { -0.96598327222227631, 2, 1.7500000000000000 },
157 { -0.73399142468765399, 2, 2.0000000000000000 },
158 { -0.57482378498626008, 2, 2.2500000000000000 },
159 { -0.45390450120366133, 2, 2.5000000000000000 },
160 { -0.35417978840293796, 2, 2.7500000000000000 },
161 { -0.26703833526449916, 2, 3.0000000000000000 },
162 { -0.18827719584715374, 2, 3.2500000000000000 },
163 { -0.11612829076848646, 2, 3.5000000000000000 },
164 { -0.050202037185754500, 2, 3.7500000000000000 },
165 { 0.0091291073823153435, 2, 4.0000000000000000 },
166 { 0.061120084680974532, 2, 4.2500000000000000 },
167 { 0.10491554511163632, 2, 4.5000000000000000 },
168 { 0.13973362282567303, 2, 4.7500000000000000 },
169 { 0.16499545760110443, 2, 5.0000000000000000 },
170 };
171
172 // Test function for n=2.
173 template<typename Tp>
174 void
test003()175 test003()
176 {
177 bool test [[gnu::unused]] = true;
178 const Tp eps = std::numeric_limits<Tp>::epsilon();
179 Tp max_abs_diff = -Tp(1);
180 Tp max_abs_frac = -Tp(1);
181 unsigned int num_datum = sizeof(data003)
182 / sizeof(testcase_sph_neumann<double>);
183 for (unsigned int i = 0; i < num_datum; ++i)
184 {
185 const Tp f = std::tr1::sph_neumann(Tp(data003[i].n), Tp(data003[i].x));
186 const Tp f0 = data003[i].f0;
187 const Tp diff = f - f0;
188 if (std::abs(diff) > max_abs_diff)
189 max_abs_diff = std::abs(diff);
190 if (std::abs(f0) > Tp(10) * eps
191 && std::abs(f) > Tp(10) * eps)
192 {
193 const Tp frac = diff / f0;
194 if (std::abs(frac) > max_abs_frac)
195 max_abs_frac = std::abs(frac);
196 }
197 }
198 VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
199 }
200
201 // Test data for n=5.
202 testcase_sph_neumann<double> data004[] = {
203 { -3884190.0626637731, 5, 0.25000000000000000 },
204 { -61327.563166980639, 5, 0.50000000000000000 },
205 { -5478.9529323190836, 5, 0.75000000000000000 },
206 { -999.44034339223640, 5, 1.0000000000000000 },
207 { -270.49720502942358, 5, 1.2500000000000000 },
208 { -94.236110085232468, 5, 1.5000000000000000 },
209 { -39.182827786584333, 5, 1.7500000000000000 },
210 { -18.591445311190984, 5, 2.0000000000000000 },
211 { -9.7821420203182274, 5, 2.2500000000000000 },
212 { -5.5991001548063233, 5, 2.5000000000000000 },
213 { -3.4400655233636823, 5, 2.7500000000000000 },
214 { -2.2470233284653904, 5, 3.0000000000000000 },
215 { -1.5491439945779160, 5, 3.2500000000000000 },
216 { -1.1205896325654248, 5, 3.5000000000000000 },
217 { -0.84592255605194844, 5, 3.7500000000000000 },
218 { -0.66280126645045878, 5, 4.0000000000000000 },
219 { -0.53589374436038528, 5, 4.2500000000000000 },
220 { -0.44430324229090551, 5, 4.5000000000000000 },
221 { -0.37520157232899892, 5, 4.7500000000000000 },
222 { -0.32046504674973919, 5, 5.0000000000000000 },
223 };
224
225 // Test function for n=5.
226 template<typename Tp>
227 void
test004()228 test004()
229 {
230 bool test [[gnu::unused]] = true;
231 const Tp eps = std::numeric_limits<Tp>::epsilon();
232 Tp max_abs_diff = -Tp(1);
233 Tp max_abs_frac = -Tp(1);
234 unsigned int num_datum = sizeof(data004)
235 / sizeof(testcase_sph_neumann<double>);
236 for (unsigned int i = 0; i < num_datum; ++i)
237 {
238 const Tp f = std::tr1::sph_neumann(Tp(data004[i].n), Tp(data004[i].x));
239 const Tp f0 = data004[i].f0;
240 const Tp diff = f - f0;
241 if (std::abs(diff) > max_abs_diff)
242 max_abs_diff = std::abs(diff);
243 if (std::abs(f0) > Tp(10) * eps
244 && std::abs(f) > Tp(10) * eps)
245 {
246 const Tp frac = diff / f0;
247 if (std::abs(frac) > max_abs_frac)
248 max_abs_frac = std::abs(frac);
249 }
250 }
251 VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
252 }
253
254 // Test data for n=10.
255 testcase_sph_neumann<double> data005[] = {
256 { -2750653598174213.5, 10, 0.25000000000000000 },
257 { -1349739281107.0554, 10, 0.50000000000000000 },
258 { -15733380424.953760, 10, 0.75000000000000000 },
259 { -672215008.25620842, 10, 1.0000000000000000 },
260 { -58607405.988679446, 10, 1.2500000000000000 },
261 { -8032728.8148234813, 10, 1.5000000000000000 },
262 { -1505955.5720640516, 10, 1.7500000000000000 },
263 { -355414.72008543846, 10, 2.0000000000000000 },
264 { -100086.80374425423, 10, 2.2500000000000000 },
265 { -32423.794085334419, 10, 2.5000000000000000 },
266 { -11772.863161809979, 10, 2.7500000000000000 },
267 { -4699.8591888113924, 10, 3.0000000000000000 },
268 { -2033.0183273853759, 10, 3.2500000000000000 },
269 { -942.19075028425493, 10, 3.5000000000000000 },
270 { -463.65206971202474, 10, 3.7500000000000000 },
271 { -240.53552987988931, 10, 4.0000000000000000 },
272 { -130.78478404631085, 10, 4.2500000000000000 },
273 { -74.170665501737531, 10, 4.5000000000000000 },
274 { -43.698249898184983, 10, 4.7500000000000000 },
275 { -26.656114405718711, 10, 5.0000000000000000 },
276 };
277
278 // Test function for n=10.
279 template<typename Tp>
280 void
test005()281 test005()
282 {
283 bool test [[gnu::unused]] = true;
284 const Tp eps = std::numeric_limits<Tp>::epsilon();
285 Tp max_abs_diff = -Tp(1);
286 Tp max_abs_frac = -Tp(1);
287 unsigned int num_datum = sizeof(data005)
288 / sizeof(testcase_sph_neumann<double>);
289 for (unsigned int i = 0; i < num_datum; ++i)
290 {
291 const Tp f = std::tr1::sph_neumann(Tp(data005[i].n), Tp(data005[i].x));
292 const Tp f0 = data005[i].f0;
293 const Tp diff = f - f0;
294 if (std::abs(diff) > max_abs_diff)
295 max_abs_diff = std::abs(diff);
296 if (std::abs(f0) > Tp(10) * eps
297 && std::abs(f) > Tp(10) * eps)
298 {
299 const Tp frac = diff / f0;
300 if (std::abs(frac) > max_abs_frac)
301 max_abs_frac = std::abs(frac);
302 }
303 }
304 VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
305 }
306
307 // Test data for n=20.
308 testcase_sph_neumann<double> data006[] = {
309 { -1.4077591402542251e+36, 20, 0.25000000000000000 },
310 { -6.7288761838234712e+29, 20, 0.50000000000000000 },
311 { -1.3544611382105945e+26, 20, 0.75000000000000000 },
312 { -3.2395922185789833e+23, 20, 1.0000000000000000 },
313 { -3.0096416715953060e+21, 20, 1.2500000000000000 },
314 { -6.5999646851668173e+19, 20, 1.5000000000000000 },
315 { -2.6193364753070735e+18, 20, 1.7500000000000000 },
316 { -1.6054364928152224e+17, 20, 2.0000000000000000 },
317 { -13719071872797762., 20, 2.2500000000000000 },
318 { -1524247248298953.8, 20, 2.5000000000000000 },
319 { -209484650509384.06, 20, 2.7500000000000000 },
320 { -34327545666696.488, 20, 3.0000000000000000 },
321 { -6522260876203.3174, 20, 3.2500000000000000 },
322 { -1406018871897.2307, 20, 3.5000000000000000 },
323 { -338025193731.78882, 20, 3.7500000000000000 },
324 { -89381690326.018677, 20, 4.0000000000000000 },
325 { -25701805899.474934, 20, 4.2500000000000000 },
326 { -7961859734.2407761, 20, 4.5000000000000000 },
327 { -2636237230.0850010, 20, 4.7500000000000000 },
328 { -926795140.30575466, 20, 5.0000000000000000 },
329 };
330
331 // Test function for n=20.
332 template<typename Tp>
333 void
test006()334 test006()
335 {
336 bool test [[gnu::unused]] = true;
337 const Tp eps = std::numeric_limits<Tp>::epsilon();
338 Tp max_abs_diff = -Tp(1);
339 Tp max_abs_frac = -Tp(1);
340 unsigned int num_datum = sizeof(data006)
341 / sizeof(testcase_sph_neumann<double>);
342 for (unsigned int i = 0; i < num_datum; ++i)
343 {
344 const Tp f = std::tr1::sph_neumann(Tp(data006[i].n), Tp(data006[i].x));
345 const Tp f0 = data006[i].f0;
346 const Tp diff = f - f0;
347 if (std::abs(diff) > max_abs_diff)
348 max_abs_diff = std::abs(diff);
349 if (std::abs(f0) > Tp(10) * eps
350 && std::abs(f) > Tp(10) * eps)
351 {
352 const Tp frac = diff / f0;
353 if (std::abs(frac) > max_abs_frac)
354 max_abs_frac = std::abs(frac);
355 }
356 }
357 VERIFY(max_abs_frac < Tp(2.5000000000000020e-13));
358 }
359
360 // Test data for n=50.
361 testcase_sph_neumann<double> data007[] = {
362 { -1.3823742808004061e+109, 50, 0.25000000000000000 },
363 { -6.1447912922121694e+93, 50, 0.50000000000000000 },
364 { -6.4348494908900529e+84, 50, 0.75000000000000000 },
365 { -2.7391922846297569e+78, 50, 1.0000000000000000 },
366 { -3.1365037573299931e+73, 50, 1.2500000000000000 },
367 { -2.8821098528635756e+69, 50, 1.5000000000000000 },
368 { -1.1148255024189452e+66, 50, 1.7500000000000000 },
369 { -1.2350219443670970e+63, 50, 2.0000000000000000 },
370 { -3.0565226939717125e+60, 50, 2.2500000000000000 },
371 { -1.4262702131152733e+58, 50, 2.5000000000000000 },
372 { -1.1118745474840939e+56, 50, 2.7500000000000000 },
373 { -1.3243260716629126e+54, 50, 3.0000000000000000 },
374 { -2.2519472094129334e+52, 50, 3.2500000000000000 },
375 { -5.1861507201100364e+50, 50, 3.5000000000000000 },
376 { -1.5513212909461383e+49, 50, 3.7500000000000000 },
377 { -5.8276471407899822e+47, 50, 4.0000000000000000 },
378 { -2.6745414086542661e+46, 50, 4.2500000000000000 },
379 { -1.4657308996352322e+45, 50, 4.5000000000000000 },
380 { -9.4102674366685358e+43, 50, 4.7500000000000000 },
381 { -6.9641091882698388e+42, 50, 5.0000000000000000 },
382 };
383
384 // Test function for n=50.
385 template<typename Tp>
386 void
test007()387 test007()
388 {
389 bool test [[gnu::unused]] = true;
390 const Tp eps = std::numeric_limits<Tp>::epsilon();
391 Tp max_abs_diff = -Tp(1);
392 Tp max_abs_frac = -Tp(1);
393 unsigned int num_datum = sizeof(data007)
394 / sizeof(testcase_sph_neumann<double>);
395 for (unsigned int i = 0; i < num_datum; ++i)
396 {
397 const Tp f = std::tr1::sph_neumann(Tp(data007[i].n), Tp(data007[i].x));
398 const Tp f0 = data007[i].f0;
399 const Tp diff = f - f0;
400 if (std::abs(diff) > max_abs_diff)
401 max_abs_diff = std::abs(diff);
402 if (std::abs(f0) > Tp(10) * eps
403 && std::abs(f) > Tp(10) * eps)
404 {
405 const Tp frac = diff / f0;
406 if (std::abs(frac) > max_abs_frac)
407 max_abs_frac = std::abs(frac);
408 }
409 }
410 VERIFY(max_abs_frac < Tp(5.0000000000000029e-12));
411 }
412
413 // Test data for n=100.
414 testcase_sph_neumann<double> data008[] = {
415 { -4.2856109460516407e+247, 100, 0.25000000000000000 },
416 { -1.6911720011753781e+217, 100, 0.50000000000000000 },
417 { -2.7753107402139484e+199, 100, 0.75000000000000000 },
418 { -6.6830794632586774e+186, 100, 1.0000000000000000 },
419 { -1.0906342369729277e+177, 100, 1.2500000000000000 },
420 { -1.0993184254131119e+169, 100, 1.5000000000000000 },
421 { -1.9071480498141315e+162, 100, 1.7500000000000000 },
422 { -2.6559558301924957e+156, 100, 2.0000000000000000 },
423 { -1.8154136926485787e+151, 100, 2.2500000000000000 },
424 { -4.3527631662111383e+146, 100, 2.5000000000000000 },
425 { -2.8809537014100589e+142, 100, 2.7500000000000000 },
426 { -4.4102229953033134e+138, 100, 3.0000000000000000 },
427 { -1.3651904154045514e+135, 100, 3.2500000000000000 },
428 { -7.6980749101080730e+131, 100, 3.5000000000000000 },
429 { -7.2790553499254927e+128, 100, 3.7500000000000000 },
430 { -1.0796647795893970e+126, 100, 4.0000000000000000 },
431 { -2.3785795774445298e+123, 100, 4.2500000000000000 },
432 { -7.4391596631955861e+120, 100, 4.5000000000000000 },
433 { -3.1802258278279400e+118, 100, 4.7500000000000000 },
434 { -1.7997139826259740e+116, 100, 5.0000000000000000 },
435 };
436
437 // Test function for n=100.
438 template<typename Tp>
439 void
test008()440 test008()
441 {
442 bool test [[gnu::unused]] = true;
443 const Tp eps = std::numeric_limits<Tp>::epsilon();
444 Tp max_abs_diff = -Tp(1);
445 Tp max_abs_frac = -Tp(1);
446 unsigned int num_datum = sizeof(data008)
447 / sizeof(testcase_sph_neumann<double>);
448 for (unsigned int i = 0; i < num_datum; ++i)
449 {
450 const Tp f = std::tr1::sph_neumann(Tp(data008[i].n), Tp(data008[i].x));
451 const Tp f0 = data008[i].f0;
452 const Tp diff = f - f0;
453 if (std::abs(diff) > max_abs_diff)
454 max_abs_diff = std::abs(diff);
455 if (std::abs(f0) > Tp(10) * eps
456 && std::abs(f) > Tp(10) * eps)
457 {
458 const Tp frac = diff / f0;
459 if (std::abs(frac) > max_abs_frac)
460 max_abs_frac = std::abs(frac);
461 }
462 }
463 VERIFY(max_abs_frac < Tp(5.0000000000000029e-12));
464 }
465 // sph_neumann
466
467 // Test data for n=0.
468 testcase_sph_neumann<double> data009[] = {
469 { -0.056732437092645263, 0, 5.0000000000000000 },
470 { 0.083907152907645249, 0, 10.000000000000000 },
471 { 0.050645860857254747, 0, 15.000000000000000 },
472 { -0.020404103090669597, 0, 20.000000000000000 },
473 { -0.039648112474538942, 0, 25.000000000000000 },
474 { -0.0051417149962528020, 0, 30.000000000000000 },
475 { 0.025819777288328762, 0, 35.000000000000000 },
476 { 0.016673451541306544, 0, 40.000000000000000 },
477 { -0.011673821973727327, 0, 45.000000000000000 },
478 { -0.019299320569842265, 0, 50.000000000000000 },
479 { -0.00040230465930828606, 0, 55.000000000000000 },
480 { 0.015873549673585938, 0, 60.000000000000000 },
481 { 0.0086531361728949541, 0, 65.000000000000000 },
482 { -0.0090474171869471404, 0, 70.000000000000000 },
483 { -0.012290016929663325, 0, 75.000000000000000 },
484 { 0.0013798405479880944, 0, 80.000000000000000 },
485 { 0.011580901686988727, 0, 85.000000000000000 },
486 { 0.0049785957347685574, 0, 90.000000000000000 },
487 { -0.0076860374841559963, 0, 95.000000000000000 },
488 { -0.0086231887228768404, 0, 100.00000000000000 },
489 };
490
491 // Test function for n=0.
492 template<typename Tp>
493 void
test009()494 test009()
495 {
496 bool test [[gnu::unused]] = true;
497 const Tp eps = std::numeric_limits<Tp>::epsilon();
498 Tp max_abs_diff = -Tp(1);
499 Tp max_abs_frac = -Tp(1);
500 unsigned int num_datum = sizeof(data009)
501 / sizeof(testcase_sph_neumann<double>);
502 for (unsigned int i = 0; i < num_datum; ++i)
503 {
504 const Tp f = std::tr1::sph_neumann(Tp(data009[i].n), Tp(data009[i].x));
505 const Tp f0 = data009[i].f0;
506 const Tp diff = f - f0;
507 if (std::abs(diff) > max_abs_diff)
508 max_abs_diff = std::abs(diff);
509 if (std::abs(f0) > Tp(10) * eps
510 && std::abs(f) > Tp(10) * eps)
511 {
512 const Tp frac = diff / f0;
513 if (std::abs(frac) > max_abs_frac)
514 max_abs_frac = std::abs(frac);
515 }
516 }
517 VERIFY(max_abs_frac < Tp(5.0000000000000028e-11));
518 }
519
520 // Test data for n=1.
521 testcase_sph_neumann<double> data010[] = {
522 { 0.18043836751409864, 1, 5.0000000000000000 },
523 { 0.062792826379701502, 1, 10.000000000000000 },
524 { -0.039976131953324147, 1, 15.000000000000000 },
525 { -0.046667467690914864, 1, 20.000000000000000 },
526 { 0.0037081455049293634, 1, 25.000000000000000 },
527 { 0.032762996969886965, 1, 30.000000000000000 },
528 { 0.012971498479556563, 1, 35.000000000000000 },
529 { -0.018210992723451058, 1, 40.000000000000000 },
530 { -0.019168385477952129, 1, 45.000000000000000 },
531 { 0.0048615106626817301, 1, 50.000000000000000 },
532 { 0.018170052158169303, 1, 55.000000000000000 },
533 { 0.0053447361795967109, 1, 60.000000000000000 },
534 { -0.012587316051033977, 1, 65.000000000000000 },
535 { -0.011184829982069090, 1, 70.000000000000000 },
536 { 0.0050065549130635621, 1, 75.000000000000000 },
537 { 0.012440856180892041, 1, 80.000000000000000 },
538 { 0.0022077237839479508, 1, 85.000000000000000 },
539 { -0.0098779785318421041, 1, 90.000000000000000 },
540 { -0.0072731342338976518, 1, 95.000000000000000 },
541 { 0.0049774245238688201, 1, 100.00000000000000 },
542 };
543
544 // Test function for n=1.
545 template<typename Tp>
546 void
test010()547 test010()
548 {
549 bool test [[gnu::unused]] = true;
550 const Tp eps = std::numeric_limits<Tp>::epsilon();
551 Tp max_abs_diff = -Tp(1);
552 Tp max_abs_frac = -Tp(1);
553 unsigned int num_datum = sizeof(data010)
554 / sizeof(testcase_sph_neumann<double>);
555 for (unsigned int i = 0; i < num_datum; ++i)
556 {
557 const Tp f = std::tr1::sph_neumann(Tp(data010[i].n), Tp(data010[i].x));
558 const Tp f0 = data010[i].f0;
559 const Tp diff = f - f0;
560 if (std::abs(diff) > max_abs_diff)
561 max_abs_diff = std::abs(diff);
562 if (std::abs(f0) > Tp(10) * eps
563 && std::abs(f) > Tp(10) * eps)
564 {
565 const Tp frac = diff / f0;
566 if (std::abs(frac) > max_abs_frac)
567 max_abs_frac = std::abs(frac);
568 }
569 }
570 VERIFY(max_abs_frac < Tp(2.5000000000000014e-11));
571 }
572
573 // Test data for n=2.
574 testcase_sph_neumann<double> data011[] = {
575 { 0.16499545760110443, 2, 5.0000000000000000 },
576 { -0.065069304993734783, 2, 10.000000000000000 },
577 { -0.058641087247919575, 2, 15.000000000000000 },
578 { 0.013403982937032370, 2, 20.000000000000000 },
579 { 0.040093089935130458, 2, 25.000000000000000 },
580 { 0.0084180146932414986, 2, 30.000000000000000 },
581 { -0.024707934561509628, 2, 35.000000000000000 },
582 { -0.018039275995565374, 2, 40.000000000000000 },
583 { 0.010395929608530518, 2, 45.000000000000000 },
584 { 0.019591011209603170, 2, 50.000000000000000 },
585 { 0.0013933984133902479, 2, 55.000000000000000 },
586 { -0.015606312864606101, 2, 60.000000000000000 },
587 { -0.0092340892214042153, 2, 65.000000000000000 },
588 { 0.0085680673305727519, 2, 70.000000000000000 },
589 { 0.012490279126185866, 2, 75.000000000000000 },
590 { -0.00091330844120464274, 2, 80.000000000000000 },
591 { -0.011502982024025860, 2, 85.000000000000000 },
592 { -0.0053078616858299611, 2, 90.000000000000000 },
593 { 0.0074563595609802797, 2, 95.000000000000000 },
594 { 0.0087725114585929052, 2, 100.00000000000000 },
595 };
596
597 // Test function for n=2.
598 template<typename Tp>
599 void
test011()600 test011()
601 {
602 bool test [[gnu::unused]] = true;
603 const Tp eps = std::numeric_limits<Tp>::epsilon();
604 Tp max_abs_diff = -Tp(1);
605 Tp max_abs_frac = -Tp(1);
606 unsigned int num_datum = sizeof(data011)
607 / sizeof(testcase_sph_neumann<double>);
608 for (unsigned int i = 0; i < num_datum; ++i)
609 {
610 const Tp f = std::tr1::sph_neumann(Tp(data011[i].n), Tp(data011[i].x));
611 const Tp f0 = data011[i].f0;
612 const Tp diff = f - f0;
613 if (std::abs(diff) > max_abs_diff)
614 max_abs_diff = std::abs(diff);
615 if (std::abs(f0) > Tp(10) * eps
616 && std::abs(f) > Tp(10) * eps)
617 {
618 const Tp frac = diff / f0;
619 if (std::abs(frac) > max_abs_frac)
620 max_abs_frac = std::abs(frac);
621 }
622 }
623 VERIFY(max_abs_frac < Tp(5.0000000000000028e-11));
624 }
625
626 // Test data for n=5.
627 testcase_sph_neumann<double> data012[] = {
628 { -0.32046504674973919, 5, 5.0000000000000000 },
629 { 0.093833541678691818, 5, 10.000000000000000 },
630 { 0.020475698281859061, 5, 15.000000000000000 },
631 { -0.048172347757372780, 5, 20.000000000000000 },
632 { -0.018309489232548347, 5, 25.000000000000000 },
633 { 0.026639390496569996, 5, 30.000000000000000 },
634 { 0.022006038985576210, 5, 35.000000000000000 },
635 { -0.011268975348057965, 5, 40.000000000000000 },
636 { -0.021770388372274858, 5, 45.000000000000000 },
637 { -0.00069711319645853701, 5, 50.000000000000000 },
638 { 0.017439589450220901, 5, 55.000000000000000 },
639 { 0.0088699170919343089, 5, 60.000000000000000 },
640 { -0.010421334444951861, 5, 65.000000000000000 },
641 { -0.012746769858008553, 5, 70.000000000000000 },
642 { 0.0026282888028967737, 5, 75.000000000000000 },
643 { 0.012477658581324189, 5, 80.000000000000000 },
644 { 0.0040771816818182642, 5, 85.000000000000000 },
645 { -0.0089777759570579818, 5, 90.000000000000000 },
646 { -0.0083184557896676149, 5, 95.000000000000000 },
647 { 0.0037206784862748965, 5, 100.00000000000000 },
648 };
649
650 // Test function for n=5.
651 template<typename Tp>
652 void
test012()653 test012()
654 {
655 bool test [[gnu::unused]] = true;
656 const Tp eps = std::numeric_limits<Tp>::epsilon();
657 Tp max_abs_diff = -Tp(1);
658 Tp max_abs_frac = -Tp(1);
659 unsigned int num_datum = sizeof(data012)
660 / sizeof(testcase_sph_neumann<double>);
661 for (unsigned int i = 0; i < num_datum; ++i)
662 {
663 const Tp f = std::tr1::sph_neumann(Tp(data012[i].n), Tp(data012[i].x));
664 const Tp f0 = data012[i].f0;
665 const Tp diff = f - f0;
666 if (std::abs(diff) > max_abs_diff)
667 max_abs_diff = std::abs(diff);
668 if (std::abs(f0) > Tp(10) * eps
669 && std::abs(f) > Tp(10) * eps)
670 {
671 const Tp frac = diff / f0;
672 if (std::abs(frac) > max_abs_frac)
673 max_abs_frac = std::abs(frac);
674 }
675 }
676 VERIFY(max_abs_frac < Tp(5.0000000000000028e-11));
677 }
678
679 // Test data for n=10.
680 testcase_sph_neumann<double> data013[] = {
681 { -26.656114405718711, 10, 5.0000000000000000 },
682 { -0.17245367208805784, 10, 10.000000000000000 },
683 { 0.078461689849642580, 10, 15.000000000000000 },
684 { -0.036843410496289961, 10, 20.000000000000000 },
685 { -0.021158339301097475, 10, 25.000000000000000 },
686 { 0.031219591064754939, 10, 30.000000000000000 },
687 { 0.012840593422414807, 10, 35.000000000000000 },
688 { -0.021803068636888072, 10, 40.000000000000000 },
689 { -0.014071636804469044, 10, 45.000000000000000 },
690 { 0.013524687511158758, 10, 50.000000000000000 },
691 { 0.015684932653180595, 10, 55.000000000000000 },
692 { -0.0056356895567262122, 10, 60.000000000000000 },
693 { -0.015364490270315362, 10, 65.000000000000000 },
694 { -0.0014525575672261295, 10, 70.000000000000000 },
695 { 0.012648951699549433, 10, 75.000000000000000 },
696 { 0.0068571608061120367, 10, 80.000000000000000 },
697 { -0.0080151152941401460, 10, 85.000000000000000 },
698 { -0.0098139742219019149, 10, 90.000000000000000 },
699 { 0.0025002854072314951, 10, 95.000000000000000 },
700 { 0.010025777373636155, 10, 100.00000000000000 },
701 };
702
703 // Test function for n=10.
704 template<typename Tp>
705 void
test013()706 test013()
707 {
708 bool test [[gnu::unused]] = true;
709 const Tp eps = std::numeric_limits<Tp>::epsilon();
710 Tp max_abs_diff = -Tp(1);
711 Tp max_abs_frac = -Tp(1);
712 unsigned int num_datum = sizeof(data013)
713 / sizeof(testcase_sph_neumann<double>);
714 for (unsigned int i = 0; i < num_datum; ++i)
715 {
716 const Tp f = std::tr1::sph_neumann(Tp(data013[i].n), Tp(data013[i].x));
717 const Tp f0 = data013[i].f0;
718 const Tp diff = f - f0;
719 if (std::abs(diff) > max_abs_diff)
720 max_abs_diff = std::abs(diff);
721 if (std::abs(f0) > Tp(10) * eps
722 && std::abs(f) > Tp(10) * eps)
723 {
724 const Tp frac = diff / f0;
725 if (std::abs(frac) > max_abs_frac)
726 max_abs_frac = std::abs(frac);
727 }
728 }
729 VERIFY(max_abs_frac < Tp(5.0000000000000028e-11));
730 }
731
732 // Test data for n=20.
733 testcase_sph_neumann<double> data014[] = {
734 { -926795140.30575466, 20, 5.0000000000000000 },
735 { -1211.2106053526036, 20, 10.000000000000000 },
736 { -1.5559965765652175, 20, 15.000000000000000 },
737 { -0.093401132250914398, 20, 20.000000000000000 },
738 { 0.044031985675276462, 20, 25.000000000000000 },
739 { -0.036078033606613907, 20, 30.000000000000000 },
740 { 0.029828405631319645, 20, 35.000000000000000 },
741 { -0.0048414810986760759, 20, 40.000000000000000 },
742 { -0.020504694681516944, 20, 45.000000000000000 },
743 { 0.013759531302541216, 20, 50.000000000000000 },
744 { 0.012783038861734196, 20, 55.000000000000000 },
745 { -0.013117009421906418, 20, 60.000000000000000 },
746 { -0.010338106075674407, 20, 65.000000000000000 },
747 { 0.010538610814111244, 20, 70.000000000000000 },
748 { 0.010200029094273744, 20, 75.000000000000000 },
749 { -0.0073123450945617122, 20, 80.000000000000000 },
750 { -0.010581510354950906, 20, 85.000000000000000 },
751 { 0.0036866374015298723, 20, 90.000000000000000 },
752 { 0.010498384318338270, 20, 95.000000000000000 },
753 { 5.6317293788334978e-05, 20, 100.00000000000000 },
754 };
755
756 // Test function for n=20.
757 template<typename Tp>
758 void
test014()759 test014()
760 {
761 bool test [[gnu::unused]] = true;
762 const Tp eps = std::numeric_limits<Tp>::epsilon();
763 Tp max_abs_diff = -Tp(1);
764 Tp max_abs_frac = -Tp(1);
765 unsigned int num_datum = sizeof(data014)
766 / sizeof(testcase_sph_neumann<double>);
767 for (unsigned int i = 0; i < num_datum; ++i)
768 {
769 const Tp f = std::tr1::sph_neumann(Tp(data014[i].n), Tp(data014[i].x));
770 const Tp f0 = data014[i].f0;
771 const Tp diff = f - f0;
772 if (std::abs(diff) > max_abs_diff)
773 max_abs_diff = std::abs(diff);
774 if (std::abs(f0) > Tp(10) * eps
775 && std::abs(f) > Tp(10) * eps)
776 {
777 const Tp frac = diff / f0;
778 if (std::abs(frac) > max_abs_frac)
779 max_abs_frac = std::abs(frac);
780 }
781 }
782 VERIFY(max_abs_frac < Tp(1.0000000000000007e-09));
783 }
784
785 // Test data for n=50.
786 testcase_sph_neumann<double> data015[] = {
787 { -6.9641091882698388e+42, 50, 5.0000000000000000 },
788 { -4.5282272723512023e+27, 50, 10.000000000000000 },
789 { -9.0004902645887037e+18, 50, 15.000000000000000 },
790 { -9542541667002.5117, 50, 20.000000000000000 },
791 { -363518140.71026671, 50, 25.000000000000000 },
792 { -152551.57233157745, 50, 30.000000000000000 },
793 { -386.26599186208625, 50, 35.000000000000000 },
794 { -4.3290507947291035, 50, 40.000000000000000 },
795 { -0.19968460851503758, 50, 45.000000000000000 },
796 { -0.041900001504607758, 50, 50.000000000000000 },
797 { 0.010696040672421902, 50, 55.000000000000000 },
798 { 0.0078198768555267188, 50, 60.000000000000000 },
799 { -0.010088474938191242, 50, 65.000000000000000 },
800 { 0.0062423671279824801, 50, 70.000000000000000 },
801 { 0.0011284242794941733, 50, 75.000000000000000 },
802 { -0.0093934266037485562, 50, 80.000000000000000 },
803 { 0.013108079602843424, 50, 85.000000000000000 },
804 { -0.0075396607225722626, 50, 90.000000000000000 },
805 { -0.0042605703552836558, 50, 95.000000000000000 },
806 { 0.010747822973682470, 50, 100.00000000000000 },
807 };
808
809 // Test function for n=50.
810 template<typename Tp>
811 void
test015()812 test015()
813 {
814 bool test [[gnu::unused]] = true;
815 const Tp eps = std::numeric_limits<Tp>::epsilon();
816 Tp max_abs_diff = -Tp(1);
817 Tp max_abs_frac = -Tp(1);
818 unsigned int num_datum = sizeof(data015)
819 / sizeof(testcase_sph_neumann<double>);
820 for (unsigned int i = 0; i < num_datum; ++i)
821 {
822 const Tp f = std::tr1::sph_neumann(Tp(data015[i].n), Tp(data015[i].x));
823 const Tp f0 = data015[i].f0;
824 const Tp diff = f - f0;
825 if (std::abs(diff) > max_abs_diff)
826 max_abs_diff = std::abs(diff);
827 if (std::abs(f0) > Tp(10) * eps
828 && std::abs(f) > Tp(10) * eps)
829 {
830 const Tp frac = diff / f0;
831 if (std::abs(frac) > max_abs_frac)
832 max_abs_frac = std::abs(frac);
833 }
834 }
835 VERIFY(max_abs_frac < Tp(2.5000000000000014e-11));
836 }
837
838 // Test data for n=100.
839 testcase_sph_neumann<double> data016[] = {
840 { -1.7997139826259740e+116, 100, 5.0000000000000000 },
841 { -8.5732263093296268e+85, 100, 10.000000000000000 },
842 { -1.9270658593711677e+68, 100, 15.000000000000000 },
843 { -7.2208893582952385e+55, 100, 20.000000000000000 },
844 { -2.0868752613007946e+46, 100, 25.000000000000000 },
845 { -4.2496124023612646e+38, 100, 30.000000000000000 },
846 { -1.7042898348910271e+32, 100, 35.000000000000000 },
847 { -6.3021565260724554e+26, 100, 40.000000000000000 },
848 { -1.3199917400494367e+22, 100, 45.000000000000000 },
849 { -1.1256928913265988e+18, 100, 50.000000000000000 },
850 { -309801083340343.25, 100, 55.000000000000000 },
851 { -232585620046.64737, 100, 60.000000000000000 },
852 { -421135935.93756074, 100, 65.000000000000000 },
853 { -1680637.4531202621, 100, 70.000000000000000 },
854 { -13868.302591128844, 100, 75.000000000000000 },
855 { -227.24385709173322, 100, 80.000000000000000 },
856 { -7.2807038787138731, 100, 85.000000000000000 },
857 { -0.46648154448250878, 100, 90.000000000000000 },
858 { -0.067270772720654556, 100, 95.000000000000000 },
859 { -0.022983850491562267, 100, 100.00000000000000 },
860 };
861
862 // Test function for n=100.
863 template<typename Tp>
864 void
test016()865 test016()
866 {
867 bool test [[gnu::unused]] = true;
868 const Tp eps = std::numeric_limits<Tp>::epsilon();
869 Tp max_abs_diff = -Tp(1);
870 Tp max_abs_frac = -Tp(1);
871 unsigned int num_datum = sizeof(data016)
872 / sizeof(testcase_sph_neumann<double>);
873 for (unsigned int i = 0; i < num_datum; ++i)
874 {
875 const Tp f = std::tr1::sph_neumann(Tp(data016[i].n), Tp(data016[i].x));
876 const Tp f0 = data016[i].f0;
877 const Tp diff = f - f0;
878 if (std::abs(diff) > max_abs_diff)
879 max_abs_diff = std::abs(diff);
880 if (std::abs(f0) > Tp(10) * eps
881 && std::abs(f) > Tp(10) * eps)
882 {
883 const Tp frac = diff / f0;
884 if (std::abs(frac) > max_abs_frac)
885 max_abs_frac = std::abs(frac);
886 }
887 }
888 VERIFY(max_abs_frac < Tp(2.5000000000000015e-12));
889 }
890
891 int
main()892 main()
893 {
894 test001<double>();
895 test002<double>();
896 test003<double>();
897 test004<double>();
898 test005<double>();
899 test006<double>();
900 test007<double>();
901 test008<double>();
902 test009<double>();
903 test010<double>();
904 test011<double>();
905 test012<double>();
906 test013<double>();
907 test014<double>();
908 test015<double>();
909 test016<double>();
910 return 0;
911 }
912