Home | History | Annotate | Line # | Download | only in intrinsics
      1      1.1  mrg /* Implementation of the IRAND, RAND, and SRAND intrinsics.
      2  1.1.1.4  mrg    Copyright (C) 2004-2024 Free Software Foundation, Inc.
      3      1.1  mrg    Contributed by Steven G. Kargl <kargls (at) comcast.net>.
      4      1.1  mrg 
      5      1.1  mrg This file is part of the GNU Fortran 95 runtime library (libgfortran).
      6      1.1  mrg 
      7      1.1  mrg Libgfortran is free software; you can redistribute it and/or
      8      1.1  mrg modify it under the terms of the GNU General Public
      9      1.1  mrg License as published by the Free Software Foundation; either
     10      1.1  mrg version 3 of the License, or (at your option) any later version.
     11      1.1  mrg 
     12      1.1  mrg Libgfortran is distributed in the hope that it will be useful,
     13      1.1  mrg but WITHOUT ANY WARRANTY; without even the implied warranty of
     14      1.1  mrg MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     15      1.1  mrg GNU General Public License for more details.
     16      1.1  mrg 
     17      1.1  mrg Under Section 7 of GPL version 3, you are granted additional
     18      1.1  mrg permissions described in the GCC Runtime Library Exception, version
     19      1.1  mrg 3.1, as published by the Free Software Foundation.
     20      1.1  mrg 
     21      1.1  mrg You should have received a copy of the GNU General Public License and
     22      1.1  mrg a copy of the GCC Runtime Library Exception along with this program;
     23      1.1  mrg see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
     24      1.1  mrg <http://www.gnu.org/licenses/>.  */
     25      1.1  mrg 
     26      1.1  mrg /* Simple multiplicative congruent algorithm.
     27      1.1  mrg    The period of this generator is approximately 2^31-1, which means that
     28      1.1  mrg    it should not be used for anything serious.  The implementation here
     29      1.1  mrg    is based of an algorithm from  S.K. Park and K.W. Miller, Comm. ACM,
     30      1.1  mrg    31, 1192-1201 (1988).  It is also provided solely for compatibility
     31      1.1  mrg    with G77.  */
     32      1.1  mrg 
     33      1.1  mrg #include "libgfortran.h"
     34      1.1  mrg #include <gthr.h>
     35      1.1  mrg 
     36      1.1  mrg #define GFC_RAND_A	16807
     37      1.1  mrg #define GFC_RAND_M	2147483647
     38      1.1  mrg #define GFC_RAND_M1	(GFC_RAND_M - 1)
     39      1.1  mrg 
     40      1.1  mrg static GFC_UINTEGER_8 rand_seed = 1;
     41      1.1  mrg #ifdef __GTHREAD_MUTEX_INIT
     42      1.1  mrg static __gthread_mutex_t rand_seed_lock = __GTHREAD_MUTEX_INIT;
     43      1.1  mrg #else
     44      1.1  mrg static __gthread_mutex_t rand_seed_lock;
     45      1.1  mrg #endif
     46      1.1  mrg 
     47      1.1  mrg 
     48      1.1  mrg /* Set the seed of the irand generator.  Note 0 is a bad seed.  */
     49      1.1  mrg 
     50      1.1  mrg static void
     51      1.1  mrg srand_internal (GFC_INTEGER_8 i)
     52      1.1  mrg {
     53      1.1  mrg   rand_seed = i ? i : 123459876;
     54      1.1  mrg }
     55      1.1  mrg 
     56      1.1  mrg extern void PREFIX(srand) (GFC_INTEGER_4 *i);
     57      1.1  mrg export_proto_np(PREFIX(srand));
     58      1.1  mrg 
     59      1.1  mrg void
     60      1.1  mrg PREFIX(srand) (GFC_INTEGER_4 *i)
     61      1.1  mrg {
     62      1.1  mrg   __gthread_mutex_lock (&rand_seed_lock);
     63      1.1  mrg   srand_internal (*i);
     64      1.1  mrg   __gthread_mutex_unlock (&rand_seed_lock);
     65      1.1  mrg }
     66      1.1  mrg 
     67      1.1  mrg /* Return an INTEGER in the range [1,GFC_RAND_M-1].  */
     68      1.1  mrg 
     69      1.1  mrg extern GFC_INTEGER_4 irand (GFC_INTEGER_4 *);
     70      1.1  mrg iexport_proto(irand);
     71      1.1  mrg 
     72      1.1  mrg GFC_INTEGER_4
     73      1.1  mrg irand (GFC_INTEGER_4 *i)
     74      1.1  mrg {
     75      1.1  mrg   GFC_INTEGER_4 j;
     76      1.1  mrg   if (i)
     77      1.1  mrg     j = *i;
     78      1.1  mrg   else
     79      1.1  mrg     j = 0;
     80      1.1  mrg 
     81      1.1  mrg   __gthread_mutex_lock (&rand_seed_lock);
     82      1.1  mrg 
     83      1.1  mrg   switch (j)
     84      1.1  mrg   {
     85      1.1  mrg     /* Return the next RN. */
     86      1.1  mrg     case 0:
     87      1.1  mrg       break;
     88      1.1  mrg 
     89      1.1  mrg     /* Reset the RN sequence to system-dependent sequence and return the
     90      1.1  mrg        first value.  */
     91      1.1  mrg     case 1:
     92      1.1  mrg       srand_internal (0);
     93      1.1  mrg       break;
     94      1.1  mrg 
     95      1.1  mrg     /* Seed the RN sequence with j and return the first value.  */
     96      1.1  mrg     default:
     97      1.1  mrg       srand_internal (j);
     98      1.1  mrg       break;
     99      1.1  mrg    }
    100      1.1  mrg 
    101      1.1  mrg    rand_seed = GFC_RAND_A * rand_seed % GFC_RAND_M;
    102      1.1  mrg    j = (GFC_INTEGER_4) rand_seed;
    103      1.1  mrg 
    104      1.1  mrg   __gthread_mutex_unlock (&rand_seed_lock);
    105      1.1  mrg 
    106      1.1  mrg    return j;
    107      1.1  mrg }
    108      1.1  mrg iexport(irand);
    109      1.1  mrg 
    110      1.1  mrg 
    111      1.1  mrg /*  Return a random REAL in the range [0,1).  */
    112      1.1  mrg 
    113      1.1  mrg extern GFC_REAL_4 PREFIX(rand) (GFC_INTEGER_4 *i);
    114      1.1  mrg export_proto_np(PREFIX(rand));
    115      1.1  mrg 
    116      1.1  mrg GFC_REAL_4
    117      1.1  mrg PREFIX(rand) (GFC_INTEGER_4 *i)
    118      1.1  mrg {
    119      1.1  mrg   GFC_UINTEGER_4 mask;
    120      1.1  mrg #if GFC_REAL_4_RADIX == 2
    121      1.1  mrg   mask = ~ (GFC_UINTEGER_4) 0u << (32 - GFC_REAL_4_DIGITS + 1);
    122      1.1  mrg #elif GFC_REAL_4_RADIX == 16
    123      1.1  mrg   mask = ~ (GFC_UINTEGER_4) 0u << ((8 - GFC_REAL_4_DIGITS) * 4 + 1);
    124      1.1  mrg #else
    125      1.1  mrg #error "GFC_REAL_4_RADIX has unknown value"
    126      1.1  mrg #endif
    127      1.1  mrg   return ((GFC_UINTEGER_4) (irand(i) -1) & mask) * (GFC_REAL_4) 0x1.p-31f;
    128      1.1  mrg }
    129      1.1  mrg 
    130      1.1  mrg #ifndef __GTHREAD_MUTEX_INIT
    131      1.1  mrg static void __attribute__((constructor))
    132      1.1  mrg init (void)
    133      1.1  mrg {
    134      1.1  mrg   __GTHREAD_MUTEX_INIT_FUNCTION (&rand_seed_lock);
    135      1.1  mrg }
    136      1.1  mrg #endif
    137