otsdaq_utilities  v2_05_02_indev
JSRootMath.js
1 
4 (function( factory ) {
5  if ( typeof define === "function" && define.amd ) {
6  // AMD. Register as an anonymous module.
7  define( ['JSRootCore'], factory );
8  } else {
9  if (typeof JSROOT == 'undefined')
10  throw new Error("This extension requires JSRootCore.js", "JSRootMath.js");
11 
12  if (typeof JSROOT.Math == "object")
13  throw new Error("This JSROOT Math already loaded", "JSRootMath.js");
14 
15  factory(JSROOT);
16  }
17 } (function(JSROOT) {
18  // math methods for Javascript ROOT
19 
20  JSROOT.Math = {};
21 
22  JSROOT.Math.lgam = function( x ) {
23  var p, q, u, w, z, i, sgngam = 1;
24  var kMAXLGM = 2.556348e305;
25  var LS2PI = 0.91893853320467274178;
26 
27  var A = [
28  8.11614167470508450300E-4,
29  -5.95061904284301438324E-4,
30  7.93650340457716943945E-4,
31  -2.77777777730099687205E-3,
32  8.33333333333331927722E-2
33  ];
34 
35  var B = [
36  -1.37825152569120859100E3,
37  -3.88016315134637840924E4,
38  -3.31612992738871184744E5,
39  -1.16237097492762307383E6,
40  -1.72173700820839662146E6,
41  -8.53555664245765465627E5
42  ];
43 
44  var C = [
45  /* 1.00000000000000000000E0, */
46  -3.51815701436523470549E2,
47  -1.70642106651881159223E4,
48  -2.20528590553854454839E5,
49  -1.13933444367982507207E6,
50  -2.53252307177582951285E6,
51  -2.01889141433532773231E6
52  ];
53 
54  if (x >= Number.POSITIVE_INFINITY)
55  return(Number.POSITIVE_INFINITY);
56 
57  if ( x < -34.0 ) {
58  q = -x;
59  w = this.lgam(q);
60  p = Math.floor(q);
61  if ( p==q )//_unur_FP_same(p,q)
62  return (Number.POSITIVE_INFINITY);
63  i = Math.round(p);
64  if ( (i & 1) == 0 )
65  sgngam = -1;
66  else
67  sgngam = 1;
68  z = q - p;
69  if ( z > 0.5 ) {
70  p += 1.0;
71  z = p - q;
72  }
73  z = q * Math.sin( Math.PI * z );
74  if ( z < 1e-300 )
75  return (Number.POSITIVE_INFINITY);
76  z = Math.log(Math.PI) - Math.log( z ) - w;
77  return( z );
78  }
79  if ( x < 13.0 ) {
80  z = 1.0;
81  p = 0.0;
82  u = x;
83  while ( u >= 3.0 ) {
84  p -= 1.0;
85  u = x + p;
86  z *= u;
87  }
88  while ( u < 2.0 ) {
89  if ( u < 1e-300 )
90  return (Number.POSITIVE_INFINITY);
91  z /= u;
92  p += 1.0;
93  u = x + p;
94  }
95  if ( z < 0.0 ) {
96  sgngam = -1;
97  z = -z;
98  }
99  else
100  sgngam = 1;
101  if ( u == 2.0 )
102  return( Math.log(z) );
103  p -= 2.0;
104  x = x + p;
105  p = x * this.Polynomialeval(x, B, 5 ) / this.Polynomial1eval( x, C, 6);
106  return( Math.log(z) + p );
107  }
108  if ( x > kMAXLGM )
109  return( sgngam * Number.POSITIVE_INFINITY );
110 
111  q = ( x - 0.5 ) * Math.log(x) - x + LS2PI;
112  if ( x > 1.0e8 )
113  return( q );
114 
115  p = 1.0/(x*x);
116  if ( x >= 1000.0 )
117  q += ((7.9365079365079365079365e-4 * p
118  - 2.7777777777777777777778e-3) *p
119  + 0.0833333333333333333333) / x;
120  else
121  q += this.Polynomialeval( p, A, 4 ) / x;
122  return( q );
123  };
124 
125  /*
126  * calculates a value of a polynomial of the form:
127  * a[0]x^N+a[1]x^(N-1) + ... + a[N]
128  */
129  JSROOT.Math.Polynomialeval = function(x, a, N) {
130  if (N==0) return a[0];
131  else {
132  var pom = a[0];
133  for (var i=1; i <= N; ++i)
134  pom = pom *x + a[i];
135  return pom;
136  }
137  };
138 
139  /*
140  * calculates a value of a polynomial of the form:
141  * x^N+a[0]x^(N-1) + ... + a[N-1]
142  */
143  JSROOT.Math.Polynomial1eval = function(x, a, N) {
144  if (N==0) return a[0];
145  else {
146  var pom = x + a[0];
147  for (var i=1; i < N; ++i)
148  pom = pom *x + a[i];
149  return pom;
150  }
151  };
152 
153  JSROOT.Math.ndtri = function( y0 ) {
154  if ( y0 <= 0.0 )
155  return( Number.NEGATIVE_INFINITY );
156  if ( y0 >= 1.0 )
157  return( Number.POSITIVE_INFINITY );
158 
159  var P0 = new Array(
160  -5.99633501014107895267E1,
161  9.80010754185999661536E1,
162  -5.66762857469070293439E1,
163  1.39312609387279679503E1,
164  -1.23916583867381258016E0
165  );
166 
167  var Q0 = new Array(
168  1.95448858338141759834E0,
169  4.67627912898881538453E0,
170  8.63602421390890590575E1,
171  -2.25462687854119370527E2,
172  2.00260212380060660359E2,
173  -8.20372256168333339912E1,
174  1.59056225126211695515E1,
175  -1.18331621121330003142E0
176  );
177 
178  var P1 = new Array(
179  4.05544892305962419923E0,
180  3.15251094599893866154E1,
181  5.71628192246421288162E1,
182  4.40805073893200834700E1,
183  1.46849561928858024014E1,
184  2.18663306850790267539E0,
185  -1.40256079171354495875E-1,
186  -3.50424626827848203418E-2,
187  -8.57456785154685413611E-4
188  );
189 
190  var Q1 = new Array(
191  1.57799883256466749731E1,
192  4.53907635128879210584E1,
193  4.13172038254672030440E1,
194  1.50425385692907503408E1,
195  2.50464946208309415979E0,
196  -1.42182922854787788574E-1,
197  -3.80806407691578277194E-2,
198  -9.33259480895457427372E-4
199  );
200 
201  var P2 = new Array(
202  3.23774891776946035970E0,
203  6.91522889068984211695E0,
204  3.93881025292474443415E0,
205  1.33303460815807542389E0,
206  2.01485389549179081538E-1,
207  1.23716634817820021358E-2,
208  3.01581553508235416007E-4,
209  2.65806974686737550832E-6,
210  6.23974539184983293730E-9
211  );
212 
213  var Q2 = new Array(
214  6.02427039364742014255E0,
215  3.67983563856160859403E0,
216  1.37702099489081330271E0,
217  2.16236993594496635890E-1,
218  1.34204006088543189037E-2,
219  3.28014464682127739104E-4,
220  2.89247864745380683936E-6,
221  6.79019408009981274425E-9
222  );
223 
224  var s2pi = 2.50662827463100050242e0;
225  var code = 1;
226  var y = y0;
227  var x, z, y2, x0, x1;
228 
229  if ( y > (1.0 - 0.13533528323661269189) ) {
230  y = 1.0 - y;
231  code = 0;
232  }
233  if ( y > 0.13533528323661269189 ) {
234  y = y - 0.5;
235  y2 = y * y;
236  x = y + y * (y2 * this.Polynomialeval( y2, P0, 4)/ this.Polynomial1eval( y2, Q0, 8 ));
237  x = x * s2pi;
238  return(x);
239  }
240  x = Math.sqrt( -2.0 * Math.log(y) );
241  x0 = x - Math.log(x)/x;
242  z = 1.0/x;
243  if ( x < 8.0 )
244  x1 = z * this.Polynomialeval( z, P1, 8 )/ this.Polynomial1eval( z, Q1, 8 );
245  else
246  x1 = z * this.Polynomialeval( z, P2, 8 )/ this.Polynomial1eval( z, Q2, 8 );
247  x = x0 - x1;
248  if ( code != 0 )
249  x = -x;
250  return( x );
251  }
252 
253  JSROOT.Math.igam = function(a,x) {
254  var kMACHEP = 1.11022302462515654042363166809e-16;
255  var kMAXLOG = 709.782712893383973096206318587;
256  var ans, ax, c, r;
257 
258  // LM: for negative values returns 1.0 instead of zero
259  // This is correct if a is a negative integer since Gamma(-n) = +/- inf
260  if (a <= 0) return 1.0;
261 
262  if (x <= 0) return 0.0;
263 
264  if( (x > 1.0) && (x > a ) )
265  return( 1.0 - this.igamc(a,x) );
266 
267  /* Compute x**a * exp(-x) / gamma(a) */
268  ax = a * Math.log(x) - x - this.lgam(a);
269  if( ax < -kMAXLOG )
270  return( 0.0 );
271 
272  ax = Math.exp(ax);
273 
274  /* power series */
275  r = a;
276  c = 1.0;
277  ans = 1.0;
278 
279  do
280  {
281  r += 1.0;
282  c *= x/r;
283  ans += c;
284  }
285  while( c/ans > kMACHEP );
286 
287  return( ans * ax/a );
288  }
289 
290  JSROOT.Math.igamc = function(a,x) {
291  var kMACHEP = 1.11022302462515654042363166809e-16;
292  var kMAXLOG = 709.782712893383973096206318587;
293 
294  var kBig = 4.503599627370496e15;
295  var kBiginv = 2.22044604925031308085e-16;
296 
297  var ans, ax, c, yc, r, t, y, z;
298  var pk, pkm1, pkm2, qk, qkm1, qkm2;
299 
300  // LM: for negative values returns 0.0
301  // This is correct if a is a negative integer since Gamma(-n) = +/- inf
302  if (a <= 0) return 0.0;
303 
304  if (x <= 0) return 1.0;
305 
306  if( (x < 1.0) || (x < a) )
307  return ( 1.0 - JSROOT.Math.igam(a,x) );
308 
309  ax = a * Math.log(x) - x - JSROOT.Math.lgam(a);
310  if( ax < -kMAXLOG )
311  return( 0.0 );
312 
313  ax = Math.exp(ax);
314 
315  /* continued fraction */
316  y = 1.0 - a;
317  z = x + y + 1.0;
318  c = 0.0;
319  pkm2 = 1.0;
320  qkm2 = x;
321  pkm1 = x + 1.0;
322  qkm1 = z * x;
323  ans = pkm1/qkm1;
324 
325  do
326  {
327  c += 1.0;
328  y += 1.0;
329  z += 2.0;
330  yc = y * c;
331  pk = pkm1 * z - pkm2 * yc;
332  qk = qkm1 * z - qkm2 * yc;
333  if(qk)
334  {
335  r = pk/qk;
336  t = Math.abs( (ans - r)/r );
337  ans = r;
338  }
339  else
340  t = 1.0;
341  pkm2 = pkm1;
342  pkm1 = pk;
343  qkm2 = qkm1;
344  qkm1 = qk;
345  if( Math.abs(pk) > kBig )
346  {
347  pkm2 *= kBiginv;
348  pkm1 *= kBiginv;
349  qkm2 *= kBiginv;
350  qkm1 *= kBiginv;
351  }
352  }
353  while( t > kMACHEP );
354 
355  return( ans * ax );
356  }
357 
358 
359  JSROOT.Math.igami = function(a, y0) {
360  var x0, x1, x, yl, yh, y, d, lgm, dithresh;
361  var i, dir;
362  var kMACHEP = 1.11022302462515654042363166809e-16;
363 
364  // check the domain
365  if (a <= 0) {
366  alert("igami : Wrong domain for parameter a (must be > 0)");
367  return 0;
368  }
369  if (y0 <= 0) {
370  return Number.POSITIVE_INFINITY;
371  }
372  if (y0 >= 1) {
373  return 0;
374  }
375  /* bound the solution */
376  var kMAXNUM = Number.MAX_VALUE;
377  x0 = kMAXNUM;
378  yl = 0;
379  x1 = 0;
380  yh = 1.0;
381  dithresh = 5.0 * kMACHEP;
382 
383  /* approximation to inverse function */
384  d = 1.0/(9.0*a);
385  y = ( 1.0 - d - this.ndtri(y0) * Math.sqrt(d) );
386  x = a * y * y * y;
387 
388  lgm = this.lgam(a);
389 
390  for( i=0; i<10; ++i ) {
391  if ( x > x0 || x < x1 )
392  break;
393  y = this.igamc(a,x);
394  if ( y < yl || y > yh )
395  break;
396  if ( y < y0 ) {
397  x0 = x;
398  yl = y;
399  }
400  else {
401  x1 = x;
402  yh = y;
403  }
404  /* compute the derivative of the function at this point */
405  d = (a - 1.0) * Math.log(x) - x - lgm;
406  if ( d < -kMAXLOG )
407  break;
408  d = -Math.exp(d);
409  /* compute the step to the next approximation of x */
410  d = (y - y0)/d;
411  if ( Math.abs(d/x) < kMACHEP )
412  return( x );
413  x = x - d;
414  }
415  /* Resort to interval halving if Newton iteration did not converge. */
416  d = 0.0625;
417  if ( x0 == kMAXNUM ) {
418  if ( x <= 0.0 )
419  x = 1.0;
420  while ( x0 == kMAXNUM ) {
421  x = (1.0 + d) * x;
422  y = this.igamc( a, x );
423  if ( y < y0 ) {
424  x0 = x;
425  yl = y;
426  break;
427  }
428  d = d + d;
429  }
430  }
431  d = 0.5;
432  dir = 0;
433 
434  for( i=0; i<400; ++i ) {
435  x = x1 + d * (x0 - x1);
436  y = this.igamc( a, x );
437  lgm = (x0 - x1)/(x1 + x0);
438  if ( Math.abs(lgm) < dithresh )
439  break;
440  lgm = (y - y0)/y0;
441  if ( Math.abs(lgm) < dithresh )
442  break;
443  if ( x <= 0.0 )
444  break;
445  if ( y >= y0 ) {
446  x1 = x;
447  yh = y;
448  if ( dir < 0 ) {
449  dir = 0;
450  d = 0.5;
451  }
452  else if ( dir > 1 )
453  d = 0.5 * d + 0.5;
454  else
455  d = (y0 - yl)/(yh - yl);
456  dir += 1;
457  }
458  else {
459  x0 = x;
460  yl = y;
461  if ( dir > 0 ) {
462  dir = 0;
463  d = 0.5;
464  }
465  else if ( dir < -1 )
466  d = 0.5 * d;
467  else
468  d = (y0 - yl)/(yh - yl);
469  dir -= 1;
470  }
471  }
472  return( x );
473  };
474 
475  JSROOT.Math.gamma_quantile_c = function(z, alpha, theta) {
476  return theta * this.igami( alpha, z);
477  };
478 
479  JSROOT.Math.gamma_quantile = function(z, alpha, theta) {
480  return theta * this.igami( alpha, 1.- z);
481  };
482 
483  JSROOT.Math.log10 = function(n) {
484  return Math.log(n) / Math.log(10);
485  };
486 
487  JSROOT.Math.landau_pdf = function(x, xi, x0) {
488  // LANDAU pdf : algorithm from CERNLIB G110 denlan
489  // same algorithm is used in GSL
490  if (xi <= 0) return 0;
491  var v = (x - x0)/xi;
492  var u, ue, us, denlan;
493  var p1 = new Array(0.4259894875,-0.1249762550, 0.03984243700, -0.006298287635, 0.001511162253);
494  var q1 = new Array(1.0 ,-0.3388260629, 0.09594393323, -0.01608042283, 0.003778942063);
495  var p2 = new Array(0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411, 0.0001283617211);
496  var q2 = new Array(1.0 , 0.7428795082, 0.3153932961, 0.06694219548, 0.008790609714);
497  var p3 = new Array(0.1788544503, 0.09359161662,0.006325387654, 0.00006611667319,-0.000002031049101);
498  var q3 = new Array(1.0 , 0.6097809921, 0.2560616665, 0.04746722384, 0.006957301675);
499  var p4 = new Array(0.9874054407, 118.6723273, 849.2794360, -743.7792444, 427.0262186);
500  var q4 = new Array(1.0 , 106.8615961, 337.6496214, 2016.712389, 1597.063511);
501  var p5 = new Array(1.003675074, 167.5702434, 4789.711289, 21217.86767, -22324.94910);
502  var q5 = new Array(1.0 , 156.9424537, 3745.310488, 9834.698876, 66924.28357);
503  var p6 = new Array(1.000827619, 664.9143136, 62972.92665, 475554.6998, -5743609.109);
504  var q6 = new Array(1.0 , 651.4101098, 56974.73333, 165917.4725, -2815759.939);
505  var a1 = new Array(0.04166666667,-0.01996527778, 0.02709538966);
506  var a2 = new Array(-1.845568670,-4.284640743);
507 
508  if (v < -5.5) {
509  u = Math.exp(v+1.0);
510  if (u < 1e-10) return 0.0;
511  ue = Math.exp(-1/u);
512  us = Math.sqrt(u);
513  denlan = 0.3989422803*(ue/us)*(1+(a1[0]+(a1[1]+a1[2]*u)*u)*u);
514  } else if(v < -1) {
515  u = Math.exp(-v-1);
516  denlan = Math.exp(-u)*Math.sqrt(u)*
517  (p1[0]+(p1[1]+(p1[2]+(p1[3]+p1[4]*v)*v)*v)*v)/
518  (q1[0]+(q1[1]+(q1[2]+(q1[3]+q1[4]*v)*v)*v)*v);
519  } else if(v < 1) {
520  denlan = (p2[0]+(p2[1]+(p2[2]+(p2[3]+p2[4]*v)*v)*v)*v)/
521  (q2[0]+(q2[1]+(q2[2]+(q2[3]+q2[4]*v)*v)*v)*v);
522  } else if(v < 5) {
523  denlan = (p3[0]+(p3[1]+(p3[2]+(p3[3]+p3[4]*v)*v)*v)*v)/
524  (q3[0]+(q3[1]+(q3[2]+(q3[3]+q3[4]*v)*v)*v)*v);
525  } else if(v < 12) {
526  u = 1/v;
527  denlan = u*u*(p4[0]+(p4[1]+(p4[2]+(p4[3]+p4[4]*u)*u)*u)*u)/
528  (q4[0]+(q4[1]+(q4[2]+(q4[3]+q4[4]*u)*u)*u)*u);
529  } else if(v < 50) {
530  u = 1/v;
531  denlan = u*u*(p5[0]+(p5[1]+(p5[2]+(p5[3]+p5[4]*u)*u)*u)*u)/
532  (q5[0]+(q5[1]+(q5[2]+(q5[3]+q5[4]*u)*u)*u)*u);
533  } else if(v < 300) {
534  u = 1/v;
535  denlan = u*u*(p6[0]+(p6[1]+(p6[2]+(p6[3]+p6[4]*u)*u)*u)*u)/
536  (q6[0]+(q6[1]+(q6[2]+(q6[3]+q6[4]*u)*u)*u)*u);
537  } else {
538  u = 1/(v-v*Math.log(v)/(v+1));
539  denlan = u*u*(1+(a2[0]+a2[1]*u)*u);
540  }
541  return denlan/xi;
542  };
543 
544  JSROOT.Math.Landau = function(x, mpv, sigma, norm) {
545  if (sigma <= 0) return 0;
546  var den = JSROOT.Math.landau_pdf((x - mpv) / sigma, 1, 0);
547  if (!norm) return den;
548  return den/sigma;
549  }
550 
551  JSROOT.Math.inc_gamma_c = function(a,x) {
552  return JSROOT.Math.igamc(a,x);
553  }
554 
555  JSROOT.Math.chisquared_cdf_c = function(x,r,x0) {
556  return JSROOT.Math.inc_gamma_c ( 0.5 * r , 0.5* (x-x0) );
557  }
558 
559  JSROOT.Math.Prob = function(chi2, ndf) {
560  if (ndf <= 0) return 0; // Set CL to zero in case ndf<=0
561 
562  if (chi2 <= 0) {
563  if (chi2 < 0) return 0;
564  else return 1;
565  }
566 
567  return JSROOT.Math.chisquared_cdf_c(chi2,ndf,0);
568  }
569 
570  JSROOT.Math.gaus = function(f, x, i) {
571  return f['fParams'][i+0] * Math.exp(-0.5 * Math.pow((x-f['fParams'][i+1]) / f['fParams'][i+2], 2));
572  };
573 
574  JSROOT.Math.gausn = function(f, x, i) {
575  return JSROOT.Math.gaus(f, x, i)/(Math.sqrt(2 * Math.PI) * f['fParams'][i+2]);
576  };
577 
578  JSROOT.Math.expo = function(f, x, i) {
579  return Math.exp(f['fParams'][i+0] + f['fParams'][i+1] * x);
580  };
581 
582  JSROOT.Math.landau = function(f, x, i) {
583  return JSROOT.Math.Landau(x, f['fParams'][i+1],f['fParams'][i+2], false);
584  }
585 
586  JSROOT.Math.landaun = function(f, x, i) {
587  return JSROOT.Math.Landau(x, f['fParams'][i+1],f['fParams'][i+2], true);
588  }
589 
590  return JSROOT;
591 
592 }));