Kernel/detail/LandauQuantile.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2018 CERN for the benefit of the Acts project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include <cmath>
0012 #include <limits>
0013
0014 static inline double landau_quantile(double z, double xi) {
0015   // LANDAU quantile : algorithm from CERNLIB G110 ranlan
0016   // with scale parameter xi
0017   // Converted by Rene Brun from CERNLIB routine ranlan(G110),
0018   // Moved and adapted to QuantFuncMathCore by B. List 29.4.2010
0019
0020   static double f[982] = {
0021       0,         0,         0,         0,         0,         -2.244733,
0022       -2.204365, -2.168163, -2.135219, -2.104898, -2.076740, -2.050397,
0023       -2.025605, -2.002150, -1.979866, -1.958612, -1.938275, -1.918760,
0024       -1.899984, -1.881879, -1.864385, -1.847451, -1.831030, -1.815083,
0025       -1.799574, -1.784473, -1.769751, -1.755383, -1.741346, -1.727620,
0026       -1.714187, -1.701029, -1.688130, -1.675477, -1.663057, -1.650858,
0027       -1.638868, -1.627078, -1.615477, -1.604058, -1.592811, -1.581729,
0028       -1.570806, -1.560034, -1.549407, -1.538919, -1.528565, -1.518339,
0029       -1.508237, -1.498254, -1.488386, -1.478628, -1.468976, -1.459428,
0030       -1.449979, -1.440626, -1.431365, -1.422195, -1.413111, -1.404112,
0031       -1.395194, -1.386356, -1.377594, -1.368906, -1.360291, -1.351746,
0032       -1.343269, -1.334859, -1.326512, -1.318229, -1.310006, -1.301843,
0033       -1.293737, -1.285688, -1.277693, -1.269752, -1.261863, -1.254024,
0034       -1.246235, -1.238494, -1.230800, -1.223153, -1.215550, -1.207990,
0035       -1.200474, -1.192999, -1.185566, -1.178172, -1.170817, -1.163500,
0036       -1.156220, -1.148977, -1.141770, -1.134598, -1.127459, -1.120354,
0037       -1.113282, -1.106242, -1.099233, -1.092255, -1.085306, -1.078388,
0038       -1.071498, -1.064636, -1.057802, -1.050996, -1.044215, -1.037461,
0039       -1.030733, -1.024029, -1.017350, -1.010695, -1.004064, -.997456,
0040       -.990871,  -.984308,  -.977767,  -.971247,  -.964749,  -.958271,
0041       -.951813,  -.945375,  -.938957,  -.932558,  -.926178,  -.919816,
0042       -.913472,  -.907146,  -.900838,  -.894547,  -.888272,  -.882014,
0043       -.875773,  -.869547,  -.863337,  -.857142,  -.850963,  -.844798,
0044       -.838648,  -.832512,  -.826390,  -.820282,  -.814187,  -.808106,
0045       -.802038,  -.795982,  -.789940,  -.783909,  -.777891,  -.771884,
0046       -.765889,  -.759906,  -.753934,  -.747973,  -.742023,  -.736084,
0047       -.730155,  -.724237,  -.718328,  -.712429,  -.706541,  -.700661,
0048       -.694791,  -.688931,  -.683079,  -.677236,  -.671402,  -.665576,
0049       -.659759,  -.653950,  -.648149,  -.642356,  -.636570,  -.630793,
0050       -.625022,  -.619259,  -.613503,  -.607754,  -.602012,  -.596276,
0051       -.590548,  -.584825,  -.579109,  -.573399,  -.567695,  -.561997,
0052       -.556305,  -.550618,  -.544937,  -.539262,  -.533592,  -.527926,
0053       -.522266,  -.516611,  -.510961,  -.505315,  -.499674,  -.494037,
0054       -.488405,  -.482777,  -.477153,  -.471533,  -.465917,  -.460305,
0055       -.454697,  -.449092,  -.443491,  -.437893,  -.432299,  -.426707,
0056       -.421119,  -.415534,  -.409951,  -.404372,  -.398795,  -.393221,
0057       -.387649,  -.382080,  -.376513,  -.370949,  -.365387,  -.359826,
0058       -.354268,  -.348712,  -.343157,  -.337604,  -.332053,  -.326503,
0059       -.320955,  -.315408,  -.309863,  -.304318,  -.298775,  -.293233,
0060       -.287692,  -.282152,  -.276613,  -.271074,  -.265536,  -.259999,
0061       -.254462,  -.248926,  -.243389,  -.237854,  -.232318,  -.226783,
0062       -.221247,  -.215712,  -.210176,  -.204641,  -.199105,  -.193568,
0063       -.188032,  -.182495,  -.176957,  -.171419,  -.165880,  -.160341,
0064       -.154800,  -.149259,  -.143717,  -.138173,  -.132629,  -.127083,
0065       -.121537,  -.115989,  -.110439,  -.104889,  -.099336,  -.093782,
0066       -.088227,  -.082670,  -.077111,  -.071550,  -.065987,  -.060423,
0067       -.054856,  -.049288,  -.043717,  -.038144,  -.032569,  -.026991,
0068       -.021411,  -.015828,  -.010243,  -.004656,  .000934,   .006527,
0069       .012123,   .017722,   .023323,   .028928,   .034535,   .040146,
0070       .045759,   .051376,   .056997,   .062620,   .068247,   .073877,
0071       .079511,   .085149,   .090790,   .096435,   .102083,   .107736,
0072       .113392,   .119052,   .124716,   .130385,   .136057,   .141734,
0073       .147414,   .153100,   .158789,   .164483,   .170181,   .175884,
0074       .181592,   .187304,   .193021,   .198743,   .204469,   .210201,
0075       .215937,   .221678,   .227425,   .233177,   .238933,   .244696,
0076       .250463,   .256236,   .262014,   .267798,   .273587,   .279382,
0077       .285183,   .290989,   .296801,   .302619,   .308443,   .314273,
0078       .320109,   .325951,   .331799,   .337654,   .343515,   .349382,
0079       .355255,   .361135,   .367022,   .372915,   .378815,   .384721,
0080       .390634,   .396554,   .402481,   .408415,   .414356,   .420304,
0081       .426260,   .432222,   .438192,   .444169,   .450153,   .456145,
0082       .462144,   .468151,   .474166,   .480188,   .486218,   .492256,
0083       .498302,   .504356,   .510418,   .516488,   .522566,   .528653,
0084       .534747,   .540850,   .546962,   .553082,   .559210,   .565347,
0085       .571493,   .577648,   .583811,   .589983,   .596164,   .602355,
0086       .608554,   .614762,   .620980,   .627207,   .633444,   .639689,
0087       .645945,   .652210,   .658484,   .664768,   .671062,   .677366,
0088       .683680,   .690004,   .696338,   .702682,   .709036,   .715400,
0089       .721775,   .728160,   .734556,   .740963,   .747379,   .753807,
0090       .760246,   .766695,   .773155,   .779627,   .786109,   .792603,
0091       .799107,   .805624,   .812151,   .818690,   .825241,   .831803,
0092       .838377,   .844962,   .851560,   .858170,   .864791,   .871425,
0093       .878071,   .884729,   .891399,   .898082,   .904778,   .911486,
0094       .918206,   .924940,   .931686,   .938446,   .945218,   .952003,
0095       .958802,   .965614,   .972439,   .979278,   .986130,   .992996,
0096       .999875,   1.006769,  1.013676,  1.020597,  1.027533,  1.034482,
0097       1.041446,  1.048424,  1.055417,  1.062424,  1.069446,  1.076482,
0098       1.083534,  1.090600,  1.097681,  1.104778,  1.111889,  1.119016,
0099       1.126159,  1.133316,  1.140490,  1.147679,  1.154884,  1.162105,
0100       1.169342,  1.176595,  1.183864,  1.191149,  1.198451,  1.205770,
0101       1.213105,  1.220457,  1.227826,  1.235211,  1.242614,  1.250034,
0102       1.257471,  1.264926,  1.272398,  1.279888,  1.287395,  1.294921,
0103       1.302464,  1.310026,  1.317605,  1.325203,  1.332819,  1.340454,
0104       1.348108,  1.355780,  1.363472,  1.371182,  1.378912,  1.386660,
0105       1.394429,  1.402216,  1.410024,  1.417851,  1.425698,  1.433565,
0106       1.441453,  1.449360,  1.457288,  1.465237,  1.473206,  1.481196,
0107       1.489208,  1.497240,  1.505293,  1.513368,  1.521465,  1.529583,
0108       1.537723,  1.545885,  1.554068,  1.562275,  1.570503,  1.578754,
0109       1.587028,  1.595325,  1.603644,  1.611987,  1.620353,  1.628743,
0110       1.637156,  1.645593,  1.654053,  1.662538,  1.671047,  1.679581,
0111       1.688139,  1.696721,  1.705329,  1.713961,  1.722619,  1.731303,
0112       1.740011,  1.748746,  1.757506,  1.766293,  1.775106,  1.783945,
0113       1.792810,  1.801703,  1.810623,  1.819569,  1.828543,  1.837545,
0114       1.846574,  1.855631,  1.864717,  1.873830,  1.882972,  1.892143,
0115       1.901343,  1.910572,  1.919830,  1.929117,  1.938434,  1.947781,
0116       1.957158,  1.966566,  1.976004,  1.985473,  1.994972,  2.004503,
0117       2.014065,  2.023659,  2.033285,  2.042943,  2.052633,  2.062355,
0118       2.072110,  2.081899,  2.091720,  2.101575,  2.111464,  2.121386,
0119       2.131343,  2.141334,  2.151360,  2.161421,  2.171517,  2.181648,
0120       2.191815,  2.202018,  2.212257,  2.222533,  2.232845,  2.243195,
0121       2.253582,  2.264006,  2.274468,  2.284968,  2.295507,  2.306084,
0122       2.316701,  2.327356,  2.338051,  2.348786,  2.359562,  2.370377,
0123       2.381234,  2.392131,  2.403070,  2.414051,  2.425073,  2.436138,
0124       2.447246,  2.458397,  2.469591,  2.480828,  2.492110,  2.503436,
0125       2.514807,  2.526222,  2.537684,  2.549190,  2.560743,  2.572343,
0126       2.583989,  2.595682,  2.607423,  2.619212,  2.631050,  2.642936,
0127       2.654871,  2.666855,  2.678890,  2.690975,  2.703110,  2.715297,
0128       2.727535,  2.739825,  2.752168,  2.764563,  2.777012,  2.789514,
0129       2.802070,  2.814681,  2.827347,  2.840069,  2.852846,  2.865680,
0130       2.878570,  2.891518,  2.904524,  2.917588,  2.930712,  2.943894,
0131       2.957136,  2.970439,  2.983802,  2.997227,  3.010714,  3.024263,
0132       3.037875,  3.051551,  3.065290,  3.079095,  3.092965,  3.106900,
0133       3.120902,  3.134971,  3.149107,  3.163312,  3.177585,  3.191928,
0134       3.206340,  3.220824,  3.235378,  3.250005,  3.264704,  3.279477,
0135       3.294323,  3.309244,  3.324240,  3.339312,  3.354461,  3.369687,
0136       3.384992,  3.400375,  3.415838,  3.431381,  3.447005,  3.462711,
0137       3.478500,  3.494372,  3.510328,  3.526370,  3.542497,  3.558711,
0138       3.575012,  3.591402,  3.607881,  3.624450,  3.641111,  3.657863,
0139       3.674708,  3.691646,  3.708680,  3.725809,  3.743034,  3.760357,
0140       3.777779,  3.795300,  3.812921,  3.830645,  3.848470,  3.866400,
0141       3.884434,  3.902574,  3.920821,  3.939176,  3.957640,  3.976215,
0142       3.994901,  4.013699,  4.032612,  4.051639,  4.070783,  4.090045,
0143       4.109425,  4.128925,  4.148547,  4.168292,  4.188160,  4.208154,
0144       4.228275,  4.248524,  4.268903,  4.289413,  4.310056,  4.330832,
0145       4.351745,  4.372794,  4.393982,  4.415310,  4.436781,  4.458395,
0146       4.480154,  4.502060,  4.524114,  4.546319,  4.568676,  4.591187,
0147       4.613854,  4.636678,  4.659662,  4.682807,  4.706116,  4.729590,
0148       4.753231,  4.777041,  4.801024,  4.825179,  4.849511,  4.874020,
0149       4.898710,  4.923582,  4.948639,  4.973883,  4.999316,  5.024942,
0150       5.050761,  5.076778,  5.102993,  5.129411,  5.156034,  5.182864,
0151       5.209903,  5.237156,  5.264625,  5.292312,  5.320220,  5.348354,
0152       5.376714,  5.405306,  5.434131,  5.463193,  5.492496,  5.522042,
0153       5.551836,  5.581880,  5.612178,  5.642734,  5.673552,  5.704634,
0154       5.735986,  5.767610,  5.799512,  5.831694,  5.864161,  5.896918,
0155       5.929968,  5.963316,  5.996967,  6.030925,  6.065194,  6.099780,
0156       6.134687,  6.169921,  6.205486,  6.241387,  6.277630,  6.314220,
0157       6.351163,  6.388465,  6.426130,  6.464166,  6.502578,  6.541371,
0158       6.580553,  6.620130,  6.660109,  6.700495,  6.741297,  6.782520,
0159       6.824173,  6.866262,  6.908795,  6.951780,  6.995225,  7.039137,
0160       7.083525,  7.128398,  7.173764,  7.219632,  7.266011,  7.312910,
0161       7.360339,  7.408308,  7.456827,  7.505905,  7.555554,  7.605785,
0162       7.656608,  7.708035,  7.760077,  7.812747,  7.866057,  7.920019,
0163       7.974647,  8.029953,  8.085952,  8.142657,  8.200083,  8.258245,
0164       8.317158,  8.376837,  8.437300,  8.498562,  8.560641,  8.623554,
0165       8.687319,  8.751955,  8.817481,  8.883916,  8.951282,  9.019600,
0166       9.088889,  9.159174,  9.230477,  9.302822,  9.376233,  9.450735,
0167       9.526355,  9.603118,  9.681054,  9.760191,  9.840558,  9.922186,
0168       10.005107, 10.089353, 10.174959, 10.261958, 10.350389, 10.440287,
0169       10.531693, 10.624646, 10.719188, 10.815362, 10.913214, 11.012789,
0170       11.114137, 11.217307, 11.322352, 11.429325, 11.538283, 11.649285,
0171       11.762390, 11.877664, 11.995170, 12.114979, 12.237161, 12.361791,
0172       12.488946, 12.618708, 12.751161, 12.886394, 13.024498, 13.165570,
0173       13.309711, 13.457026, 13.607625, 13.761625, 13.919145, 14.080314,
0174       14.245263, 14.414134, 14.587072, 14.764233, 14.945778, 15.131877,
0175       15.322712, 15.518470, 15.719353, 15.925570, 16.137345, 16.354912,
0176       16.578520, 16.808433, 17.044929, 17.288305, 17.538873, 17.796967,
0177       18.062943, 18.337176, 18.620068, 18.912049, 19.213574, 19.525133,
0178       19.847249, 20.180480, 20.525429, 20.882738, 21.253102, 21.637266,
0179       22.036036, 22.450278, 22.880933, 23.329017, 23.795634, 24.281981,
0180       24.789364, 25.319207, 25.873062, 26.452634, 27.059789, 27.696581,
0181       28.365274, 29.068370, 29.808638, 30.589157, 31.413354, 32.285060,
0182       33.208568, 34.188705, 35.230920, 36.341388, 37.527131, 38.796172,
0183       40.157721, 41.622399, 43.202525, 44.912465, 46.769077, 48.792279,
0184       51.005773, 53.437996, 56.123356, 59.103894};
0185
0186   if (xi <= 0)
0187     return 0;
0188   if (z <= 0)
0189     return -std::numeric_limits<double>::infinity();
0190   if (z >= 1)
0191     return std::numeric_limits<double>::infinity();
0192
0193   double ranlan, u, v;
0194   u = 1000 * z;
0195   int i = int(u);
0196   u -= i;
0197   if (i >= 70 && i < 800) {
0198     ranlan = f[i - 1] + u * (f[i] - f[i - 1]);
0199   } else if (i >= 7 && i <= 980) {
0200     ranlan = f[i - 1] +
0201              u * (f[i] - f[i - 1] -
0202                   0.25 * (1 - u) * (f[i + 1] - f[i] - f[i - 1] + f[i - 2]));
0203   } else if (i < 7) {
0204     v = std::log(z);
0205     u = 1 / v;
0206     ranlan = ((0.99858950 + (3.45213058E1 + 1.70854528E1 * u) * u) /
0207               (1 + (3.41760202E1 + 4.01244582 * u) * u)) *
0208              (-std::log(-0.91893853 - v) - 1);
0209   } else {
0210     u = 1 - z;
0211     v = u * u;
0212     if (z <= 0.999) {
0213       ranlan = (1.00060006 + 2.63991156E2 * u + 4.37320068E3 * v) /
0214                ((1 + 2.57368075E2 * u + 3.41448018E3 * v) * u);
0215     } else {
0216       ranlan = (1.00001538 + 6.07514119E3 * u + 7.34266409E5 * v) /
0217                ((1 + 6.06511919E3 * u + 6.94021044E5 * v) * u);
0218     }
0219   }
0220   return xi * ranlan;
0221 }