gcc/ada/a-nuflra.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT RUNTIME COMPONENTS                          --
   4 --                                                                          --
   5 --            A D A . N U M E R I C S . F L O A T _ R A N D O M             --
   6 --                                                                          --
   7 --                                 B o d y                                  --
   8 --                                                                          --
   9 --                            $Revision: 1.17 $                             --
  10 --                                                                          --
  11 --          Copyright (C) 1992-1998, Free Software Foundation, Inc.         --
  12 --                                                                          --
  13 -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  14 -- terms of the  GNU General Public License as published  by the Free Soft- --
  15 -- ware  Foundation;  either version 2,  or (at your option) any later ver- --
  16 -- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
  17 -- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
  18 -- or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License --
  19 -- for  more details.  You should have  received  a copy of the GNU General --
  20 -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  21 -- to  the Free Software Foundation,  59 Temple Place - Suite 330,  Boston, --
  22 -- MA 02111-1307, USA.                                                      --
  23 --                                                                          --
  24 -- As a special exception,  if other files  instantiate  generics from this --
  25 -- unit, or you link  this unit with other files  to produce an executable, --
  26 -- this  unit  does not  by itself cause  the resulting  executable  to  be --
  27 -- covered  by the  GNU  General  Public  License.  This exception does not --
  28 -- however invalidate  any other reasons why  the executable file  might be --
  29 -- covered by the  GNU Public License.                                      --
  30 --                                                                          --
  31 -- GNAT was originally developed  by the GNAT team at  New York University. --
  32 -- It is now maintained by Ada Core Technologies Inc (http://www.gnat.com). --
  33 --                                                                          --
  34 ------------------------------------------------------------------------------
  35
  36 with Ada.Calendar;
  37
  38 package body Ada.Numerics.Float_Random is
  39
  40    -------------------------
  41    -- Implementation Note --
  42    -------------------------
  43
  44    --  The design of this spec is very awkward, as a result of Ada 95 not
  45    --  permitting in-out parameters for function formals (most naturally
  46    --  Generator values would be passed this way). In pure Ada 95, the only
  47    --  solution is to use the heap and pointers, and, to avoid memory leaks,
  48    --  controlled types.
  49
  50    --  This is awfully heavy, so what we do is to use Unrestricted_Access to
  51    --  get a pointer to the state in the passed Generator. This works because
  52    --  Generator is a limited type and will thus always be passed by reference.
  53
  54    type Pointer is access all State;
  55
  56    -----------------------
  57    -- Local Subprograms --
  58    -----------------------
  59
  60    procedure Euclid (P, Q : in Int; X, Y : out Int; GCD : out Int);
  61
  62    function  Euclid (P, Q : Int) return Int;
  63
  64    function Square_Mod_N (X, N : Int) return Int;
  65
  66    ------------
  67    -- Euclid --
  68    ------------
  69
  70    procedure Euclid (P, Q : in Int; X, Y : out Int; GCD : out Int) is
  71
  72       XT : Int := 1;
  73       YT : Int := 0;
  74
  75       procedure Recur
  76         (P,  Q  : in Int;                 --  a (i-1), a (i)
  77          X,  Y  : in Int;                 --  x (i),   y (i)
  78          XP, YP : in out Int;             --  x (i-1), y (i-1)
  79          GCD    : out Int);
  80
  81       procedure Recur
  82         (P,  Q  : in Int;
  83          X,  Y  : in Int;
  84          XP, YP : in out Int;
  85          GCD    : out Int)
  86       is
  87          Quo : Int  := P / Q;             --  q <-- |_ a (i-1) / a (i) _|
  88          XT  : Int := X;                  --  x (i)
  89          YT  : Int := Y;                  --  y (i)
  90
  91       begin
  92          if P rem Q = 0 then                 --  while does not divide
  93             GCD := Q;
  94             XP  := X;
  95             YP  := Y;
  96          else
  97             Recur (Q, P - Q * Quo, XP - Quo * X, YP - Quo * Y, XT, YT, Quo);
  98
  99             --  a (i) <== a (i)
 100             --  a (i+1) <-- a (i-1) - q*a (i)
 101             --  x (i+1) <-- x (i-1) - q*x (i)
 102             --  y (i+1) <-- y (i-1) - q*y (i)
 103             --  x (i) <== x (i)
 104             --  y (i) <== y (i)
 105
 106             XP  := XT;
 107             YP  := YT;
 108             GCD := Quo;
 109          end if;
 110       end Recur;
 111
 112    --  Start of processing for Euclid
 113
 114    begin
 115       Recur (P, Q, 0, 1, XT, YT, GCD);
 116       X := XT;
 117       Y := YT;
 118    end Euclid;
 119
 120    function Euclid (P, Q : Int) return Int is
 121       X, Y, GCD : Int;
 122
 123    begin
 124       Euclid (P, Q, X, Y, GCD);
 125       return X;
 126    end Euclid;
 127
 128    -----------
 129    -- Image --
 130    -----------
 131
 132    function Image (Of_State : State) return String is
 133    begin
 134       return Int'Image (Of_State.X1) & ',' & Int'Image (Of_State.X2)
 135              & ',' &
 136              Int'Image (Of_State.P)  & ',' & Int'Image (Of_State.Q);
 137    end Image;
 138
 139    ------------
 140    -- Random --
 141    ------------
 142
 143    function Random  (Gen : Generator) return Uniformly_Distributed is
 144       Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access;
 145
 146    begin
 147       Genp.X1 := Square_Mod_N (Genp.X1,  Genp.P);
 148       Genp.X2 := Square_Mod_N (Genp.X2,  Genp.Q);
 149       return
 150         Float ((Flt (((Genp.X2 - Genp.X1) * Genp.X)
 151                   mod Genp.Q) * Flt (Genp.P)
 152           + Flt (Genp.X1)) * Genp.Scl);
 153    end Random;
 154
 155    -----------
 156    -- Reset --
 157    -----------
 158
 159    --  Version that works from given initiator value
 160
 161    procedure Reset (Gen : in Generator; Initiator : in Integer) is
 162       Genp   : constant Pointer := Gen.Gen_State'Unrestricted_Access;
 163       X1, X2 : Int;
 164
 165    begin
 166       X1 := 2 + Int (Initiator) mod (K1 - 3);
 167       X2 := 2 + Int (Initiator) mod (K2 - 3);
 168
 169       --  Eliminate effects of small Initiators.
 170
 171       for J in 1 .. 5 loop
 172          X1 := Square_Mod_N (X1, K1);
 173          X2 := Square_Mod_N (X2, K2);
 174       end loop;
 175
 176       Genp.all :=
 177         (X1  => X1,
 178          X2  => X2,
 179          P   => K1,
 180          Q   => K2,
 181          X   => 1,
 182          Scl => Scal);
 183    end Reset;
 184
 185    --  Version that works from specific saved state
 186
 187    procedure Reset (Gen : Generator; From_State : State) is
 188       Genp : constant Pointer := Gen.Gen_State'Unrestricted_Access;
 189
 190    begin
 191       Genp.all := From_State;
 192    end Reset;
 193
 194    --  Version that works from calendar
 195
 196    procedure Reset (Gen : Generator) is
 197       Genp   : constant Pointer       := Gen.Gen_State'Unrestricted_Access;
 198       Now    : constant Calendar.Time := Calendar.Clock;
 199       X1, X2 : Int;
 200
 201    begin
 202       X1 := Int (Calendar.Year  (Now)) * 12 * 31 +
 203             Int (Calendar.Month (Now)) * 31 +
 204             Int (Calendar.Day   (Now));
 205
 206       X2 := Int (Calendar.Seconds (Now) * Duration (1000.0));
 207
 208       X1 := 2 + X1 mod (K1 - 3);
 209       X2 := 2 + X2 mod (K2 - 3);
 210
 211       --  Eliminate visible effects of same day starts
 212
 213       for J in 1 .. 5 loop
 214          X1 := Square_Mod_N (X1, K1);
 215          X2 := Square_Mod_N (X2, K2);
 216       end loop;
 217
 218
 219       Genp.all :=
 220         (X1  => X1,
 221          X2  => X2,
 222          P   => K1,
 223          Q   => K2,
 224          X   => 1,
 225          Scl => Scal);
 226
 227    end Reset;
 228
 229    ----------
 230    -- Save --
 231    ----------
 232
 233    procedure Save (Gen : in Generator; To_State : out State) is
 234    begin
 235       To_State := Gen.Gen_State;
 236    end Save;
 237
 238    ------------------
 239    -- Square_Mod_N --
 240    ------------------
 241
 242    function Square_Mod_N (X, N : Int) return Int is
 243       Temp : Flt := Flt (X) * Flt (X);
 244       Div  : Int := Int (Temp / Flt (N));
 245
 246    begin
 247       Div := Int (Temp - Flt (Div) * Flt (N));
 248
 249       if Div < 0 then
 250          return Div + N;
 251       else
 252          return Div;
 253       end if;
 254    end Square_Mod_N;
 255
 256    -----------
 257    -- Value --
 258    -----------
 259
 260    function Value (Coded_State : String) return State is
 261       Start : Positive := Coded_State'First;
 262       Stop  : Positive := Coded_State'First;
 263       Outs  : State;
 264
 265    begin
 266       while Coded_State (Stop) /= ',' loop
 267          Stop := Stop + 1;
 268       end loop;
 269
 270       Outs.X1 := Int'Value (Coded_State (Start .. Stop - 1));
 271       Start := Stop + 1;
 272
 273       loop
 274          Stop := Stop + 1;
 275          exit when Coded_State (Stop) = ',';
 276       end loop;
 277
 278       Outs.X2 := Int'Value (Coded_State (Start .. Stop - 1));
 279       Start := Stop + 1;
 280
 281       loop
 282          Stop := Stop + 1;
 283          exit when Coded_State (Stop) = ',';
 284       end loop;
 285
 286       Outs.P   := Int'Value (Coded_State (Start .. Stop - 1));
 287       Outs.Q   := Int'Value (Coded_State (Stop + 1 .. Coded_State'Last));
 288       Outs.X   := Euclid (Outs.P, Outs.Q);
 289       Outs.Scl := 1.0 / (Flt (Outs.P) * Flt (Outs.Q));
 290
 291       --  Now do *some* sanity checks.
 292
 293       if Outs.Q < 31 or else Outs.P < 31
 294         or else Outs.X1 not in 2 .. Outs.P - 1
 295         or else Outs.X2 not in 2 .. Outs.Q - 1
 296       then
 297          raise Constraint_Error;
 298       end if;
 299
 300       return Outs;
 301    end Value;
 302 end Ada.Numerics.Float_Random;