--- capa/capa51/pProj/ranlib.c 1999/09/28 21:26:21 1.1 +++ capa/capa51/pProj/ranlib.c 2001/10/29 19:36:53 1.6 @@ -1,3 +1,26 @@ +/* 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 @@ -39,8 +62,8 @@ static long qsame; qsame = olda == aa && oldb == bb; if(qsame) goto S20; if(!(aa <= 0.0 || bb <= 0.0)) goto S10; - puts(" AA or BB <= 0 in GENBET - Abort!"); - printf(" AA: %16.6E BB %16.6E\n",aa,bb); + fprintf(stderr," AA or BB <= 0 in GENBET - Abort!\n"); + fprintf(stderr," AA: %16.6E BB %16.6E\n",aa,bb); exit(1); S10: olda = aa; @@ -48,7 +71,7 @@ S10: S20: if(!(min(aa,bb) > 1.0)) goto S100; /* - Alborithm BB + Algorithm BB Initialize */ if(qsame) goto S30; @@ -190,8 +213,8 @@ float genchi(float df) static float genchi; if(!(df <= 0.0)) goto S10; - puts("DF <= 0 in GENCHI - ABORT"); - printf("Value of DF: %16.6E\n",df); + fprintf(stderr,"DF <= 0 in GENCHI - ABORT\n"); + fprintf(stderr,"Value of DF: %16.6E\n",df); exit(1); S10: genchi = 2.0*gengam(1.0,df/2.0); @@ -246,8 +269,8 @@ float genf(float dfn,float dfd) static float genf,xden,xnum; if(!(dfn <= 0.0 || dfd <= 0.0)) goto S10; - puts("Degrees of freedom nonpositive in GENF - abort!"); - printf("DFN value: %16.6EDFD value: %16.6E\n",dfn,dfd); + fprintf(stderr,"Degrees of freedom nonpositive in GENF - abort!\n"); + fprintf(stderr,"DFN value: %16.6EDFD value: %16.6E\n",dfn,dfd); exit(1); S10: xnum = genchi(dfn)/dfn; @@ -256,9 +279,9 @@ S10: */ xden = genchi(dfd)/dfd; if(!(xden <= 9.999999999998E-39*xnum)) goto S20; - puts(" GENF - generated numbers would cause overflow"); - printf(" Numerator %16.6E Denominator %16.6E\n",xnum,xden); - puts(" GENF returning 1.0E38"); + fprintf(stderr," GENF - generated numbers would cause overflow\n"); + fprintf(stderr," Numerator %16.6E Denominator %16.6E\n",xnum,xden); + fprintf(stderr," GENF returning 1.0E38\n"); genf = 1.0E38; goto S30; S20: @@ -378,8 +401,8 @@ float gennch(float df,float xnonc) static float gennch; if(!(df <= 1.0 || xnonc < 0.0)) goto S10; - puts("DF <= 1 or XNONC < 0 in GENNCH - ABORT"); - printf("Value of DF: %16.6E Value of XNONC%16.6E\n",df,xnonc); + fprintf(stderr,"DF <= 1 or XNONC < 0 in GENNCH - ABORT\n"); + fprintf(stderr,"Value of DF: %16.6E Value of XNONC%16.6E\n",df,xnonc); exit(1); S10: gennch = genchi(df-1.0)+pow(gennor(sqrt(xnonc),1.0),2.0); @@ -413,10 +436,10 @@ static long qcond; qcond = dfn <= 1.0 || dfd <= 0.0 || xnonc < 0.0; if(!qcond) goto S10; - puts("In GENNF - Either (1) Numerator DF <= 1.0 or"); - puts("(2) Denominator DF < 0.0 or "); - puts("(3) Noncentrality parameter < 0.0"); - printf("DFN value: %16.6EDFD value: %16.6EXNONC value: \n%16.6E\n",dfn,dfd, + fprintf(stderr,"In GENNF - Either (1) Numerator DF <= 1.0 or\n"); + fprintf(stderr,"(2) Denominator DF < 0.0 or \n"); + fprintf(stderr,"(3) Noncentrality parameter < 0.0\n"); + fprintf(stderr,"DFN value: %16.6EDFD value: %16.6EXNONC value: \n%16.6E\n",dfn,dfd, xnonc); exit(1); S10: @@ -426,9 +449,9 @@ S10: */ xden = genchi(dfd)/dfd; if(!(xden <= 9.999999999998E-39*xnum)) goto S20; - puts(" GENNF - generated numbers would cause overflow"); - printf(" Numerator %16.6E Denominator %16.6E\n",xnum,xden); - puts(" GENNF returning 1.0E38"); + fprintf(stderr," GENNF - generated numbers would cause overflow\n"); + fprintf(stderr," Numerator %16.6E Denominator %16.6E\n",xnum,xden); + fprintf(stderr," GENNF returning 1.0E38\n"); gennf = 1.0E38; goto S30; S20: @@ -436,6 +459,7 @@ S20: S30: return gennf; } + float gennor(float av,float sd) /* ********************************************************************** @@ -458,11 +482,51 @@ 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; - return gennor; + tmp_f = snorm(); + + 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) /* ********************************************************************** @@ -500,8 +564,8 @@ float genunf(float low,float high) static float genunf; if(!(low > high)) goto S10; - printf("LOW > HIGH in GENUNF: LOW %16.6E HIGH: %16.6E\n",low,high); - puts("Abort"); + fprintf(stderr,"LOW > HIGH in GENUNF: LOW %16.6E HIGH: %16.6E\n",low,high); + fprintf(stderr,"Abort\n"); exit(1); S10: genunf = low+(high-low)*ranf(); @@ -525,7 +589,7 @@ static long curntg = 1; if(getset == 0) *g = curntg; else { 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); } curntg = *g; @@ -1083,13 +1147,13 @@ long ignuin(long low,long high) static long ignuin,ign,maxnow,range,ranp1; if(!(low > high)) goto S10; - puts(" low > high in ignuin - ABORT"); + fprintf(stderr," low > high in ignuin - ABORT\n"); exit(1); S10: range = high-low; 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); S20: @@ -1150,9 +1214,9 @@ static long mltmod,a0,a1,k,p,q,qh,rh; machine. On a different machine recompute H */ if(!(a <= 0 || a >= m || s <= 0 || s >= m)) goto S10; - puts(" a, m, s out of order in mltmod - ABORT!"); - printf(" a = %12ld s = %12ld m = %12ld\n",a,s,m); - puts(" mltmod requires: 0 < a < m; 0 < s < m"); + fprintf(stderr," a, m, s out of order in mltmod - ABORT!\n"); + fprintf(stderr," a = %12ld s = %12ld m = %12ld\n",a,s,m); + fprintf(stderr," mltmod requires: 0 < a < m; 0 < s < m\n"); exit(1); S10: if(!(a < h)) goto S20; @@ -1300,13 +1364,51 @@ float ranf(void) */ { static float ranf; +long tmp_l; +double tmp_d; /* 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. */ - 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; } +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) /* ********************************************************************** @@ -1339,8 +1441,8 @@ static long i,icount,info,j,D2,D3,D4,D5; TEST THE INPUT */ if(!(p <= 0)) goto S10; - puts("P nonpositive in SETGMN"); - printf("Value of P: %12ld\n",p); + fprintf(stderr,"P nonpositive in SETGMN\n"); + fprintf(stderr,"Value of P: %12ld\n",p); exit(1); S10: *parm = p; @@ -1353,7 +1455,7 @@ S10: */ spofa(covm,p,p,&info); if(!(info != 0)) goto S30; - puts(" COVM not positive definite in SETGMN"); + fprintf(stderr," COVM not positive definite in SETGMN\n"); exit(1); S30: icount = p+1;