/*
 *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 *  Copyright (C) 2008 - INRIA - Arnaud TORSET
 *
 *  This file must be used under the terms of the CeCILL.
 *  This source file is licensed as described in the file COPYING, which
 *  you should have received as part of this distribution.  The terms
 *  are also available at
 *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 *
 */
 
 
/*
   URAND, A UNIVERSAL RANDOM NUMBER GENERATOR
   BY, MICHAEL A. MALCOLM, CLEVE B. MOLER,
   STAN-CS-73-334, JANUARY 1973,
   COMPUTER SCIENCE  DEPARTMENT,
   School of Humanities and Sciences, STANFORD UNIVERSITY,
   ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/73/334/CS-TR-73-334.pdf

*/



#include "rand.h"
#include <stdio.h>

double	drands(void) {
	int m=1;
	const int itwo=2;
	static int m2=0,halfm=0,ia=0, ic=0, mic=0, iy=0;
	static double s=0.0;
  
      if (m2==0){
	/*  if first entry, compute machine integer word length */
               while (m>m2){
		m2=m;
		m=itwo*m2;
               
		}
               
		halfm = m2;
               
	        /* compute multiplier and increment for linear congruential method */
		ia = 8*(int)(halfm*atan(1.0)/8.0) + 5;
               
                ic = 2*(int)(halfm*(0.5-sqrt(3.0)/6.0)) + 1;
		
                mic = (m2 - ic) + m2;
               
                /* s is the scale factor for converting to floating point */
		s = 0.5/halfm;
                

	}
	
	/* compute next random number	*/
	iy = iy*ia;
      


	/* the following statement is for computers which do not allow
	   integer overflow on addition  */
      if (iy > mic) 
       {
         iy = (iy - m2) - m2;
        
       }  

      iy = iy + ic;

	/* the following statement is for computers where the
	   word length for addition is greater than for multiplication */
      if (iy/2 > m2)
       { 
         iy = (iy - m2) - m2;
        
       }

	/* the following statement is for computers where integer
	   overflow affects the sign bit */
      if (iy < 0) 
       {
        iy = (iy + m2) + m2;
       

       }
      

      return (double)iy*s;
}