View Javadoc

1   /*
2    * RandLandau.java
3    *
4    * Created on June 6, 2005, 4:05 PM
5    *
6    * To change this template, choose Tools | Options and locate the template under
7    * the Source Creation and Management node. Right-click the template and choose
8    * Open. You can then make changes to the template in the Source Editor.
9    */
10  
11  package org.lcsim.math.distribution;
12  
13  
14  import static java.lang.Math.log;
15  import java.util.Random;
16  /**
17   *
18   * @author ngraf
19   */
20  public class RandLandau
21  {
22      private Random _ran = new Random();
23      private double _sigma;
24      private double _mean;
25      
26      /** Creates a new instance of RandLandau */
27      public RandLandau()
28      {
29          _mean = 0.;
30          _sigma = 1.0;
31      }
32      
33      public RandLandau(double mean, double sigma)
34      {
35          _mean = mean;
36          _sigma = sigma;
37      }
38      
39      
40      public double shoot()
41      {
42          return _mean+_sigma*transform(_ran.nextDouble());
43      }
44      
45      private double transform(double r)
46      {
47          
48          double  u = r * TABLE_MULTIPLIER;
49          int index = (int) u;
50          double du = u - index;
51          
52          // du is scaled such that the we dont have to multiply by TABLE_INTERVAL
53          // when interpolating.
54          
55          // Five cases:
56          // A) Between .070 and .800 the function is so smooth, straight
57          //	linear interpolation is adequate.
58          // B) Between .007 and .070, and between .800 and .980, quadratic
59          //    interpolation is used.  This requires the same 4 points as
60          //	a cubic spline (thus we need .006 and .981 and .982) but
61          //	the quadratic interpolation is accurate enough and quicker.
62          // C) Below .007 an asymptotic expansion for low negative lambda
63          //    (involving two logs) is used; there is a pade-style correction
64          //	factor.
65          // D) Above .980, a simple pade approximation is made (asymptotic to
66          //    1/(1-r)), but...
67          // E) the coefficients in that pade are different above r=.999.
68          
69          if ( index >= 70 && index <= 800 )
70          {		// (A)
71              
72              double f0 = inverseLandau [index];
73              double f1 = inverseLandau [index+1];
74              return f0 + du * (f1 - f0);
75              
76          }
77          else if ( index >= 7 && index <= 980 )
78          {	// (B)
79              
80              double f_1 = inverseLandau [index-1];
81              double f0  = inverseLandau [index];
82              double f1  = inverseLandau [index+1];
83              double f2  = inverseLandau [index+2];
84              
85              return f0 + du * (f1 - f0 - .25*(1-du)* (f2 -f1 - f0 + f_1) );
86              
87          }
88          else if ( index < 7 )
89          {			// (C)
90              
91              final double n0 =  0.99858950;
92              final double n1 = 34.5213058;	final double d1 = 34.1760202;
93              final double n2 = 17.0854528;	final double d2 =  4.01244582;
94              
95              double logr = log(r);
96              double x    = 1/logr;
97              double x2   = x*x;
98              
99              double pade = (n0 + n1*x + n2*x2) / (1.0 + d1*x + d2*x2);
100             
101             return ( - log( -.91893853 - logr ) -1 ) * pade;
102             
103         }
104         else if ( index <= 999 )
105         {			// (D)
106             
107             final double n0 =    1.00060006;
108             final double n1 =  263.991156;	final double d1 =  257.368075;
109             final double n2 = 4373.20068;	final double d2 = 3414.48018;
110             
111             double x = 1-r;
112             double x2   = x*x;
113             
114             return (n0 + n1*x + n2*x2) / (x * (1.0 + d1*x + d2*x2));
115             
116         }
117         else
118         { 					// (E)
119             
120             final double n0 =      1.00001538;
121             final double n1 =   6075.14119;	final double d1 =   6065.11919;
122             final double n2 = 734266.409;	final double d2 = 694021.044;
123             
124             double x = 1-r;
125             double x2   = x*x;
126             
127             return (n0 + n1*x + n2*x2) / (x * (1.0 + d1*x + d2*x2));
128             
129         }
130         
131     } // transform()
132     
133     
134     
135 //
136 // Table of values of inverse Landau, from r = .060 to .982
137 //
138     
139 // Since all these are this is static to this compilation unit only, the
140 // info is establised a priori and not at each invocation.
141     
142     static final float TABLE_INTERVAL   = .001f;
143     static final int   TABLE_END        =  982;
144     static final float TABLE_MULTIPLIER = 1.0f/TABLE_INTERVAL;
145     
146 // Here comes the big (4K bytes) table ---
147 //
148 // inverseLandau[ n ] = the inverse cdf at r = n*TABLE_INTERVAL = n/1000.
149 //
150 // Credit CERNLIB for these computations.
151 //
152 // This data is float because the algortihm does not benefit from further
153 // data accuracy.  The numbers below .006 or above .982 are moot since
154 // non-table-based methods are used below r=.007 and above .980.
155     
156     static final float[] inverseLandau = {
157         
158         0.0f,							     // .000
159                 0.0f, 	    0.0f, 	0.0f, 	    0.0f, 	0.0f, 	     // .001 - .005
160                 -2.244733f, -2.204365f, -2.168163f, -2.135219f, -2.104898f,  // .006 - .010
161                 -2.076740f, -2.050397f, -2.025605f, -2.002150f, -1.979866f,
162                 -1.958612f, -1.938275f, -1.918760f, -1.899984f, -1.881879f,  // .020
163                 -1.864385f, -1.847451f, -1.831030f, -1.815083f, -1.799574f,
164                 -1.784473f, -1.769751f, -1.755383f, -1.741346f, -1.727620f,  // .030
165                 -1.714187f, -1.701029f, -1.688130f, -1.675477f, -1.663057f,
166                 -1.650858f, -1.638868f, -1.627078f, -1.615477f, -1.604058f,  // .040
167                 -1.592811f, -1.581729f, -1.570806f, -1.560034f, -1.549407f,
168                 -1.538919f, -1.528565f, -1.518339f, -1.508237f, -1.498254f,  // .050
169                 -1.488386f, -1.478628f, -1.468976f, -1.459428f, -1.449979f,
170                 -1.440626f, -1.431365f, -1.422195f, -1.413111f, -1.404112f,  // .060
171                 -1.395194f, -1.386356f, -1.377594f, -1.368906f, -1.360291f,
172                 -1.351746f, -1.343269f, -1.334859f, -1.326512f, -1.318229f,  // .070
173                 -1.310006f, -1.301843f, -1.293737f, -1.285688f, -1.277693f,
174                 -1.269752f, -1.261863f, -1.254024f, -1.246235f, -1.238494f,  // .080
175                 -1.230800f, -1.223153f, -1.215550f, -1.207990f, -1.200474f,
176                 -1.192999f, -1.185566f, -1.178172f, -1.170817f, -1.163500f,  // .090
177                 -1.156220f, -1.148977f, -1.141770f, -1.134598f, -1.127459f,
178                 -1.120354f, -1.113282f, -1.106242f, -1.099233f, -1.092255f,  // .100
179                 
180                 -1.085306f, -1.078388f, -1.071498f, -1.064636f, -1.057802f,
181                 -1.050996f, -1.044215f, -1.037461f, -1.030733f, -1.024029f,
182                 -1.017350f, -1.010695f, -1.004064f, -.997456f,  -.990871f,
183                 -.984308f, -.977767f, -.971247f, -.964749f, -.958271f,
184                 -.951813f, -.945375f, -.938957f, -.932558f, -.926178f,
185                 -.919816f, -.913472f, -.907146f, -.900838f, -.894547f,
186                 -.888272f, -.882014f, -.875773f, -.869547f, -.863337f,
187                 -.857142f, -.850963f, -.844798f, -.838648f, -.832512f,
188                 -.826390f, -.820282f, -.814187f, -.808106f, -.802038f,
189                 -.795982f, -.789940f, -.783909f, -.777891f, -.771884f, 	// .150
190                 -.765889f, -.759906f, -.753934f, -.747973f, -.742023f,
191                 -.736084f, -.730155f, -.724237f, -.718328f, -.712429f,
192                 -.706541f, -.700661f, -.694791f, -.688931f, -.683079f,
193                 -.677236f, -.671402f, -.665576f, -.659759f, -.653950f,
194                 -.648149f, -.642356f, -.636570f, -.630793f, -.625022f,
195                 -.619259f, -.613503f, -.607754f, -.602012f, -.596276f,
196                 -.590548f, -.584825f, -.579109f, -.573399f, -.567695f,
197                 -.561997f, -.556305f, -.550618f, -.544937f, -.539262f,
198                 -.533592f, -.527926f, -.522266f, -.516611f, -.510961f,
199                 -.505315f, -.499674f, -.494037f, -.488405f, -.482777f,	// .200
200                 
201                 -.477153f, -.471533f, -.465917f, -.460305f, -.454697f,
202                 -.449092f, -.443491f, -.437893f, -.432299f, -.426707f,
203                 -.421119f, -.415534f, -.409951f, -.404372f, -.398795f,
204                 -.393221f, -.387649f, -.382080f, -.376513f, -.370949f,
205                 -.365387f, -.359826f, -.354268f, -.348712f, -.343157f,
206                 -.337604f, -.332053f, -.326503f, -.320955f, -.315408f,
207                 -.309863f, -.304318f, -.298775f, -.293233f, -.287692f,
208                 -.282152f, -.276613f, -.271074f, -.265536f, -.259999f,
209                 -.254462f, -.248926f, -.243389f, -.237854f, -.232318f,
210                 -.226783f, -.221247f, -.215712f, -.210176f, -.204641f, 	// .250
211                 -.199105f, -.193568f, -.188032f, -.182495f, -.176957f,
212                 -.171419f, -.165880f, -.160341f, -.154800f, -.149259f,
213                 -.143717f, -.138173f, -.132629f, -.127083f, -.121537f,
214                 -.115989f, -.110439f, -.104889f, -.099336f, -.093782f,
215                 -.088227f, -.082670f, -.077111f, -.071550f, -.065987f,
216                 -.060423f, -.054856f, -.049288f, -.043717f, -.038144f,
217                 -.032569f, -.026991f, -.021411f, -.015828f, -.010243f,
218                 -.004656f,  .000934f,  .006527f,  .012123f,  .017722f,
219                 .023323f, .028928f,  .034535f,  .040146f,  .045759f,
220                 .051376f, .056997f,  .062620f,  .068247f,  .073877f,	// .300
221                 
222                 .079511f,  .085149f,  .090790f,  .096435f,  .102083f,
223                 .107736f,  .113392f,  .119052f,  .124716f,  .130385f,
224                 .136057f,  .141734f,  .147414f,  .153100f,  .158789f,
225                 .164483f,  .170181f,  .175884f,  .181592f,  .187304f,
226                 .193021f,  .198743f,  .204469f,  .210201f,  .215937f,
227                 .221678f,  .227425f,  .233177f,  .238933f,  .244696f,
228                 .250463f,  .256236f,  .262014f,  .267798f,  .273587f,
229                 .279382f,  .285183f,  .290989f,  .296801f,  .302619f,
230                 .308443f,  .314273f,  .320109f,  .325951f,  .331799f,
231                 .337654f,  .343515f,  .349382f,  .355255f,  .361135f,	// .350
232                 .367022f,  .372915f,  .378815f,  .384721f,  .390634f,
233                 .396554f,  .402481f,  .408415f,  .414356f,  .420304f,
234                 .426260f,  .432222f,  .438192f,  .444169f,  .450153f,
235                 .456145f,  .462144f,  .468151f,  .474166f,  .480188f,
236                 .486218f,  .492256f,  .498302f,  .504356f,  .510418f,
237                 .516488f,  .522566f,  .528653f,  .534747f,  .540850f,
238                 .546962f,  .553082f,  .559210f,  .565347f,  .571493f,
239                 .577648f,  .583811f,  .589983f,  .596164f,  .602355f,
240                 .608554f,  .614762f,  .620980f,  .627207f,  .633444f,
241                 .639689f,  .645945f,  .652210f,  .658484f,  .664768f,	// .400
242                 
243                 .671062f,  .677366f,  .683680f,  .690004f,  .696338f,
244                 .702682f,  .709036f,  .715400f,  .721775f,  .728160f,
245                 .734556f,  .740963f,  .747379f,  .753807f,  .760246f,
246                 .766695f,  .773155f,  .779627f,  .786109f,  .792603f,
247                 .799107f,  .805624f,  .812151f,  .818690f,  .825241f,
248                 .831803f,  .838377f,  .844962f,  .851560f,  .858170f,
249                 .864791f,  .871425f,  .878071f,  .884729f,  .891399f,
250                 .898082f,  .904778f,  .911486f,  .918206f,  .924940f,
251                 .931686f,  .938446f,  .945218f,  .952003f,  .958802f,
252                 .965614f,  .972439f,  .979278f,  .986130f,  .992996f, 	// .450
253                 .999875f,  1.006769f, 1.013676f, 1.020597f, 1.027533f,
254                 1.034482f, 1.041446f, 1.048424f, 1.055417f, 1.062424f,
255                 1.069446f, 1.076482f, 1.083534f, 1.090600f, 1.097681f,
256                 1.104778f, 1.111889f, 1.119016f, 1.126159f, 1.133316f,
257                 1.140490f, 1.147679f, 1.154884f, 1.162105f, 1.169342f,
258                 1.176595f, 1.183864f, 1.191149f, 1.198451f, 1.205770f,
259                 1.213105f, 1.220457f, 1.227826f, 1.235211f, 1.242614f,
260                 1.250034f, 1.257471f, 1.264926f, 1.272398f, 1.279888f,
261                 1.287395f, 1.294921f, 1.302464f, 1.310026f, 1.317605f,
262                 1.325203f, 1.332819f, 1.340454f, 1.348108f, 1.355780f,	// .500
263                 
264                 1.363472f, 1.371182f, 1.378912f, 1.386660f, 1.394429f,
265                 1.402216f, 1.410024f, 1.417851f, 1.425698f, 1.433565f,
266                 1.441453f, 1.449360f, 1.457288f, 1.465237f, 1.473206f,
267                 1.481196f, 1.489208f, 1.497240f, 1.505293f, 1.513368f,
268                 1.521465f, 1.529583f, 1.537723f, 1.545885f, 1.554068f,
269                 1.562275f, 1.570503f, 1.578754f, 1.587028f, 1.595325f,
270                 1.603644f, 1.611987f, 1.620353f, 1.628743f, 1.637156f,
271                 1.645593f, 1.654053f, 1.662538f, 1.671047f, 1.679581f,
272                 1.688139f, 1.696721f, 1.705329f, 1.713961f, 1.722619f,
273                 1.731303f, 1.740011f, 1.748746f, 1.757506f, 1.766293f, 	// .550
274                 1.775106f, 1.783945f, 1.792810f, 1.801703f, 1.810623f,
275                 1.819569f, 1.828543f, 1.837545f, 1.846574f, 1.855631f,
276                 1.864717f, 1.873830f, 1.882972f, 1.892143f, 1.901343f,
277                 1.910572f, 1.919830f, 1.929117f, 1.938434f, 1.947781f,
278                 1.957158f, 1.966566f, 1.976004f, 1.985473f, 1.994972f,
279                 2.004503f, 2.014065f, 2.023659f, 2.033285f, 2.042943f,
280                 2.052633f, 2.062355f, 2.072110f, 2.081899f, 2.091720f,
281                 2.101575f, 2.111464f, 2.121386f, 2.131343f, 2.141334f,
282                 2.151360f, 2.161421f, 2.171517f, 2.181648f, 2.191815f,
283                 2.202018f, 2.212257f, 2.222533f, 2.232845f, 2.243195f,	// .600
284                 
285                 2.253582f, 2.264006f, 2.274468f, 2.284968f, 2.295507f,
286                 2.306084f, 2.316701f, 2.327356f, 2.338051f, 2.348786f,
287                 2.359562f, 2.370377f, 2.381234f, 2.392131f, 2.403070f,
288                 2.414051f, 2.425073f, 2.436138f, 2.447246f, 2.458397f,
289                 2.469591f, 2.480828f, 2.492110f, 2.503436f, 2.514807f,
290                 2.526222f, 2.537684f, 2.549190f, 2.560743f, 2.572343f,
291                 2.583989f, 2.595682f, 2.607423f, 2.619212f, 2.631050f,
292                 2.642936f, 2.654871f, 2.666855f, 2.678890f, 2.690975f,
293                 2.703110f, 2.715297f, 2.727535f, 2.739825f, 2.752168f,
294                 2.764563f, 2.777012f, 2.789514f, 2.802070f, 2.814681f,	// .650
295                 2.827347f, 2.840069f, 2.852846f, 2.865680f, 2.878570f,
296                 2.891518f, 2.904524f, 2.917588f, 2.930712f, 2.943894f,
297                 2.957136f, 2.970439f, 2.983802f, 2.997227f, 3.010714f,
298                 3.024263f, 3.037875f, 3.051551f, 3.065290f, 3.079095f,
299                 3.092965f, 3.106900f, 3.120902f, 3.134971f, 3.149107f,
300                 3.163312f, 3.177585f, 3.191928f, 3.206340f, 3.220824f,
301                 3.235378f, 3.250005f, 3.264704f, 3.279477f, 3.294323f,
302                 3.309244f, 3.324240f, 3.339312f, 3.354461f, 3.369687f,
303                 3.384992f, 3.400375f, 3.415838f, 3.431381f, 3.447005f,
304                 3.462711f, 3.478500f, 3.494372f, 3.510328f, 3.526370f,	// .700
305                 
306                 3.542497f, 3.558711f, 3.575012f, 3.591402f, 3.607881f,
307                 3.624450f, 3.641111f, 3.657863f, 3.674708f, 3.691646f,
308                 3.708680f, 3.725809f, 3.743034f, 3.760357f, 3.777779f,
309                 3.795300f, 3.812921f, 3.830645f, 3.848470f, 3.866400f,
310                 3.884434f, 3.902574f, 3.920821f, 3.939176f, 3.957640f,
311                 3.976215f, 3.994901f, 4.013699f, 4.032612f, 4.051639f,
312                 4.070783f, 4.090045f, 4.109425f, 4.128925f, 4.148547f,
313                 4.168292f, 4.188160f, 4.208154f, 4.228275f, 4.248524f,
314                 4.268903f, 4.289413f, 4.310056f, 4.330832f, 4.351745f,
315                 4.372794f, 4.393982f, 4.415310f, 4.436781f, 4.458395f,
316                 4.480154f, 4.502060f, 4.524114f, 4.546319f, 4.568676f,	// .750
317                 4.591187f, 4.613854f, 4.636678f, 4.659662f, 4.682807f,
318                 4.706116f, 4.729590f, 4.753231f, 4.777041f, 4.801024f,
319                 4.825179f, 4.849511f, 4.874020f, 4.898710f, 4.923582f,
320                 4.948639f, 4.973883f, 4.999316f, 5.024942f, 5.050761f,
321                 5.076778f, 5.102993f, 5.129411f, 5.156034f, 5.182864f,
322                 5.209903f, 5.237156f, 5.264625f, 5.292312f, 5.320220f,
323                 5.348354f, 5.376714f, 5.405306f, 5.434131f, 5.463193f,
324                 5.492496f, 5.522042f, 5.551836f, 5.581880f, 5.612178f,
325                 5.642734f, 5.673552f, 5.704634f, 5.735986f, 5.767610f,	// .800
326                 
327                 5.799512f, 5.831694f, 5.864161f, 5.896918f, 5.929968f,
328                 5.963316f, 5.996967f, 6.030925f, 6.065194f, 6.099780f,
329                 6.134687f, 6.169921f, 6.205486f, 6.241387f, 6.277630f,
330                 6.314220f, 6.351163f, 6.388465f, 6.426130f, 6.464166f,
331                 6.502578f, 6.541371f, 6.580553f, 6.620130f, 6.660109f,
332                 6.700495f, 6.741297f, 6.782520f, 6.824173f, 6.866262f,
333                 6.908795f, 6.951780f, 6.995225f, 7.039137f, 7.083525f,
334                 7.128398f, 7.173764f, 7.219632f, 7.266011f, 7.312910f,
335                 7.360339f, 7.408308f, 7.456827f, 7.505905f, 7.555554f,
336                 7.605785f, 7.656608f, 7.708035f, 7.760077f, 7.812747f, 	// .850
337                 7.866057f, 7.920019f, 7.974647f, 8.029953f, 8.085952f,
338                 8.142657f, 8.200083f, 8.258245f, 8.317158f, 8.376837f,
339                 8.437300f, 8.498562f, 8.560641f, 8.623554f, 8.687319f,
340                 8.751955f, 8.817481f, 8.883916f, 8.951282f, 9.019600f,
341                 9.088889f, 9.159174f, 9.230477f, 9.302822f, 9.376233f,
342                 9.450735f, 9.526355f, 9.603118f, 9.681054f, 9.760191f,
343                 9.840558f,  9.922186f, 10.005107f, 10.089353f, 10.174959f,
344                 10.261958f, 10.350389f, 10.440287f, 10.531693f, 10.624646f,
345                 10.719188f, 10.815362f, 10.913214f, 11.012789f, 11.114137f,
346                 11.217307f, 11.322352f, 11.429325f, 11.538283f, 11.649285f,	// .900
347                 
348                 11.762390f, 11.877664f, 11.995170f, 12.114979f, 12.237161f,
349                 12.361791f, 12.488946f, 12.618708f, 12.751161f, 12.886394f,
350                 13.024498f, 13.165570f, 13.309711f, 13.457026f, 13.607625f,
351                 13.761625f, 13.919145f, 14.080314f, 14.245263f, 14.414134f,
352                 14.587072f, 14.764233f, 14.945778f, 15.131877f, 15.322712f,
353                 15.518470f, 15.719353f, 15.925570f, 16.137345f, 16.354912f,
354                 16.578520f, 16.808433f, 17.044929f, 17.288305f, 17.538873f,
355                 17.796967f, 18.062943f, 18.337176f, 18.620068f, 18.912049f,
356                 19.213574f, 19.525133f, 19.847249f, 20.180480f, 20.525429f,
357                 20.882738f, 21.253102f, 21.637266f, 22.036036f, 22.450278f, 	// .950
358                 22.880933f, 23.329017f, 23.795634f, 24.281981f, 24.789364f,
359                 25.319207f, 25.873062f, 26.452634f, 27.059789f, 27.696581f,     // .960
360                 28.365274f, 29.068370f, 29.808638f, 30.589157f, 31.413354f,
361                 32.285060f, 33.208568f, 34.188705f, 35.230920f, 36.341388f,     // .970
362                 37.527131f, 38.796172f, 40.157721f, 41.622399f, 43.202525f,
363                 44.912465f, 46.769077f, 48.792279f, 51.005773f, 53.437996f,     // .980
364                 56.123356f, 59.103894f, 					// .982
365                 
366     };  // End of the inverseLandau table
367     
368     
369 // root method
370     
371     static final double[] f = {
372         0       , 0       , 0       ,0        ,0        ,-2.244733,
373                 -2.204365,-2.168163,-2.135219,-2.104898,-2.076740,-2.050397,
374                 -2.025605,-2.002150,-1.979866,-1.958612,-1.938275,-1.918760,
375                 -1.899984,-1.881879,-1.864385,-1.847451,-1.831030,-1.815083,
376                 -1.799574,-1.784473,-1.769751,-1.755383,-1.741346,-1.727620,
377                 -1.714187,-1.701029,-1.688130,-1.675477,-1.663057,-1.650858,
378                 -1.638868,-1.627078,-1.615477,-1.604058,-1.592811,-1.581729,
379                 -1.570806,-1.560034,-1.549407,-1.538919,-1.528565,-1.518339,
380                 -1.508237,-1.498254,-1.488386,-1.478628,-1.468976,-1.459428,
381                 -1.449979,-1.440626,-1.431365,-1.422195,-1.413111,-1.404112,
382                 -1.395194,-1.386356,-1.377594,-1.368906,-1.360291,-1.351746,
383                 -1.343269,-1.334859,-1.326512,-1.318229,-1.310006,-1.301843,
384                 -1.293737,-1.285688,-1.277693,-1.269752,-1.261863,-1.254024,
385                 -1.246235,-1.238494,-1.230800,-1.223153,-1.215550,-1.207990,
386                 -1.200474,-1.192999,-1.185566,-1.178172,-1.170817,-1.163500,
387                 -1.156220,-1.148977,-1.141770,-1.134598,-1.127459,-1.120354,
388                 -1.113282,-1.106242,-1.099233,-1.092255,
389                 -1.085306,-1.078388,-1.071498,-1.064636,-1.057802,-1.050996,
390                 -1.044215,-1.037461,-1.030733,-1.024029,-1.017350,-1.010695,
391                 -1.004064, -.997456, -.990871, -.984308, -.977767, -.971247,
392                 -.964749, -.958271, -.951813, -.945375, -.938957, -.932558,
393                 -.926178, -.919816, -.913472, -.907146, -.900838, -.894547,
394                 -.888272, -.882014, -.875773, -.869547, -.863337, -.857142,
395                 -.850963, -.844798, -.838648, -.832512, -.826390, -.820282,
396                 -.814187, -.808106, -.802038, -.795982, -.789940, -.783909,
397                 -.777891, -.771884, -.765889, -.759906, -.753934, -.747973,
398                 -.742023, -.736084, -.730155, -.724237, -.718328, -.712429,
399                 -.706541, -.700661, -.694791, -.688931, -.683079, -.677236,
400                 -.671402, -.665576, -.659759, -.653950, -.648149, -.642356,
401                 -.636570, -.630793, -.625022, -.619259, -.613503, -.607754,
402                 -.602012, -.596276, -.590548, -.584825, -.579109, -.573399,
403                 -.567695, -.561997, -.556305, -.550618, -.544937, -.539262,
404                 -.533592, -.527926, -.522266, -.516611, -.510961, -.505315,
405                 -.499674, -.494037, -.488405, -.482777,
406                 -.477153, -.471533, -.465917, -.460305, -.454697, -.449092,
407                 -.443491, -.437893, -.432299, -.426707, -.421119, -.415534,
408                 -.409951, -.404372, -.398795, -.393221, -.387649, -.382080,
409                 -.376513, -.370949, -.365387, -.359826, -.354268, -.348712,
410                 -.343157, -.337604, -.332053, -.326503, -.320955, -.315408,
411                 -.309863, -.304318, -.298775, -.293233, -.287692, -.282152,
412                 -.276613, -.271074, -.265536, -.259999, -.254462, -.248926,
413                 -.243389, -.237854, -.232318, -.226783, -.221247, -.215712,
414                 -.210176, -.204641, -.199105, -.193568, -.188032, -.182495,
415                 -.176957, -.171419, -.165880, -.160341, -.154800, -.149259,
416                 -.143717, -.138173, -.132629, -.127083, -.121537, -.115989,
417                 -.110439, -.104889, -.099336, -.093782, -.088227, -.082670,
418                 -.077111, -.071550, -.065987, -.060423, -.054856, -.049288,
419                 -.043717, -.038144, -.032569, -.026991, -.021411, -.015828,
420                 -.010243, -.004656,  .000934,  .006527,  .012123,  .017722,
421                 .023323,  .028928,  .034535,  .040146,  .045759,  .051376,
422                 .056997,  .062620,  .068247,  .073877,
423                 .079511,  .085149,  .090790,  .096435,  .102083,  .107736,
424                 .113392,  .119052,  .124716,  .130385,  .136057,  .141734,
425                 .147414,  .153100,  .158789,  .164483,  .170181,  .175884,
426                 .181592,  .187304,  .193021,  .198743,  .204469,  .210201,
427                 .215937,  .221678,  .227425,  .233177,  .238933,  .244696,
428                 .250463,  .256236,  .262014,  .267798,  .273587,  .279382,
429                 .285183,  .290989,  .296801,  .302619,  .308443,  .314273,
430                 .320109,  .325951,  .331799,  .337654,  .343515,  .349382,
431                 .355255,  .361135,  .367022,  .372915,  .378815,  .384721,
432                 .390634,  .396554,  .402481,  .408415,  .414356,  .420304,
433                 .426260,  .432222,  .438192,  .444169,  .450153,  .456145,
434                 .462144,  .468151,  .474166,  .480188,  .486218,  .492256,
435                 .498302,  .504356,  .510418,  .516488,  .522566,  .528653,
436                 .534747,  .540850,  .546962,  .553082,  .559210,  .565347,
437                 .571493,  .577648,  .583811,  .589983,  .596164,  .602355,
438                 .608554,  .614762,  .620980,  .627207,  .633444,  .639689,
439                 .645945,  .652210,  .658484,  .664768,
440                 .671062,  .677366,  .683680,  .690004,  .696338,  .702682,
441                 .709036,  .715400,  .721775,  .728160,  .734556,  .740963,
442                 .747379,  .753807,  .760246,  .766695,  .773155,  .779627,
443                 .786109,  .792603,  .799107,  .805624,  .812151,  .818690,
444                 .825241,  .831803,  .838377,  .844962,  .851560,  .858170,
445                 .864791,  .871425,  .878071,  .884729,  .891399,  .898082,
446                 .904778,  .911486,  .918206,  .924940,  .931686,  .938446,
447                 .945218,  .952003,  .958802,  .965614,  .972439,  .979278,
448                 .986130,  .992996,  .999875, 1.006769, 1.013676, 1.020597,
449                 1.027533, 1.034482, 1.041446, 1.048424, 1.055417, 1.062424,
450                 1.069446, 1.076482, 1.083534, 1.090600, 1.097681, 1.104778,
451                 1.111889, 1.119016, 1.126159, 1.133316, 1.140490, 1.147679,
452                 1.154884, 1.162105, 1.169342, 1.176595, 1.183864, 1.191149,
453                 1.198451, 1.205770, 1.213105, 1.220457, 1.227826, 1.235211,
454                 1.242614, 1.250034, 1.257471, 1.264926, 1.272398, 1.279888,
455                 1.287395, 1.294921, 1.302464, 1.310026, 1.317605, 1.325203,
456                 1.332819, 1.340454, 1.348108, 1.355780,
457                 1.363472, 1.371182, 1.378912, 1.386660, 1.394429, 1.402216,
458                 1.410024, 1.417851, 1.425698, 1.433565, 1.441453, 1.449360,
459                 1.457288, 1.465237, 1.473206, 1.481196, 1.489208, 1.497240,
460                 1.505293, 1.513368, 1.521465, 1.529583, 1.537723, 1.545885,
461                 1.554068, 1.562275, 1.570503, 1.578754, 1.587028, 1.595325,
462                 1.603644, 1.611987, 1.620353, 1.628743, 1.637156, 1.645593,
463                 1.654053, 1.662538, 1.671047, 1.679581, 1.688139, 1.696721,
464                 1.705329, 1.713961, 1.722619, 1.731303, 1.740011, 1.748746,
465                 1.757506, 1.766293, 1.775106, 1.783945, 1.792810, 1.801703,
466                 1.810623, 1.819569, 1.828543, 1.837545, 1.846574, 1.855631,
467                 1.864717, 1.873830, 1.882972, 1.892143, 1.901343, 1.910572,
468                 1.919830, 1.929117, 1.938434, 1.947781, 1.957158, 1.966566,
469                 1.976004, 1.985473, 1.994972, 2.004503, 2.014065, 2.023659,
470                 2.033285, 2.042943, 2.052633, 2.062355, 2.072110, 2.081899,
471                 2.091720, 2.101575, 2.111464, 2.121386, 2.131343, 2.141334,
472                 2.151360, 2.161421, 2.171517, 2.181648, 2.191815, 2.202018,
473                 2.212257, 2.222533, 2.232845, 2.243195,
474                 2.253582, 2.264006, 2.274468, 2.284968, 2.295507, 2.306084,
475                 2.316701, 2.327356, 2.338051, 2.348786, 2.359562, 2.370377,
476                 2.381234, 2.392131, 2.403070, 2.414051, 2.425073, 2.436138,
477                 2.447246, 2.458397, 2.469591, 2.480828, 2.492110, 2.503436,
478                 2.514807, 2.526222, 2.537684, 2.549190, 2.560743, 2.572343,
479                 2.583989, 2.595682, 2.607423, 2.619212, 2.631050, 2.642936,
480                 2.654871, 2.666855, 2.678890, 2.690975, 2.703110, 2.715297,
481                 2.727535, 2.739825, 2.752168, 2.764563, 2.777012, 2.789514,
482                 2.802070, 2.814681, 2.827347, 2.840069, 2.852846, 2.865680,
483                 2.878570, 2.891518, 2.904524, 2.917588, 2.930712, 2.943894,
484                 2.957136, 2.970439, 2.983802, 2.997227, 3.010714, 3.024263,
485                 3.037875, 3.051551, 3.065290, 3.079095, 3.092965, 3.106900,
486                 3.120902, 3.134971, 3.149107, 3.163312, 3.177585, 3.191928,
487                 3.206340, 3.220824, 3.235378, 3.250005, 3.264704, 3.279477,
488                 3.294323, 3.309244, 3.324240, 3.339312, 3.354461, 3.369687,
489                 3.384992, 3.400375, 3.415838, 3.431381, 3.447005, 3.462711,
490                 3.478500, 3.494372, 3.510328, 3.526370,
491                 3.542497, 3.558711, 3.575012, 3.591402, 3.607881, 3.624450,
492                 3.641111, 3.657863, 3.674708, 3.691646, 3.708680, 3.725809,
493                 3.743034, 3.760357, 3.777779, 3.795300, 3.812921, 3.830645,
494                 3.848470, 3.866400, 3.884434, 3.902574, 3.920821, 3.939176,
495                 3.957640, 3.976215, 3.994901, 4.013699, 4.032612, 4.051639,
496                 4.070783, 4.090045, 4.109425, 4.128925, 4.148547, 4.168292,
497                 4.188160, 4.208154, 4.228275, 4.248524, 4.268903, 4.289413,
498                 4.310056, 4.330832, 4.351745, 4.372794, 4.393982, 4.415310,
499                 4.436781, 4.458395, 4.480154, 4.502060, 4.524114, 4.546319,
500                 4.568676, 4.591187, 4.613854, 4.636678, 4.659662, 4.682807,
501                 4.706116, 4.729590, 4.753231, 4.777041, 4.801024, 4.825179,
502                 4.849511, 4.874020, 4.898710, 4.923582, 4.948639, 4.973883,
503                 4.999316, 5.024942, 5.050761, 5.076778, 5.102993, 5.129411,
504                 5.156034, 5.182864, 5.209903, 5.237156, 5.264625, 5.292312,
505                 5.320220, 5.348354, 5.376714, 5.405306, 5.434131, 5.463193,
506                 5.492496, 5.522042, 5.551836, 5.581880, 5.612178, 5.642734,
507                 5.673552, 5.704634, 5.735986, 5.767610,
508                 5.799512, 5.831694, 5.864161, 5.896918, 5.929968, 5.963316,
509                 5.996967, 6.030925, 6.065194, 6.099780, 6.134687, 6.169921,
510                 6.205486, 6.241387, 6.277630, 6.314220, 6.351163, 6.388465,
511                 6.426130, 6.464166, 6.502578, 6.541371, 6.580553, 6.620130,
512                 6.660109, 6.700495, 6.741297, 6.782520, 6.824173, 6.866262,
513                 6.908795, 6.951780, 6.995225, 7.039137, 7.083525, 7.128398,
514                 7.173764, 7.219632, 7.266011, 7.312910, 7.360339, 7.408308,
515                 7.456827, 7.505905, 7.555554, 7.605785, 7.656608, 7.708035,
516                 7.760077, 7.812747, 7.866057, 7.920019, 7.974647, 8.029953,
517                 8.085952, 8.142657, 8.200083, 8.258245, 8.317158, 8.376837,
518                 8.437300, 8.498562, 8.560641, 8.623554, 8.687319, 8.751955,
519                 8.817481, 8.883916, 8.951282, 9.019600, 9.088889, 9.159174,
520                 9.230477, 9.302822, 9.376233, 9.450735, 9.526355, 9.603118,
521                 9.681054, 9.760191, 9.840558, 9.922186,10.005107,10.089353,
522                 10.174959,10.261958,10.350389,10.440287,10.531693,10.624646,
523                 10.719188,10.815362,10.913214,11.012789,11.114137,11.217307,
524                 11.322352,11.429325,11.538283,11.649285,
525                 11.762390,11.877664,11.995170,12.114979,12.237161,12.361791,
526                 12.488946,12.618708,12.751161,12.886394,13.024498,13.165570,
527                 13.309711,13.457026,13.607625,13.761625,13.919145,14.080314,
528                 14.245263,14.414134,14.587072,14.764233,14.945778,15.131877,
529                 15.322712,15.518470,15.719353,15.925570,16.137345,16.354912,
530                 16.578520,16.808433,17.044929,17.288305,17.538873,17.796967,
531                 18.062943,18.337176,18.620068,18.912049,19.213574,19.525133,
532                 19.847249,20.180480,20.525429,20.882738,21.253102,21.637266,
533                 22.036036,22.450278,22.880933,23.329017,23.795634,24.281981,
534                 24.789364,25.319207,25.873062,26.452634,27.059789,27.696581,
535                 28.365274,29.068370,29.808638,30.589157,31.413354,32.285060,
536                 33.208568,34.188705,35.230920,36.341388,37.527131,38.796172,
537                 40.157721,41.622399,43.202525,44.912465,46.769077,48.792279,
538                 51.005773,53.437996,56.123356,59.103894 };
539                 
540                 public double nextLandau()
541                 {
542                     if (_sigma <= 0) return 0;
543                     double ranlan, x, u, v;
544                     x = _ran.nextDouble();
545                     u = 1000.*x;
546                     int i = (int) u;
547                     u -= i;
548                     if (i >= 70 && i < 800)
549                     {
550                         ranlan = f[i-1] + u*(f[i] - f[i-1]);
551                     }
552                     else if (i >= 7 && i <= 980)
553                     {
554                         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]));
555                     }
556                     else if (i < 7)
557                     {
558                         v = log(x);
559                         u = 1/v;
560                         ranlan = ((0.99858950+(3.45213058E1+1.70854528E1*u)*u)/
561                                 (1         +(3.41760202E1+4.01244582  *u)*u))*
562                                 (log(-0.91893853-v)-1);
563                     }
564                     else
565                     {
566                         u = 1-x;
567                         v = u*u;
568                         if (x <= 0.999)
569                         {
570                             ranlan = (1.00060006+2.63991156E2*u+4.37320068E3*v)/
571                                     ((1         +2.57368075E2*u+3.41448018E3*v)*u);
572                         }
573                         else
574                         {
575                             ranlan = (1.00001538+6.07514119E3*u+7.34266409E5*v)/
576                                     ((1         +6.06511919E3*u+6.94021044E5*v)*u);
577                         }
578                     }
579                     return _mean + _sigma*ranlan;
580                 }
581                              
582 }