kolmogorov.c
Go to the documentation of this file.
1 /* File altered for inclusion in cephes module for Python:
2  * Main loop commented out.... */
3 /* Travis Oliphant Nov. 1998 */
4 
5 /* Re Kolmogorov statistics, here is Birnbaum and Tingey's (actually it was already present
6  * in Smirnov's paper) formula for the
7  * distribution of D+, the maximum of all positive deviations between a
8  * theoretical distribution function P(x) and an empirical one Sn(x)
9  * from n samples.
10  *
11  * +
12  * D = sup [P(x) - S (x)]
13  * n -inf < x < inf n
14  *
15  *
16  * [n(1-d)]
17  * + - v-1 n-v
18  * Pr{D > d} = > C d (d + v/n) (1 - d - v/n)
19  * n - n v
20  * v=0
21  *
22  * (also equals the following sum, but note the terms may be large and alternating in sign)
23  * See Smirnov 1944, Dwass 1959
24  * n
25  * - v-1 n-v
26  * = 1 - > C d (d + v/n) (1 - d - v/n)
27  * - n v
28  * v=[n(1-d)]+1
29  *
30  * [n(1-d)] is the largest integer not exceeding n(1-d).
31  * nCv is the number of combinations of n things taken v at a time.
32 
33  * Sources:
34  * [1] Smirnov, N.V. "Approximate laws of distribution of random variables from empirical data"
35  * Usp. Mat. Nauk, 1944. http://mi.mathnet.ru/umn8798
36  * [2] Birnbaum, Z. W. and Tingey, Fred H.
37  * "One-Sided Confidence Contours for Probability Distribution Functions",
38  * Ann. Math. Statist. 1951. https://doi.org/10.1214/aoms/1177729550
39  * [3] Dwass, Meyer, "The Distribution of a Generalized $\mathrm{D}^+_n$ Statistic",
40  * Ann. Math. Statist., 1959. https://doi.org/10.1214/aoms/1177706085
41  * [4] van Mulbregt, Paul, "Computing the Cumulative Distribution Function and Quantiles of the One-sided Kolmogorov-Smirnov Statistic"
42  * http://arxiv.org/abs/1802.06966
43  * [5] van Mulbregt, Paul, "Computing the Cumulative Distribution Function and Quantiles of the limit of the Two-sided Kolmogorov-Smirnov Statistic"
44  * https://arxiv.org/abs/1803.00426
45  *
46  */
47 
48 #include "mconf.h"
49 #include <float.h>
50 #include <math.h>
51 #include <assert.h>
52 
53 
54 /* ************************************************************************ */
55 /* Algorithm Configuration */
56 
57 /*
58  * Kolmogorov Two-sided:
59  * Switchover between the two series to compute K(x)
60  * 0 <= x <= KOLMOG_CUTOVER and
61  * KOLMOG_CUTOVER < x < infty
62  */
63 #define KOLMOG_CUTOVER 0.82
64 
65 
66 /*
67  * Smirnov One-sided:
68  * n larger than SMIRNOV_MAX_COMPUTE_N will result in an approximation
69  */
70 const int SMIRNOV_MAX_COMPUTE_N = 1000000;
71 
72 /*
73  * Use the upper sum formula, if the number of terms is at most SM_UPPER_MAX_TERMS,
74  * and n is at least SM_UPPERSUM_MIN_N
75  * Don't use the upper sum if lots of terms are involved as the series alternates
76  * sign and the terms get much bigger than 1.
77  */
78 #define SM_UPPER_MAX_TERMS 3
79 #define SM_UPPERSUM_MIN_N 10
80 
81 /* ************************************************************************ */
82 /* ************************************************************************ */
83 
84 /* Assuming LOW and HIGH are constants. */
85 #define CLIP(X, LOW, HIGH) ((X) < LOW ? LOW : MIN(X, HIGH))
86 #ifndef MIN
87 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
88 #endif
89 #ifndef MAX
90 #define MAX(a,b) (((a) < (b)) ? (b) : (a))
91 #endif
92 
93 /* from cephes constants */
94 extern double MINLOG;
95 
96 /* exp() of anything below this returns 0 */
97 static const int MIN_EXPABLE = (-708 - 38);
98 
99 #ifndef LOGSQRT2PI
100 #define LOGSQRT2PI 0.91893853320467274178032973640561764
101 #endif
102 
103 /* Struct to hold the CDF, SF and PDF, which are computed simultaneously */
104 typedef struct ThreeProbs {
105  double sf;
106  double cdf;
107  double pdf;
108 } ThreeProbs;
109 
110 #define RETURN_3PROBS(PSF, PCDF, PDF) \
111  ret.cdf = (PCDF); \
112  ret.sf = (PSF); \
113  ret.pdf = (PDF); \
114  return ret;
115 
116 static const double _xtol = DBL_EPSILON;
117 static const double _rtol = 2*DBL_EPSILON;
118 
119 static int
120 _within_tol(double x, double y, double atol, double rtol)
121 {
122  double diff = fabs(x-y);
123  int result = (diff <= (atol + rtol * fabs(y)));
124  return result;
125 }
126 
127 #include "dd_real.h"
128 
129 /* Shorten some of the double-double names for readibility */
130 #define valueD dd_to_double
131 #define add_dd dd_add_d_d
132 #define sub_dd dd_sub_d_d
133 #define mul_dd dd_mul_d_d
134 #define neg_D dd_neg
135 #define div_dd dd_div_d_d
136 #define add_DD dd_add
137 #define sub_DD dd_sub
138 #define mul_DD dd_mul
139 #define div_DD dd_div
140 #define add_Dd dd_add_dd_d
141 #define add_dD dd_add_d_dd
142 #define sub_Dd dd_sub_dd_d
143 #define sub_dD dd_sub_d_dd
144 #define mul_Dd dd_mul_dd_d
145 #define mul_dD dd_mul_d_dd
146 #define div_Dd dd_div_dd_d
147 #define div_dD dd_div_d_dd
148 #define frexpD dd_frexp
149 #define ldexpD dd_ldexp
150 #define logD dd_log
151 #define log1pD dd_log1p
152 
153 
154 /* ************************************************************************ */
155 /* Kolmogorov : Two-sided **************************** */
156 /* ************************************************************************ */
157 
158 static ThreeProbs
159 _kolmogorov(double x)
160 {
161  double P = 1.0;
162  double D = 0;
163  double sf, cdf, pdf;
164  ThreeProbs ret;
165 
166  if (isnan(x)) {
167  RETURN_3PROBS(NAN, NAN, NAN);
168  }
169  if (x <= 0) {
170  RETURN_3PROBS(1.0, 0.0, 0);
171  }
172  /* x <= 0.040611972203751713 */
173  if (x <= (double)M_PI/sqrt(-MIN_EXPABLE * 8)) {
174  RETURN_3PROBS(1.0, 0.0, 0);
175  }
176 
177  P = 1.0;
178  if (x <= KOLMOG_CUTOVER) {
179  /*
180  * u = e^(-pi^2/(8x^2))
181  * w = sqrt(2pi)/x
182  * P = w*u * (1 + u^8 + u^24 + u^48 + ...)
183  */
184  double w = sqrt(2 * M_PI)/x;
185  double logu8 = -M_PI * M_PI/(x * x); /* log(u^8) */
186  double u = exp(logu8/8);
187  if (u == 0) {
188  /*
189  * P = w*u, but u < 1e-308, and w > 1,
190  * so compute as logs, then exponentiate
191  */
192  double logP = logu8/8 + log(w);
193  P = exp(logP);
194  } else {
195  /* Just unroll the loop, 3 iterations */
196  double u8 = exp(logu8);
197  double u8cub = pow(u8, 3);
198  P = 1 + u8cub * P;
199  D = 5*5 + u8cub * D;
200  P = 1 + u8*u8 * P;
201  D = 3*3 + u8*u8 * D;
202  P = 1 + u8 * P;
203  D = 1*1 + u8 * D;
204 
205  D = M_PI * M_PI/4/(x*x) * D - P;
206  D *= w * u/x;
207  P = w * u * P;
208  }
209  cdf = P;
210  sf = 1-P;
211  pdf = D;
212  }
213  else {
214  /*
215  * v = e^(-2x^2)
216  * P = 2 (v - v^4 + v^9 - v^16 + ...)
217  * = 2v(1 - v^3*(1 - v^5*(1 - v^7*(1 - ...)))
218  */
219  double logv = -2*x*x;
220  double v = exp(logv);
221  /*
222  * Want q^((2k-1)^2)(1-q^(4k-1)) / q(1-q^3) < epsilon to break out of loop.
223  * With KOLMOG_CUTOVER ~ 0.82, k <= 4. Just unroll the loop, 4 iterations
224  */
225  double vsq = v*v;
226  double v3 = pow(v, 3);
227  double vpwr;
228 
229  vpwr = v3*v3*v; /* v**7 */
230  P = 1 - vpwr * P; /* P <- 1 - (1-v**(2k-1)) * P */
231  D = 3*3 - vpwr * D;
232 
233  vpwr = v3*vsq;
234  P = 1 - vpwr * P;
235  D = 2*2 - vpwr * D;
236 
237  vpwr = v3;
238  P = 1 - vpwr * P;
239  D = 1*1 - vpwr * D;
240 
241  P = 2 * v * P;
242  D = 8 * v * x * D;
243  sf = P;
244  cdf = 1 - sf;
245  pdf = D;
246  }
247  pdf = MAX(0, pdf);
248  cdf = CLIP(cdf, 0, 1);
249  sf = CLIP(sf, 0, 1);
250  RETURN_3PROBS(sf, cdf, pdf);
251 }
252 
253 
254 /* Find x such kolmogorov(x)=psf, kolmogc(x)=pcdf */
255 static double
256 _kolmogi(double psf, double pcdf)
257 {
258  double x, t;
259  double xmin = 0;
260  double xmax = INFINITY;
261  int iterations;
262  double a = xmin, b = xmax;
263 
264  if (!(psf >= 0.0 && pcdf >= 0.0 && pcdf <= 1.0 && psf <= 1.0)) {
265  sf_error("kolmogi", SF_ERROR_DOMAIN, NULL);
266  return (NAN);
267  }
268  if (fabs(1.0 - pcdf - psf) > 4* DBL_EPSILON) {
269  sf_error("kolmogi", SF_ERROR_DOMAIN, NULL);
270  return (NAN);
271  }
272  if (pcdf == 0.0) {
273  return 0.0;
274  }
275  if (psf == 0.0) {
276  return INFINITY;
277  }
278 
279  if (pcdf <= 0.5) {
280  /* p ~ (sqrt(2pi)/x) *exp(-pi^2/8x^2). Generate lower and upper bounds */
281  double logpcdf = log(pcdf);
282  const double SQRT2 = M_SQRT2;
283  /* Now that 1 >= x >= sqrt(p) */
284  /* Iterate twice: x <- pi/(sqrt(8) sqrt(log(sqrt(2pi)) - log(x) - log(pdf))) */
285  a = M_PI / (2 * SQRT2 * sqrt(-(logpcdf + logpcdf/2 - LOGSQRT2PI)));
286  b = M_PI / (2 * SQRT2 * sqrt(-(logpcdf + 0 - LOGSQRT2PI)));
287  a = M_PI / (2 * SQRT2 * sqrt(-(logpcdf + log(a) - LOGSQRT2PI)));
288  b = M_PI / (2 * SQRT2 * sqrt(-(logpcdf + log(b) - LOGSQRT2PI)));
289  x = (a + b) / 2.0;
290  }
291  else {
292  /*
293  * Based on the approximation p ~ 2 exp(-2x^2)
294  * Found that needed to replace psf with a slightly smaller number in the second element
295  * as otherwise _kolmogorov(b) came back as a very small number but with
296  * the same sign as _kolmogorov(a)
297  * kolmogi(0.5) = 0.82757355518990772
298  * so (1-q^(-(4-1)*2*x^2)) = (1-exp(-6*0.8275^2) ~ (1-exp(-4.1)
299  */
300  const double jiggerb = 256 * DBL_EPSILON;
301  double pba = psf/(1.0 - exp(-4))/2, pbb = psf * (1 - jiggerb)/2;
302  double q0;
303  a = sqrt(-0.5 * log(pba));
304  b = sqrt(-0.5 * log(pbb));
305  /*
306  * Use inversion of
307  * p = q - q^4 + q^9 - q^16 + ...:
308  * q = p + p^4 + 4p^7 - p^9 + 22p^10 - 13p^12 + 140*p^13 ...
309  */
310  {
311  double p = psf/2.0;
312  double p2 = p*p;
313  double p3 = p*p*p;
314  q0 = 1 + p3 * (1 + p3 * (4 + p2 *(-1 + p*(22 + p2* (-13 + 140 * p)))));
315  q0 *= p;
316  }
317  x = sqrt(-log(q0) / 2);
318  if (x < a || x > b) {
319  x = (a+b)/2;
320  }
321  }
322  assert(a <= b);
323 
324  iterations = 0;
325  do {
326  double x0 = x;
327  ThreeProbs probs = _kolmogorov(x0);
328  double df = ((pcdf < 0.5) ? (pcdf - probs.cdf) : (probs.sf - psf));
329  double dfdx;
330 
331  if (fabs(df) == 0) {
332  break;
333  }
334  /* Update the bracketing interval */
335  if (df > 0 && x > a) {
336  a = x;
337  } else if (df < 0 && x < b) {
338  b = x;
339  }
340 
341  dfdx = -probs.pdf;
342  if (fabs(dfdx) <= 0.0) {
343  x = (a+b)/2;
344  t = x0 - x;
345  } else {
346  t = df/dfdx;
347  x = x0 - t;
348  }
349 
350  /*
351  * Check out-of-bounds.
352  * Not expecting this to happen often --- kolmogorov is convex near x=infinity and
353  * concave near x=0, and we should be approaching from the correct side.
354  * If out-of-bounds, replace x with a midpoint of the bracket.
355  */
356  if (x >= a && x <= b) {
357  if (_within_tol(x, x0, _xtol, _rtol)) {
358  break;
359  }
360  if ((x == a) || (x == b)) {
361  x = (a + b) / 2.0;
362  /* If the bracket is already so small ... */
363  if (x == a || x == b) {
364  break;
365  }
366  }
367  } else {
368  x = (a + b) / 2.0;
369  if (_within_tol(x, x0, _xtol, _rtol)) {
370  break;
371  }
372  }
373 
374  if (++iterations > MAXITER) {
375  sf_error("kolmogi", SF_ERROR_SLOW, NULL);
376  break;
377  }
378  } while(1);
379  return (x);
380 }
381 
382 
383 double
384 kolmogorov(double x)
385 {
386  if (isnan(x)) {
387  return NAN;
388  }
389  return _kolmogorov(x).sf;
390 }
391 
392 double
393 kolmogc(double x)
394 {
395  if (isnan(x)) {
396  return NAN;
397  }
398  return _kolmogorov(x).cdf;
399 }
400 
401 double
402 kolmogp(double x)
403 {
404  if (isnan(x)) {
405  return NAN;
406  }
407  if (x <= 0) {
408  return -0.0;
409  }
410  return -_kolmogorov(x).pdf;
411 }
412 
413 /* Functional inverse of Kolmogorov survival statistic for two-sided test.
414  * Finds x such that kolmogorov(x) = p.
415  */
416 double
417 kolmogi(double p)
418 {
419  if (isnan(p)) {
420  return NAN;
421  }
422  return _kolmogi(p, 1-p);
423 }
424 
425 /* Functional inverse of Kolmogorov cumulative statistic for two-sided test.
426  * Finds x such that kolmogc(x) = p = (or kolmogorov(x) = 1-p).
427  */
428 double
429 kolmogci(double p)
430 {
431  if (isnan(p)) {
432  return NAN;
433  }
434  return _kolmogi(1-p, p);
435 }
436 
437 
438 
439 /* ************************************************************************ */
440 /* ********** Smirnov : One-sided ***************************************** */
441 /* ************************************************************************ */
442 
443 static double
445 {
446  double q = ldexp(x, 1-DBL_MANT_DIG);
447  double L = fabs(q+x);
448  if (L == 0) {
449  L = fabs(x);
450  } else {
451  int Lint = (int)(L);
452  if (Lint == L) {
453  L = Lint;
454  }
455  }
456  return L;
457 }
458 
459 static double
460 modNX(int n, double x, int *pNXFloor, double *pNX)
461 {
462  /*
463  * Compute floor(n*x) and remainder *exactly*.
464  * If remainder is too close to 1 (E.g. (1, -DBL_EPSILON/2))
465  * round up and adjust */
466  double2 alphaD, nxD, nxfloorD;
467  int nxfloor;
468  double alpha;
469 
470  nxD = mul_dd(n, x);
471  nxfloorD = dd_floor(nxD);
472  alphaD = sub_DD(nxD, nxfloorD);
473  alpha = dd_hi(alphaD);
474  nxfloor = dd_to_int(nxfloorD);
475  assert(alpha >= 0);
476  assert(alpha <= 1);
477  if (alpha == 1) {
478  nxfloor += 1;
479  alpha = 0;
480  }
481  assert(alpha < 1.0);
482  *pNX = dd_to_double(nxD);
483  *pNXFloor = nxfloor;
484  return alpha;
485 }
486 
487 /*
488  * The binomial coefficient C overflows a 64 bit double, as the 11-bit
489  * exponent is too small.
490  * Store C as (Cman:double2, Cexpt:int).
491  * I.e a Mantissa/significand, and an exponent.
492  * Cman lies between 0.5 and 1, and the exponent has >=32-bit.
493  */
494 static void
495 updateBinomial(double2 *Cman, int *Cexpt, int n, int j)
496 {
497  int expt;
498  double2 rat = div_dd(n - j, j + 1.0);
499  double2 man2 = mul_DD(*Cman, rat);
500  man2 = frexpD(man2, &expt);
501  assert (!dd_is_zero(man2));
502  *Cexpt += expt;
503  *Cman = man2;
504 }
505 
506 
507 static double2
509 {
510  /*
511  * Using dd_npwr() here would be quite time-consuming.
512  * Tradeoff accuracy-time by using pow().
513  */
514  double ans, r, adj;
515  if (m <= 0) {
516  if (m == 0) {
517  return DD_C_ONE;
518  }
519  return dd_inv(pow_D(a, -m));
520  }
521  if (dd_is_zero(a)) {
522  return DD_C_ZERO;
523  }
524  ans = pow(a.x[0], m);
525  r = a.x[1]/a.x[0];
526  adj = m*r;
527  if (fabs(adj) > 1e-8) {
528  if (fabs(adj) < 1e-4) {
529  /* Take 1st two terms of Taylor Series for (1+r)^m */
530  adj += (m*r) * ((m-1)/2.0 * r);
531  } else {
532  /* Take exp of scaled log */
533  adj = expm1(m*log1p(r));
534  }
535  }
536  return dd_add_d_d(ans, ans*adj);
537 }
538 
539 static double
540 pow2(double a, double b, int m)
541 {
542  return dd_to_double(pow_D(add_dd(a, b), m));
543 }
544 
545 /*
546  * Not 1024 as too big. Want _MAX_EXPONENT < 1023-52 so as to keep both
547  * elements of the double2 normalized
548  */
549 #define _MAX_EXPONENT 960
550 
551 #define RETURN_M_E(MAND, EXPT) \
552  *pExponent = EXPT;\
553  return MAND;
554 
555 
556 static double2
557 pow2Scaled_D(double2 a, int m, int *pExponent)
558 {
559  /* Compute a^m = significand*2^expt and return as (significand, expt) */
560  double2 ans, y;
561  int ansE, yE;
562  int maxExpt = _MAX_EXPONENT;
563  int q, r, y2mE, y2rE, y2mqE;
564  double2 y2r, y2m, y2mq;
565 
566  if (m <= 0)
567  {
568  int aE1, aE2;
569  if (m == 0) {
570  RETURN_M_E(DD_C_ONE, 0);
571  }
572  ans = pow2Scaled_D(a, -m, &aE1);
573  ans = frexpD(dd_inv(ans), &aE2);
574  ansE = -aE1 + aE2;
575  RETURN_M_E(ans, ansE);
576  }
577  y = frexpD(a, &yE);
578  if (m == 1) {
579  RETURN_M_E(y, yE);
580  }
581  /*
582  * y ^ maxExpt >= 2^{-960}
583  * => maxExpt = 960 / log2(y.x[0]) = 708 / log(y.x[0])
584  * = 665/((1-y.x[0] + y.x[0]^2/2 - ...)
585  * <= 665/(1-y.x[0])
586  * Quick check to see if we might need to break up the exponentiation
587  */
588  if (m*(y.x[0]-1) / y.x[0] < -_MAX_EXPONENT * M_LN2) {
589  /* Now do it carefully, calling log() */
590  double lg2y = log(y.x[0]) / M_LN2;
591  double lgAns = m * lg2y;
592  if (lgAns <= -_MAX_EXPONENT) {
593  maxExpt = (int)(nextPowerOf2(-_MAX_EXPONENT / lg2y + 1)/2);
594  }
595  }
596  if (m <= maxExpt)
597  {
598  double2 ans1 = pow_D(y, m);
599  ans = frexpD(ans1, &ansE);
600  ansE += m * yE;
601  RETURN_M_E(ans, ansE);
602  }
603 
604  q = m / maxExpt;
605  r = m % maxExpt;
606  /* y^m = (y^maxExpt)^q * y^r */
607  y2r = pow2Scaled_D(y, r, &y2rE);
608  y2m = pow2Scaled_D(y, maxExpt, &y2mE);
609  y2mq = pow2Scaled_D(y2m, q, &y2mqE);
610  ans = frexpD(mul_DD(y2r, y2mq), &ansE);
611  y2mqE += y2mE * q;
612  ansE += y2mqE + y2rE;
613  ansE += m * yE;
614  RETURN_M_E(ans, ansE);
615 }
616 
617 
618 static double2
619 pow4_D(double a, double b, double c, double d, int m)
620 {
621  /* Compute ((a+b)/(c+d)) ^ m */
622  double2 A, C, X;
623  if (m <= 0){
624  if (m == 0) {
625  return DD_C_ONE;
626  }
627  return pow4_D(c, d, a, b, -m);
628  }
629  A = add_dd(a, b);
630  C = add_dd(c, d);
631  if (dd_is_zero(A)) {
632  return (dd_is_zero(C) ? DD_C_NAN : DD_C_ZERO);
633  }
634  if (dd_is_zero(C)) {
635  return (dd_is_negative(A) ? DD_C_NEGINF : DD_C_INF);
636  }
637  X = div_DD(A, C);
638  return pow_D(X, m);
639 }
640 
641 static double
642 pow4(double a, double b, double c, double d, int m)
643 {
644  double2 ret = pow4_D(a, b, c, d, m);
645  return dd_to_double(ret);
646 }
647 
648 
649 static double2
650 logpow4_D(double a, double b, double c, double d, int m)
651 {
652  /*
653  * Compute log(((a+b)/(c+d)) ^ m)
654  * == m * log((a+b)/(c+d))
655  * == m * log( 1 + (a+b-c-d)/(c+d))
656  */
657  double2 ans;
658  double2 A, C, X;
659  if (m == 0) {
660  return DD_C_ZERO;
661  }
662  A = add_dd(a, b);
663  C = add_dd(c, d);
664  if (dd_is_zero(A)) {
665  return (dd_is_zero(C) ? DD_C_ZERO : DD_C_NEGINF);
666  }
667  if (dd_is_zero(C)) {
668  return DD_C_INF;
669  }
670  X = div_DD(A, C);
671  assert(X.x[0] >= 0);
672  if (0.5 <= X.x[0] && X.x[0] <= 1.5) {
673  double2 A1 = sub_DD(A, C);
674  double2 X1 = div_DD(A1, C);
675  ans = log1pD(X1);
676  } else {
677  ans = logD(X);
678  }
679  ans = mul_dD(m, ans);
680  return ans;
681 }
682 
683 static double
684 logpow4(double a, double b, double c, double d, int m)
685 {
686  double2 ans = logpow4_D(a, b, c, d, m);
687  return dd_to_double(ans);
688 }
689 
690 /*
691  * Compute a single term in the summation, A_v(n, x):
692  * A_v(n, x) = Binomial(n,v) * (1-x-v/n)^(n-v) * (x+v/n)^(v-1)
693  */
694 static void
695 computeAv(int n, double x, int v, double2 Cman, int Cexpt,
696  double2 *pt1, double2 *pt2, double2 *pAv)
697 {
698  int t1E, t2E, ansE;
699  double2 Av;
700  double2 t2x = sub_Dd(div_dd(n - v, n), x); /* 1 - x - v/n */
701  double2 t2 = pow2Scaled_D(t2x, n-v, &t2E);
702  double2 t1x = add_Dd(div_dd(v, n), x); /* x + v/n */
703  double2 t1 = pow2Scaled_D(t1x, v-1, &t1E);
704  double2 ans = mul_DD(t1, t2);
705  ans = mul_DD(ans, Cman);
706  ansE = Cexpt + t1E + t2E;
707  Av = ldexpD(ans, ansE);
708  *pAv = Av;
709  *pt1 = t1;
710  *pt2 = t2;
711 }
712 
713 
714 static ThreeProbs
715 _smirnov(int n, double x)
716 {
717  double nx, alpha;
718  double2 AjSum = DD_C_ZERO;
719  double2 dAjSum = DD_C_ZERO;
720  double cdf, sf, pdf;
721 
722  int bUseUpperSum;
723  int nxfl, n1mxfl, n1mxceil;
724  ThreeProbs ret;
725 
726  if (!(n > 0 && x >= 0.0 && x <= 1.0)) {
727  RETURN_3PROBS(NAN, NAN, NAN);
728  }
729  if (n == 1) {
730  RETURN_3PROBS(1-x, x, 1.0);
731  }
732  if (x == 0.0) {
733  RETURN_3PROBS(1.0, 0.0, 1.0);
734  }
735  if (x == 1.0) {
736  RETURN_3PROBS(0.0, 1.0, 0.0);
737  }
738 
739  alpha = modNX(n, x, &nxfl, &nx);
740  n1mxfl = n - nxfl - (alpha == 0 ? 0 : 1);
741  n1mxceil = n - nxfl;
742  /*
743  * If alpha is 0, don't actually want to include the last term
744  * in either the lower or upper summations.
745  */
746  if (alpha == 0) {
747  n1mxfl -= 1;
748  n1mxceil += 1;
749  }
750 
751  /* Special case: x <= 1/n */
752  if (nxfl == 0 || (nxfl == 1 && alpha == 0)) {
753  double t = pow2(1, x, n-1);
754  pdf = (nx + 1) * t / (1+x);
755  cdf = x * t;
756  sf = 1 - cdf;
757  /* Adjust if x=1/n *exactly* */
758  if (nxfl == 1) {
759  assert(alpha == 0);
760  pdf -= 0.5;
761  }
762  RETURN_3PROBS(sf, cdf, pdf);
763  }
764  /* Special case: x is so big, the sf underflows double64 */
765  if (-2 * n * x*x < MINLOG) {
766  RETURN_3PROBS(0, 1, 0);
767  }
768  /* Special case: x >= 1 - 1/n */
769  if (nxfl >= n-1) {
770  sf = pow2(1, -x, n);
771  cdf = 1 - sf;
772  pdf = n * sf/(1-x);
773  RETURN_3PROBS(sf, cdf, pdf);
774  }
775  /* Special case: n is so big, take too long to compute */
776  if (n > SMIRNOV_MAX_COMPUTE_N) {
777  /* p ~ e^(-(6nx+1)^2 / 18n) */
778  double logp = -pow(6.0*n*x+1, 2)/18.0/n;
779  /* Maximise precision for small p-value. */
780  if (logp < -M_LN2) {
781  sf = exp(logp);
782  cdf = 1 - sf;
783  } else {
784  cdf = -expm1(logp);
785  sf = 1 - cdf;
786  }
787  pdf = (6.0*n*x+1) * 2 * sf/3;
788  RETURN_3PROBS(sf, cdf, pdf);
789  }
790  {
791  /*
792  * Use the upper sum if n is large enough, and x is small enough and
793  * the number of terms is going to be small enough.
794  * Otherwise it just drops accuracy, about 1.6bits * nUpperTerms
795  */
796  int nUpperTerms = n - n1mxceil + 1;
797  bUseUpperSum = (nUpperTerms <= 1 && x < 0.5);
798  bUseUpperSum = (bUseUpperSum ||
799  ((n >= SM_UPPERSUM_MIN_N)
800  && (nUpperTerms <= SM_UPPER_MAX_TERMS)
801  && (x <= 0.5 / sqrt(n))));
802  }
803 
804  {
805  int start=0, step=1, nTerms=n1mxfl+1;
806  int j, firstJ = 0;
807  int vmid = n/2;
808  double2 Cman = DD_C_ONE;
809  int Cexpt = 0;
810  double2 Aj, dAj, t1, t2, dAjCoeff;
811  double2 oneOverX = div_dd(1, x);
812 
813  if (bUseUpperSum) {
814  start = n;
815  step = -1;
816  nTerms = n - n1mxceil + 1;
817 
818  t1 = pow4_D(1, x, 1, 0, n - 1);
819  t2 = DD_C_ONE;
820  Aj = t1;
821 
822  dAjCoeff = div_dD(n - 1, add_dd(1, x));
823  dAjCoeff = add_DD(dAjCoeff, oneOverX);
824  } else {
825  t1 = oneOverX;
826  t2 = pow4_D(1, -x, 1, 0, n);
827  Aj = div_Dd(t2, x);
828 
829  dAjCoeff = div_DD(sub_dD(-1, mul_dd(n - 1, x)), sub_dd(1, x));
830  dAjCoeff = div_Dd(dAjCoeff, x);
831  dAjCoeff = add_DD(dAjCoeff, oneOverX);
832  }
833 
834  dAj = mul_DD(Aj, dAjCoeff);
835  AjSum = add_DD(AjSum, Aj);
836  dAjSum = add_DD(dAjSum, dAj);
837 
838  updateBinomial(&Cman, &Cexpt, n, 0);
839  firstJ ++;
840 
841  for (j = firstJ; j < nTerms; j += 1) {
842  int v = start + j * step;
843 
844  computeAv(n, x, v, Cman, Cexpt, &t1, &t2, &Aj);
845 
846  if (dd_isfinite(Aj) && !dd_is_zero(Aj)) {
847  /* coeff = 1/x + (j-1)/(x+j/n) - (n-j)/(1-x-j/n) */
848  dAjCoeff = sub_DD(div_dD((n * (v - 1)), add_dd(nxfl + v, alpha)),
849  div_dD(((n - v) * n), sub_dd(n - nxfl - v, alpha)));
850  dAjCoeff = add_DD(dAjCoeff, oneOverX);
851  dAj = mul_DD(Aj, dAjCoeff);
852 
853  assert(dd_isfinite(Aj));
854  AjSum = add_DD(AjSum, Aj);
855  dAjSum = add_DD(dAjSum, dAj);
856  }
857  /* Safe to terminate early? */
858  if (!dd_is_zero(Aj)) {
859  if ((4*(nTerms-j) * fabs(dd_to_double(Aj)) < DBL_EPSILON * dd_to_double(AjSum))
860  && (j != nTerms - 1)) {
861  break;
862  }
863  }
864  else if (j > vmid) {
865  assert(dd_is_zero(Aj));
866  break;
867  }
868 
869  updateBinomial(&Cman, &Cexpt, n, j);
870  }
871  assert(dd_isfinite(AjSum));
872  assert(dd_isfinite(dAjSum));
873  {
874  double2 derivD = mul_dD(x, dAjSum);
875  double2 probD = mul_dD(x, AjSum);
876  double deriv = dd_to_double(derivD);
877  double prob = dd_to_double(probD);
878 
879  assert (nx != 1 || alpha > 0);
880  if (step < 0) {
881  cdf = prob;
882  sf = 1-prob;
883  pdf = deriv;
884  } else {
885  cdf = 1-prob;
886  sf = prob;
887  pdf = -deriv;
888  }
889  }
890  }
891 
892  pdf = MAX(0, pdf);
893  cdf = CLIP(cdf, 0, 1);
894  sf = CLIP(sf, 0, 1);
895  RETURN_3PROBS(sf, cdf, pdf);
896 }
897 
898 /*
899  * Functional inverse of Smirnov distribution
900  * finds x such that smirnov(n, x) = psf; smirnovc(n, x) = pcdf).
901  */
902 static double
903 _smirnovi(int n, double psf, double pcdf)
904 {
905  /*
906  * Need to use a bracketing NR algorithm here and be very careful
907  * about the starting point.
908  */
909  double x, logpcdf;
910  int iterations = 0;
911  int function_calls = 0;
912  double a=0, b=1;
913  double maxlogpcdf, psfrootn;
914  double dx, dxold;
915 
916  if (!(n > 0 && psf >= 0.0 && pcdf >= 0.0 && pcdf <= 1.0 && psf <= 1.0)) {
917  sf_error("smirnovi", SF_ERROR_DOMAIN, NULL);
918  return (NAN);
919  }
920  if (fabs(1.0 - pcdf - psf) > 4* DBL_EPSILON) {
921  sf_error("smirnovi", SF_ERROR_DOMAIN, NULL);
922  return (NAN);
923  }
924  /* STEP 1: Handle psf==0, or pcdf == 0 */
925  if (pcdf == 0.0) {
926  return 0.0;
927  }
928  if (psf == 0.0) {
929  return 1.0;
930  }
931  /* STEP 2: Handle n=1 */
932  if (n == 1) {
933  return pcdf;
934  }
935 
936  /* STEP 3 Handle psf *very* close to 0. Correspond to (n-1)/n < x < 1 */
937  psfrootn = pow(psf, 1.0 / n);
938  /* xmin > 1 - 1.0 / n */
939  if (n < 150 && n*psfrootn <= 1) {
940  /* Solve exactly. */
941  x = 1 - psfrootn;
942  return x;
943  }
944 
945  logpcdf = (pcdf < 0.5 ? log(pcdf) : log1p(-psf));
946 
947  /*
948  * STEP 4 Find bracket and initial estimate for use in N-R
949  * 4(a) Handle 0 < x <= 1/n: pcdf = x * (1+x)^*(n-1)
950  */
951  maxlogpcdf = logpow4(1, 0.0, n, 0, 1) + logpow4(n, 1, n, 0, n - 1);
952  if (logpcdf <= maxlogpcdf) {
953  double xmin = pcdf / SCIPY_El;
954  double xmax = pcdf;
955  double P1 = pow4(n, 1, n, 0, n - 1) / n;
956  double R = pcdf/P1;
957  double z0 = R;
958  /*
959  * Do one iteration of N-R solving: z*e^(z-1) = R, with z0=pcdf/P1
960  * z <- z - (z exp(z-1) - pcdf)/((z+1)exp(z-1))
961  * If z_0 = R, z_1 = R(1-exp(1-R))/(R+1)
962  */
963  if (R >= 1) {
964  /*
965  * R=1 is OK;
966  * R>1 can happen due to truncation error for x = (1-1/n)+-eps
967  */
968  R = 1;
969  x = R/n;
970  return x;
971  }
972  z0 = (z0*z0 + R * exp(1-z0))/(1+z0);
973  x = z0/n;
974  a = xmin*(1 - 4 * DBL_EPSILON);
975  a = MAX(a, 0);
976  b = xmax * (1 + 4 * DBL_EPSILON);
977  b = MIN(b, 1.0/n);
978  x = CLIP(x, a, b);
979  }
980  else
981  {
982  /* 4(b) : 1/n < x < (n-1)/n */
983  double xmin = 1 - psfrootn;
984  double logpsf = (psf < 0.5 ? log(psf) : log1p(-pcdf));
985  double xmax = sqrt(-logpsf / (2.0L * n));
986  double xmax6 = xmax - 1.0L / (6 * n);
987  a = xmin;
988  b = xmax;
989  /* Allow for a little rounding error */
990  a *= 1 - 4 * DBL_EPSILON;
991  b *= 1 + 4 * DBL_EPSILON;
992  a = MAX(xmin, 1.0/n);
993  b = MIN(xmax, 1-1.0/n);
994  x = xmax6;
995  }
996  if (x < a || x > b) {
997  x = (a + b)/2;
998  }
999  assert (x < 1);
1000 
1001  /*
1002  * Skip computing fa, fb as that takes cycles and the exact values
1003  * are not needed.
1004  */
1005 
1006  /* STEP 5 Run N-R.
1007  * smirnov should be well-enough behaved for NR starting at this location.
1008  * Use smirnov(n, x)-psf, or pcdf - smirnovc(n, x), whichever has smaller p.
1009  */
1010  dxold = b - a;
1011  dx = dxold;
1012  do {
1013  double dfdx, x0 = x, deltax, df;
1014  assert(x < 1);
1015  assert(x > 0);
1016  {
1017  ThreeProbs probs = _smirnov(n, x0);
1018  ++function_calls;
1019  df = ((pcdf < 0.5) ? (pcdf - probs.cdf) : (probs.sf - psf));
1020  dfdx = -probs.pdf;
1021  }
1022  if (df == 0) {
1023  return x;
1024  }
1025  /* Update the bracketing interval */
1026  if (df > 0 && x > a) {
1027  a = x;
1028  } else if (df < 0 && x < b) {
1029  b = x;
1030  }
1031 
1032  if (dfdx == 0) {
1033  /*
1034  * x was not within tolerance, but now we hit a 0 derivative.
1035  * This implies that x >> 1/sqrt(n), and even then |smirnovp| >= |smirnov|
1036  * so this condition is unexpected. Do a bisection step.
1037  */
1038  x = (a+b)/2;
1039  deltax = x0 - x;
1040  } else {
1041  deltax = df / dfdx;
1042  x = x0 - deltax;
1043  }
1044  /*
1045  * Check out-of-bounds.
1046  * Not expecting this to happen ofen --- smirnov is convex near x=1 and
1047  * concave near x=0, and we should be approaching from the correct side.
1048  * If out-of-bounds, replace x with a midpoint of the bracket.
1049  * Also check fast enough convergence.
1050  */
1051  if ((a <= x) && (x <= b) && (fabs(2 * deltax) <= fabs(dxold) || fabs(dxold) < 256 * DBL_EPSILON)) {
1052  dxold = dx;
1053  dx = deltax;
1054  } else {
1055  dxold = dx;
1056  dx = dx / 2;
1057  x = (a + b) / 2;
1058  deltax = x0 - x;
1059  }
1060  /*
1061  * Note that if psf is close to 1, f(x) -> 1, f'(x) -> -1.
1062  * => abs difference |x-x0| is approx |f(x)-p| >= DBL_EPSILON,
1063  * => |x-x0|/x >= DBL_EPSILON/x.
1064  * => cannot use a purely relative criteria as it will fail for x close to 0.
1065  */
1066  if (_within_tol(x, x0, (psf < 0.5 ? 0 : _xtol), _rtol)) {
1067  break;
1068  }
1069  if (++iterations > MAXITER) {
1070  sf_error("smirnovi", SF_ERROR_SLOW, NULL);
1071  return (x);
1072  }
1073  } while (1);
1074  return x;
1075 }
1076 
1077 
1078 double
1079 smirnov(int n, double d)
1080 {
1081  ThreeProbs probs;
1082  if (isnan(d)) {
1083  return NAN;
1084  }
1085  probs = _smirnov(n, d);
1086  return probs.sf;
1087 }
1088 
1089 double
1090 smirnovc(int n, double d)
1091 {
1092  ThreeProbs probs;
1093  if (isnan(d)) {
1094  return NAN;
1095  }
1096  probs = _smirnov(n, d);
1097  return probs.cdf;
1098 }
1099 
1100 
1101 /*
1102  * Derivative of smirnov(n, d)
1103  * One interior point of discontinuity at d=1/n.
1104 */
1105 double
1106 smirnovp(int n, double d)
1107 {
1108  ThreeProbs probs;
1109  if (!(n > 0 && d >= 0.0 && d <= 1.0)) {
1110  return (NAN);
1111  }
1112  if (n == 1) {
1113  /* Slope is always -1 for n=1, even at d = 1.0 */
1114  return -1.0;
1115  }
1116  if (d == 1.0) {
1117  return -0.0;
1118  }
1119  /*
1120  * If d is 0, the derivative is discontinuous, but approaching
1121  * from the right the limit is -1
1122  */
1123  if (d == 0.0) {
1124  return -1.0;
1125  }
1126  probs = _smirnov(n, d);
1127  return -probs.pdf;
1128 }
1129 
1130 
1131 double
1132 smirnovi(int n, double p)
1133 {
1134  if (isnan(p)) {
1135  return NAN;
1136  }
1137  return _smirnovi(n, p, 1-p);
1138 }
1139 
1140 double
1141 smirnovci(int n, double p)
1142 {
1143  if (isnan(p)) {
1144  return NAN;
1145  }
1146  return _smirnovi(n, 1-p, p);
1147 }
DD_C_NAN
#define DD_C_NAN
Definition: dd_real.h:105
dd_floor
static double2 dd_floor(const double2 a)
Definition: dd_real_idefs.h:207
dd_to_int
static int dd_to_int(const double2 a)
Definition: dd_real_idefs.h:97
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
gtsam.examples.DogLegOptimizerExample.int
int
Definition: DogLegOptimizerExample.py:111
DD_C_ONE
const double2 DD_C_ONE
Definition: dd_real.c:47
ldexpD
#define ldexpD
Definition: kolmogorov.c:149
add_Dd
#define add_Dd
Definition: kolmogorov.c:140
modNX
static double modNX(int n, double x, int *pNXFloor, double *pNX)
Definition: kolmogorov.c:460
div_dD
#define div_dD
Definition: kolmogorov.c:147
dd_real.h
pow4
static double pow4(double a, double b, double c, double d, int m)
Definition: kolmogorov.c:642
D
MatrixXcd D
Definition: EigenSolver_EigenSolver_MatrixType.cpp:14
double2
Definition: dd_real.h:74
ThreeProbs
struct ThreeProbs ThreeProbs
MAXITER
#define MAXITER
Definition: igam.c:110
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:147
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
d
static const double d[K][N]
Definition: igam.h:11
CLIP
#define CLIP(X, LOW, HIGH)
Definition: kolmogorov.c:85
_smirnovi
static double _smirnovi(int n, double psf, double pcdf)
Definition: kolmogorov.c:903
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
b
Scalar * b
Definition: benchVecAdd.cpp:17
RETURN_M_E
#define RETURN_M_E(MAND, EXPT)
Definition: kolmogorov.c:551
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
smirnovi
double smirnovi(int n, double p)
Definition: kolmogorov.c:1132
ret
DenseIndex ret
Definition: level1_cplx_impl.h:44
ThreeProbs::pdf
double pdf
Definition: kolmogorov.c:107
updateBinomial
static void updateBinomial(double2 *Cman, int *Cexpt, int n, int j)
Definition: kolmogorov.c:495
MIN
#define MIN(a, b)
Definition: kolmogorov.c:87
_MAX_EXPONENT
#define _MAX_EXPONENT
Definition: kolmogorov.c:549
mul_DD
#define mul_DD
Definition: kolmogorov.c:138
_within_tol
static int _within_tol(double x, double y, double atol, double rtol)
Definition: kolmogorov.c:120
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
X
#define X
Definition: icosphere.cpp:20
isnan
#define isnan(X)
Definition: main.h:93
RETURN_3PROBS
#define RETURN_3PROBS(PSF, PCDF, PDF)
Definition: kolmogorov.c:110
sub_DD
#define sub_DD
Definition: kolmogorov.c:137
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
dd_inv
static double2 dd_inv(const double2 a)
Definition: dd_real_idefs.h:503
result
Values result
Definition: OdometryOptimize.cpp:8
kolmogp
double kolmogp(double x)
Definition: kolmogorov.c:402
DD_C_ZERO
const double2 DD_C_ZERO
Definition: dd_real.c:46
dd_add_d_d
static double2 dd_add_d_d(double a, double b)
Definition: dd_real_idefs.h:286
smirnovci
double smirnovci(int n, double p)
Definition: kolmogorov.c:1141
_xtol
static const double _xtol
Definition: kolmogorov.c:116
frexpD
#define frexpD
Definition: kolmogorov.c:148
A
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
SF_ERROR_SLOW
@ SF_ERROR_SLOW
Definition: sf_error.h:13
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
computeAv
static void computeAv(int n, double x, int v, double2 Cman, int Cexpt, double2 *pt1, double2 *pt2, double2 *pAv)
Definition: kolmogorov.c:695
n
int n
Definition: BiCGSTAB_simple.cpp:1
simple::p2
static Point3 p2
Definition: testInitializePose3.cpp:51
MAX
#define MAX(a, b)
Definition: kolmogorov.c:90
smirnovc
double smirnovc(int n, double d)
Definition: kolmogorov.c:1090
smirnovp
double smirnovp(int n, double d)
Definition: kolmogorov.c:1106
logpow4
static double logpow4(double a, double b, double c, double d, int m)
Definition: kolmogorov.c:684
SCIPY_El
#define SCIPY_El
Definition: mconf.h:129
div_dd
#define div_dd
Definition: kolmogorov.c:135
expm1
double expm1(double x)
Definition: unity.c:106
A
Definition: test_numpy_dtypes.cpp:298
smirnov
double smirnov(int n, double d)
Definition: kolmogorov.c:1079
kolmogci
double kolmogci(double p)
Definition: kolmogorov.c:429
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:1984
P1
static double P1[]
Definition: jv.c:552
add_DD
#define add_DD
Definition: kolmogorov.c:136
pow4_D
static double2 pow4_D(double a, double b, double c, double d, int m)
Definition: kolmogorov.c:619
dd_to_double
static double dd_to_double(const double2 a)
Definition: dd_real_idefs.h:90
log1pD
#define log1pD
Definition: kolmogorov.c:151
LOGSQRT2PI
#define LOGSQRT2PI
Definition: kolmogorov.c:100
SMIRNOV_MAX_COMPUTE_N
const int SMIRNOV_MAX_COMPUTE_N
Definition: kolmogorov.c:70
log1p
double log1p(double x)
Definition: unity.c:49
ThreeProbs::sf
double sf
Definition: kolmogorov.c:105
x0
static Symbol x0('x', 0)
pow2Scaled_D
static double2 pow2Scaled_D(double2 a, int m, int *pExponent)
Definition: kolmogorov.c:557
gtsam::utils.visual_isam.step
def step(data, isam, result, truth, currPoseIndex, isamArgs=())
Definition: visual_isam.py:82
dd_hi
static double dd_hi(const double2 a)
Definition: dd_real_idefs.h:41
dd_isfinite
static int dd_isfinite(const double2 a)
Definition: dd_real_idefs.h:53
L
MatrixXd L
Definition: LLT_example.cpp:6
simple::p3
static Point3 p3
Definition: testInitializePose3.cpp:53
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
dd_is_zero
static int dd_is_zero(const double2 a)
Definition: dd_real_idefs.h:65
kolmogi
double kolmogi(double p)
Definition: kolmogorov.c:417
_rtol
static const double _rtol
Definition: kolmogorov.c:117
MINLOG
double MINLOG
Definition: const.c:60
MIN_EXPABLE
static const int MIN_EXPABLE
Definition: kolmogorov.c:97
pow2
static double pow2(double a, double b, int m)
Definition: kolmogorov.c:540
_kolmogi
static double _kolmogi(double psf, double pcdf)
Definition: kolmogorov.c:256
y
Scalar * y
Definition: level1_cplx_impl.h:124
nextPowerOf2
static double nextPowerOf2(double x)
Definition: kolmogorov.c:444
gtsam::utils.numerical_derivative.X1
X1
Definition: numerical_derivative.py:25
kolmogc
double kolmogc(double x)
Definition: kolmogorov.c:393
Eigen::ArrayBase::pow
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_pow_op< typename Derived::Scalar, typename ExponentDerived::Scalar >, const Derived, const ExponentDerived > pow(const Eigen::ArrayBase< Derived > &x, const Eigen::ArrayBase< ExponentDerived > &exponents)
Definition: GlobalFunctions.h:143
sub_Dd
#define sub_Dd
Definition: kolmogorov.c:142
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
_kolmogorov
static ThreeProbs _kolmogorov(double x)
Definition: kolmogorov.c:159
mconf.h
DD_C_NEGINF
#define DD_C_NEGINF
Definition: dd_real.h:107
C
Matrix< Scalar, Dynamic, Dynamic > C
Definition: bench_gemm.cpp:50
DD_C_INF
#define DD_C_INF
Definition: dd_real.h:106
_smirnov
static ThreeProbs _smirnov(int n, double x)
Definition: kolmogorov.c:715
sub_dD
#define sub_dD
Definition: kolmogorov.c:143
SM_UPPER_MAX_TERMS
#define SM_UPPER_MAX_TERMS
Definition: kolmogorov.c:78
p
float * p
Definition: Tutorial_Map_using.cpp:9
A1
static const double A1[]
Definition: expn.h:6
sub_dd
#define sub_dd
Definition: kolmogorov.c:132
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
M_LN2
#define M_LN2
Definition: mconf.h:115
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
P
static double P[]
Definition: ellpe.c:68
kolmogorov
double kolmogorov(double x)
Definition: kolmogorov.c:384
div_DD
#define div_DD
Definition: kolmogorov.c:139
KOLMOG_CUTOVER
#define KOLMOG_CUTOVER
Definition: kolmogorov.c:63
mul_dD
#define mul_dD
Definition: kolmogorov.c:145
NULL
#define NULL
Definition: ccolamd.c:609
v3
Vector v3
Definition: testSerializationBase.cpp:40
M_PI
#define M_PI
Definition: mconf.h:117
ThreeProbs
Definition: kolmogorov.c:104
logpow4_D
static double2 logpow4_D(double a, double b, double c, double d, int m)
Definition: kolmogorov.c:650
div_Dd
#define div_Dd
Definition: kolmogorov.c:146
align_3::t
Point2 t(10, 10)
SM_UPPERSUM_MIN_N
#define SM_UPPERSUM_MIN_N
Definition: kolmogorov.c:79
SF_ERROR_DOMAIN
@ SF_ERROR_DOMAIN
Definition: sf_error.h:16
pow_D
static double2 pow_D(double2 a, int m)
Definition: kolmogorov.c:508
ThreeProbs::cdf
double cdf
Definition: kolmogorov.c:106
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
R
Rot2 R(Rot2::fromAngle(0.1))
logD
#define logD
Definition: kolmogorov.c:150
dd_is_negative
static int dd_is_negative(const double2 a)
Definition: dd_real_idefs.h:83
M_SQRT2
#define M_SQRT2
Definition: mconf.h:123
mul_dd
#define mul_dd
Definition: kolmogorov.c:133
add_dd
#define add_dd
Definition: kolmogorov.c:131


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:33:01