[BACK]Return to ranlib.c CVS log [TXT] [DIR] Up to [LON-CAPA] / capa / capa51 / pProj

Diff for /capa/capa51/pProj/ranlib.c between versions 1.1 and 1.6

version 1.1, 1999/09/28 21:26:21 version 1.6, 2001/10/29 19:36:53
Line 1 Line 1
   /* another library of funtions to support the rand number generator
      Copyright (C) 1992-2000 Michigan State University
   
      The CAPA system is free software; you can redistribute it and/or
      modify it under the terms of the GNU General Public License as
      published by the Free Software Foundation; either version 2 of the
      License, or (at your option) any later version.
   
      The CAPA system is distributed in the hope that it will be useful,
      but WITHOUT ANY WARRANTY; without even the implied warranty of
      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      General Public License for more details.
   
      You should have received a copy of the GNU General Public
      License along with the CAPA system; see the file COPYING.  If not,
      write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
      Boston, MA 02111-1307, USA.
   
      As a special exception, you have permission to link this program
      with the TtH/TtM library and distribute executables, as long as you
      follow the requirements of the GNU GPL in regard to all of the
      software in the executable aside from TtH/TtM.
   */
   
 #include "ranlib.h"  #include "ranlib.h"
 #include <stdio.h>  #include <stdio.h>
Line 39  static long qsame; Line 62  static long qsame;
     qsame = olda == aa && oldb == bb;      qsame = olda == aa && oldb == bb;
     if(qsame) goto S20;      if(qsame) goto S20;
     if(!(aa <= 0.0 || bb <= 0.0)) goto S10;      if(!(aa <= 0.0 || bb <= 0.0)) goto S10;
     puts(" AA or BB <= 0 in GENBET - Abort!");      fprintf(stderr," AA or BB <= 0 in GENBET - Abort!\n");
     printf(" AA: %16.6E BB %16.6E\n",aa,bb);      fprintf(stderr," AA: %16.6E BB %16.6E\n",aa,bb);
     exit(1);      exit(1);
 S10:  S10:
     olda = aa;      olda = aa;
Line 48  S10: Line 71  S10:
 S20:  S20:
     if(!(min(aa,bb) > 1.0)) goto S100;      if(!(min(aa,bb) > 1.0)) goto S100;
 /*  /*
      Alborithm BB       Algorithm BB
      Initialize       Initialize
 */  */
     if(qsame) goto S30;      if(qsame) goto S30;
Line 190  float genchi(float df) Line 213  float genchi(float df)
 static float genchi;  static float genchi;
   
     if(!(df <= 0.0)) goto S10;      if(!(df <= 0.0)) goto S10;
     puts("DF <= 0 in GENCHI - ABORT");      fprintf(stderr,"DF <= 0 in GENCHI - ABORT\n");
     printf("Value of DF: %16.6E\n",df);      fprintf(stderr,"Value of DF: %16.6E\n",df);
     exit(1);      exit(1);
 S10:  S10:
     genchi = 2.0*gengam(1.0,df/2.0);      genchi = 2.0*gengam(1.0,df/2.0);
Line 246  float genf(float dfn,float dfd) Line 269  float genf(float dfn,float dfd)
 static float genf,xden,xnum;  static float genf,xden,xnum;
   
     if(!(dfn <= 0.0 || dfd <= 0.0)) goto S10;      if(!(dfn <= 0.0 || dfd <= 0.0)) goto S10;
     puts("Degrees of freedom nonpositive in GENF - abort!");      fprintf(stderr,"Degrees of freedom nonpositive in GENF - abort!\n");
     printf("DFN value: %16.6EDFD value: %16.6E\n",dfn,dfd);      fprintf(stderr,"DFN value: %16.6EDFD value: %16.6E\n",dfn,dfd);
     exit(1);      exit(1);
 S10:  S10:
     xnum = genchi(dfn)/dfn;      xnum = genchi(dfn)/dfn;
