home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 41 / IOPROG_41.ISO / soft / c++ / NUMCPP11.ZIP / deviate.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1999-04-17  |  11.2 KB  |  445 lines
  1. //===================================================================
  2. // deviate.cpp
  3. //
  4. // Version 1.1
  5. //
  6. // Written by:
  7. //   Brent Worden
  8. //   WordenWare
  9. //   email:  Brent@Worden.org
  10. //
  11. // Copyright (c) 1998-1999 WordenWare
  12. //
  13. // Created:  August 28, 1998
  14. // Revised:  April 10, 1999
  15. //===================================================================
  16.  
  17. #include <cmath>
  18.  
  19. #include "contdist.h"
  20. #include "deviate.h"
  21. #include "mathx.h"
  22.  
  23. NUM_BEGIN
  24.  
  25. static const long   IA   = 16087L;
  26. static const long   IM   = 2147483647L;
  27. static const double AM   = 1.0 / (double)IM;
  28. static const long   IQ   = 127773L;
  29. static const long   IR   = 2836L;
  30. static const long   MASK = 123459876L;
  31.  
  32. NUMERICS_EXPORT double ran0(long *idum)
  33. {
  34.     long k;
  35.     double ans;
  36.     
  37.     *idum ^= MASK;
  38.     k = (*idum)/IQ;
  39.     *idum = IA*(*idum-k*IQ)-IR*k;
  40.     if(*idum < 0) *idum += IM;
  41.     ans = AM*(*idum);
  42.     *idum ^= MASK;
  43.     
  44.     return ans;
  45. }
  46.  
  47. static const long   NTAB = 32L;
  48. static const double NDIV = 1.0 + (double)(IM-1.0) / (double)NTAB;
  49. static const double RNMX = 1.0 - NUMERICS_MAX_ERROR;
  50.  
  51. NUMERICS_EXPORT double ran1(long *idum)
  52. {
  53.     int j;
  54.     long k;
  55.     static long iy=0;
  56.     static long iv[NTAB];
  57.     double temp;
  58.     
  59.     if(*idum <= 0 || !iy){
  60.         if(-(*idum) < 1) *idum = 1;
  61.         else *idum = -(*idum);
  62.         for(j = NTAB+7; j >= 0; j--){
  63.             k = (*idum) / IQ;
  64.             *idum = IA * (*idum - k * IQ) - IR * k;
  65.             if(*idum < 0) *idum += IM;
  66.             if(j < NTAB) iv[j] = *idum;
  67.         }
  68.         iy=iv[0];
  69.     }
  70.     k=(*idum)/IQ;
  71.     *idum=IA*(*idum-k*IQ)-IR*k;
  72.     if(*idum<0) *idum += IM;
  73.     j=(int)(iy/NDIV);
  74.     iy=iv[j];
  75.     iv[j]= *idum;
  76.     if((temp=AM*iy) > RNMX) return RNMX;
  77.     else return temp;
  78. }
  79.  
  80. static const long   IA1   = 40014L;
  81. static const long   IA2   = 40692L;
  82. static const long   IM1   = 2147483563L;
  83. static const double AM1   = 1.0 / (double)IM1;
  84. static const long   IM2   = 2147483399L;
  85. static const long   IMM1  = IM1 - 1;
  86. static const long   IQ1   = 53668L;
  87. static const long   IQ2   = 52774L;
  88. static const long   IR1   = 12211L;
  89. static const long   IR2   = 3791;
  90. static const double NDIV1 = 1.0+ (double)IMM1 / (double)NTAB;
  91.  
  92. NUMERICS_EXPORT double ran2(long *idum)
  93. {
  94.     int j;
  95.     long k;
  96.     static long idum2 = 123456789L;
  97.     static long iy=0;
  98.     static long iv[NTAB];
  99.     double temp;
  100.     
  101.     if(*idum <= 0){
  102.         if(-(*idum) < 1) *idum = 1;
  103.         else *idum = -(*idum);
  104.         idum2=*idum;
  105.         for(j=NTAB+7;j>=0;j--){
  106.             k=(*idum)/IQ1;
  107.             *idum =IA1*(*idum-k*IQ1)-k*IR1;
  108.             if(*idum < 0) *idum += IM1;
  109.             if(j < NTAB) iv[j] = *idum;
  110.         }
  111.         iy=iv[0];
  112.     }
  113.     k=(*idum)/IQ1;
  114.     *idum=IA1*(*idum-k*IQ1)-k*IR1;
  115.     if(*idum < 0) *idum += IM1;
  116.     k = idum2/IQ2;
  117.     idum2=IA2*(idum2-k*IQ2)-k*IR2;
  118.     if(idum2 < 0) idum2 += IM2;
  119.     j = (int)(iy/NDIV1);
  120.     iy=iv[j]-idum2;
  121.     iv[j]=*idum;
  122.     if(iy < 1)iy +=IMM1;
  123.     if((temp=AM1*iy) > RNMX) return RNMX;
  124.     else return temp;
  125. }
  126.  
  127. static const long   MBIG  = 1000000000L;
  128. static const double FAC   = 1.0 / (double)MBIG;
  129. static const long   MSEED = 161803398L;
  130. static const long   MZ    = 0L;
  131.  
  132. NUMERICS_EXPORT double ran3(long *idum)
  133. {
  134.     static int inext, inextp;
  135.     static long ma[56];
  136.     static int iff = 0;
  137.     long mj, mk;
  138.     int i, ii, k;
  139.     
  140.     if(*idum < 0 || iff == 0){
  141.         iff = 1;
  142.         mj = MSEED - (*idum < 0 ? -*idum : *idum);
  143.         mj %= MBIG;
  144.         ma[55] = mj;
  145.         mk = 1;
  146.         for(i = 1; i <= 54; i++){
  147.             ii = (21 * i) % 55;
  148.             ma[ii] = mk;
  149.             mk = mj - mk;
  150.             if(mk < MZ) mk += MBIG;
  151.             mj = ma[ii];
  152.         }
  153.         for(k = 1; k <= 4; k++){
  154.             for(i = 1; i <= 55; i++){
  155.                 ma[i] -= ma[1+(i+30) % 55];
  156.                 if(ma[i] < MZ) ma[i] += MBIG;
  157.             }
  158.         }
  159.         inext = 0;
  160.         inextp = 31;
  161.         *idum = 1;
  162.     }
  163.     if(++inext == 56) inext = 1;
  164.     if(++inextp == 56) inextp = 1;
  165.     mj = ma[inext] - ma[inextp];
  166.     if(mj < MZ) mj += MBIG;
  167.     ma[inext] = mj;
  168.     
  169.     return mj * FAC;
  170. }
  171.  
  172. NUMERICS_EXPORT double berdev(double p, long *idum)
  173. {
  174.     if(RAND(idum) <= p) return 1.0;
  175.  
  176.     return 0.0;
  177. }
  178.  
  179. NUMERICS_EXPORT double betadev(double a, double b, long *idum)
  180. {
  181.     double aa, bb, c, l, m, s, y, x;
  182.     
  183.     if(a <= 1.0 || b <= 1.0) return betav(RAND(idum), a, b);
  184.     aa = a-1.0; bb = b-1.0; c = aa + bb; l = c * log(c); m = aa / c; s = .5 / sqrt(c);
  185.     do {
  186.         do {
  187.             y = gasdev(0.0, 1.0, idum);
  188.             x = s * y + m;
  189.         } while(x < 0.0 || x > 1.0);
  190.     } while(log(RAND(idum)) > aa * log(x/aa) + bb*log((1.0-x)/bb) + l + .5*y*y);
  191.     
  192.     return x;
  193. }
  194.  
  195. NUMERICS_EXPORT double bnldev(int n, double pp, long *idum)
  196. {
  197.     int i;
  198.     static int nold = -1;
  199.     double am, em, g, angle, p, bnl, sq, t, y;
  200.     static double pold = -1.0, pc, plog, pclog, en, oldg;
  201.     
  202.     p = (pp <= .5 ? pp : 1.0 - pp);
  203.     am = n * p;
  204.     if(n < 25){
  205.         bnl = 0.0;
  206.         for(i = 0; i < n; i++) bnl += berdev(p, idum);
  207.     } else if(am < 1.0){
  208.         g = exp(-am);
  209.         t = 1.0;
  210.         for(i = 0; i <= n; i++){
  211.             t *= RAND(idum);
  212.             if(t < g) break;
  213.         }
  214.         bnl = (i <= n ? i : n);
  215.     } else {
  216.         if(n != nold){
  217.             en = n;
  218.             oldg = gammln(en + 1.0);
  219.             nold = n;
  220.         }
  221.         if(p != pold){
  222.             pc = 1.0 - p;
  223.             plog = log(p);
  224.             pclog = log(pc);
  225.             pold = p;
  226.         }
  227.         sq = sqrt(2.0 * am * pc);
  228.         do {
  229.             do {
  230.                 angle = NUMERICS_PI * RAND(idum);
  231.                 y = tan(angle);
  232.                 em = sq * y + am;
  233.             } while(em < 0.0 || em >= (en + 1.0));
  234.             em = floor(em);
  235.             t = 1.2 * sq * (1.0 + y * y) * exp(oldg - gammln(em + 1.0) - gammln(en - em + 1.0) + em * plog + (en - em) * pclog);
  236.         } while(RAND(idum) > t);
  237.         bnl = em;
  238.     }
  239.     if(p != pp) bnl = n - bnl;
  240.     
  241.     return bnl;
  242. }
  243.  
  244. NUMERICS_EXPORT double caudev(double m, double s, long *idum)
  245. {
  246.     return m + s*tan(NUMERICS_PI*(RAND(idum)-.5));
  247. }
  248.  
  249. NUMERICS_EXPORT double chidev(double v, long *idum)
  250. {
  251.     return gamdev(v * .5, 2.0, idum);
  252. }
  253.  
  254. NUMERICS_EXPORT double dexpdev(double m, double s, long *idum)
  255. {
  256.     double u = RAND(idum), sn = sign(2.0 * u - 1.0);
  257.     
  258.     return m - s * sn * log(1.0-sn*(2.0*u-1.0));
  259. }
  260.  
  261. NUMERICS_EXPORT double expdev(double b, long *idum)
  262. {
  263.     double dum;
  264.     
  265.     do
  266.         dum=RAND(idum);
  267.     while(dum == 0.0);
  268.     return -log(dum)*b;
  269. }
  270.  
  271. NUMERICS_EXPORT double fdev(double dfn, double dfd, long *idum)
  272. {
  273.     double xnum = chidev(dfn, idum) / dfn, xden;
  274.     
  275.     do {
  276.         xden = chidev(dfd, idum) / dfd;
  277.     } while(xden <= 9.999999999998E-39*xnum);
  278.     
  279.     return (xnum / xden);
  280. }
  281.  
  282. NUMERICS_EXPORT double gamdev(double a, double b, long *idum)
  283. {
  284.     double u, bb, p, x, y;
  285.     
  286.     if(a < 1.0){
  287.         do {
  288.             u = RAND(idum);
  289.             bb = (NUMERICS_E + a) / NUMERICS_E;
  290.             if((p = bb * u) > .1){
  291.                 x = -log((bb-p) / a);
  292.                 if(RAND(idum) <= pow(x, a - 1.0)) return b * x;
  293.             } else {
  294.                 x = pow(p, 1.0/a);
  295.                 if(RAND(idum) <= exp(-x)) return b * x;
  296.             }
  297.         } while(1);
  298.     } else {
  299.         do {
  300.             y = expdev(1.0, idum);
  301.         } while(RAND(idum) > pow(y / exp(y-1.0), a - 1.0));
  302.         
  303.         return a * b * y;
  304.     }
  305. }
  306.  
  307. NUMERICS_EXPORT double gasdev(double m, double v, long *idum)
  308. {
  309.     static int iset=0;
  310.     static double gset;
  311.     double fac, rsq, v1, v2;
  312.     
  313.     if(iset==0){
  314.         do{
  315.             v1=2.0*RAND(idum)-1.0;
  316.             v2=2.0*RAND(idum)-1.0;
  317.             rsq=v1*v1+v2*v2;
  318.         } while(rsq>=1.0 || rsq==0.0);
  319.         fac=sqrt(-2.0*log(rsq)/rsq);
  320.         gset=v1*fac;
  321.         iset=1;
  322.         return v2*fac*sqrt(v)+m;
  323.     } else {
  324.         iset = 0;
  325.         return gset*sqrt(v)+m;
  326.     }
  327. }
  328.  
  329. NUMERICS_EXPORT double logdev(double m, double b, long *idum)
  330. {
  331.     return m-b*log(1.0/RAND(idum)-1.0);
  332. }
  333.  
  334. NUMERICS_EXPORT double nchidev(double df, double xnonc, long *idum)
  335. {
  336.     double gas = gasdev(sqrt(xnonc), 1.0, idum);
  337.     
  338.     return (chidev(df-1.0, idum) + gas * gas);
  339. }
  340.  
  341. NUMERICS_EXPORT double nfdev(double dfn, double dfd, double xnonc, long *idum)
  342. {
  343.     double xden, xnum;
  344.     
  345.     do {
  346.         xnum = nchidev(dfn, xnonc, idum) / dfn;
  347.         xden = chidev(dfd, idum) / dfd;
  348.     } while(xden < 1.0 && xnum > NUMERICS_FLOAT_MAX * xden);
  349.     
  350.     return (xnum/xden);
  351. }
  352.  
  353. NUMERICS_EXPORT double poisdev(double xm, long *idum)
  354. {
  355.     static double sq, alxm, g, oldm = -1.0;
  356.     double em, t, y;
  357.     
  358.     if(xm < 12.0){
  359.         if(xm != oldm){
  360.             oldm = xm;
  361.             g = exp(-xm);
  362.         }
  363.         em = -1;
  364.         t = 1.0;
  365.         do {
  366.             ++em;
  367.             t *= RAND(idum);
  368.         } while(t > g);
  369.     } else {
  370.         if(xm != oldm){
  371.             oldm = xm;
  372.             sq = sqrt(2.0 * xm);
  373.             alxm = log(xm);
  374.             g = xm * alxm - gammln(xm + 1.0);
  375.         }
  376.         do {
  377.             do {
  378.                 y = tan(NUMERICS_PI * RAND(idum));
  379.                 em = sq * y + xm;
  380.             } while (em < 0.0);
  381.             em = floor(em);
  382.             t = .9 * (1.0 + y*y) * exp(em *alxm-gammln(em+1.0)-g);
  383.         } while(RAND(idum) > t);
  384.     }
  385.     
  386.     return em;
  387. }
  388.  
  389. double funt(double x, double v)
  390. // internal function needed for tdev(v, idum)
  391. {
  392.     return pow(1.0 + x * x / v, -(v + 1.0) * .5);
  393. }
  394.  
  395. NUMERICS_EXPORT double tdev(double v, long *idum)
  396. {
  397.     double c, u1, u2, u3, x, s;
  398.     
  399.     c = (v + 1.0) * .5 / (sqrt(NUMERICS_PI * v) * v * .5);
  400.     u1 = RAND(idum);
  401.     if(u1 >= sqrt(2.0 / NUMERICS_PI) || u1 >= 2.0 * c){
  402.         do {
  403.             u2 = RAND(idum);
  404.             if(u2 <= .3622520694){
  405.                 if(u2 <= .05300969080){
  406.                     s = sign(7.840088159 * u2 - .2078);
  407.                     x = s * (fabs(7.840088159 * u2 - .2078) + 1.7922);
  408.                     u3 = RAND(idum);
  409.                     if(.2 * u3 <= funt(x, v) - 1.0 + fabs(x) * .5) return x;
  410.                 } else {
  411.                     x = 11.5909050257 * u2 - 2.406629332;
  412.                     u3 = RAND(idum);
  413.                     if(.13528 * u3 <= funt(x, v) - 1.0 + fabs(x) * .5) return x;
  414.                 }
  415.             } else {
  416.                 x = 1.0 / (1.0680176321 - 1.5680176321 * u2);
  417.                 u3 = RAND(idum);
  418.                 if(u3 <= x * x * funt(x, v)) return x;
  419.             }
  420.         } while(1);
  421.     }
  422.     
  423.     return 2.0 * (RAND(idum) + RAND(idum) - 1.0);
  424. }
  425.  
  426. NUMERICS_EXPORT double unifdev(double lo, double hi, long *idum)
  427. {
  428.     return (lo + RAND(idum) * (hi - lo));
  429. }
  430.  
  431. NUMERICS_EXPORT double weibdev(double a, double b, long *idum)
  432. {
  433.     return b*pow(log(1.0/RAND(idum)),1.0/a);
  434. }
  435.  
  436. NUM_END
  437.  
  438. //===================================================================
  439. // Revision History
  440. //
  441. // Version 1.0 - 08/28/1998 - New.
  442. // Version 1.1 - 04/10/1999 - Added Numerics namespace.
  443. //===================================================================
  444.  
  445.