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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
|
// miscfp.cc -- misc floating-point functions
//
// jtutil::signum - sign of a floating point number
// jtutil::hypot3 - 3D Euclidean distance
// jtutil::arctan_xy - 4-quadrant arc tangent
// jtutil::modulo_reduce - reduce an angle modulo 2*pi radians (360 degrees)
//
// jtutil::array_zero - zero an array
// jtutil::array_copy - copy one array to another
// jtutil::array_scale - scale an array by a scalar
//
// ***** template instantiations *****
//
#include <math.h>
#include "stdc.h"
#include "util.hh"
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// This function computes the floating point "signum" function (as in APL),
// signum(x) = -1.0, 0.0, or +1.0, according to the sign of x.
//
namespace jtutil
{
double signum(double x)
{
if (x == 0.0)
then return 0.0;
else return (x > 0.0) ? 1.0 : -1.0;
}
} // namespace jtutil::
//******************************************************************************
//
// This function computes the 3-dimensional Euclidean distance,
// analogously to the standard math library function hypot(2) .
//
// Arguments:
// (x,y,z) = (in) The rectangular coordinates of a point in $\Re^3$.
//
// Result:
// The function returns the Euclidean distance of (x,y,z) from the origin.
//
// Bugs:
// Unlike hypot(2), this function takes no special care to avoid
// unwarranted IEEE exceptions if any of |x|, |y|, or |z| is close to
// the overflow and/or underflow threshold.
//
namespace jtutil
{
double hypot3(double x, double y, double z)
{
return std::sqrt(x*x + y*y + z*z);
}
} // namespace jtutil::
//******************************************************************************
//
// This function computes a four-quadrant arctangent, using what I think
// is the "right" ordering of the arguments, and returning 0.0 if both
// arguments are 0.0.
//
// Arguments:
// (x,y) = (in) The rectangular coordinates of a point in $\Re^2$.
//
// Result:
// The function returns a value $\theta$ such that $x + iy = R e^{i\theta}$ for
// some real $R$, i.e. it returns the angle between the positive $x$ axis and
// the line joining the origin and the point $(x,y)$.
//
namespace jtutil
{
double arctan_xy(double x, double y)
{
return ((x == 0.0) && (y == 0.0))
?
: std::atan2(y, x); // note reversed argument order (y,x)!!
}
} // namespace jtutil::
//******************************************************************************
//
// This function reduces x modulo xmod to be (fuzzily) in the range
// [xmin, xmax] , or does an error_exit() if no such value exists.
//
namespace jtutil
{
double modulo_reduce(double x, double xmod, double xmin, double xmax)
{
double xx = x;
while (fuzzy<double>::LT(xx, xmin))
{
xx += xmod;
}
while (fuzzy<double>::GT(xx, xmax))
{
xx -= xmod;
}
if (! (fuzzy<double>::GE(xx, xmin) && fuzzy<double>::LE(xx, xmax)) )
then jtutil::error_exit(ERROR_EXIT,
"***** modulo_reduce(): no modulo value is fuzzily within specified range!\n"
" x = %g xmod = %g\n"
" [xmin,xmax] = [%g,%g]\n"
" ==> xx = %g\n"
,
x, xmod,
xmin, xmax,
xx); /*NOTREACHED*/
return xx;
}
} // namespace jtutil::
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// This function zeros the size-N array dst[].
//
namespace jtutil
{
template <typename fp>
void array_zero(int N, fp dst[])
{
for (int i = 0 ; i < N ; ++i)
{
dst[i] = 0.0;
}
}
} // namespace jtutil::
//******************************************************************************
//
// This function copies the size-N array src[] to the size-N array dst[],
// analogously to the BLAS routines xCOPY().
//
namespace jtutil
{
template <typename fp>
void array_copy(int N, const fp src[], fp dst[])
{
for (int i = 0 ; i < N ; ++i)
{
dst[i] = src[i];
}
}
} // namespace jtutil::
//******************************************************************************
//
// This function multiplies the size-N array dst[] by the scalar alpha,
// analogously to the BLAS routines xSCAL().
//
namespace jtutil
{
template <typename fp>
void array_scale(int N, fp alpha, fp dst[])
{
for (int i = 0 ; i < N ; ++i)
{
dst[i] *= alpha;
}
}
} // namespace jtutil::
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// ***** template instantiations *****
//
namespace jtutil
{
template void array_zero<float >(int, float []);
template void array_zero<double>(int, double[]);
template void array_copy<float >(int, const float [], float []);
template void array_copy<double>(int, const double[], double[]);
template void array_scale<float >(int, float , float []);
template void array_scale<double>(int, double, double[]);
} // namespace jtutil::
|