Line 256  S10: Line 279  S10:
 */  */
     xden = genchi(dfd)/dfd;      xden = genchi(dfd)/dfd;
     if(!(xden <= 9.999999999998E-39*xnum)) goto S20;      if(!(xden <= 9.999999999998E-39*xnum)) goto S20;
     puts(" GENF - generated numbers would cause overflow");      fprintf(stderr," GENF - generated numbers would cause overflow\n");
     printf(" Numerator %16.6E Denominator %16.6E\n",xnum,xden);      fprintf(stderr," Numerator %16.6E Denominator %16.6E\n",xnum,xden);
     puts(" GENF returning 1.0E38");      fprintf(stderr," GENF returning 1.0E38\n");
     genf = 1.0E38;      genf = 1.0E38;
     goto S30;      goto S30;
 S20:  S20:
Line 378  float gennch(float df,float xnonc) Line 401  float gennch(float df,float xnonc)
 static float gennch;  static float gennch;
   
     if(!(df <= 1.0 || xnonc < 0.0)) goto S10;      if(!(df <= 1.0 || xnonc < 0.0)) goto S10;
     puts("DF <= 1 or XNONC < 0 in GENNCH - ABORT");      fprintf(stderr,"DF <= 1 or XNONC < 0 in GENNCH - ABORT\n");
     printf("Value of DF: %16.6E Value of XNONC%16.6E\n",df,xnonc);      fprintf(stderr,"Value of DF: %16.6E Value of XNONC%16.6E\n",df,xnonc);
     exit(1);      exit(1);
 S10:  S10:
     gennch = genchi(df-1.0)+pow(gennor(sqrt(xnonc),1.0),2.0);      gennch = genchi(df-1.0)+pow(gennor(sqrt(xnonc),1.0),2.0);
Line 413  static long qcond; Line 436  static long qcond;
   
     qcond = dfn <= 1.0 || dfd <= 0.0 || xnonc < 0.0;      qcond = dfn <= 1.0 || dfd <= 0.0 || xnonc < 0.0;
     if(!qcond) goto S10;      if(!qcond) goto S10;
     puts("In GENNF - Either (1) Numerator DF <= 1.0 or");      fprintf(stderr,"In GENNF - Either (1) Numerator DF <= 1.0 or\n");
     puts("(2) Denominator DF < 0.0 or ");      fprintf(stderr,"(2) Denominator DF < 0.0 or \n");
     puts("(3) Noncentrality parameter < 0.0");      fprintf(stderr,"(3) Noncentrality parameter < 0.0\n");
     printf("DFN value: %16.6EDFD value: %16.6EXNONC value: \n%16.6E\n",dfn,dfd,      fprintf(stderr,"DFN value: %16.6EDFD value: %16.6EXNONC value: \n%16.6E\n",dfn,dfd,
       xnonc);        xnonc);
     exit(1);      exit(1);
 S10:  S10:
Line 426  S10: Line 449  S10:
 */  */
     xden = genchi(dfd)/dfd;      xden = genchi(dfd)/dfd;
     if(!(xden <= 9.999999999998E-39*xnum)) goto S20;      if(!(xden <= 9.999999999998E-39*xnum)) goto S20;
     puts(" GENNF - generated numbers would cause overflow");      fprintf(stderr," GENNF - generated numbers would cause overflow\n");
     printf(" Numerator %16.6E Denominator %16.6E\n",xnum,xden);      fprintf(stderr," Numerator %16.6E Denominator %16.6E\n",xnum,xden);
     puts(" GENNF returning 1.0E38");      fprintf(stderr," GENNF returning 1.0E38\n");
     gennf = 1.0E38;      gennf = 1.0E38;
     goto S30;      goto S30;
 S20:  S20:
