/*****

*Property of the University of British Columbia (UBC),

*Copyright 2001, by UBC.

*

*By receiving this code, you are agreeing to the following terms:

*1. You will use this code for academic purposes only.

*2. For academic use only, you may distribute the binary or executable code

*   to persons at UBC or the Univ. of Western Australia who have previously 

*   read and agreed to these terms, but you must distribute the SOURCE code

*   with it. 

*3. Each file of source code so distributed must have this header attached.

*4. If the code is revised, the programmer's name and revision date must be added

*   to the Revision List below, as well as the revisions identified in the code.

*5. You will not make this code more widely available via any method such as 

*   publishing in print, email mail-list, usenet posting, website etc.

*6. UBC reserves all rights to this work and all derivative works.

*

*For other proposed purposes please contact:

*The University-Industry Liaison Office 

*IRC Room 331 - 2194 Health Sciences Mall 

*University of British Columbia 

*Vancouver, BC, Canada V6T 1Z3 

*Tel: (604) 822-8580 

*Fax: (604) 822-8589 

*

*or contact: 

*Peter D. Lawrence, Professor at peterl@ece.ubc.ca or 

*Greg Z. Grudic, Assistant Professor, at grudic@cs.colorado.edu

*

*Revision List: 

*Greg Grudic, August 28, 1998.

*****/



/*

%

% File:		rand.c

% Program: 	Random Number functions for Functional Approximation code

%

% Author:	Greg Grudic 

%

*/



#include "b_pc.h"



double Greg_Rand(unsigned short xsubi[3])

{

  float ran4(long *idum);



  double t1;



  static int Initialize = 0;



  static long my_idum;



  if ( Initialize == 0 )

    {

      my_idum = (long)xsubi[0] * (long)xsubi[1] * (long)xsubi[2];



      Initialize = 1;

    }



  t1 = (double)ran4(&my_idum);



  return(t1);

}









#define NITER 4



void psdes(unsigned long *lword, unsigned long *irword)

{

	unsigned long i,ia,ib,iswap,itmph=0,itmpl=0;

	static unsigned long c1[NITER]={

		0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L};

	static unsigned long c2[NITER]={

		0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L};



	for (i=0;i<NITER;i++) {

		ia=(iswap=(*irword)) ^ c1[i];

		itmpl = ia & 0xffff;

		itmph = ia >> 16;

		ib=itmpl*itmpl+ ~(itmph*itmph);

		*irword=(*lword) ^ (((ia = (ib >> 16) |

			((ib & 0xffff) << 16)) ^ c2[i])+itmpl*itmph);

		*lword=iswap;

	}

}

#undef NITER







/* (C) Copr. 1986-92 Numerical Recipes Software 'j''25E#. */

float ran4(long *idum)

{

	void psdes(unsigned long *lword, unsigned long *irword);

	unsigned long irword,itemp,lword;

	static long idums = 0;

#if defined(vax) || defined(_vax_) || defined(__vax__) || defined(VAX)

	static unsigned long jflone = 0x00004080;

	static unsigned long jflmsk = 0xffff007f;

#else

	static unsigned long jflone = 0x3f800000;

	static unsigned long jflmsk = 0x007fffff;

#endif



	if (*idum < 0) {

		idums = -(*idum);

		*idum=1;

	}

	irword=(*idum);

	lword=idums;

	psdes(&lword,&irword);

	itemp=jflone | (jflmsk & irword);

	++(*idum);

	return (*(float *)&itemp)-(float)1.0;

}

/* (C) Copr. 1986-92 Numerical Recipes Software 'j''25E#. */