OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libgfortran/] [intrinsics/] [random.c] - Blame information for rev 775

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 733 jeremybenn
/* Implementation of the RANDOM intrinsics
2
   Copyright 2002, 2004, 2005, 2006, 2007, 2009, 2010
3
   Free Software Foundation, Inc.
4
   Contributed by Lars Segerlund <seger@linuxmail.org>
5
   and Steve Kargl.
6
 
7
This file is part of the GNU Fortran 95 runtime library (libgfortran).
8
 
9
Libgfortran is free software; you can redistribute it and/or
10
modify it under the terms of the GNU General Public
11
License as published by the Free Software Foundation; either
12
version 3 of the License, or (at your option) any later version.
13
 
14
Ligbfortran is distributed in the hope that it will be useful,
15
but WITHOUT ANY WARRANTY; without even the implied warranty of
16
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17
GNU General Public License for more details.
18
 
19
Under Section 7 of GPL version 3, you are granted additional
20
permissions described in the GCC Runtime Library Exception, version
21
3.1, as published by the Free Software Foundation.
22
 
23
You should have received a copy of the GNU General Public License and
24
a copy of the GCC Runtime Library Exception along with this program;
25
see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
26
<http://www.gnu.org/licenses/>.  */
27
 
28
#include "libgfortran.h"
29
#include <gthr.h>
30
#include <string.h>
31
 
32
extern void random_r4 (GFC_REAL_4 *);
33
iexport_proto(random_r4);
34
 
35
extern void random_r8 (GFC_REAL_8 *);
36
iexport_proto(random_r8);
37
 
38
extern void arandom_r4 (gfc_array_r4 *);
39
export_proto(arandom_r4);
40
 
41
extern void arandom_r8 (gfc_array_r8 *);
42
export_proto(arandom_r8);
43
 
44
#ifdef HAVE_GFC_REAL_10
45
 
46
extern void random_r10 (GFC_REAL_10 *);
47
iexport_proto(random_r10);
48
 
49
extern void arandom_r10 (gfc_array_r10 *);
50
export_proto(arandom_r10);
51
 
52
#endif
53
 
54
#ifdef HAVE_GFC_REAL_16
55
 
56
extern void random_r16 (GFC_REAL_16 *);
57
iexport_proto(random_r16);
58
 
59
extern void arandom_r16 (gfc_array_r16 *);
60
export_proto(arandom_r16);
61
 
62
#endif
63
 
64
#ifdef __GTHREAD_MUTEX_INIT
65
static __gthread_mutex_t random_lock = __GTHREAD_MUTEX_INIT;
66
#else
67
static __gthread_mutex_t random_lock;
68
#endif
69
 
70
/* Helper routines to map a GFC_UINTEGER_* to the corresponding
71
   GFC_REAL_* types in the range of [0,1).  If GFC_REAL_*_RADIX are 2
72
   or 16, respectively, we mask off the bits that don't fit into the
73
   correct GFC_REAL_*, convert to the real type, then multiply by the
74
   correct offset.  */
75
 
76
 
77
static void
78
rnumber_4 (GFC_REAL_4 *f, GFC_UINTEGER_4 v)
79
{
80
  GFC_UINTEGER_4 mask;
81
#if GFC_REAL_4_RADIX == 2
82
  mask = ~ (GFC_UINTEGER_4) 0u << (32 - GFC_REAL_4_DIGITS);
83
#elif GFC_REAL_4_RADIX == 16
84
  mask = ~ (GFC_UINTEGER_4) 0u << ((8 - GFC_REAL_4_DIGITS) * 4);
85
#else
86
#error "GFC_REAL_4_RADIX has unknown value"
87
#endif
88
  v = v & mask;
89
  *f = (GFC_REAL_4) v * GFC_REAL_4_LITERAL(0x1.p-32);
90
}
91
 
92
static void
93
rnumber_8 (GFC_REAL_8 *f, GFC_UINTEGER_8 v)
94
{
95
  GFC_UINTEGER_8 mask;
96
#if GFC_REAL_8_RADIX == 2
97
  mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_8_DIGITS);
98
#elif GFC_REAL_8_RADIX == 16
99
  mask = ~ (GFC_UINTEGER_8) 0u << (16 - GFC_REAL_8_DIGITS) * 4);
