Print

Print
Commit in lcsim/src/org/lcsim/math on MAIN
chisq/ChisqProb.java+247added 1.1
distribution/LandauDistribution.java+137added 1.1
            /MoyalDistribution.java+23added 1.1
            /RandLandau.java+582added 1.1
+989
4 added files
First release of some basic math distributions. May move to freehep.
lcsim/src/org/lcsim/math/chisq
ChisqProb.java added at 1.1
diff -N ChisqProb.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ ChisqProb.java	18 Jul 2005 07:12:46 -0000	1.1
@@ -0,0 +1,247 @@

+package org.lcsim.math.chisq;
+/** A utility class to return the probability that a value of chi-squared <b> x </b> measured
+ * for a system with <b> n </b>degrees of freedom would be exceeded by chance. Based on routines
+ * in the book <em> "Numerical Recipes: The Art of Scientific Computing"</em>.
+ * @author Norman Graf
+ * @version 1.0
+ *
+ */
+
+public class ChisqProb
+{
+    
+    /**Returns the incomplete gamma function P(a, x).
+     *
+     * <br clear="all" /><table border="0" width="100%"><tr><td>
+     * <table align="center"><tr><td nowrap="nowrap" align="center">
+     * <i>P</i>(<i>a</i>,<i>x</i>)  <font face="symbol">º</font
+     * > </td><td nowrap="nowrap" align="center">
+     * <font face="symbol">g</font
+     * >(<i>a</i>,<i>x</i>)<hr noshade="noshade" /><font face="symbol">G</font
+     * >(<i>a</i>)<br /></td><td nowrap="nowrap" align="center">
+     * <font face="symbol">º</font
+     * > </td><td nowrap="nowrap" align="center">
+     * 1<hr noshade="noshade" /><font face="symbol">G</font
+     * >(<i>a</i>)<br /></td><td nowrap="nowrap" align="center">
+     * </td><td nowrap="nowrap" align="center">
+     * <font size="-1"><i>x</i></font><!--sup
+     * --><br /><font face="symbol">ó<br />õ<br /></font><font size="-1">0</font>&nbsp;<br /></td><td nowrap="nowrap" align="center">
+     * <i>e</i><sup> <font face="symbol">-</font
+     * > <i>t</i></sup> <i>t</i><sup><i>a</i> <font face="symbol">-</font
+     * > 1</sup> <i>dt</i> </td></tr></table>
+     * </td></tr></table>
+     *
+     * This is the probability that the observed chi-square for a correct
+     * model should be less than a value of chi-square.
+     *  @param int a   : the number of degrees of freedom
+     *  @param double x: the value of chi-squared
+     *  @return The incomplete gamma function P(a, x).
+     */
+    public static double gammp(int a, double x)
+    {
+        if (x < 0.0 || a <= 0.0) System.out.println("Invalid arguments in routine gammp");
+        if (x < (a+1.0))
+        { //Use the series representation.
+            return gser(a/2.,x/2.);
+        }
+        else
+        { //Use the continued fraction representation and take its complement.
+            return 1.0-gcf(a/2.,x/2.);
+        }
+    }
+    
+    /**Returns the incomplete gamma function Q(a, x)= 1 - P(a, x).
+     *
+     * <br clear="all" /><table border="0" width="100%"><tr><td>
+     * <table align="center"><tr><td nowrap="nowrap" align="center">
+     * <i>Q</i>(<i>a</i>,<i>x</i>)  <font face="symbol">º</font
+     * > 1 <font face="symbol">-</font
+     * > <i>P</i>(<i>a</i>,<i>x</i>)  <font face="symbol">º</font
+     * > </td><td nowrap="nowrap" align="center">
+     * <font face="symbol">G</font
+     * >(<i>a</i>,<i>x</i>)<hr noshade="noshade" /><font face="symbol">G</font
+     * >(<i>a</i>)<br /></td><td nowrap="nowrap" align="center">
+     * </td><td nowrap="nowrap" align="center">
+     * <font size="-1"><font face="symbol">¥</font
+     * ></font><!--sup
+     * --><br /><font face="symbol">ó<br />õ<br /></font><font size="-1"><i>x</i></font>&nbsp;<br /></td><td nowrap="nowrap" align="center">
+     * <i>e</i><sup> <font face="symbol">-</font
+     * > <i>t</i></sup> <i>t</i><sup><i>a</i> <font face="symbol">-</font
+     * > 1</sup> <i>dt</i> </td></tr></table>
+     * </td></tr></table>
+     *
+     * This is the probability that the observed chi-square
+     * will exceed the value of chi-square by chance even for a correct model.
+     *  @param int a   : the number of degrees of freedom
+     *  @param double x: the value of chi-squared
+     *  @return The incomplete gamma function Q(a, x)= 1 - P(a, x).
+     */
+    public static double gammq(int a, double x)
+    {
+        
+        if (x < 0.0 || a <= 0.0) System.out.println("Invalid arguments in routine gammq");
+        if (x < (a+1.0))
+        { //Use the series representation and take its complement.
+            return 1.0-gser(a/2.,x/2.);
+        }
+        else
+        { //Use the continued fraction representation.
+            return gcf(a/2.,x/2.);
+        }
+    }
+    
+    /**Returns the incomplete gamma function P(a, x) evaluated by its series representation.
+     *
+     * <br clear="all" /><table border="0" width="100%"><tr><td>
+     * <table align="center"><tr><td nowrap="nowrap" align="center">
+     * <font face="symbol">g</font
+     * >(<i>a</i>,<i>x</i>) = <i>e</i><sup> <font face="symbol">-</font
+     * > <i>x</i></sup> <i>x</i><sup><i>a</i></sup> </td><td nowrap="nowrap" align="center">
+     * <font size="-1"><font face="symbol">¥</font
+     * ></font><!--sup
+     * --><br /><font size="+3"><font face="symbol">å<br />
+     * </font></font><font size="-1"><i>n</i> = 0</font>&nbsp;<br /></td><td nowrap="nowrap" align="center">
+     * </td><td nowrap="nowrap" align="center">
+     * &nbsp;<font face="symbol">G</font
+     * >(<i>a</i>)
+     * <div class="hrcomp"><hr noshade="noshade" size="1"/></div><font face="symbol">G</font
+     * >(<i>a</i> + 1 + <i>n</i>)<br /></td><td nowrap="nowrap" align="center">
+     * <i>x</i><sup><i>n</i></sup></td></tr></table>
+     * </td></tr></table>
+     *
+     *  @param double a   : the number of degrees of freedom/2
+     *  @param double x: the value of chi-squared/2
+     *  @return The incomplete gamma P(a, x) evaluated by its series representation.
+     */
+    public static double gser(double a, double x)
+    {
+        int ITMAX=100;
+        double EPS=3.0e-7;
+        
+        int n;
+        double sum,del,ap;
+        double gln=gammln(a);
+        if (x <= 0.0)
+        {
+            if (x < 0.0) System.out.println("x less than 0 in routine gser");
+            return 0.0;
+        }
+        else
+        {
+            ap=a;
+            del=sum=1.0/a;
+            for (n=1;n<=ITMAX;n++)
+            {
+                ++ap;
+                del *= x/ap;
+                sum += del;
+                if (Math.abs(del) < Math.abs(sum)*EPS)
+                {
+                    return sum*Math.exp(-x+a*Math.log(x)-(gln));
+                }
+            }
+            System.out.println("a too large, ITMAX too small in routine gser");
+            return 0.;
+        }
+    }
+    
+    /**Returns the incomplete gamma function Q(a, x) evaluated by its continued fraction representation.
+     *
+     * <br clear="all" /><table border="0" width="100%"><tr><td>
+     * <table align="center"><tr><td nowrap="nowrap" align="center">
+     * <font face="symbol">G</font
+     * >(<i>a</i>,<i>x</i>) = <i>e</i><sup> <font face="symbol">-</font
+     * > <i>x</i></sup> <i>x</i><sup><i>a</i></sup> </td><td align="left" class="cl"><font face="symbol">
+     * æ<br />è
+     * </font> </td><td nowrap="nowrap" align="center">
+     * &nbsp;1
+     * <div class="hrcomp"><hr noshade="noshade" size="1"/></div><i>x</i> + <br /></td><td nowrap="nowrap" align="center">
+     * </td><td nowrap="nowrap" align="center">
+     * &nbsp;1 <font face="symbol">-</font
+     * > <i>a</i>
+     * <div class="hrcomp"><hr noshade="noshade" size="1"/></div>1 + <br /></td><td nowrap="nowrap" align="center">
+     * </td><td nowrap="nowrap" align="center">
+     * &nbsp;1
+     * <div class="hrcomp"><hr noshade="noshade" size="1"/></div><i>x</i> + <br /></td><td nowrap="nowrap" align="center">
+     * </td><td nowrap="nowrap" align="center">
+     * &nbsp;2 <font face="symbol">-</font
+     * > <i>a</i>
+     * <div class="hrcomp"><hr noshade="noshade" size="1"/></div>1 + <br /></td><td nowrap="nowrap" align="center">
+     * </td><td nowrap="nowrap" align="center">
+     * &nbsp;2
+     * <div class="hrcomp"><hr noshade="noshade" size="1"/></div><i>x</i> + <br /></td><td nowrap="nowrap" align="center">
+     * <font face="symbol">¼</font
+     * > </td><td align="left" class="cl"><font face="symbol">
+     * ö<br />ø
+     * </font></td><td nowrap="nowrap" align="center">
+     * &nbsp;&nbsp;&nbsp;(<i>x</i>  &gt;  0)</td></tr></table>
+     * </td></tr></table>
+     *
+     *  @param double a   : the number of degrees of freedom/2
+     *  @param double x: the value of chi-squared/2
+     *  @return The incomplete gamma Q(a, x) evaluated by its continued fraction representation.
+     */
+    public static double gcf(double a, double x)
+    {
+        int ITMAX=100; //Maximum allowed number of iteration
+        double EPS=3.0e-7; //Relative accuracy.
+        double FPMIN=1.0e-30; //Number near the smallest representable doubleing-point number.
+        {
+            int i;
+            double an,b,c,d,del,h;
+            double gln=gammln(a);
+            b=x+1.0-a; //Set up for evaluating continued fraction by modified Lentz's method with b0 = 0.
+            c=1.0/FPMIN;
+            d=1.0/b;
+            h=d;
+            for (i=1;i<=ITMAX;i++)
+            { //Iterate to convergence.
+                an = -i*(i-a);
+                b += 2.0;
+                d=an*d+b;
+                if (Math.abs(d) < FPMIN) d=FPMIN;
+                c=b+an/c;
+                if (Math.abs(c) < FPMIN) c=FPMIN;
+                d=1.0/d;
+                del=d*c;
+                h *= del;
+                if (Math.abs(del-1.0) < EPS) break;
+            }
+            if (i > ITMAX) System.out.println("a too large, ITMAX too small in gcf");
+            return Math.exp(-x+a*Math.log(x)-(gln))*h; //Put factors in front.
+        }
+    }
+    
+    /**Returns the logarithm of the Gamma function of x, for x>0.
+     *
+     * <br clear="all" /><table border="0" width="100%"><tr><td>
+     * <table align="center"><tr><td nowrap="nowrap" align="center">
+     * <font face="symbol">G</font
+     * >(<i>z</i>) = </td><td nowrap="nowrap" align="center">
+     * <font size="-1"><font face="symbol">¥</font
+     * ></font><!--sup
+     * --><br /><font face="symbol">ó<br />õ<br /></font><font size="-1">0</font>&nbsp;<br /></td><td nowrap="nowrap" align="center">
+     * <i>t</i><sup><i>z</i> <font face="symbol">-</font
+     * > 1</sup> <i>e</i><sup> <font face="symbol">-</font
+     * > <i>t</i></sup> <i>dt</i> </td></tr></table>
+     * </td></tr></table>
+     *
+     *@param double x
+     *@return The value ln(Gamma(x))
+     */
+    public static double gammln(double x)
+    {
+        double y,tmp,ser;
+        double[] cof=
+        {76.18009172947146,-86.50532032941677,
+         24.01409824083091,-1.231739572450155,
+         0.1208650973866179e-2,-0.5395239384953e-5};
+         int j;
+         y=x;
+         tmp=x+5.5;
+         tmp -= (x+0.5)*Math.log(tmp);
+         ser=1.000000000190015;
+         for (j=0;j<=5;j++) ser += cof[j]/++y;
+         return -tmp+Math.log(2.5066282746310005*ser/x);
+    }
+}

