CVS difference for ai12s/ai12-0208-1.txt

Differences between 1.4 and version 1.5
Log of other versions for file ai12s/ai12-0208-1.txt

--- ai12s/ai12-0208-1.txt	2018/02/27 04:29:33	1.4
+++ ai12s/ai12-0208-1.txt	2018/02/27 07:03:53	1.5
@@ -1180,3 +1180,1488 @@
 sense.
 
 ****************************************************************
+
+From: John Barnes
+Sent: Tuesday, January 30, 2018  4:22 AM
+
+I am feverishly giving lectures on numbers at Oxford at the moment but I am 
+trying to keep an eye on what the ARG is up to.
+
+Did you know that a new Mersenne prime was discovered on Boxing Day (26
+December) 2017. It is 2**77232917 - 1 and has only 23,249,425 digits. Will
+the Bignum package cope with it?
+
+****************************************************************
+
+From: Tucker Taft
+Sent: Tuesday, January 30, 2018  3:02 PM
+
+Yes, I noticed that new Mersenne prime as well.  And 23 mega-digit is nothing
+for a modern iPhone. ;-)  Just be sure to set aside a bit of extra time to
+print it out using the Image function.  Except in base 2, of course, which
+I could do right now.  Ready:  1111111111 [... 77,232,900 1's] 1111111!
+
+****************************************************************
+
+From: Jeff Cousins
+Sent: Wednesday, January 31, 2018  6:31 AM
+
+[This is John Barnes' Bignum package and some test programs - Editor.]
+
+-- file books\fun\progs\numbers.ada
+
+-- Restructured using children
+
+-- Types and No_Of_Places in parent package 
+
+-- 20-10-06
+
+package Numbers is
+
+   Max_Index: constant := 1000;
+   subtype Index is Integer range 0 .. Max_Index;
+   type Number(Max_Digits: Index := 1) is private;
+
+   Zero, One: constant Number;
+
+   Number_Error : exception;
+
+private
+   Base_Exp: constant := 4;
+   Base: constant := 10 ** Base_Exp;
+
+   type Base_Digit is range -Base .. 2 * Base - 1;
+
+   type Base_Digit_Array is
+      array(Index range <>) of Base_Digit;
+
+   type Number(Max_Digits: Index := 1) is
+      record
+         Sign: Integer := +1;
+         Length: Index := 0;
+         D: Base_Digit_Array(1..Max_Digits);
+      end record;
+
+   Zero: constant Number := (0, +1, 0, (others => 0));
+   One: constant Number := (1, +1, 1, (1 => 1));
+
+   function No_Of_Places(N: Number) return Integer;
+
+end Numbers;
+
+
+package body Numbers is
+
+   function No_Of_Places(N: Number) return Integer is
+   begin
+      if N.Length = 0 then 
+         return 1;
+      else
+         return N.Length * Base_Exp;
+      end if;
+   end No_Of_Places;
+
+end Numbers;
+
+--------------------------------------------------------------------
+
+package Numbers.Proc is
+
+   subtype Index is Numbers.Index;
+   subtype Number is Numbers.Number;
+   Zero: Number renames Numbers.Zero;
+   One: Number renames Numbers.One;
+   Number_Error: exception renames Numbers.Number_Error;
+
+   procedure Add(X, Y: Number; Z: out Number);
+   procedure Sub(X, Y: Number; Z: out Number);
+   procedure Mul(X, Y: Number; Z: out Number);
+   procedure Div(X, Y: Number; Quotient,
+                               Remainder: out Number);
+
+   procedure Neg(X: in out Number);
+
+   function Compare(X, Y: Number) return Integer;
+
+   function Length(N: Number) return Index;
+
+   procedure To_Number(S: String; N: out Number);
+   procedure To_Number(I: Integer; N: out Number);
+   procedure To_Text(N: Number; S: out String);
+   procedure To_Integer(N: Number; I: out Integer);
+
+end Numbers.Proc;
+
+--------------------------------------------------------------------
+
+package Numbers.IO is
+
+   Default_Width: Natural := 0;
+   procedure Put(Item: Number;
+                Width: Natural := Default_Width);
+   procedure Get(Item: out Number);
+
+end Numbers.IO;
+
+--------------------------------------------------------------------
+
+package Numbers.Func is
+
+   subtype Number is Numbers.Number;
+   Zero: Number renames Numbers.Zero;
+   One: Number renames Numbers.One;
+   Number_Error: exception renames Numbers.Number_Error;
+
+   function "+" (X: Number) return Number;
+   function "-" (X: Number) return Number;
+   function "abs" (X: Number) return Number;
+   function "+" (X, Y: Number) return Number;
+   function "-" (X, Y: Number) return Number;
+   function "*" (X, Y: Number) return Number;
+   function "/" (X, Y: Number) return Number;
+   function "rem" (X, Y: Number) return Number;
+   function "mod" (X, Y: Number) return Number;
+   function "**" (X: Number; N: Natural) return Number;
+   function "<" (X, Y: Number) return Boolean;
+   function "<=" (X, Y: Number) return Boolean;
+   function ">" (X, Y: Number) return Boolean;
+   function ">=" (X, Y: Number) return Boolean;
+
+   function "=" (X, Y: Number) return Boolean;
+
+   function Make_Number(S: String) return Number;
+   function Make_Number(I: Integer) return Number;
+   function String_Of(N: Number) return String;
+   function Integer_Of(N: Number) return Integer;
+
+end Numbers.Func;
+
+--------------------------------------------------------------------
+
+package body Numbers.Proc is
+
+   Base_Squared: constant := Base**2;
+   subtype Single is Base_Digit;
+   type Double is range -Base_Squared .. 2*Base_Squared - 1;
+
+
+   function Unsigned_Compare(X, Y: Number) return Integer is
+   -- ignoring signs
+   -- returns +1, 0 or -1 according as X >, = or < Y
+   begin
+      if X.Length > Y.Length then return +1; end if;
+      if X.Length < Y.Length then return -1; end if;
+      for I in reverse 1 .. X.Length loop
+         if X.D(I) > Y.D(I) then return +1; end if;
+         if X.D(I) < Y.D(I) then return -1; end if;
+      end loop;
+      return 0;  -- the numbers are equal
+   end Unsigned_Compare;
+
+   function Compare(X, Y: Number) return Integer is
+   -- returns +1, 0 or -1 according as X >, = or < Y
+   begin
+      if X.Sign /= Y.Sign then return X.Sign; end if;
+      return Unsigned_Compare(X, Y) * X.Sign;
+   end Compare;
+
+   procedure Raw_Add(X, Y: Number; Z: out Number) is
+   -- assumes X not smaller than Y
+      Carry: Single := 0;
+      Digit: Single;
+      ZL: Index := X.Length;  -- length of answer
+   begin
+      if Z.Max_Digits < ZL then
+         raise Number_Error;  -- Z not big enough to hold X
+      end if;
+      for I in 1 .. ZL loop
+         Digit := X.D(I) + Carry;
+         if I <= Y.Length then
+            Digit := Digit + Y.D(I);
+         end if;
+         if Digit >= Base then
+            Carry := 1;  Digit := Digit - Base;
+         else
+            Carry := 0;
+         end if;
+         Z.D(I) := Digit;
+      end loop;
+      if Carry /= 0 then
+         if ZL = Z.Max_Digits then
+            raise Number_Error;  --  too big to fit in Z
+         end if;
+         ZL := ZL + 1;
+         Z.D(ZL) := Carry;
+      end if;
+      Z.Length := ZL;
+   end Raw_Add;
+
+   procedure Raw_Sub(X, Y: Number; Z: out Number) is
+   -- assumes X not smaller than Y
+      Carry: Single := 0;
+      Digit: Single;
+      ZL: Index := X.Length;  -- length of answer
+   begin
+      if Z.Max_Digits < ZL then
+         raise Number_Error;  -- Z not big enough to hold X
+      end if;
+      for I in 1 .. ZL loop
+         Digit := X.D(I) - Carry;
+         if I <= Y.Length then
+            Digit := Digit - Y.D(I);
+         end if;
+         if Digit < 0 then
+            Carry := 1;  Digit := Digit + Base;
+         else
+            Carry := 0;
+         end if;
+         Z.D(I) := Digit;
+      end loop;
+      while Z.D(ZL) = 0 loop  -- SHOULD THIS NOT FAIL???
+         ZL := ZL - 1;         -- remove leading zeroes
+      end loop;
+      Z.Length := ZL;
+   end Raw_Sub;
+
+   procedure Add(X, Y: Number; Z: out Number) is
+      UCMPXY: Integer := Unsigned_Compare(X, Y);
+   begin
+      if X.Sign = Y.Sign then
+         Z.Sign := X.Sign;
+         if UCMPXY >= 0 then
+            Raw_Add(X, Y, Z);
+         else
+            Raw_Add(Y, X, Z);  -- reverse if Y larger
+         end if;
+      else
+         if UCMPXY > 0 then
+            Raw_Sub(X, Y, Z);
+            Z.Sign := X.Sign;
+         elsif UCMPXY < 0 then
+            Raw_Sub(Y, X, Z);
+            Z.Sign := -X.Sign;
+         else  --  answer is zero
+            Z.Sign := +1; Z.Length := 0;
+         end if;
+      end if;
+   end Add;
+
+   procedure Sub(X, Y: Number; Z: out Number) is
+      UCMPXY: Integer := Unsigned_Compare(X, Y);
+   begin
+      if X.Sign /= Y.Sign then
+         Z.Sign := X.Sign;
+         if UCMPXY >= 0 then
+            Raw_Add(X, Y, Z);
+         else
+            Raw_Add(Y, X, Z);  -- reverse if Y larger
+         end if;
+      else
+         if UCMPXY > 0 then
+            Raw_Sub(X, Y, Z);
+            Z.Sign := X.Sign;
+         elsif UCMPXY < 0 then
+            Raw_Sub(Y, X, Z);
+            Z.Sign := -X.Sign;
+         else  --  answer is zero
+            Z.Sign := +1; Z.Length := 0;
+         end if;
+      end if;
+   end Sub;
+
+   procedure Neg(X: in out Number) is
+   begin
+      -- do nothing in zero case
+      if X.Length = 0 then return; end if;
+      X.Sign := -X.Sign;
+   end Neg;
+
+   function Length(N: Number) return Index is
+   begin
+      return N.Length;
+   end Length;
+
+   procedure Mul(X, Y: Number; Z: out Number) is
+
+      Carry: Double;
+      Digit: Double;
+      ZL: Index;
+   begin
+      if Z.Max_Digits < X.Length + Y.Length then
+         raise Number_Error;
+      end if;
+      if X.Length = 0 or Y.Length = 0 then  -- zero case
+         Z.Sign := +1;  Z.Length := 0;
+         return;
+      end if;
+
+      ZL := X.Length + Y.Length - 1;
+                   -- lower possible length of answer
+
+      -- copy X to top of Z; so X and Z can be same array
+      for I in reverse 1 .. X.Length loop
+         Z.D(I + Y.Length) := X.D(I);
+      end loop;
+      declare -- initialise limits and length of cycle
+         Z_Index: Index;
+         Y_Index: Index;
+         Initial_Z_Index: Index := Y.Length + 1; 
+         Initial_Y_Index: Index := 1;            
+         Cycle_Length: Index := 1;               
+      begin
+         Carry := 0;
+         for I in 1 .. ZL loop
+            Digit := Carry;
+            Carry := 0;
+            Z_Index := Initial_Z_Index;          
+            Y_Index := Initial_Y_Index;          
+            for J in 1 .. Cycle_Length loop
+               if Digit > Base_Squared then
+                  Digit := Digit - Base_Squared;
+                  Carry := Carry + Base;
+               end if;
+               Digit := Digit + Double(Z.D(Z_Index))
+                              * Double(Y.D(Y_Index));
+               Z_Index := Z_Index + 1;
+               Y_Index := Y_Index - 1;
+            end loop;
+            -- now adjust limits and length of cycle
+            if I < Y.Length then
+               Cycle_Length := Cycle_Length + 1;
+               Initial_Y_Index := Initial_Y_Index + 1;
+            else
+               Initial_Z_Index := Initial_Z_Index + 1;
+            end if;
+            if I < X.Length then
+               Cycle_Length := Cycle_Length + 1;
+            end if;
+            Cycle_Length := Cycle_length - 1;
+            Carry := Carry + Digit / Base;
+            Z.D(I) := Single(Digit mod Base);
+         end loop;
+      end;
+      if Carry /= 0 then -- one more digit in answer
+         ZL := ZL + 1;
+         Z.D(ZL) := Single(Carry);
+      end if;
+      Z.Length := ZL;
+      Z.Sign := X.Sign * Y.Sign;
+   end Mul;
+
+   procedure Div(X, Y: Number; Quotient,
+                               Remainder: out Number) is
+      U: Number renames Quotient;
+      V: Number renames Remainder;
+      Digit, Scale, Carry: Double;
+      U0, U1, U2: Double;
+      V1, V2: Double;
+      QD: Double;
+      LOQ: constant Index := Y.Length;
+      HIQ: constant Index := X.Length;
+      QL: Index;
+      RL: Index;
+      QStart: Index;
+      J : Index;
+   begin
+      if Y.Length = 0 then
+         raise Number_Error;
+      end if;
+      if Quotient.Max_Digits < X.Length or
+         Remainder.Max_Digits < Y.Length then
+         raise Number_Error;
+      end if;
+      if X.Length < Y.Length then  -- Quotient is definitely zero
+         Quotient.Sign := +1;
+         Quotient.Length := 0;
+         Remainder.Sign := X.Sign;
+         Remainder.Length := X.Length;
+         for I in 1 .. X.Length loop
+            Remainder.D(I) := X.D(I);
+         end loop;
+         return;
+      end if;
+
+      QL := X.Length - Y.Length + 1;
+      RL := Y.Length;
+      QStart := QL;
+
+      -- compute normalizing factor
+      Scale := Base/Double(Y.D(Y.Length)+1);
+
+      -- scale X and copy to U
+      Carry := 0;
+      for I in 1 .. X.Length loop
+         Digit := Double(X.D(I)) * Scale + Carry;
+         Carry := Digit / Base;
+         U.D(I) := Single(Digit mod Base);
+      end loop;
+      U0 := Carry;  -- leading digit of dividend
+
+      -- scale Y and copy to V
+      Carry := 0;
+      for I in 1 .. Y.Length loop
+         Digit := Double(Y.D(I)) * Scale + Carry;
+         Carry := Digit / Base;
+         V.D(I) := Single(Digit mod Base);
+      end loop;
+      -- no further carry
+
+      -- set V1 and V2 to first two digits of divisor
+      V1 := Double(V.D(Y.Length));
+      if Y.Length > 1 then
+         V2 := Double(V.D(Y.Length-1));
+      else
+         V2 := 0;
+      end if;
+
+      -- now iterate over digits in answer
+      -- with U0, U1 and U2 being first three digits of dividend
+      for I in reverse LOQ .. HIQ loop
+         U1 := Double(U.D(I));
+         if Y.Length > 1 then
+            U2 := Double(U.D(I-1));
+         else
+            U2 := 0;
+         end if;
+
+         -- now set initial estimate of digit in quotient
+         if U0 = V1 then
+            QD := Base - 1;
+         else
+            QD := (U0 * Base + U1) / V1;
+         end if;
+
+         -- now refine estimate by considering U2 also
+         while V2*QD > (U0*Base+U1-QD*V1)*Base + U2 loop
+            QD := QD - 1;
+         end loop;
+
+         -- QD is now correct digit or possibly one too big
+         -- subtract QD times V from U
+         Carry := 0;
+         J := QStart;
+         for I in 1 .. Y.Length loop
+            Digit := Double(U.D(J)) - Carry
+                        - QD * Double(V.D(I));
+            if Digit < 0 then
+               Carry := (-1-Digit) / Base + 1;
+               Digit := Digit + Carry * Base;
+            else
+               Carry := 0;
+            end if;
+            U.D(J) := Single(Digit);
+            J := J + 1;
+         end loop;
+         if Carry > U0 then -- estimate was too large
+            declare
+               Carry, Digit: Single;
+            begin
+               QD := QD - 1;
+               Carry := 0;  J := QStart;
+               for I in 1 .. Y.Length loop
+                  Digit := U.D(J) + Carry + V.D(I);
+                  if Digit >= Base then
+                     Carry := 1;
+                     Digit := Digit - Base;
+                  else
+                     Carry := 0;
+                  end if;
+                  U.D(J) := Digit;
+                  J := J + 1;
+               end loop;
+            end;
+         end if;
+         -- QD is now the required digit
+         U0 := Double(U.D(I));  U.D(I) := Single(QD);
+         QStart := QStart - 1;
+      end loop;
+      
+      -- delete possible leading zero in quotient
+      if U.D(HIQ) = 0 then
+         QL := QL - 1;
+      end if;
+   
+      -- copy remainder into place and scale 
+      -- top digit is in U0 still
+      Digit := U0;
+      for I in reverse 2 .. RL loop
+         Remainder.D(I) := Single(Digit/Scale);
+         Carry := Digit mod Scale;
+         Digit := Double(U.D(I-1)) + Carry * Base;
+      end loop;
+      Remainder.D(1) := Single(Digit/Scale);
+      -- delete leading zeroes in remainder
+      while RL > 0 and then Remainder.D(RL) = 0 loop
+         RL := RL - 1;
+      end loop;
+      Remainder.Length := RL;
+      if Remainder.Length = 0 then
+         Remainder.Sign := +1;
+      else
+         Remainder.Sign := X.Sign;
+      end if;
+
+      -- slide quotient into place
+      -- Quotient.D(1 .. QL) := U.D(LOQ .. HIQ);
+      for I in 1 .. QL loop
+         Quotient.D(I) := U.D(I + LOQ - 1);
+      end loop;
+      Quotient.Length := QL;
+      if Quotient.Length = 0 then
+         Quotient.Sign := +1;
+      else
+         Quotient.Sign := X.Sign * Y.Sign;
+      end if;
+
+   end Div;
+
+   procedure To_Number(S: String; N: out Number) is
+      NL: Index := 0;
+      Place: Integer := 0;
+      Is_A_Number: Boolean := False;
+      Digit: Single := 0;
+      Ch: Character;
+      Last_I: Positive;
+      Dig_Of: constant array (Character range '0' .. '9') of
+                   Single := (0, 1, 2, 3, 4, 5, 6, 7, 8, 9);
+   begin
+      N.Sign := +1;   -- set default sign
+      -- scan string from end
+      for I in reverse S'Range loop
+         Last_I := I;   -- note how far we have got
+         Ch := S(I);
+         case Ch is
+            when '0' .. '9' =>
+               -- add digit to number so far
+               if Place = 0 then
+                  NL := NL + 1;
+                  if NL > N.Max_Digits then
+                     raise Number_Error;
+                  end if;
+               end if;
+               Digit := Digit + Dig_Of(Ch) * 10**Place;
+               Place := Place + 1;
+               if Place = Base_Exp then
+                  N.D(NL) := Digit;
+                  Digit := 0;
+                  Place := 0;
+               end if;
+               Is_A_Number := True;
+            when '_' =>
+               -- underscore must be embedded in digits
+               if not Is_A_Number then
+                  raise Number_Error;
+               end if;
+               Is_A_Number := False;
+            when '+' | '-' | ' ' =>
+               -- lump so far must be a valid number
+               if not Is_A_Number then
+                  raise Number_Error;
+               end if;
+               if Ch ='-' then N.Sign := -1; end if;
+               exit;         -- leave loop
+            when others =>
+               raise Number_Error;
+         end case;
+      end loop;
+      -- check we had a number
+      if not Is_A_Number then
+         raise Number_Error;
+      end if;
+      -- add the last digit if necessary
+      if Place /= 0 then
+         N.D(NL) := Digit;
+      end if;
+      -- check that any other characters are leading spaces
+      for I in S'First .. Last_I - 1 loop
+         if S(I) /= ' ' then
+            raise Number_Error;
+         end if;
+      end loop;
+      -- remove leading zeroes if any, beware zero case
+      while NL > 0 and then N.D(NL) = 0 loop
+         NL := NL - 1;
+      end loop;
+      N.Length := NL;
+   end To_Number;
+
+   procedure To_Number(I: Integer; N: out Number) is
+      NL: Index := 0;
+      II: Integer;
+   begin
+      if I = 0 then
+         N.Sign := +1;  N.Length := 0;
+         return;
+      end if;
+      if I > 0 then
+         II := I;  N.Sign := +1;
+      else
+         II := -I;  N.Sign := -1;
+      end if;
+      while II /= 0 loop
+         NL := NL + 1;
+         if NL > N.Max_Digits then
+            raise Number_Error;
+         end if;
+         N.D(NL) := Single(II mod Base);
+         II := II / Base;
+      end loop;
+      N.Length := NL;
+   end To_Number;
+
+   procedure To_Text(N: Number; S: out String) is
+      SI: Natural := S'Last;
+      Digit: Single;
+      Char_Of: constant array (Single range 0 .. 9) of
+                              Character := "0123456789";
+   begin
+      if N.Length = 0 then   -- zero case
+         if SI < 2 then
+            raise Number_Error;
+         end if;
+         S(SI) := Char_Of(0);
+         S(SI-1) := '+';
+         for I in 1 .. SI-2 loop
+            S(I) := ' ';
+         end loop;
+         return;
+      end if;
+      if SI < Base_Exp * N.Length + 1 then
+         raise Number_Error;
+      end if;
+      for I in 1 .. N.Length loop
+         Digit := N.D(I);
+         for J in 1 .. Base_Exp loop
+            S(SI) := Char_Of(Digit mod 10);
+            Digit := Digit / 10;
+            SI := SI - 1;
+         end loop;
+      end loop;
+      while S(SI + 1) = '0' loop
+         SI := SI + 1;  -- delete leading zeroes
+      end loop;
+      if N.Sign = +1 then
+         S(SI) := '+';
+      else
+         S(SI) := '-';
+      end if;
+      for I in 1 .. SI - 1 loop
+         S(I) := ' ';
+      end loop;
+   end To_Text;
+
+   procedure To_Integer(N: Number; I: out Integer) is
+      II: Integer := 0;
+   begin
+      for I in reverse 1 .. N.Length loop
+         II := II * Base + Integer(N.D(I));
+      end loop;
+      if N.Sign = -1 then II := -II; end if;
+      I := II;
+   end To_Integer;
+
+end Numbers.Proc;
+
+--------------------------------------------------------------------
+
+with Ada.Text_IO; use Ada;
+with Numbers.Proc;
+package body Numbers.IO is
+
+   use Proc;
+
+   procedure Put(Item: Number;
+                 Width: Natural := Default_Width) is
+      Block_Size: constant := 3;
+      Places: Integer := No_Of_Places(Item);
+      S: String(1 .. Places + 1);
+      SP: Positive := 1;
+      Before_Break: Integer;
+   begin
+      To_Text(Item, S);
+      -- allow for leading spaces in S
+      while S(SP) = ' ' loop
+         SP := SP + 1;  Places := Places - 1;
+      end loop;
+      -- now output leading spaces for padding if any
+      for I in 1 .. Width -
+            (Places + 1 + (Places - 1) / Block_Size) loop
+         Text_IO.Put(' ');
+      end loop;
+      if S(SP) = '+' then S(SP) := ' '; end if;
+      Text_IO.Put(S(SP));  -- output minus or space
+      -- output digits with underscores every "Blocksize"
+      Before_Break := (Places - 1) rem Block_Size + 1;
+      for I in SP + 1 .. S'Last loop
+         if Before_Break = 0 then
+            Text_IO.Put('_');
+            Before_Break := Block_Size;
+         end if;
+         Text_IO.Put(S(I));
+         Before_Break := Before_Break - 1;
+      end loop;
+   end Put;
+
+   procedure Get(Item: out Number) is
+      -- declare string large enough to hold maximum value
+      -- allows every other character to be an underscore!
+      S: String(1 .. Base_Exp * Max_Index * 2);
+      SP: Positive := 1;
+      Places: Integer := 0;
+      Ch: Character;
+      EOL: Boolean;       -- end of line
+   begin
+      -- loop for first digit or sign, skipping spaces
+      loop
+         Text_IO.Get(Ch);
+         case Ch is
+            when ' ' =>
+               null;
+            when '+'| '-' =>
+               S(SP) := Ch;  SP := SP + 1;
+               exit;
+            when '0' .. '9' =>
+               S(SP) := Ch;  SP := SP + 1; Places := 1;
+               exit;
+            when others =>
+               raise Number_Error;
+         end case;
+      end loop;
+      -- now accept only digits and underscores
+      -- count the digits in Places
+      -- stop on end of line or other character
+      loop
+         Text_IO.Look_Ahead(Ch, EOL);
+         exit when EOL;
+         Text_IO.Get(Ch);
+         case Ch is
+            when '0' .. '9' =>
+               S(SP) := Ch;  SP := SP + 1;
+               Places := Places + 1;
+            when '_' =>
+               S(SP) := Ch;  SP := SP + 1;
+            when others =>
+               exit;
+         end case;
+      end loop;
+      -- now declare a Number big enough
+      -- note Item assumed unconstrained
+      declare
+         Result: Number((Places - 1)/Base_Exp + 1);
+      begin
+         To_Number(S(1 .. SP - 1), Result);
+         Item := Result;
+      end;
+   end Get;
+
+end Numbers.IO;
+
+--------------------------------------------------------------------
+
+with Numbers.Proc;
+package body Numbers.Func is
+
+   use Proc;
+
+   function "+" (X: Number) return Number is
+   begin
+      return X;
+   end "+";
+
+   function "-" (X: Number) return Number is
+      N: Number(X.Max_Digits);
+   begin
+      N := X;  Neg(N);
+      return N;
+   end "-";
+
+   function "abs" (X: Number) return Number is
+   begin
+      if X < Zero then return -X; else return X; end if;
+   end "abs";
+
+   function "+" (X, Y: Number) return Number is
+      Z: Number(Index'Max(Length(X), Length(Y)) + 1);
+   begin
+      Add(X, Y, Z);
+      return Z;
+   end "+";
+
+   function "-" (X, Y: Number) return Number is
+      Z: Number(Index'Max(Length(X), Length(Y)) + 1);
+   begin
+      Sub(X, Y, Z);
+      return Z;
+   end "-";
+
+   function "*" (X, Y: Number) return Number is
+      Z: Number(Length(X) + Length(Y));
+   begin
+      Mul(X, Y, Z);
+      return Z;
+   end "*";
+
+   function "/" (X, Y: Number) return Number is
+      Q: Number(Length(X));
+      R: Number(Length(Y));
+   begin
+      Div(X, Y, Q, R);
+      return Q;
+   end "/";
+
+   function "rem" (X, Y: Number) return Number is
+      Q: Number(Length(X));
+      R: Number(Length(Y));
+   begin
+      Div(X, Y, Q, R);
+      return R;
+   end "rem";
+
+   function "mod" (X, Y: Number) return Number is
+      Q: Number(Length(X));
+      R: Number(Length(Y));
+   begin
+      Div(X, Y, Q, R);
+      if (X < Zero and Y > Zero) or (X > Zero and Y < Zero) then
+         R := R + Y; 
+      end if;
+      return R;
+   end "mod";
+
+   function "**" (X: Number; N: Natural) return Number is
+      Result: Number := One;
+      Term: Number := X;
+      M: Natural := N;
+   begin
+      loop
+         if M rem 2 /= 0 then
+            Result := Term * Result;
+         end if;
+         M := M / 2;
+         exit when M = 0;
+         Term := Term * Term;
+      end loop;
+      return Result;
+   end "**";
+
+   function "<" (X, Y: Number) return Boolean is
+   begin
+      return Compare(X, Y) < 0;
+   end "<";
+
+   function "<=" (X, Y: Number) return Boolean is
+   begin
+      return Compare(X, Y) <= 0;
+   end "<=";
+
+   function ">" (X, Y: Number) return Boolean is
+   begin
+      return Compare(X, Y) > 0;
+   end ">";
+
+   function ">=" (X, Y: Number) return Boolean is
+   begin
+      return Compare(X, Y) >= 0;
+   end ">=";
+
+   function "=" (X, Y: Number) return Boolean is
+   begin
+      return Compare(X, Y) = 0;
+   end "=";
+
+   function Make_Number(S: String) return Number is
+      Result: Number((S'Length - 1) / Base_Exp + 1);
+   begin
+      To_Number(S, Result);
+      return Result;
+   end Make_Number;
+
+   function Make_Number(I: Integer) return Number is
+      Base_Digits: Index := 0;
+      II: Integer := abs I;
+      Base: constant := 10 ** Base_Exp;
+   begin
+      -- loop to determine discriminant for result
+      while II /= 0 loop
+         Base_Digits := Base_Digits + 1;
+         II := II / Base;
+      end loop;
+      declare
+         Result: Number(Base_Digits);
+      begin
+         To_Number(I, Result);
+         return Result;
+      end;
+   end Make_Number;
+
+   function String_Of(N: Number) return String is
+      Places: Integer := No_Of_Places(N);
+      S: String(1 .. Places + 1);
+      SP: Positive := 1;
+   begin
+      To_Text(N, S);
+      -- allow for leading spaces in S
+      while S(SP) = ' ' loop
+         SP := SP + 1;
+      end loop;
+      return S(SP .. S'Last);
+   end String_Of;
+
+   function Integer_Of(N: Number) return Integer is
+      Result: Integer;
+   begin
+      To_Integer(N, Result);
+      return Result;
+   end Integer_Of;
+      
+end Numbers.Func;
+
+
+--- Numbers calculator ---------------------------------------------
+
+--with ID, Numbers.Func; use Numbers.Func;
+--package Numbers_ID is new ID(Number, Zero, One);
+
+--with ID.IO, Numbers.IO; use Numbers.IO;
+--package Numbers_ID.The_IO is new Numbers_ID.IO;
+
+--with Numbers_ID.The_IO;
+--with Numbers.Func;
+--with Calculator;
+--procedure Numcalc is
+  -- new Calculator.Run(Numbers_ID,
+                     -- Numbers_ID.The_IO,
+                     -- Numbers.Func."+");
+
+
+--- test 1 - powers of 11 and 99 -----------------------------------
+
+with Numbers.Proc; use Numbers.Proc;
+with Ada.Text_IO; use Ada;
+procedure Test_11_99 is
+
+   U: Number(2);
+
+   procedure P(X: Number) is
+      S: String(1 .. 150);
+   begin
+      To_Text(X, S);
+      Text_IO.Put(S); Text_IO.New_Line;
+   end P;
+
+   procedure Power(U: Number; V: Integer) is
+      W: Number(50);
+   begin
+      To_Number(1, W);
+      for I in 1 .. V loop
+         Mul(W, U, W);
+         P(W);
+      end loop;
+   end Power;
+
+begin
+   To_Number(11, U);
+   Power(U, 50);
+   To_Number(99, U);
+   Power(U, 50);
+end;
+
+-- test 2 - Mersenne using procedural forms ------------------------
+
+with Ada.Calendar; use Ada.Calendar;
+with Numbers.Proc; use Numbers.Proc;
+with Ada.Text_IO, Ada.Integer_Text_IO;
+use Ada.Text_IO, Ada.Integer_Text_IO;
+procedure Test2 is
+
+   Nop: constant := 30;
+
+   Loop_Start, Loop_End: Integer;
+   T_Start, T_End: Time;
+
+   Is_Prime: Boolean;
+   MM: Number(50);
+
+   Primes: array(1 .. Nop)  of Integer := 
+   (3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,
+    61,67,71,73,79,83,89,97,101,103,107,109,113,127);
+
+   package Duration_IO is new Fixed_IO(Duration);
+   use Duration_IO;
+
+
+procedure Lucas_Lehmer (P: Integer;
+                        Mersenne: out Number;
+                        Is_Prime: out Boolean) is
+
+   Two: Number(1);
+   M: Number(Mersenne.Max_Digits); 
+   L, W, Quotient: Number(M.Max_Digits*2);
+begin
+   To_Number(2, Two);
+   To_Number(4, L);
+   To_Number(1, M);
+   for I in 1 .. P loop
+      Mul(M, Two, M);
+   end loop;
+   Sub(M, One, M);
+   for I in 1 .. P-2 loop
+      Mul(L, L, W);
+      Sub(W, Two, W);
+      Div(W, M, Quotient, L);
+    --  L := (L**2 - Two) mod M;
+   end loop;
+   Is_Prime := Compare(L, Zero) = 0;
+   Mersenne := M;
+end Lucas_Lehmer;
+
+   procedure Put(X: Number) is
+      S: String(1 .. 45);
+   begin
+      To_Text(X, S);
+      Put(S); 
+   end Put;
+
+begin
+   Put_Line("Start loop? ");  Get(Loop_Start);
+   Put_Line("End loop? ");  Get(Loop_End);
+
+   Put_Line("Mersenne Primes");
+   Put_Line("    Time     P                           2**P-1");
+   for I in Loop_Start .. Loop_End loop
+      New_Line;
+      T_Start := Clock;
+      Lucas_Lehmer(Primes(I), MM, Is_Prime);
+      T_End := Clock;
+      New_Line;  Put(T_End - T_Start, 2, 1);
+      Put(Primes(I), 4);  Put(" : ");
+      Put(MM);
+      if Is_Prime then
+         Put(" is prime");
+      else
+         Put(" is not prime");
+      end if;
+   end loop;
+end Test2;
+
+
+--- test 5 - Mersenne using functional forms -----------------------
+
+with Ada.Calendar; use Ada.Calendar;
+with Numbers.Func; use Numbers.Func;
+with Numbers.IO; use Numbers.IO;
+with Ada.Text_IO, Ada.Integer_Text_IO;
+use Ada.Text_IO, Ada.Integer_Text_IO;
+procedure Test5 is
+
+   Nop: constant := 30;
+
+   Loop_Start, Loop_End: Integer;
+   T_Start, T_End: Time;
+
+   Is_Prime: Boolean;
+   LL, MM: Number;
+
+   Primes: array(1 .. Nop)  of Integer := 
+   (3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,
+    61,67,71,73,79,83,89,97,101,103,107,109,113,127);
+
+   package Duration_IO is new Fixed_IO(Duration);
+   use Duration_IO;
+
+
+procedure Lucas_Lehmer (Q: Integer;
+                        Mersenne: out Number;
+						Lout: out Number;
+                        Is_Prime: out Boolean) is
+
+   Two: constant Number := Make_Number(2);
+   M: constant Number := Two**Q - One; 
+   L: Number := Make_Number(4);
+begin
+   for I in 1 .. Q-2 loop
+      L := (L**2 - Two); -- mod M; -- mod M here is optional;
+	  Put(L);  New_line(2);
+   end loop;
+   Is_Prime := L mod M = Zero;
+   Lout := l;
+   Mersenne := M;
+end Lucas_Lehmer;
+
+begin
+   Put_Line("Start loop? ");  Get(Loop_Start);
+   Put_Line("End loop? ");  Get(Loop_End);
+
+   Put_Line("Mersenne Primes");
+   Put_Line("   Time      P                          2**P-1");
+   for I in Loop_Start .. Loop_End loop
+      New_Line;
+      T_Start := Clock;
+      Lucas_Lehmer(Primes(I), MM, LL, Is_Prime);
+      T_End := Clock;
+      New_Line;  Put(T_End - T_Start, 2, 1);
+      Put(Primes(I), 4);  Put(" : ");
+      Put(MM, 20);
+      if Is_Prime then
+         Put(" is prime");
+		New_line;
+		 put((MM+One)/Make_number(2)*MM, 30);  Put(" is perfect");
+		 New_Line;
+      else
+         Put(" is not prime");
+      end if;
+      New_line;
+	  -- comment out next four lines to avoid detail
+	  Put("L is "); Put(LL); Put(" equals "); Put(MM); Put(" times "); Put(LL/MM);
+	  if not Is_prime then
+	  	New_Line; Put(" remainder = "); Put(LL mod MM);
+	  end if;
+	-- end of comment
+   end loop;
+   skip_Line(2);
+end Test5;
+
+--- test 6 - GCD and Mersenne --------------------------------------
+
+with Ada.Text_IO, Ada.Integer_Text_IO;
+use Ada.Text_IO, Ada.Integer_Text_IO;
+with Numbers.Proc; use Numbers.Proc;
+procedure Test6 is
+   subtype CNumber is Number(50);
+   M1, M2, M3: CNumber;
+   P1, P2, P3: Integer;
+   G, H: CNumber;
+   XX, YY, QQ, ZZ: CNumber;
+   Two: CNumber;
+   N: Cnumber; 
+   Start, Stop: Integer;
+
+procedure GCD(X, Y: CNumber; Z: out CNumber) is
+begin
+   XX := X;  YY := Y;
+   while Compare(YY, Zero) /= 0 loop
+      Div(XX, YY, QQ, ZZ);
+      XX := YY;
+      YY := ZZ;
+   end loop;
+   Z := XX;
+end GCD;
+
+procedure Mersenne(P: Integer; M:  in out CNumber) is
+begin
+   To_Number(2, Two);
+   To_Number(1, N);
+   for I in 1 .. P loop
+      Mul(N, Two, N);
+   end loop;
+   Sub(N, One, N);
+   M := N;
+end Mersenne;
+
+begin
+   To_Number(2, Two);
+   Put_Line("Start? "); Get(Start);
+   Put_Line("Stop?  "); Get(Stop);
+   for I in Start .. Stop loop
+      P1 := 2*I + 1;
+      Mersenne(P1, M1);
+      for J in 1 .. I -1 loop
+         P2 := 2*J + 1;
+         Mersenne(P2, M2);
+         GCD(M1, M2, G);
+         for K in 1 .. I loop
+            P3 := 2*K + 1;
+            Mersenne(P3, M3);
+            if Compare(M3, G) > 0 then
+               exit;
+            end if;
+            if Compare(M3, G) =0 then
+               New_Line;
+               Put(P1); Put(P2); Put(P3);
+               exit;
+            end if;
+         end loop;
+      end loop;
+   end loop;
+end Test6;
+
+------- test 7 multipication
+
+with Ada.Text_IO, Ada.Integer_Text_IO;
+use Ada.Text_IO, Ada.Integer_Text_IO;
+with Numbers.Proc; use Numbers.Proc;
+with Numbers.IO; use Numbers.IO;
+procedure Test7 is
+	X: Number;
+	Y: Number;
+ 	Z: Number(50);
+begin
+	Put_line("Multiplier test");
+	Put("X = ");  Get(X);
+	Put("Y = ");  Get(Y);
+	Mul(X, Y, Z);
+	Put_Line("product is ");
+	Put(Z);
+	New_Line(2);
+end Test7;
+  
+****************************************************************
+
+From: Jeff Cousins
+Sent: Wednesday, January 31, 2018  6:17 AM
+
+[This some additional test programs for John Barnes' Bignum package; the
+Numbers package was duplicated at the front, which I removed. - Editor.]
+
+
+--package Monitor is
+--	C1, C2, C3, C4: Integer;
+--end;
+
+package Primes is
+	Max: constant := 20000;
+	Prime: array (1..Max) of Integer;
+	pragma Elaborate_Body;
+end;
+
+
+package body Primes is
+	N: Integer:= 2;
+	Index: Integer := Prime'First;
+	Found: Boolean := False;
+begin
+	-- initialization part to build prime table
+	Prime(Prime'First) := 2;
+	loop
+		N := N+1; Found := True;
+		for I in Prime'First .. Index loop
+			if n rem Prime(I) = 0 then -- divides by existing prime
+				Found := False; exit;
+			end if;
+		end loop;
+		if Found then
+
+			-- found a new prime
+			Index := Index+1; Prime(Index) := n;
+			exit when Index = Max;
+		end if;
+	end loop;
+end Primes;
+
+Package Squares is
+	Square_Digits: Integer := 4;
+	Square_Ten_Power: Integer := 10**Square_Digits;
+	Poss_Last_Digits: array(0..Square_Ten_Power-1) of Boolean := (others => False);
+	pragma Elaborate_Body;
+end;
+
+package body Squares is
+begin
+	
+	-- make Poss_Last_Digits array
+	for I in 1 .. Square_Ten_Power loop
+		Poss_Last_Digits(I*I mod Square_Ten_Power) := True;
+	end loop;
+end Squares;		
+	
+
+
+with Numbers.Func; use Numbers.Func;
+-- with Monitor;
+procedure Square_root(XX: in Number; Try: Number; X: out Number; R: out Number) is
+	-- X is the largest X such that X*X + R equals XX with R >= zero
+	-- Try is initial guess 
+	K, KK: Number;
+	DK: Number;
+	Three: Number := Make_Number(3);
+begin
+	k := Try;
+	loop
+		KK := K*K;
+		DK := (XX - KK)/(K+K);
+		-- iterate until nearly there
+		if abs DK < Three then -- nearly there
+			if XX < KK then
+				K := XX/K;		-- ensure K*K is less than XX
+			end if; 
+			exit;
+		end if;
+		-- do another iterate
+		K := K+DK;
+	--	Monitor.c2 := Monitor.C2 + 1;
+	end loop;
+
+	-- now loop from below
+
+	loop
+		KK := K*K;
+	--	Monitor.C2 := Monitor.C2 + 1;
+		if KK >= XX then 
+			if KK = XX then
+				X := K; R := Zero; return;
+			end if;
+			X := K - One; R := XX - X*X; return;
+		end if;
+		K := K+One;
+	end loop;
+
+
+end Square_Root;
+
+
+with Square_Root;
+with Numbers.Func; use Numbers.Func;
+-- with Monitor;
+with Squares; use Squares;
+procedure Fermat_Long(N: in Number; Min_Prime: in Number; P, Q: out Number) is
+	-- we know that factors up to Min_Prime have been removed, so max square to try is
+	-- roughly (N/Min_Prime + Min_Prime)/2; we add 2 for luck
+	Two: Number := Make_Number(2);
+	Number_Square_Digits: Number := Make_Number(Square_Digits);
+	Number_Square_Ten_Power: Number := Make_Number(Square_Ten_Power);
+	X: Number;
+	Y: Number;
+	R: Number;
+	K: Number;
+	DK: Number;
+	N2, X2, K2: Integer;
+	Last_Digits: Integer;
+	Try: Number := One;
+	Max_square : Number := (N/Min_Prime + Min_Prime)/ Two + Two;
+
+	
+begin
+
+	Square_Root(N, One, X, R);
+	if R = Zero then
+		-- N was a perfect square
+		P := X; Q := X;
+		return;
+	end if;
+	--	Monitor.C1 := Monitor.C2;
+	--	Monitor.C2 := 0;
+
+		K := X*X-n;
+		DK := X+X+One;
+		N2 := Integer_Of(N rem Number_Square_Ten_Power);
+		X2 := Integer_Of(X rem Number_Square_Ten_Power);
+
+	loop
+	--	Monitor.C3 := Monitor.C3 + 1;
+		X := X + One;
+		if X > Max_Square then
+			-- must be prime
+			P := N; Q := One;
+			return;
+		end if;
+
+		X2 := (X2 + 1) rem Square_Ten_Power;
+		K := K + DK;	-- omit if DK not used
+		DK := DK + Two;
+	--	K := X*X-N;  -- omit if DK used
+		K2 := (X2*X2-N2) mod Square_Ten_Power;
+	--	Last_Digits := Integer_Of(K rem Number_Square_Ten_Power);
+		Last_Digits := K2;
+		if Poss_Last_Digits(Last_Digits) then
+	--		Monitor.C4 := Monitor.C4 + 1;
+
+			Square_Root(K, Try, Y, R);
+			if R = Zero then
+				-- X*X-N was a perfect square
+				P := X+Y; Q := X-Y;
+				return;
+			end if;
+			Try := Y;
+		end if;
+	end loop;
+exception
+	when others => p := Zero; Q := Zero;
+end Fermat_long;
+
+with Numbers.IO; use Numbers.IO;
+with Numbers.Func; use Numbers.Func;
+with Ada.Text_IO; use Ada.Text_IO;
+with Ada.Integer_Text_IO; use Ada.Integer_Text_IO;
+-- with Monitor;
+with Primes; use Primes;
+with Fermat_Long;
+with Ada.Calendar; use Ada.Calendar;
+procedure Do_Fermat_Long is
+	NN, PP, QQ: Number;
+	T_Start, T_End: Time;
+	package Duration_IO is new Fixed_IO(Duration);
+	use Duration_IO;
+begin	
+	Put_Line("Welcome to Fermat's method (multilength)");
+	loop
+<<Again>>
+		Begin	
+			New_Line;
+			Put("Insert number N = "); Get(NN);
+		exception
+			when others =>
+				Put_Line("Not a number or too big");  Skip_Line;
+				goto Again;
+		end;
+
+		exit when NN = Zero;
+
+			-- check to see if N divides by a known prime
+
+		for i in Prime'Range loop
+			loop
+				if NN rem make_Number(Prime(I)) = Zero then
+					Put("N divides by "); Put(Prime(I), 0); Put_Line(" so removing factor");
+					NN := NN / Make_Number(Prime(I));
+				else
+					exit;
+				end if;
+			end loop;
+		end loop;
+
+
+		if nN rem Make_number(4) = Make_Number(2) then
+			Put_Line("Algorithm fails on odd multiples of 2, so halving N"); 
+			NN := NN/Make_Number(2);
+			Put("New N = "); Put(NN, 0); New_Line;
+		end if;
+		if NN = One then
+			Put_Line("all factors removed");
+		else
+			
+	--	Monitor.C1 := 0;
+	--	Monitor.C2 := 0;
+	--	Monitor.C3 := 0;
+	--	Monitor.C4 := 0;
+		T_Start := Clock;
+		Fermat_Long(NN, Make_Number(Prime(Max)), PP, QQ);
+		T_End := Clock;
+		if PP = Zero and QQ = zero then
+			Put_Line("Failed internally");
+			goto Again;
+		end if;
+		Put("Two factors are   "); Put(PP, 0); Put("   "); Put(QQ, 0); New_Line;
+		Put(T_End-T_Start, 2, 1);  New_Line;
+	 --	Put(Monitor.C1, 9);  Put(Monitor.C2, 9); Put(Monitor.C3, 9); Put(Monitor.C4, 9); New_Line;
+	end if;
+	end loop;
+	Put_line("Goodbye");  Skip_Line(2);
+end Do_Fermat_Long;
+
+****************************************************************
+
+From: Randy Brukardt
+Sent: Wednesday, January 31, 2018  5:03 PM
+
+>Sending in three parts as the original message was too big for the mail
+>server.
+
+Just so the rest of you know, "Part 3" was an executable Windows program which
+we've decided not to distribute at all (it's too big for the list, and
+executables for specific computers seem out-of-bounds for this list anyway,
+given all of us know how to operate our favorite Ada compiler and probably
+several non-favorite compilers as well).
+
+So don't look for the non-existent part 3.
+
+****************************************************************

Questions? Ask the ACAA Technical Agent