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
111
112
113
114
115
116
117
118
|
#include <ctime>
#include <random>
#include "cubicEqn.H"
#include "IOstreams.H"
#include "stringList.H"
using namespace Foam;
scalar randomScalar(const scalar min, const scalar max)
{
static_assert
(
sizeof(long) == sizeof(scalar),
"Scalar and long are not the same size"
);
static std::default_random_engine generator(std::time(0));
static std::uniform_int_distribution<long>
distribution
(
std::numeric_limits<long>::min(),
std::numeric_limits<long>::max()
);
scalar x;
do
{
long i = distribution(generator);
x = reinterpret_cast<scalar&>(i);
}
while (min > mag(x) || mag(x) > max || !std::isfinite(x));
return x;
};
template <class Type>
void test(const Type& polynomialEqn, const scalar tol)
{
Roots<Type::nComponents - 1> r = polynomialEqn.roots();
const scalar nan = std::numeric_limits<scalar>::quiet_NaN();
const scalar inf = std::numeric_limits<scalar>::infinity();
FixedList<label, Type::nComponents - 1> t;
FixedList<scalar, Type::nComponents - 1> v(nan);
FixedList<scalar, Type::nComponents - 1> e(nan);
bool ok = true;
forAll(r, i)
{
t[i] = r.type(i);
switch (t[i])
{
case roots::real:
v[i] = polynomialEqn.value(r[i]);
e[i] = polynomialEqn.error(r[i]);
ok = ok && mag(v[i]) <= tol*mag(e[i]);
break;
case roots::posInf:
v[i] = + inf;
e[i] = nan;
break;
case roots::negInf:
v[i] = - inf;
e[i] = nan;
break;
default:
v[i] = e[i] = nan;
break;
}
}
if (!ok)
{
Info<< "Coeffs: " << polynomialEqn << endl
<< " Types: " << t << endl
<< " Roots: " << r << endl
<< "Values: " << v << endl
<< "Errors: " << e << endl << endl;
}
}
int main()
{
const scalar tol = 5;
const label nTests = 1000000;
for (label t = 0; t < nTests; ++ t)
{
test
(
cubicEqn
(
randomScalar(1e-50, 1e+50),
randomScalar(1e-50, 1e+50),
randomScalar(1e-50, 1e+50),
randomScalar(1e-50, 1e+50)
),
tol
);
}
Info << nTests << " random cubics tested" << endl;
const label coeffMin = -9, coeffMax = 10, nCoeff = coeffMax - coeffMin;
for (label a = coeffMin; a < coeffMax; ++ a)
{
for (label b = coeffMin; b < coeffMax; ++ b)
{
for (label c = coeffMin; c < coeffMax; ++ c)
{
for (label d = coeffMin; d < coeffMax; ++ d)
{
test(cubicEqn(a, b, c, d), tol);
}
}
}
}
Info<< nCoeff*nCoeff*nCoeff*nCoeff << " integer cubics tested" << endl;
return 0;
}
|
