mapm_rnd.c

任意精度的数学库

/* 
 *  M_APM  -  mapm_rnd.c
 *
 *  Copyright (C) 1999   Michael C. Ring
 *
 *  Permission to use, copy, and distribute this software and its
 *  documentation for any purpose with or without fee is hereby granted, 
 *  provided that the above copyright notice appear in all copies and 
 *  that both that copyright notice and this permission notice appear 
 *  in supporting documentation.
 *
 *  Permission to modify the software is granted, but not the right to
 *  distribute the modified code.  Modifications are to be distributed 
 *  as patches to released version.
 *  
 *  This software is provided "as is" without express or implied warranty.
 */

/*
 *      $Id: mapm_rnd.c,v 1.5 1999/09/19 23:32:14 mike Exp $
 *
 *      This file contains the Random Number Generator function.
 *
 *      $Log: mapm_rnd.c,v $
 *      Revision 1.5  1999/09/19 23:32:14  mike
 *      added comments
 *
 *      Revision 1.4  1999/09/18 03:49:25  mike
 *      *** empty log message ***
 *
 *      Revision 1.3  1999/09/18 03:35:36  mike
 *      only prototype get_microsec for non-DOS
 *
 *      Revision 1.2  1999/09/18 02:35:36  mike
 *      delete debug printf's
 *
 *      Revision 1.1  1999/09/18 02:26:52  mike
 *      Initial revision
 */

#include "m_apm_lc.h"

#ifdef MSDOS
#include <time.h>
#include <sys/timeb.h>
#else
#include <sys/time.h>
extern  void	M_get_microsec(unsigned long *, long *);
#endif

extern  void	M_reverse_string(char *);
extern  void    M_get_rnd_seed(M_APM);

static	M_APM   M_rnd_aa;
static  M_APM   M_rnd_mm;
static  M_APM   M_rnd_XX;
static  M_APM   M_rtmp0;
static  M_APM   M_rtmp1;

/*
        Used Knuth's The Art of Computer Programming, Volume 2 as
        the basis. Assuming the random number is X, compute
	(where all the math is performed on integers) :

	X = (a * X + c) MOD m

	From Knuth:

	'm' should be large, at least 2^30 : we use 1.0E+15

	'a' should be between .01m and .99m and not have a simple
	pattern. 'a' should not have any large factors in common 
	with 'm' and (since 'm' is a power of 10) if 'a' MOD 200 
	= 21 then all 'm' different possible values will be 
	generated before 'X' starts to repeat.

	We use 'a' = 649378126517621.

	This is a prime number and also meets 'a' MOD 200 = 21. 
	Commented out below are many potential multipliers that 
	are all prime and meet 'a' MOD 200 = 21.

	There are few restrictions on 'c' except 'c' can have no
	factor in common with 'm', hence we set 'c' = 'a'.

	On the first call, the system time is used to initialize X.
*/

/*
 *  the following constants are all potential multipliers. they are
 *  all prime numbers that also meet the criteria of NUM mod 200 = 21.
 */

/*
439682071525421   439682071528421   439682071529221   439682071529821
439682071530421   439682071532021   439682071538821   439682071539421
439682071540021   439682071547021   439682071551221   439682071553821
439682071555421   439682071557221   439682071558021   439682071558621
439682071559821   439652381461621   439652381465221   439652381465621
439652381466421   439652381467421   439652381468621   439652381470021
439652381471221   439652381477021   439652381484221   439652381488421
439652381491021   439652381492021   439652381494021   439652381496821
617294387035621   617294387038621   617294387039221   617294387044421
617294387045221   617294387048621   617294387051621   617294387051821
617294387053621   617294387058421   617294387064221   617294387065621
617294387068621   617294387069221   617294387069821   617294387070421
617294387072021   617294387072621   617294387073821   617294387076821
649378126517621   649378126517821   649378126518221   649378126520821
649378126523821   649378126525621   649378126526621   649378126528421
649378126529621   649378126530821   649378126532221   649378126533221
649378126535221   649378126539421   649378126543621   649378126546021
649378126546421   649378126549421   649378126550821   649378126555021
*/