Line 436  S20: Line 459  S20:
 S30:  S30:
     return gennf;      return gennf;
 }  }
   
 float gennor(float av,float sd)  float gennor(float av,float sd)
 /*  /*
 **********************************************************************  **********************************************************************
Line 458  float gennor(float av,float sd) Line 482  float gennor(float av,float sd)
 **********************************************************************  **********************************************************************
 */  */
 {  {
 static float gennor;  float  gennor;
   float  tmp_f;
   
       tmp_f = snorm();
       
       gennor = sd*tmp_f+av;
       return (gennor);
   }
   
   float capa_gennor(double *num_d, float av,float sd)
   /*
   **********************************************************************
        float gennor(float av,float sd)
            GENerate random deviate from a NORmal distribution
                                 Function
        Generates a single random deviate from a normal distribution
        with mean, AV, and standard deviation, SD.
                                 Arguments
        av --> Mean of the normal distribution.
        sd --> Standard deviation of the normal distribution.
                                 Method
        Renames SNORM from TOMS as slightly modified by BWB to use RANF
        instead of SUNIF.
        For details see:
                  Ahrens, J.H. and Dieter, U.
                  Extensions of Forsythe's Method for Random
                  Sampling from the Normal Distribution.
                  Math. Comput., 27,124 (Oct. 1973), 927 - 937.
   **********************************************************************
   */
   {
   float  gen_num;
   float  tmp_f;
   
     gennor = sd*snorm()+av;      tmp_f = snorm();
     return gennor;      
       gen_num = sd*tmp_f+av;
       /* fprintf(stderr,"SNORM()=%f,GENNOR()=%f,%f*%f+%f\n",tmp_f,gen_num,sd,tmp_f,av); */
       *num_d = (double)gen_num;
       
       gen_num = (float)37.358341;
       return (gen_num);
 }  }
   
   
 void genprm(long *iarray,int larray)  void genprm(long *iarray,int larray)
 /*  /*
 **********************************************************************  **********************************************************************
Line 500  float genunf(float low,float high) Line 564  float genunf(float low,float high)
 static float genunf;  static float genunf;
   
     if(!(low > high)) goto S10;      if(!(low > high)) goto S10;
     printf("LOW > HIGH in GENUNF: LOW %16.6E HIGH: %16.6E\n",low,high);      fprintf(stderr,"LOW > HIGH in GENUNF: LOW %16.6E HIGH: %16.6E\n",low,high);
     puts("Abort");      fprintf(stderr,"Abort\n");
     exit(1);      exit(1);
 S10:  S10:
     genunf = low+(high-low)*ranf();      genunf = low+(high-low)*ranf();
Line 525  static long curntg = 1; Line 589  static long curntg = 1;
     if(getset == 0) *g = curntg;      if(getset == 0) *g = curntg;
     else  {      else  {
         if(*g < 0 || *g > numg) {          if(*g < 0 || *g > numg) {
             puts(" Generator number out of range in GSCGN");              fprintf(stderr," Generator number out of range in GSCGN\n");
             exit(0);              exit(0);
         }          }
         curntg = *g;          curntg = *g;
Line 1083  long ignuin(long low,long high) Line 1147  long ignuin(long low,long high)
 static long ignuin,ign,maxnow,range,ranp1;  static long ignuin,ign,maxnow,range,ranp1;
   
     if(!(low > high)) goto S10;      if(!(low > high)) goto S10;
     puts(" low > high in ignuin - ABORT");      fprintf(stderr," low > high in ignuin - ABORT\n");
     exit(1);      exit(1);
   
 S10:  S10:
     range = high-low;      range = high-low;
     if(!(range > maxnum)) goto S20;      if(!(range > maxnum)) goto S20;
     puts(" high - low too large in ignuin - ABORT");      fprintf(stderr," high - low too large in ignuin - ABORT\n");
     exit(1);      exit(1);
   
 S20:  S20:
Line 1150  static long mltmod,a0,a1,k,p,q,qh,rh; Line 1214  static long mltmod,a0,a1,k,p,q,qh,rh;
       machine. On a different machine recompute H        machine. On a different machine recompute H
 */  */
     if(!(a <= 0 || a >= m || s <= 0 || s >= m)) goto S10;      if(!(a <= 0 || a >= m || s <= 0 || s >= m)) goto S10;
     puts(" a, m, s out of order in mltmod - ABORT!");      fprintf(stderr," a, m, s out of order in mltmod - ABORT!\n");
     printf(" a = %12ld s = %12ld m = %12ld\n",a,s,m);      fprintf(stderr," a = %12ld s = %12ld m = %12ld\n",a,s,m);
     puts(" mltmod requires: 0 < a < m; 0 < s < m");      fprintf(stderr," mltmod requires: 0 < a < m; 0 < s < m\n");
     exit(1);      exit(1);
 S10:  S10:
     if(!(a < h)) goto S20;      if(!(a < h)) goto S20;
