Article 4853 of comp.lang.perl:
Xref: feenix.metronet.com comp.lang.perl:4853 alt.sources:1678
Newsgroups: comp.lang.perl,alt.sources
Path: feenix.metronet.com!news.ecn.bgu.edu!usenet.ins.cwru.edu!agate!doc.ic.ac.uk!uknet!pipex!uunet!stephsf.com!wengland
From: wengland@stephsf.com (Bill England)
Subject: LCG Random number generator for Perl
References: <1993Jul25.200610.22409@rtsg.mot.com> 
    <TSHIONO.93Jul26200050@cxn22.cv.sony.co.jp> <2324csINNiv@gap.caltech.edu>
Message-ID: <1993Aug7.104053.29725@stephsf.com>
Date: Sat, 7 Aug 1993 15:40:53 GMT
Organization: Stephen Software Systems Inc., Chicago, +1 708 776 0286
X-Bytes: 19234
Lines: 733

In article <2324csINNiv@gap.caltech.edu> ccwf@gg.caltech.edu (Charles Fu) writes:
>(3) Write a better number generator, perhaps based upon perl's rand
>    (see _Numerical Recipes_, for example).
>
>-ccwf
>Charles C. Fu
>Consultant, System Administrator, and Programmer

  Thats a good idea ;-). Fortunately someone has already done it and it
  is really not as easy to get a good set of LCG numbers as one might
  hope.  

  This LCG method has the advantage that it can provide 32 bit numbers
  on 16 bit machines (I've had to hack in in assembler on some systems
  where the 16 bit C could not cut it.)

  Only the perl code is GNU Copyleft the random.c module is Public 
  Domain.

  Read the 'Perl' file below for instructions on how to install this
  into Larry Walls perl.

  Correcting the modulus value in eval.c as described in the 'Perl'
  instructions is required.

#! /bin/sh

# Make a new directory for the random sources, cd to it, and run kits 1
# thru 1 through sh.  When all 1 kits have been run, read README.

