--- capa/capa51/pProj/ranlib.c 1999/09/28 21:26:21 1.1.1.1 +++ capa/capa51/pProj/ranlib.c 2000/07/07 18:33:03 1.4 @@ -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 Library 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 + Library General Public License for more details. + + You should have received a copy of the GNU Library 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 @@ -48,7 +71,7 @@ S10: S20: if(!(min(aa,bb) > 1.0)) goto S100; /* - Alborithm BB + Algorithm BB Initialize */ if(qsame) goto S30; @@ -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; + /* printf("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) /* ********************************************************************** @@ -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(); + /* printf("MY_ignlgi=%ld -- first time\n",my_ran); */ + /* ran_f = my_ran * 4.656613057E-10; */ + + my_doub = (double)my_ran * (double)4.656613057E-10; + printf("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) /* **********************************************************************