Line 1300  float ranf(void) Line 1364  float ranf(void)
 */  */
 {  {
 static float ranf;  static float ranf;
   long    tmp_l;
   double  tmp_d;
 /*  /*
      4.656613057E-10 is 1/M1  M1 is set in a data statement in IGNLGI       4.656613057E-10 is 1/M1  M1 is set in a data statement in IGNLGI
       and is currently 2147483563. If M1 changes, change this also.        and is currently 2147483563. If M1 changes, change this also.
 */  */
     ranf = ignlgi()*4.656613057E-10;      tmp_l = ignlgi();
       tmp_d = (double)tmp_l * (double)4.656613057E-10;
       ranf = (float)tmp_d;
       /* printf("RANF()=%f\n",ranf); */
     return ranf;      return ranf;
 }  }
   float capa_ranf(void)
   /*
   **********************************************************************
        float ranf(void)
                   RANDom number generator as a Function
        Returns a random floating point number from a uniform distribution
        over 0 - 1 (endpoints of this interval are not returned) using the
        current generator
        This is a transcription from Pascal to Fortran of routine
        Uniform_01 from the paper
        L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
        with Splitting Facilities." ACM Transactions on Mathematical
        Software, 17:98-111 (1991)
   **********************************************************************
   */
   {
     float ran_f;
     long  my_ran;
     double my_doub;
   /*
        4.656613057E-10 is 1/M1  M1 is set in a data statement in IGNLGI
         and is currently 2147483563. If M1 changes, change this also.
   */
       my_ran = ignlgi();
       /* fprintf(stderr,"MY_ignlgi=%ld -- first time\n",my_ran); */
       /* ran_f = my_ran * 4.656613057E-10; */
       
       my_doub = (double)my_ran * (double)4.656613057E-10;
       fprintf(stderr,"MY_ranf in double=%.15g -- first time\n",my_doub);
       ran_f = (float)my_doub;
       return (ran_f);
   }
   
 void setgmn(float *meanv,float *covm,long p,float *parm)  void setgmn(float *meanv,float *covm,long p,float *parm)
 /*  /*
 **********************************************************************  **********************************************************************
Line 1339  static long i,icount,info,j,D2,D3,D4,D5; Line 1441  static long i,icount,info,j,D2,D3,D4,D5;
      TEST THE INPUT       TEST THE INPUT
 */  */
     if(!(p <= 0)) goto S10;      if(!(p <= 0)) goto S10;
     puts("P nonpositive in SETGMN");      fprintf(stderr,"P nonpositive in SETGMN\n");
     printf("Value of P: %12ld\n",p);      fprintf(stderr,"Value of P: %12ld\n",p);
     exit(1);      exit(1);
 S10:  S10:
     *parm = p;      *parm = p;
Line 1353  S10: Line 1455  S10:
 */  */
     spofa(covm,p,p,&info);      spofa(covm,p,p,&info);
     if(!(info != 0)) goto S30;      if(!(info != 0)) goto S30;
     puts(" COVM not positive definite in SETGMN");      fprintf(stderr," COVM not positive definite in SETGMN\n");
     exit(1);      exit(1);
 S30:  S30:
     icount = p+1;      icount = p+1;
Removed from v.1.1  
changed lines
  Added in v.1.6

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>