100
#else
101
#error "GFC_REAL_8_RADIX has unknown value"
102
#endif
103
  v = v & mask;
104
  *f = (GFC_REAL_8) v * GFC_REAL_8_LITERAL(0x1.p-64);
105
}
106
 
107
#ifdef HAVE_GFC_REAL_10
108
 
109
static void
110
rnumber_10 (GFC_REAL_10 *f, GFC_UINTEGER_8 v)
111
{
112
  GFC_UINTEGER_8 mask;
113
#if GFC_REAL_10_RADIX == 2
114
  mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_10_DIGITS);
115
#elif GFC_REAL_10_RADIX == 16
116
  mask = ~ (GFC_UINTEGER_10) 0u << ((16 - GFC_REAL_10_DIGITS) * 4);
117
#else
118
#error "GFC_REAL_10_RADIX has unknown value"
119
#endif
120
  v = v & mask;
121
  *f = (GFC_REAL_10) v * GFC_REAL_10_LITERAL(0x1.p-64);
122
}
123
#endif
124
 
125
#ifdef HAVE_GFC_REAL_16
126
 
127
/* For REAL(KIND=16), we only need to mask off the lower bits.  */
128
 
129
static void
130
rnumber_16 (GFC_REAL_16 *f, GFC_UINTEGER_8 v1, GFC_UINTEGER_8 v2)
131
{
132
  GFC_UINTEGER_8 mask;
133
#if GFC_REAL_16_RADIX == 2
134
  mask = ~ (GFC_UINTEGER_8) 0u << (128 - GFC_REAL_16_DIGITS);
135
#elif GFC_REAL_16_RADIX == 16
136
  mask = ~ (GFC_UINTEGER_8) 0u << ((32 - GFC_REAL_16_DIGITS) * 4);
137
#else
138
#error "GFC_REAL_16_RADIX has unknown value"
139
#endif
140
  v2 = v2 & mask;
141
  *f = (GFC_REAL_16) v1 * GFC_REAL_16_LITERAL(0x1.p-64)
142
    + (GFC_REAL_16) v2 * GFC_REAL_16_LITERAL(0x1.p-128);
143
}
144
#endif
145
/* libgfortran previously had a Mersenne Twister, taken from the paper:
146
 
147
        Mersenne Twister:       623-dimensionally equidistributed
148
                                uniform pseudorandom generator.
149
 
150
        by Makoto Matsumoto & Takuji Nishimura
151
        which appeared in the: ACM Transactions on Modelling and Computer
152
        Simulations: Special Issue on Uniform Random Number
153
        Generation. ( Early in 1998 ).
154
 
155
   The Mersenne Twister code was replaced due to
156
 
157
    (1) Simple user specified seeds lead to really bad sequences for
158
        nearly 100000 random numbers.
159
    (2) open(), read(), and close() were not properly declared via header
160
        files.
161
    (3) The global index i was abused and caused unexpected behavior with
162
        GET and PUT.
163
    (4) See PR 15619.
164
 
165
 
166
   libgfortran currently uses George Marsaglia's KISS (Keep It Simple Stupid)
167
   random number generator.  This PRNG combines:
168
 
169
   (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
170
       of 2^32,
171
   (2) A 3-shift shift-register generator with a period of 2^32-1,
172
   (3) Two 16-bit multiply-with-carry generators with a period of
173
       597273182964842497 > 2^59.
174
 
175
   The overall period exceeds 2^123.
176
 
177
   http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
178
 
179
   The above web site has an archive of a newsgroup posting from George
180
   Marsaglia with the statement:
181
 
182
   Subject: Random numbers for C: Improvements.
183
   Date: Fri, 15 Jan 1999 11:41:47 -0500
184
   From: George Marsaglia <geo@stat.fsu.edu>
185
   Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
186
   References: <369B5E30.65A55FD1@stat.fsu.edu>
187
   Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
188
   Lines: 93
189
 
190
   As I hoped, several suggestions have led to
191
   improvements in the code for RNG's I proposed for
192
   use in C. (See the thread "Random numbers for C: Some
193
   suggestions" in previous postings.) The improved code
194
   is listed below.
195
 
196
   A question of copyright has also been raised.  Unlike
197
   DIEHARD, there is no copyright on the code below. You
198
   are free to use it in any way you want, but you may
199
   wish to acknowledge the source, as a courtesy.
200
 
201
"There is no copyright on the code below." included the original
202
KISS algorithm.  */
203
 
