home *** CD-ROM | disk | FTP | other *** search
/ Frozen Fish 1: Amiga / FrozenFish-Apr94.iso / bbs / alib / d3xx / d386 / xlispstat.lha / XLispStat / src2.lzh / XLisp-Stat / distribs.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-10-04  |  17.7 KB  |  651 lines
  1. /* distributions - Basic continuous probability distributions          */
  2. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  3. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  4. /* You may give out copies of this software; for conditions see the    */
  5. /* file COPYING included with this distribution.                       */
  6.  
  7. #include "xlisp.h"
  8. #include "osdef.h"
  9. #ifdef ANSI
  10. #include "xlproto.h"
  11. #include "xlsproto.h"
  12. #else
  13. #include "xlfun.h"
  14. #include "xlsfun.h"
  15. #endif ANSI
  16.  
  17. /* forward declarations */
  18. #ifdef ANSI
  19. LVAL normalcdf(void),betacdf(void),gammacdf(void),chisqcdf(void),tcdf(void),
  20.      fcdf(void),cauchycdf(void),loggamma(void),bnormcdf(void),quant(int),
  21.      normalquant(void),cauchyquant(void),betaquant(void),gammaquant(void),
  22.      chisqquant(void),tquant(void),fquant(void),dens(int),normal_dens(void),
  23.      cauchy_dens(void),beta_dens(void),gamma_dens(void),chisq_dens(void),
  24.      t_dens(void),f_dens(void),contrand(int),xsuniformrand(void),
  25.      xsnormalrand(void),xscauchyrand(void),xsgammarand(void),xschisqrand(void),
  26.      xstrand(void),xsbetarand(void),xsfrand(void);
  27. double getXarg(void),inverse_cdf(double,double,double,int),logbeta(double,double),
  28.        betadens(double,double,double),gammadens(double,double),
  29.        tdens(double,double),gammainit(double,double),
  30.        betainit(double,double,double),getposdouble(void),unirand(void),
  31.        normrand(void),cauchyrand(void),gammarand(double),chisqrand(double),
  32.        trand(double),betarand(double,double),frand(double,double);
  33. void getbetaargs(double *,double *,int *,int *),getgxtarg(double *),
  34.      getfargs(double *,double *,double *,int *,int *),check_one(LVAL,double),
  35.      checkstrict(double);
  36. #else
  37. LVAL normalcdf(),betacdf(),gammacdf(),chisqcdf(),tcdf(),
  38.      fcdf(),cauchycdf(),loggamma(),bnormcdf(),quant(),
  39.      normalquant(),cauchyquant(),betaquant(),gammaquant(),
  40.      chisqquant(),tquant(),fquant(),dens(int),normal_dens(),
  41.      cauchy_dens(),beta_dens(),gamma_dens(),chisq_dens(),
  42.      t_dens(),f_dens(),contrand(),xsuniformrand(),
  43.      xsnormalrand(),xscauchyrand(),xsgammarand(),xschisqrand(),
  44.      xstrand(),xsbetarand(),xsfrand();
  45. double getXarg(),inverse_cdf(),logbeta(),
  46.        betadens(),gammadens(),
  47.        tdens(),gammainit(),betainit(),
  48.        getposdouble(),unirand(),
  49.        normrand(),cauchyrand(),gammarand(),chisqrand(),
  50.        trand(),betarand(),frand();
  51. void getbetaargs(),getgxtarg(),
  52.      getfargs(),check_one(),
  53.      checkstrict();
  54. #endif ANSI
  55.  
  56. # ifndef PI
  57. # define PI 3.14159265358979323846
  58. # endif PI
  59. # ifndef nil
  60. # define nil 0L
  61. # endif
  62.  
  63. /***************************************************************************/
  64. /**                                                                       **/
  65. /**                         Argument Readers                              **/
  66. /**                                                                       **/
  67. /***************************************************************************/
  68.  
  69. static void getbetaargs(pa, pb, pia, pib)
  70.      double *pa, *pb;
  71.      int *pia, *pib;
  72. {
  73.   LVAL La, Lb;
  74.   double da, db;
  75.   
  76.   La = xlgetarg(); da = makedouble(La);
  77.   Lb = xlgetarg(); db = makedouble(Lb);
  78.   xllastarg();
  79.   if (da <= 0.0) xlerror("alpha is too small", La);
  80.   if (db <= 0.0) xlerror("beta is too small", Lb);
  81.  
  82.   if (pa != nil) *pa = da; 
  83.   if (pb != nil) *pb = db;
  84.   if (pia != nil) *pia = floor(da);
  85.   if (pib != nil) *pib = floor(db);
  86. }
  87.  
  88. static void getgxtarg(pa)
  89.      double *pa;
  90. {
  91.   LVAL La;
  92.   double da;
  93.   
  94.   La = xlgetarg(); da = makedouble(La);
  95.   xllastarg();
  96.   if (da <= 0.0) xlerror("alpha is too small", La);
  97.   if (pa != nil) *pa = da; 
  98. }
  99.  
  100. static void getfargs(px, pa, pb, pia, pib)
  101.      double *px, *pa, *pb;
  102.      int *pia, *pib;
  103. {
  104.   LVAL La, Lb;
  105.   double da, db;
  106.   
  107.   La = xlgetarg(); da = makedouble(La);
  108.   Lb = xlgetarg(); db = makedouble(Lb);
  109.   xllastarg();
  110.   if (da <= 0.0) xlerror("alpha is too small", La);
  111.   if (db <= 0.0) xlerror("beta is too small", Lb);
  112.   da = 0.5 * da; db = 0.5 * db; 
  113.   
  114.   if (px != nil) *px = db / (db + da * *px);
  115.   if (pa != nil) *pa = da; 
  116.   if (pb != nil) *pb = db;
  117.   if (pia != nil) *pia = floor(da);
  118.   if (pib != nil) *pib = floor(db);
  119. }
  120.  
  121. static double getXarg() { return(makedouble(xlgetarg())); }
  122.  
  123. static void check_one(p, dp)
  124.      LVAL p;
  125.      double dp;
  126. {
  127.   if (dp < 0.0 || dp >= 1.0)
  128.     xlerror("probability not between 0 and 1", p);
  129. }
  130.  
  131. /***************************************************************************/
  132. /**                                                                       **/
  133. /**                          Numerical Cdf's                              **/
  134. /**                                                                       **/
  135. /***************************************************************************/
  136.  
  137. LOCAL LVAL normalcdf()
  138. {
  139.   double dx, dp;
  140.   
  141.   dx = getXarg();
  142.   normbase(&dx, &dp);
  143.   return(cvflonum((FLOTYPE) dp));
  144. }
  145.  
  146. LOCAL LVAL betacdf()
  147. {
  148.   double dx, da, db, dp;
  149.   int ia, ib;
  150.   
  151.   dx = getXarg();
  152.   getbetaargs(&da, &db, &ia, &ib);
  153.   betabase(&dx, &da, &db, &ia, &ib, &dp);
  154.   return(cvflonum((FLOTYPE) dp));
  155. }
  156.  
  157. LOCAL LVAL gammacdf()
  158. {
  159.   double dx, da, dp;
  160.  
  161.   dx = getXarg();
  162.   getgxtarg(&da);
  163.   gammabase(&dx, &da, &dp);
  164.   return(cvflonum((FLOTYPE) dp));
  165. }
  166.  
  167. LOCAL LVAL chisqcdf()
  168. {
  169.   double dx, da, dp;
  170.  
  171.   dx = getXarg();
  172.   getgxtarg(&da);
  173.   da = 0.5 * da; dx = 0.5 * dx;
  174.   gammabase(&dx, &da, &dp);
  175.   return(cvflonum((FLOTYPE) dp));
  176. }
  177.  
  178. LOCAL LVAL tcdf()
  179. {
  180.   double dx, da, dp;
  181.  
  182.   dx = getXarg();
  183.   getgxtarg(&da);
  184.   studentbase(&dx, &da, &dp);
  185.   return(cvflonum((FLOTYPE) dp));
  186. }
  187.  
  188. LOCAL LVAL fcdf()
  189. {
  190.   double dx, da, db, dp;
  191.   int ia, ib;
  192.  
  193.   dx = getXarg();
  194.   getfargs(&dx, &da, &db, &ia, &ib);
  195.   betabase(&dx, &db, &da, &ib, &ia, &dp);
  196.   dp = 1.0 - dp;
  197.   return(cvflonum((FLOTYPE) dp));
  198. }
  199.  
  200. LOCAL LVAL cauchycdf()
  201. {
  202.   double dx, dp;
  203.   
  204.   dx = getXarg();
  205.   dp = (atan(dx) + PI / 2) / PI;
  206.   return(cvflonum((FLOTYPE) dp));
  207. }
  208.  
  209. /* log-gamma function does not really belong, but... */
  210. LOCAL LVAL loggamma()
  211. {
  212.   LVAL x;
  213.   double dx, dp;
  214.  
  215.   x = xlgetarg();
  216.   dx = makedouble(x);
  217.   if (dx <= 0) xlerror("non positive argument", x);
  218.   dp = gamma(dx);
  219.   return(cvflonum((FLOTYPE) dp));
  220. }
  221.  
  222. /* bivariate normal cdf */
  223. LOCAL LVAL bnormcdf()
  224. {
  225.   LVAL R;
  226.   double x, y, r;
  227.   x = makedouble(xlgetarg());
  228.   y = makedouble(xlgetarg());
  229.   R = xlgetarg(); r = makedouble(R);
  230.   xllastarg();
  231.   
  232.   if (r < -1 || r > 1) xlerror("correlation out of range", R);
  233.   return(cvflonum((FLOTYPE) bivnor(-x, -y, r)));
  234. }
  235.  
  236. /* recursive distribution functions */
  237. LVAL xsrnormalcdf()
  238.     { return(recursive_subr_map_elements(normalcdf, xsrnormalcdf)); }
  239. LVAL xsrbetacdf() 
  240.     { return(recursive_subr_map_elements(betacdf, xsrbetacdf));   }
  241. LVAL xsrgammacdf()
  242.     { return(recursive_subr_map_elements(gammacdf, xsrgammacdf));  }
  243. LVAL xsrchisqcdf()
  244.     { return(recursive_subr_map_elements(chisqcdf, xsrchisqcdf));  }
  245. LVAL xsrtcdf()
  246.     { return(recursive_subr_map_elements(tcdf, xsrtcdf));      }
  247. LVAL xsrfcdf()
  248.     { return(recursive_subr_map_elements(fcdf, xsrfcdf));      }
  249. LVAL xsrcauchycdf()
  250.     { return(recursive_subr_map_elements(cauchycdf, xsrcauchycdf)); }
  251. LVAL xsrloggamma() 
  252.     { return(recursive_subr_map_elements(loggamma, xsrloggamma));  }
  253. LVAL xsrbnormcdf() 
  254.     { return(recursive_subr_map_elements(bnormcdf, xsrbnormcdf));  }
  255.  
  256. /***************************************************************************/
  257. /**                                                                       **/
  258. /**                    Numerical Quantile Functions                       **/
  259. /**                                                                       **/
  260. /***************************************************************************/
  261.  
  262. LOCAL LVAL quant(dist)
  263.      int dist;
  264. {
  265.   LVAL p;
  266.   double dp, da, db, dx;
  267.   int ia, ib;
  268.  
  269.   p = xlgetarg(); dp = makedouble(p);
  270.   if (dp < 0.0 || dp > 1.0) xlerror("probability out of range", p);
  271.  
  272.   switch (dist) {
  273.   case 'N': xllastarg(); checkstrict(dp); dx = ppnd(dp, &ia); break;
  274.   case 'C': xllastarg(); checkstrict(dp); dx = tan(PI * (dp - 0.5)); break;
  275.   case 'B': getbetaargs(&da, &db, &ia, &ib);
  276.             check_one(p, dp);
  277.             dx = ppbeta(dp, da, db, &ia);
  278.             break;
  279.   case 'G': getgxtarg(&da);
  280.             db = 0.0;
  281.             check_one(p, dp);
  282.             dx = ppgamma(dp, da, &ia);
  283.             break;
  284.   case 'X': getgxtarg(&da);
  285.             da = 0.5 * da; db = 0.0;
  286.             check_one(p, dp);
  287.             dx = 2.0 * ppgamma(dp, da, &ia);
  288.             break;
  289.   case 'T': getgxtarg(&da);
  290.             db = 0.0;
  291.             checkstrict(dp); 
  292.             dx = ppstudent(dp, da, &ia);
  293.             break;
  294.   case 'F': getfargs(nil, &da, &db, &ia, &ib);
  295.             check_one(p, dp);
  296.             if (dp == 0.0) dx = 0.0;
  297.             else {           
  298.               dp = 1.0 - dp;
  299.               dx = ppbeta(dp, db, da, &ia);
  300.               dx = db * (1.0 / dx - 1.0) / da;
  301.             }
  302.             break;
  303.   default:  xlfail("unknown distribution");
  304.   }
  305.   return(cvflonum((FLOTYPE) dx));
  306. }
  307.  
  308. LOCAL LVAL normalquant() { return(quant('N')); }
  309. LOCAL LVAL cauchyquant() { return(quant('C')); }
  310. LOCAL LVAL betaquant()   { return(quant('B')); }
  311. LOCAL LVAL gammaquant()  { return(quant('G')); }
  312. LOCAL LVAL chisqquant()  { return(quant('X')); }
  313. LOCAL LVAL tquant()      { return(quant('T')); }
  314. LOCAL LVAL fquant()      { return(quant('F')); }
  315.  
  316. /* recursive quantile functions */
  317. LVAL xsrnormalquant()  
  318.     { return(recursive_subr_map_elements(normalquant, xsrnormalquant)); }
  319. LVAL xsrcauchyquant()
  320.     { return(recursive_subr_map_elements(cauchyquant, xsrcauchyquant)); }
  321. LVAL xsrbetaquant()
  322.     { return(recursive_subr_map_elements(betaquant, xsrbetaquant)); }
  323. LVAL xsrgammaquant()
  324.     { return(recursive_subr_map_elements(gammaquant, xsrgammaquant)); }
  325. LVAL xsrchisqquant()
  326.     { return(recursive_subr_map_elements(chisqquant, xsrchisqquant)); }
  327. LVAL xsrtquant()
  328.     { return(recursive_subr_map_elements(tquant, xsrtquant)); }
  329. LVAL xsrfquant()
  330.     { return(recursive_subr_map_elements(fquant, xsrfquant)); }
  331.  
  332. /***************************************************************************/
  333. /**                                                                       **/
  334. /**                    Numerical Density Functions                       **/
  335. /**                                                                       **/
  336. /***************************************************************************/
  337.  
  338. LOCAL LVAL dens(dist)
  339.      int dist;
  340. {
  341.   LVAL x;
  342.   double dx, da, db, dens;
  343.  
  344.   x = xlgetarg(); dx = makedouble(x);
  345.  
  346.   switch (dist) {
  347.   case 'N': xllastarg(); dens = exp(- 0.5 * dx * dx) / sqrt(2.0 * PI); break;
  348.   case 'B': getbetaargs(&da, &db, nil, nil);
  349.             dens = betadens(dx, da, db);
  350.             break;
  351.   case 'G': getgxtarg(&da);
  352.             dens = gammadens(dx, da);
  353.             break;
  354.   case 'X': getgxtarg(&da);
  355.             da = 0.5 * da; dx = 0.5 * dx;
  356.             dens = 0.5 * gammadens(dx, da);
  357.             break;
  358.   case 'T': getgxtarg(&da);
  359.             dens = tdens(dx, da);
  360.             break;
  361.   case 'F': getbetaargs(&da, &db, nil, nil);
  362.             if (dx <= 0.0) dens = 0.0;
  363.             else {
  364.               dens = exp(0.5 * da * log(da) + 0.5 * db *log(db)
  365.                          + (0.5 * da - 1.0) * log(dx)
  366.                          - logbeta(0.5 * da, 0.5 * db)
  367.                          - 0.5 * (da + db) * log(db + da * dx));
  368.             }
  369.             break;
  370.   case 'C': xllastarg(); dens = tdens(dx, 1.0); break;
  371.   default:  xlfail(" unknown distribution");
  372.   }
  373.   
  374.   return(cvflonum((FLOTYPE) dens));
  375. }
  376.  
  377. /* density functions */
  378. LOCAL LVAL normal_dens() { return(dens('N')); }
  379. LOCAL LVAL cauchy_dens() { return(dens('C')); }
  380. LOCAL LVAL beta_dens()   { return(dens('B')); }
  381. LOCAL LVAL gamma_dens()  { return(dens('G')); }
  382. LOCAL LVAL chisq_dens()  { return(dens('X')); }
  383. LOCAL LVAL t_dens()      { return(dens('T')); }
  384. LOCAL LVAL f_dens()      { return(dens('F')); }
  385.  
  386. /* recursive density functions */
  387. LVAL xsrnormaldens()  
  388.     { return(recursive_subr_map_elements(normal_dens, xsrnormaldens)); }
  389. LVAL xsrcauchydens()
  390.     { return(recursive_subr_map_elements(cauchy_dens, xsrcauchydens)); }
  391. LVAL xsrbetadens()
  392.     { return(recursive_subr_map_elements(beta_dens, xsrbetadens)); }
  393. LVAL xsrgammadens()
  394.     { return(recursive_subr_map_elements(gamma_dens, xsrgammadens)); }
  395. LVAL xsrchisqdens()
  396.     { return(recursive_subr_map_elements(chisq_dens, xsrchisqdens)); }
  397. LVAL xsrtdens()
  398.     { return(recursive_subr_map_elements(t_dens, xsrtdens)); }
  399. LVAL xsrfdens()
  400.     { return(recursive_subr_map_elements(f_dens, xsrfdens)); }
  401.  
  402. LOCAL double logbeta(a, b)
  403.         double a, b;
  404. {
  405.   static double da = 0.0, db = 0.0, lbeta = 0.0;
  406.   
  407.   if (a != da || b != db) { /* cache most recent call */
  408.     da = a; db = b;
  409.     lbeta = gamma(da) + gamma(db) - gamma(da + db);
  410.   }
  411.   return(lbeta);
  412. }
  413.  
  414. LOCAL double betadens(x, a, b)
  415.         double x, a, b;
  416. {
  417.   double dens;
  418.   
  419.   if (x <= 0.0 || x >= 1.0) dens = 0.0;
  420.   else {
  421.     dens = exp(log(x) * (a - 1) + log(1 - x) * (b - 1) - logbeta(a, b));
  422.   }
  423.   return(dens);
  424. }
  425.  
  426. LOCAL double gammadens(x, a)
  427.         double x, a;
  428. {
  429.   double dens;
  430.   if (x <= 0.0) dens = 0.0;
  431.   else {
  432.     dens = exp(log(x) * (a - 1) - x - gamma(a));
  433.   }
  434.   return(dens);
  435. }
  436.  
  437. LOCAL double tdens(x, a)
  438.         double x, a;
  439. {
  440.   double dens;
  441.   
  442.   dens = (1.0 / sqrt(a * PI)) 
  443.        * exp(gamma(0.5 * (a + 1)) - gamma(0.5 * a) 
  444.              - 0.5 * (a + 1) * log(1.0 + x * x / a));
  445.   return(dens);
  446. }
  447.  
  448. LOCAL void checkstrict(dp)
  449.      double dp;
  450. {
  451.   if (dp <= 0.0 || dp >= 1.0)
  452.     xlfail("probbility not strictly between 0 and 1");
  453. }
  454.  
  455. LOCAL double getposdouble()
  456. {
  457.   LVAL x;
  458.   double dx;
  459.   
  460.   x = xlgetarg();
  461.   dx = makedouble(x);
  462.   if (dx <= 0.0) xlerror("not a positive number", x);
  463.   return(dx);
  464. }
  465.  
  466. LOCAL double unirand()
  467. {
  468.   double u;
  469.   do {
  470.     u = uni();
  471.   } while ((u <= 0.0) || (u >= 1.0));
  472.   return(u);
  473. }
  474.  
  475. LOCAL double normrand()
  476. {
  477.   double x, y, u, u1, v;
  478.   static double c = -1.0;
  479.   
  480.   if (c < 0.0) c = sqrt(2.0 / exp(1.0));
  481.   
  482.   /* ratio of uniforms with linear pretest */
  483.   do {
  484.     u = unirand();
  485.     u1 = unirand();
  486.     v = c * (2 * u1 - 1);
  487.     x = v / u;
  488.     y = x * x / 4.0;
  489.   } while(y > (1 - u) && y > - log(u));
  490.   return(x);
  491. }
  492.  
  493. LOCAL double cauchyrand()
  494. {
  495.   double u1, u2, v1, v2;
  496.   
  497.   /* ratio of uniforms on half disk */
  498.   do {
  499.     u1 = unirand();
  500.     u2 = unirand();
  501.     v1 = 2.0 * u1 - 1.0;
  502.     v2 = u2;
  503.   } while(v1 * v1 + v2 * v2 > 1.0);
  504.   return(v1 / v2);
  505. }
  506.  
  507. LOCAL double gammarand(a)
  508.      double a;
  509. {
  510.   double x, u0, u1, u2, v, w, c, c1, c2, c3, c4, c5;
  511.   static double e = -1.0;
  512.   int done;
  513.   
  514.   if (e < 0.0) e = exp(1.0);
  515.   
  516.   if (a < 1.0) {
  517.     /* Ahrens and Dieter algorithm */
  518.     done = FALSE;
  519.     c = (a + e) / e;
  520.     do {
  521.       u0 = unirand();
  522.       u1 = unirand();
  523.       v = c * u0;
  524.       if (v <= 1.0) {
  525.         x = exp(log(v) / a);
  526.         if (u1 <= exp(-x)) done = TRUE;
  527.       }
  528.       else {
  529.         x = -log((c - v) / a);
  530.         if (x > 0.0 && u1 < exp((a - 1.0) * log(x))) done = TRUE;
  531.       }
  532.     } while(! done);
  533.   }
  534.   else if (a == 1.0) x = -log(unirand());
  535.   else {
  536.     /* Cheng and Feast algorithm */
  537.     c1 = a - 1.0;
  538.     c2 = (a - 1.0 / (6.0 * a)) / c1;
  539.     c3 = 2.0 / c1;
  540.     c4 = 2.0 / (a - 1.0) + 2.0;
  541.     c5 = 1.0 / sqrt(a);
  542.     do {
  543.       do {
  544.         u1 = unirand();
  545.         u2 = unirand();
  546.         if (a > 2.5) u1 = u2 + c5 * (1.0 - 1.86 * u1);
  547.       } while (u1 <= 0.0 || u1 >= 1.0);
  548.       w = c2 * u2 / u1;
  549.     } while ((c3 * u1 + w + 1.0/w) > c4 && (c3 * log(u1) - log(w) + w) > 1.0);
  550.     x = c1 * w;
  551.   }
  552.   return(x);
  553. }
  554.  
  555. LOCAL double chisqrand(df)
  556.         double df;
  557. {
  558.   return(2.0 * gammarand(df / 2.0));
  559. }
  560.  
  561. LOCAL double trand(df)
  562.      double df;
  563. {
  564.   return(normrand() / sqrt(chisqrand(df) / df));
  565. }
  566.  
  567. LOCAL double betarand(a, b)
  568.      double a, b;
  569. {
  570.   double x, y;
  571.   
  572.   x = gammarand(a);
  573.   y = gammarand(b);
  574.   return(x / (x + y));
  575. }
  576.  
  577. LOCAL double frand(ndf, ddf)
  578.      double ndf, ddf;
  579. {
  580.   return((ddf * chisqrand(ndf)) / (ndf * chisqrand(ddf)));
  581. }
  582.  
  583. LOCAL LVAL contrand(which)
  584.      int which;
  585. {
  586.   LVAL next, result;
  587.   int n;
  588.   double dx, da, db;
  589.   
  590.   n = getfixnum(xlgafixnum());
  591.   switch (which) {
  592.   case 'G':
  593.   case 'X': 
  594.   case 'T': da = getposdouble(); break;
  595.   case 'B': 
  596.   case 'F': da = getposdouble(); db = getposdouble(); break;
  597.   }
  598.   xllastarg();
  599.   
  600.   if (n <= 0) return(NIL);
  601.   
  602.   /* protect result pointer */
  603.   xlsave1(result);
  604.   
  605.   result = mklist(n, NIL);
  606.   for (next = result; consp(next); next = cdr(next)) {
  607.     switch (which) {
  608.     case 'U': dx = unirand();         break;
  609.     case 'N': dx = normrand();        break;
  610.     case 'C': dx = cauchyrand();      break;
  611.     case 'G': dx = gammarand(da);     break;
  612.     case 'X': dx = chisqrand(da);     break;
  613.     case 'T': dx = trand(da);         break;
  614.     case 'B': dx = betarand(da, db);  break;
  615.     case 'F': dx = frand(da, db);     break;
  616.     }
  617.     rplaca(next, cvflonum((FLOTYPE) dx));
  618.   }
  619.   
  620.   /* restore the stack frame */
  621.   xlpop();
  622.   
  623.   return(result);
  624. }
  625.  
  626. LOCAL LVAL xsuniformrand() { return(contrand('U')); }
  627. LOCAL LVAL xsnormalrand()  { return(contrand('N')); }
  628. LOCAL LVAL xscauchyrand()  { return(contrand('C')); }
  629. LOCAL LVAL xsgammarand()   { return(contrand('G')); }
  630. LOCAL LVAL xschisqrand()   { return(contrand('X')); }
  631. LOCAL LVAL xstrand()       { return(contrand('T')); }
  632. LOCAL LVAL xsbetarand()    { return(contrand('B')); }
  633. LOCAL LVAL xsfrand()       { return(contrand('F')); }
  634.  
  635. LVAL xsruniformrand() 
  636.     { return(recursive_subr_map_elements(xsuniformrand, xsruniformrand)); }
  637. LVAL xsrnormalrand() 
  638.     { return(recursive_subr_map_elements(xsnormalrand, xsrnormalrand)); }
  639. LVAL xsrcauchyrand() 
  640.     { return(recursive_subr_map_elements(xscauchyrand, xsrcauchyrand)); }
  641. LVAL xsrgammarand() 
  642.     { return(recursive_subr_map_elements(xsgammarand, xsrgammarand)); }
  643. LVAL xsrchisqrand() 
  644.     { return(recursive_subr_map_elements(xschisqrand, xsrchisqrand)); }
  645. LVAL xsrtrand() 
  646.     { return(recursive_subr_map_elements(xstrand, xsrtrand)); }
  647. LVAL xsrbetarand() 
  648.     { return(recursive_subr_map_elements(xsbetarand, xsrbetarand)); }
  649. LVAL xsrfrand() 
  650.     { return(recursive_subr_map_elements(xsfrand, xsrfrand)); }
  651.