1 /*****************************************************************************
2  *
3  *  Elmer, A Finite Element Software for Multiphysical Problems
4  *
5  *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public
9  * License as published by the Free Software Foundation; either
10  * version 2.1 of the License, or (at your option) any later version.
11  *
12  * This library is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15  * Lesser General Public License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with this library (in file ../LGPL-2.1); if not, write
19  * to the Free Software Foundation, Inc., 51 Franklin Street,
20  * Fifth Floor, Boston, MA  02110-1301  USA
21  *
22  *****************************************************************************/
23 
24 /*******************************************************************************
25  *
26  *     Random number generator.
27  *
28  *******************************************************************************
29  *
30  *                     Author:       Juha Ruokolainen
31  *
32  *                    Address: CSC - IT Center for Science Ltd.
33  *                                Keilaranta 14, P.O. BOX 405
34  *                                  02101 Espoo, Finland
35  *                                  Tel. +358 0 457 2723
36  *                                Telefax: +358 0 457 2302
37  *                              EMail: Juha.Ruokolainen@csc.fi
38  *
39  *                       Date: 30 May 1996
40  *
41  *                Modified by:
42  *
43  *       Date of modification:
44  *
45  ******************************************************************************/
46 /***********************************************************************
47 |
48 |  URAND.C - Last Edited 6. 8. 1988
49 |
50 ***********************************************************************/
51 
52 /*======================================================================
53 |Syntax of the manual pages:
54 |
55 |FUNCTION NAME(...) params ...
56 |
57 $  usage of the function and type of the parameters
58 ?  explane the effects of the function
59 =  return value and the type of value if not of type int
60 @  globals effected directly by this routine
61 !  current known bugs or limitations
62 &  functions called by this function
63 ~  these functions may interest you as an alternative function or
64 |  because they control this function somehow
65 ^=====================================================================*/
66 
67 
68 /*
69  * $Id: urand.c,v 1.4 2006/02/07 10:24:44 jpr Exp $
70  *
71  * $Log: urand.c,v $
72  * Revision 1.4  2006/02/07 10:24:44  jpr
73  * *** empty log message ***
74  *
75  * Revision 1.3  2006/02/07 10:21:42  jpr
76  * Changed visibility of some variables to local scope.
77  *
78  * Revision 1.2  2005/05/27 12:26:22  vierinen
79  * changed header install location
80  *
81  * Revision 1.1.1.1  2005/04/14 13:29:14  vierinen
82  * initial matc automake package
83  *
84  * Revision 1.2  1998/08/01 12:34:57  jpr
85  *
86  * Added Id, started Log.
87  *
88  *
89  */
90 
91 #include "elmer/matc.h"
92 
urand(int * iy)93 double urand(int *iy)
94 /*======================================================================
95 ?  urand is a uniform random number generator based  on  theory  and
96 |  suggestions  given  in  d.e. knuth (1969),  vol  2.   the integer  iy
97 |  should be initialized to an arbitrary integer prior to the first call
98 |  to urand.  the calling program should  not  alter  the  value  of  iy
99 |  between  subsequent calls to urand.  values of urand will be returned
100 |  in the interval (0,1).
101 |  see forsythe, malcolm and moler (1977).
102 ^=====================================================================*/
103 {
104    static double s, halfm;
105    static int  ia, ic, m, mic;
106    static int m2 = 0, itwo = 2;
107 #pragma omp threadprivate (s, halfm, ia, ic, m, mic, m2, itwo)
108 
109   if (m2 == 0)
110   {
111 
112     /* if first entry, compute machine integer word length */
113 
114     m = 1;
115     do
116     {
117       m2 = m;
118       m = itwo * m2;
119     } while(m > m2);
120     halfm = m2;
121 
122     /* compute multiplier and increment for linear congruential method */
123 
124     ia = 8 * (int)(halfm * atan(1.0) / 8.00) + 5;
125     ic = 2*(int)(halfm * (0.50 - sqrt(3.0) / 6.0)) + 1;
126     mic = (m2 - ic) + m2;
127 
128     /* s is the scale factor for converting to floating point */
129 
130     s = 0.5 / halfm;
131 
132   }
133 
134   /* compute next random number */
135 
136   *iy = *iy * ia;
137 
138   /*
139     the following statement is for computers which do not allow
140     integer overflow on addition
141     */
142 
143   if (*iy > mic) *iy = (*iy - m2) - m2;
144 
145   *iy = *iy + ic;
146 
147   /*
148     the following statement is for computers where the
149     word length for addition is greater than for multiplication
150     */
151 
152   if (*iy / 2 > m2) *iy = (*iy - m2) - m2;
153 
154   /*
155     the following statement is for computers where integer
156     overflow affects the sign bit
157     */
158 
159   if (*iy < 0) *iy = (*iy + m2) + m2;
160 
161   return *iy * s;
162 }
163