204
/* We use three KISS random number generators, with different
205
   seeds.
206
   As a matter of Quality of Implementation, the random numbers
207
   we generate for different REAL kinds, starting from the same
208
   seed, are always the same up to the precision of these types.
209
   We do this by using three generators with different seeds, the
210
   first one always for the most significant bits, the second one
211
   for bits 33..64 (if present in the REAL kind), and the third one
212
   (called twice) for REAL(16).  */
213
 
214
#define GFC_SL(k, n)    ((k)^((k)<<(n)))
215
#define GFC_SR(k, n)    ((k)^((k)>>(n)))
216
 
217
/*   Reference for the seed:
218
   From: "George Marsaglia" <g...@stat.fsu.edu>
219
   Newsgroups: sci.math
220
   Message-ID: <e7CcnWxczriWssCjXTWc3A@comcast.com>
221
 
222
   The KISS RNG uses four seeds, x, y, z, c,
223
   with 0<=x<2^32, 0<y<2^32, 0<=z<2^32, 0<=c<698769069
224
   except that the two pairs
225
   z=0,c=0 and z=2^32-1,c=698769068
226
   should be avoided.  */
227
 
228
/* Any modifications to the seeds that change kiss_size below need to be
229
   reflected in check.c (gfc_check_random_seed) to enable correct
230
   compile-time checking of PUT size for the RANDOM_SEED intrinsic.  */
231
 
232
#define KISS_DEFAULT_SEED_1 123456789, 362436069, 521288629, 316191069
233
#define KISS_DEFAULT_SEED_2 987654321, 458629013, 582859209, 438195021
234
#ifdef HAVE_GFC_REAL_16
235
#define KISS_DEFAULT_SEED_3 573658661, 185639104, 582619469, 296736107
236
#endif
237
 
238
static GFC_UINTEGER_4 kiss_seed[] = {
239
  KISS_DEFAULT_SEED_1,
240
  KISS_DEFAULT_SEED_2,
241
#ifdef HAVE_GFC_REAL_16
242
  KISS_DEFAULT_SEED_3
243
#endif
244
};
245
 
246
static GFC_UINTEGER_4 kiss_default_seed[] = {
247
  KISS_DEFAULT_SEED_1,
248
  KISS_DEFAULT_SEED_2,
249
#ifdef HAVE_GFC_REAL_16
250
  KISS_DEFAULT_SEED_3
251
#endif
252
};
253
 
254
static const GFC_INTEGER_4 kiss_size = sizeof(kiss_seed)/sizeof(kiss_seed[0]);
255
 
256
static GFC_UINTEGER_4 * const kiss_seed_1 = kiss_seed;
257
static GFC_UINTEGER_4 * const kiss_seed_2 = kiss_seed + 4;
258
 
259
#ifdef HAVE_GFC_REAL_16
260
static GFC_UINTEGER_4 * const kiss_seed_3 = kiss_seed + 8;
261
#endif
262
 
263
/* kiss_random_kernel() returns an integer value in the range of
264
   (0, GFC_UINTEGER_4_HUGE].  The distribution of pseudorandom numbers
265
   should be uniform.  */
266
 
267
static GFC_UINTEGER_4
268
kiss_random_kernel(GFC_UINTEGER_4 * seed)
269
{
270
  GFC_UINTEGER_4 kiss;
271
 
272
  seed[0] = 69069 * seed[0] + 1327217885;
273
  seed[1] = GFC_SL(GFC_SR(GFC_SL(seed[1],13),17),5);
274
  seed[2] = 18000 * (seed[2] & 65535) + (seed[2] >> 16);
275
  seed[3] = 30903 * (seed[3] & 65535) + (seed[3] >> 16);
276
  kiss = seed[0] + seed[1] + (seed[2] << 16) + seed[3];
277
 
278
  return kiss;
279
}
280
 
281
/*  This function produces a REAL(4) value from the uniform distribution
282
    with range [0,1).  */
283
 
