123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566 |
- /* randist/landau.c
- *
- * Copyright (C) 2001, 2004 David Morrison
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 3 of the License, or (at
- * your option) any later version.
- *
- * This program is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- * General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- */
- /* Adapted from the CERN library routines DENLAN, RANLAN, and DISLAN
- * as described in http://consult.cern.ch/shortwrups/g110/top.html.
- * Original author: K.S. K\"olbig.
- *
- * The distribution is given by the complex path integral,
- *
- * p(x) = (1/(2 pi i)) \int_{c-i\inf}^{c+i\inf} ds exp(s log(s) + x s)
- *
- * which can be converted into a real integral over [0,+\inf]
- *
- * p(x) = (1/pi) \int_0^\inf dt \exp(-t log(t) - x t) sin(pi t)
- *
- */
- #include "gsl__config.h"
- #include <math.h>
- #include "gsl_rng.h"
- #include "gsl_randist.h"
- double
- gsl_ran_landau_pdf(const double x)
- {
- static double P1[5] =
- {
- 0.4259894875E0, -0.1249762550E0, 0.3984243700E-1,
- -0.6298287635E-2, 0.1511162253E-2
- };
- static double P2[5] =
- {
- 0.1788541609E0, 0.1173957403E0, 0.1488850518E-1,
- -0.1394989411E-2, 0.1283617211E-3
- };
- static double P3[5] =
- {
- 0.1788544503E0, 0.9359161662E-1, 0.6325387654E-2,
- 0.6611667319E-4, -0.2031049101E-5
- };
- static double P4[5] =
- {
- 0.9874054407E0, 0.1186723273E3, 0.8492794360E3,
- -0.7437792444E3, 0.4270262186E3
- };
- static double P5[5] =
- {
- 0.1003675074E1, 0.1675702434E3, 0.4789711289E4,
- 0.2121786767E5, -0.2232494910E5
- };
- static double P6[5] =
- {
- 0.1000827619E1, 0.6649143136E3, 0.6297292665E5,
- 0.4755546998E6, -0.5743609109E7
- };
- static double Q1[5] =
- {
- 1.0, -0.3388260629E0, 0.9594393323E-1,
- -0.1608042283E-1, 0.3778942063E-2
- };
- static double Q2[5] =
- {
- 1.0, 0.7428795082E0, 0.3153932961E0,
- 0.6694219548E-1, 0.8790609714E-2
- };
- static double Q3[5] =
- {
- 1.0, 0.6097809921E0, 0.2560616665E0,
- 0.4746722384E-1, 0.6957301675E-2
- };
- static double Q4[5] =
- {
- 1.0, 0.1068615961E3, 0.3376496214E3,
- 0.2016712389E4, 0.1597063511E4
- };
- static double Q5[5] =
- {
- 1.0, 0.1569424537E3, 0.3745310488E4,
- 0.9834698876E4, 0.6692428357E5
- };
- static double Q6[5] =
- {
- 1.0, 0.6514101098E3, 0.5697473333E5,
- 0.1659174725E6, -0.2815759939E7
- };
- static double A1[3] =
- {
- 0.4166666667E-1, -0.1996527778E-1, 0.2709538966E-1
- };
- static double A2[2] =
- {
- -0.1845568670E1, -0.4284640743E1
- };
- double U, V, DENLAN;
- V = x;
- if (V < -5.5)
- {
- U = exp(V + 1.0);
- DENLAN = 0.3989422803 * (exp( -1 / U) / sqrt(U)) *
- (1 + (A1[0] + (A1[1] + A1[2] * U) * U) * U);
- }
- else if (V < -1)
- {
- U = exp( -V - 1);
- DENLAN = exp( -U) * sqrt(U) *
- (P1[0] + (P1[1] + (P1[2] + (P1[3] + P1[4] * V) * V) * V) * V) /
- (Q1[0] + (Q1[1] + (Q1[2] + (Q1[3] + Q1[4] * V) * V) * V) * V);
- }
- else if (V < 1)
- {
- DENLAN = (P2[0] + (P2[1] + (P2[2] + (P2[3] + P2[4] * V) * V) * V) * V) /
- (Q2[0] + (Q2[1] + (Q2[2] + (Q2[3] + Q2[4] * V) * V) * V) * V);
- }
- else if (V < 5)
- {
- DENLAN = (P3[0] + (P3[1] + (P3[2] + (P3[3] + P3[4] * V) * V) * V) * V) /
- (Q3[0] + (Q3[1] + (Q3[2] + (Q3[3] + Q3[4] * V) * V) * V) * V);
- }
- else if (V < 12)
- {
- U = 1 / V;
- DENLAN = U * U *
- (P4[0] + (P4[1] + (P4[2] + (P4[3] + P4[4] * U) * U) * U) * U) /
- (Q4[0] + (Q4[1] + (Q4[2] + (Q4[3] + Q4[4] * U) * U) * U) * U);
- }
- else if (V < 50)
- {
- U = 1 / V;
- DENLAN = U * U *
- (P5[0] + (P5[1] + (P5[2] + (P5[3] + P5[4] * U) * U) * U) * U) /
- (Q5[0] + (Q5[1] + (Q5[2] + (Q5[3] + Q5[4] * U) * U) * U) * U);
- }
- else if (V < 300)
- {
- U = 1 / V;
- DENLAN = U * U *
- (P6[0] + (P6[1] + (P6[2] + (P6[3] + P6[4] * U) * U) * U) * U) /
- (Q6[0] + (Q6[1] + (Q6[2] + (Q6[3] + Q6[4] * U) * U) * U) * U);
- }
- else
- {
- U = 1 / (V - V * log(V) / (V + 1));
- DENLAN = U * U * (1 + (A2[0] + A2[1] * U) * U);
- }
- return DENLAN;
- }
- #if 0 /* Not needed yet */
- /* This function is a translation from the original Fortran of the
- * CERN library routine DISLAN, the integral from -inf to x of the
- * Landau p.d.f.
- */
- static
- double
- gsl_ran_landau_dislan(const double x)
- {
- static double P1[5] =
- {
- 0.2514091491E0, -0.6250580444E-1,
- 0.1458381230E-1, -0.2108817737E-2,
- 0.7411247290E-3
- };
- static double P2[4] =
- {
- 0.2868328584E0, 0.3564363231E0,
- 0.1523518695E0, 0.2251304883E-1
- };
- static double P3[4] =
- {
- 0.2868329066E0, 0.3003828436E0,
- 0.9950951941E-1, 0.8733827185E-2
- };
- static double P4[4] =
- {
- 0.1000351630E1, 0.4503592498E1,
- 0.1085883880E2, 0.7536052269E1
- };
- static double P5[4] =
- {
- 0.1000006517E1, 0.4909414111E2,
- 0.8505544753E2, 0.1532153455E3
- };
- static double P6[4] =
- {
- 0.1000000983E1, 0.1329868456E3,
- 0.9162149244E3, -0.9605054274E3
- };
- static double Q1[5] =
- {
- 1.0, -0.5571175625E-2,
- 0.6225310236E-1, -0.3137378427E-2,
- 0.1931496439E-2
- };
- static double Q2[4] =
- {
- 1.0, 0.6191136137E0,
- 0.1720721448E0, 0.2278594771E-1
- };
- static double Q3[4] =
- {
- 1.0, 0.4237190502E0,
- 0.1095631512E0, 0.8693851567E-2
- };
- static double Q4[4] =
- {
- 1.0, 0.5539969678E1,
- 0.1933581111E2, 0.2721321508E2
- };
- static double Q5[4] =
- {
- 1.0, 0.5009928881E2,
- 0.1399819104E3, 0.4200002909E3
- };
- static double Q6[4] =
- {
- 1.0, 0.1339887843E3,
- 0.1055990413E4, 0.5532224619E3
- };
- static double A1[3] =
- {
- -0.4583333333E0, 0.6675347222E0, -0.1641741416E1
- };
- static double A2[3] =
- {
- 1.0, -0.4227843351E0, -0.2043403138E1
- };
- double U, V, DISLAN;
- V = x;
- if (V < -5.5)
- {
- U = exp(V + 1);
- DISLAN = 0.3989422803 * exp( -1 / U) * sqrt(U) *
- (1 + (A1[0] + (A1[1] + A1[2] * U) * U) * U);
- }
- else if (V < -1)
- {
- U = exp( -V - 1);
- DISLAN = (exp( -U) / sqrt(U)) *
- (P1[0] + (P1[1] + (P1[2] + (P1[3] + P1[4] * V) * V) * V) * V) /
- (Q1[0] + (Q1[1] + (Q1[2] + (Q1[3] + Q1[4] * V) * V) * V) * V);
- }
- else if (V < 1)
- {
- DISLAN = (P2[0] + (P2[1] + (P2[2] + P2[3] * V) * V) * V) /
- (Q2[0] + (Q2[1] + (Q2[2] + Q2[3] * V) * V) * V);
- }
- else if (V < 4)
- {
- DISLAN = (P3[0] + (P3[1] + (P3[2] + P3[3] * V) * V) * V) /
- (Q3[0] + (Q3[1] + (Q3[2] + Q3[3] * V) * V) * V);
- }
- else if (V < 12)
- {
- U = 1 / V;
- DISLAN = (P4[0] + (P4[1] + (P4[2] + P4[3] * U) * U) * U) /
- (Q4[0] + (Q4[1] + (Q4[2] + Q4[3] * U) * U) * U);
- }
- else if (V < 50)
- {
- U = 1 / V;
- DISLAN = (P5[0] + (P5[1] + (P5[2] + P5[3] * U) * U) * U) /
- (Q5[0] + (Q5[1] + (Q5[2] + Q5[3] * U) * U) * U);
- }
- else if (V < 300)
- {
- U = 1 / V;
- DISLAN = (P6[0] + (P6[1] + (P6[2] + P6[3] * U) * U) * U) /
- (Q6[0] + (Q6[1] + (Q6[2] + Q6[3] * U) * U) * U);
- }
- else
- {
- U = 1 / (V - V * log(V) / (V + 1));
- DISLAN = 1 - (A2[0] + (A2[1] + A2[2] * U) * U) * U;
- }
- return DISLAN;
- }
- #endif
- double
- gsl_ran_landau(const gsl_rng * r)
- {
- static double F[983] =
- {
- 0.0000000, /* Add empty element [0] to account for difference
- between C and Fortran convention for lower bound. */
- 00.000000, 00.000000, 00.000000, 00.000000, 00.000000,
- -2.244733, -2.204365, -2.168163, -2.135219, -2.104898,
- -2.076740, -2.050397, -2.025605, -2.002150, -1.979866,
- -1.958612, -1.938275, -1.918760, -1.899984, -1.881879,
- -1.864385, -1.847451, -1.831030, -1.815083, -1.799574,
- -1.784473, -1.769751, -1.755383, -1.741346, -1.727620,
- -1.714187, -1.701029, -1.688130, -1.675477, -1.663057,
- -1.650858, -1.638868, -1.627078, -1.615477, -1.604058,
- -1.592811, -1.581729, -1.570806, -1.560034, -1.549407,
- -1.538919, -1.528565, -1.518339, -1.508237, -1.498254,
- -1.488386, -1.478628, -1.468976, -1.459428, -1.449979,
- -1.440626, -1.431365, -1.422195, -1.413111, -1.404112,
- -1.395194, -1.386356, -1.377594, -1.368906, -1.360291,
- -1.351746, -1.343269, -1.334859, -1.326512, -1.318229,
- -1.310006, -1.301843, -1.293737, -1.285688, -1.277693,
- -1.269752, -1.261863, -1.254024, -1.246235, -1.238494,
- -1.230800, -1.223153, -1.215550, -1.207990, -1.200474,
- -1.192999, -1.185566, -1.178172, -1.170817, -1.163500,
- -1.156220, -1.148977, -1.141770, -1.134598, -1.127459,
- -1.120354, -1.113282, -1.106242, -1.099233, -1.092255,
- -1.085306, -1.078388, -1.071498, -1.064636, -1.057802,
- -1.050996, -1.044215, -1.037461, -1.030733, -1.024029,
- -1.017350, -1.010695, -1.004064, -0.997456, -0.990871,
- -0.984308, -0.977767, -0.971247, -0.964749, -0.958271,
- -0.951813, -0.945375, -0.938957, -0.932558, -0.926178,
- -0.919816, -0.913472, -0.907146, -0.900838, -0.894547,
- -0.888272, -0.882014, -0.875773, -0.869547, -0.863337,
- -0.857142, -0.850963, -0.844798, -0.838648, -0.832512,
- -0.826390, -0.820282, -0.814187, -0.808106, -0.802038,
- -0.795982, -0.789940, -0.783909, -0.777891, -0.771884,
- -0.765889, -0.759906, -0.753934, -0.747973, -0.742023,
- -0.736084, -0.730155, -0.724237, -0.718328, -0.712429,
- -0.706541, -0.700661, -0.694791, -0.688931, -0.683079,
- -0.677236, -0.671402, -0.665576, -0.659759, -0.653950,
- -0.648149, -0.642356, -0.636570, -0.630793, -0.625022,
- -0.619259, -0.613503, -0.607754, -0.602012, -0.596276,
- -0.590548, -0.584825, -0.579109, -0.573399, -0.567695,
- -0.561997, -0.556305, -0.550618, -0.544937, -0.539262,
- -0.533592, -0.527926, -0.522266, -0.516611, -0.510961,
- -0.505315, -0.499674, -0.494037, -0.488405, -0.482777,
- -0.477153, -0.471533, -0.465917, -0.460305, -0.454697,
- -0.449092, -0.443491, -0.437893, -0.432299, -0.426707,
- -0.421119, -0.415534, -0.409951, -0.404372, -0.398795,
- -0.393221, -0.387649, -0.382080, -0.376513, -0.370949,
- -0.365387, -0.359826, -0.354268, -0.348712, -0.343157,
- -0.337604, -0.332053, -0.326503, -0.320955, -0.315408,
- -0.309863, -0.304318, -0.298775, -0.293233, -0.287692,
- -0.282152, -0.276613, -0.271074, -0.265536, -0.259999,
- -0.254462, -0.248926, -0.243389, -0.237854, -0.232318,
- -0.226783, -0.221247, -0.215712, -0.210176, -0.204641,
- -0.199105, -0.193568, -0.188032, -0.182495, -0.176957,
- -0.171419, -0.165880, -0.160341, -0.154800, -0.149259,
- -0.143717, -0.138173, -0.132629, -0.127083, -0.121537,
- -0.115989, -0.110439, -0.104889, -0.099336, -0.093782,
- -0.088227, -0.082670, -0.077111, -0.071550, -0.065987,
- -0.060423, -0.054856, -0.049288, -0.043717, -0.038144,
- -0.032569, -0.026991, -0.021411, -0.015828, -0.010243,
- -0.004656, 00.000934, 00.006527, 00.012123, 00.017722,
- 00.023323, 00.028928, 00.034535, 00.040146, 00.045759,
- 00.051376, 00.056997, 00.062620, 00.068247, 00.073877,
- 00.079511, 00.085149, 00.090790, 00.096435, 00.102083,
- 00.107736, 00.113392, 00.119052, 00.124716, 00.130385,
- 00.136057, 00.141734, 00.147414, 00.153100, 00.158789,
- 00.164483, 00.170181, 00.175884, 00.181592, 00.187304,
- 00.193021, 00.198743, 00.204469, 00.210201, 00.215937,
- 00.221678, 00.227425, 00.233177, 00.238933, 00.244696,
- 00.250463, 00.256236, 00.262014, 00.267798, 00.273587,
- 00.279382, 00.285183, 00.290989, 00.296801, 00.302619,
- 00.308443, 00.314273, 00.320109, 00.325951, 00.331799,
- 00.337654, 00.343515, 00.349382, 00.355255, 00.361135,
- 00.367022, 00.372915, 00.378815, 00.384721, 00.390634,
- 00.396554, 00.402481, 00.408415, 00.414356, 00.420304,
- 00.426260, 00.432222, 00.438192, 00.444169, 00.450153,
- 00.456145, 00.462144, 00.468151, 00.474166, 00.480188,
- 00.486218, 00.492256, 00.498302, 00.504356, 00.510418,
- 00.516488, 00.522566, 00.528653, 00.534747, 00.540850,
- 00.546962, 00.553082, 00.559210, 00.565347, 00.571493,
- 00.577648, 00.583811, 00.589983, 00.596164, 00.602355,
- 00.608554, 00.614762, 00.620980, 00.627207, 00.633444,
- 00.639689, 00.645945, 00.652210, 00.658484, 00.664768,
- 00.671062, 00.677366, 00.683680, 00.690004, 00.696338,
- 00.702682, 00.709036, 00.715400, 00.721775, 00.728160,
- 00.734556, 00.740963, 00.747379, 00.753807, 00.760246,
- 00.766695, 00.773155, 00.779627, 00.786109, 00.792603,
- 00.799107, 00.805624, 00.812151, 00.818690, 00.825241,
- 00.831803, 00.838377, 00.844962, 00.851560, 00.858170,
- 00.864791, 00.871425, 00.878071, 00.884729, 00.891399,
- 00.898082, 00.904778, 00.911486, 00.918206, 00.924940,
- 00.931686, 00.938446, 00.945218, 00.952003, 00.958802,
- 00.965614, 00.972439, 00.979278, 00.986130, 00.992996,
- 00.999875, 01.006769, 01.013676, 01.020597, 01.027533,
- 01.034482, 01.041446, 01.048424, 01.055417, 01.062424,
- 01.069446, 01.076482, 01.083534, 01.090600, 01.097681,
- 01.104778, 01.111889, 01.119016, 01.126159, 01.133316,
- 01.140490, 01.147679, 01.154884, 01.162105, 01.169342,
- 01.176595, 01.183864, 01.191149, 01.198451, 01.205770,
- 01.213105, 01.220457, 01.227826, 01.235211, 01.242614,
- 01.250034, 01.257471, 01.264926, 01.272398, 01.279888,
- 01.287395, 01.294921, 01.302464, 01.310026, 01.317605,
- 01.325203, 01.332819, 01.340454, 01.348108, 01.355780,
- 01.363472, 01.371182, 01.378912, 01.386660, 01.394429,
- 01.402216, 01.410024, 01.417851, 01.425698, 01.433565,
- 01.441453, 01.449360, 01.457288, 01.465237, 01.473206,
- 01.481196, 01.489208, 01.497240, 01.505293, 01.513368,
- 01.521465, 01.529583, 01.537723, 01.545885, 01.554068,
- 01.562275, 01.570503, 01.578754, 01.587028, 01.595325,
- 01.603644, 01.611987, 01.620353, 01.628743, 01.637156,
- 01.645593, 01.654053, 01.662538, 01.671047, 01.679581,
- 01.688139, 01.696721, 01.705329, 01.713961, 01.722619,
- 01.731303, 01.740011, 01.748746, 01.757506, 01.766293,
- 01.775106, 01.783945, 01.792810, 01.801703, 01.810623,
- 01.819569, 01.828543, 01.837545, 01.846574, 01.855631,
- 01.864717, 01.873830, 01.882972, 01.892143, 01.901343,
- 01.910572, 01.919830, 01.929117, 01.938434, 01.947781,
- 01.957158, 01.966566, 01.976004, 01.985473, 01.994972,
- 02.004503, 02.014065, 02.023659, 02.033285, 02.042943,
- 02.052633, 02.062355, 02.072110, 02.081899, 02.091720,
- 02.101575, 02.111464, 02.121386, 02.131343, 02.141334,
- 02.151360, 02.161421, 02.171517, 02.181648, 02.191815,
- 02.202018, 02.212257, 02.222533, 02.232845, 02.243195,
- 02.253582, 02.264006, 02.274468, 02.284968, 02.295507,
- 02.306084, 02.316701, 02.327356, 02.338051, 02.348786,
- 02.359562, 02.370377, 02.381234, 02.392131, 02.403070,
- 02.414051, 02.425073, 02.436138, 02.447246, 02.458397,
- 02.469591, 02.480828, 02.492110, 02.503436, 02.514807,
- 02.526222, 02.537684, 02.549190, 02.560743, 02.572343,
- 02.583989, 02.595682, 02.607423, 02.619212, 02.631050,
- 02.642936, 02.654871, 02.666855, 02.678890, 02.690975,
- 02.703110, 02.715297, 02.727535, 02.739825, 02.752168,
- 02.764563, 02.777012, 02.789514, 02.802070, 02.814681,
- 02.827347, 02.840069, 02.852846, 02.865680, 02.878570,
- 02.891518, 02.904524, 02.917588, 02.930712, 02.943894,
- 02.957136, 02.970439, 02.983802, 02.997227, 03.010714,
- 03.024263, 03.037875, 03.051551, 03.065290, 03.079095,
- 03.092965, 03.106900, 03.120902, 03.134971, 03.149107,
- 03.163312, 03.177585, 03.191928, 03.206340, 03.220824,
- 03.235378, 03.250005, 03.264704, 03.279477, 03.294323,
- 03.309244, 03.324240, 03.339312, 03.354461, 03.369687,
- 03.384992, 03.400375, 03.415838, 03.431381, 03.447005,
- 03.462711, 03.478500, 03.494372, 03.510328, 03.526370,
- 03.542497, 03.558711, 03.575012, 03.591402, 03.607881,
- 03.624450, 03.641111, 03.657863, 03.674708, 03.691646,
- 03.708680, 03.725809, 03.743034, 03.760357, 03.777779,
- 03.795300, 03.812921, 03.830645, 03.848470, 03.866400,
- 03.884434, 03.902574, 03.920821, 03.939176, 03.957640,
- 03.976215, 03.994901, 04.013699, 04.032612, 04.051639,
- 04.070783, 04.090045, 04.109425, 04.128925, 04.148547,
- 04.168292, 04.188160, 04.208154, 04.228275, 04.248524,
- 04.268903, 04.289413, 04.310056, 04.330832, 04.351745,
- 04.372794, 04.393982, 04.415310, 04.436781, 04.458395,
- 04.480154, 04.502060, 04.524114, 04.546319, 04.568676,
- 04.591187, 04.613854, 04.636678, 04.659662, 04.682807,
- 04.706116, 04.729590, 04.753231, 04.777041, 04.801024,
- 04.825179, 04.849511, 04.874020, 04.898710, 04.923582,
- 04.948639, 04.973883, 04.999316, 05.024942, 05.050761,
- 05.076778, 05.102993, 05.129411, 05.156034, 05.182864,
- 05.209903, 05.237156, 05.264625, 05.292312, 05.320220,
- 05.348354, 05.376714, 05.405306, 05.434131, 05.463193,
- 05.492496, 05.522042, 05.551836, 05.581880, 05.612178,
- 05.642734, 05.673552, 05.704634, 05.735986, 05.767610,
- 05.799512, 05.831694, 05.864161, 05.896918, 05.929968,
- 05.963316, 05.996967, 06.030925, 06.065194, 06.099780,
- 06.134687, 06.169921, 06.205486, 06.241387, 06.277630,
- 06.314220, 06.351163, 06.388465, 06.426130, 06.464166,
- 06.502578, 06.541371, 06.580553, 06.620130, 06.660109,
- 06.700495, 06.741297, 06.782520, 06.824173, 06.866262,
- 06.908795, 06.951780, 06.995225, 07.039137, 07.083525,
- 07.128398, 07.173764, 07.219632, 07.266011, 07.312910,
- 07.360339, 07.408308, 07.456827, 07.505905, 07.555554,
- 07.605785, 07.656608, 07.708035, 07.760077, 07.812747,
- 07.866057, 07.920019, 07.974647, 08.029953, 08.085952,
- 08.142657, 08.200083, 08.258245, 08.317158, 08.376837,
- 08.437300, 08.498562, 08.560641, 08.623554, 08.687319,
- 08.751955, 08.817481, 08.883916, 08.951282, 09.019600,
- 09.088889, 09.159174, 09.230477, 09.302822, 09.376233,
- 09.450735, 09.526355, 09.603118, 09.681054, 09.760191,
- 09.840558, 09.922186, 10.005107, 10.089353, 10.174959,
- 10.261958, 10.350389, 10.440287, 10.531693, 10.624646,
- 10.719188, 10.815362, 10.913214, 11.012789, 11.114137,
- 11.217307, 11.322352, 11.429325, 11.538283, 11.649285,
- 11.762390, 11.877664, 11.995170, 12.114979, 12.237161,
- 12.361791, 12.488946, 12.618708, 12.751161, 12.886394,
- 13.024498, 13.165570, 13.309711, 13.457026, 13.607625,
- 13.761625, 13.919145, 14.080314, 14.245263, 14.414134,
- 14.587072, 14.764233, 14.945778, 15.131877, 15.322712,
- 15.518470, 15.719353, 15.925570, 16.137345, 16.354912,
- 16.578520, 16.808433, 17.044929, 17.288305, 17.538873,
- 17.796967, 18.062943, 18.337176, 18.620068, 18.912049,
- 19.213574, 19.525133, 19.847249, 20.180480, 20.525429,
- 20.882738, 21.253102, 21.637266, 22.036036, 22.450278,
- 22.880933, 23.329017, 23.795634, 24.281981, 24.789364,
- 25.319207, 25.873062, 26.452634, 27.059789, 27.696581,
- 28.365274, 29.068370, 29.808638, 30.589157, 31.413354,
- 32.285060, 33.208568, 34.188705, 35.230920, 36.341388,
- 37.527131, 38.796172, 40.157721, 41.622399, 43.202525,
- 44.912465, 46.769077, 48.792279, 51.005773, 53.437996,
- 56.123356, 59.103894
- };
- double X, U, V, RANLAN;
- int I;
- X = gsl_rng_uniform_pos(r);
- U = 1000.0 * X;
- I = U;
- U = U - I;
- if (I >= 70 && I <= 800)
- {
- RANLAN = F[I] + U * (F[I + 1] - F[I]);
- }
- else if (I >= 7 && I <= 980)
- {
- RANLAN = F[I]
- + U * (F[I + 1] - F[I]
- - 0.25 * (1 - U) * (F[I + 2] - F[I + 1] - F[I] + F[I - 1]));
- }
- else if (I < 7)
- {
- V = log(X);
- U = 1 / V;
- RANLAN = ((0.99858950 + (3.45213058E1 + 1.70854528E1 * U) * U) /
- (1 + (3.41760202E1 + 4.01244582 * U) * U)) *
- ( -log( -0.91893853 - V) - 1);
- }
- else
- {
- U = 1 - X;
- V = U * U;
- if (X <= 0.999)
- {
- RANLAN = (1.00060006 + 2.63991156E2 * U + 4.37320068E3 * V) /
- ((1 + 2.57368075E2 * U + 3.41448018E3 * V) * U);
- }
- else
- {
- RANLAN = (1.00001538 + 6.07514119E3 * U + 7.34266409E5 * V) /
- ((1 + 6.06511919E3 * U + 6.94021044E5 * V) * U);
- }
- }
- return RANLAN;
- }
|