00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include <stdio.h>
00015 #include <stdlib.h>
00016 #include <time.h>
00017 #include <math.h>
00018
00019 #ifndef __make_fast__
00020 #include "randy.h"
00021 #endif
00022
00023 # ifndef M_PI
00024 # define M_PI 3.14159265358979323846
00025 # endif
00026
00027 # ifndef RAND_MAX
00028 # define RAND_MAX 32767
00029 # endif
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040 double ran2 ( long &idum )
00041 {
00042 const long IM1=2147483563 ;
00043 const long IM2=2147483399 ;
00044 const long IA1=40014 ;
00045 const long IA2=40692 ;
00046 const long IQ1=53668 ;
00047 const long IQ2=52774 ;
00048 const long IR1=12211 ;
00049 const long IR2=3791 ;
00050 const long NTAB=32 ;
00051 const long IMM1=IM1-1 ;
00052 const long NDIV=1+IMM1/NTAB ;
00053 const double EPS=3.0e-16 ;
00054 const double RNMX=1.0-EPS ;
00055 const double AM=1.0/(double)IM1 ;
00056 static long idum2=123456789 ;
00057 static long iy=0 ;
00058 static long iv[NTAB] ;
00059
00060 long j,k ;
00061 double temp;
00062
00063 if (idum <= 0) {
00064 idum=(idum==0 ? 1 : -idum) ;
00065 idum2=idum ;
00066 for (j=NTAB+7;j>=0; j--) {
00067 k=idum/IQ1 ;
00068 idum=IA1*(idum-k*IQ1)-k*IR1 ;
00069 if (idum < 0) idum += IM1 ;
00070 if (j < NTAB) iv[j] = idum ;
00071 }
00072 iy=iv[0] ;
00073 }
00074 k=idum/IQ1;
00075 idum=IA1*(idum-k*IQ1)-k*IR1;
00076 if (idum < 0) idum += IM1;
00077 k=idum2/IQ2;
00078 idum2=IA2*(idum2-k*IQ2)-k*IR2;
00079 if (idum2 < 0) idum2 += IM2;
00080 j=iy/NDIV;
00081 iy=iv[j]-idum2;
00082 iv[j] = idum;
00083 if (iy < 1) iy += IMM1;
00084 if ((temp=AM*iy) > RNMX) return RNMX;
00085 else return temp;
00086 }
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096 long rnd ( void )
00097 {
00098 #ifdef __TEST_BOUNDS__
00099
00100 const long pole[4]={ 0L, 1L, RAND_MAX-1L, RAND_MAX } ;
00101 static long i=-1 ;
00102
00103 i++ ;
00104
00105 return ( pole[i%4] ) ;
00106
00107 #else
00108
00109 return ( rand() ) ;
00110
00111 #endif
00112 }
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122 void randy::init ( long orseed )
00123 {
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133 if ( orseed==-1 )
00134 {
00135 time_t t;
00136 t = time(&t);
00137 rseed=t ;
00138 }
00139 else
00140 {
00141 rseed=orseed ;
00142 }
00143
00144 srand( rseed ) ;
00145 }
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 long randy::give_long ( long omax )
00160 {
00161 if ( omax==1 )
00162 return ( 0 ) ;
00163 else if ( omax>0 )
00164 {
00165
00166 return ( (long)floor((double)omax * ((double)(rnd() - 1)/(double)( RAND_MAX ))) ) ;
00167 }
00168 else
00169 {
00170 fprintf ( stderr, "\n\n Higher value is lower than zero. \n" ) ;
00171 fprintf ( stderr, " randy (%s, line %d)\n",__FILE__,__LINE__) ;
00172 exit ( 1 ) ;
00173 }
00174 return( 0 ) ;
00175 }
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 long randy::give_long ( long olow, long ohigh )
00193 {
00194 long div = ohigh-olow ;
00195 if ( div>0 )
00196 return ( this->give_long( div+1 )+olow ) ;
00197 else if ( div==0 )
00198 return ( olow ) ;
00199 else
00200 {
00201 fprintf ( stderr, "\n\n Higher value is smaller than lower value. \n" ) ;
00202 fprintf ( stderr, " randy (%s, line %d)\n",__FILE__,__LINE__) ;
00203 exit ( 1 ) ;
00204 }
00205 return( 0 ) ;
00206 }
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225 long randy::give_long ( long olow, long ohigh, long oprecision )
00226 {
00227 return( this->round( this->give_long( olow, ohigh ), oprecision )) ;
00228 }
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247 double randy::give_double ( long oa, long ob )
00248 {
00249 if (( oa==0 || oa==1 ) && ( ob==0 || ob==1 ))
00250 {
00251
00252 double tmp = ((double)( rnd() ) + (double)( oa )) / ((double)( RAND_MAX ) + (double)( oa + ob )) ;
00253 return ( tmp ) ;
00254 }
00255 else
00256 {
00257 fprintf ( stderr, "\n\n Values oa and ob must be zeros or ones.. \n" ) ;
00258 fprintf ( stderr, " randy (%s, line %d)\n",__FILE__,__LINE__) ;
00259 exit ( 1 ) ;
00260 }
00261
00262 return( 0. ) ;
00263 }
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285 double randy::give_double ( double olow, double ohigh, long oa, long ob )
00286 {
00287 double div = ohigh-olow ;
00288 if ( div>NearZero )
00289 return ( this->give_double(oa, ob)*div+olow ) ;
00290 else if ( div>=0. )
00291 return ( olow ) ;
00292 else
00293 {
00294 fprintf ( stderr, "\n\n Higher value is smaller than lower value. \n" ) ;
00295 fprintf ( stderr, " randy (%s, line %d)\n",__FILE__,__LINE__) ;
00296 exit ( 1 ) ;
00297 }
00298 return( 0. ) ;
00299 }
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322 double randy::give_double ( double olow, double ohigh, double oprecision, long oa, long ob )
00323 {
00324 return( this->round( this->give_double( olow, ohigh, oa, ob ), oprecision )) ;
00325 }
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340 long randy::is_less ( double oa )
00341 {
00342 return( this->give_double(0L,1L)<oa ) ;
00343 }
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358 long randy::round ( long oa, long oprecision )
00359 {
00360 if ( oprecision==0 )
00361 {
00362 return ( oa ) ;
00363 }
00364 else if ( oprecision > 0 )
00365 {
00366 return( (long)floor( this->round( (double) oa, (double) oprecision )) ) ;
00367 }
00368 else
00369 {
00370 fprintf ( stderr, "\n\n Precision has to be non-negative. \n" ) ;
00371 fprintf ( stderr, " randy (%s, line %d)\n",__FILE__,__LINE__) ;
00372 exit ( 1 ) ;
00373 }
00374 return oa;
00375 }
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 long randy::move_to_bounds ( long oa, long olow, long ohigh )
00389 {
00390 if( oa<olow )
00391 return( olow ) ;
00392 else if( oa>ohigh )
00393 return( ohigh ) ;
00394 else
00395 return( oa ) ;
00396 }
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411 double randy::round ( double oa, double oprecision )
00412 {
00413 if ( fabs(oprecision) < NearZero )
00414 {
00415 return ( oa ) ;
00416 }
00417 else if ( oprecision > 0 )
00418 {
00419 return( ( floor(( oa/oprecision )+0.5 ))*oprecision ) ;
00420 }
00421 else
00422 {
00423 fprintf ( stderr, "\n\n Precision has to be non-negative. \n" ) ;
00424 fprintf ( stderr, " randy (%s, line %d)\n",__FILE__,__LINE__) ;
00425 exit ( 1 ) ;
00426 }
00427 return oa;
00428 }
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441 double randy::move_to_bounds ( double oa, double olow, double ohigh )
00442 {
00443 if( oa<olow )
00444 return( olow ) ;
00445 else if( oa>ohigh )
00446 return( ohigh ) ;
00447 else
00448 return( oa ) ;
00449 }
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465 long randy::give_gauss_long ( double omi, double oss )
00466 {
00467 return( (long)round( give_gauss_double( omi,oss ), 1.0 )) ;
00468 }
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484 double randy::give_gauss_double ( double omi, double oss )
00485 {
00486
00487 static bool value_stored=false ;
00488 static double value=0.0 ;
00489 double fac, v1, v2, rsq ;
00490
00491 if ( value_stored )
00492 {
00493 value_stored=false ;
00494 return value;
00495 }
00496 else
00497 {
00498
00499 do
00500 {
00501 v1=2.*give_double( 0L,1L )-1. ;
00502 v2=2.*give_double( 0L,1L )-1. ;
00503 rsq=v1*v1 + v2*v2 ;
00504 }
00505 while ( rsq>=1. || rsq==0. ) ;
00506
00507 fac=sqrt(-2. * log(rsq)/rsq * oss) ;
00508 value=v1*fac + omi ;
00509 value_stored=true ;
00510 return ( v2*fac + omi ) ;
00511 }
00512 }
00513