284
void
285
random_r4 (GFC_REAL_4 *x)
286
{
287
  GFC_UINTEGER_4 kiss;
288
 
289
  __gthread_mutex_lock (&random_lock);
290
  kiss = kiss_random_kernel (kiss_seed_1);
291
  rnumber_4 (x, kiss);
292
  __gthread_mutex_unlock (&random_lock);
293
}
294
iexport(random_r4);
295
 
296
/*  This function produces a REAL(8) value from the uniform distribution
297
    with range [0,1).  */
298
 
299
void
300
random_r8 (GFC_REAL_8 *x)
301
{
302
  GFC_UINTEGER_8 kiss;
303
 
304
  __gthread_mutex_lock (&random_lock);
305
  kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
306
  kiss += kiss_random_kernel (kiss_seed_2);
307
  rnumber_8 (x, kiss);
308
  __gthread_mutex_unlock (&random_lock);
309
}
310
iexport(random_r8);
311
 
312
#ifdef HAVE_GFC_REAL_10
313
 
314
/*  This function produces a REAL(10) value from the uniform distribution
315
    with range [0,1).  */
316
 
317
void
318
random_r10 (GFC_REAL_10 *x)
319
{
320
  GFC_UINTEGER_8 kiss;
321
 
322
  __gthread_mutex_lock (&random_lock);
323
  kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
324
  kiss += kiss_random_kernel (kiss_seed_2);
325
  rnumber_10 (x, kiss);
326
  __gthread_mutex_unlock (&random_lock);
327
}
328
iexport(random_r10);
329
 
330
#endif
331
 
332
/*  This function produces a REAL(16) value from the uniform distribution
333
    with range [0,1).  */
334
 
335
#ifdef HAVE_GFC_REAL_16
336
 
337
void
338
random_r16 (GFC_REAL_16 *x)
339
{
340
  GFC_UINTEGER_8 kiss1, kiss2;
341
 
342
  __gthread_mutex_lock (&random_lock);
343
  kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
344
  kiss1 += kiss_random_kernel (kiss_seed_2);
345
 
346
  kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32;
347
  kiss2 += kiss_random_kernel (kiss_seed_3);
348
 
349
  rnumber_16 (x, kiss1, kiss2);
350
  __gthread_mutex_unlock (&random_lock);
351
}
352
iexport(random_r16);
353
 
354
 
355
#endif
356
/*  This function fills a REAL(4) array with values from the uniform
357
    distribution with range [0,1).  */
358
 
359
void
360
arandom_r4 (gfc_array_r4 *x)
361
{
362
  index_type count[GFC_MAX_DIMENSIONS];
363
  index_type extent[GFC_MAX_DIMENSIONS];
364
  index_type stride[GFC_MAX_DIMENSIONS];
365
  index_type stride0;
366
  index_type dim;
367
  GFC_REAL_4 *dest;
368
  GFC_UINTEGER_4 kiss;
369
  int n;
370
 
371
  dest = x->data;
372
 
373
  dim = GFC_DESCRIPTOR_RANK (x);
374
 
375
  for (n = 0; n < dim; n++)
376
    {
377
      count[n] = 0;
378
      stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
379
      extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
380
      if (extent[n] <= 0)
381
        return;
382
    }
383
 
384
  stride0 = stride[0];
385
 
386
  __gthread_mutex_lock (&random_lock);
387
 
388
  while (dest)
389
    {
390
      /* random_r4 (dest);  */
391
      kiss = kiss_random_kernel (kiss_seed_1);
392
      rnumber_4 (dest, kiss);
393
 
394
      /* Advance to the next element.  */
395
      dest += stride0;
396
      count[0]++;
397
      /* Advance to the next source element.  */
398
      n = 0;
399
      while (count[n] == extent[n])
400
        {
401
          /* When we get to the end of a dimension, reset it and increment
402
             the next dimension.  */
403
          count[n] = 0;
404
          /* We could precalculate these products, but this is a less
405
             frequently used path so probably not worth it.  */
406
          dest -= stride[n] * extent[n];
407
          n++;
408
          if (n == dim)
409
            {
410
              dest = NULL;
411
              break;
412
            }
413
          else
414
            {
415
              count[n]++;
416
              dest += stride[n];
417
            }
418
        }
419
    }
420
  __gthread_mutex_unlock (&random_lock);
421
}
422
 
