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>
72extern  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
82extern  void	M_reverse_string(char *);
83extern  void    M_get_rnd_seed(M_APM);
84
85static	M_APM   M_rnd_aa;
86static  M_APM   M_rnd_mm;
87static  M_APM   M_rnd_XX;
88static  M_APM   M_rtmp0;
89static  M_APM   M_rtmp1;
90
91static  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/*
128439682071525421   439682071528421   439682071529221   439682071529821
129439682071530421   439682071532021   439682071538821   439682071539421
130439682071540021   439682071547021   439682071551221   439682071553821
131439682071555421   439682071557221   439682071558021   439682071558621
132439682071559821   439652381461621   439652381465221   439652381465621
133439652381466421   439652381467421   439652381468621   439652381470021
134439652381471221   439652381477021   439652381484221   439652381488421
135439652381491021   439652381492021   439652381494021   439652381496821
136617294387035621   617294387038621   617294387039221   617294387044421
137617294387045221   617294387048621   617294387051621   617294387051821
138617294387053621   617294387058421   617294387064221   617294387065621
139617294387068621   617294387069221   617294387069821   617294387070421
140617294387072021   617294387072621   617294387073821   617294387076821
141649378126517621   649378126517821   649378126518221   649378126520821
142649378126523821   649378126525621   649378126526621   649378126528421
143649378126529621   649378126530821   649378126532221   649378126533221
144649378126535221   649378126539421   649378126543621   649378126546021
145649378126546421   649378126549421   649378126550821   649378126555021
146649378126557421   649378126560221   649378126561621   649378126562021
147649378126564621   649378126565821   672091582360421   672091582364221
148672091582364621   672091582367021   672091582368421   672091582369021
149672091582370821   672091582371421   672091582376821   672091582380821
150716805243983221   716805243984821   716805947623621   716805947624621
151716805947629021   716805947629621   716805947630621   716805947633621
152716805947634221   716805947635021   716805947635621   716805947642221
153*/
154
155/****************************************************************************/
156void	M_free_all_rnd()
157{
158if (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/****************************************************************************/
170void	m_apm_set_random_seed(char *ss)
171{
172M_APM   btmp;
173
174if (M_firsttime2)
175  {
176   btmp = M_get_stack_var();
177   m_apm_get_random(btmp);
178   M_restore_stack(1);
179  }
180
181m_apm_set_string(M_rnd_XX, ss);
182}
183/****************************************************************************/
184/*
185 *  compute X = (a * X + c) MOD m       where c = a
186 */
187void	m_apm_get_random(M_APM mrnd)
188{
189
190if (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
208m_apm_multiply(M_rtmp0, M_rnd_XX, M_rnd_aa);
209m_apm_add(M_rtmp1, M_rtmp0, M_rnd_aa);
210m_apm_integer_div_rem(M_rtmp0, M_rnd_XX, M_rtmp1, M_rnd_mm);
211m_apm_copy(mrnd, M_rnd_XX);
212mrnd->m_apm_exponent -= 15;
213}
214/****************************************************************************/
215void	M_reverse_string(char *s)
216{
217int	ct;
218char	ch, *p1, *p2;
219
220if ((ct = strlen(s)) <= 1)
221  return;
222
223p1 = s;
224p2 = s + ct - 1;
225ct /= 2;
226
227while (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 */
247void	M_get_rnd_seed(M_APM mm)
248{
249int              millisec;
250time_t 		 timestamp;
251unsigned long    ul;
252char             ss[32], buf1[48], buf2[32];
253struct timeb     timebuffer;
254M_APM		 atmp;
255
256atmp = M_get_stack_var();
257
258ftime(&timebuffer);
259
260millisec  = (int)timebuffer.millitm;
261timestamp = timebuffer.time;
262ul        = (unsigned long)(timestamp / 7);
263ul       += timestamp + 537;
264strcpy(ss,ctime(&timestamp));        /* convert to string and copy to ss */
265
266sprintf(buf1,"%d",(millisec / 10));
267sprintf(buf2,"%lu",ul);
268
269ss[0] = ss[18];
270ss[1] = ss[17];
271ss[2] = ss[15];
272ss[3] = ss[14];
273ss[4] = ss[12];
274ss[5] = ss[11];
275ss[6] = ss[9];
276ss[7] = ss[23];
277ss[8] = ss[20];
278ss[9] = '\0';
279
280M_reverse_string(buf2);
281strcat(buf1,buf2);
282strcat(buf1,ss);
283
284m_apm_set_string(atmp, buf1);
285atmp->m_apm_exponent = 15;
286m_apm_integer_divide(mm, atmp, MM_One);
287
288M_restore_stack(1);
289}
290/****************************************************************************/
291
292#else
293
294/****************************************************************************/
295/*
296 *  for unix systems : use 'gettimeofday'
297 */
298void	M_get_rnd_seed(M_APM mm)
299{
300unsigned long    sec3;
301long             usec3;
302char             buf1[32], buf2[32];
303M_APM		 atmp;
304
305atmp = M_get_stack_var();
306M_get_microsec(&sec3,&usec3);
307
308sprintf(buf1,"%ld",usec3);
309sprintf(buf2,"%lu",sec3);
310M_reverse_string(buf2);
311strcat(buf1,buf2);
312
313m_apm_set_string(atmp, buf1);
314atmp->m_apm_exponent = 15;
315m_apm_integer_divide(mm, atmp, MM_One);
316
317M_restore_stack(1);
318}
319/****************************************************************************/
320void	 M_get_microsec(unsigned long *sec, long *usec)
321{
322struct timeval time_now;           /* current time for elapsed time check */
323struct timezone time_zone;         /* time zone for gettimeofday call     */
324
325gettimeofday(&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
342void	M_get_rnd_seed(M_APM mm)
343{
344double		 timer0;
345int		 millisec;
346char             *cvi_time, *cvi_date, buf1[64], buf2[32];
347M_APM		 atmp;
348
349atmp = M_get_stack_var();
350
351cvi_date = DateStr();
352cvi_time = TimeStr();
353timer0   = 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
361millisec = (int)(0.01 + 1000.0 * (timer0 - floor(timer0)));
362
363sprintf(buf1, "%d", millisec);
364
365buf2[0]  = cvi_time[6];	/* time format: "HH:MM:SS" */
366buf2[1]  = cvi_time[7];
367buf2[2]  = cvi_time[3];
368buf2[3]  = cvi_time[4];
369buf2[4]  = cvi_time[0];
370buf2[5]  = cvi_time[1];
371
372buf2[6]  = cvi_date[3];	/* date format: "MM-DD-YYYY" */
373buf2[7]  = cvi_date[4];
374buf2[8]  = cvi_date[0];
375buf2[9]  = cvi_date[1];
376buf2[10] = cvi_date[8];
377buf2[11] = cvi_date[9];
378buf2[12] = cvi_date[7];
379
380buf2[13] = '4';
381buf2[14] = '7';
382buf2[15] = '\0';
383
384strcat(buf1, buf2);
385
386m_apm_set_string(atmp, buf1);
387atmp->m_apm_exponent = 15;
388m_apm_integer_divide(mm, atmp, MM_One);
389
390M_restore_stack(1);
391}
392
393#endif
394
395/****************************************************************************/
396
397