ADMB Documentation  11.1.2439
 All Classes Files Functions Variables Typedefs Friends Defines
ranfill.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: ranfill.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2009-2012 ADMB Foundation
00006  */
00011 #include <fvar.hpp>
00012 
00013 
00014 #ifdef __TURBOC__
00015   #pragma hdrstop
00016   #include <iostream.h>
00017 #endif
00018 
00019 #ifdef __ZTC__
00020   #include <iostream.hpp>
00021 #endif
00022 
00023 #include <math.h>
00024 
00025 #define M1 714025
00026 #define IA1 1366
00027 #define IC1 150889
00028 #define RM1 (1.0/M1)
00029 #define M3 134456
00030 #define IA3 8121
00031 #define IC3 28411
00032 #define M2 243000
00033 #define IA2 4561
00034 #define IC2 51349
00035 #define RM2 (1.0/M2)
00036 
00037 double auto_rand(long int& idum, int reset);
00038 void reinitialize_auto_rand();
00039 double randn(long int& n);
00040 
00045 void reinitialize_auto_rand()
00046 {
00047   long int One=1;
00048   auto_rand(One,-One);
00049 }
00050 
00057 double auto_rand(long int& idum, int reset)
00058 {
00059   static long ix1,ix2,ix3;
00060   static float r[108];
00061   double temp;
00062   static int iff=0;
00063   int j;
00064 
00065   if (reset < 0)
00066   {
00067     iff =0;
00068     return .5;
00069   }
00070 
00071   if (idum < 0 || iff == 0)
00072   {
00073     iff=2;
00074     ix1=(IC1-(idum))%M1;
00075     ix1=ix1 % M1;
00076 
00077 
00078     ix1=(IA1*ix1+IC1);
00079     ix1=ix1 % M1;
00080 
00081     ix2=ix1 % M2;
00082     ix1=(IA1*ix1+IC1);
00083     ix1=ix1 % M1;
00084     ix3=ix1 % M3;
00085     for (j=1;j<=107;j++)
00086     {
00087       ix2=(IA2*ix2+IC2) % M2;
00088       ix1=(IA1*ix1+IC1);
00089               ix1=ix1 % M1;
00090 
00091       long int iu = (long int)(ix2 * RM2);
00092       r[j]=(ix1+iu)*RM1;
00093     }
00094     idum =6;
00095   }
00096 
00097   ix3=(IA3*ix3+IC3) % M3;
00098   ix3=ix3 % M3;
00099   ix1=(IA1*ix1+IC1) % M1;
00100   ix1=ix1 % M1;
00101   ix2=(IA2*ix2+IC2) % M2;
00102   ix2=ix2 % M2;
00103   j=1 + ((107*ix3)/M3);
00104   if (j > 107 || j < 1) cerr << " Error in random number generator\n";
00105   temp=r[j];
00106   r[j]=ix2*RM2;
00107   r[j]=(ix1+r[j]);
00108   r[j]=r[j]*RM1;
00109   return temp;
00110 }
00111 
00112 #undef M1
00113 #undef IA1
00114 #undef IC1
00115 #undef RM1
00116 #undef M2
00117 #undef IA2
00118 #undef IC2
00119 #undef RM2
00120 #undef M3
00121 #undef IA3
00122 #undef IC3
00123 
00128 double randn(long int& n)
00129 {
00130   long int nn;
00131   nn=n;
00132   double x,y;
00133   x=auto_rand(nn,1);
00134   y=auto_rand(nn,1);
00135   double u=sqrt(-2*log(x))*cos(2*PI*y);
00136   return(u);
00137 }
00138 
00143   void dvector::fill_randbi(long int& n, double p)
00144   {
00145     if ( p<0 || p>1)
00146     {
00147       cerr << "Error in dvar_vector::fill_randbi proportions of"
00148        " successes must lie between 0 and 1\n";
00149       ad_exit(1);
00150     }
00151     long int nn;
00152     nn=n;
00153     for (int i=indexmin(); i<=indexmax(); i++)
00154     {
00155       if (auto_rand(nn,1)<=p)
00156       {
00157         elem(i)=1;
00158       }
00159       else
00160       {
00161         elem(i)=0;
00162       }
00163     }
00164     reinitialize_auto_rand();
00165   }
00166 
00171   void dvector::fill_randu(long int& n)
00172   {
00173     long int nn;
00174     nn=n;
00175     for (int i=indexmin(); i<=indexmax(); i++)
00176     {
00177       elem(i)=auto_rand(nn,1);
00178     }
00179     reinitialize_auto_rand();
00180   }
00181 
00186 void dmatrix::colfill_randu(const int &j, long int &n)
00187   {
00188     long int nn;
00189     nn=n;
00190     for (int i=rowmin(); i<=rowmax(); i++)
00191     {
00192       elem(i,j)=auto_rand(nn,1);
00193     }
00194     reinitialize_auto_rand();
00195   }
00196 
00201 void dmatrix::rowfill_randu(const int& i, long int& n)
00202   {
00203     long int nn;
00204     nn=n;
00205     for (int j=colmin(); j<=colmax(); j++)
00206     {
00207       elem(i,j)=auto_rand(nn,1);
00208     }
00209     reinitialize_auto_rand();
00210   }
00211 
00216   void dvector::fill_randn(long int& n)
00217   {
00218     long int nn;
00219     nn=n;
00220     for (int i=indexmin(); i<=indexmax(); i++)
00221     {
00222       (*this)(i)=randn(nn);
00223     }
00224     reinitialize_auto_rand();
00225   }
00226 
00231   void dmatrix::fill_randn(long int& n)
00232   {
00233     long int nn=n;
00234     for (int i=rowmin(); i<=rowmax(); i++)
00235     {
00236       elem(i).fill_randn_ni(nn);
00237       nn+=2;
00238     }
00239     reinitialize_auto_rand();
00240   }
00241 
00246   void d3_array::fill_randn(long int& n)
00247   {
00248     long int nn;
00249     nn=n;
00250     for (int i=slicemin(); i<=slicemax(); i++)
00251     {
00252       elem(i).fill_randn_ni(nn);
00253       nn+=2;
00254     }
00255     reinitialize_auto_rand();
00256   }
00257 
00262   void d3_array::fill_randu(long int& n)
00263   {
00264     long int nn;
00265     nn=n;
00266     for (int i=slicemin(); i<=slicemax(); i++)
00267     {
00268       elem(i).fill_randu_ni(nn);
00269       nn+=2;
00270     }
00271     reinitialize_auto_rand();
00272   }
00273 
00278   void dmatrix::fill_randu(long int& n)
00279   {
00280     long int nn;
00281     nn=n;
00282     for (int i=rowmin(); i<=rowmax(); i++)
00283     {
00284       elem(i).fill_randu_ni(nn);
00285       nn+=2;
00286     }
00287     reinitialize_auto_rand();
00288   }
00289 
00294 void dmatrix::colfill_randn(const int &j,long int &n)
00295   {
00296     long int nn;
00297     nn=n;
00298     for (int i=rowmin(); i<=rowmax(); i++)
00299     {
00300       elem(i,j)=randn(nn);
00301     }
00302     reinitialize_auto_rand();
00303   }
00304 
00309 void dmatrix::rowfill_randn(const int& i, long int& n)
00310   {
00311     long int nn;
00312     nn=n;
00313     for (int j=colmin(); j<=colmax(); j++)
00314     {
00315       elem(i,j)=randn(nn);
00316     }
00317     reinitialize_auto_rand();
00318   }