423
/*  This function fills a REAL(8) array with values from the uniform
424
    distribution with range [0,1).  */
425
 
426
void
427
arandom_r8 (gfc_array_r8 *x)
428
{
429
  index_type count[GFC_MAX_DIMENSIONS];
430
  index_type extent[GFC_MAX_DIMENSIONS];
431
  index_type stride[GFC_MAX_DIMENSIONS];
432
  index_type stride0;
433
  index_type dim;
434
  GFC_REAL_8 *dest;
435
  GFC_UINTEGER_8 kiss;
436
  int n;
437
 
438
  dest = x->data;
439
 
440
  dim = GFC_DESCRIPTOR_RANK (x);
441
 
442
  for (n = 0; n < dim; n++)
443
    {
444
      count[n] = 0;
445
      stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
446
      extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
447
      if (extent[n] <= 0)
448
        return;
449
    }
450
 
451
  stride0 = stride[0];
452
 
453
  __gthread_mutex_lock (&random_lock);
454
 
455
  while (dest)
456
    {
457
      /* random_r8 (dest);  */
458
      kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
459
      kiss += kiss_random_kernel (kiss_seed_2);
460
      rnumber_8 (dest, kiss);
461
 
462
      /* Advance to the next element.  */
463
      dest += stride0;
464
      count[0]++;
465
      /* Advance to the next source element.  */
466
      n = 0;
467
      while (count[n] == extent[n])
468
        {
469
          /* When we get to the end of a dimension, reset it and increment
470
             the next dimension.  */
471
          count[n] = 0;
472
          /* We could precalculate these products, but this is a less
473
             frequently used path so probably not worth it.  */
474
          dest -= stride[n] * extent[n];
475
          n++;
476
          if (n == dim)
477
            {
478
              dest = NULL;
479
              break;
480
            }
481
          else
482
            {
483
              count[n]++;
484
              dest += stride[n];
485
            }
486
        }
487
    }
488
  __gthread_mutex_unlock (&random_lock);
489
}
490
 
491
#ifdef HAVE_GFC_REAL_10
492
 
493
/*  This function fills a REAL(10) array with values from the uniform
494
    distribution with range [0,1).  */
495
 
496
void
497
arandom_r10 (gfc_array_r10 *x)
498
{
499
  index_type count[GFC_MAX_DIMENSIONS];
500
  index_type extent[GFC_MAX_DIMENSIONS];
501
  index_type stride[GFC_MAX_DIMENSIONS];
502
  index_type stride0;
503
  index_type dim;
504
  GFC_REAL_10 *dest;
505
  GFC_UINTEGER_8 kiss;
506
  int n;
507
 
508
  dest = x->data;
509
 
510
  dim = GFC_DESCRIPTOR_RANK (x);
511
 
512
  for (n = 0; n < dim; n++)
513
    {
514
      count[n] = 0;
515
      stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
516
      extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
517
      if (extent[n] <= 0)
518
        return;
519
    }
520
 
521
  stride0 = stride[0];
522
 
523
  __gthread_mutex_lock (&random_lock);
524
 
525
  while (dest)
526
    {
527
      /* random_r10 (dest);  */
528
      kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
529
      kiss += kiss_random_kernel (kiss_seed_2);
530
      rnumber_10 (dest, kiss);
531
 
532
      /* Advance to the next element.  */
533
      dest += stride0;
534
      count[0]++;
535
      /* Advance to the next source element.  */
536
      n = 0;
537
      while (count[n] == extent[n])
538
        {
539
          /* When we get to the end of a dimension, reset it and increment
540
             the next dimension.  */
541
          count[n] = 0;
542
          /* We could precalculate these products, but this is a less
543
             frequently used path so probably not worth it.  */
544
          dest -= stride[n] * extent[n];
545
          n++;
546
          if (n == dim)
547
            {
548
              dest = NULL;
549
              break;
550
            }
551
          else
552
            {
553
              count[n]++;
554
              dest += stride[n];
555
            }
556
        }
557
    }
558
  __gthread_mutex_unlock (&random_lock);
559
}
560
 
561
#endif
562
 
