owens_t.c
Go to the documentation of this file.
1 /* Copyright Benjamin Sobotta 2012
2  *
3  * Use, modification and distribution are subject to the
4  * Boost Software License, Version 1.0. (See accompanying file
5  * LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt)
6  */
7 
8 /*
9  * Reference:
10  * Mike Patefield, David Tandy
11  * FAST AND ACCURATE CALCULATION OF OWEN'S T-FUNCTION
12  * Journal of Statistical Software, 5 (5), 1-25
13  */
14 #include "mconf.h"
15 
16 static const int SELECT_METHOD[] = {
17  0, 0, 1, 12, 12, 12, 12, 12, 12, 12, 12, 15, 15, 15, 8,
18  0, 1, 1, 2, 2, 4, 4, 13, 13, 14, 14, 15, 15, 15, 8,
19  1, 1, 2, 2, 2, 4, 4, 14, 14, 14, 14, 15, 15, 15, 9,
20  1, 1, 2, 4, 4, 4, 4, 6, 6, 15, 15, 15, 15, 15, 9,
21  1, 2 , 2, 4, 4, 5 , 5, 7, 7, 16 ,16, 16, 11, 11, 10,
22  1, 2 , 4, 4 , 4, 5 , 5, 7, 7, 16, 16, 16, 11, 11, 11,
23  1, 2 , 3, 3, 5, 5 , 7, 7, 16, 16, 16, 16, 16, 11, 11,
24  1, 2 , 3 , 3 , 5, 5, 17, 17, 17, 17, 16, 16, 16, 11, 11
25 };
26 
27 static const double HRANGE[] = {0.02, 0.06, 0.09, 0.125, 0.26, 0.4, 0.6, 1.6,
28  1.7, 2.33, 2.4, 3.36, 3.4, 4.8};
29 
30 static const double ARANGE[] = {0.025, 0.09, 0.15, 0.36, 0.5, 0.9, 0.99999};
31 
32 static const double ORD[] = {2, 3, 4, 5, 7, 10, 12, 18, 10, 20, 30, 0, 4, 7,
33  8, 20, 0, 0};
34 
35 static const int METHODS[] = {1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4,
36  5, 6};
37 
38 static const double C[] = {
39  0.99999999999999999999999729978162447266851932041876728736094298092917625009873,
40  -0.99999999999999999999467056379678391810626533251885323416799874878563998732905968,
41  0.99999999999999999824849349313270659391127814689133077036298754586814091034842536,
42  -0.9999999999999997703859616213643405880166422891953033591551179153879839440241685,
43  0.99999999999998394883415238173334565554173013941245103172035286759201504179038147,
44  -0.9999999999993063616095509371081203145247992197457263066869044528823599399470977,
45  0.9999999999797336340409464429599229870590160411238245275855903767652432017766116267,
46  -0.999999999574958412069046680119051639753412378037565521359444170241346845522403274,
47  0.9999999933226234193375324943920160947158239076786103108097456617750134812033362048,
48  -0.9999999188923242461073033481053037468263536806742737922476636768006622772762168467,
49  0.9999992195143483674402853783549420883055129680082932629160081128947764415749728967,
50  -0.999993935137206712830997921913316971472227199741857386575097250553105958772041501,
51  0.99996135597690552745362392866517133091672395614263398912807169603795088421057688716,
52  -0.99979556366513946026406788969630293820987757758641211293079784585126692672425362469,
53  0.999092789629617100153486251423850590051366661947344315423226082520411961968929483,
54  -0.996593837411918202119308620432614600338157335862888580671450938858935084316004769854,
55  0.98910017138386127038463510314625339359073956513420458166238478926511821146316469589567,
56  -0.970078558040693314521331982203762771512160168582494513347846407314584943870399016019,
57  0.92911438683263187495758525500033707204091967947532160289872782771388170647150321633673,
58  -0.8542058695956156057286980736842905011429254735181323743367879525470479126968822863,
59  0.73796526033030091233118357742803709382964420335559408722681794195743240930748630755,
60  -0.58523469882837394570128599003785154144164680587615878645171632791404210655891158,
61  0.415997776145676306165661663581868460503874205343014196580122174949645271353372263,
62  -0.2588210875241943574388730510317252236407805082485246378222935376279663808416534365,
63  0.1375535825163892648504646951500265585055789019410617565727090346559210218472356689,
64  -0.0607952766325955730493900985022020434830339794955745989150270485056436844239206648,
65  0.0216337683299871528059836483840390514275488679530797294557060229266785853764115,
66  -0.00593405693455186729876995814181203900550014220428843483927218267309209471516256,
67  0.0011743414818332946510474576182739210553333860106811865963485870668929503649964142,
68  -1.489155613350368934073453260689881330166342484405529981510694514036264969925132E-4,
69  9.072354320794357587710929507988814669454281514268844884841547607134260303118208E-6
70 };
71 
72 static const double PTS[] = {
73  0.35082039676451715489E-02, 0.31279042338030753740E-01,
74  0.85266826283219451090E-01, 0.16245071730812277011E+00,
75  0.25851196049125434828E+00, 0.36807553840697533536E+00,
76  0.48501092905604697475E+00, 0.60277514152618576821E+00,
77  0.71477884217753226516E+00, 0.81475510988760098605E+00,
78  0.89711029755948965867E+00, 0.95723808085944261843E+00,
79  0.99178832974629703586E+00
80 };
81 
82 static const double WTS[] = {
83  0.18831438115323502887E-01, 0.18567086243977649478E-01,
84  0.18042093461223385584E-01, 0.17263829606398753364E-01,
85  0.16243219975989856730E-01, 0.14994592034116704829E-01,
86  0.13535474469662088392E-01, 0.11886351605820165233E-01,
87  0.10070377242777431897E-01, 0.81130545742299586629E-02,
88  0.60419009528470238773E-02, 0.38862217010742057883E-02,
89  0.16793031084546090448E-02
90 };
91 
92 
93 static int get_method(double h, double a) {
94  int ihint, iaint, i;
95 
96  ihint = 14;
97  iaint = 7;
98 
99  for (i = 0; i < 14; i++) {
100  if (h <= HRANGE[i]) {
101  ihint = i;
102  break;
103  }
104  }
105 
106  for (i = 0; i < 7; i++) {
107  if (a <= ARANGE[i]) {
108  iaint = i;
109  break;
110  }
111  }
112  return SELECT_METHOD[iaint * 15 + ihint];
113 }
114 
115 
116 static double owens_t_norm1(double x) {
117  return erf(x / sqrt(2)) / 2;
118 }
119 
120 
121 static double owens_t_norm2(double x) {
122  return erfc(x / sqrt(2)) / 2;
123 }
124 
125 
126 static double owensT1(double h, double a, double m) {
127  int j = 1;
128  int jj = 1;
129 
130  double hs = -0.5 * h * h;
131  double dhs = exp(hs);
132  double as = a * a;
133  double aj = a / (2 * M_PI);
134  double dj = expm1(hs);
135  double gj = hs * dhs;
136 
137  double val = atan(a) / (2 * M_PI);
138 
139  while (1) {
140  val += dj*aj / jj;
141 
142  if (m <= j) {
143  break;
144  }
145  j++;
146  jj += 2;
147  aj *= as;
148  dj = gj - dj;
149  gj *= hs / j;
150  }
151 
152  return val;
153 }
154 
155 
156 static double owensT2(double h, double a, double ah, double m) {
157  int i = 1;
158  int maxi = 2 * m + 1;
159  double hs = h * h;
160  double as = -a * a;
161  double y = 1.0 / hs;
162  double val = 0.0;
163  double vi = a*exp(-0.5 * ah * ah) / sqrt(2 * M_PI);
164  double z = (ndtr(ah) - 0.5) / h;
165 
166  while (1) {
167  val += z;
168  if (maxi <= i) {
169  break;
170  }
171  z = y * (vi - i * z);
172  vi *= as;
173  i += 2;
174  }
175  val *= exp(-0.5 * hs) / sqrt(2 * M_PI);
176 
177  return val;
178 }
179 
180 
181 static double owensT3(double h, double a, double ah) {
182  double aa, hh, y, vi, zi, result;
183  int i;
184 
185  aa = a * a;
186  hh = h * h;
187  y = 1 / hh;
188 
189  vi = a * exp(-ah * ah/ 2) / sqrt(2 * M_PI);
190  zi = owens_t_norm1(ah) / h;
191  result = 0;
192 
193  for(i = 0; i<= 30; i++) {
194  result += zi * C[i];
195  zi = y * ((2 * i + 1) * zi - vi);
196  vi *= aa;
197  }
198 
199  result *= exp(-hh / 2) / sqrt(2 * M_PI);
200 
201  return result;
202 }
203 
204 
205 static double owensT4(double h, double a, double m) {
206  double maxi, hh, naa, ai, yi, result;
207  int i;
208 
209  maxi = 2 * m + 1;
210  hh = h * h;
211  naa = -a * a;
212 
213  i = 1;
214  ai = a * exp(-hh * (1 - naa) / 2) / (2 * M_PI);
215  yi = 1;
216  result = 0;
217 
218  while (1) {
219  result += ai * yi;
220 
221  if (maxi <= i) {
222  break;
223  }
224 
225  i += 2;
226  yi = (1 - hh * yi) / i;
227  ai *= naa;
228  }
229 
230  return result;
231 }
232 
233 
234 static double owensT5(double h, double a) {
235  double result, r, aa, nhh;
236  int i;
237 
238  result = 0;
239  r = 0;
240  aa = a * a;
241  nhh = -0.5 * h * h;
242 
243  for (i = 1; i < 14; i++) {
244  r = 1 + aa * PTS[i - 1];
245  result += WTS[i - 1] * exp(nhh * r) / r;
246  }
247 
248  result *= a;
249 
250  return result;
251 }
252 
253 
254 static double owensT6(double h, double a) {
255  double normh, y, r, result;
256 
257  normh = owens_t_norm2(h);
258  y = 1 - a;
259  r = atan2(y, (1 + a));
260  result = normh * (1 - normh) / 2;
261 
262  if (r != 0) {
263  result -= r * exp(-y * h * h / (2 * r)) / (2 * M_PI);
264  }
265 
266  return result;
267 }
268 
269 
270 static double owens_t_dispatch(double h, double a, double ah) {
271  int index, meth_code;
272  double m, result;
273 
274  if (h == 0) {
275  return atan(a) / (2 * M_PI);
276  }
277  if (a == 0) {
278  return 0;
279  }
280  if (a == 1) {
281  return owens_t_norm2(-h) * owens_t_norm2(h) / 2;
282  }
283 
284  index = get_method(h, a);
285  m = ORD[index];
286  meth_code = METHODS[index];
287 
288  switch(meth_code) {
289  case 1:
290  result = owensT1(h, a, m);
291  break;
292  case 2:
293  result = owensT2(h, a, ah, m);
294  break;
295  case 3:
296  result = owensT3(h, a, ah);
297  break;
298  case 4:
299  result = owensT4(h, a, m);
300  break;
301  case 5:
302  result = owensT5(h, a);
303  break;
304  case 6:
305  result = owensT6(h, a);
306  break;
307  default:
308  result = NAN;
309  }
310 
311  return result;
312 }
313 
314 
315 double owens_t(double h, double a) {
316  double result, fabs_a, fabs_ah, normh, normah;
317 
318  if (cephes_isnan(h) || cephes_isnan(a)) {
319  return NAN;
320  }
321 
322  /* exploit that T(-h,a) == T(h,a) */
323  h = fabs(h);
324 
325  /*
326  * Use equation (2) in the paper to remap the arguments such that
327  * h >= 0 and 0 <= a <= 1 for the call of the actual computation
328  * routine.
329  */
330  fabs_a = fabs(a);
331  fabs_ah = fabs_a * h;
332 
333  if (fabs_a == INFINITY) {
334  /* See page 13 in the paper */
335  result = 0.5 * owens_t_norm2(h);
336  }
337  else if (h == INFINITY) {
338  result = 0;
339  }
340  else if (fabs_a <= 1) {
341  result = owens_t_dispatch(h, fabs_a, fabs_ah);
342  }
343  else {
344  if (fabs_ah <= 0.67) {
345  normh = owens_t_norm1(h);
346  normah = owens_t_norm1(fabs_ah);
347  result = 0.25 - normh * normah -
348  owens_t_dispatch(fabs_ah, (1 / fabs_a), h);
349  }
350  else {
351  normh = owens_t_norm2(h);
352  normah = owens_t_norm2(fabs_ah);
353  result = (normh + normah) / 2 - normh * normah -
354  owens_t_dispatch(fabs_ah, (1 / fabs_a), h);
355  }
356  }
357 
358  if (a < 0) {
359  /* exploit that T(h,-a) == -T(h,a) */
360  return -result;
361  }
362 
363  return result;
364 }
owens_t_norm2
static double owens_t_norm2(double x)
Definition: owens_t.c:121
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
WTS
static const double WTS[]
Definition: owens_t.c:82
owensT3
static double owensT3(double h, double a, double ah)
Definition: owens_t.c:181
cephes_isnan
#define cephes_isnan(x)
Definition: mconf.h:99
atan
const EIGEN_DEVICE_FUNC AtanReturnType atan() const
Definition: ArrayCwiseUnaryOps.h:283
PTS
static const double PTS[]
Definition: owens_t.c:72
C
static const double C[]
Definition: owens_t.c:38
METHODS
static const int METHODS[]
Definition: owens_t.c:35
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
owensT5
static double owensT5(double h, double a)
Definition: owens_t.c:234
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
h
const double h
Definition: testSimpleHelicopter.cpp:19
result
Values result
Definition: OdometryOptimize.cpp:8
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
owensT2
static double owensT2(double h, double a, double ah, double m)
Definition: owens_t.c:156
expm1
double expm1(double x)
Definition: unity.c:106
owensT4
static double owensT4(double h, double a, double m)
Definition: owens_t.c:205
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
owens_t_norm1
static double owens_t_norm1(double x)
Definition: owens_t.c:116
owensT1
static double owensT1(double h, double a, double m)
Definition: owens_t.c:126
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
ARANGE
static const double ARANGE[]
Definition: owens_t.c:30
owens_t_dispatch
static double owens_t_dispatch(double h, double a, double ah)
Definition: owens_t.c:270
ndtr
double ndtr(double a)
Definition: ndtr.c:201
y
Scalar * y
Definition: level1_cplx_impl.h:124
atan2
AnnoyingScalar atan2(const AnnoyingScalar &y, const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:110
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
SELECT_METHOD
static const int SELECT_METHOD[]
Definition: owens_t.c:16
mconf.h
erf
double erf(double x)
Definition: ndtr.c:285
internal
Definition: BandTriangularSolver.h:13
Eigen::numext::maxi
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
Definition: Eigen/src/Core/MathFunctions.h:1093
M_PI
#define M_PI
Definition: mconf.h:117
erfc
double erfc(double a)
Definition: ndtr.c:227
HRANGE
static const double HRANGE[]
Definition: owens_t.c:27
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
owens_t
double owens_t(double h, double a)
Definition: owens_t.c:315
owensT6
static double owensT6(double h, double a)
Definition: owens_t.c:254
get_method
static int get_method(double h, double a)
Definition: owens_t.c:93
ORD
static const double ORD[]
Definition: owens_t.c:32


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:03:14