+2019-12-16 Arnaud Charlet <charlet@adacore.com>
+
+ * impunit.adb: Add a-nbnbin, a-nbnbre, a-nubinu to Ada 2020
+ units.
+ * Makefile.rtl: Enable new file.
+ * libgnat/a-nbnbin.adb, libgnat/a-nbnbin.ads,
+ libgnat/a-nbnbre.adb, libgnat/a-nbnbre.ads,
+ libgnat/a-nubinu.ads: New files. Provide default standalone
+ implementation of Ada.Numerics.Big_Numbers.Big_* based on
+ System.Generic_Bignum.
+ * libgnat/a-nbnbin__gmp.adb: Alternate implementation of
+ Ada.Numerics.Big_Numbers.Big_Integers based on GMP. Not enabled
+ for now.
+ * libgnat/s-bignum.ads, libgnat/s-bignum.adb: Now a simple
+ wrapper on top of s-genbig.ads.
+ * libgnat/s-genbig.ads, libgnat/s-genbig.adb: New files, making
+ s-bignum generic for reuse in Ada.Numerics.Big_Numbers.
+
2019-12-16 Bob Duff <duff@adacore.com>
* doc/gnat_ugn/building_executable_programs_with_gnat.rst:
a-lliwti$(objext) \
a-llizti$(objext) \
a-locale$(objext) \
+ a-nbnbin$(objext) \
+ a-nbnbre$(objext) \
a-ncelfu$(objext) \
a-ngcefu$(objext) \
a-ngcoar$(objext) \
a-nscefu$(objext) \
a-nscoty$(objext) \
a-nselfu$(objext) \
+ a-nubinu$(objext) \
a-nucoar$(objext) \
a-nucoty$(objext) \
a-nudira$(objext) \
s-flocon$(objext) \
s-fore$(objext) \
s-gearop$(objext) \
+ s-genbig$(objext) \
s-geveop$(objext) \
s-gloloc$(objext) \
s-htable$(objext) \
-- The following units should be used only in Ada 202X mode
Non_Imp_File_Names_2X : constant File_List := (
- 0 => ("a-stteou", T) -- Ada.Strings.Text_Output
- -- ???We use named notation, because there is only one of these so far.
- -- When we add more, we should switch to positional notation, and erase
- -- the "0 =>".
- );
+ ("a-stteou", T), -- Ada.Strings.Text_Output
+ ("a-nubinu", T), -- Ada.Numerics.Big_Numbers
+ ("a-nbnbin", T), -- Ada.Numerics.Big_Numbers.Big_Integers
+ ("a-nbnbre", T)); -- Ada.Numerics.Big_Numbers.Big_Reals
-----------------------
-- Alternative Units --
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- GNAT RUN-TIME COMPONENTS --
+-- --
+-- ADA.NUMERICS.BIG_NUMBERS.BIG_INTEGERS --
+-- --
+-- B o d y --
+-- --
+-- Copyright (C) 2019, Free Software Foundation, Inc. --
+-- --
+-- GNAT is free software; you can redistribute it and/or modify it under --
+-- terms of the GNU General Public License as published by the Free Soft- --
+-- ware Foundation; either version 3, or (at your option) any later ver- --
+-- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
+-- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
+-- or FITNESS FOR A PARTICULAR PURPOSE. --
+-- --
+-- As a special exception under Section 7 of GPL version 3, you are granted --
+-- additional permissions described in the GCC Runtime Library Exception, --
+-- version 3.1, as published by the Free Software Foundation. --
+-- --
+-- You should have received a copy of the GNU General Public License and --
+-- a copy of the GCC Runtime Library Exception along with this program; --
+-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
+-- <http://www.gnu.org/licenses/>. --
+-- --
+-- GNAT was originally developed by the GNAT team at New York University. --
+-- Extensive contributions were provided by Ada Core Technologies Inc. --
+-- --
+------------------------------------------------------------------------------
+
+with Ada.Unchecked_Deallocation;
+with Ada.Characters.Conversions; use Ada.Characters.Conversions;
+
+with Interfaces; use Interfaces;
+
+with System.Generic_Bignums;
+
+package body Ada.Numerics.Big_Numbers.Big_Integers is
+
+ package Bignums is new
+ System.Generic_Bignums (Use_Secondary_Stack => False);
+ use Bignums, System;
+
+ procedure Free is new Ada.Unchecked_Deallocation (Bignum_Data, Bignum);
+
+ function Get_Bignum (Arg : Optional_Big_Integer) return Bignum is
+ (To_Bignum (Arg.Value.C));
+ -- Return the Bignum value stored in Arg
+
+ procedure Set_Bignum (Arg : in out Optional_Big_Integer; Value : Bignum)
+ with Inline;
+ -- Set the Bignum value stored in Arg to Value
+
+ ----------------
+ -- Set_Bignum --
+ ----------------
+
+ procedure Set_Bignum (Arg : in out Optional_Big_Integer; Value : Bignum) is
+ begin
+ Arg.Value.C := To_Address (Value);
+ end Set_Bignum;
+
+ --------------
+ -- Is_Valid --
+ --------------
+
+ function Is_Valid (Arg : Optional_Big_Integer) return Boolean is
+ (Arg.Value.C /= System.Null_Address);
+
+ --------------------------
+ -- Invalid_Big_Integer --
+ --------------------------
+
+ function Invalid_Big_Integer return Optional_Big_Integer is
+ (Value => (Ada.Finalization.Controlled with C => System.Null_Address));
+
+ ---------
+ -- "=" --
+ ---------
+
+ function "=" (L, R : Big_Integer) return Boolean is
+ begin
+ return Big_EQ (Get_Bignum (L), Get_Bignum (R));
+ end "=";
+
+ ---------
+ -- "<" --
+ ---------
+
+ function "<" (L, R : Big_Integer) return Boolean is
+ begin
+ return Big_LT (Get_Bignum (L), Get_Bignum (R));
+ end "<";
+
+ ----------
+ -- "<=" --
+ ----------
+
+ function "<=" (L, R : Big_Integer) return Boolean is
+ begin
+ return Big_LE (Get_Bignum (L), Get_Bignum (R));
+ end "<=";
+
+ ---------
+ -- ">" --
+ ---------
+
+ function ">" (L, R : Big_Integer) return Boolean is
+ begin
+ return Big_GT (Get_Bignum (L), Get_Bignum (R));
+ end ">";
+
+ ----------
+ -- ">=" --
+ ----------
+
+ function ">=" (L, R : Big_Integer) return Boolean is
+ begin
+ return Big_GE (Get_Bignum (L), Get_Bignum (R));
+ end ">=";
+
+ --------------------
+ -- To_Big_Integer --
+ --------------------
+
+ function To_Big_Integer (Arg : Integer) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Set_Bignum (Result, To_Bignum (Long_Long_Integer (Arg)));
+ return Result;
+ end To_Big_Integer;
+
+ ----------------
+ -- To_Integer --
+ ----------------
+
+ function To_Integer (Arg : Big_Integer) return Integer is
+ begin
+ return Integer (From_Bignum (Get_Bignum (Arg)));
+ end To_Integer;
+
+ ------------------------
+ -- Signed_Conversions --
+ ------------------------
+
+ package body Signed_Conversions is
+
+ --------------------
+ -- To_Big_Integer --
+ --------------------
+
+ function To_Big_Integer (Arg : Int) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Set_Bignum (Result, To_Bignum (Long_Long_Integer (Arg)));
+ return Result;
+ end To_Big_Integer;
+
+ ----------------------
+ -- From_Big_Integer --
+ ----------------------
+
+ function From_Big_Integer (Arg : Big_Integer) return Int is
+ begin
+ return Int (From_Bignum (Get_Bignum (Arg)));
+ end From_Big_Integer;
+
+ end Signed_Conversions;
+
+ --------------------------
+ -- Unsigned_Conversions --
+ --------------------------
+
+ package body Unsigned_Conversions is
+
+ --------------------
+ -- To_Big_Integer --
+ --------------------
+
+ function To_Big_Integer (Arg : Int) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Set_Bignum (Result, To_Bignum (Unsigned_64 (Arg)));
+ return Result;
+ end To_Big_Integer;
+
+ ----------------------
+ -- From_Big_Integer --
+ ----------------------
+
+ function From_Big_Integer (Arg : Big_Integer) return Int is
+ begin
+ return Int (From_Bignum (Get_Bignum (Arg)));
+ end From_Big_Integer;
+
+ end Unsigned_Conversions;
+
+ ---------------
+ -- To_String --
+ ---------------
+
+ Hex_Chars : constant array (0 .. 15) of Character := "0123456789ABCDEF";
+
+ function To_String
+ (Arg : Big_Integer; Width : Field := 0; Base : Number_Base := 10)
+ return String
+ is
+ Big_Base : constant Big_Integer := To_Big_Integer (Integer (Base));
+
+ function Add_Base (S : String) return String;
+ -- Add base information if Base /= 10
+
+ function Leading_Padding
+ (Str : String;
+ Min_Length : Field;
+ Char : Character := ' ') return String;
+ -- Return padding of Char concatenated with Str so that the resulting
+ -- string is at least Min_Length long.
+
+ function Image (Arg : Big_Integer) return String;
+ -- Return image of Arg, assuming Arg is positive.
+
+ function Image (N : Natural) return String;
+ -- Return image of N, with no leading space.
+
+ --------------
+ -- Add_Base --
+ --------------
+
+ function Add_Base (S : String) return String is
+ begin
+ if Base = 10 then
+ return S;
+ else
+ return Image (Base) & "#" & S & "#";
+ end if;
+ end Add_Base;
+
+ -----------
+ -- Image --
+ -----------
+
+ function Image (N : Natural) return String is
+ S : constant String := Natural'Image (N);
+ begin
+ return S (2 .. S'Last);
+ end Image;
+
+ function Image (Arg : Big_Integer) return String is
+ begin
+ if Arg < Big_Base then
+ return (1 => Hex_Chars (To_Integer (Arg)));
+ else
+ return Image (Arg / Big_Base)
+ & Hex_Chars (To_Integer (Arg rem Big_Base));
+ end if;
+ end Image;
+
+ ---------------------
+ -- Leading_Padding --
+ ---------------------
+
+ function Leading_Padding
+ (Str : String;
+ Min_Length : Field;
+ Char : Character := ' ') return String is
+ begin
+ return (1 .. Integer'Max (Integer (Min_Length) - Str'Length, 0)
+ => Char) & Str;
+ end Leading_Padding;
+
+ begin
+ if Arg < To_Big_Integer (0) then
+ return Leading_Padding ("-" & Add_Base (Image (-Arg)), Width);
+ else
+ return Leading_Padding (" " & Add_Base (Image (Arg)), Width);
+ end if;
+ end To_String;
+
+ -----------------
+ -- From_String --
+ -----------------
+
+ function From_String (Arg : String) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ -- ??? only support Long_Long_Integer, good enough for now
+ Set_Bignum (Result, To_Bignum (Long_Long_Integer'Value (Arg)));
+ return Result;
+ end From_String;
+
+ ---------------
+ -- Put_Image --
+ ---------------
+
+ procedure Put_Image
+ (Stream : not null access Ada.Streams.Root_Stream_Type'Class;
+ Arg : Big_Integer) is
+ begin
+ Wide_Wide_String'Write (Stream, To_Wide_Wide_String (To_String (Arg)));
+ end Put_Image;
+
+ ---------
+ -- "+" --
+ ---------
+
+ function "+" (L : Big_Integer) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Set_Bignum (Result, new Bignum_Data'(Get_Bignum (L).all));
+ return Result;
+ end "+";
+
+ ---------
+ -- "-" --
+ ---------
+
+ function "-" (L : Big_Integer) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Set_Bignum (Result, Big_Neg (Get_Bignum (L)));
+ return Result;
+ end "-";
+
+ -----------
+ -- "abs" --
+ -----------
+
+ function "abs" (L : Big_Integer) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Set_Bignum (Result, Big_Abs (Get_Bignum (L)));
+ return Result;
+ end "abs";
+
+ ---------
+ -- "+" --
+ ---------
+
+ function "+" (L, R : Big_Integer) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Set_Bignum (Result, Big_Add (Get_Bignum (L), Get_Bignum (R)));
+ return Result;
+ end "+";
+
+ ---------
+ -- "-" --
+ ---------
+
+ function "-" (L, R : Big_Integer) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Set_Bignum (Result, Big_Sub (Get_Bignum (L), Get_Bignum (R)));
+ return Result;
+ end "-";
+
+ ---------
+ -- "*" --
+ ---------
+
+ function "*" (L, R : Big_Integer) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Set_Bignum (Result, Big_Mul (Get_Bignum (L), Get_Bignum (R)));
+ return Result;
+ end "*";
+
+ ---------
+ -- "/" --
+ ---------
+
+ function "/" (L, R : Big_Integer) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Set_Bignum (Result, Big_Div (Get_Bignum (L), Get_Bignum (R)));
+ return Result;
+ end "/";
+
+ -----------
+ -- "mod" --
+ -----------
+
+ function "mod" (L, R : Big_Integer) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Set_Bignum (Result, Big_Mod (Get_Bignum (L), Get_Bignum (R)));
+ return Result;
+ end "mod";
+
+ -----------
+ -- "rem" --
+ -----------
+
+ function "rem" (L, R : Big_Integer) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Set_Bignum (Result, Big_Rem (Get_Bignum (L), Get_Bignum (R)));
+ return Result;
+ end "rem";
+
+ ----------
+ -- "**" --
+ ----------
+
+ function "**" (L : Big_Integer; R : Natural) return Big_Integer is
+ Exp : Bignum := To_Bignum (Long_Long_Integer (R));
+ Result : Optional_Big_Integer;
+ begin
+ Set_Bignum (Result, Big_Exp (Get_Bignum (L), Exp));
+ Free (Exp);
+ return Result;
+ end "**";
+
+ ---------
+ -- Min --
+ ---------
+
+ function Min (L, R : Big_Integer) return Big_Integer is
+ (if L < R then L else R);
+
+ ---------
+ -- Max --
+ ---------
+
+ function Max (L, R : Big_Integer) return Big_Integer is
+ (if L > R then L else R);
+
+ -----------------------------
+ -- Greatest_Common_Divisor --
+ -----------------------------
+
+ function Greatest_Common_Divisor (L, R : Big_Integer) return Big_Positive is
+ function GCD (A, B : Big_Integer) return Big_Integer;
+ -- Recursive internal version
+
+ ---------
+ -- GCD --
+ ---------
+
+ function GCD (A, B : Big_Integer) return Big_Integer is
+ begin
+ if Is_Zero (Get_Bignum (B)) then
+ return A;
+ else
+ return GCD (B, A rem B);
+ end if;
+ end GCD;
+
+ begin
+ return GCD (abs L, abs R);
+ end Greatest_Common_Divisor;
+
+ ------------
+ -- Adjust --
+ ------------
+
+ procedure Adjust (This : in out Controlled_Bignum) is
+ begin
+ if This.C /= System.Null_Address then
+ This.C := To_Address (new Bignum_Data'(To_Bignum (This.C).all));
+ end if;
+ end Adjust;
+
+ --------------
+ -- Finalize --
+ --------------
+
+ procedure Finalize (This : in out Controlled_Bignum) is
+ Tmp : Bignum := To_Bignum (This.C);
+ begin
+ Free (Tmp);
+ This.C := System.Null_Address;
+ end Finalize;
+
+end Ada.Numerics.Big_Numbers.Big_Integers;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- GNAT RUN-TIME COMPONENTS --
+-- --
+-- ADA.NUMERICS.BIG_NUMBERS.BIG_INTEGERS --
+-- --
+-- S p e c --
+-- --
+-- This specification is derived from the Ada Reference Manual for use with --
+-- GNAT. In accordance with the copyright of that document, you can freely --
+-- copy and modify this specification, provided that if you redistribute a --
+-- modified version, any changes that you have made are clearly indicated. --
+-- --
+------------------------------------------------------------------------------
+
+with Ada.Finalization;
+with Ada.Streams;
+
+private with System;
+
+-- Note that some Ada 2020 aspects are commented out since they are not
+-- supported yet.
+
+package Ada.Numerics.Big_Numbers.Big_Integers
+ with Preelaborate
+-- Nonblocking
+is
+ type Optional_Big_Integer is private
+ with Default_Initial_Condition => not Is_Valid (Optional_Big_Integer);
+ -- Integer_Literal => From_String,
+ -- Put_Image => Put_Image;
+
+ function Is_Valid (Arg : Optional_Big_Integer) return Boolean;
+
+ subtype Big_Integer is Optional_Big_Integer
+ with Dynamic_Predicate => Is_Valid (Big_Integer),
+ Predicate_Failure => (raise Constraint_Error);
+
+ function Invalid_Big_Integer return Optional_Big_Integer
+ with Post => not Is_Valid (Invalid_Big_Integer'Result);
+
+ function "=" (L, R : Big_Integer) return Boolean;
+
+ function "<" (L, R : Big_Integer) return Boolean;
+
+ function "<=" (L, R : Big_Integer) return Boolean;
+
+ function ">" (L, R : Big_Integer) return Boolean;
+
+ function ">=" (L, R : Big_Integer) return Boolean;
+
+ function To_Big_Integer (Arg : Integer) return Big_Integer;
+
+ subtype Optional_Big_Positive is Optional_Big_Integer
+ with Dynamic_Predicate =>
+ (not Is_Valid (Optional_Big_Positive))
+ or else (Optional_Big_Positive > To_Big_Integer (0)),
+ Predicate_Failure => (raise Constraint_Error);
+
+ subtype Optional_Big_Natural is Optional_Big_Integer
+ with Dynamic_Predicate =>
+ (not Is_Valid (Optional_Big_Natural))
+ or else (Optional_Big_Natural >= To_Big_Integer (0)),
+ Predicate_Failure => (raise Constraint_Error);
+
+ subtype Big_Positive is Big_Integer
+ with Dynamic_Predicate => Big_Positive > To_Big_Integer (0),
+ Predicate_Failure => (raise Constraint_Error);
+
+ subtype Big_Natural is Big_Integer
+ with Dynamic_Predicate => Big_Natural >= To_Big_Integer (0),
+ Predicate_Failure => (raise Constraint_Error);
+
+ function In_Range (Arg, Low, High : Big_Integer) return Boolean is
+ ((Low <= Arg) and (Arg <= High));
+
+ function To_Integer (Arg : Big_Integer) return Integer
+ with Pre => In_Range (Arg,
+ Low => To_Big_Integer (Integer'First),
+ High => To_Big_Integer (Integer'Last))
+ or else (raise Constraint_Error);
+
+ generic
+ type Int is range <>;
+ package Signed_Conversions is
+
+ function To_Big_Integer (Arg : Int) return Big_Integer;
+
+ function From_Big_Integer (Arg : Big_Integer) return Int
+ with Pre => In_Range (Arg,
+ Low => To_Big_Integer (Int'First),
+ High => To_Big_Integer (Int'Last))
+ or else (raise Constraint_Error);
+
+ end Signed_Conversions;
+
+ generic
+ type Int is mod <>;
+ package Unsigned_Conversions is
+
+ function To_Big_Integer (Arg : Int) return Big_Integer;
+
+ function From_Big_Integer (Arg : Big_Integer) return Int
+ with Pre => In_Range (Arg,
+ Low => To_Big_Integer (Int'First),
+ High => To_Big_Integer (Int'Last))
+ or else (raise Constraint_Error);
+
+ end Unsigned_Conversions;
+
+ function To_String (Arg : Big_Integer;
+ Width : Field := 0;
+ Base : Number_Base := 10) return String
+ with Post => To_String'Result'First = 1;
+
+ function From_String (Arg : String) return Big_Integer;
+
+ procedure Put_Image
+ (Stream : not null access Ada.Streams.Root_Stream_Type'Class;
+ Arg : Big_Integer);
+
+ function "+" (L : Big_Integer) return Big_Integer;
+
+ function "-" (L : Big_Integer) return Big_Integer;
+
+ function "abs" (L : Big_Integer) return Big_Integer;
+
+ function "+" (L, R : Big_Integer) return Big_Integer;
+
+ function "-" (L, R : Big_Integer) return Big_Integer;
+
+ function "*" (L, R : Big_Integer) return Big_Integer;
+
+ function "/" (L, R : Big_Integer) return Big_Integer;
+
+ function "mod" (L, R : Big_Integer) return Big_Integer;
+
+ function "rem" (L, R : Big_Integer) return Big_Integer;
+
+ function "**" (L : Big_Integer; R : Natural) return Big_Integer;
+
+ function Min (L, R : Big_Integer) return Big_Integer;
+
+ function Max (L, R : Big_Integer) return Big_Integer;
+
+ function Greatest_Common_Divisor
+ (L, R : Big_Integer) return Big_Positive
+ with Pre => (L /= To_Big_Integer (0) and R /= To_Big_Integer (0))
+ or else (raise Constraint_Error);
+
+private
+
+ type Controlled_Bignum is new Ada.Finalization.Controlled with record
+ C : System.Address := System.Null_Address;
+ end record;
+
+ procedure Adjust (This : in out Controlled_Bignum);
+ procedure Finalize (This : in out Controlled_Bignum);
+
+ type Optional_Big_Integer is record
+ Value : Controlled_Bignum;
+ end record;
+
+end Ada.Numerics.Big_Numbers.Big_Integers;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- GNAT RUN-TIME COMPONENTS --
+-- --
+-- ADA.NUMERICS.BIG_NUMBERS.BIG_INTEGERS --
+-- --
+-- B o d y --
+-- --
+-- Copyright (C) 2019, Free Software Foundation, Inc. --
+-- --
+-- GNAT is free software; you can redistribute it and/or modify it under --
+-- terms of the GNU General Public License as published by the Free Soft- --
+-- ware Foundation; either version 3, or (at your option) any later ver- --
+-- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
+-- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
+-- or FITNESS FOR A PARTICULAR PURPOSE. --
+-- --
+-- As a special exception under Section 7 of GPL version 3, you are granted --
+-- additional permissions described in the GCC Runtime Library Exception, --
+-- version 3.1, as published by the Free Software Foundation. --
+-- --
+-- You should have received a copy of the GNU General Public License and --
+-- a copy of the GCC Runtime Library Exception along with this program; --
+-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
+-- <http://www.gnu.org/licenses/>. --
+-- --
+-- GNAT was originally developed by the GNAT team at New York University. --
+-- Extensive contributions were provided by Ada Core Technologies Inc. --
+-- --
+------------------------------------------------------------------------------
+
+-- This is the GMP version of this package
+
+with Ada.Unchecked_Conversion;
+with Ada.Unchecked_Deallocation;
+with Interfaces.C; use Interfaces.C;
+with Interfaces.C.Strings; use Interfaces.C.Strings;
+with Ada.Characters.Conversions; use Ada.Characters.Conversions;
+with Ada.Characters.Handling; use Ada.Characters.Handling;
+
+package body Ada.Numerics.Big_Numbers.Big_Integers is
+
+ use System;
+
+ pragma Linker_Options ("-lgmp");
+
+ type mpz_t is record
+ mp_alloc : Integer;
+ mp_size : Integer;
+ mp_d : System.Address;
+ end record;
+ pragma Convention (C, mpz_t);
+ type mpz_t_ptr is access all mpz_t;
+
+ function To_Mpz is new Ada.Unchecked_Conversion (System.Address, mpz_t_ptr);
+ function To_Address is new
+ Ada.Unchecked_Conversion (mpz_t_ptr, System.Address);
+
+ function Get_Mpz (Arg : Optional_Big_Integer) return mpz_t_ptr is
+ (To_Mpz (Arg.Value.C));
+ -- Return the mpz_t value stored in Arg
+
+ procedure Set_Mpz (Arg : in out Optional_Big_Integer; Value : mpz_t_ptr)
+ with Inline;
+ -- Set the mpz_t value stored in Arg to Value
+
+ procedure Allocate (This : in out Optional_Big_Integer) with Inline;
+ -- Allocate an Optional_Big_Integer, including the underlying mpz
+
+ procedure mpz_init_set (ROP : access mpz_t; OP : access constant mpz_t);
+ pragma Import (C, mpz_init_set, "__gmpz_init_set");
+
+ procedure mpz_set (ROP : access mpz_t; OP : access constant mpz_t);
+ pragma Import (C, mpz_set, "__gmpz_set");
+
+ function mpz_cmp (OP1, OP2 : access constant mpz_t) return Integer;
+ pragma Import (C, mpz_cmp, "__gmpz_cmp");
+
+ function mpz_cmp_ui
+ (OP1 : access constant mpz_t; OP2 : unsigned_long) return Integer;
+ pragma Import (C, mpz_cmp_ui, "__gmpz_cmp_ui");
+
+ procedure mpz_set_si (ROP : access mpz_t; OP : long);
+ pragma Import (C, mpz_set_si, "__gmpz_set_si");
+
+ procedure mpz_set_ui (ROP : access mpz_t; OP : unsigned_long);
+ pragma Import (C, mpz_set_ui, "__gmpz_set_ui");
+
+ function mpz_get_si (OP : access constant mpz_t) return long;
+ pragma Import (C, mpz_get_si, "__gmpz_get_si");
+
+ function mpz_get_ui (OP : access constant mpz_t) return unsigned_long;
+ pragma Import (C, mpz_get_ui, "__gmpz_get_ui");
+
+ procedure mpz_neg (ROP : access mpz_t; OP : access constant mpz_t);
+ pragma Import (C, mpz_neg, "__gmpz_neg");
+
+ procedure mpz_sub (ROP : access mpz_t; OP1, OP2 : access constant mpz_t);
+ pragma Import (C, mpz_sub, "__gmpz_sub");
+
+ -------------
+ -- Set_Mpz --
+ -------------
+
+ procedure Set_Mpz (Arg : in out Optional_Big_Integer; Value : mpz_t_ptr) is
+ begin
+ Arg.Value.C := To_Address (Value);
+ end Set_Mpz;
+
+ --------------
+ -- Is_Valid --
+ --------------
+
+ function Is_Valid (Arg : Optional_Big_Integer) return Boolean is
+ (Arg.Value.C /= System.Null_Address);
+
+ --------------------------
+ -- Invalid_Big_Integer --
+ --------------------------
+
+ function Invalid_Big_Integer return Optional_Big_Integer is
+ (Value => (Ada.Finalization.Controlled with C => System.Null_Address));
+
+ ---------
+ -- "=" --
+ ---------
+
+ function "=" (L, R : Big_Integer) return Boolean is
+ begin
+ return mpz_cmp (Get_Mpz (L), Get_Mpz (R)) = 0;
+ end "=";
+
+ ---------
+ -- "<" --
+ ---------
+
+ function "<" (L, R : Big_Integer) return Boolean is
+ begin
+ return mpz_cmp (Get_Mpz (L), Get_Mpz (R)) < 0;
+ end "<";
+
+ ----------
+ -- "<=" --
+ ----------
+
+ function "<=" (L, R : Big_Integer) return Boolean is
+ begin
+ return mpz_cmp (Get_Mpz (L), Get_Mpz (R)) <= 0;
+ end "<=";
+
+ ---------
+ -- ">" --
+ ---------
+
+ function ">" (L, R : Big_Integer) return Boolean is
+ begin
+ return mpz_cmp (Get_Mpz (L), Get_Mpz (R)) > 0;
+ end ">";
+
+ ----------
+ -- ">=" --
+ ----------
+
+ function ">=" (L, R : Big_Integer) return Boolean is
+ begin
+ return mpz_cmp (Get_Mpz (L), Get_Mpz (R)) >= 0;
+ end ">=";
+
+ --------------------
+ -- To_Big_Integer --
+ --------------------
+
+ function To_Big_Integer (Arg : Integer) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Allocate (Result);
+ mpz_set_si (Get_Mpz (Result), long (Arg));
+ return Result;
+ end To_Big_Integer;
+
+ ----------------
+ -- To_Integer --
+ ----------------
+
+ function To_Integer (Arg : Big_Integer) return Integer is
+ begin
+ return Integer (mpz_get_si (Get_Mpz (Arg)));
+ end To_Integer;
+
+ ------------------------
+ -- Signed_Conversions --
+ ------------------------
+
+ package body Signed_Conversions is
+
+ --------------------
+ -- To_Big_Integer --
+ --------------------
+
+ function To_Big_Integer (Arg : Int) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Allocate (Result);
+ mpz_set_si (Get_Mpz (Result), long (Arg));
+ return Result;
+ end To_Big_Integer;
+
+ ----------------------
+ -- From_Big_Integer --
+ ----------------------
+
+ function From_Big_Integer (Arg : Big_Integer) return Int is
+ begin
+ return Int (mpz_get_si (Get_Mpz (Arg)));
+ end From_Big_Integer;
+
+ end Signed_Conversions;
+
+ --------------------------
+ -- Unsigned_Conversions --
+ --------------------------
+
+ package body Unsigned_Conversions is
+
+ --------------------
+ -- To_Big_Integer --
+ --------------------
+
+ function To_Big_Integer (Arg : Int) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Allocate (Result);
+ mpz_set_ui (Get_Mpz (Result), unsigned_long (Arg));
+ return Result;
+ end To_Big_Integer;
+
+ ----------------------
+ -- From_Big_Integer --
+ ----------------------
+
+ function From_Big_Integer (Arg : Big_Integer) return Int is
+ begin
+ return Int (mpz_get_ui (Get_Mpz (Arg)));
+ end From_Big_Integer;
+
+ end Unsigned_Conversions;
+
+ ---------------
+ -- To_String --
+ ---------------
+
+ function To_String
+ (Arg : Big_Integer; Width : Field := 0; Base : Number_Base := 10)
+ return String
+ is
+ function mpz_get_str
+ (STR : System.Address;
+ BASE : Integer;
+ OP : access constant mpz_t) return chars_ptr;
+ pragma Import (C, mpz_get_str, "__gmpz_get_str");
+
+ function mpz_sizeinbase
+ (this : access constant mpz_t; base : Integer) return size_t;
+ pragma Import (C, mpz_sizeinbase, "__gmpz_sizeinbase");
+
+ function Add_Base (S : String) return String;
+ -- Add base information if Base /= 10
+
+ function Leading_Padding
+ (Str : String;
+ Min_Length : Field;
+ Char : Character := ' ') return String;
+ -- Return padding of Char concatenated with Str so that the resulting
+ -- string is at least Min_Length long.
+
+ function Image (N : Natural) return String;
+ -- Return image of N, with no leading space.
+
+ --------------
+ -- Add_Base --
+ --------------
+
+ function Add_Base (S : String) return String is
+ begin
+ if Base = 10 then
+ return S;
+ else
+ return Image (Base) & "#" & To_Upper (S) & "#";
+ end if;
+ end Add_Base;
+
+ -----------
+ -- Image --
+ -----------
+
+ function Image (N : Natural) return String is
+ S : constant String := Natural'Image (N);
+ begin
+ return S (2 .. S'Last);
+ end Image;
+
+ ---------------------
+ -- Leading_Padding --
+ ---------------------
+
+ function Leading_Padding
+ (Str : String;
+ Min_Length : Field;
+ Char : Character := ' ') return String is
+ begin
+ return (1 .. Integer'Max (Integer (Min_Length) - Str'Length, 0)
+ => Char) & Str;
+ end Leading_Padding;
+
+ Number_Digits : constant Integer :=
+ Integer (mpz_sizeinbase (Get_Mpz (Arg), Integer (abs Base)));
+
+ Buffer : aliased String (1 .. Number_Digits + 2);
+ -- The correct number to allocate is 2 more than Number_Digits in order
+ -- to handle a possible minus sign and the null-terminator.
+
+ Result : constant chars_ptr :=
+ mpz_get_str (Buffer'Address, Integer (Base), Get_Mpz (Arg));
+ S : constant String := Value (Result);
+
+ begin
+ if S (1) = '-' then
+ return Leading_Padding ("-" & Add_Base (S (2 .. S'Last)), Width);
+ else
+ return Leading_Padding (" " & Add_Base (S), Width);
+ end if;
+ end To_String;
+
+ -----------------
+ -- From_String --
+ -----------------
+
+ function From_String (Arg : String) return Big_Integer is
+ function mpz_set_str
+ (this : access mpz_t;
+ str : System.Address;
+ base : Integer := 10) return Integer;
+ pragma Import (C, mpz_set_str, "__gmpz_set_str");
+
+ Result : Optional_Big_Integer;
+ First : Natural;
+ Last : Natural;
+ Base : Natural;
+
+ begin
+ Allocate (Result);
+
+ if Arg (Arg'Last) /= '#' then
+
+ -- Base 10 number
+
+ First := Arg'First;
+ Last := Arg'Last;
+ Base := 10;
+ else
+ -- Compute the xx base in a xx#yyyyy# number
+
+ if Arg'Length < 4 then
+ raise Constraint_Error;
+ end if;
+
+ First := 0;
+ Last := Arg'Last - 1;
+
+ for J in Arg'First + 1 .. Last loop
+ if Arg (J) = '#' then
+ First := J;
+ exit;
+ end if;
+ end loop;
+
+ if First = 0 then
+ raise Constraint_Error;
+ end if;
+
+ Base := Natural'Value (Arg (Arg'First .. First - 1));
+ First := First + 1;
+ end if;
+
+ declare
+ Str : aliased String (1 .. Last - First + 2);
+ Index : Natural := 0;
+ begin
+ -- Strip underscores
+
+ for J in First .. Last loop
+ if Arg (J) /= '_' then
+ Index := Index + 1;
+ Str (Index) := Arg (J);
+ end if;
+ end loop;
+
+ Index := Index + 1;
+ Str (Index) := ASCII.NUL;
+
+ if mpz_set_str (Get_Mpz (Result), Str'Address, Base) /= 0 then
+ raise Constraint_Error;
+ end if;
+ end;
+
+ return Result;
+ end From_String;
+
+ ---------------
+ -- Put_Image --
+ ---------------
+
+ procedure Put_Image
+ (Stream : not null access Ada.Streams.Root_Stream_Type'Class;
+ Arg : Big_Integer) is
+ begin
+ Wide_Wide_String'Write (Stream, To_Wide_Wide_String (To_String (Arg)));
+ end Put_Image;
+
+ ---------
+ -- "+" --
+ ---------
+
+ function "+" (L : Big_Integer) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Set_Mpz (Result, new mpz_t);
+ mpz_init_set (Get_Mpz (Result), Get_Mpz (L));
+ return Result;
+ end "+";
+
+ ---------
+ -- "-" --
+ ---------
+
+ function "-" (L : Big_Integer) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Allocate (Result);
+ mpz_neg (Get_Mpz (Result), Get_Mpz (L));
+ return Result;
+ end "-";
+
+ -----------
+ -- "abs" --
+ -----------
+
+ function "abs" (L : Big_Integer) return Big_Integer is
+ procedure mpz_abs (ROP : access mpz_t; OP : access constant mpz_t);
+ pragma Import (C, mpz_abs, "__gmpz_abs");
+
+ Result : Optional_Big_Integer;
+ begin
+ Allocate (Result);
+ mpz_abs (Get_Mpz (Result), Get_Mpz (L));
+ return Result;
+ end "abs";
+
+ ---------
+ -- "+" --
+ ---------
+
+ function "+" (L, R : Big_Integer) return Big_Integer is
+ procedure mpz_add
+ (ROP : access mpz_t; OP1, OP2 : access constant mpz_t);
+ pragma Import (C, mpz_add, "__gmpz_add");
+
+ Result : Optional_Big_Integer;
+
+ begin
+ Allocate (Result);
+ mpz_add (Get_Mpz (Result), Get_Mpz (L), Get_Mpz (R));
+ return Result;
+ end "+";
+
+ ---------
+ -- "-" --
+ ---------
+
+ function "-" (L, R : Big_Integer) return Big_Integer is
+ Result : Optional_Big_Integer;
+ begin
+ Allocate (Result);
+ mpz_sub (Get_Mpz (Result), Get_Mpz (L), Get_Mpz (R));
+ return Result;
+ end "-";
+
+ ---------
+ -- "*" --
+ ---------
+
+ function "*" (L, R : Big_Integer) return Big_Integer is
+ procedure mpz_mul
+ (ROP : access mpz_t; OP1, OP2 : access constant mpz_t);
+ pragma Import (C, mpz_mul, "__gmpz_mul");
+
+ Result : Optional_Big_Integer;
+
+ begin
+ Allocate (Result);
+ mpz_mul (Get_Mpz (Result), Get_Mpz (L), Get_Mpz (R));
+ return Result;
+ end "*";
+
+ ---------
+ -- "/" --
+ ---------
+
+ function "/" (L, R : Big_Integer) return Big_Integer is
+ procedure mpz_tdiv_q (Q : access mpz_t; N, D : access constant mpz_t);
+ pragma Import (C, mpz_tdiv_q, "__gmpz_tdiv_q");
+ begin
+ if mpz_cmp_ui (Get_Mpz (R), 0) = 0 then
+ raise Constraint_Error;
+ end if;
+
+ declare
+ Result : Optional_Big_Integer;
+ begin
+ Allocate (Result);
+ mpz_tdiv_q (Get_Mpz (Result), Get_Mpz (L), Get_Mpz (R));
+ return Result;
+ end;
+ end "/";
+
+ -----------
+ -- "mod" --
+ -----------
+
+ function "mod" (L, R : Big_Integer) return Big_Integer is
+ procedure mpz_mod (R : access mpz_t; N, D : access constant mpz_t);
+ pragma Import (C, mpz_mod, "__gmpz_mod");
+ -- result is always non-negative
+
+ L_Negative, R_Negative : Boolean;
+
+ begin
+ if mpz_cmp_ui (Get_Mpz (R), 0) = 0 then
+ raise Constraint_Error;
+ end if;
+
+ declare
+ Result : Optional_Big_Integer;
+ begin
+ Allocate (Result);
+ L_Negative := mpz_cmp_ui (Get_Mpz (L), 0) < 0;
+ R_Negative := mpz_cmp_ui (Get_Mpz (R), 0) < 0;
+
+ if not (L_Negative or R_Negative) then
+ mpz_mod (Get_Mpz (Result), Get_Mpz (L), Get_Mpz (R));
+ else
+ -- The GMP library provides operators defined by C semantics, but
+ -- the semantics of Ada's mod operator are not the same as C's
+ -- when negative values are involved. We do the following to
+ -- implement the required Ada semantics.
+
+ declare
+ Temp_Left : Big_Integer;
+ Temp_Right : Big_Integer;
+ Temp_Result : Big_Integer;
+
+ begin
+ Allocate (Temp_Result);
+ Set_Mpz (Temp_Left, new mpz_t);
+ Set_Mpz (Temp_Right, new mpz_t);
+ mpz_init_set (Get_Mpz (Temp_Left), Get_Mpz (L));
+ mpz_init_set (Get_Mpz (Temp_Right), Get_Mpz (R));
+
+ if L_Negative then
+ mpz_neg (Get_Mpz (Temp_Left), Get_Mpz (Temp_Left));
+ end if;
+
+ if R_Negative then
+ mpz_neg (Get_Mpz (Temp_Right), Get_Mpz (Temp_Right));
+ end if;
+
+ -- now both Temp_Left and Temp_Right are nonnegative
+
+ mpz_mod (Get_Mpz (Temp_Result),
+ Get_Mpz (Temp_Left),
+ Get_Mpz (Temp_Right));
+
+ if mpz_cmp_ui (Get_Mpz (Temp_Result), 0) = 0 then
+ -- if Temp_Result is zero we are done
+ mpz_set (Get_Mpz (Result), Get_Mpz (Temp_Result));
+
+ elsif L_Negative then
+ if R_Negative then
+ mpz_neg (Get_Mpz (Result), Get_Mpz (Temp_Result));
+ else -- L is negative but R is not
+ mpz_sub (Get_Mpz (Result),
+ Get_Mpz (Temp_Right),
+ Get_Mpz (Temp_Result));
+ end if;
+ else
+ pragma Assert (R_Negative);
+ mpz_sub (Get_Mpz (Result),
+ Get_Mpz (Temp_Result),
+ Get_Mpz (Temp_Right));
+ end if;
+ end;
+ end if;
+
+ return Result;
+ end;
+ end "mod";
+
+ -----------
+ -- "rem" --
+ -----------
+
+ function "rem" (L, R : Big_Integer) return Big_Integer is
+ procedure mpz_tdiv_r (R : access mpz_t; N, D : access constant mpz_t);
+ pragma Import (C, mpz_tdiv_r, "__gmpz_tdiv_r");
+ -- R will have the same sign as N.
+
+ begin
+ if mpz_cmp_ui (Get_Mpz (R), 0) = 0 then
+ raise Constraint_Error;
+ end if;
+
+ declare
+ Result : Optional_Big_Integer;
+ begin
+ Allocate (Result);
+ mpz_tdiv_r (R => Get_Mpz (Result),
+ N => Get_Mpz (L),
+ D => Get_Mpz (R));
+ -- the result takes the sign of N, as required by the RM
+
+ return Result;
+ end;
+ end "rem";
+
+ ----------
+ -- "**" --
+ ----------
+
+ function "**" (L : Big_Integer; R : Natural) return Big_Integer is
+ procedure mpz_pow_ui (ROP : access mpz_t;
+ BASE : access constant mpz_t;
+ EXP : unsigned_long);
+ pragma Import (C, mpz_pow_ui, "__gmpz_pow_ui");
+
+ Result : Optional_Big_Integer;
+
+ begin
+ Allocate (Result);
+ mpz_pow_ui (Get_Mpz (Result), Get_Mpz (L), unsigned_long (R));
+ return Result;
+ end "**";
+
+ ---------
+ -- Min --
+ ---------
+
+ function Min (L, R : Big_Integer) return Big_Integer is
+ (if L < R then L else R);
+
+ ---------
+ -- Max --
+ ---------
+
+ function Max (L, R : Big_Integer) return Big_Integer is
+ (if L > R then L else R);
+
+ -----------------------------
+ -- Greatest_Common_Divisor --
+ -----------------------------
+
+ function Greatest_Common_Divisor (L, R : Big_Integer) return Big_Integer is
+ procedure mpz_gcd
+ (ROP : access mpz_t; Op1, Op2 : access constant mpz_t);
+ pragma Import (C, mpz_gcd, "__gmpz_gcd");
+
+ Result : Optional_Big_Integer;
+
+ begin
+ Allocate (Result);
+ mpz_gcd (Get_Mpz (Result), Get_Mpz (L), Get_Mpz (R));
+ return Result;
+ end Greatest_Common_Divisor;
+
+ --------------
+ -- Allocate --
+ --------------
+
+ procedure Allocate (This : in out Optional_Big_Integer) is
+ procedure mpz_init (this : access mpz_t);
+ pragma Import (C, mpz_init, "__gmpz_init");
+ begin
+ Set_Mpz (This, new mpz_t);
+ mpz_init (Get_Mpz (This));
+ end Allocate;
+
+ ------------
+ -- Adjust --
+ ------------
+
+ procedure Adjust (This : in out Controlled_Bignum) is
+ Value : constant mpz_t_ptr := To_Mpz (This.C);
+ begin
+ if Value /= null then
+ This.C := To_Address (new mpz_t);
+ mpz_init_set (To_Mpz (This.C), Value);
+ end if;
+ end Adjust;
+
+ --------------
+ -- Finalize --
+ --------------
+
+ procedure Finalize (This : in out Controlled_Bignum) is
+ procedure Free is new Ada.Unchecked_Deallocation (mpz_t, mpz_t_ptr);
+
+ procedure mpz_clear (this : access mpz_t);
+ pragma Import (C, mpz_clear, "__gmpz_clear");
+
+ Mpz : mpz_t_ptr;
+
+ begin
+ if This.C /= System.Null_Address then
+ Mpz := To_Mpz (This.C);
+ mpz_clear (Mpz);
+ Free (Mpz);
+ This.C := System.Null_Address;
+ end if;
+ end Finalize;
+
+end Ada.Numerics.Big_Numbers.Big_Integers;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- GNAT RUN-TIME COMPONENTS --
+-- --
+-- ADA.NUMERICS.BIG_NUMBERS.BIG_REALS --
+-- --
+-- B o d y --
+-- --
+-- Copyright (C) 2019, Free Software Foundation, Inc. --
+-- --
+-- GNAT is free software; you can redistribute it and/or modify it under --
+-- terms of the GNU General Public License as published by the Free Soft- --
+-- ware Foundation; either version 3, or (at your option) any later ver- --
+-- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
+-- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
+-- or FITNESS FOR A PARTICULAR PURPOSE. --
+-- --
+-- As a special exception under Section 7 of GPL version 3, you are granted --
+-- additional permissions described in the GCC Runtime Library Exception, --
+-- version 3.1, as published by the Free Software Foundation. --
+-- --
+-- You should have received a copy of the GNU General Public License and --
+-- a copy of the GCC Runtime Library Exception along with this program; --
+-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
+-- <http://www.gnu.org/licenses/>. --
+-- --
+-- GNAT was originally developed by the GNAT team at New York University. --
+-- Extensive contributions were provided by Ada Core Technologies Inc. --
+-- --
+------------------------------------------------------------------------------
+
+-- This is the default version of this package, based on Big_Integers only.
+
+with Ada.Characters.Conversions; use Ada.Characters.Conversions;
+
+package body Ada.Numerics.Big_Numbers.Big_Reals is
+
+ use Big_Integers;
+
+ procedure Normalize (Arg : in out Big_Real);
+ -- Normalize Arg by ensuring that Arg.Den is always positive and that
+ -- Arg.Num and Arg.Den always have a GCD of 1.
+
+ --------------
+ -- Is_Valid --
+ --------------
+
+ function Is_Valid (Arg : Optional_Big_Real) return Boolean is
+ (Is_Valid (Arg.Num) and then Is_Valid (Arg.Den));
+
+ -----------------
+ -- No_Big_Real --
+ -----------------
+
+ function No_Big_Real return Optional_Big_Real is
+ (Num => Invalid_Big_Integer, Den => Invalid_Big_Integer);
+
+ ---------
+ -- "/" --
+ ---------
+
+ function "/" (Num, Den : Big_Integer) return Big_Real is
+ Result : Optional_Big_Real;
+ begin
+ if Den = To_Big_Integer (0) then
+ raise Constraint_Error with "divide by zero";
+ end if;
+
+ Result.Num := Num;
+ Result.Den := Den;
+ Normalize (Result);
+ return Result;
+ end "/";
+
+ ---------------
+ -- Numerator --
+ ---------------
+
+ function Numerator (Arg : Big_Real) return Big_Integer is (Arg.Num);
+
+ -----------------
+ -- Denominator --
+ -----------------
+
+ function Denominator (Arg : Big_Real) return Big_Positive is (Arg.Den);
+
+ ---------
+ -- "=" --
+ ---------
+
+ function "=" (L, R : Big_Real) return Boolean is
+ (abs L.Num = abs R.Num and then L.Den = R.Den);
+
+ ---------
+ -- "<" --
+ ---------
+
+ function "<" (L, R : Big_Real) return Boolean is
+ (abs L.Num * R.Den < abs R.Num * L.Den);
+
+ ----------
+ -- "<=" --
+ ----------
+
+ function "<=" (L, R : Big_Real) return Boolean is (not (R < L));
+
+ ---------
+ -- ">" --
+ ---------
+
+ function ">" (L, R : Big_Real) return Boolean is (R < L);
+
+ ----------
+ -- ">=" --
+ ----------
+
+ function ">=" (L, R : Big_Real) return Boolean is (not (L < R));
+
+ -----------------------
+ -- Float_Conversions --
+ -----------------------
+
+ package body Float_Conversions is
+
+ -----------------
+ -- To_Big_Real --
+ -----------------
+
+ function To_Big_Real (Arg : Num) return Big_Real is
+ begin
+ return From_String (Arg'Image);
+ end To_Big_Real;
+
+ -------------------
+ -- From_Big_Real --
+ -------------------
+
+ function From_Big_Real (Arg : Big_Real) return Num is
+ begin
+ return Num'Value (To_String (Arg));
+ end From_Big_Real;
+
+ end Float_Conversions;
+
+ -----------------------
+ -- Fixed_Conversions --
+ -----------------------
+
+ package body Fixed_Conversions is
+
+ -----------------
+ -- To_Big_Real --
+ -----------------
+
+ function To_Big_Real (Arg : Num) return Big_Real is
+ begin
+ return From_String (Arg'Image);
+ end To_Big_Real;
+
+ -------------------
+ -- From_Big_Real --
+ -------------------
+
+ function From_Big_Real (Arg : Big_Real) return Num is
+ begin
+ return Num'Value (To_String (Arg));
+ end From_Big_Real;
+
+ end Fixed_Conversions;
+
+ ---------------
+ -- To_String --
+ ---------------
+
+ function To_String
+ (Arg : Big_Real; Fore : Field := 2; Aft : Field := 3; Exp : Field := 0)
+ return String
+ is
+ Zero : constant Big_Integer := To_Big_Integer (0);
+ Ten : constant Big_Integer := To_Big_Integer (10);
+
+ function Leading_Padding
+ (Str : String;
+ Min_Length : Field;
+ Char : Character := ' ') return String;
+ -- Return padding of Char concatenated with Str so that the resulting
+ -- string is at least Min_Length long.
+
+ function Trailing_Padding
+ (Str : String;
+ Length : Field;
+ Char : Character := '0') return String;
+ -- Return Str with trailing Char removed, and if needed either
+ -- truncated or concatenated with padding of Char so that the resulting
+ -- string is Length long.
+
+ function Image (N : Natural) return String;
+ -- Return image of N, with no leading space.
+
+ function Numerator_Image
+ (Num : Big_Integer;
+ After : Natural) return String;
+ -- Return image of Num as a float value with After digits after the "."
+ -- and taking Fore, Aft, Exp into account.
+
+ -----------
+ -- Image --
+ -----------
+
+ function Image (N : Natural) return String is
+ S : constant String := Natural'Image (N);
+ begin
+ return S (2 .. S'Last);
+ end Image;
+
+ ---------------------
+ -- Leading_Padding --
+ ---------------------
+
+ function Leading_Padding
+ (Str : String;
+ Min_Length : Field;
+ Char : Character := ' ') return String is
+ begin
+ if Str = "" then
+ return Leading_Padding ("0", Min_Length, Char);
+ else
+ return (1 .. Integer'Max (Integer (Min_Length) - Str'Length, 0)
+ => Char) & Str;
+ end if;
+ end Leading_Padding;
+
+ ----------------------
+ -- Trailing_Padding --
+ ----------------------
+
+ function Trailing_Padding
+ (Str : String;
+ Length : Field;
+ Char : Character := '0') return String is
+ begin
+ if Str'Length > 0 and then Str (Str'Last) = Char then
+ for J in reverse Str'Range loop
+ if Str (J) /= '0' then
+ return Trailing_Padding
+ (Str (Str'First .. J), Length, Char);
+ end if;
+ end loop;
+ end if;
+
+ if Str'Length >= Length then
+ return Str (Str'First .. Str'First + Length - 1);
+ else
+ return Str &
+ (1 .. Integer'Max (Integer (Length) - Str'Length, 0)
+ => Char);
+ end if;
+ end Trailing_Padding;
+
+ ---------------------
+ -- Numerator_Image --
+ ---------------------
+
+ function Numerator_Image
+ (Num : Big_Integer;
+ After : Natural) return String
+ is
+ Tmp : constant String := To_String (Num);
+ Str : constant String (1 .. Tmp'Last - 1) := Tmp (2 .. Tmp'Last);
+ Index : Integer;
+
+ begin
+ if After = 0 then
+ return Leading_Padding (Str, Fore) & "."
+ & Trailing_Padding ("0", Aft);
+ else
+ Index := Str'Last - After;
+
+ if Index < 0 then
+ return Leading_Padding ("0", Fore)
+ & "."
+ & Trailing_Padding ((1 .. -Index => '0') & Str, Aft)
+ & (if Exp = 0 then "" else "E+" & Image (Natural (Exp)));
+ else
+ return Leading_Padding (Str (Str'First .. Index), Fore)
+ & "."
+ & Trailing_Padding (Str (Index + 1 .. Str'Last), Aft)
+ & (if Exp = 0 then "" else "E+" & Image (Natural (Exp)));
+ end if;
+ end if;
+ end Numerator_Image;
+
+ begin
+ if Arg.Num < Zero then
+ declare
+ Str : String := To_String (-Arg, Fore, Aft, Exp);
+ begin
+ if Str (1) = ' ' then
+ for J in 1 .. Str'Last - 1 loop
+ if Str (J + 1) /= ' ' then
+ Str (J) := '-';
+ exit;
+ end if;
+ end loop;
+
+ return Str;
+ else
+ return '-' & Str;
+ end if;
+ end;
+ else
+ -- Compute Num * 10^Aft so that we get Aft significant digits
+ -- in the integer part (rounded) to display.
+
+ return Numerator_Image
+ ((Arg.Num * Ten ** Aft) / Arg.Den, After => Exp + Aft);
+ end if;
+ end To_String;
+
+ -----------------
+ -- From_String --
+ -----------------
+
+ function From_String (Arg : String) return Big_Real is
+ Ten : constant Big_Integer := To_Big_Integer (10);
+ Frac : Optional_Big_Integer;
+ Exp : Integer := 0;
+ Pow : Natural := 0;
+ Index : Natural := 0;
+ Last : Natural := Arg'Last;
+
+ begin
+ for J in reverse Arg'Range loop
+ if Arg (J) in 'e' | 'E' then
+ if Last /= Arg'Last then
+ raise Constraint_Error with "multiple exponents specified";
+ end if;
+
+ Last := J - 1;
+ Exp := Integer'Value (Arg (J + 1 .. Arg'Last));
+ Pow := 0;
+
+ elsif Arg (J) = '.' then
+ Index := J - 1;
+ exit;
+ else
+ Pow := Pow + 1;
+ end if;
+ end loop;
+
+ if Index = 0 then
+ raise Constraint_Error with "invalid real value";
+ end if;
+
+ declare
+ Result : Optional_Big_Real;
+ begin
+ Result.Den := Ten ** Pow;
+ Result.Num := From_String (Arg (Arg'First .. Index)) * Result.Den;
+ Frac := From_String (Arg (Index + 2 .. Last));
+
+ if Result.Num < To_Big_Integer (0) then
+ Result.Num := Result.Num - Frac;
+ else
+ Result.Num := Result.Num + Frac;
+ end if;
+
+ if Exp > 0 then
+ Result.Num := Result.Num * Ten ** Exp;
+ elsif Exp < 0 then
+ Result.Den := Result.Den * Ten ** (-Exp);
+ end if;
+
+ Normalize (Result);
+ return Result;
+ end;
+ end From_String;
+
+ --------------------------
+ -- From_Quotient_String --
+ --------------------------
+
+ function From_Quotient_String (Arg : String) return Big_Real is
+ Index : Natural := 0;
+ begin
+ for J in Arg'First + 1 .. Arg'Last - 1 loop
+ if Arg (J) = '/' then
+ Index := J;
+ exit;
+ end if;
+ end loop;
+
+ if Index = 0 then
+ raise Constraint_Error with "no quotient found";
+ end if;
+
+ return Big_Integers.From_String (Arg (Arg'First .. Index - 1)) /
+ Big_Integers.From_String (Arg (Index + 1 .. Arg'Last));
+ end From_Quotient_String;
+
+ ---------------
+ -- Put_Image --
+ ---------------
+
+ procedure Put_Image
+ (Stream : not null access Ada.Streams.Root_Stream_Type'Class;
+ Arg : Big_Real) is
+ begin
+ Wide_Wide_String'Write (Stream, To_Wide_Wide_String (To_String (Arg)));
+ end Put_Image;
+
+ ---------
+ -- "+" --
+ ---------
+
+ function "+" (L : Big_Real) return Big_Real is
+ Result : Optional_Big_Real;
+ begin
+ Result.Num := L.Num;
+ Result.Den := L.Den;
+ return Result;
+ end "+";
+
+ ---------
+ -- "-" --
+ ---------
+
+ function "-" (L : Big_Real) return Big_Real is
+ (Num => -L.Num, Den => L.Den);
+
+ -----------
+ -- "abs" --
+ -----------
+
+ function "abs" (L : Big_Real) return Big_Real is
+ (Num => abs L.Num, Den => L.Den);
+
+ ---------
+ -- "+" --
+ ---------
+
+ function "+" (L, R : Big_Real) return Big_Real is
+ Result : Optional_Big_Real;
+ begin
+ Result.Num := L.Num * R.Den + R.Num * L.Den;
+ Result.Den := L.Den * R.Den;
+ Normalize (Result);
+ return Result;
+ end "+";
+
+ ---------
+ -- "-" --
+ ---------
+
+ function "-" (L, R : Big_Real) return Big_Real is
+ Result : Optional_Big_Real;
+ begin
+ Result.Num := L.Num * R.Den - R.Num * L.Den;
+ Result.Den := L.Den * R.Den;
+ Normalize (Result);
+ return Result;
+ end "-";
+
+ ---------
+ -- "*" --
+ ---------
+
+ function "*" (L, R : Big_Real) return Big_Real is
+ Result : Optional_Big_Real;
+ begin
+ Result.Num := L.Num * R.Num;
+ Result.Den := L.Den * R.Den;
+ Normalize (Result);
+ return Result;
+ end "*";
+
+ ---------
+ -- "/" --
+ ---------
+
+ function "/" (L, R : Big_Real) return Big_Real is
+ Result : Optional_Big_Real;
+ begin
+ Result.Num := L.Num * R.Den;
+ Result.Den := L.Den * R.Num;
+ Normalize (Result);
+ return Result;
+ end "/";
+
+ ----------
+ -- "**" --
+ ----------
+
+ function "**" (L : Big_Real; R : Integer) return Big_Real is
+ Result : Optional_Big_Real;
+ begin
+ if R = 0 then
+ Result.Num := To_Big_Integer (1);
+ Result.Den := To_Big_Integer (1);
+ else
+ if R < 0 then
+ Result.Num := L.Den ** (-R);
+ Result.Den := L.Num ** (-R);
+ else
+ Result.Num := L.Num ** R;
+ Result.Den := L.Den ** R;
+ end if;
+
+ Normalize (Result);
+ end if;
+
+ return Result;
+ end "**";
+
+ ---------
+ -- Min --
+ ---------
+
+ function Min (L, R : Big_Real) return Big_Real is (if L < R then L else R);
+
+ ---------
+ -- Max --
+ ---------
+
+ function Max (L, R : Big_Real) return Big_Real is (if L > R then L else R);
+
+ ---------------
+ -- Normalize --
+ ---------------
+
+ procedure Normalize (Arg : in out Big_Real) is
+ begin
+ if Arg.Den < To_Big_Integer (0) then
+ Arg.Num := -Arg.Num;
+ Arg.Den := -Arg.Den;
+ end if;
+
+ declare
+ GCD : constant Big_Integer :=
+ Greatest_Common_Divisor (Arg.Num, Arg.Den);
+ begin
+ Arg.Num := Arg.Num / GCD;
+ Arg.Den := Arg.Den / GCD;
+ end;
+ end Normalize;
+
+end Ada.Numerics.Big_Numbers.Big_Reals;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- GNAT RUN-TIME COMPONENTS --
+-- --
+-- ADA.NUMERICS.BIG_NUMBERS.BIG_REALS --
+-- --
+-- S p e c --
+-- --
+-- This specification is derived from the Ada Reference Manual for use with --
+-- GNAT. In accordance with the copyright of that document, you can freely --
+-- copy and modify this specification, provided that if you redistribute a --
+-- modified version, any changes that you have made are clearly indicated. --
+-- --
+------------------------------------------------------------------------------
+
+with Ada.Numerics.Big_Numbers.Big_Integers;
+with Ada.Streams;
+
+-- Note that some Ada 2020 aspects are commented out since they are not
+-- supported yet.
+
+package Ada.Numerics.Big_Numbers.Big_Reals
+ with Preelaborate
+-- Nonblocking, Global => in out synchronized Big_Reals
+is
+ type Optional_Big_Real is private with
+ Default_Initial_Condition => not Is_Valid (Optional_Big_Real);
+-- Real_Literal => From_String,
+-- Put_Image => Put_Image;
+
+ function Is_Valid (Arg : Optional_Big_Real) return Boolean;
+
+ function No_Big_Real return Optional_Big_Real
+ with Post => not Is_Valid (No_Big_Real'Result);
+
+ subtype Big_Real is Optional_Big_Real
+ with Dynamic_Predicate => Is_Valid (Big_Real),
+ Predicate_Failure => (raise Constraint_Error);
+
+ function "/" (Num, Den : Big_Integers.Big_Integer) return Big_Real;
+-- with Pre => (if Big_Integers."=" (Den, Big_Integers.To_Big_Integer (0))
+-- then raise Constraint_Error);
+
+ function Numerator (Arg : Big_Real) return Big_Integers.Big_Integer;
+
+ function Denominator (Arg : Big_Real) return Big_Integers.Big_Positive
+ with Post =>
+ (Arg = To_Real (0)) or else
+ (Big_Integers."="
+ (Big_Integers.Greatest_Common_Divisor
+ (Numerator (Arg), Denominator'Result),
+ Big_Integers.To_Big_Integer (1)));
+
+ function To_Big_Real
+ (Arg : Big_Integers.Big_Integer)
+ return Big_Real is (Arg / Big_Integers.To_Big_Integer (1));
+
+ function To_Real (Arg : Integer) return Big_Real is
+ (Big_Integers.To_Big_Integer (Arg) / Big_Integers.To_Big_Integer (1));
+
+ function "=" (L, R : Big_Real) return Boolean;
+
+ function "<" (L, R : Big_Real) return Boolean;
+
+ function "<=" (L, R : Big_Real) return Boolean;
+
+ function ">" (L, R : Big_Real) return Boolean;
+
+ function ">=" (L, R : Big_Real) return Boolean;
+
+ function In_Range (Arg, Low, High : Big_Real) return Boolean is
+ (Low <= Arg and then Arg <= High);
+
+ generic
+ type Num is digits <>;
+ package Float_Conversions is
+
+ function To_Big_Real (Arg : Num) return Big_Real;
+
+ function From_Big_Real (Arg : Big_Real) return Num
+ with Pre => In_Range (Arg,
+ Low => To_Big_Real (Num'First),
+ High => To_Big_Real (Num'Last))
+ or else (raise Constraint_Error);
+
+ end Float_Conversions;
+
+ generic
+ type Num is delta <>;
+ package Fixed_Conversions is
+
+ function To_Big_Real (Arg : Num) return Big_Real;
+
+ function From_Big_Real (Arg : Big_Real) return Num
+ with Pre => In_Range (Arg,
+ Low => To_Big_Real (Num'First),
+ High => To_Big_Real (Num'Last))
+ or else (raise Constraint_Error);
+
+ end Fixed_Conversions;
+
+ function To_String (Arg : Big_Real;
+ Fore : Field := 2;
+ Aft : Field := 3;
+ Exp : Field := 0) return String
+ with Post => To_String'Result'First = 1;
+
+ function From_String (Arg : String) return Big_Real;
+
+ function To_Quotient_String (Arg : Big_Real) return String is
+ (Big_Integers.To_String (Numerator (Arg)) & " / "
+ & Big_Integers.To_String (Denominator (Arg)));
+
+ function From_Quotient_String (Arg : String) return Big_Real;
+
+ procedure Put_Image
+ (Stream : not null access Ada.Streams.Root_Stream_Type'Class;
+ Arg : Big_Real);
+
+ function "+" (L : Big_Real) return Big_Real;
+
+ function "-" (L : Big_Real) return Big_Real;
+
+ function "abs" (L : Big_Real) return Big_Real;
+
+ function "+" (L, R : Big_Real) return Big_Real;
+
+ function "-" (L, R : Big_Real) return Big_Real;
+
+ function "*" (L, R : Big_Real) return Big_Real;
+
+ function "/" (L, R : Big_Real) return Big_Real;
+
+ function "**" (L : Big_Real; R : Integer) return Big_Real;
+
+ function Min (L, R : Big_Real) return Big_Real;
+
+ function Max (L, R : Big_Real) return Big_Real;
+
+private
+
+ type Optional_Big_Real is record
+ Num, Den : Big_Integers.Optional_Big_Integer;
+ end record;
+
+end Ada.Numerics.Big_Numbers.Big_Reals;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- GNAT RUN-TIME COMPONENTS --
+-- --
+-- A D A . N U M E R I C S --
+-- --
+-- S p e c --
+-- --
+-- This specification is derived from the Ada Reference Manual for use with --
+-- GNAT. In accordance with the copyright of that document, you can freely --
+-- copy and modify this specification, provided that if you redistribute a --
+-- modified version, any changes that you have made are clearly indicated. --
+-- --
+------------------------------------------------------------------------------
+
+-- Note that some Ada 2020 aspects are commented out since they are not
+-- supported yet.
+
+package Ada.Numerics.Big_Numbers
+ -- with Pure, Nonblocking, Global => null
+ with Pure
+is
+ subtype Field is Integer range 0 .. 255;
+ subtype Number_Base is Integer range 2 .. 16;
+end Ada.Numerics.Big_Numbers;
-- --
------------------------------------------------------------------------------
--- This package provides arbitrary precision signed integer arithmetic for
--- use in computing intermediate values in expressions for the case where
--- pragma Overflow_Check (Eliminate) is in effect.
-
-with System; use System;
-with System.Secondary_Stack; use System.Secondary_Stack;
-with System.Storage_Elements; use System.Storage_Elements;
+with System.Generic_Bignums;
+with Ada.Unchecked_Conversion;
package body System.Bignums is
- use Interfaces;
- -- So that operations on Unsigned_32 are available
-
- type DD is mod Base ** 2;
- -- Double length digit used for intermediate computations
-
- function MSD (X : DD) return SD is (SD (X / Base));
- function LSD (X : DD) return SD is (SD (X mod Base));
- -- Most significant and least significant digit of double digit value
-
- function "&" (X, Y : SD) return DD is (DD (X) * Base + DD (Y));
- -- Compose double digit value from two single digit values
-
- subtype LLI is Long_Long_Integer;
-
- One_Data : constant Digit_Vector (1 .. 1) := (1 => 1);
- -- Constant one
-
- Zero_Data : constant Digit_Vector (1 .. 0) := (1 .. 0 => 0);
- -- Constant zero
-
- -----------------------
- -- Local Subprograms --
- -----------------------
-
- function Add
- (X, Y : Digit_Vector;
- X_Neg : Boolean;
- Y_Neg : Boolean) return Bignum
- with
- Pre => X'First = 1 and then Y'First = 1;
- -- This procedure adds two signed numbers returning the Sum, it is used
- -- for both addition and subtraction. The value computed is X + Y, with
- -- X_Neg and Y_Neg giving the signs of the operands.
-
- function Allocate_Bignum (Len : Length) return Bignum with
- Post => Allocate_Bignum'Result.Len = Len;
- -- Allocate Bignum value of indicated length on secondary stack. On return
- -- the Neg and D fields are left uninitialized.
-
- type Compare_Result is (LT, EQ, GT);
- -- Indicates result of comparison in following call
-
- function Compare
- (X, Y : Digit_Vector;
- X_Neg, Y_Neg : Boolean) return Compare_Result
- with
- Pre => X'First = 1 and then Y'First = 1;
- -- Compare (X with sign X_Neg) with (Y with sign Y_Neg), and return the
- -- result of the signed comparison.
-
- procedure Div_Rem
- (X, Y : Bignum;
- Quotient : out Bignum;
- Remainder : out Bignum;
- Discard_Quotient : Boolean := False;
- Discard_Remainder : Boolean := False);
- -- Returns the Quotient and Remainder from dividing abs (X) by abs (Y). The
- -- values of X and Y are not modified. If Discard_Quotient is True, then
- -- Quotient is undefined on return, and if Discard_Remainder is True, then
- -- Remainder is undefined on return. Service routine for Big_Div/Rem/Mod.
-
- procedure Free_Bignum (X : Bignum) is null;
- -- Called to free a Bignum value used in intermediate computations. In
- -- this implementation using the secondary stack, it does nothing at all,
- -- because we rely on Mark/Release, but it may be of use for some
- -- alternative implementation.
-
- function Normalize
- (X : Digit_Vector;
- Neg : Boolean := False) return Bignum;
- -- Given a digit vector and sign, allocate and construct a Bignum value.
- -- Note that X may have leading zeroes which must be removed, and if the
- -- result is zero, the sign is forced positive.
-
- ---------
- -- Add --
- ---------
-
- function Add
- (X, Y : Digit_Vector;
- X_Neg : Boolean;
- Y_Neg : Boolean) return Bignum
- is
- begin
- -- If signs are the same, we are doing an addition, it is convenient to
- -- ensure that the first operand is the longer of the two.
-
- if X_Neg = Y_Neg then
- if X'Last < Y'Last then
- return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg);
-
- -- Here signs are the same, and the first operand is the longer
-
- else
- pragma Assert (X_Neg = Y_Neg and then X'Last >= Y'Last);
-
- -- Do addition, putting result in Sum (allowing for carry)
-
- declare
- Sum : Digit_Vector (0 .. X'Last);
- RD : DD;
-
- begin
- RD := 0;
- for J in reverse 1 .. X'Last loop
- RD := RD + DD (X (J));
-
- if J >= 1 + (X'Last - Y'Last) then
- RD := RD + DD (Y (J - (X'Last - Y'Last)));
- end if;
-
- Sum (J) := LSD (RD);
- RD := RD / Base;
- end loop;
-
- Sum (0) := SD (RD);
- return Normalize (Sum, X_Neg);
- end;
- end if;
-
- -- Signs are different so really this is a subtraction, we want to make
- -- sure that the largest magnitude operand is the first one, and then
- -- the result will have the sign of the first operand.
-
- else
- declare
- CR : constant Compare_Result := Compare (X, Y, False, False);
-
- begin
- if CR = EQ then
- return Normalize (Zero_Data);
-
- elsif CR = LT then
- return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg);
-
- else
- pragma Assert (X_Neg /= Y_Neg and then CR = GT);
-
- -- Do subtraction, putting result in Diff
-
- declare
- Diff : Digit_Vector (1 .. X'Length);
- RD : DD;
-
- begin
- RD := 0;
- for J in reverse 1 .. X'Last loop
- RD := RD + DD (X (J));
-
- if J >= 1 + (X'Last - Y'Last) then
- RD := RD - DD (Y (J - (X'Last - Y'Last)));
- end if;
-
- Diff (J) := LSD (RD);
- RD := (if RD < Base then 0 else -1);
- end loop;
-
- return Normalize (Diff, X_Neg);
- end;
- end if;
- end;
- end if;
- end Add;
-
- ---------------------
- -- Allocate_Bignum --
- ---------------------
-
- function Allocate_Bignum (Len : Length) return Bignum is
- Addr : Address;
-
- begin
- -- Change the if False here to if True to get allocation on the heap
- -- instead of the secondary stack, which is convenient for debugging
- -- System.Bignum itself.
-
- if False then
- declare
- B : Bignum;
- begin
- B := new Bignum_Data'(Len, False, (others => 0));
- return B;
- end;
-
- -- Normal case of allocation on the secondary stack
-
- else
- -- Note: The approach used here is designed to avoid strict aliasing
- -- warnings that appeared previously using unchecked conversion.
-
- SS_Allocate (Addr, Storage_Offset (4 + 4 * Len));
-
- declare
- B : Bignum;
- for B'Address use Addr'Address;
- pragma Import (Ada, B);
-
- BD : Bignum_Data (Len);
- for BD'Address use Addr;
- pragma Import (Ada, BD);
-
- -- Expose a writable view of discriminant BD.Len so that we can
- -- initialize it. We need to use the exact layout of the record
- -- to ensure that the Length field has 24 bits as expected.
-
- type Bignum_Data_Header is record
- Len : Length;
- Neg : Boolean;
- end record;
-
- for Bignum_Data_Header use record
- Len at 0 range 0 .. 23;
- Neg at 3 range 0 .. 7;
- end record;
-
- BDH : Bignum_Data_Header;
- for BDH'Address use BD'Address;
- pragma Import (Ada, BDH);
-
- pragma Assert (BDH.Len'Size = BD.Len'Size);
-
- begin
- BDH.Len := Len;
- return B;
- end;
- end if;
- end Allocate_Bignum;
-
- -------------
- -- Big_Abs --
- -------------
-
- function Big_Abs (X : Bignum) return Bignum is
- begin
- return Normalize (X.D);
- end Big_Abs;
-
- -------------
- -- Big_Add --
- -------------
-
- function Big_Add (X, Y : Bignum) return Bignum is
- begin
- return Add (X.D, Y.D, X.Neg, Y.Neg);
- end Big_Add;
-
- -------------
- -- Big_Div --
- -------------
-
- -- This table is excerpted from RM 4.5.5(28-30) and shows how the result
- -- varies with the signs of the operands.
-
- -- A B A/B A B A/B
- --
- -- 10 5 2 -10 5 -2
- -- 11 5 2 -11 5 -2
- -- 12 5 2 -12 5 -2
- -- 13 5 2 -13 5 -2
- -- 14 5 2 -14 5 -2
- --
- -- A B A/B A B A/B
- --
- -- 10 -5 -2 -10 -5 2
- -- 11 -5 -2 -11 -5 2
- -- 12 -5 -2 -12 -5 2
- -- 13 -5 -2 -13 -5 2
- -- 14 -5 -2 -14 -5 2
-
- function Big_Div (X, Y : Bignum) return Bignum is
- Q, R : Bignum;
- begin
- Div_Rem (X, Y, Q, R, Discard_Remainder => True);
- Q.Neg := Q.Len > 0 and then (X.Neg xor Y.Neg);
- return Q;
- end Big_Div;
-
- -------------
- -- Big_Exp --
- -------------
-
- function Big_Exp (X, Y : Bignum) return Bignum is
+ package Sec_Stack_Bignums is new
+ System.Generic_Bignums (Use_Secondary_Stack => True);
+ use Sec_Stack_Bignums;
- function "**" (X : Bignum; Y : SD) return Bignum;
- -- Internal routine where we know right operand is one word
+ function "+" is new Ada.Unchecked_Conversion
+ (Bignum, Sec_Stack_Bignums.Bignum);
- ----------
- -- "**" --
- ----------
+ function "-" is new Ada.Unchecked_Conversion
+ (Sec_Stack_Bignums.Bignum, Bignum);
- function "**" (X : Bignum; Y : SD) return Bignum is
- begin
- case Y is
+ function Big_Add (X, Y : Bignum) return Bignum is
+ (-Sec_Stack_Bignums.Big_Add (+X, +Y));
- -- X ** 0 is 1
-
- when 0 =>
- return Normalize (One_Data);
-
- -- X ** 1 is X
-
- when 1 =>
- return Normalize (X.D);
-
- -- X ** 2 is X * X
-
- when 2 =>
- return Big_Mul (X, X);
-
- -- For X greater than 2, use the recursion
-
- -- X even, X ** Y = (X ** (Y/2)) ** 2;
- -- X odd, X ** Y = (X ** (Y/2)) ** 2 * X;
-
- when others =>
- declare
- XY2 : constant Bignum := X ** (Y / 2);
- XY2S : constant Bignum := Big_Mul (XY2, XY2);
- Res : Bignum;
-
- begin
- Free_Bignum (XY2);
-
- -- Raise storage error if intermediate value is getting too
- -- large, which we arbitrarily define as 200 words for now.
-
- if XY2S.Len > 200 then
- Free_Bignum (XY2S);
- raise Storage_Error with
- "exponentiation result is too large";
- end if;
-
- -- Otherwise take care of even/odd cases
-
- if (Y and 1) = 0 then
- return XY2S;
-
- else
- Res := Big_Mul (XY2S, X);
- Free_Bignum (XY2S);
- return Res;
- end if;
- end;
- end case;
- end "**";
-
- -- Start of processing for Big_Exp
-
- begin
- -- Error if right operand negative
-
- if Y.Neg then
- raise Constraint_Error with "exponentiation to negative power";
-
- -- X ** 0 is always 1 (including 0 ** 0, so do this test first)
-
- elsif Y.Len = 0 then
- return Normalize (One_Data);
-
- -- 0 ** X is always 0 (for X non-zero)
-
- elsif X.Len = 0 then
- return Normalize (Zero_Data);
-
- -- (+1) ** Y = 1
- -- (-1) ** Y = +/-1 depending on whether Y is even or odd
-
- elsif X.Len = 1 and then X.D (1) = 1 then
- return Normalize
- (X.D, Neg => X.Neg and then ((Y.D (Y.Len) and 1) = 1));
-
- -- If the absolute value of the base is greater than 1, then the
- -- exponent must not be bigger than one word, otherwise the result
- -- is ludicrously large, and we just signal Storage_Error right away.
-
- elsif Y.Len > 1 then
- raise Storage_Error with "exponentiation result is too large";
-
- -- Special case (+/-)2 ** K, where K is 1 .. 31 using a shift
-
- elsif X.Len = 1 and then X.D (1) = 2 and then Y.D (1) < 32 then
- declare
- D : constant Digit_Vector (1 .. 1) :=
- (1 => Shift_Left (SD'(1), Natural (Y.D (1))));
- begin
- return Normalize (D, X.Neg);
- end;
-
- -- Remaining cases have right operand of one word
-
- else
- return X ** Y.D (1);
- end if;
- end Big_Exp;
-
- ------------
- -- Big_EQ --
- ------------
-
- function Big_EQ (X, Y : Bignum) return Boolean is
- begin
- return Compare (X.D, Y.D, X.Neg, Y.Neg) = EQ;
- end Big_EQ;
-
- ------------
- -- Big_GE --
- ------------
-
- function Big_GE (X, Y : Bignum) return Boolean is
- begin
- return Compare (X.D, Y.D, X.Neg, Y.Neg) /= LT;
- end Big_GE;
-
- ------------
- -- Big_GT --
- ------------
-
- function Big_GT (X, Y : Bignum) return Boolean is
- begin
- return Compare (X.D, Y.D, X.Neg, Y.Neg) = GT;
- end Big_GT;
-
- ------------
- -- Big_LE --
- ------------
-
- function Big_LE (X, Y : Bignum) return Boolean is
- begin
- return Compare (X.D, Y.D, X.Neg, Y.Neg) /= GT;
- end Big_LE;
-
- ------------
- -- Big_LT --
- ------------
-
- function Big_LT (X, Y : Bignum) return Boolean is
- begin
- return Compare (X.D, Y.D, X.Neg, Y.Neg) = LT;
- end Big_LT;
-
- -------------
- -- Big_Mod --
- -------------
-
- -- This table is excerpted from RM 4.5.5(28-30) and shows how the result
- -- of Rem and Mod vary with the signs of the operands.
-
- -- A B A mod B A rem B A B A mod B A rem B
-
- -- 10 5 0 0 -10 5 0 0
- -- 11 5 1 1 -11 5 4 -1
- -- 12 5 2 2 -12 5 3 -2
- -- 13 5 3 3 -13 5 2 -3
- -- 14 5 4 4 -14 5 1 -4
-
- -- A B A mod B A rem B A B A mod B A rem B
-
- -- 10 -5 0 0 -10 -5 0 0
- -- 11 -5 -4 1 -11 -5 -1 -1
- -- 12 -5 -3 2 -12 -5 -2 -2
- -- 13 -5 -2 3 -13 -5 -3 -3
- -- 14 -5 -1 4 -14 -5 -4 -4
-
- function Big_Mod (X, Y : Bignum) return Bignum is
- Q, R : Bignum;
-
- begin
- -- If signs are same, result is same as Rem
-
- if X.Neg = Y.Neg then
- return Big_Rem (X, Y);
-
- -- Case where Mod is different
-
- else
- -- Do division
-
- Div_Rem (X, Y, Q, R, Discard_Quotient => True);
-
- -- Zero result is unchanged
-
- if R.Len = 0 then
- return R;
-
- -- Otherwise adjust result
-
- else
- declare
- T1 : constant Bignum := Big_Sub (Y, R);
- begin
- T1.Neg := Y.Neg;
- Free_Bignum (R);
- return T1;
- end;
- end if;
- end if;
- end Big_Mod;
-
- -------------
- -- Big_Mul --
- -------------
+ function Big_Sub (X, Y : Bignum) return Bignum is
+ (-Sec_Stack_Bignums.Big_Sub (+X, +Y));
function Big_Mul (X, Y : Bignum) return Bignum is
- Result : Digit_Vector (1 .. X.Len + Y.Len) := (others => 0);
- -- Accumulate result (max length of result is sum of operand lengths)
-
- L : Length;
- -- Current result digit
-
- D : DD;
- -- Result digit
-
- begin
- for J in 1 .. X.Len loop
- for K in 1 .. Y.Len loop
- L := Result'Last - (X.Len - J) - (Y.Len - K);
- D := DD (X.D (J)) * DD (Y.D (K)) + DD (Result (L));
- Result (L) := LSD (D);
- D := D / Base;
-
- -- D is carry which must be propagated
+ (-Sec_Stack_Bignums.Big_Mul (+X, +Y));
- while D /= 0 and then L >= 1 loop
- L := L - 1;
- D := D + DD (Result (L));
- Result (L) := LSD (D);
- D := D / Base;
- end loop;
+ function Big_Div (X, Y : Bignum) return Bignum is
+ (-Sec_Stack_Bignums.Big_Div (+X, +Y));
- -- Must not have a carry trying to extend max length
+ function Big_Exp (X, Y : Bignum) return Bignum is
+ (-Sec_Stack_Bignums.Big_Exp (+X, +Y));
- pragma Assert (D = 0);
- end loop;
- end loop;
-
- -- Return result
-
- return Normalize (Result, X.Neg xor Y.Neg);
- end Big_Mul;
-
- ------------
- -- Big_NE --
- ------------
-
- function Big_NE (X, Y : Bignum) return Boolean is
- begin
- return Compare (X.D, Y.D, X.Neg, Y.Neg) /= EQ;
- end Big_NE;
-
- -------------
- -- Big_Neg --
- -------------
-
- function Big_Neg (X : Bignum) return Bignum is
- begin
- return Normalize (X.D, not X.Neg);
- end Big_Neg;
-
- -------------
- -- Big_Rem --
- -------------
-
- -- This table is excerpted from RM 4.5.5(28-30) and shows how the result
- -- varies with the signs of the operands.
-
- -- A B A rem B A B A rem B
-
- -- 10 5 0 -10 5 0
- -- 11 5 1 -11 5 -1
- -- 12 5 2 -12 5 -2
- -- 13 5 3 -13 5 -3
- -- 14 5 4 -14 5 -4
-
- -- A B A rem B A B A rem B
-
- -- 10 -5 0 -10 -5 0
- -- 11 -5 1 -11 -5 -1
- -- 12 -5 2 -12 -5 -2
- -- 13 -5 3 -13 -5 -3
- -- 14 -5 4 -14 -5 -4
+ function Big_Mod (X, Y : Bignum) return Bignum is
+ (-Sec_Stack_Bignums.Big_Mod (+X, +Y));
function Big_Rem (X, Y : Bignum) return Bignum is
- Q, R : Bignum;
- begin
- Div_Rem (X, Y, Q, R, Discard_Quotient => True);
- R.Neg := R.Len > 0 and then X.Neg;
- return R;
- end Big_Rem;
-
- -------------
- -- Big_Sub --
- -------------
-
- function Big_Sub (X, Y : Bignum) return Bignum is
- begin
- -- If right operand zero, return left operand (avoiding sharing)
-
- if Y.Len = 0 then
- return Normalize (X.D, X.Neg);
-
- -- Otherwise add negative of right operand
-
- else
- return Add (X.D, Y.D, X.Neg, not Y.Neg);
- end if;
- end Big_Sub;
-
- -------------
- -- Compare --
- -------------
-
- function Compare
- (X, Y : Digit_Vector;
- X_Neg, Y_Neg : Boolean) return Compare_Result
- is
- begin
- -- Signs are different, that's decisive, since 0 is always plus
-
- if X_Neg /= Y_Neg then
- return (if X_Neg then LT else GT);
-
- -- Lengths are different, that's decisive since no leading zeroes
-
- elsif X'Last /= Y'Last then
- return (if (X'Last > Y'Last) xor X_Neg then GT else LT);
-
- -- Need to compare data
-
- else
- for J in X'Range loop
- if X (J) /= Y (J) then
- return (if (X (J) > Y (J)) xor X_Neg then GT else LT);
- end if;
- end loop;
-
- return EQ;
- end if;
- end Compare;
-
- -------------
- -- Div_Rem --
- -------------
-
- procedure Div_Rem
- (X, Y : Bignum;
- Quotient : out Bignum;
- Remainder : out Bignum;
- Discard_Quotient : Boolean := False;
- Discard_Remainder : Boolean := False)
- is
- begin
- -- Error if division by zero
-
- if Y.Len = 0 then
- raise Constraint_Error with "division by zero";
- end if;
-
- -- Handle simple cases with special tests
-
- -- If X < Y then quotient is zero and remainder is X
-
- if Compare (X.D, Y.D, False, False) = LT then
- Remainder := Normalize (X.D);
- Quotient := Normalize (Zero_Data);
- return;
-
- -- If both X and Y are less than 2**63-1, we can use Long_Long_Integer
- -- arithmetic. Note it is good not to do an accurate range check against
- -- Long_Long_Integer since -2**63 / -1 overflows.
-
- elsif (X.Len <= 1 or else (X.Len = 2 and then X.D (1) < 2**31))
- and then
- (Y.Len <= 1 or else (Y.Len = 2 and then Y.D (1) < 2**31))
- then
- declare
- A : constant LLI := abs (From_Bignum (X));
- B : constant LLI := abs (From_Bignum (Y));
- begin
- Quotient := To_Bignum (A / B);
- Remainder := To_Bignum (A rem B);
- return;
- end;
-
- -- Easy case if divisor is one digit
-
- elsif Y.Len = 1 then
- declare
- ND : DD;
- Div : constant DD := DD (Y.D (1));
-
- Result : Digit_Vector (1 .. X.Len);
- Remdr : Digit_Vector (1 .. 1);
-
- begin
- ND := 0;
- for J in 1 .. X.Len loop
- ND := Base * ND + DD (X.D (J));
- Result (J) := SD (ND / Div);
- ND := ND rem Div;
- end loop;
-
- Quotient := Normalize (Result);
- Remdr (1) := SD (ND);
- Remainder := Normalize (Remdr);
- return;
- end;
- end if;
-
- -- The complex full multi-precision case. We will employ algorithm
- -- D defined in the section "The Classical Algorithms" (sec. 4.3.1)
- -- of Donald Knuth's "The Art of Computer Programming", Vol. 2, 2nd
- -- edition. The terminology is adjusted for this section to match that
- -- reference.
-
- -- We are dividing X.Len digits of X (called u here) by Y.Len digits
- -- of Y (called v here), developing the quotient and remainder. The
- -- numbers are represented using Base, which was chosen so that we have
- -- the operations of multiplying to single digits (SD) to form a double
- -- digit (DD), and dividing a double digit (DD) by a single digit (SD)
- -- to give a single digit quotient and a single digit remainder.
-
- -- Algorithm D from Knuth
-
- -- Comments here with square brackets are directly from Knuth
-
- Algorithm_D : declare
-
- -- The following lower case variables correspond exactly to the
- -- terminology used in algorithm D.
-
- m : constant Length := X.Len - Y.Len;
- n : constant Length := Y.Len;
- b : constant DD := Base;
-
- u : Digit_Vector (0 .. m + n);
- v : Digit_Vector (1 .. n);
- q : Digit_Vector (0 .. m);
- r : Digit_Vector (1 .. n);
-
- u0 : SD renames u (0);
- v1 : SD renames v (1);
- v2 : SD renames v (2);
-
- d : DD;
- j : Length;
- qhat : DD;
- rhat : DD;
- temp : DD;
-
- begin
- -- Initialize data of left and right operands
-
- for J in 1 .. m + n loop
- u (J) := X.D (J);
- end loop;
-
- for J in 1 .. n loop
- v (J) := Y.D (J);
- end loop;
-
- -- [Division of nonnegative integers.] Given nonnegative integers u
- -- = (ul,u2..um+n) and v = (v1,v2..vn), where v1 /= 0 and n > 1, we
- -- form the quotient u / v = (q0,ql..qm) and the remainder u mod v =
- -- (r1,r2..rn).
-
- pragma Assert (v1 /= 0);
- pragma Assert (n > 1);
-
- -- Dl. [Normalize.] Set d = b/(vl + 1). Then set (u0,u1,u2..um+n)
- -- equal to (u1,u2..um+n) times d, and set (v1,v2..vn) equal to
- -- (v1,v2..vn) times d. Note the introduction of a new digit position
- -- u0 at the left of u1; if d = 1 all we need to do in this step is
- -- to set u0 = 0.
-
- d := b / (DD (v1) + 1);
-
- if d = 1 then
- u0 := 0;
-
- else
- declare
- Carry : DD;
- Tmp : DD;
-
- begin
- -- Multiply Dividend (u) by d
-
- Carry := 0;
- for J in reverse 1 .. m + n loop
- Tmp := DD (u (J)) * d + Carry;
- u (J) := LSD (Tmp);
- Carry := Tmp / Base;
- end loop;
-
- u0 := SD (Carry);
-
- -- Multiply Divisor (v) by d
-
- Carry := 0;
- for J in reverse 1 .. n loop
- Tmp := DD (v (J)) * d + Carry;
- v (J) := LSD (Tmp);
- Carry := Tmp / Base;
- end loop;
-
- pragma Assert (Carry = 0);
- end;
- end if;
-
- -- D2. [Initialize j.] Set j = 0. The loop on j, steps D2 through D7,
- -- will be essentially a division of (uj, uj+1..uj+n) by (v1,v2..vn)
- -- to get a single quotient digit qj.
-
- j := 0;
-
- -- Loop through digits
-
- loop
- -- Note: In the original printing, step D3 was as follows:
-
- -- D3. [Calculate qhat.] If uj = v1, set qhat to b-l; otherwise
- -- set qhat to (uj,uj+1)/v1. Now test if v2 * qhat is greater than
- -- (uj*b + uj+1 - qhat*v1)*b + uj+2. If so, decrease qhat by 1 and
- -- repeat this test
-
- -- This had a bug not discovered till 1995, see Vol 2 errata:
- -- http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz. Under
- -- rare circumstances the expression in the test could overflow.
- -- This version was further corrected in 2005, see Vol 2 errata:
- -- http://www-cs-faculty.stanford.edu/~uno/all2-pre.ps.gz.
- -- The code below is the fixed version of this step.
-
- -- D3. [Calculate qhat.] Set qhat to (uj,uj+1)/v1 and rhat to
- -- to (uj,uj+1) mod v1.
-
- temp := u (j) & u (j + 1);
- qhat := temp / DD (v1);
- rhat := temp mod DD (v1);
-
- -- D3 (continued). Now test if qhat >= b or v2*qhat > (rhat,uj+2):
- -- if so, decrease qhat by 1, increase rhat by v1, and repeat this
- -- test if rhat < b. [The test on v2 determines at high speed
- -- most of the cases in which the trial value qhat is one too
- -- large, and eliminates all cases where qhat is two too large.]
-
- while qhat >= b
- or else DD (v2) * qhat > LSD (rhat) & u (j + 2)
- loop
- qhat := qhat - 1;
- rhat := rhat + DD (v1);
- exit when rhat >= b;
- end loop;
-
- -- D4. [Multiply and subtract.] Replace (uj,uj+1..uj+n) by
- -- (uj,uj+1..uj+n) minus qhat times (v1,v2..vn). This step
- -- consists of a simple multiplication by a one-place number,
- -- combined with a subtraction.
-
- -- The digits (uj,uj+1..uj+n) are always kept positive; if the
- -- result of this step is actually negative then (uj,uj+1..uj+n)
- -- is left as the true value plus b**(n+1), i.e. as the b's
- -- complement of the true value, and a "borrow" to the left is
- -- remembered.
-
- declare
- Borrow : SD;
- Carry : DD;
- Temp : DD;
-
- Negative : Boolean;
- -- Records if subtraction causes a negative result, requiring
- -- an add back (case where qhat turned out to be 1 too large).
-
- begin
- Borrow := 0;
- for K in reverse 1 .. n loop
- Temp := qhat * DD (v (K)) + DD (Borrow);
- Borrow := MSD (Temp);
-
- if LSD (Temp) > u (j + K) then
- Borrow := Borrow + 1;
- end if;
-
- u (j + K) := u (j + K) - LSD (Temp);
- end loop;
-
- Negative := u (j) < Borrow;
- u (j) := u (j) - Borrow;
-
- -- D5. [Test remainder.] Set qj = qhat. If the result of step
- -- D4 was negative, we will do the add back step (step D6).
-
- q (j) := LSD (qhat);
-
- if Negative then
-
- -- D6. [Add back.] Decrease qj by 1, and add (0,v1,v2..vn)
- -- to (uj,uj+1,uj+2..uj+n). (A carry will occur to the left
- -- of uj, and it is be ignored since it cancels with the
- -- borrow that occurred in D4.)
-
- q (j) := q (j) - 1;
-
- Carry := 0;
- for K in reverse 1 .. n loop
- Temp := DD (v (K)) + DD (u (j + K)) + Carry;
- u (j + K) := LSD (Temp);
- Carry := Temp / Base;
- end loop;
-
- u (j) := u (j) + SD (Carry);
- end if;
- end;
-
- -- D7. [Loop on j.] Increase j by one. Now if j <= m, go back to
- -- D3 (the start of the loop on j).
-
- j := j + 1;
- exit when not (j <= m);
- end loop;
-
- -- D8. [Unnormalize.] Now (qo,ql..qm) is the desired quotient, and
- -- the desired remainder may be obtained by dividing (um+1..um+n)
- -- by d.
-
- if not Discard_Quotient then
- Quotient := Normalize (q);
- end if;
-
- if not Discard_Remainder then
- declare
- Remdr : DD;
-
- begin
- Remdr := 0;
- for K in 1 .. n loop
- Remdr := Base * Remdr + DD (u (m + K));
- r (K) := SD (Remdr / d);
- Remdr := Remdr rem d;
- end loop;
-
- pragma Assert (Remdr = 0);
- end;
-
- Remainder := Normalize (r);
- end if;
- end Algorithm_D;
- end Div_Rem;
-
- -----------------
- -- From_Bignum --
- -----------------
-
- function From_Bignum (X : Bignum) return Long_Long_Integer is
- begin
- if X.Len = 0 then
- return 0;
-
- elsif X.Len = 1 then
- return (if X.Neg then -LLI (X.D (1)) else LLI (X.D (1)));
-
- elsif X.Len = 2 then
- declare
- Mag : constant DD := X.D (1) & X.D (2);
- begin
- if X.Neg and then Mag <= 2 ** 63 then
- return -LLI (Mag);
- elsif Mag < 2 ** 63 then
- return LLI (Mag);
- end if;
- end;
- end if;
-
- raise Constraint_Error with "expression value out of range";
- end From_Bignum;
-
- -------------------------
- -- Bignum_In_LLI_Range --
- -------------------------
+ (-Sec_Stack_Bignums.Big_Rem (+X, +Y));
+
+ function Big_Neg (X : Bignum) return Bignum is
+ (-Sec_Stack_Bignums.Big_Neg (+X));
+
+ function Big_Abs (X : Bignum) return Bignum is
+ (-Sec_Stack_Bignums.Big_Abs (+X));
+
+ function Big_EQ (X, Y : Bignum) return Boolean is
+ (Sec_Stack_Bignums.Big_EQ (+X, +Y));
+ function Big_NE (X, Y : Bignum) return Boolean is
+ (Sec_Stack_Bignums.Big_NE (+X, +Y));
+ function Big_GE (X, Y : Bignum) return Boolean is
+ (Sec_Stack_Bignums.Big_GE (+X, +Y));
+ function Big_LE (X, Y : Bignum) return Boolean is
+ (Sec_Stack_Bignums.Big_LE (+X, +Y));
+ function Big_GT (X, Y : Bignum) return Boolean is
+ (Sec_Stack_Bignums.Big_GT (+X, +Y));
+ function Big_LT (X, Y : Bignum) return Boolean is
+ (Sec_Stack_Bignums.Big_LT (+X, +Y));
function Bignum_In_LLI_Range (X : Bignum) return Boolean is
- begin
- -- If length is 0 or 1, definitely fits
-
- if X.Len <= 1 then
- return True;
-
- -- If length is greater than 2, definitely does not fit
-
- elsif X.Len > 2 then
- return False;
-
- -- Length is 2, more tests needed
-
- else
- declare
- Mag : constant DD := X.D (1) & X.D (2);
- begin
- return Mag < 2 ** 63 or else (X.Neg and then Mag = 2 ** 63);
- end;
- end if;
- end Bignum_In_LLI_Range;
-
- ---------------
- -- Normalize --
- ---------------
-
- function Normalize
- (X : Digit_Vector;
- Neg : Boolean := False) return Bignum
- is
- B : Bignum;
- J : Length;
-
- begin
- J := X'First;
- while J <= X'Last and then X (J) = 0 loop
- J := J + 1;
- end loop;
-
- B := Allocate_Bignum (X'Last - J + 1);
- B.Neg := B.Len > 0 and then Neg;
- B.D := X (J .. X'Last);
- return B;
- end Normalize;
-
- ---------------
- -- To_Bignum --
- ---------------
+ (Sec_Stack_Bignums.Bignum_In_LLI_Range (+X));
function To_Bignum (X : Long_Long_Integer) return Bignum is
- R : Bignum;
-
- begin
- if X = 0 then
- R := Allocate_Bignum (0);
+ (-Sec_Stack_Bignums.To_Bignum (X));
- -- One word result
-
- elsif X in -(2 ** 32 - 1) .. +(2 ** 32 - 1) then
- R := Allocate_Bignum (1);
- R.D (1) := SD (abs (X));
-
- -- Largest negative number annoyance
-
- elsif X = Long_Long_Integer'First then
- R := Allocate_Bignum (2);
- R.D (1) := 2 ** 31;
- R.D (2) := 0;
-
- -- Normal two word case
-
- else
- R := Allocate_Bignum (2);
- R.D (2) := SD (abs (X) mod Base);
- R.D (1) := SD (abs (X) / Base);
- end if;
-
- R.Neg := X < 0;
- return R;
- end To_Bignum;
+ function From_Bignum (X : Bignum) return Long_Long_Integer is
+ (Sec_Stack_Bignums.From_Bignum (+X));
end System.Bignums;
-- use in computing intermediate values in expressions for the case where
-- pragma Overflow_Check (Eliminated) is in effect.
-with Interfaces;
+-- Note that we cannot use a straight instantiation of System.Generic_Bignums
+-- because the rtsfind mechanism is not ready to handle instantiations.
package System.Bignums is
+ pragma Preelaborate;
- pragma Assert (Long_Long_Integer'Size = 64);
- -- This package assumes that Long_Long_Integer size is 64 bit (i.e. that it
- -- has a range of -2**63 to 2**63-1). The front end ensures that the mode
- -- ELIMINATED is not allowed for overflow checking if this is not the case.
-
- subtype Length is Natural range 0 .. 2 ** 23 - 1;
- -- Represent number of words in Digit_Vector
-
- Base : constant := 2 ** 32;
- -- Digit vectors use this base
-
- subtype SD is Interfaces.Unsigned_32;
- -- Single length digit
-
- type Digit_Vector is array (Length range <>) of SD;
- -- Represent digits of a number (most significant digit first)
-
- type Bignum_Data (Len : Length) is record
- Neg : Boolean;
- -- Set if value is negative, never set for zero
-
- D : Digit_Vector (1 .. Len);
- -- Digits of number, most significant first, represented in base
- -- 2**Base. No leading zeroes are stored, and the value of zero is
- -- represented using an empty vector for D.
- end record;
-
- for Bignum_Data use record
- Len at 0 range 0 .. 23;
- Neg at 3 range 0 .. 7;
- end record;
-
- type Bignum is access all Bignum_Data;
- -- This is the type that is used externally. Possibly this could be a
- -- private type, but we leave the structure exposed for now. For one
- -- thing it helps with debugging. Note that this package never shares
- -- an allocated Bignum value, so for example for X + 0, a copy of X is
- -- returned, not X itself.
-
- -- Note: none of the subprograms in this package modify the Bignum_Data
- -- records referenced by Bignum arguments of mode IN.
+ type Bignum is private;
function Big_Add (X, Y : Bignum) return Bignum; -- "+"
function Big_Sub (X, Y : Bignum) return Bignum; -- "-"
-- Convert Bignum to Long_Long_Integer. Constraint_Error raised with
-- appropriate message if value is out of range of Long_Long_Integer.
+private
+
+ type Bignum is new System.Address;
+
+ pragma Inline (Big_Add);
+ pragma Inline (Big_Sub);
+ pragma Inline (Big_Mul);
+ pragma Inline (Big_Div);
+ pragma Inline (Big_Exp);
+ pragma Inline (Big_Mod);
+ pragma Inline (Big_Rem);
+ pragma Inline (Big_Neg);
+ pragma Inline (Big_Abs);
+ pragma Inline (Big_EQ);
+ pragma Inline (Big_NE);
+ pragma Inline (Big_GE);
+ pragma Inline (Big_LE);
+ pragma Inline (Big_GT);
+ pragma Inline (Big_LT);
+ pragma Inline (Bignum_In_LLI_Range);
+ pragma Inline (To_Bignum);
+ pragma Inline (From_Bignum);
+
end System.Bignums;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- GNAT COMPILER COMPONENTS --
+-- --
+-- S Y S T E M . G E N E R I C _ B I G N U M S --
+-- --
+-- B o d y --
+-- --
+-- Copyright (C) 2012-2019, Free Software Foundation, Inc. --
+-- --
+-- GNAT is free software; you can redistribute it and/or modify it under --
+-- terms of the GNU General Public License as published by the Free Soft- --
+-- ware Foundation; either version 3, or (at your option) any later ver- --
+-- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
+-- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
+-- or FITNESS FOR A PARTICULAR PURPOSE. --
+-- --
+-- As a special exception under Section 7 of GPL version 3, you are granted --
+-- additional permissions described in the GCC Runtime Library Exception, --
+-- version 3.1, as published by the Free Software Foundation. --
+-- --
+-- You should have received a copy of the GNU General Public License and --
+-- a copy of the GCC Runtime Library Exception along with this program; --
+-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
+-- <http://www.gnu.org/licenses/>. --
+-- --
+-- GNAT was originally developed by the GNAT team at New York University. --
+-- Extensive contributions were provided by Ada Core Technologies Inc. --
+-- --
+------------------------------------------------------------------------------
+
+-- This package provides arbitrary precision signed integer arithmetic.
+
+with System; use System;
+with System.Secondary_Stack; use System.Secondary_Stack;
+with System.Storage_Elements; use System.Storage_Elements;
+
+package body System.Generic_Bignums is
+
+ use Interfaces;
+ -- So that operations on Unsigned_32/Unsigned_64 are available
+
+ type DD is mod Base ** 2;
+ -- Double length digit used for intermediate computations
+
+ function MSD (X : DD) return SD is (SD (X / Base));
+ function LSD (X : DD) return SD is (SD (X mod Base));
+ -- Most significant and least significant digit of double digit value
+
+ function "&" (X, Y : SD) return DD is (DD (X) * Base + DD (Y));
+ -- Compose double digit value from two single digit values
+
+ subtype LLI is Long_Long_Integer;
+
+ One_Data : constant Digit_Vector (1 .. 1) := (1 => 1);
+ -- Constant one
+
+ Zero_Data : constant Digit_Vector (1 .. 0) := (1 .. 0 => 0);
+ -- Constant zero
+
+ -----------------------
+ -- Local Subprograms --
+ -----------------------
+
+ function Add
+ (X, Y : Digit_Vector;
+ X_Neg : Boolean;
+ Y_Neg : Boolean) return Bignum
+ with
+ Pre => X'First = 1 and then Y'First = 1;
+ -- This procedure adds two signed numbers returning the Sum, it is used
+ -- for both addition and subtraction. The value computed is X + Y, with
+ -- X_Neg and Y_Neg giving the signs of the operands.
+
+ function Allocate_Bignum (Len : Length) return Bignum with
+ Post => Allocate_Bignum'Result.Len = Len;
+ -- Allocate Bignum value of indicated length on secondary stack. On return
+ -- the Neg and D fields are left uninitialized.
+
+ type Compare_Result is (LT, EQ, GT);
+ -- Indicates result of comparison in following call
+
+ function Compare
+ (X, Y : Digit_Vector;
+ X_Neg, Y_Neg : Boolean) return Compare_Result
+ with
+ Pre => X'First = 1 and then Y'First = 1;
+ -- Compare (X with sign X_Neg) with (Y with sign Y_Neg), and return the
+ -- result of the signed comparison.
+
+ procedure Div_Rem
+ (X, Y : Bignum;
+ Quotient : out Bignum;
+ Remainder : out Bignum;
+ Discard_Quotient : Boolean := False;
+ Discard_Remainder : Boolean := False);
+ -- Returns the Quotient and Remainder from dividing abs (X) by abs (Y). The
+ -- values of X and Y are not modified. If Discard_Quotient is True, then
+ -- Quotient is undefined on return, and if Discard_Remainder is True, then
+ -- Remainder is undefined on return. Service routine for Big_Div/Rem/Mod.
+
+ procedure Free_Bignum (X : Bignum) is null;
+ -- Called to free a Bignum value used in intermediate computations. In
+ -- this implementation using the secondary stack, it does nothing at all,
+ -- because we rely on Mark/Release, but it may be of use for some
+ -- alternative implementation.
+
+ function Normalize
+ (X : Digit_Vector;
+ Neg : Boolean := False) return Bignum;
+ -- Given a digit vector and sign, allocate and construct a Bignum value.
+ -- Note that X may have leading zeroes which must be removed, and if the
+ -- result is zero, the sign is forced positive.
+
+ ---------
+ -- Add --
+ ---------
+
+ function Add
+ (X, Y : Digit_Vector;
+ X_Neg : Boolean;
+ Y_Neg : Boolean) return Bignum
+ is
+ begin
+ -- If signs are the same, we are doing an addition, it is convenient to
+ -- ensure that the first operand is the longer of the two.
+
+ if X_Neg = Y_Neg then
+ if X'Last < Y'Last then
+ return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg);
+
+ -- Here signs are the same, and the first operand is the longer
+
+ else
+ pragma Assert (X_Neg = Y_Neg and then X'Last >= Y'Last);
+
+ -- Do addition, putting result in Sum (allowing for carry)
+
+ declare
+ Sum : Digit_Vector (0 .. X'Last);
+ RD : DD;
+
+ begin
+ RD := 0;
+ for J in reverse 1 .. X'Last loop
+ RD := RD + DD (X (J));
+
+ if J >= 1 + (X'Last - Y'Last) then
+ RD := RD + DD (Y (J - (X'Last - Y'Last)));
+ end if;
+
+ Sum (J) := LSD (RD);
+ RD := RD / Base;
+ end loop;
+
+ Sum (0) := SD (RD);
+ return Normalize (Sum, X_Neg);
+ end;
+ end if;
+
+ -- Signs are different so really this is a subtraction, we want to make
+ -- sure that the largest magnitude operand is the first one, and then
+ -- the result will have the sign of the first operand.
+
+ else
+ declare
+ CR : constant Compare_Result := Compare (X, Y, False, False);
+
+ begin
+ if CR = EQ then
+ return Normalize (Zero_Data);
+
+ elsif CR = LT then
+ return Add (X => Y, Y => X, X_Neg => Y_Neg, Y_Neg => X_Neg);
+
+ else
+ pragma Assert (X_Neg /= Y_Neg and then CR = GT);
+
+ -- Do subtraction, putting result in Diff
+
+ declare
+ Diff : Digit_Vector (1 .. X'Length);
+ RD : DD;
+
+ begin
+ RD := 0;
+ for J in reverse 1 .. X'Last loop
+ RD := RD + DD (X (J));
+
+ if J >= 1 + (X'Last - Y'Last) then
+ RD := RD - DD (Y (J - (X'Last - Y'Last)));
+ end if;
+
+ Diff (J) := LSD (RD);
+ RD := (if RD < Base then 0 else -1);
+ end loop;
+
+ return Normalize (Diff, X_Neg);
+ end;
+ end if;
+ end;
+ end if;
+ end Add;
+
+ ---------------------
+ -- Allocate_Bignum --
+ ---------------------
+
+ function Allocate_Bignum (Len : Length) return Bignum is
+ Addr : Address;
+
+ begin
+ -- Allocation on the heap
+
+ if not Use_Secondary_Stack then
+ declare
+ B : Bignum;
+ begin
+ B := new Bignum_Data'(Len, False, (others => 0));
+ return B;
+ end;
+
+ -- Allocation on the secondary stack
+
+ else
+ -- Note: The approach used here is designed to avoid strict aliasing
+ -- warnings that appeared previously using unchecked conversion.
+
+ SS_Allocate (Addr, Storage_Offset (4 + 4 * Len));
+
+ declare
+ B : Bignum;
+ for B'Address use Addr'Address;
+ pragma Import (Ada, B);
+
+ BD : Bignum_Data (Len);
+ for BD'Address use Addr;
+ pragma Import (Ada, BD);
+
+ -- Expose a writable view of discriminant BD.Len so that we can
+ -- initialize it. We need to use the exact layout of the record
+ -- to ensure that the Length field has 24 bits as expected.
+
+ type Bignum_Data_Header is record
+ Len : Length;
+ Neg : Boolean;
+ end record;
+
+ for Bignum_Data_Header use record
+ Len at 0 range 0 .. 23;
+ Neg at 3 range 0 .. 7;
+ end record;
+
+ BDH : Bignum_Data_Header;
+ for BDH'Address use BD'Address;
+ pragma Import (Ada, BDH);
+
+ pragma Assert (BDH.Len'Size = BD.Len'Size);
+
+ begin
+ BDH.Len := Len;
+ return B;
+ end;
+ end if;
+ end Allocate_Bignum;
+
+ -------------
+ -- Big_Abs --
+ -------------
+
+ function Big_Abs (X : Bignum) return Bignum is
+ begin
+ return Normalize (X.D);
+ end Big_Abs;
+
+ -------------
+ -- Big_Add --
+ -------------
+
+ function Big_Add (X, Y : Bignum) return Bignum is
+ begin
+ return Add (X.D, Y.D, X.Neg, Y.Neg);
+ end Big_Add;
+
+ -------------
+ -- Big_Div --
+ -------------
+
+ -- This table is excerpted from RM 4.5.5(28-30) and shows how the result
+ -- varies with the signs of the operands.
+
+ -- A B A/B A B A/B
+ --
+ -- 10 5 2 -10 5 -2
+ -- 11 5 2 -11 5 -2
+ -- 12 5 2 -12 5 -2
+ -- 13 5 2 -13 5 -2
+ -- 14 5 2 -14 5 -2
+ --
+ -- A B A/B A B A/B
+ --
+ -- 10 -5 -2 -10 -5 2
+ -- 11 -5 -2 -11 -5 2
+ -- 12 -5 -2 -12 -5 2
+ -- 13 -5 -2 -13 -5 2
+ -- 14 -5 -2 -14 -5 2
+
+ function Big_Div (X, Y : Bignum) return Bignum is
+ Q, R : Bignum;
+ begin
+ Div_Rem (X, Y, Q, R, Discard_Remainder => True);
+ Q.Neg := Q.Len > 0 and then (X.Neg xor Y.Neg);
+ return Q;
+ end Big_Div;
+
+ -------------
+ -- Big_Exp --
+ -------------
+
+ function Big_Exp (X, Y : Bignum) return Bignum is
+
+ function "**" (X : Bignum; Y : SD) return Bignum;
+ -- Internal routine where we know right operand is one word
+
+ ----------
+ -- "**" --
+ ----------
+
+ function "**" (X : Bignum; Y : SD) return Bignum is
+ begin
+ case Y is
+
+ -- X ** 0 is 1
+
+ when 0 =>
+ return Normalize (One_Data);
+
+ -- X ** 1 is X
+
+ when 1 =>
+ return Normalize (X.D);
+
+ -- X ** 2 is X * X
+
+ when 2 =>
+ return Big_Mul (X, X);
+
+ -- For X greater than 2, use the recursion
+
+ -- X even, X ** Y = (X ** (Y/2)) ** 2;
+ -- X odd, X ** Y = (X ** (Y/2)) ** 2 * X;
+
+ when others =>
+ declare
+ XY2 : constant Bignum := X ** (Y / 2);
+ XY2S : constant Bignum := Big_Mul (XY2, XY2);
+ Res : Bignum;
+
+ begin
+ Free_Bignum (XY2);
+
+ -- Raise storage error if intermediate value is getting too
+ -- large, which we arbitrarily define as 200 words for now.
+
+ if XY2S.Len > 200 then
+ Free_Bignum (XY2S);
+ raise Storage_Error with
+ "exponentiation result is too large";
+ end if;
+
+ -- Otherwise take care of even/odd cases
+
+ if (Y and 1) = 0 then
+ return XY2S;
+
+ else
+ Res := Big_Mul (XY2S, X);
+ Free_Bignum (XY2S);
+ return Res;
+ end if;
+ end;
+ end case;
+ end "**";
+
+ -- Start of processing for Big_Exp
+
+ begin
+ -- Error if right operand negative
+
+ if Y.Neg then
+ raise Constraint_Error with "exponentiation to negative power";
+
+ -- X ** 0 is always 1 (including 0 ** 0, so do this test first)
+
+ elsif Y.Len = 0 then
+ return Normalize (One_Data);
+
+ -- 0 ** X is always 0 (for X non-zero)
+
+ elsif X.Len = 0 then
+ return Normalize (Zero_Data);
+
+ -- (+1) ** Y = 1
+ -- (-1) ** Y = +/-1 depending on whether Y is even or odd
+
+ elsif X.Len = 1 and then X.D (1) = 1 then
+ return Normalize
+ (X.D, Neg => X.Neg and then ((Y.D (Y.Len) and 1) = 1));
+
+ -- If the absolute value of the base is greater than 1, then the
+ -- exponent must not be bigger than one word, otherwise the result
+ -- is ludicrously large, and we just signal Storage_Error right away.
+
+ elsif Y.Len > 1 then
+ raise Storage_Error with "exponentiation result is too large";
+
+ -- Special case (+/-)2 ** K, where K is 1 .. 31 using a shift
+
+ elsif X.Len = 1 and then X.D (1) = 2 and then Y.D (1) < 32 then
+ declare
+ D : constant Digit_Vector (1 .. 1) :=
+ (1 => Shift_Left (SD'(1), Natural (Y.D (1))));
+ begin
+ return Normalize (D, X.Neg);
+ end;
+
+ -- Remaining cases have right operand of one word
+
+ else
+ return X ** Y.D (1);
+ end if;
+ end Big_Exp;
+
+ ------------
+ -- Big_EQ --
+ ------------
+
+ function Big_EQ (X, Y : Bignum) return Boolean is
+ begin
+ return Compare (X.D, Y.D, X.Neg, Y.Neg) = EQ;
+ end Big_EQ;
+
+ ------------
+ -- Big_GE --
+ ------------
+
+ function Big_GE (X, Y : Bignum) return Boolean is
+ begin
+ return Compare (X.D, Y.D, X.Neg, Y.Neg) /= LT;
+ end Big_GE;
+
+ ------------
+ -- Big_GT --
+ ------------
+
+ function Big_GT (X, Y : Bignum) return Boolean is
+ begin
+ return Compare (X.D, Y.D, X.Neg, Y.Neg) = GT;
+ end Big_GT;
+
+ ------------
+ -- Big_LE --
+ ------------
+
+ function Big_LE (X, Y : Bignum) return Boolean is
+ begin
+ return Compare (X.D, Y.D, X.Neg, Y.Neg) /= GT;
+ end Big_LE;
+
+ ------------
+ -- Big_LT --
+ ------------
+
+ function Big_LT (X, Y : Bignum) return Boolean is
+ begin
+ return Compare (X.D, Y.D, X.Neg, Y.Neg) = LT;
+ end Big_LT;
+
+ -------------
+ -- Big_Mod --
+ -------------
+
+ -- This table is excerpted from RM 4.5.5(28-30) and shows how the result
+ -- of Rem and Mod vary with the signs of the operands.
+
+ -- A B A mod B A rem B A B A mod B A rem B
+
+ -- 10 5 0 0 -10 5 0 0
+ -- 11 5 1 1 -11 5 4 -1
+ -- 12 5 2 2 -12 5 3 -2
+ -- 13 5 3 3 -13 5 2 -3
+ -- 14 5 4 4 -14 5 1 -4
+
+ -- A B A mod B A rem B A B A mod B A rem B
+
+ -- 10 -5 0 0 -10 -5 0 0
+ -- 11 -5 -4 1 -11 -5 -1 -1
+ -- 12 -5 -3 2 -12 -5 -2 -2
+ -- 13 -5 -2 3 -13 -5 -3 -3
+ -- 14 -5 -1 4 -14 -5 -4 -4
+
+ function Big_Mod (X, Y : Bignum) return Bignum is
+ Q, R : Bignum;
+
+ begin
+ -- If signs are same, result is same as Rem
+
+ if X.Neg = Y.Neg then
+ return Big_Rem (X, Y);
+
+ -- Case where Mod is different
+
+ else
+ -- Do division
+
+ Div_Rem (X, Y, Q, R, Discard_Quotient => True);
+
+ -- Zero result is unchanged
+
+ if R.Len = 0 then
+ return R;
+
+ -- Otherwise adjust result
+
+ else
+ declare
+ T1 : constant Bignum := Big_Sub (Y, R);
+ begin
+ T1.Neg := Y.Neg;
+ Free_Bignum (R);
+ return T1;
+ end;
+ end if;
+ end if;
+ end Big_Mod;
+
+ -------------
+ -- Big_Mul --
+ -------------
+
+ function Big_Mul (X, Y : Bignum) return Bignum is
+ Result : Digit_Vector (1 .. X.Len + Y.Len) := (others => 0);
+ -- Accumulate result (max length of result is sum of operand lengths)
+
+ L : Length;
+ -- Current result digit
+
+ D : DD;
+ -- Result digit
+
+ begin
+ for J in 1 .. X.Len loop
+ for K in 1 .. Y.Len loop
+ L := Result'Last - (X.Len - J) - (Y.Len - K);
+ D := DD (X.D (J)) * DD (Y.D (K)) + DD (Result (L));
+ Result (L) := LSD (D);
+ D := D / Base;
+
+ -- D is carry which must be propagated
+
+ while D /= 0 and then L >= 1 loop
+ L := L - 1;
+ D := D + DD (Result (L));
+ Result (L) := LSD (D);
+ D := D / Base;
+ end loop;
+
+ -- Must not have a carry trying to extend max length
+
+ pragma Assert (D = 0);
+ end loop;
+ end loop;
+
+ -- Return result
+
+ return Normalize (Result, X.Neg xor Y.Neg);
+ end Big_Mul;
+
+ ------------
+ -- Big_NE --
+ ------------
+
+ function Big_NE (X, Y : Bignum) return Boolean is
+ begin
+ return Compare (X.D, Y.D, X.Neg, Y.Neg) /= EQ;
+ end Big_NE;
+
+ -------------
+ -- Big_Neg --
+ -------------
+
+ function Big_Neg (X : Bignum) return Bignum is
+ begin
+ return Normalize (X.D, not X.Neg);
+ end Big_Neg;
+
+ -------------
+ -- Big_Rem --
+ -------------
+
+ -- This table is excerpted from RM 4.5.5(28-30) and shows how the result
+ -- varies with the signs of the operands.
+
+ -- A B A rem B A B A rem B
+
+ -- 10 5 0 -10 5 0
+ -- 11 5 1 -11 5 -1
+ -- 12 5 2 -12 5 -2
+ -- 13 5 3 -13 5 -3
+ -- 14 5 4 -14 5 -4
+
+ -- A B A rem B A B A rem B
+
+ -- 10 -5 0 -10 -5 0
+ -- 11 -5 1 -11 -5 -1
+ -- 12 -5 2 -12 -5 -2
+ -- 13 -5 3 -13 -5 -3
+ -- 14 -5 4 -14 -5 -4
+
+ function Big_Rem (X, Y : Bignum) return Bignum is
+ Q, R : Bignum;
+ begin
+ Div_Rem (X, Y, Q, R, Discard_Quotient => True);
+ R.Neg := R.Len > 0 and then X.Neg;
+ return R;
+ end Big_Rem;
+
+ -------------
+ -- Big_Sub --
+ -------------
+
+ function Big_Sub (X, Y : Bignum) return Bignum is
+ begin
+ -- If right operand zero, return left operand (avoiding sharing)
+
+ if Y.Len = 0 then
+ return Normalize (X.D, X.Neg);
+
+ -- Otherwise add negative of right operand
+
+ else
+ return Add (X.D, Y.D, X.Neg, not Y.Neg);
+ end if;
+ end Big_Sub;
+
+ -------------
+ -- Compare --
+ -------------
+
+ function Compare
+ (X, Y : Digit_Vector;
+ X_Neg, Y_Neg : Boolean) return Compare_Result
+ is
+ begin
+ -- Signs are different, that's decisive, since 0 is always plus
+
+ if X_Neg /= Y_Neg then
+ return (if X_Neg then LT else GT);
+
+ -- Lengths are different, that's decisive since no leading zeroes
+
+ elsif X'Last /= Y'Last then
+ return (if (X'Last > Y'Last) xor X_Neg then GT else LT);
+
+ -- Need to compare data
+
+ else
+ for J in X'Range loop
+ if X (J) /= Y (J) then
+ return (if (X (J) > Y (J)) xor X_Neg then GT else LT);
+ end if;
+ end loop;
+
+ return EQ;
+ end if;
+ end Compare;
+
+ -------------
+ -- Div_Rem --
+ -------------
+
+ procedure Div_Rem
+ (X, Y : Bignum;
+ Quotient : out Bignum;
+ Remainder : out Bignum;
+ Discard_Quotient : Boolean := False;
+ Discard_Remainder : Boolean := False)
+ is
+ begin
+ -- Error if division by zero
+
+ if Y.Len = 0 then
+ raise Constraint_Error with "division by zero";
+ end if;
+
+ -- Handle simple cases with special tests
+
+ -- If X < Y then quotient is zero and remainder is X
+
+ if Compare (X.D, Y.D, False, False) = LT then
+ Remainder := Normalize (X.D);
+ Quotient := Normalize (Zero_Data);
+ return;
+
+ -- If both X and Y are less than 2**63-1, we can use Long_Long_Integer
+ -- arithmetic. Note it is good not to do an accurate range check against
+ -- Long_Long_Integer since -2**63 / -1 overflows.
+
+ elsif (X.Len <= 1 or else (X.Len = 2 and then X.D (1) < 2**31))
+ and then
+ (Y.Len <= 1 or else (Y.Len = 2 and then Y.D (1) < 2**31))
+ then
+ declare
+ A : constant LLI := abs (From_Bignum (X));
+ B : constant LLI := abs (From_Bignum (Y));
+ begin
+ Quotient := To_Bignum (A / B);
+ Remainder := To_Bignum (A rem B);
+ return;
+ end;
+
+ -- Easy case if divisor is one digit
+
+ elsif Y.Len = 1 then
+ declare
+ ND : DD;
+ Div : constant DD := DD (Y.D (1));
+
+ Result : Digit_Vector (1 .. X.Len);
+ Remdr : Digit_Vector (1 .. 1);
+
+ begin
+ ND := 0;
+ for J in 1 .. X.Len loop
+ ND := Base * ND + DD (X.D (J));
+ Result (J) := SD (ND / Div);
+ ND := ND rem Div;
+ end loop;
+
+ Quotient := Normalize (Result);
+ Remdr (1) := SD (ND);
+ Remainder := Normalize (Remdr);
+ return;
+ end;
+ end if;
+
+ -- The complex full multi-precision case. We will employ algorithm
+ -- D defined in the section "The Classical Algorithms" (sec. 4.3.1)
+ -- of Donald Knuth's "The Art of Computer Programming", Vol. 2, 2nd
+ -- edition. The terminology is adjusted for this section to match that
+ -- reference.
+
+ -- We are dividing X.Len digits of X (called u here) by Y.Len digits
+ -- of Y (called v here), developing the quotient and remainder. The
+ -- numbers are represented using Base, which was chosen so that we have
+ -- the operations of multiplying to single digits (SD) to form a double
+ -- digit (DD), and dividing a double digit (DD) by a single digit (SD)
+ -- to give a single digit quotient and a single digit remainder.
+
+ -- Algorithm D from Knuth
+
+ -- Comments here with square brackets are directly from Knuth
+
+ Algorithm_D : declare
+
+ -- The following lower case variables correspond exactly to the
+ -- terminology used in algorithm D.
+
+ m : constant Length := X.Len - Y.Len;
+ n : constant Length := Y.Len;
+ b : constant DD := Base;
+
+ u : Digit_Vector (0 .. m + n);
+ v : Digit_Vector (1 .. n);
+ q : Digit_Vector (0 .. m);
+ r : Digit_Vector (1 .. n);
+
+ u0 : SD renames u (0);
+ v1 : SD renames v (1);
+ v2 : SD renames v (2);
+
+ d : DD;
+ j : Length;
+ qhat : DD;
+ rhat : DD;
+ temp : DD;
+
+ begin
+ -- Initialize data of left and right operands
+
+ for J in 1 .. m + n loop
+ u (J) := X.D (J);
+ end loop;
+
+ for J in 1 .. n loop
+ v (J) := Y.D (J);
+ end loop;
+
+ -- [Division of nonnegative integers.] Given nonnegative integers u
+ -- = (ul,u2..um+n) and v = (v1,v2..vn), where v1 /= 0 and n > 1, we
+ -- form the quotient u / v = (q0,ql..qm) and the remainder u mod v =
+ -- (r1,r2..rn).
+
+ pragma Assert (v1 /= 0);
+ pragma Assert (n > 1);
+
+ -- Dl. [Normalize.] Set d = b/(vl + 1). Then set (u0,u1,u2..um+n)
+ -- equal to (u1,u2..um+n) times d, and set (v1,v2..vn) equal to
+ -- (v1,v2..vn) times d. Note the introduction of a new digit position
+ -- u0 at the left of u1; if d = 1 all we need to do in this step is
+ -- to set u0 = 0.
+
+ d := b / (DD (v1) + 1);
+
+ if d = 1 then
+ u0 := 0;
+
+ else
+ declare
+ Carry : DD;
+ Tmp : DD;
+
+ begin
+ -- Multiply Dividend (u) by d
+
+ Carry := 0;
+ for J in reverse 1 .. m + n loop
+ Tmp := DD (u (J)) * d + Carry;
+ u (J) := LSD (Tmp);
+ Carry := Tmp / Base;
+ end loop;
+
+ u0 := SD (Carry);
+
+ -- Multiply Divisor (v) by d
+
+ Carry := 0;
+ for J in reverse 1 .. n loop
+ Tmp := DD (v (J)) * d + Carry;
+ v (J) := LSD (Tmp);
+ Carry := Tmp / Base;
+ end loop;
+
+ pragma Assert (Carry = 0);
+ end;
+ end if;
+
+ -- D2. [Initialize j.] Set j = 0. The loop on j, steps D2 through D7,
+ -- will be essentially a division of (uj, uj+1..uj+n) by (v1,v2..vn)
+ -- to get a single quotient digit qj.
+
+ j := 0;
+
+ -- Loop through digits
+
+ loop
+ -- Note: In the original printing, step D3 was as follows:
+
+ -- D3. [Calculate qhat.] If uj = v1, set qhat to b-l; otherwise
+ -- set qhat to (uj,uj+1)/v1. Now test if v2 * qhat is greater than
+ -- (uj*b + uj+1 - qhat*v1)*b + uj+2. If so, decrease qhat by 1 and
+ -- repeat this test
+
+ -- This had a bug not discovered till 1995, see Vol 2 errata:
+ -- http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz. Under
+ -- rare circumstances the expression in the test could overflow.
+ -- This version was further corrected in 2005, see Vol 2 errata:
+ -- http://www-cs-faculty.stanford.edu/~uno/all2-pre.ps.gz.
+ -- The code below is the fixed version of this step.
+
+ -- D3. [Calculate qhat.] Set qhat to (uj,uj+1)/v1 and rhat to
+ -- to (uj,uj+1) mod v1.
+
+ temp := u (j) & u (j + 1);
+ qhat := temp / DD (v1);
+ rhat := temp mod DD (v1);
+
+ -- D3 (continued). Now test if qhat >= b or v2*qhat > (rhat,uj+2):
+ -- if so, decrease qhat by 1, increase rhat by v1, and repeat this
+ -- test if rhat < b. [The test on v2 determines at high speed
+ -- most of the cases in which the trial value qhat is one too
+ -- large, and eliminates all cases where qhat is two too large.]
+
+ while qhat >= b
+ or else DD (v2) * qhat > LSD (rhat) & u (j + 2)
+ loop
+ qhat := qhat - 1;
+ rhat := rhat + DD (v1);
+ exit when rhat >= b;
+ end loop;
+
+ -- D4. [Multiply and subtract.] Replace (uj,uj+1..uj+n) by
+ -- (uj,uj+1..uj+n) minus qhat times (v1,v2..vn). This step
+ -- consists of a simple multiplication by a one-place number,
+ -- combined with a subtraction.
+
+ -- The digits (uj,uj+1..uj+n) are always kept positive; if the
+ -- result of this step is actually negative then (uj,uj+1..uj+n)
+ -- is left as the true value plus b**(n+1), i.e. as the b's
+ -- complement of the true value, and a "borrow" to the left is
+ -- remembered.
+
+ declare
+ Borrow : SD;
+ Carry : DD;
+ Temp : DD;
+
+ Negative : Boolean;
+ -- Records if subtraction causes a negative result, requiring
+ -- an add back (case where qhat turned out to be 1 too large).
+
+ begin
+ Borrow := 0;
+ for K in reverse 1 .. n loop
+ Temp := qhat * DD (v (K)) + DD (Borrow);
+ Borrow := MSD (Temp);
+
+ if LSD (Temp) > u (j + K) then
+ Borrow := Borrow + 1;
+ end if;
+
+ u (j + K) := u (j + K) - LSD (Temp);
+ end loop;
+
+ Negative := u (j) < Borrow;
+ u (j) := u (j) - Borrow;
+
+ -- D5. [Test remainder.] Set qj = qhat. If the result of step
+ -- D4 was negative, we will do the add back step (step D6).
+
+ q (j) := LSD (qhat);
+
+ if Negative then
+
+ -- D6. [Add back.] Decrease qj by 1, and add (0,v1,v2..vn)
+ -- to (uj,uj+1,uj+2..uj+n). (A carry will occur to the left
+ -- of uj, and it is be ignored since it cancels with the
+ -- borrow that occurred in D4.)
+
+ q (j) := q (j) - 1;
+
+ Carry := 0;
+ for K in reverse 1 .. n loop
+ Temp := DD (v (K)) + DD (u (j + K)) + Carry;
+ u (j + K) := LSD (Temp);
+ Carry := Temp / Base;
+ end loop;
+
+ u (j) := u (j) + SD (Carry);
+ end if;
+ end;
+
+ -- D7. [Loop on j.] Increase j by one. Now if j <= m, go back to
+ -- D3 (the start of the loop on j).
+
+ j := j + 1;
+ exit when not (j <= m);
+ end loop;
+
+ -- D8. [Unnormalize.] Now (qo,ql..qm) is the desired quotient, and
+ -- the desired remainder may be obtained by dividing (um+1..um+n)
+ -- by d.
+
+ if not Discard_Quotient then
+ Quotient := Normalize (q);
+ end if;
+
+ if not Discard_Remainder then
+ declare
+ Remdr : DD;
+
+ begin
+ Remdr := 0;
+ for K in 1 .. n loop
+ Remdr := Base * Remdr + DD (u (m + K));
+ r (K) := SD (Remdr / d);
+ Remdr := Remdr rem d;
+ end loop;
+
+ pragma Assert (Remdr = 0);
+ end;
+
+ Remainder := Normalize (r);
+ end if;
+ end Algorithm_D;
+ end Div_Rem;
+
+ -----------------
+ -- From_Bignum --
+ -----------------
+
+ function From_Bignum (X : Bignum) return Long_Long_Integer is
+ begin
+ if X.Len = 0 then
+ return 0;
+
+ elsif X.Len = 1 then
+ return (if X.Neg then -LLI (X.D (1)) else LLI (X.D (1)));
+
+ elsif X.Len = 2 then
+ declare
+ Mag : constant DD := X.D (1) & X.D (2);
+ begin
+ if X.Neg and then Mag <= 2 ** 63 then
+ return -LLI (Mag);
+ elsif Mag < 2 ** 63 then
+ return LLI (Mag);
+ end if;
+ end;
+ end if;
+
+ raise Constraint_Error with "expression value out of range";
+ end From_Bignum;
+
+ -------------------------
+ -- Bignum_In_LLI_Range --
+ -------------------------
+
+ function Bignum_In_LLI_Range (X : Bignum) return Boolean is
+ begin
+ -- If length is 0 or 1, definitely fits
+
+ if X.Len <= 1 then
+ return True;
+
+ -- If length is greater than 2, definitely does not fit
+
+ elsif X.Len > 2 then
+ return False;
+
+ -- Length is 2, more tests needed
+
+ else
+ declare
+ Mag : constant DD := X.D (1) & X.D (2);
+ begin
+ return Mag < 2 ** 63 or else (X.Neg and then Mag = 2 ** 63);
+ end;
+ end if;
+ end Bignum_In_LLI_Range;
+
+ ---------------
+ -- Normalize --
+ ---------------
+
+ function Normalize
+ (X : Digit_Vector;
+ Neg : Boolean := False) return Bignum
+ is
+ B : Bignum;
+ J : Length;
+
+ begin
+ J := X'First;
+ while J <= X'Last and then X (J) = 0 loop
+ J := J + 1;
+ end loop;
+
+ B := Allocate_Bignum (X'Last - J + 1);
+ B.Neg := B.Len > 0 and then Neg;
+ B.D := X (J .. X'Last);
+ return B;
+ end Normalize;
+
+ ---------------
+ -- To_Bignum --
+ ---------------
+
+ function To_Bignum (X : Long_Long_Integer) return Bignum is
+ R : Bignum;
+
+ begin
+ if X = 0 then
+ R := Allocate_Bignum (0);
+
+ -- One word result
+
+ elsif X in -(2 ** 32 - 1) .. +(2 ** 32 - 1) then
+ R := Allocate_Bignum (1);
+ R.D (1) := SD (abs (X));
+
+ -- Largest negative number annoyance
+
+ elsif X = Long_Long_Integer'First then
+ R := Allocate_Bignum (2);
+ R.D (1) := 2 ** 31;
+ R.D (2) := 0;
+
+ -- Normal two word case
+
+ else
+ R := Allocate_Bignum (2);
+ R.D (2) := SD (abs (X) mod Base);
+ R.D (1) := SD (abs (X) / Base);
+ end if;
+
+ R.Neg := X < 0;
+ return R;
+ end To_Bignum;
+
+ function To_Bignum (X : Unsigned_64) return Bignum is
+ R : Bignum;
+
+ begin
+ if X = 0 then
+ R := Allocate_Bignum (0);
+
+ -- One word result
+
+ elsif X < 2 ** 32 then
+ R := Allocate_Bignum (1);
+ R.D (1) := SD (X);
+
+ -- Two word result
+
+ else
+ R := Allocate_Bignum (2);
+ R.D (2) := SD (X mod Base);
+ R.D (1) := SD (X / Base);
+ end if;
+
+ R.Neg := False;
+ return R;
+ end To_Bignum;
+
+ -------------
+ -- Is_Zero --
+ -------------
+
+ function Is_Zero (X : Bignum) return Boolean is
+ (X /= null and then X.D = Zero_Data);
+
+end System.Generic_Bignums;
--- /dev/null
+------------------------------------------------------------------------------
+-- --
+-- GNAT COMPILER COMPONENTS --
+-- --
+-- S Y S T E M . G E N E R I C _ B I G N U M S --
+-- --
+-- S p e c --
+-- --
+-- Copyright (C) 2012-2019, Free Software Foundation, Inc. --
+-- --
+-- GNAT is free software; you can redistribute it and/or modify it under --
+-- terms of the GNU General Public License as published by the Free Soft- --
+-- ware Foundation; either version 3, or (at your option) any later ver- --
+-- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
+-- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
+-- or FITNESS FOR A PARTICULAR PURPOSE. --
+-- --
+-- As a special exception under Section 7 of GPL version 3, you are granted --
+-- additional permissions described in the GCC Runtime Library Exception, --
+-- version 3.1, as published by the Free Software Foundation. --
+-- --
+-- You should have received a copy of the GNU General Public License and --
+-- a copy of the GCC Runtime Library Exception along with this program; --
+-- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
+-- <http://www.gnu.org/licenses/>. --
+-- --
+-- GNAT was originally developed by the GNAT team at New York University. --
+-- Extensive contributions were provided by Ada Core Technologies Inc. --
+-- --
+------------------------------------------------------------------------------
+
+-- This package provides arbitrary precision signed integer arithmetic
+-- and can be used either built into the compiler via System.Bignums or to
+-- implement a default version of Ada.Numerics.Big_Numbers.Big_Integers.
+
+-- If Use_Secondary_Stack is True then all Bignum values are allocated on the
+-- secondary stack. If False, the heap is used and the caller is responsible
+-- for memory management.
+
+with Ada.Unchecked_Conversion;
+with Interfaces;
+
+generic
+ Use_Secondary_Stack : Boolean;
+package System.Generic_Bignums is
+ pragma Preelaborate;
+
+ pragma Assert (Long_Long_Integer'Size = 64);
+ -- This package assumes that Long_Long_Integer size is 64 bit (i.e. that it
+ -- has a range of -2**63 to 2**63-1). The front end ensures that the mode
+ -- ELIMINATED is not allowed for overflow checking if this is not the case.
+
+ subtype Length is Natural range 0 .. 2 ** 23 - 1;
+ -- Represent number of words in Digit_Vector
+
+ Base : constant := 2 ** 32;
+ -- Digit vectors use this base
+
+ subtype SD is Interfaces.Unsigned_32;
+ -- Single length digit
+
+ type Digit_Vector is array (Length range <>) of SD;
+ -- Represent digits of a number (most significant digit first)
+
+ type Bignum_Data (Len : Length) is record
+ Neg : Boolean;
+ -- Set if value is negative, never set for zero
+
+ D : Digit_Vector (1 .. Len);
+ -- Digits of number, most significant first, represented in base
+ -- 2**Base. No leading zeroes are stored, and the value of zero is
+ -- represented using an empty vector for D.
+ end record;
+
+ for Bignum_Data use record
+ Len at 0 range 0 .. 23;
+ Neg at 3 range 0 .. 7;
+ end record;
+
+ type Bignum is access all Bignum_Data;
+ -- This is the type that is used externally. Possibly this could be a
+ -- private type, but we leave the structure exposed for now. For one
+ -- thing it helps with debugging. Note that this package never shares
+ -- an allocated Bignum value, so for example for X + 0, a copy of X is
+ -- returned, not X itself.
+
+ function To_Bignum is new Ada.Unchecked_Conversion (System.Address, Bignum);
+ function To_Address is new
+ Ada.Unchecked_Conversion (Bignum, System.Address);
+
+ -- Note: none of the subprograms in this package modify the Bignum_Data
+ -- records referenced by Bignum arguments of mode IN.
+
+ function Big_Add (X, Y : Bignum) return Bignum; -- "+"
+ function Big_Sub (X, Y : Bignum) return Bignum; -- "-"
+ function Big_Mul (X, Y : Bignum) return Bignum; -- "*"
+ function Big_Div (X, Y : Bignum) return Bignum; -- "/"
+ function Big_Exp (X, Y : Bignum) return Bignum; -- "**"
+ function Big_Mod (X, Y : Bignum) return Bignum; -- "mod"
+ function Big_Rem (X, Y : Bignum) return Bignum; -- "rem"
+ function Big_Neg (X : Bignum) return Bignum; -- "-"
+ function Big_Abs (X : Bignum) return Bignum; -- "abs"
+ -- Perform indicated arithmetic operation on bignum values. No exception
+ -- raised except for Div/Mod/Rem by 0 which raises Constraint_Error with
+ -- an appropriate message.
+
+ function Big_EQ (X, Y : Bignum) return Boolean; -- "="
+ function Big_NE (X, Y : Bignum) return Boolean; -- "/="
+ function Big_GE (X, Y : Bignum) return Boolean; -- ">="
+ function Big_LE (X, Y : Bignum) return Boolean; -- "<="
+ function Big_GT (X, Y : Bignum) return Boolean; -- ">"
+ function Big_LT (X, Y : Bignum) return Boolean; -- "<"
+ -- Perform indicated comparison on bignums, returning result as Boolean.
+ -- No exception raised for any input arguments.
+
+ function Bignum_In_LLI_Range (X : Bignum) return Boolean;
+ -- Returns True if the Bignum value is in the range of Long_Long_Integer,
+ -- so that a call to From_Bignum is guaranteed not to raise an exception.
+
+ function To_Bignum (X : Long_Long_Integer) return Bignum;
+ -- Convert Long_Long_Integer to Bignum. No exception can be raised for any
+ -- input argument.
+
+ function To_Bignum (X : Interfaces.Unsigned_64) return Bignum;
+ -- Convert Unsigned_64 to Bignum. No exception can be raised for any
+ -- input argument.
+
+ function From_Bignum (X : Bignum) return Long_Long_Integer;
+ -- Convert Bignum to Long_Long_Integer. Constraint_Error raised with
+ -- appropriate message if value is out of range of Long_Long_Integer.
+
+ function Is_Zero (X : Bignum) return Boolean;
+ -- Return True if X = 0
+
+end System.Generic_Bignums;