563
#ifdef HAVE_GFC_REAL_16
564
 
565
/*  This function fills a REAL(16) array with values from the uniform
566
    distribution with range [0,1).  */
567
 
568
void
569
arandom_r16 (gfc_array_r16 *x)
570
{
571
  index_type count[GFC_MAX_DIMENSIONS];
572
  index_type extent[GFC_MAX_DIMENSIONS];
573
  index_type stride[GFC_MAX_DIMENSIONS];
574
  index_type stride0;
575
  index_type dim;
576
  GFC_REAL_16 *dest;
577
  GFC_UINTEGER_8 kiss1, kiss2;
578
  int n;
579
 
580
  dest = x->data;
581
 
582
  dim = GFC_DESCRIPTOR_RANK (x);
583
 
584
  for (n = 0; n < dim; n++)
585
    {
586
      count[n] = 0;
587
      stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
588
      extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
589
      if (extent[n] <= 0)
590
        return;
591
    }
592
 
593
  stride0 = stride[0];
594
 
595
  __gthread_mutex_lock (&random_lock);
596
 
597
  while (dest)
598
    {
599
      /* random_r16 (dest);  */
600
      kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32;
601
      kiss1 += kiss_random_kernel (kiss_seed_2);
602
 
603
      kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32;
604
      kiss2 += kiss_random_kernel (kiss_seed_3);
605
 
606
      rnumber_16 (dest, kiss1, kiss2);
607
 
608
      /* Advance to the next element.  */
609
      dest += stride0;
610
      count[0]++;
611
      /* Advance to the next source element.  */
612
      n = 0;
613
      while (count[n] == extent[n])
614
        {
615
          /* When we get to the end of a dimension, reset it and increment
616
             the next dimension.  */
617
          count[n] = 0;
618
          /* We could precalculate these products, but this is a less
619
             frequently used path so probably not worth it.  */
620
          dest -= stride[n] * extent[n];
621
          n++;
622
          if (n == dim)
623
            {
624
              dest = NULL;
625
              break;
626
            }
627
          else
628
            {
629
              count[n]++;
630
              dest += stride[n];
631
            }
632
        }
633
    }
634
  __gthread_mutex_unlock (&random_lock);
635
}
636
 
637
#endif
638
 
639
 
640
 
641
static void
642
scramble_seed (unsigned char *dest, unsigned char *src, int size)
643
{
644
  int i;
645
 
646
  for (i = 0; i < size; i++)
647
    dest[(i % 2) * (size / 2) + i / 2] = src[i];
648
}
649
 
650
 
651
static void
652
unscramble_seed (unsigned char *dest, unsigned char *src, int size)
653
{
654
  int i;
655
 
656
  for (i = 0; i < size; i++)
657
    dest[i] = src[(i % 2) * (size / 2) + i / 2];
658
}
659
 
660
 
661
 
662
/* random_seed is used to seed the PRNG with either a default
663
   set of seeds or user specified set of seeds.  random_seed
664
   must be called with no argument or exactly one argument.  */
665
 
666
void
667
random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
668
{
669
  int i;
670
  unsigned char seed[4*kiss_size];
671
 
672
  __gthread_mutex_lock (&random_lock);
673
 
674
  /* Check that we only have one argument present.  */
675
  if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1)
676
    runtime_error ("RANDOM_SEED should have at most one argument present.");
677
 
678
  /* From the standard: "If no argument is present, the processor assigns
679
     a processor-dependent value to the seed."  */
680
  if (size == NULL && put == NULL && get == NULL)
681
      for (i = 0; i < kiss_size; i++)
682
        kiss_seed[i] = kiss_default_seed[i];
683
 
684
  if (size != NULL)
685
    *size = kiss_size;
686
 
687
  if (put != NULL)