\ No newline at end of file
lcsim/src/org/lcsim/math/distribution
LandauDistribution.java added at 1.1
diff -N LandauDistribution.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ LandauDistribution.java	18 Jul 2005 07:12:46 -0000	1.1
@@ -0,0 +1,137 @@

+/*
+ * LandauDistribution.java
+ * This corresponds to the CLHEP version, which
+ * is based on the CERNLIB routine denlan (G110).
+ * Note that there is some variation in whether
+ * to divide by sigma, and whether the peak
+ * coincides with the most probable value.
+ *
+ * Created on May 23, 2005, 4:21 PM
+ */
+
+package org.lcsim.math.distribution;
+
+import static java.lang.Math.sqrt;
+import static java.lang.Math.exp;
+import static java.lang.Math.log;
+/**
+ *
+ * @author ngraf
+ */
+public class LandauDistribution
+{
+    private double _peak;
+    private double _width;
+    
+    /** Creates a new instance of LandauDistribution */
+    public LandauDistribution(double peak, double width)
+    {
+        _peak = peak;
+        _width = width;
+    }
+    
+    public double peak()
+    {
+        return _peak;
+    }
+    
+    public double width()
+    {
+        return _width;
+    }
+    
+    public double value(double x)
+    {
+        double xs  = _peak + 0.222782*_width;
+        return denlan((x-xs)/_width)/_width;
+    }
+    
+    private double denlan(double x)
+    {
+        /* Initialized data */
+        
+        /* System generated locals */
+        double ret_val, r__1;
+        
+        /* Local variables */
+        double u, v;
+        v = x;
+        if (v < -5.5)
+        {
+            u = exp(v + 1.);
+            ret_val = exp(-1 / u) / sqrt(u) * .3989422803 * ((a1[0] + (a1[
+                    1] + a1[2] * u) * u) * u + 1);
+        }
+        else if (v < -1.)
+        {
+            u = exp(-v - 1);
+            ret_val = 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.)
+        {
+            ret_val = (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.)
+        {
+            ret_val = (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;
+            /* Computing 2nd power */
+            r__1 = u;
+            ret_val = r__1 * r__1 * (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;
+            /* Computing 2nd power */
+            r__1 = u;
+            ret_val = r__1 * r__1 * (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;
+            /* Computing 2nd power */
+            r__1 = u;
+            ret_val = r__1 * r__1 * (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));
+            /* Computing 2nd power */
+            r__1 = u;
+            ret_val = r__1 * r__1 * ((a2[0] + a2[1] * u) * u + 1);
+        }
+        return ret_val;
+        
+    }
+    
+    static final double[] p1 = { .4259894875,-.124976255,.039842437,-.006298287635,.001511162253 };
+    static final double[] q5 = { 1.,156.9424537,3745.310488,9834.698876,66924.28357 };
+    static final double[] p6 = { 1.000827619,664.9143136,62972.92665,475554.6998,-5743609.109 };
+    static final double[] q6 = { 1.,651.4101098,56974.73333,165917.4725,-2815759.939 };
+    static final double[] a1 = { .04166666667,-.01996527778,.02709538966 };
+    static final double[] a2 = { -1.84556867,-4.284640743 };
+    static final double[] q1 = { 1.,-.3388260629,.09594393323,-.01608042283,.003778942063 };
+    static final double[] p2 = { .1788541609,.1173957403,.01488850518,-.001394989411,1.283617211e-4 };
+    static final double[] q2 = { 1.,.7428795082,.3153932961,.06694219548,.008790609714 };
+    static final double[] p3 = { .1788544503,.09359161662,.006325387654,6.611667319e-5,-2.031049101e-6 };
+    static final double[] q3 = { 1.,.6097809921,.2560616665,.04746722384,.006957301675 };
+    static final double[] p4 = { .9874054407,118.6723273,849.279436,-743.7792444,427.0262186 };
+    static final double[] q4 = { 1.,106.8615961,337.6496214,2016.712389,1597.063511 };
+    static final double[] p5 = { 1.003675074,167.5702434,4789.711289,21217.86767,-22324.9491 };
+    
+}
lcsim/src/org/lcsim/math/distribution
MoyalDistribution.java added at 1.1
diff -N MoyalDistribution.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ MoyalDistribution.java	18 Jul 2005 07:12:46 -0000	1.1
@@ -0,0 +1,23 @@

+import java.util.Random;
+
+public class MoyalDistribution
+{
+    
+    private String[] parameters = {"N", "x0", "sigma"};
+    private double[] parValues = {1, 1, 1};
+    private String title = "MoyalDistribution";
+    private String[] variables = {"x"};
+    
+    public double value(double x)
+    {
+        double var = x;
+        double n = parValues[0];
+        double x0 = parValues[1];
+        double s  = parValues[2];
+        double lambda = (var-x0)/s;
+        double arg = lambda + Math.exp( -1.*lambda );
+        double num = Math.exp( -1.*arg );
+        return n*Math.sqrt( num/(2*Math.PI) );
+    }
+    
+}
lcsim/src/org/lcsim/math/distribution
RandLandau.java added at 1.1
diff -N RandLandau.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ RandLandau.java	18 Jul 2005 07:12:46 -0000	1.1
@@ -0,0 +1,582 @@

+/*
+ * RandLandau.java
+ *
+ * Created on June 6, 2005, 4:05 PM
+ *
+ * To change this template, choose Tools | Options and locate the template under
+ * the Source Creation and Management node. Right-click the template and choose
+ * Open. You can then make changes to the template in the Source Editor.
+ */
+
+package org.lcsim.math.distribution;
+
+
+import static java.lang.Math.log;
+import java.util.Random;
+/**
+ *
+ * @author ngraf
+ */
+public class RandLandau
+{
+    private Random _ran = new Random();
+    private double _sigma;
+    private double _mean;
+    
+    /** Creates a new instance of RandLandau */
+    public RandLandau()
+    {
+        _mean = 0.;
+        _sigma = 1.0;
+    }
+    
+    public RandLandau(double mean, double sigma)
+    {
+        _mean = mean;
+        _sigma = sigma;
+    }
+    
+    
+    public double shoot()
+    {
+        return _mean+_sigma*transform(_ran.nextDouble());
+    }
+    
+    private double transform(double r)
+    {
+        
+        double  u = r * TABLE_MULTIPLIER;
+        int index = (int) u;
+        double du = u - index;
+        
+        // du is scaled such that the we dont have to multiply by TABLE_INTERVAL
+        // when interpolating.
+        
+        // Five cases:
+        // A) Between .070 and .800 the function is so smooth, straight
+        //	linear interpolation is adequate.
+        // B) Between .007 and .070, and between .800 and .980, quadratic
+        //    interpolation is used.  This requires the same 4 points as
+        //	a cubic spline (thus we need .006 and .981 and .982) but
+        //	the quadratic interpolation is accurate enough and quicker.
+        // C) Below .007 an asymptotic expansion for low negative lambda
+        //    (involving two logs) is used; there is a pade-style correction
+        //	factor.
+        // D) Above .980, a simple pade approximation is made (asymptotic to
+        //    1/(1-r)), but...
+        // E) the coefficients in that pade are different above r=.999.
+        
+        if ( index >= 70 && index <= 800 )
+        {		// (A)
+            
+            double f0 = inverseLandau [index];
+            double f1 = inverseLandau [index+1];
+            return f0 + du * (f1 - f0);
+            
+        }
+        else if ( index >= 7 && index <= 980 )
+        {	// (B)
+            
+            double f_1 = inverseLandau [index-1];
+            double f0  = inverseLandau [index];
+            double f1  = inverseLandau [index+1];
+            double f2  = inverseLandau [index+2];
+            
+            return f0 + du * (f1 - f0 - .25*(1-du)* (f2 -f1 - f0 + f_1) );
+            
+        }
+        else if ( index < 7 )
+        {			// (C)
+            
+            final double n0 =  0.99858950;
+            final double n1 = 34.5213058;	final double d1 = 34.1760202;
+            final double n2 = 17.0854528;	final double d2 =  4.01244582;
+            
+            double logr = log(r);
+            double x    = 1/logr;
+            double x2   = x*x;
+            
+            double pade = (n0 + n1*x + n2*x2) / (1.0 + d1*x + d2*x2);
+            
+            return ( - log( -.91893853 - logr ) -1 ) * pade;
+            
+        }
+        else if ( index <= 999 )
+        {			// (D)
+            
+            final double n0 =    1.00060006;
+            final double n1 =  263.991156;	final double d1 =  257.368075;
+            final double n2 = 4373.20068;	final double d2 = 3414.48018;
+            
+            double x = 1-r;
+            double x2   = x*x;
+            
+            return (n0 + n1*x + n2*x2) / (x * (1.0 + d1*x + d2*x2));
+            
+        }
+        else
+        { 					// (E)
+            
+            final double n0 =      1.00001538;
+            final double n1 =   6075.14119;	final double d1 =   6065.11919;
+            final double n2 = 734266.409;	final double d2 = 694021.044;
+            
+            double x = 1-r;
+            double x2   = x*x;
+            
+            return (n0 + n1*x + n2*x2) / (x * (1.0 + d1*x + d2*x2));
+            
+        }
+        
+    } // transform()
+    
+    
+    
+//
+// Table of values of inverse Landau, from r = .060 to .982
+//
+    
+// Since all these are this is static to this compilation unit only, the
+// info is establised a priori and not at each invocation.
+    
+    static final float TABLE_INTERVAL   = .001f;
+    static final int   TABLE_END        =  982;
+    static final float TABLE_MULTIPLIER = 1.0f/TABLE_INTERVAL;
+    
+// Here comes the big (4K bytes) table ---
+//
+// inverseLandau[ n ] = the inverse cdf at r = n*TABLE_INTERVAL = n/1000.
+//
+// Credit CERNLIB for these computations.
+//
+// This data is float because the algortihm does not benefit from further
+// data accuracy.  The numbers below .006 or above .982 are moot since
+// non-table-based methods are used below r=.007 and above .980.
+    
+    static final float[] inverseLandau = {
+        
+        0.0f,							     // .000
+                0.0f, 	    0.0f, 	0.0f, 	    0.0f, 	0.0f, 	     // .001 - .005
+                -2.244733f, -2.204365f, -2.168163f, -2.135219f, -2.104898f,  // .006 - .010
+                -2.076740f, -2.050397f, -2.025605f, -2.002150f, -1.979866f,
+                -1.958612f, -1.938275f, -1.918760f, -1.899984f, -1.881879f,  // .020
+                -1.864385f, -1.847451f, -1.831030f, -1.815083f, -1.799574f,
+                -1.784473f, -1.769751f, -1.755383f, -1.741346f, -1.727620f,  // .030
+                -1.714187f, -1.701029f, -1.688130f, -1.675477f, -1.663057f,
+                -1.650858f, -1.638868f, -1.627078f, -1.615477f, -1.604058f,  // .040
+                -1.592811f, -1.581729f, -1.570806f, -1.560034f, -1.549407f,
+                -1.538919f, -1.528565f, -1.518339f, -1.508237f, -1.498254f,  // .050
+                -1.488386f, -1.478628f, -1.468976f, -1.459428f, -1.449979f,
+                -1.440626f, -1.431365f, -1.422195f, -1.413111f, -1.404112f,  // .060
+                -1.395194f, -1.386356f, -1.377594f, -1.368906f, -1.360291f,
+                -1.351746f, -1.343269f, -1.334859f, -1.326512f, -1.318229f,  // .070
+                -1.310006f, -1.301843f, -1.293737f, -1.285688f, -1.277693f,
+                -1.269752f, -1.261863f, -1.254024f, -1.246235f, -1.238494f,  // .080
+                -1.230800f, -1.223153f, -1.215550f, -1.207990f, -1.200474f,
+                -1.192999f, -1.185566f, -1.178172f, -1.170817f, -1.163500f,  // .090
+                -1.156220f, -1.148977f, -1.141770f, -1.134598f, -1.127459f,
+                -1.120354f, -1.113282f, -1.106242f, -1.099233f, -1.092255f,  // .100
+                
+                -1.085306f, -1.078388f, -1.071498f, -1.064636f, -1.057802f,
+                -1.050996f, -1.044215f, -1.037461f, -1.030733f, -1.024029f,
+                -1.017350f, -1.010695f, -1.004064f, -.997456f,  -.990871f,
+                -.984308f, -.977767f, -.971247f, -.964749f, -.958271f,
+                -.951813f, -.945375f, -.938957f, -.932558f, -.926178f,
+                -.919816f, -.913472f, -.907146f, -.900838f, -.894547f,
+                -.888272f, -.882014f, -.875773f, -.869547f, -.863337f,
+                -.857142f, -.850963f, -.844798f, -.838648f, -.832512f,
+                -.826390f, -.820282f, -.814187f, -.808106f, -.802038f,
+                -.795982f, -.789940f, -.783909f, -.777891f, -.771884f, 	// .150
+                -.765889f, -.759906f, -.753934f, -.747973f, -.742023f,
+                -.736084f, -.730155f, -.724237f, -.718328f, -.712429f,
+                -.706541f, -.700661f, -.694791f, -.688931f, -.683079f,
+                -.677236f, -.671402f, -.665576f, -.659759f, -.653950f,
+                -.648149f, -.642356f, -.636570f, -.630793f, -.625022f,
+                -.619259f, -.613503f, -.607754f, -.602012f, -.596276f,
+                -.590548f, -.584825f, -.579109f, -.573399f, -.567695f,
+                -.561997f, -.556305f, -.550618f, -.544937f, -.539262f,
+                -.533592f, -.527926f, -.522266f, -.516611f, -.510961f,
+                -.505315f, -.499674f, -.494037f, -.488405f, -.482777f,	// .200
+                
+                -.477153f, -.471533f, -.465917f, -.460305f, -.454697f,
+                -.449092f, -.443491f, -.437893f, -.432299f, -.426707f,
+                -.421119f, -.415534f, -.409951f, -.404372f, -.398795f,
+                -.393221f, -.387649f, -.382080f, -.376513f, -.370949f,
+                -.365387f, -.359826f, -.354268f, -.348712f, -.343157f,
+                -.337604f, -.332053f, -.326503f, -.320955f, -.315408f,
+                -.309863f, -.304318f, -.298775f, -.293233f, -.287692f,
+                -.282152f, -.276613f, -.271074f, -.265536f, -.259999f,
+                -.254462f, -.248926f, -.243389f, -.237854f, -.232318f,
+                -.226783f, -.221247f, -.215712f, -.210176f, -.204641f, 	// .250
+                -.199105f, -.193568f, -.188032f, -.182495f, -.176957f,
+                -.171419f, -.165880f, -.160341f, -.154800f, -.149259f,
+                -.143717f, -.138173f, -.132629f, -.127083f, -.121537f,
+                -.115989f, -.110439f, -.104889f, -.099336f, -.093782f,
+                -.088227f, -.082670f, -.077111f, -.071550f, -.065987f,
+                -.060423f, -.054856f, -.049288f, -.043717f, -.038144f,
+                -.032569f, -.026991f, -.021411f, -.015828f, -.010243f,
+                -.004656f,  .000934f,  .006527f,  .012123f,  .017722f,
+                .023323f, .028928f,  .034535f,  .040146f,  .045759f,
+                .051376f, .056997f,  .062620f,  .068247f,  .073877f,	// .300
+                
+                .079511f,  .085149f,  .090790f,  .096435f,  .102083f,
+                .107736f,  .113392f,  .119052f,  .124716f,  .130385f,
+                .136057f,  .141734f,  .147414f,  .153100f,  .158789f,
+                .164483f,  .170181f,  .175884f,  .181592f,  .187304f,
+                .193021f,  .198743f,  .204469f,  .210201f,  .215937f,
+                .221678f,  .227425f,  .233177f,  .238933f,  .244696f,
+                .250463f,  .256236f,  .262014f,  .267798f,  .273587f,
+                .279382f,  .285183f,  .290989f,  .296801f,  .302619f,
+                .308443f,  .314273f,  .320109f,  .325951f,  .331799f,
+                .337654f,  .343515f,  .349382f,  .355255f,  .361135f,	// .350
+                .367022f,  .372915f,  .378815f,  .384721f,  .390634f,
+                .396554f,  .402481f,  .408415f,  .414356f,  .420304f,
+                .426260f,  .432222f,  .438192f,  .444169f,  .450153f,
+                .456145f,  .462144f,  .468151f,  .474166f,  .480188f,
+                .486218f,  .492256f,  .498302f,  .504356f,  .510418f,
+                .516488f,  .522566f,  .528653f,  .534747f,  .540850f,
+                .546962f,  .553082f,  .559210f,  .565347f,  .571493f,
+                .577648f,  .583811f,  .589983f,  .596164f,  .602355f,
+                .608554f,  .614762f,  .620980f,  .627207f,  .633444f,
+                .639689f,  .645945f,  .652210f,  .658484f,  .664768f,	// .400
+                
+                .671062f,  .677366f,  .683680f,  .690004f,  .696338f,
+                .702682f,  .709036f,  .715400f,  .721775f,  .728160f,
+                .734556f,  .740963f,  .747379f,  .753807f,  .760246f,
+                .766695f,  .773155f,  .779627f,  .786109f,  .792603f,
+                .799107f,  .805624f,  .812151f,  .818690f,  .825241f,
+                .831803f,  .838377f,  .844962f,  .851560f,  .858170f,
+                .864791f,  .871425f,  .878071f,  .884729f,  .891399f,
+                .898082f,  .904778f,  .911486f,  .918206f,  .924940f,
+                .931686f,  .938446f,  .945218f,  .952003f,  .958802f,
+                .965614f,  .972439f,  .979278f,  .986130f,  .992996f, 	// .450
+                .999875f,  1.006769f, 1.013676f, 1.020597f, 1.027533f,
+                1.034482f, 1.041446f, 1.048424f, 1.055417f, 1.062424f,
+                1.069446f, 1.076482f, 1.083534f, 1.090600f, 1.097681f,
+                1.104778f, 1.111889f, 1.119016f, 1.126159f, 1.133316f,
+                1.140490f, 1.147679f, 1.154884f, 1.162105f, 1.169342f,
+                1.176595f, 1.183864f, 1.191149f, 1.198451f, 1.205770f,
+                1.213105f, 1.220457f, 1.227826f, 1.235211f, 1.242614f,
+                1.250034f, 1.257471f, 1.264926f, 1.272398f, 1.279888f,
+                1.287395f, 1.294921f, 1.302464f, 1.310026f, 1.317605f,
+                1.325203f, 1.332819f, 1.340454f, 1.348108f, 1.355780f,	// .500
+                
+                1.363472f, 1.371182f, 1.378912f, 1.386660f, 1.394429f,
+                1.402216f, 1.410024f, 1.417851f, 1.425698f, 1.433565f,
+                1.441453f, 1.449360f, 1.457288f, 1.465237f, 1.473206f,
+                1.481196f, 1.489208f, 1.497240f, 1.505293f, 1.513368f,
+                1.521465f, 1.529583f, 1.537723f, 1.545885f, 1.554068f,
+                1.562275f, 1.570503f, 1.578754f, 1.587028f, 1.595325f,
+                1.603644f, 1.611987f, 1.620353f, 1.628743f, 1.637156f,
+                1.645593f, 1.654053f, 1.662538f, 1.671047f, 1.679581f,
+                1.688139f, 1.696721f, 1.705329f, 1.713961f, 1.722619f,
+                1.731303f, 1.740011f, 1.748746f, 1.757506f, 1.766293f, 	// .550
+                1.775106f, 1.783945f, 1.792810f, 1.801703f, 1.810623f,
+                1.819569f, 1.828543f, 1.837545f, 1.846574f, 1.855631f,
+                1.864717f, 1.873830f, 1.882972f, 1.892143f, 1.901343f,
+                1.910572f, 1.919830f, 1.929117f, 1.938434f, 1.947781f,
+                1.957158f, 1.966566f, 1.976004f, 1.985473f, 1.994972f,
+                2.004503f, 2.014065f, 2.023659f, 2.033285f, 2.042943f,
+                2.052633f, 2.062355f, 2.072110f, 2.081899f, 2.091720f,
+                2.101575f, 2.111464f, 2.121386f, 2.131343f, 2.141334f,
+                2.151360f, 2.161421f, 2.171517f, 2.181648f, 2.191815f,
+                2.202018f, 2.212257f, 2.222533f, 2.232845f, 2.243195f,	// .600
+                
+                2.253582f, 2.264006f, 2.274468f, 2.284968f, 2.295507f,
+                2.306084f, 2.316701f, 2.327356f, 2.338051f, 2.348786f,
+                2.359562f, 2.370377f, 2.381234f, 2.392131f, 2.403070f,
+                2.414051f, 2.425073f, 2.436138f, 2.447246f, 2.458397f,
+                2.469591f, 2.480828f, 2.492110f, 2.503436f, 2.514807f,
+                2.526222f, 2.537684f, 2.549190f, 2.560743f, 2.572343f,
+                2.583989f, 2.595682f, 2.607423f, 2.619212f, 2.631050f,
+                2.642936f, 2.654871f, 2.666855f, 2.678890f, 2.690975f,
+                2.703110f, 2.715297f, 2.727535f, 2.739825f, 2.752168f,
+                2.764563f, 2.777012f, 2.789514f, 2.802070f, 2.814681f,	// .650
+                2.827347f, 2.840069f, 2.852846f, 2.865680f, 2.878570f,
+                2.891518f, 2.904524f, 2.917588f, 2.930712f, 2.943894f,
+                2.957136f, 2.970439f, 2.983802f, 2.997227f, 3.010714f,
+                3.024263f, 3.037875f, 3.051551f, 3.065290f, 3.079095f,
+                3.092965f, 3.106900f, 3.120902f, 3.134971f, 3.149107f,
+                3.163312f, 3.177585f, 3.191928f, 3.206340f, 3.220824f,
+                3.235378f, 3.250005f, 3.264704f, 3.279477f, 3.294323f,
+                3.309244f, 3.324240f, 3.339312f, 3.354461f, 3.369687f,
+                3.384992f, 3.400375f, 3.415838f, 3.431381f, 3.447005f,
+                3.462711f, 3.478500f, 3.494372f, 3.510328f, 3.526370f,	// .700
+                
+                3.542497f, 3.558711f, 3.575012f, 3.591402f, 3.607881f,
+                3.624450f, 3.641111f, 3.657863f, 3.674708f, 3.691646f,
+                3.708680f, 3.725809f, 3.743034f, 3.760357f, 3.777779f,
+                3.795300f, 3.812921f, 3.830645f, 3.848470f, 3.866400f,
+                3.884434f, 3.902574f, 3.920821f, 3.939176f, 3.957640f,
+                3.976215f, 3.994901f, 4.013699f, 4.032612f, 4.051639f,
+                4.070783f, 4.090045f, 4.109425f, 4.128925f, 4.148547f,
+                4.168292f, 4.188160f, 4.208154f, 4.228275f, 4.248524f,
+                4.268903f, 4.289413f, 4.310056f, 4.330832f, 4.351745f,
+                4.372794f, 4.393982f, 4.415310f, 4.436781f, 4.458395f,
+                4.480154f, 4.502060f, 4.524114f, 4.546319f, 4.568676f,	// .750
+                4.591187f, 4.613854f, 4.636678f, 4.659662f, 4.682807f,
+                4.706116f, 4.729590f, 4.753231f, 4.777041f, 4.801024f,
+                4.825179f, 4.849511f, 4.874020f, 4.898710f, 4.923582f,
+                4.948639f, 4.973883f, 4.999316f, 5.024942f, 5.050761f,
+                5.076778f, 5.102993f, 5.129411f, 5.156034f, 5.182864f,
+                5.209903f, 5.237156f, 5.264625f, 5.292312f, 5.320220f,
+                5.348354f, 5.376714f, 5.405306f, 5.434131f, 5.463193f,
+                5.492496f, 5.522042f, 5.551836f, 5.581880f, 5.612178f,
+                5.642734f, 5.673552f, 5.704634f, 5.735986f, 5.767610f,	// .800
+                
+                5.799512f, 5.831694f, 5.864161f, 5.896918f, 5.929968f,
+                5.963316f, 5.996967f, 6.030925f, 6.065194f, 6.099780f,
+                6.134687f, 6.169921f, 6.205486f, 6.241387f, 6.277630f,
+                6.314220f, 6.351163f, 6.388465f, 6.426130f, 6.464166f,
+                6.502578f, 6.541371f, 6.580553f, 6.620130f, 6.660109f,
+                6.700495f, 6.741297f, 6.782520f, 6.824173f, 6.866262f,
+                6.908795f, 6.951780f, 6.995225f, 7.039137f, 7.083525f,
+                7.128398f, 7.173764f, 7.219632f, 7.266011f, 7.312910f,
+                7.360339f, 7.408308f, 7.456827f, 7.505905f, 7.555554f,
+                7.605785f, 7.656608f, 7.708035f, 7.760077f, 7.812747f, 	// .850
+                7.866057f, 7.920019f, 7.974647f, 8.029953f, 8.085952f,
+                8.142657f, 8.200083f, 8.258245f, 8.317158f, 8.376837f,
+                8.437300f, 8.498562f, 8.560641f, 8.623554f, 8.687319f,
+                8.751955f, 8.817481f, 8.883916f, 8.951282f, 9.019600f,
+                9.088889f, 9.159174f, 9.230477f, 9.302822f, 9.376233f,
+                9.450735f, 9.526355f, 9.603118f, 9.681054f, 9.760191f,
+                9.840558f,  9.922186f, 10.005107f, 10.089353f, 10.174959f,
+                10.261958f, 10.350389f, 10.440287f, 10.531693f, 10.624646f,
+                10.719188f, 10.815362f, 10.913214f, 11.012789f, 11.114137f,
+                11.217307f, 11.322352f, 11.429325f, 11.538283f, 11.649285f,	// .900
+                
+                11.762390f, 11.877664f, 11.995170f, 12.114979f, 12.237161f,
+                12.361791f, 12.488946f, 12.618708f, 12.751161f, 12.886394f,
+                13.024498f, 13.165570f, 13.309711f, 13.457026f, 13.607625f,
+                13.761625f, 13.919145f, 14.080314f, 14.245263f, 14.414134f,
+                14.587072f, 14.764233f, 14.945778f, 15.131877f, 15.322712f,
+                15.518470f, 15.719353f, 15.925570f, 16.137345f, 16.354912f,
+                16.578520f, 16.808433f, 17.044929f, 17.288305f, 17.538873f,
+                17.796967f, 18.062943f, 18.337176f, 18.620068f, 18.912049f,
+                19.213574f, 19.525133f, 19.847249f, 20.180480f, 20.525429f,
+                20.882738f, 21.253102f, 21.637266f, 22.036036f, 22.450278f, 	// .950
+                22.880933f, 23.329017f, 23.795634f, 24.281981f, 24.789364f,
+                25.319207f, 25.873062f, 26.452634f, 27.059789f, 27.696581f,     // .960
+                28.365274f, 29.068370f, 29.808638f, 30.589157f, 31.413354f,
+                32.285060f, 33.208568f, 34.188705f, 35.230920f, 36.341388f,     // .970
+                37.527131f, 38.796172f, 40.157721f, 41.622399f, 43.202525f,
+                44.912465f, 46.769077f, 48.792279f, 51.005773f, 53.437996f,     // .980
+                56.123356f, 59.103894f, 					// .982
+                
+    };  // End of the inverseLandau table
+    
+    
+// root method
+    
+    static final double[] f = {
+        0       , 0       , 0       ,0        ,0        ,-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, -.997456, -.990871, -.984308, -.977767, -.971247,
+                -.964749, -.958271, -.951813, -.945375, -.938957, -.932558,
+                -.926178, -.919816, -.913472, -.907146, -.900838, -.894547,
+                -.888272, -.882014, -.875773, -.869547, -.863337, -.857142,
+                -.850963, -.844798, -.838648, -.832512, -.826390, -.820282,
+                -.814187, -.808106, -.802038, -.795982, -.789940, -.783909,
+                -.777891, -.771884, -.765889, -.759906, -.753934, -.747973,
+                -.742023, -.736084, -.730155, -.724237, -.718328, -.712429,
+                -.706541, -.700661, -.694791, -.688931, -.683079, -.677236,
+                -.671402, -.665576, -.659759, -.653950, -.648149, -.642356,
+                -.636570, -.630793, -.625022, -.619259, -.613503, -.607754,
+                -.602012, -.596276, -.590548, -.584825, -.579109, -.573399,
+                -.567695, -.561997, -.556305, -.550618, -.544937, -.539262,
+                -.533592, -.527926, -.522266, -.516611, -.510961, -.505315,
+                -.499674, -.494037, -.488405, -.482777,
+                -.477153, -.471533, -.465917, -.460305, -.454697, -.449092,
+                -.443491, -.437893, -.432299, -.426707, -.421119, -.415534,
+                -.409951, -.404372, -.398795, -.393221, -.387649, -.382080,
+                -.376513, -.370949, -.365387, -.359826, -.354268, -.348712,
+                -.343157, -.337604, -.332053, -.326503, -.320955, -.315408,
+                -.309863, -.304318, -.298775, -.293233, -.287692, -.282152,
+                -.276613, -.271074, -.265536, -.259999, -.254462, -.248926,
+                -.243389, -.237854, -.232318, -.226783, -.221247, -.215712,
+                -.210176, -.204641, -.199105, -.193568, -.188032, -.182495,
+                -.176957, -.171419, -.165880, -.160341, -.154800, -.149259,
+                -.143717, -.138173, -.132629, -.127083, -.121537, -.115989,
+                -.110439, -.104889, -.099336, -.093782, -.088227, -.082670,
+                -.077111, -.071550, -.065987, -.060423, -.054856, -.049288,
+                -.043717, -.038144, -.032569, -.026991, -.021411, -.015828,
+                -.010243, -.004656,  .000934,  .006527,  .012123,  .017722,
+                .023323,  .028928,  .034535,  .040146,  .045759,  .051376,
+                .056997,  .062620,  .068247,  .073877,
+                .079511,  .085149,  .090790,  .096435,  .102083,  .107736,
+                .113392,  .119052,  .124716,  .130385,  .136057,  .141734,
+                .147414,  .153100,  .158789,  .164483,  .170181,  .175884,
+                .181592,  .187304,  .193021,  .198743,  .204469,  .210201,
+                .215937,  .221678,  .227425,  .233177,  .238933,  .244696,
+                .250463,  .256236,  .262014,  .267798,  .273587,  .279382,
+                .285183,  .290989,  .296801,  .302619,  .308443,  .314273,
+                .320109,  .325951,  .331799,  .337654,  .343515,  .349382,
+                .355255,  .361135,  .367022,  .372915,  .378815,  .384721,
+                .390634,  .396554,  .402481,  .408415,  .414356,  .420304,
+                .426260,  .432222,  .438192,  .444169,  .450153,  .456145,
+                .462144,  .468151,  .474166,  .480188,  .486218,  .492256,
+                .498302,  .504356,  .510418,  .516488,  .522566,  .528653,
+                .534747,  .540850,  .546962,  .553082,  .559210,  .565347,
+                .571493,  .577648,  .583811,  .589983,  .596164,  .602355,
+                .608554,  .614762,  .620980,  .627207,  .633444,  .639689,
+                .645945,  .652210,  .658484,  .664768,
+                .671062,  .677366,  .683680,  .690004,  .696338,  .702682,
+                .709036,  .715400,  .721775,  .728160,  .734556,  .740963,
+                .747379,  .753807,  .760246,  .766695,  .773155,  .779627,
+                .786109,  .792603,  .799107,  .805624,  .812151,  .818690,
+                .825241,  .831803,  .838377,  .844962,  .851560,  .858170,
+                .864791,  .871425,  .878071,  .884729,  .891399,  .898082,
+                .904778,  .911486,  .918206,  .924940,  .931686,  .938446,
+                .945218,  .952003,  .958802,  .965614,  .972439,  .979278,
+                .986130,  .992996,  .999875, 1.006769, 1.013676, 1.020597,
+                1.027533, 1.034482, 1.041446, 1.048424, 1.055417, 1.062424,
+                1.069446, 1.076482, 1.083534, 1.090600, 1.097681, 1.104778,
+                1.111889, 1.119016, 1.126159, 1.133316, 1.140490, 1.147679,
+                1.154884, 1.162105, 1.169342, 1.176595, 1.183864, 1.191149,
+                1.198451, 1.205770, 1.213105, 1.220457, 1.227826, 1.235211,
+                1.242614, 1.250034, 1.257471, 1.264926, 1.272398, 1.279888,
+                1.287395, 1.294921, 1.302464, 1.310026, 1.317605, 1.325203,
+                1.332819, 1.340454, 1.348108, 1.355780,
+                1.363472, 1.371182, 1.378912, 1.386660, 1.394429, 1.402216,
+                1.410024, 1.417851, 1.425698, 1.433565, 1.441453, 1.449360,
+                1.457288, 1.465237, 1.473206, 1.481196, 1.489208, 1.497240,
+                1.505293, 1.513368, 1.521465, 1.529583, 1.537723, 1.545885,
+                1.554068, 1.562275, 1.570503, 1.578754, 1.587028, 1.595325,
+                1.603644, 1.611987, 1.620353, 1.628743, 1.637156, 1.645593,
+                1.654053, 1.662538, 1.671047, 1.679581, 1.688139, 1.696721,
+                1.705329, 1.713961, 1.722619, 1.731303, 1.740011, 1.748746,
+                1.757506, 1.766293, 1.775106, 1.783945, 1.792810, 1.801703,
+                1.810623, 1.819569, 1.828543, 1.837545, 1.846574, 1.855631,
+                1.864717, 1.873830, 1.882972, 1.892143, 1.901343, 1.910572,
+                1.919830, 1.929117, 1.938434, 1.947781, 1.957158, 1.966566,
+                1.976004, 1.985473, 1.994972, 2.004503, 2.014065, 2.023659,
+                2.033285, 2.042943, 2.052633, 2.062355, 2.072110, 2.081899,
+                2.091720, 2.101575, 2.111464, 2.121386, 2.131343, 2.141334,
+                2.151360, 2.161421, 2.171517, 2.181648, 2.191815, 2.202018,
+                2.212257, 2.222533, 2.232845, 2.243195,
+                2.253582, 2.264006, 2.274468, 2.284968, 2.295507, 2.306084,
+                2.316701, 2.327356, 2.338051, 2.348786, 2.359562, 2.370377,
+                2.381234, 2.392131, 2.403070, 2.414051, 2.425073, 2.436138,
+                2.447246, 2.458397, 2.469591, 2.480828, 2.492110, 2.503436,
+                2.514807, 2.526222, 2.537684, 2.549190, 2.560743, 2.572343,
+                2.583989, 2.595682, 2.607423, 2.619212, 2.631050, 2.642936,
+                2.654871, 2.666855, 2.678890, 2.690975, 2.703110, 2.715297,
+                2.727535, 2.739825, 2.752168, 2.764563, 2.777012, 2.789514,
+                2.802070, 2.814681, 2.827347, 2.840069, 2.852846, 2.865680,
+                2.878570, 2.891518, 2.904524, 2.917588, 2.930712, 2.943894,
+                2.957136, 2.970439, 2.983802, 2.997227, 3.010714, 3.024263,
+                3.037875, 3.051551, 3.065290, 3.079095, 3.092965, 3.106900,
+                3.120902, 3.134971, 3.149107, 3.163312, 3.177585, 3.191928,
+                3.206340, 3.220824, 3.235378, 3.250005, 3.264704, 3.279477,
+                3.294323, 3.309244, 3.324240, 3.339312, 3.354461, 3.369687,
+                3.384992, 3.400375, 3.415838, 3.431381, 3.447005, 3.462711,
+                3.478500, 3.494372, 3.510328, 3.526370,
+                3.542497, 3.558711, 3.575012, 3.591402, 3.607881, 3.624450,
+                3.641111, 3.657863, 3.674708, 3.691646, 3.708680, 3.725809,
+                3.743034, 3.760357, 3.777779, 3.795300, 3.812921, 3.830645,
+                3.848470, 3.866400, 3.884434, 3.902574, 3.920821, 3.939176,
+                3.957640, 3.976215, 3.994901, 4.013699, 4.032612, 4.051639,
+                4.070783, 4.090045, 4.109425, 4.128925, 4.148547, 4.168292,
+                4.188160, 4.208154, 4.228275, 4.248524, 4.268903, 4.289413,
+                4.310056, 4.330832, 4.351745, 4.372794, 4.393982, 4.415310,
+                4.436781, 4.458395, 4.480154, 4.502060, 4.524114, 4.546319,
+                4.568676, 4.591187, 4.613854, 4.636678, 4.659662, 4.682807,
+                4.706116, 4.729590, 4.753231, 4.777041, 4.801024, 4.825179,
+                4.849511, 4.874020, 4.898710, 4.923582, 4.948639, 4.973883,
+                4.999316, 5.024942, 5.050761, 5.076778, 5.102993, 5.129411,
+                5.156034, 5.182864, 5.209903, 5.237156, 5.264625, 5.292312,
+                5.320220, 5.348354, 5.376714, 5.405306, 5.434131, 5.463193,
+                5.492496, 5.522042, 5.551836, 5.581880, 5.612178, 5.642734,
+                5.673552, 5.704634, 5.735986, 5.767610,
+                5.799512, 5.831694, 5.864161, 5.896918, 5.929968, 5.963316,
+                5.996967, 6.030925, 6.065194, 6.099780, 6.134687, 6.169921,
+                6.205486, 6.241387, 6.277630, 6.314220, 6.351163, 6.388465,
+                6.426130, 6.464166, 6.502578, 6.541371, 6.580553, 6.620130,
+                6.660109, 6.700495, 6.741297, 6.782520, 6.824173, 6.866262,
+                6.908795, 6.951780, 6.995225, 7.039137, 7.083525, 7.128398,
+                7.173764, 7.219632, 7.266011, 7.312910, 7.360339, 7.408308,
+                7.456827, 7.505905, 7.555554, 7.605785, 7.656608, 7.708035,
+                7.760077, 7.812747, 7.866057, 7.920019, 7.974647, 8.029953,
+                8.085952, 8.142657, 8.200083, 8.258245, 8.317158, 8.376837,
+                8.437300, 8.498562, 8.560641, 8.623554, 8.687319, 8.751955,
+                8.817481, 8.883916, 8.951282, 9.019600, 9.088889, 9.159174,
+                9.230477, 9.302822, 9.376233, 9.450735, 9.526355, 9.603118,
+                9.681054, 9.760191, 9.840558, 9.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 };
+                
+                public double nextLandau()
+                {
+                    if (_sigma <= 0) return 0;
+                    double ranlan, x, u, v;
+                    x = _ran.nextDouble();
+                    u = 1000.*x;
+                    int i = (int) u;
+                    u -= i;
+                    if (i >= 70 && i < 800)
+                    {
+                        ranlan = f[i-1] + u*(f[i] - f[i-1]);
+                    }
+                    else if (i >= 7 && i <= 980)
+                    {
+                        ranlan =  f[i-1] + u*(f[i]-f[i-1]-0.25*(1-u)*(f[i+1]-f[i]-f[i-1]+f[i-2]));
+                    }
+                    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 _mean + _sigma*ranlan;
+                }
+                             
+}
CVSspam 0.2.8