/****************************************************************************/
/*
 *  compute X = (a * X + c) MOD m       where c = a
 */
void	m_apm_get_random(mrnd)
M_APM   mrnd;
{
static  int  firsttimer = TRUE;

if (firsttimer)           /* use the system time as the initial seed value */
  {
   firsttimer = FALSE;
   
   M_rnd_aa = m_apm_init();
   M_rnd_XX = m_apm_init();
   M_rnd_mm = m_apm_init();
   M_rtmp0  = m_apm_init();
   M_rtmp1  = m_apm_init();

   /* set the multiplier M_rnd_aa and M_rnd_mm */
   
   m_apm_set_string(M_rnd_aa, "649378126517621");
   m_apm_set_string(M_rnd_mm, "1.0E15");

   M_get_rnd_seed(M_rnd_XX);
  }

m_apm_multiply(M_rtmp0, M_rnd_XX, M_rnd_aa);
m_apm_add(M_rtmp1, M_rtmp0, M_rnd_aa);
m_apm_integer_div_rem(M_rtmp0, M_rnd_XX, M_rtmp1, M_rnd_mm);
m_apm_copy(mrnd, M_rnd_XX);
mrnd->m_apm_exponent -= 15;
}
/****************************************************************************/
void	M_reverse_string(s)
char    *s;
{
int	ct;
char	ch, *p1, *p2;

if ((ct = strlen(s)) <= 1)
  return;

p1 = s;
p2 = s + ct - 1;
ct /= 2;

while (TRUE)
  {
   ch    = *p1;
   *p1++ = *p2;
   *p2-- = ch;

   if (--ct == 0)
     break;
  }
}
/****************************************************************************/

#ifdef MSDOS

/****************************************************************************/
/*
 *  for DOS systems : use 'ftime' 
 */
void	M_get_rnd_seed(mm)
M_APM   mm;
{
int              millisec;
time_t 		 timestamp;
unsigned long    ul;
char             ss[32], buf1[48], buf2[32];
struct timeb     timebuffer;
M_APM		 atmp;

atmp = M_get_stack_var();

ftime(&timebuffer);

millisec  = (int)timebuffer.millitm;    
timestamp = timebuffer.time;
ul        = (unsigned long)(timestamp / 7);
ul       += timestamp + 537;
strcpy(ss,ctime(&timestamp));        /* convert to string and copy to ss */

sprintf(buf1,"%d",(millisec / 10));
sprintf(buf2,"%lu",ul);

ss[0] = ss[18];
ss[1] = ss[17];
ss[2] = ss[15];
ss[3] = ss[14];
ss[4] = ss[12];
ss[5] = ss[11];
ss[6] = ss[9];
ss[7] = ss[23];
ss[8] = ss[20];
ss[9] = '\0';

M_reverse_string(buf2);
strcat(buf1,buf2);
strcat(buf1,ss);

m_apm_set_string(atmp, buf1);
atmp->m_apm_exponent = 15;
m_apm_integer_divide(mm, atmp, MM_One);

M_restore_stack(1);
}
/****************************************************************************/

#else

/****************************************************************************/
/*
 *  for unix systems : use 'gettimeofday'
 */
void	M_get_rnd_seed(mm)
M_APM   mm;
{
unsigned long    sec3;
long             usec3;
char             buf1[32], buf2[32];
M_APM		 atmp;

atmp = M_get_stack_var();
M_get_microsec(&sec3,&usec3);

sprintf(buf1,"%ld",usec3);
sprintf(buf2,"%lu",sec3);
M_reverse_string(buf2);
strcat(buf1,buf2);

m_apm_set_string(atmp, buf1);
atmp->m_apm_exponent = 15;
m_apm_integer_divide(mm, atmp, MM_One);

M_restore_stack(1);
}
/****************************************************************************/
void	 M_get_microsec(sec,usec)
unsigned long *sec;
long     *usec;
{
struct timeval time_now;           /* current time for elapsed time check */
struct timezone time_zone;         /* time zone for gettimeofday call     */

gettimeofday(&time_now, &time_zone);                  /* get current time */

*sec  = time_now.tv_sec;
*usec = time_now.tv_usec;
}
/****************************************************************************/

#endif