688
    {
689
      /* If the rank of the array is not 1, abort.  */
690
      if (GFC_DESCRIPTOR_RANK (put) != 1)
691
        runtime_error ("Array rank of PUT is not 1.");
692
 
693
      /* If the array is too small, abort.  */
694
      if (GFC_DESCRIPTOR_EXTENT(put,0) < kiss_size)
695
        runtime_error ("Array size of PUT is too small.");
696
 
697
      /*  We copy the seed given by the user.  */
698
      for (i = 0; i < kiss_size; i++)
699
        memcpy (seed + i * sizeof(GFC_UINTEGER_4),
700
                &(put->data[(kiss_size - 1 - i) * GFC_DESCRIPTOR_STRIDE(put,0)]),
701
                sizeof(GFC_UINTEGER_4));
702
 
703
      /* We put it after scrambling the bytes, to paper around users who
704
         provide seeds with quality only in the lower or upper part.  */
705
      scramble_seed ((unsigned char *) kiss_seed, seed, 4*kiss_size);
706
    }
707
 
708
  /* Return the seed to GET data.  */
709
  if (get != NULL)
710
    {
711
      /* If the rank of the array is not 1, abort.  */
712
      if (GFC_DESCRIPTOR_RANK (get) != 1)
713
        runtime_error ("Array rank of GET is not 1.");
714
 
715
      /* If the array is too small, abort.  */
716
      if (GFC_DESCRIPTOR_EXTENT(get,0) < kiss_size)
717
        runtime_error ("Array size of GET is too small.");
718
 
719
      /* Unscramble the seed.  */
720
      unscramble_seed (seed, (unsigned char *) kiss_seed, 4*kiss_size);
721
 
722
      /*  Then copy it back to the user variable.  */
723
      for (i = 0; i < kiss_size; i++)
724
        memcpy (&(get->data[(kiss_size - 1 - i) * GFC_DESCRIPTOR_STRIDE(get,0)]),
725
               seed + i * sizeof(GFC_UINTEGER_4),
726
               sizeof(GFC_UINTEGER_4));
727
    }
728
 
729
  __gthread_mutex_unlock (&random_lock);
730
}
731
iexport(random_seed_i4);
732
 
733
 
734
void
735
random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get)
736
{
737
  int i;
738
 
739
  __gthread_mutex_lock (&random_lock);
740
 
741
  /* Check that we only have one argument present.  */
742
  if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1)
743
    runtime_error ("RANDOM_SEED should have at most one argument present.");
744
 
745
  /* From the standard: "If no argument is present, the processor assigns
746
     a processor-dependent value to the seed."  */
747
  if (size == NULL && put == NULL && get == NULL)
748
      for (i = 0; i < kiss_size; i++)
749
        kiss_seed[i] = kiss_default_seed[i];
750
 
751
  if (size != NULL)
752
    *size = kiss_size / 2;
753
 
754
  if (put != NULL)
755
    {
756
      /* If the rank of the array is not 1, abort.  */
757
      if (GFC_DESCRIPTOR_RANK (put) != 1)
758
        runtime_error ("Array rank of PUT is not 1.");
759
 
760
      /* If the array is too small, abort.  */
761
      if (GFC_DESCRIPTOR_EXTENT(put,0) < kiss_size / 2)
762
        runtime_error ("Array size of PUT is too small.");
763
 
764
      /*  This code now should do correct strides.  */
765
      for (i = 0; i < kiss_size / 2; i++)
766
        memcpy (&kiss_seed[2*i], &(put->data[i * GFC_DESCRIPTOR_STRIDE(put,0)]),
767
                sizeof (GFC_UINTEGER_8));
768
    }
769
 
770
  /* Return the seed to GET data.  */
771
  if (get != NULL)
772
    {
773
      /* If the rank of the array is not 1, abort.  */
774
      if (GFC_DESCRIPTOR_RANK (get) != 1)
775
        runtime_error ("Array rank of GET is not 1.");
776
 
777
      /* If the array is too small, abort.  */
778
      if (GFC_DESCRIPTOR_EXTENT(get,0) < kiss_size / 2)
779
        runtime_error ("Array size of GET is too small.");
780
 
781
      /*  This code now should do correct strides.  */
782
      for (i = 0; i < kiss_size / 2; i++)
783
        memcpy (&(get->data[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &kiss_seed[2*i],
784
                sizeof (GFC_UINTEGER_8));
785
    }
786
 
787
  __gthread_mutex_unlock (&random_lock);
788
}
789
iexport(random_seed_i8);
790
 
791
 
792
#ifndef __GTHREAD_MUTEX_INIT
793
static void __attribute__((constructor))
794
init (void)
795
{
796
  __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock);
797
}
798
#endif

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.