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