htswanalysis/MACS/lib/gsl/gsl-1.11/randist/landau.c

   1 /* randist/landau.c
   2  *
   3  * Copyright (C) 2001, 2004 David Morrison
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Adapted from the CERN library routines DENLAN, RANLAN, and DISLAN
  21  * as described in http://consult.cern.ch/shortwrups/g110/top.html.
  22  * Original author: K.S. K\"olbig.
  23  *
  24  * The distribution is given by the complex path integral,
  25  *
  26  *  p(x) = (1/(2 pi i)) \int_{c-i\inf}^{c+i\inf} ds exp(s log(s) + x s)
  27  *
  28  * which can be converted into a real integral over [0,+\inf]
  29  *
  30  *  p(x) = (1/pi) \int_0^\inf dt \exp(-t log(t) - x t) sin(pi t)
  31  *
  32  */
  33
  34 #include <config.h>
  35 #include <math.h>
  36 #include <gsl/gsl_rng.h>
  37 #include <gsl/gsl_randist.h>
  38
  39 double
  40 gsl_ran_landau_pdf(const double x)
  41 {
  42   static double P1[5] =
  43     {
  44       0.4259894875E0, -0.1249762550E0, 0.3984243700E-1,
  45       -0.6298287635E-2, 0.1511162253E-2
  46     };
  47   static double P2[5] =
  48     {
  49       0.1788541609E0, 0.1173957403E0, 0.1488850518E-1,
  50       -0.1394989411E-2, 0.1283617211E-3
  51     };
  52   static double P3[5] =
  53     {
  54       0.1788544503E0, 0.9359161662E-1, 0.6325387654E-2,
  55       0.6611667319E-4, -0.2031049101E-5
  56     };
  57   static double P4[5] =
  58     {
  59       0.9874054407E0, 0.1186723273E3, 0.8492794360E3,
  60       -0.7437792444E3, 0.4270262186E3
  61     };
  62   static double P5[5] =
  63     {
  64       0.1003675074E1, 0.1675702434E3, 0.4789711289E4,
  65       0.2121786767E5, -0.2232494910E5
  66     };
  67   static double P6[5] =
  68     {
  69       0.1000827619E1, 0.6649143136E3, 0.6297292665E5,
  70       0.4755546998E6, -0.5743609109E7
  71     };
  72
  73   static double Q1[5] =
  74     {
  75       1.0, -0.3388260629E0, 0.9594393323E-1,
  76       -0.1608042283E-1, 0.3778942063E-2
  77     };
  78   static double Q2[5] =
  79     {
  80       1.0, 0.7428795082E0, 0.3153932961E0,
  81       0.6694219548E-1, 0.8790609714E-2
  82     };
  83   static double Q3[5] =
  84     {
  85       1.0, 0.6097809921E0, 0.2560616665E0,
  86       0.4746722384E-1, 0.6957301675E-2
  87     };
  88   static double Q4[5] =
  89     {
  90       1.0, 0.1068615961E3, 0.3376496214E3,
  91       0.2016712389E4, 0.1597063511E4
  92     };
  93   static double Q5[5] =
  94     {
  95       1.0, 0.1569424537E3, 0.3745310488E4,
  96       0.9834698876E4, 0.6692428357E5
  97     };
  98   static double Q6[5] =
  99     {
 100       1.0, 0.6514101098E3, 0.5697473333E5,
 101       0.1659174725E6, -0.2815759939E7
 102     };
 103
 104   static double A1[3] =
 105     {
 106       0.4166666667E-1, -0.1996527778E-1, 0.2709538966E-1
 107     };
 108   static double A2[2] =
 109     {
 110       -0.1845568670E1, -0.4284640743E1
 111     };
 112
 113   double U, V, DENLAN;
 114
 115   V = x;
 116   if (V < -5.5)
 117     {
 118       U = exp(V + 1.0);
 119       DENLAN = 0.3989422803 * (exp( -1 / U) / sqrt(U)) *
 120         (1 + (A1[0] + (A1[1] + A1[2] * U) * U) * U);
 121     }
 122   else if (V < -1)
 123     {
 124       U = exp( -V - 1);
 125       DENLAN = exp( -U) * sqrt(U) *
 126         (P1[0] + (P1[1] + (P1[2] + (P1[3] + P1[4] * V) * V) * V) * V) /
 127         (Q1[0] + (Q1[1] + (Q1[2] + (Q1[3] + Q1[4] * V) * V) * V) * V);
 128     }
 129   else if (V < 1)
 130     {
 131       DENLAN = (P2[0] + (P2[1] + (P2[2] + (P2[3] + P2[4] * V) * V) * V) * V) /
 132         (Q2[0] + (Q2[1] + (Q2[2] + (Q2[3] + Q2[4] * V) * V) * V) * V);
 133     }
 134   else if (V < 5)
 135     {
 136       DENLAN = (P3[0] + (P3[1] + (P3[2] + (P3[3] + P3[4] * V) * V) * V) * V) /
 137         (Q3[0] + (Q3[1] + (Q3[2] + (Q3[3] + Q3[4] * V) * V) * V) * V);
 138     }
 139   else if (V < 12)
 140     {
 141       U = 1 / V;
 142       DENLAN = U * U *
 143         (P4[0] + (P4[1] + (P4[2] + (P4[3] + P4[4] * U) * U) * U) * U) /
 144         (Q4[0] + (Q4[1] + (Q4[2] + (Q4[3] + Q4[4] * U) * U) * U) * U);
 145     }
 146   else if (V < 50)
 147     {
 148       U = 1 / V;
 149       DENLAN = U * U *
 150         (P5[0] + (P5[1] + (P5[2] + (P5[3] + P5[4] * U) * U) * U) * U) /
 151         (Q5[0] + (Q5[1] + (Q5[2] + (Q5[3] + Q5[4] * U) * U) * U) * U);
 152     }
 153   else if (V < 300)
 154     {
 155       U = 1 / V;
 156       DENLAN = U * U *
 157         (P6[0] + (P6[1] + (P6[2] + (P6[3] + P6[4] * U) * U) * U) * U) /
 158         (Q6[0] + (Q6[1] + (Q6[2] + (Q6[3] + Q6[4] * U) * U) * U) * U);
 159     }
 160   else
 161     {
 162       U = 1 / (V - V * log(V) / (V + 1));
 163       DENLAN = U * U * (1 + (A2[0] + A2[1] * U) * U);
 164     }
 165
 166   return DENLAN;
 167 }
 168
 169 #if 0 /* Not needed yet */
 170 /* This function is a translation from the original Fortran of the
 171  * CERN library routine DISLAN, the integral from -inf to x of the
 172  * Landau p.d.f.
 173  */
 174 static
 175 double
 176 gsl_ran_landau_dislan(const double x)
 177 {
 178   static double P1[5] =
 179     {
 180       0.2514091491E0, -0.6250580444E-1,
 181       0.1458381230E-1, -0.2108817737E-2,
 182       0.7411247290E-3
 183     };
 184
 185   static double P2[4] =
 186     {
 187       0.2868328584E0, 0.3564363231E0,
 188       0.1523518695E0, 0.2251304883E-1
 189     };
 190
 191   static double P3[4] =
 192     {
 193       0.2868329066E0, 0.3003828436E0,
 194       0.9950951941E-1, 0.8733827185E-2
 195     };
 196
 197   static double P4[4] =
 198     {
 199       0.1000351630E1, 0.4503592498E1,
 200       0.1085883880E2, 0.7536052269E1
 201     };
 202
 203   static double P5[4] =
 204     {
 205       0.1000006517E1, 0.4909414111E2,
 206       0.8505544753E2, 0.1532153455E3
 207     };
 208
 209   static double P6[4] =
 210     {
 211       0.1000000983E1, 0.1329868456E3,
 212       0.9162149244E3, -0.9605054274E3
 213     };
 214
 215   static double Q1[5] =
 216     {
 217       1.0, -0.5571175625E-2,
 218       0.6225310236E-1, -0.3137378427E-2,
 219       0.1931496439E-2
 220     };
 221
 222   static double Q2[4] =
 223     {
 224       1.0, 0.6191136137E0,
 225       0.1720721448E0, 0.2278594771E-1
 226     };
 227
 228   static double Q3[4] =
 229     {
 230       1.0, 0.4237190502E0,
 231       0.1095631512E0, 0.8693851567E-2
 232     };
 233
 234   static double Q4[4] =
 235     {
 236       1.0, 0.5539969678E1,
 237       0.1933581111E2, 0.2721321508E2
 238     };
 239
 240   static double Q5[4] =
 241     {
 242       1.0, 0.5009928881E2,
 243       0.1399819104E3, 0.4200002909E3
 244     };
 245
 246   static double Q6[4] =
 247     {
 248       1.0, 0.1339887843E3,
 249       0.1055990413E4, 0.5532224619E3
 250     };
 251
 252   static double A1[3] =
 253     {
 254       -0.4583333333E0, 0.6675347222E0, -0.1641741416E1
 255     };
 256
 257   static double A2[3] =
 258     {
 259       1.0, -0.4227843351E0, -0.2043403138E1
 260     };
 261
 262   double U, V, DISLAN;
 263
 264   V = x;
 265   if (V < -5.5)
 266     {
 267       U = exp(V + 1);
 268       DISLAN = 0.3989422803 * exp( -1 / U) * sqrt(U) *
 269                (1 + (A1[0] + (A1[1] + A1[2] * U) * U) * U);
 270     }
 271   else if (V < -1)
 272     {
 273       U = exp( -V - 1);
 274       DISLAN = (exp( -U) / sqrt(U)) *
 275                (P1[0] + (P1[1] + (P1[2] + (P1[3] + P1[4] * V) * V) * V) * V) /
 276                (Q1[0] + (Q1[1] + (Q1[2] + (Q1[3] + Q1[4] * V) * V) * V) * V);
 277     }
 278   else if (V < 1)
 279     {
 280       DISLAN = (P2[0] + (P2[1] + (P2[2] + P2[3] * V) * V) * V) /
 281                (Q2[0] + (Q2[1] + (Q2[2] + Q2[3] * V) * V) * V);
 282     }
 283   else if (V < 4)
 284     {
 285       DISLAN = (P3[0] + (P3[1] + (P3[2] + P3[3] * V) * V) * V) /
 286                (Q3[0] + (Q3[1] + (Q3[2] + Q3[3] * V) * V) * V);
 287     }
 288   else if (V < 12)
 289     {
 290       U = 1 / V;
 291       DISLAN = (P4[0] + (P4[1] + (P4[2] + P4[3] * U) * U) * U) /
 292                (Q4[0] + (Q4[1] + (Q4[2] + Q4[3] * U) * U) * U);
 293     }
 294   else if (V < 50)
 295     {
 296       U = 1 / V;
 297       DISLAN = (P5[0] + (P5[1] + (P5[2] + P5[3] * U) * U) * U) /
 298                (Q5[0] + (Q5[1] + (Q5[2] + Q5[3] * U) * U) * U);
 299     }
 300   else if (V < 300)
 301     {
 302       U = 1 / V;
 303       DISLAN = (P6[0] + (P6[1] + (P6[2] + P6[3] * U) * U) * U) /
 304                (Q6[0] + (Q6[1] + (Q6[2] + Q6[3] * U) * U) * U);
 305     }
 306   else
 307     {
 308       U = 1 / (V - V * log(V) / (V + 1));
 309       DISLAN = 1 - (A2[0] + (A2[1] + A2[2] * U) * U) * U;
 310     }
 311
 312   return DISLAN;
 313 }
 314 #endif
 315
 316 double
 317 gsl_ran_landau(const gsl_rng * r)
 318 {
 319   static double F[983] =
 320     {
 321       0.0000000,   /* Add empty element [0] to account for difference
 322                       between C and Fortran convention for lower bound. */
 323       00.000000, 00.000000, 00.000000, 00.000000, 00.000000,
 324       -2.244733, -2.204365, -2.168163, -2.135219, -2.104898,
 325       -2.076740, -2.050397, -2.025605, -2.002150, -1.979866,
 326       -1.958612, -1.938275, -1.918760, -1.899984, -1.881879,
 327       -1.864385, -1.847451, -1.831030, -1.815083, -1.799574,
 328       -1.784473, -1.769751, -1.755383, -1.741346, -1.727620,
 329       -1.714187, -1.701029, -1.688130, -1.675477, -1.663057,
 330       -1.650858, -1.638868, -1.627078, -1.615477, -1.604058,
 331       -1.592811, -1.581729, -1.570806, -1.560034, -1.549407,
 332       -1.538919, -1.528565, -1.518339, -1.508237, -1.498254,
 333       -1.488386, -1.478628, -1.468976, -1.459428, -1.449979,
 334       -1.440626, -1.431365, -1.422195, -1.413111, -1.404112,
 335       -1.395194, -1.386356, -1.377594, -1.368906, -1.360291,
 336       -1.351746, -1.343269, -1.334859, -1.326512, -1.318229,
 337       -1.310006, -1.301843, -1.293737, -1.285688, -1.277693,
 338       -1.269752, -1.261863, -1.254024, -1.246235, -1.238494,
 339       -1.230800, -1.223153, -1.215550, -1.207990, -1.200474,
 340       -1.192999, -1.185566, -1.178172, -1.170817, -1.163500,
 341       -1.156220, -1.148977, -1.141770, -1.134598, -1.127459,
 342       -1.120354, -1.113282, -1.106242, -1.099233, -1.092255,
 343       -1.085306, -1.078388, -1.071498, -1.064636, -1.057802,
 344       -1.050996, -1.044215, -1.037461, -1.030733, -1.024029,
 345       -1.017350, -1.010695, -1.004064, -0.997456, -0.990871,
 346       -0.984308, -0.977767, -0.971247, -0.964749, -0.958271,
 347       -0.951813, -0.945375, -0.938957, -0.932558, -0.926178,
 348       -0.919816, -0.913472, -0.907146, -0.900838, -0.894547,
 349       -0.888272, -0.882014, -0.875773, -0.869547, -0.863337,
 350       -0.857142, -0.850963, -0.844798, -0.838648, -0.832512,
 351       -0.826390, -0.820282, -0.814187, -0.808106, -0.802038,
 352       -0.795982, -0.789940, -0.783909, -0.777891, -0.771884,
 353       -0.765889, -0.759906, -0.753934, -0.747973, -0.742023,
 354       -0.736084, -0.730155, -0.724237, -0.718328, -0.712429,
 355       -0.706541, -0.700661, -0.694791, -0.688931, -0.683079,
 356       -0.677236, -0.671402, -0.665576, -0.659759, -0.653950,
 357       -0.648149, -0.642356, -0.636570, -0.630793, -0.625022,
 358       -0.619259, -0.613503, -0.607754, -0.602012, -0.596276,
 359       -0.590548, -0.584825, -0.579109, -0.573399, -0.567695,
 360       -0.561997, -0.556305, -0.550618, -0.544937, -0.539262,
 361       -0.533592, -0.527926, -0.522266, -0.516611, -0.510961,
 362       -0.505315, -0.499674, -0.494037, -0.488405, -0.482777,
 363       -0.477153, -0.471533, -0.465917, -0.460305, -0.454697,
 364       -0.449092, -0.443491, -0.437893, -0.432299, -0.426707,
 365       -0.421119, -0.415534, -0.409951, -0.404372, -0.398795,
 366       -0.393221, -0.387649, -0.382080, -0.376513, -0.370949,
 367       -0.365387, -0.359826, -0.354268, -0.348712, -0.343157,
 368       -0.337604, -0.332053, -0.326503, -0.320955, -0.315408,
 369       -0.309863, -0.304318, -0.298775, -0.293233, -0.287692,
 370       -0.282152, -0.276613, -0.271074, -0.265536, -0.259999,
 371       -0.254462, -0.248926, -0.243389, -0.237854, -0.232318,
 372       -0.226783, -0.221247, -0.215712, -0.210176, -0.204641,
 373       -0.199105, -0.193568, -0.188032, -0.182495, -0.176957,
 374       -0.171419, -0.165880, -0.160341, -0.154800, -0.149259,
 375       -0.143717, -0.138173, -0.132629, -0.127083, -0.121537,
 376       -0.115989, -0.110439, -0.104889, -0.099336, -0.093782,
 377       -0.088227, -0.082670, -0.077111, -0.071550, -0.065987,
 378       -0.060423, -0.054856, -0.049288, -0.043717, -0.038144,
 379       -0.032569, -0.026991, -0.021411, -0.015828, -0.010243,
 380       -0.004656, 00.000934, 00.006527, 00.012123, 00.017722,
 381       00.023323, 00.028928, 00.034535, 00.040146, 00.045759,
 382       00.051376, 00.056997, 00.062620, 00.068247, 00.073877,
 383       00.079511, 00.085149, 00.090790, 00.096435, 00.102083,
 384       00.107736, 00.113392, 00.119052, 00.124716, 00.130385,
 385       00.136057, 00.141734, 00.147414, 00.153100, 00.158789,
 386       00.164483, 00.170181, 00.175884, 00.181592, 00.187304,
 387       00.193021, 00.198743, 00.204469, 00.210201, 00.215937,
 388       00.221678, 00.227425, 00.233177, 00.238933, 00.244696,
 389       00.250463, 00.256236, 00.262014, 00.267798, 00.273587,
 390       00.279382, 00.285183, 00.290989, 00.296801, 00.302619,
 391       00.308443, 00.314273, 00.320109, 00.325951, 00.331799,
 392       00.337654, 00.343515, 00.349382, 00.355255, 00.361135,
 393       00.367022, 00.372915, 00.378815, 00.384721, 00.390634,
 394       00.396554, 00.402481, 00.408415, 00.414356, 00.420304,
 395       00.426260, 00.432222, 00.438192, 00.444169, 00.450153,
 396       00.456145, 00.462144, 00.468151, 00.474166, 00.480188,
 397       00.486218, 00.492256, 00.498302, 00.504356, 00.510418,
 398       00.516488, 00.522566, 00.528653, 00.534747, 00.540850,
 399       00.546962, 00.553082, 00.559210, 00.565347, 00.571493,
 400       00.577648, 00.583811, 00.589983, 00.596164, 00.602355,
 401       00.608554, 00.614762, 00.620980, 00.627207, 00.633444,
 402       00.639689, 00.645945, 00.652210, 00.658484, 00.664768,
 403       00.671062, 00.677366, 00.683680, 00.690004, 00.696338,
 404       00.702682, 00.709036, 00.715400, 00.721775, 00.728160,
 405       00.734556, 00.740963, 00.747379, 00.753807, 00.760246,
 406       00.766695, 00.773155, 00.779627, 00.786109, 00.792603,
 407       00.799107, 00.805624, 00.812151, 00.818690, 00.825241,
 408       00.831803, 00.838377, 00.844962, 00.851560, 00.858170,
 409       00.864791, 00.871425, 00.878071, 00.884729, 00.891399,
 410       00.898082, 00.904778, 00.911486, 00.918206, 00.924940,
 411       00.931686, 00.938446, 00.945218, 00.952003, 00.958802,
 412       00.965614, 00.972439, 00.979278, 00.986130, 00.992996,
 413       00.999875, 01.006769, 01.013676, 01.020597, 01.027533,
 414       01.034482, 01.041446, 01.048424, 01.055417, 01.062424,
 415       01.069446, 01.076482, 01.083534, 01.090600, 01.097681,
 416       01.104778, 01.111889, 01.119016, 01.126159, 01.133316,
 417       01.140490, 01.147679, 01.154884, 01.162105, 01.169342,
 418       01.176595, 01.183864, 01.191149, 01.198451, 01.205770,
 419       01.213105, 01.220457, 01.227826, 01.235211, 01.242614,
 420       01.250034, 01.257471, 01.264926, 01.272398, 01.279888,
 421       01.287395, 01.294921, 01.302464, 01.310026, 01.317605,
 422       01.325203, 01.332819, 01.340454, 01.348108, 01.355780,
 423       01.363472, 01.371182, 01.378912, 01.386660, 01.394429,
 424       01.402216, 01.410024, 01.417851, 01.425698, 01.433565,
 425       01.441453, 01.449360, 01.457288, 01.465237, 01.473206,
 426       01.481196, 01.489208, 01.497240, 01.505293, 01.513368,
 427       01.521465, 01.529583, 01.537723, 01.545885, 01.554068,
 428       01.562275, 01.570503, 01.578754, 01.587028, 01.595325,
 429       01.603644, 01.611987, 01.620353, 01.628743, 01.637156,
 430       01.645593, 01.654053, 01.662538, 01.671047, 01.679581,
 431       01.688139, 01.696721, 01.705329, 01.713961, 01.722619,
 432       01.731303, 01.740011, 01.748746, 01.757506, 01.766293,
 433       01.775106, 01.783945, 01.792810, 01.801703, 01.810623,
 434       01.819569, 01.828543, 01.837545, 01.846574, 01.855631,
 435       01.864717, 01.873830, 01.882972, 01.892143, 01.901343,
 436       01.910572, 01.919830, 01.929117, 01.938434, 01.947781,
 437       01.957158, 01.966566, 01.976004, 01.985473, 01.994972,
 438       02.004503, 02.014065, 02.023659, 02.033285, 02.042943,
 439       02.052633, 02.062355, 02.072110, 02.081899, 02.091720,
 440       02.101575, 02.111464, 02.121386, 02.131343, 02.141334,
 441       02.151360, 02.161421, 02.171517, 02.181648, 02.191815,
 442       02.202018, 02.212257, 02.222533, 02.232845, 02.243195,
 443       02.253582, 02.264006, 02.274468, 02.284968, 02.295507,
 444       02.306084, 02.316701, 02.327356, 02.338051, 02.348786,
 445       02.359562, 02.370377, 02.381234, 02.392131, 02.403070,
 446       02.414051, 02.425073, 02.436138, 02.447246, 02.458397,
 447       02.469591, 02.480828, 02.492110, 02.503436, 02.514807,
 448       02.526222, 02.537684, 02.549190, 02.560743, 02.572343,
 449       02.583989, 02.595682, 02.607423, 02.619212, 02.631050,
 450       02.642936, 02.654871, 02.666855, 02.678890, 02.690975,
 451       02.703110, 02.715297, 02.727535, 02.739825, 02.752168,
 452       02.764563, 02.777012, 02.789514, 02.802070, 02.814681,
 453       02.827347, 02.840069, 02.852846, 02.865680, 02.878570,
 454       02.891518, 02.904524, 02.917588, 02.930712, 02.943894,
 455       02.957136, 02.970439, 02.983802, 02.997227, 03.010714,
 456       03.024263, 03.037875, 03.051551, 03.065290, 03.079095,
 457       03.092965, 03.106900, 03.120902, 03.134971, 03.149107,
 458       03.163312, 03.177585, 03.191928, 03.206340, 03.220824,
 459       03.235378, 03.250005, 03.264704, 03.279477, 03.294323,
 460       03.309244, 03.324240, 03.339312, 03.354461, 03.369687,
 461       03.384992, 03.400375, 03.415838, 03.431381, 03.447005,
 462       03.462711, 03.478500, 03.494372, 03.510328, 03.526370,
 463       03.542497, 03.558711, 03.575012, 03.591402, 03.607881,
 464       03.624450, 03.641111, 03.657863, 03.674708, 03.691646,
 465       03.708680, 03.725809, 03.743034, 03.760357, 03.777779,
 466       03.795300, 03.812921, 03.830645, 03.848470, 03.866400,
 467       03.884434, 03.902574, 03.920821, 03.939176, 03.957640,
 468       03.976215, 03.994901, 04.013699, 04.032612, 04.051639,
 469       04.070783, 04.090045, 04.109425, 04.128925, 04.148547,
 470       04.168292, 04.188160, 04.208154, 04.228275, 04.248524,
 471       04.268903, 04.289413, 04.310056, 04.330832, 04.351745,
 472       04.372794, 04.393982, 04.415310, 04.436781, 04.458395,
 473       04.480154, 04.502060, 04.524114, 04.546319, 04.568676,
 474       04.591187, 04.613854, 04.636678, 04.659662, 04.682807,
 475       04.706116, 04.729590, 04.753231, 04.777041, 04.801024,
 476       04.825179, 04.849511, 04.874020, 04.898710, 04.923582,
 477       04.948639, 04.973883, 04.999316, 05.024942, 05.050761,
 478       05.076778, 05.102993, 05.129411, 05.156034, 05.182864,
 479       05.209903, 05.237156, 05.264625, 05.292312, 05.320220,
 480       05.348354, 05.376714, 05.405306, 05.434131, 05.463193,
 481       05.492496, 05.522042, 05.551836, 05.581880, 05.612178,
 482       05.642734, 05.673552, 05.704634, 05.735986, 05.767610,
 483       05.799512, 05.831694, 05.864161, 05.896918, 05.929968,
 484       05.963316, 05.996967, 06.030925, 06.065194, 06.099780,
 485       06.134687, 06.169921, 06.205486, 06.241387, 06.277630,
 486       06.314220, 06.351163, 06.388465, 06.426130, 06.464166,
 487       06.502578, 06.541371, 06.580553, 06.620130, 06.660109,
 488       06.700495, 06.741297, 06.782520, 06.824173, 06.866262,
 489       06.908795, 06.951780, 06.995225, 07.039137, 07.083525,
 490       07.128398, 07.173764, 07.219632, 07.266011, 07.312910,
 491       07.360339, 07.408308, 07.456827, 07.505905, 07.555554,
 492       07.605785, 07.656608, 07.708035, 07.760077, 07.812747,
 493       07.866057, 07.920019, 07.974647, 08.029953, 08.085952,
 494       08.142657, 08.200083, 08.258245, 08.317158, 08.376837,
 495       08.437300, 08.498562, 08.560641, 08.623554, 08.687319,
 496       08.751955, 08.817481, 08.883916, 08.951282, 09.019600,
 497       09.088889, 09.159174, 09.230477, 09.302822, 09.376233,
 498       09.450735, 09.526355, 09.603118, 09.681054, 09.760191,
 499       09.840558, 09.922186, 10.005107, 10.089353, 10.174959,
 500       10.261958, 10.350389, 10.440287, 10.531693, 10.624646,
 501       10.719188, 10.815362, 10.913214, 11.012789, 11.114137,
 502       11.217307, 11.322352, 11.429325, 11.538283, 11.649285,
 503       11.762390, 11.877664, 11.995170, 12.114979, 12.237161,
 504       12.361791, 12.488946, 12.618708, 12.751161, 12.886394,
 505       13.024498, 13.165570, 13.309711, 13.457026, 13.607625,
 506       13.761625, 13.919145, 14.080314, 14.245263, 14.414134,
 507       14.587072, 14.764233, 14.945778, 15.131877, 15.322712,
 508       15.518470, 15.719353, 15.925570, 16.137345, 16.354912,
 509       16.578520, 16.808433, 17.044929, 17.288305, 17.538873,
 510       17.796967, 18.062943, 18.337176, 18.620068, 18.912049,
 511       19.213574, 19.525133, 19.847249, 20.180480, 20.525429,
 512       20.882738, 21.253102, 21.637266, 22.036036, 22.450278,
 513       22.880933, 23.329017, 23.795634, 24.281981, 24.789364,
 514       25.319207, 25.873062, 26.452634, 27.059789, 27.696581,
 515       28.365274, 29.068370, 29.808638, 30.589157, 31.413354,
 516       32.285060, 33.208568, 34.188705, 35.230920, 36.341388,
 517       37.527131, 38.796172, 40.157721, 41.622399, 43.202525,
 518       44.912465, 46.769077, 48.792279, 51.005773, 53.437996,
 519       56.123356, 59.103894
 520     };
 521   double X, U, V, RANLAN;
 522   int I;
 523
 524   X = gsl_rng_uniform_pos(r);
 525   U = 1000.0 * X;
 526   I = U;
 527   U = U - I;
 528
 529   if (I >= 70 && I <= 800)
 530     {
 531       RANLAN = F[I] + U * (F[I + 1] - F[I]);
 532     }
 533   else if (I >= 7 && I <= 980)
 534     {
 535       RANLAN = F[I]
 536         + U * (F[I + 1] - F[I]
 537                - 0.25 * (1 - U) * (F[I + 2] - F[I + 1] - F[I] + F[I - 1]));
 538     }
 539   else if (I < 7)
 540     {
 541       V = log(X);
 542       U = 1 / V;
 543       RANLAN = ((0.99858950 + (3.45213058E1 + 1.70854528E1 * U) * U) /
 544                 (1 + (3.41760202E1 + 4.01244582 * U) * U)) *
 545                ( -log( -0.91893853 - V) - 1);
 546     }
 547   else
 548     {
 549       U = 1 - X;
 550       V = U * U;
 551       if (X <= 0.999)
 552         {
 553           RANLAN = (1.00060006 + 2.63991156E2 * U + 4.37320068E3 * V) /
 554                    ((1 + 2.57368075E2 * U + 3.41448018E3 * V) * U);
 555         }
 556       else
 557         {
 558           RANLAN = (1.00001538 + 6.07514119E3 * U + 7.34266409E5 * V) /
 559                    ((1 + 6.06511919E3 * U + 6.94021044E5 * V) * U);
 560         }
 561     }
 562
 563   return RANLAN;
 564 }
 565