src/libs/mapm/mapm_rnd.c

   1
   2 /*
   3  *  M_APM  -  mapm_rnd.c
   4  *
   5  *  Copyright (C) 1999 - 2007   Michael C. Ring
   6  *
   7  *  Permission to use, copy, and distribute this software and its
   8  *  documentation for any purpose with or without fee is hereby granted,
   9  *  provided that the above copyright notice appear in all copies and
  10  *  that both that copyright notice and this permission notice appear
  11  *  in supporting documentation.
  12  *
  13  *  Permission to modify the software is granted. Permission to distribute
  14  *  the modified code is granted. Modifications are to be distributed by
  15  *  using the file 'license.txt' as a template to modify the file header.
  16  *  'license.txt' is available in the official MAPM distribution.
  17  *
  18  *  This software is provided "as is" without express or implied warranty.
  19  */
  20
  21 /*
  22  *      $Id: mapm_rnd.c,v 1.12 2007/12/03 01:47:17 mike Exp $
  23  *
  24  *      This file contains the Random Number Generator function.
  25  *
  26  *      $Log: mapm_rnd.c,v $
  27  *      Revision 1.12  2007/12/03 01:47:17  mike
  28  *      Update license
  29  *
  30  *      Revision 1.11  2003/10/25 22:55:43  mike
  31  *      add support for National Instruments LabWindows CVI
  32  *
  33  *      Revision 1.10  2002/11/03 22:41:03  mike
  34  *      Updated function parameters to use the modern style
  35  *
  36  *      Revision 1.9  2002/02/14 21:50:45  mike
  37  *      add _set_random_seed
  38  *
  39  *      Revision 1.8  2001/07/16 19:30:32  mike
  40  *      add function M_free_all_rnd
  41  *
  42  *      Revision 1.7  2001/03/20 17:19:45  mike
  43  *      use a new multiplier
  44  *
  45  *      Revision 1.6  2000/08/20 23:46:07  mike
  46  *      add more possible multupliers (no code changes)
  47  *
  48  *      Revision 1.5  1999/09/19 23:32:14  mike
  49  *      added comments
  50  *
  51  *      Revision 1.4  1999/09/18 03:49:25  mike
  52  *      *** empty log message ***
  53  *
  54  *      Revision 1.3  1999/09/18 03:35:36  mike
  55  *      only prototype get_microsec for non-DOS
  56  *
  57  *      Revision 1.2  1999/09/18 02:35:36  mike
  58  *      delete debug printf's
  59  *
  60  *      Revision 1.1  1999/09/18 02:26:52  mike
  61  *      Initial revision
  62  */
  63
  64 #include "m_apm_lc.h"
  65
  66 #ifndef _HAVE_NI_LABWIN_CVI_
  67 #ifdef MSDOS
  68 #include <time.h>
  69 #include <sys/timeb.h>
  70 #else
  71 #include <sys/time.h>
  72 extern  void    M_get_microsec(unsigned long *, long *);
  73 #endif
  74 #endif
  75
  76 #ifdef _HAVE_NI_LABWIN_CVI_
  77 #include <time.h>
  78 #include <utility.h>
  79 #include <ansi/math.h>
  80 #endif
  81
  82 extern  void    M_reverse_string(char *);
  83 extern  void    M_get_rnd_seed(M_APM);
  84
  85 static  M_APM   M_rnd_aa;
  86 static  M_APM   M_rnd_mm;
  87 static  M_APM   M_rnd_XX;
  88 static  M_APM   M_rtmp0;
  89 static  M_APM   M_rtmp1;
  90
  91 static  int     M_firsttime2 = TRUE;
  92
  93 /*
  94         Used Knuth's The Art of Computer Programming, Volume 2 as
  95         the basis. Assuming the random number is X, compute
  96         (where all the math is performed on integers) :
  97
  98         X = (a * X + c) MOD m
  99
 100         From Knuth:
 101
 102         'm' should be large, at least 2^30 : we use 1.0E+15
 103
 104         'a' should be between .01m and .99m and not have a simple
 105         pattern. 'a' should not have any large factors in common
 106         with 'm' and (since 'm' is a power of 10) if 'a' MOD 200
 107         = 21 then all 'm' different possible values will be
 108         generated before 'X' starts to repeat.
 109
 110         We use 'a' = 716805947629621.
 111
 112         This is a prime number and also meets 'a' MOD 200 = 21.
 113         Commented out below are many potential multipliers that
 114         are all prime and meet 'a' MOD 200 = 21.
 115
 116         There are few restrictions on 'c' except 'c' can have no
 117         factor in common with 'm', hence we set 'c' = 'a'.
 118
 119         On the first call, the system time is used to initialize X.
 120 */
 121
 122 /*
 123  *  the following constants are all potential multipliers. they are
 124  *  all prime numbers that also meet the criteria of NUM mod 200 = 21.
 125  */
 126
 127 /*
 128 439682071525421   439682071528421   439682071529221   439682071529821
 129 439682071530421   439682071532021   439682071538821   439682071539421
 130 439682071540021   439682071547021   439682071551221   439682071553821
 131 439682071555421   439682071557221   439682071558021   439682071558621
 132 439682071559821   439652381461621   439652381465221   439652381465621
 133 439652381466421   439652381467421   439652381468621   439652381470021
 134 439652381471221   439652381477021   439652381484221   439652381488421
 135 439652381491021   439652381492021   439652381494021   439652381496821
 136 617294387035621   617294387038621   617294387039221   617294387044421
 137 617294387045221   617294387048621   617294387051621   617294387051821
 138 617294387053621   617294387058421   617294387064221   617294387065621
 139 617294387068621   617294387069221   617294387069821   617294387070421
 140 617294387072021   617294387072621   617294387073821   617294387076821
 141 649378126517621   649378126517821   649378126518221   649378126520821
 142 649378126523821   649378126525621   649378126526621   649378126528421
 143 649378126529621   649378126530821   649378126532221   649378126533221
 144 649378126535221   649378126539421   649378126543621   649378126546021
 145 649378126546421   649378126549421   649378126550821   649378126555021
 146 649378126557421   649378126560221   649378126561621   649378126562021
 147 649378126564621   649378126565821   672091582360421   672091582364221
 148 672091582364621   672091582367021   672091582368421   672091582369021
 149 672091582370821   672091582371421   672091582376821   672091582380821
 150 716805243983221   716805243984821   716805947623621   716805947624621
 151 716805947629021   716805947629621   716805947630621   716805947633621
 152 716805947634221   716805947635021   716805947635621   716805947642221
 153 */
 154
 155 /****************************************************************************/
 156 void    M_free_all_rnd()
 157 {
 158 if (M_firsttime2 == FALSE)
 159   {
 160    m_apm_free(M_rnd_aa);
 161    m_apm_free(M_rnd_mm);
 162    m_apm_free(M_rnd_XX);
 163    m_apm_free(M_rtmp0);
 164    m_apm_free(M_rtmp1);
 165
 166    M_firsttime2 = TRUE;
 167   }
 168 }
 169 /****************************************************************************/
 170 void    m_apm_set_random_seed(char *ss)
 171 {
 172 M_APM   btmp;
 173
 174 if (M_firsttime2)
 175   {
 176    btmp = M_get_stack_var();
 177    m_apm_get_random(btmp);
 178    M_restore_stack(1);
 179   }
 180
 181 m_apm_set_string(M_rnd_XX, ss);
 182 }
 183 /****************************************************************************/
 184 /*
 185  *  compute X = (a * X + c) MOD m       where c = a
 186  */
 187 void    m_apm_get_random(M_APM mrnd)
 188 {
 189
 190 if (M_firsttime2)         /* use the system time as the initial seed value */
 191   {
 192    M_firsttime2 = FALSE;
 193
 194    M_rnd_aa = m_apm_init();
 195    M_rnd_XX = m_apm_init();
 196    M_rnd_mm = m_apm_init();
 197    M_rtmp0  = m_apm_init();
 198    M_rtmp1  = m_apm_init();
 199
 200    /* set the multiplier M_rnd_aa and M_rnd_mm */
 201
 202    m_apm_set_string(M_rnd_aa, "716805947629621");
 203    m_apm_set_string(M_rnd_mm, "1.0E15");
 204
 205    M_get_rnd_seed(M_rnd_XX);
 206   }
 207
 208 m_apm_multiply(M_rtmp0, M_rnd_XX, M_rnd_aa);
 209 m_apm_add(M_rtmp1, M_rtmp0, M_rnd_aa);
 210 m_apm_integer_div_rem(M_rtmp0, M_rnd_XX, M_rtmp1, M_rnd_mm);
 211 m_apm_copy(mrnd, M_rnd_XX);
 212 mrnd->m_apm_exponent -= 15;
 213 }
 214 /****************************************************************************/
 215 void    M_reverse_string(char *s)
 216 {
 217 int     ct;
 218 char    ch, *p1, *p2;
 219
 220 if ((ct = strlen(s)) <= 1)
 221   return;
 222
 223 p1 = s;
 224 p2 = s + ct - 1;
 225 ct /= 2;
 226
 227 while (TRUE)
 228   {
 229    ch    = *p1;
 230    *p1++ = *p2;
 231    *p2-- = ch;
 232
 233    if (--ct == 0)
 234      break;
 235   }
 236 }
 237 /****************************************************************************/
 238
 239 #ifndef _HAVE_NI_LABWIN_CVI_
 240
 241 #ifdef MSDOS
 242
 243 /****************************************************************************/
 244 /*
 245  *  for DOS / Win 9x/NT systems : use 'ftime'
 246  */
 247 void    M_get_rnd_seed(M_APM mm)
 248 {
 249 int              millisec;
 250 time_t           timestamp;
 251 unsigned long    ul;
 252 char             ss[32], buf1[48], buf2[32];
 253 struct timeb     timebuffer;
 254 M_APM            atmp;
 255
 256 atmp = M_get_stack_var();
 257
 258 ftime(&timebuffer);
 259
 260 millisec  = (int)timebuffer.millitm;
 261 timestamp = timebuffer.time;
 262 ul        = (unsigned long)(timestamp / 7);
 263 ul       += timestamp + 537;
 264 strcpy(ss,ctime(&timestamp));        /* convert to string and copy to ss */
 265
 266 sprintf(buf1,"%d",(millisec / 10));
 267 sprintf(buf2,"%lu",ul);
 268
 269 ss[0] = ss[18];
 270 ss[1] = ss[17];
 271 ss[2] = ss[15];
 272 ss[3] = ss[14];
 273 ss[4] = ss[12];
 274 ss[5] = ss[11];
 275 ss[6] = ss[9];
 276 ss[7] = ss[23];
 277 ss[8] = ss[20];
 278 ss[9] = '\0';
 279
 280 M_reverse_string(buf2);
 281 strcat(buf1,buf2);
 282 strcat(buf1,ss);
 283
 284 m_apm_set_string(atmp, buf1);
 285 atmp->m_apm_exponent = 15;
 286 m_apm_integer_divide(mm, atmp, MM_One);
 287
 288 M_restore_stack(1);
 289 }
 290 /****************************************************************************/
 291
 292 #else
 293
 294 /****************************************************************************/
 295 /*
 296  *  for unix systems : use 'gettimeofday'
 297  */
 298 void    M_get_rnd_seed(M_APM mm)
 299 {
 300 unsigned long    sec3;
 301 long             usec3;
 302 char             buf1[32], buf2[32];
 303 M_APM            atmp;
 304
 305 atmp = M_get_stack_var();
 306 M_get_microsec(&sec3,&usec3);
 307
 308 sprintf(buf1,"%ld",usec3);
 309 sprintf(buf2,"%lu",sec3);
 310 M_reverse_string(buf2);
 311 strcat(buf1,buf2);
 312
 313 m_apm_set_string(atmp, buf1);
 314 atmp->m_apm_exponent = 15;
 315 m_apm_integer_divide(mm, atmp, MM_One);
 316
 317 M_restore_stack(1);
 318 }
 319 /****************************************************************************/
 320 void     M_get_microsec(unsigned long *sec, long *usec)
 321 {
 322 struct timeval time_now;           /* current time for elapsed time check */
 323 struct timezone time_zone;         /* time zone for gettimeofday call     */
 324
 325 gettimeofday(&time_now, &time_zone);                  /* get current time */
 326
 327 *sec  = time_now.tv_sec;
 328 *usec = time_now.tv_usec;
 329 }
 330 /****************************************************************************/
 331
 332 #endif
 333 #endif
 334
 335 #ifdef _HAVE_NI_LABWIN_CVI_
 336
 337 /****************************************************************************/
 338 /*
 339  *  for National Instruments LabWindows CVI
 340  */
 341
 342 void    M_get_rnd_seed(M_APM mm)
 343 {
 344 double           timer0;
 345 int              millisec;
 346 char             *cvi_time, *cvi_date, buf1[64], buf2[32];
 347 M_APM            atmp;
 348
 349 atmp = M_get_stack_var();
 350
 351 cvi_date = DateStr();
 352 cvi_time = TimeStr();
 353 timer0   = Timer();
 354
 355 /*
 356  *  note that Timer() is not syncronized to TimeStr(),
 357  *  but we don't care here since we are just looking
 358  *  for a random source of digits.
 359  */
 360
 361 millisec = (int)(0.01 + 1000.0 * (timer0 - floor(timer0)));
 362
 363 sprintf(buf1, "%d", millisec);
 364
 365 buf2[0]  = cvi_time[6]; /* time format: "HH:MM:SS" */
 366 buf2[1]  = cvi_time[7];
 367 buf2[2]  = cvi_time[3];
 368 buf2[3]  = cvi_time[4];
 369 buf2[4]  = cvi_time[0];
 370 buf2[5]  = cvi_time[1];
 371
 372 buf2[6]  = cvi_date[3]; /* date format: "MM-DD-YYYY" */
 373 buf2[7]  = cvi_date[4];
 374 buf2[8]  = cvi_date[0];
 375 buf2[9]  = cvi_date[1];
 376 buf2[10] = cvi_date[8];
 377 buf2[11] = cvi_date[9];
 378 buf2[12] = cvi_date[7];
 379
 380 buf2[13] = '4';
 381 buf2[14] = '7';
 382 buf2[15] = '\0';
 383
 384 strcat(buf1, buf2);
 385
 386 m_apm_set_string(atmp, buf1);
 387 atmp->m_apm_exponent = 15;
 388 m_apm_integer_divide(mm, atmp, MM_One);
 389
 390 M_restore_stack(1);
 391 }
 392
 393 #endif
 394
 395 /****************************************************************************/
 396