/* 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;
}