echo "This is random 2.1 kit 1 (of 1).  If kit 1 is complete, the line"
echo '"'"End of kit 1 (of 1)"'" will echo at the end.'
echo ""
export PATH || (echo "You didn't use sh, you clunch." ; kill $$)
mkdir  2>/dev/null
echo Extracting README
sed >README <<'!STUFFY!FUNK!' -e 's/X//'
X$Id: README,v 1.1 1992/01/08 05:04:15 wengland Exp wengland $
X
X---------
X
XThis package contains source for a portable 
XLinear Congruential Generator (LCG) generator based on 
Xthe work of Park and Miller 1988.
X
XSEE:
X
X Stephen K. Park and Keith W. Miller. RANDOM NUMBER GENERATORS:
X GOOD ONES ARE HARD TO FIND. Communications of the ACM,
X New York, NY.,October 1988 p.1192
X
XINSTRUCTIONS:
X
X To compile and test random.c under Unix run 'make'.
X
X
XHELP SAVE THE WORLD: 
X
X If you have bugs, fixes, etc.  Please send them to
X support@stephsf.com.  We don't guarantee to fix them 
X but, we'll certainly look at them.
X
X
XKNOWN BUGS:
X
X  It does not work using the Metaware C compiler on a 386
X  CTOS (Unisys) system. I gave up and hacked it in assembler.
X
X
XTHEORY:
X
X  The authors of the original ACM article are continuing research 
X  and work on LCG's and are available via email.
X
X    Stephen K. Park   park@cs.wm.edu
X    Keith W. Miller   miller@cs.wm.edu
X    John Burton       jcburt@cs.wm.edu
X
X  The (random.c) software provided here is in no way sponsored  by
X  them or their university.  If you have questions about the code
X  please send them to me.  If you have questions about the theory
X  or test results please ask them.
X
X
XHave Fun!
X
X  Bill England, wengland@stephsf.com
X  Stephen Software Systems Inc., Chicago/Palatine, +1 708 776 0226
!STUFFY!FUNK!
echo Extracting lotto6
sed >lotto6 <<'!STUFFY!FUNK!' -e 's/X//'
X#!perl
X#
X# $Id: lotto6,v 1.3 1991/12/09 22:36:48 wengland Exp $
X#
X# $Log: lotto6,v $
X# Revision 1.3  1991/12/09  22:36:48  wengland
X# Modified terminology to match WSL terms and descriptions.
X# Added better statistics collection and a standard deviation
X# of percent selected to provied usefull distribution information.
X#
X# Revision 1.2  1991/12/07  06:36:30  wengland
X# Set it up to print out initial auto randomization key.
X# Print date and moved cost information around.
X#
X#
Xeval "exec perl -S $0 $*"
X	if $running_under_some_shell;
X
X## Random Washington State Loto pannel drawing program.
X##
X
X##
X #  (c) Copyright 1991, Stephen Software Systems, Inc.
X #
X #  This software is furnished under the terms of the GNU
X #  public license and is freely distributable.
X #
X #  Using this software will in no way change your odds of 
X #  winning the Washington State Loto or any other game 
X #  of chance.
X #
X #  This software is supplied FREE of charge and comes with
X #  no warrenties expressed or implied.  The author is not
X #  responsible for any damages, monetary or otherwise.
X #
X #               USE AT YOUR OWN RISK.
X##
X
X ## The following number is system dependant.  You must change it
X  # to match the bit size of your systems random number generator.
X  #
X  # (Actually it only needs to be on the same scale  a minor error
X  #   will not make a differance.)
X  # 
X  #   8 bit, 2^7 -1  or 127
X  #  16 bit, 2^15 -1 or 32,767
X  #  32 bit, 2^31 -1 or 2,147,483,647
X  #
X $rand_mod = 2147483647;
X  ##
X   # IT should be noted however that a random number generator
X   # with a Modulus of less than 14M will result in a quite
X   # non random selection.  The only way to fix this is by
X   # building your perl package with a 'good' random number
X   # generator built into it.
X   #
X   # However your odds are still so low that you really 
X   # don't need to worry about it.  :-)
X  ##
X
X## Parameters. 
X # -k Random number key to start with. (Pick your favorite lucky number.
X #     or the prior output of this program labled Next Key:  )
X #    (default, something really wierd and not predictable)
X # -g number of groups to select. (default 1)
X # -n number of numbers per group. (default 10.)
X # -p number of pannels per group. (default 10, one $5 pannel.)
X #
X
Xrequire "getopts.pl";
X    if (!&Getopts('k:g:p:n:')){
X     print "\nThe following options are available:\n\n".
X  "    -k Random number key to start with. (Pick your favorite lucky number.\n".
X  "        or the prior output of this program labled Next Key:  )\n".
X  "       (default, something really wierd and not too predictable)\n".
X  "    -g number of groups to select. (default 1)\n".
X  "    -n number of numbers per group. (default 10.)\n".
X  "    -p number of pannels per group. (default 10, one $5 pannel.)\n"; 
X	exit;
X    }
X
X    print `date`;
X
X    if ($opt_k){
X	print "User defined key=$opt_k\n\n";
X	srand($opt_k);
X    } else {
X	## Pick a fairly random seed based very on 
X	 # the time of day.
X	&randomize();
X    }
X
X    if ($opt_g){
X    }else{
X       $opt_g = 1;
X    }
X    if ($opt_p){
X    }else{
X       $opt_p = 10;
X    }
X    if ($opt_n){
X    }else{
X       $opt_n = 10;
X    }
X
X   printf("%s Groups, %s Numbers per Group, %s Pannels per Group\n",
X		$opt_g, $opt_n, $opt_p );
X
X   printf( "Cost per group: \$%.2f. Cost of all groups: \$%.2f \n\n", 
X       $opt_p * .5, $opt_p * .5 * $opt_g);
X
X## This program uses a group method of selecting
X # lotto pannel.  First 10 numbers are selected 
X # from the 40 available,  then  5 combinations of
X # the 10 numbers are chosen.
X #
X # The idea is that it is a lot more likely to pick
X # 10 possible numbers than to pick only 6 and 
X # chosing combinations of the 10 gives an illusion
X # that you are somehow more likely to win. 
X #
X # The odds are unchanged of course and your chances 
X # of winning are still 14M to one.
X #
X
X    ## Number of members in group must be greator than
X     # 6 and less than 50.
X     #
X    if ($opt_n <7 || $opt_n > 49){
X	print STDERR 
X	"Number of members in group must be greator than 6 and less than 50.\n";
X       exit 1;
X   }
X
X    ## Initialize.
X     #
X    @selection=split(/ /,"01 02 03 04 05 06 07 08 09 10 ".
X			"11 12 13 14 15 16 17 18 19 20 ".
X			"21 22 23 24 25 26 27 28 29 30 ".
X			"31 32 33 34 35 36 37 38 39 40 ".
X			"41 42 43 44 45 46 47 48 49" );
X
X    %Group=();
X    while( (@ttl_Groups=keys(%Group)) < $opt_g ){
X	## Draw one number at random and put it in line of group.
X	 #
X	@r_group=();
X	while( @selection ){
X	    push (@r_group, splice(@selection, rand @selection, 1 ));
X	}
X	@selection = @r_group;
X	@r_group = ();
X
X	$ctr = $opt_n;  # Total numbers in group.
X	while( $ctr--){
X	    $r_group[$ctr] = $selection[$ctr];
X	}
X
X	$group =  join(" ", sort @r_group);
X
X       ## Check for duplicates.
X	#
X	 if( defined $Group{$group} ){
X		printf("Dup Group: %s \n\n", $group);
X		next;
X	    }
X       $Group{$group} = join(" ", @r_group);
X    }
X
X     ## Print out all groups.
X      #
X     foreach $group ( sort(keys %Group )){
X 	&do_a_group( split(/ /, $Group{$group}) );
X    }
X
X    print "Next Key: ", int(rand($rand_mod)), "\n";
X
X##
X #
X##
Xsub do_a_group{
X    local(@r_group)= @_;
X    local($ctr,$ctr1);
X    local(%Group_stats,@stat_group)=(); 
X
X    printf("Group: %s \n", join(" ", sort @r_group));
X
X    foreach $group_member( @r_group){
X       $Group_stats{$group_member} = 0;
X   }
X
X
X    ## Pick a pannel.
X     #
X    %Pannel  = ();
X    @Dup_Pannel=();
X    while( (@ttl_Pannel=keys(%Pannel)) < $opt_p ){
X
X	@r_pannel =();
X	while( @r_group ){
X	    push (@r_pannel, splice(@r_group, rand @r_group, 1 ));
X	}
X	@r_group = @r_pannel;
X	@r_pannel = ();
X
X	$ctr = 6;  # Total numbers in pannel.
X	while( $ctr--){
X	    $r_pannel[$ctr] = $r_group[$ctr];
X	}
X	
X       $pannel =  join(" ", sort @r_pannel);
X
X       ## Check for duplicates.
X	#
X	 if( defined $Pannel{$pannel} ){
X		push(@Dup_Pannel, $pannel);
X		next;
X	    }
X
X       ## Collect group statitistic.
X	#
X       foreach $group_member( @r_pannel ){
X	   $Group_stats{$group_member}++;
X       }
X
X       $Pannel{$pannel} = 0;
X    }
X
X    ##  Print calculate group statistics.
X     #
X    foreach $group_member( sort(keys %Group_stats )){
X       push( @stat_group, 
X	   sprintf("%2s", int($Group_stats{$group_member}/($opt_p*6)*100) ));
X    }
X
X    ## Calculate  and print  s.
X     #
X     # s =   E(x-x')^2
X     #     ---- n-1 ----
X     $xx = 0;
X     foreach $group_stat( @stat_group){ 
X	 $xx += $group_stat; 
X     }
X     $xx /= @stat_group;  # Calculate mean.
X
X     $s = 0;
X     foreach $group_stat( @stat_group ){ 
X	 local($a, $b) = (0,0);
X	 $a = $group_stat - $xx;
X	 $s += $a*$a;
X     }
X
X     $s /= ($#stat_group);  # Calculate s.
X
X    printf(" Gstt: %s  sd=%.2f\n", join(" ", @stat_group), sqrt($s));
X
X
X     ## Print all Dup pannel.
X      #
X     foreach $dup_pannel ( sort @Dup_Pannel ){
X         printf("  Dup: %s \n", $dup_pannel);
X    }
X
X
X     ## Print out all ten pannel.
X      #
X     foreach $pannel ( sort(keys %Pannel )){
X	printf("    P: %s \n", $pannel );
X    }
X    print "\n";
X}
X##
X #
X##
Xsub randomize {
X   $time = (time % 100)*10000;
X   srand($time);
X   local($cycle) = int(rand(1000));
X
X   ## We can still do better.
X    #
X   rand() while $cycle-- > 0;
X
X   srand($cycle = rand($rand_mod));
X   printf("Auto randomization key=%s\n\n", $cycle);
X}
!STUFFY!FUNK!
echo Extracting random.c
sed >random.c <<'!STUFFY!FUNK!' -e 's/X//'
X/*
X $Id: random.c,v $
X
X This program is public domain and was written by William S. England
X (Oct 1988).  It is based on an article by:
X
X Stephen K. Park and Keith W. Miller. RANDOM NUMBER GENERATORS:
X GOOD ONES ARE HARD TO FIND. Communications of the ACM,
X New York, NY.,October 1988 p.1192
X
X Modifications;
X
X $Log: random.c,v $
X
X######
X
X The following is a portable c program for generating random numbers.
X The modulus and multiplier have been extensively tested and should
X not be changed except by someone who is a professional Lehmer generator
X writer.  THIS GENERATOR REPRESENTS THE MINIMUM STANDARD AGAINST WHICH
X OTHER GENERATORS SHOULD BE JUDGED. ("Quote from the referenced article's
X authors. WSE" )
X*/
X
X/*
X**  These are pre-calculated  below to compensate for c 
X**  compilers that may overflow when building the code.
X**
X**  q = (m / a)
X**  r = (m % a)
X*/
X
X/*
X** To build the generator with the original ACM
X** article's numbers use -DORIGINAL_NUMBERS when compiling.
X**
X** Original_numbers are the original published m and q in the
X** ACM article above.  John Burton has furnished numbers for
X** a reportedly better generator.  The new numbers are now
X** used in this program by default.
X*/
X
X#ifndef ORIGINAL_NUMBERS
X#define	m  (unsigned long)2147483647
X#define	q  (unsigned long)44488
X
X#define	a (unsigned int)48271
X#define	r (unsigned int)3399
X
X#define successfulltest 399268537
X#endif
X
X#ifdef ORIGINAL_NUMBERS
X#define	m  (unsigned long)2147483647
X#define	q  (unsigned long)127773
X
X#define	a (unsigned int)16807
X#define	r (unsigned int)2836
X
X#define successfulltest 1043618065
X#endif
X
X/*
X** F(z)	= (az)%m
X**	= az-m(az/m)
X**
X** F(z)  = G(z)+mT(z)
X** G(z)  = a(z%q)- r(z/q)
X** T(z)  = (z/q) - (az/m)
X**
X** F(z)  = a(z%q)- rz/q+ m((z/q) - a(z/m))
X** 	 = a(z%q)- rz/q+ m(z/q) - az
X*/
X
Xunsigned long seed;
X
Xvoid srand( /* unsigned long*/ initial_seed)
Xunsigned long initial_seed;
X{
X    seed = initial_seed; 
X}
X/*
X**
X*/
Xunsigned long rand(/*void*/){
X
Xregister
Xint 	lo, hi, test;
X
X    hi   = seed/q;
X    lo   = seed%q;
X
X    test = a*lo - r*hi;
X
X    if (test > 0)
X	seed = test;
X    else
X	seed = test+ m;
X
X    return seed;
X}
X
X#ifdef TESTRAND
X#include <stdio.h> 
X/*  
X**   The result of running this program should be
X**   "successfulltest".  If this program does not yield this
X**   value then your compiler has not implemented this
X**   program correctly.
X**   
X**   To compile with test option under unix use; 'cc -DTESTRAND random.c'
X**
X**   Be sure to compile without test option for use in applications.
X**   ( Now why did I have to say that ??? )
X*/
X
Xmain(/*void*/)
X{
Xunsigned 
Xlong	n_rand;
X
Xregister int 	i;
Xint	success = 0;
X
X    srand(1);
X
X    for( i = 1; i <= 10001; i++){
X        n_rand = rand();
X
X        if( i> 9998)  
X	    printf("Sequence %5i, Seed= %10i\n", i, seed ); 
X
X	if( i == 10000) 
X	    if( seed == successfulltest ) 
X		success = 1;
X    }
X
X    if (success){
X	printf("The random number generator works correctly.\n\n");
X	exit(0);
X    }else{
X	printf("The random number generator DOES NOT WORK!\n\n");
X	exit(1);
X    }
X}
X#endif
!STUFFY!FUNK!
echo Extracting Perl
sed >Perl <<'!STUFFY!FUNK!' -e 's/X//'
X$Id: Perl,v 1.1 1992/01/08 05:04:15 wengland Exp wengland $
X
X--------
XTips for linking 'random.o' into Larry Walls "Perl".
X
XPerl version 4.19.
X
X  Perl is a useful tool and frequently found in Unix environments.  
X  Many of the functions that Perl uses come directly from the 
X  development environment that Perl was built in.
X
X  Since Perl pulls most of it's functions from existing development
X  libraries it will also get your systems brain-dead random number
X  generator. You should, in my opinion, use the generator provided
X  here or a better one if it's available.
X
X  To use this generator in Perl instead of your systems generator:
X
X  0)  Build and test Perl to make certain that it is working
X      properly.
X
X  1)  Test and create the random.o module in a separate 
X      directory.
X
X  2)  Copy random.o to your perl building directory.
X  
X  3)  Modify eval.c.  Find the line that states #if RANDBITS == 31
X       and change it to read:
X
X       #if RANDBITS == 31
X	       value = rand() * value / 2147483647.0;
X       #else
X
X  4)  Modify config.h.  Find the line that says #define RANDBITS
X      and change it to  so that it says:
X
X      #define RANDBITS 31             /**/
X
X   5) Modify perl.c.   Find case 'v':
X      Add in a line that states you have added random.o to your
X      Perl source.  The case statement should then resemble:
X
X	case 'v':
X	fputs("\nThis is perl, version 4.0\n\n",stdout);
X	fputs(rcsid,stdout);
X	fputs("\nCopyright (c) 1989, 1990, 1991, Larry Wall\n",stdout);
X	fputs("\nRandom number patch, Better LCD installed.\n",stdout);
X    #ifdef MSDOS
X
X   6) Edit makefile.  (Warning there are two of them.  Be sure to
X      edit the correct one which is 'makefile' on  most systems. )
X
X      Add the following line just before 'mallocsrc =': 
X
X      randomobj = random.o
X
X      Add randomobj everywhere there is a mallocobj.
X
X   7) Recompile Perl and test it.
X
X   8) Test perl -v for the version comment.  ( Be sure to run 
X      the version of perl you just created by using ./perl 
X      or similar construct.)
X
X   9) Test the new random number generator by running the following 
X      perl program. ( Be sure to run the version of perl you just
X      created by using ./perl or similar construct.)
X
X----- Cut here -----
X#!perl
Xeval "exec perl -S $0 $*"
X        if $running_under_some_shell;
X
X##
X #
X#$successfulltest 1043618065;  # Optional test number for original 
X			       # generator in ACM.
X
X$successfulltest = 399268537;
X
Xsrand(1);
X
X$seq = 0;
Xwhile ( $seq++ < 10000){
X
X    $seed = rand(2147483647);
X
X    if( $seq == 10000){
X        print $seed, "\n";
X        if( $seed == $successfulltest ){
X            print"random.c is successfully integrated into Perl.\n";
X        }else{
X            print"random.c DID NOT install into Perl correctly! \n";
X        }
X    }
X}
X----- Cut here -----
X
X    10)  Install Perl as per Perl's instructions.
X
XEnd
!STUFFY!FUNK!
echo Extracting Makefile
sed >Makefile <<'!STUFFY!FUNK!' -e 's/X//'
X# $Id: Makefile,v 1.1 1992/01/08 05:04:15 wengland Exp wengland $
X#
X# $Log: Makefile,v $
X# Revision 1.1  1992/01/08  05:04:15  wengland
X# Initial revision
X#
X#
X#  A one program release does not seem to need a Makefile
X#  however, here is one anyway.
X#
X#  Bill England,
X#  Stephen Software Systems Inc., Tacoma/Seattle, +1 800 829 1684
X#
X
X# c compiler options for gcc.
X# comment out if using regular compiler cc.
X
XCC_OPT = -pipe 
XCC=gcc -O
X
XOBJECT_NAME=random.o
X
Xrandom.o:	random.c
X	echo "Testing random number generator:"
X	$(CC) $(CC_OPT) -DTESTRAND random.c
X	./a.out
X	rm a.out
X	$(CC) $(CC_OPT) -c random.c 
X
Xclean:
X	rm -f random.o a.out
!STUFFY!FUNK!
echo Extracting License
sed >License <<'!STUFFY!FUNK!' -e 's/X//'
X$Id: License,v 1.1 1992/01/08 05:04:15 wengland Exp wengland $
X
X--------
X
X  This software is public domain unless otherwise specified 
X  in the program's headers and may be used without
X  licensing or other compensation. 
X
X  THERE ARE NO WARRANTIES EXPRESSED OR IMPLIED WITH REGARDS
X  TO THIS SOFTWARE OR RELATED PACKAGES.
X
X  USE AT YOUR OWN RISK.
X
!STUFFY!FUNK!
echo Extracting MANIFEST
sed >MANIFEST <<'!STUFFY!FUNK!' -e 's/X//'
XMANIFEST                This list of files.
XREADME			General instructions.
XLicense			License agreement.
XMakefile		Makefile to create and test random.c.
XPerl			Instructions for putting random.o into Perl.
Xrandom.c		LCG sequence generator.
Xpatchlevel.h		Patchlevel.
Xlotto6			A Perl program that uses random numbers.
!STUFFY!FUNK!
echo Extracting patchlevel.h
sed >patchlevel.h <<'!STUFFY!FUNK!' -e 's/X//'
X#define PATCHLEVEL 0
!STUFFY!FUNK!
echo " "
echo "End of kit 1 (of 1)"
cat /dev/null >kit1isdone
run=''
config=''
for iskit in 1; do
    if test -f kit${iskit}isdone; then
	run="$run $iskit"
    else
	todo="$todo $iskit"
    fi
done
case $todo in
    '')
	echo "You have run all your kits.  Please read README."
	for combo in *:AA; do
	    if test -f "$combo"; then
		realfile=`basename $combo :AA`
		cat $realfile:[A-Z][A-Z] >$realfile
		rm -rf $realfile:[A-Z][A-Z]
	    fi
	done
	rm -rf kit*isdone
	;;
    *)  echo "You have run$run."
	echo "You still need to run$todo."
	;;
esac
: Someone might mail this, so...
exit

-- 
 +-  Bill England,  wengland@stephsf.COM ------------------------------+
 |   * *      H -> He +24MeV                                           |
 |  * * * ... Oooo, we're having so much fun making itty bitty suns *  |
 |__ * * ______________________________________________________________|