modules/Math-Random/helper.c

   1 /* NOTE: RETURN CODES HAVE BEEN CHANGED TO MATCH PERL, I.E.
   2    1 - NOW MEANS OK
   3    0 - NOW MEANS ERROR
   4    */
   5
   6 #include "randlib.h"
   7 #include <stdio.h>
   8 #include <stdlib.h>
   9 #include "helper.h"
  10
  11 static long   *iwork = NULL; /* perl long array, alloc. in 'rspriw'  */
  12 static double *fwork = NULL; /* perl float array, alloc. in 'rsprfw' */
  13 static double *parm  = NULL; /* maintained by 'psetmn' for 'pgenmn'  */
  14
  15 /****************************************************************************
  16                 Perl <-> C (Long) Integer Helper Functions
  17  (these pass single values back and forth, to load/read/manage working array)
  18  ****************************************************************************/
  19
  20 long gvpriw(long index) {
  21   /* Gets the Value at index of the PeRl (long) Integer Working array */
  22   extern long *iwork;
  23
  24   return *(iwork + index);
  25 }
  26
  27 int rspriw(long size) {
  28   /* Request Size for PeRl's (long) int Working array
  29    * returns:
  30    * 1 if successful
  31    * 0 if out of memory
  32    */
  33   extern long *iwork;
  34   static long siwork = 0L;
  35
  36   if (size <= siwork) return 1;
  37   /* else reset array */
  38   if (iwork != NULL) free(iwork);
  39   iwork = (long *) malloc(sizeof(long) * size);
  40   if (iwork != NULL) {
  41     siwork = size;
  42     return 1;
  43   }
  44   fputs(" Unable to allocate randlib (long) int working array:\n",stderr);
  45   fprintf(stderr," Requested number of entries = %ld\n",size);
  46   fputs(" Out of memory in RSPRIW - ABORT\n",stderr);
  47   siwork = 0L;
  48   return 0;
  49 }
  50
  51 /****************************************************************************
  52                 Perl <-> C Float Helper Functions
  53  (these pass single values back and forth, to load/read/manage working array)
  54  ****************************************************************************/
  55
  56 double gvprfw(long index) {
  57   /* Gets the Value at index of the PeRl Float Working array */
  58   extern double *fwork;
  59
  60   return *(fwork + index);
  61 }
  62
  63 void svprfw(long index, double value) {
  64   /* Sets Value in PeRl's Float Working array */
  65   extern double *fwork;
  66
  67   *(fwork + index) = value;
  68 }
  69
  70 int rsprfw(long size) {
  71   /* Request Size for PeRl's Float Working array
  72    * returns:
  73    * 1 if successful
  74    * 0 if out of memory
  75    */
  76   extern double *fwork;
  77   static long sfwork = 0L;
  78
  79   if (size <= sfwork) return 1;
  80   /* else reset array */
  81   if (fwork != NULL) free(fwork);
  82   fwork = (double*) malloc(sizeof(double) * size);
  83   if (fwork != NULL) {
  84     sfwork = size;
  85     return 1;
  86   }
  87   fputs(" Unable to allocate randlib float working array:\n",stderr);
  88   fprintf(stderr," Requested number of entries = %ld\n",size);
  89   fputs(" Out of memory in RSPRFW - ABORT\n",stderr);
  90   sfwork = 0L;
  91   return 0;
  92 }
  93
  94 /*****************************************************************************
  95                            Randlib Helper Functions
  96       These routines call those randlib routines which depend on pointers
  97               (typically those with array input and/or output)
  98  *****************************************************************************/
  99 void pgnprm(long n) {
 100   /* Perl's GeNerate PeRMutation
 101    * Fills perl's (long) integer working array with 0, ... ,n-1
 102    * and randomly permutes it.
 103    * Note: if n <= 0, it does what you'd expect:
 104    * N == 1: array of 0 of length 1
 105    * N <  1: array of length 0
 106    */
 107
 108   /* NOTE: EITHER HERE OR IN PERL IWORK MUST HAVE SIZE CHECKED */
 109
 110   extern long *iwork;
 111   long i;
 112
 113   /* Fills working array ... */
 114   for (i=0L;i<n;i++)
 115     *(iwork + i) = i;
 116
 117   /* ... and randomly permutes it */
 118   genprm(iwork,i);
 119 }
 120
 121 void pgnmul (long n, long ncat) {
 122   /* Perl's GeNerate MULtinomial observation.
 123    * Method: uses void genmul(long n,double *p,long ncat,long *ix) in 'randlib.c'
 124    * Arguments:
 125    * n    - number of events to be classified.
 126    * ncat - number of categories into which the events are classified.
 127    * Notes:
 128    * *p - must be set up first in perl's double working array *fwork.
 129    *      must have at least ncat-1 categories and otherwise make sense.
 130    * *ix - (results) will be perl's (long) integer working array *iwork.
 131    */
 132
 133   /* NOTE: FROM PERL, FWORK MUST HAVE SIZE CHECKED AND BE FILLED */
 134   /* ALSO, HERE OR IN PERL IWORK MUST HAVE SIZE CHECKED */
 135
 136   extern long   *iwork;
 137   extern double *fwork;
 138
 139   /* since all is OK so far, get the obs */
 140   genmul(n, fwork, ncat, iwork);
 141 }
 142
 143 int psetmn(long p) {
 144   /*
 145    * Perl's SET Multivariate Normal
 146    * p - dimension of multivariate normal deviate
 147    *
 148    * Input:
 149    * fwork must be loaded as follows prior to call:
 150    *    Origin = 0 indexing           Origin = 1 indexing
 151    *    (reverse odometer)
 152    *       fwork[0]                 <-> mean[1]
 153    *       fwork[1]                 <-> mean[2]
 154    *        ...                          ...
 155    *       fwork[p - 1]             <-> mean[p]
 156    *       fwork[0 + 0*p + p]       <-> covm[1,1]
 157    *       fwork[1 + 0*p + p]       <-> covm[2,1]
 158    *        ...                          ...
 159    *       fwork[i-1 + (j-1)*p + p] <-> covm[i,j]
 160    *        ...                          ...
 161    *       fwork[p-1 + (p-1)*p + p] <-> covm[p,p]
 162    * Tot:  p*p + p elements                  p*p + p elements
 163    * This should all be done by the Perl calling routine.
 164    *
 165    * Side Effects:
 166    * parm[p*(p+3)/2 + 1] is a file static array which contains all the
 167    * information needed to generate the deviates.
 168    * fwork is essentially destroyed (but not reallocated).
 169    *
 170    * Returns:
 171    * 1 if initialization succeeded
 172    * 0 if out of memory
 173    *
 174    * Method:
 175    * Calls 'setgmn' in "randlib.c":
 176    * void setgmn(double *meanv,double *covm,long p,double *parm)
 177    */
 178
 179   extern double *fwork, *parm;
 180   static long oldp = 0L; /* p from last reallocate of parm */
 181
 182   if (p > oldp) { /* pmn_param is too small; reallocate */
 183     if (parm != NULL) free(parm);
 184     parm = (double *) malloc(sizeof(double)*(p*(p+3L)/2L + 1L));
 185     if (parm == NULL) {
 186       fputs("Out of memory in PSETMN - ABORT",stderr);
 187       fprintf(stderr,
 188               "P = %ld; Requested # of doubles %ld\n",p,p*(p+3L)/2L + 1L);
 189       oldp = 0L;
 190       return 0;
 191     } else {
 192       oldp = p; /* keep track of last reallocation */
 193     }
 194   }
 195   /* initialize parm */
 196   setgmn(fwork, fwork + p, p, parm);
 197   return 1;
 198 }
 199
 200 int pgenmn(void) {
 201   /*
 202    * Perl's GENerate Multivariate Normal
 203    *
 204    * Input: (None)
 205    *
 206    * p - dimension of multivariate normal deviate - gotten from parm[].
 207    * 'psetmn' must be called successfully before this routine is called.
 208    * If that be so, then fwork[] has enough space for the deviate
 209    * and scratch space used by the routine, and parm[] has the
 210    * parameters needed.
 211    *
 212    * Output:
 213    * 0 - generation failed
 214    * 1 - generation succeeded
 215    *
 216    * Side Effects:
 217    * fwork[0] ... fwork[p-1] will contain the deviate.
 218    *
 219    * Method:
 220    * Calls 'genmn' in "randlib.c":
 221    * void genmn(double *parm,double *x,double *work)
 222    */
 223
 224   extern double *fwork, *parm;
 225
 226   /* NOTE: CHECK OF PARM ONLY NEEDED IF PERL SET/GENERATE IS SPLIT */
 227
 228   if (parm != NULL) { /* initialized OK */
 229     long p = (long) *(parm);
 230     genmn(parm,fwork,fwork+p); /* put deviate in fwork */
 231     return 1;
 232
 233   } else { /* not initialized - ABORT */
 234     fputs("PGENMN called before PSETMN called successfully - ABORT\n",
 235           stderr);
 236     fputs("parm not properly initialized in PGENMN - ABORT\n",stderr);
 237     return 0;
 238   }
 239 }
 240
 241 void salfph(char* phrase)
 242 {
 243 /*
 244 **********************************************************************
 245      void salfph(char* phrase)
 246                Set ALl From PHrase
 247
 248                               Function
 249
 250      Uses a phrase (character string) to generate two seeds for the RGN
 251      random number generator, then sets the initial seed of generator 1
 252      to the results.  The initial seeds of the other generators are set
 253      accordingly, and all generators' states are set to these seeds.
 254
 255                               Arguments
 256      phrase --> Phrase to be used for random number generation
 257
 258                               Method
 259      Calls 'setall' (from com.c) with the results of 'phrtsd' (here in
 260      randlib.c).  Please see those functions' comments for details.
 261 **********************************************************************
 262 */
 263 extern void phrtsd(char* phrase,long *seed1,long *seed2);
 264 extern void setall(long iseed1,long iseed2);
 265 static long iseed1, iseed2;
 266
 267 phrtsd(phrase,&iseed1,&iseed2);
 268 setall(iseed1,iseed2);
 269 }