/* 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 <config.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/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;
}