--  Program 4   May 2022



--  Super Sieve 

-- This program is identical to that in Programming in Ada 2012
-- except for the added comments and layout of use and with clauses.


generic

   -- The type Element represents the algebra concerned (such as
   -- type Integer for normal prime numbers).  Unit is the unit
   -- value of the type (such as 1 for type Integer).  The function
   -- Succ returns the successor of the given value such that every
   -- value follows at least one of its prime divisors (in the case
   -- of the type Integer the usual Succ function is appropriate). 
   -- The function Is_Factor checks whether the value E is exactly
   -- divisible by P.  The procedure Put outputs a value of the
   -- type in a simple format.
   --    The exported procedure Do_It tries a sequence of values of
   -- the type Element starting at Unit and as given by calls of
   -- Succ and checks to see whether they are prime.  Each prime
   -- number is output via the child package Eratosthenes.Frame.
   --    A total of To_Try values are tried and the number of primes
   -- found is returned in Found.

   type Element is private;
   Unit: in Element;
   with function Succ(E: Element) return Element;
   with function Is_Factor(E, P: Element) return Boolean;
   with procedure Put(E: in Element);
package Eratosthenes is
   procedure Do_It(To_Try: in Integer; Found: out Integer);
end;

private generic
package Eratosthenes.Frame is

   -- This interface is as described in Chapter 20.  In this
   -- implementation the procedures Clear_Frame and Write_To_Frame
   -- do nothing and Make_Frame simply outputs the value on 
   -- a new line.

   type Position is private;
   procedure Make_Frame(Prime: in Element; Where: out Position);
   procedure Write_To_Frame(Value: in Element; Where: in Position);
   procedure Clear_Frame(Where: in Position);
private
   type Position is null record;  -- a null type
end;

with Ada.Text_IO;
package body Eratosthenes.Frame is

   procedure Make_Frame(Prime: in Element; Where: out Position) is
   begin
      Ada.Text_IO.New_Line;
      Put(Prime);
      Where := (null record);   -- null aggregate
   end Make_Frame;

   procedure Write_To_Frame(Value: in Element; Where: in Position) is
   begin null; end;

   procedure Clear_Frame(Where: in Position) is
   begin null; end;

end Eratosthenes.Frame;

with Eratosthenes.Frame;
package body Eratosthenes is
   package Inner_Frame is new Frame;
   use Inner_Frame;

   Primes_Found: Integer := 0;

   protected Finished is
      entry Wait;
      procedure Signal;
   private
      Occurred: Boolean := False;
   end;

   protected body Finished is
      entry Wait when Occurred is
      begin
         null;
      end Wait;

      procedure Signal is
      begin
         Occurred := True;
      end Signal;
   end Finished;

   task type Filter is
   
      -- The first call of Input gives the new trial divisor
      -- which is a newly found prime number.  Succeeding calls 
      -- give the numbers to be tested.
      --    The entry Stop is called when the task is to close down.

      entry Input(Number: in Element);
      entry Stop;
   end Filter;

   function Make_Filter return access Filter is
   begin
      return new Filter;
   end Make_Filter;

   task body Filter is
      P: Element;      -- prime divisor
      N: Element;      -- trial element
      Here: Position;
      Next: access Filter;
   begin
      accept Input(Number: in Element) do
         P := Number;
      end;

      -- The call of Make_Frame causes the new prime number
      -- to be printed.

      Make_Frame(P, Here);
      Primes_Found := Primes_Found + 1;
      loop
         select
            accept Input(Number: in Element) do
               N := Number;
            end;
         or
            accept Stop;
            if Next = null then
               Finished.Signal;
            else
               Next.Stop;
            end if;
            exit;
         end select;
         Write_To_Frame(N, Here);
         if not Is_Factor(N, P) then
            if Next = null then
               Next := Make_Filter;
            end if;
            Next.Input(N);
         end if;
         Clear_Frame(Here);
      end loop;
   end Filter;

   procedure Do_It(To_Try: in Integer; Found: out Integer) is
      First: access Filter;
      E: Element:= Unit;
   begin
      First := new Filter;
      for I in 1 .. To_Try loop
         E := Succ(E);
         First.Input(E);
      end loop;
      First.Stop;
      Finished.Wait;
      Found := Primes_Found;
   end Do_It;

end Eratosthenes;

package Integer_Stuff is
   function Is_Factor(N, P: Integer) return Boolean;
   procedure Put(N: in Integer);
end;

with Ada.Integer_Text_IO;
package body Integer_Stuff is

   function Is_Factor(N, P: Integer) return Boolean is
   begin
      return N mod P = 0;
   end Is_Factor;

   procedure Put(N: in Integer) is
   begin
      Ada.Integer_Text_IO.Put(N, 0);
   end Put;

end Integer_Stuff;

package Complex_Stuff is
   type Complex is
      record
         Rl, Im: Integer;
      end record;

   function Succ(C: Complex) return Complex;
   function Is_Factor(N, P: Complex) return Boolean;
   procedure Put(C: Complex);
end;

with Ada.Text_IO;  use Ada.Text_IO;
with Ada.Integer_Text_IO;  use Ada.Integer_Text_IO;
package body Complex_Stuff is

   function Succ(C: Complex) return Complex is
      -- The successor is the complex number whose real part
      -- is one more and imaginary part one less than the given 
      -- number unless the given number is on the real axis in which
      -- case it is the number with real part 1 on the next parallel
      -- diagonal line further out from the origin.

      Rl: Integer := C.Rl;
      Im: Integer := C.Im;
   begin
      if Im = 0 then
         return (1, Rl);
      end if;
      return (Rl+1, Im-1);
   end Succ;

   function Is_Factor(N, P: Complex) return Boolean is
      -- The approach is to multiply both N and P by the 
      -- conjugate of P and then to check for divisibility
      -- of the real and imaginary parts separately.

      P_Mod: Integer := P.Rl*P.Rl + P.Im*P.Im;
   begin
      return (N.Rl*P.Rl + N.Im*P.Im) mod P_Mod = 0 and
             (N.Im*P.Rl - N.Rl*P.Im) mod P_Mod = 0;
   end Is_Factor;

   procedure Put(C: in Complex) is
      -- Outputs a complex number in the form a + bi but omits
      -- the imaginary part if it is zero.
   begin
      Put(C.Rl, 0);
      if C.Im /= 0 then
         Put(" + ");  Put(C.Im, 0);  Put('i');
      end if;
   end Put;

end Complex_Stuff;

package Poly_Stuff is
   Max: constant Integer := 100;
   type Coeff is mod 2;
   subtype Index is Integer range 0 .. Max;
   type Coeff_Vector is array (Integer range <>) of Coeff;
   type Polynomial(N: Index := 0) is
      record
         A: Coeff_Vector(0 .. N);
      end record;

   function Succ(P: Polynomial) return Polynomial;
   function Is_Factor(N, P: Polynomial) return Boolean;
   procedure Put(P: in Polynomial);
end;

with Ada.Text_IO;  use Ada.Text_IO;
with Ada.Integer_Text_IO;  use Ada.Integer_Text_IO;
package body Poly_Stuff is

   function Succ(P: Polynomial) return Polynomial is
      -- Consider the binary coefficients of the polynomials as
      -- representing an integer.  The successor is then that
      -- polynomial represented by the next integer.  The
      -- successor polynomial is therefore of the next order when
      -- the next integer is an exact power of 2. When this happens
      -- the return statement following the loop is taken.

      N: Integer := P.N;
      A: Coeff_Vector := P.A;
   begin
      for I in 0 .. N-1 loop
         A(I) := A(I) + 1;
         if A(I) = 1 then
            return (N, A);
         end if;
      end loop;
      return (N+1, (0 .. N => 0) & (N+1 => 1));
   end Succ;

   function Is_Factor(N, P: Polynomial) return Boolean is
      -- This just uses the normal long division algorithm
      -- simplified because each multiplier can only be 0 or 1
      -- since the coefficients can only be 0 or 1.

      Q: Polynomial := N;
   begin
      while Q.N >= P.N loop
         if Q.A(Q.N) = 1 then
            for I in 1 .. P.N loop
               Q.A(Q.N-I) := Q.A(Q.N-I) - P.A(P.N-I);
            end loop;
         end if;
         Q := (Q.N-1, Q.A(0 .. Q.N-1));
      end loop;
      return Q.A = (0 .. Q.N => 0);
   end Is_Factor;

   procedure Put(P: in Polynomial) is
      -- Outputs a polynomial giving just those terms with a 
      -- non-zero coefficient.  The coefficient of the highest
      -- term is always 1.

   begin
      Put('x');  Put(P.N, 0);
      for I in reverse 0 .. P.N-1 loop
         if P.A(I) /= 0 then
            Put(" + ");  Put('x');  Put(I, 0);
         end if;
      end loop;
   end Put;

end Poly_Stuff;

with Eratosthenes;
with Integer_Stuff;  use Integer_Stuff;
package Integer_Sieve is
   new Eratosthenes(Integer, 1, Integer'Succ, Is_Factor, Put);

with Eratosthenes;
with Complex_Stuff;  use Complex_Stuff;
package Complex_Sieve is
   new Eratosthenes(Complex, (1, 0), Succ, Is_Factor, Put);

with Eratosthenes;
with Poly_Stuff;  use Poly_Stuff;
package Poly_Sieve is
   new Eratosthenes(Polynomial, (0, (0 => 1)), Succ, Is_Factor, Put);

with Integer_Sieve;
with Complex_Sieve;
with Poly_Sieve;
with Ada.Text_IO;  use Ada.Text_IO;
with Ada.Integer_Text_IO;  use Ada.Integer_Text_IO;
procedure Super_Sieve is
   Char: Character;
   To_Try: Integer;
   Found: Integer := 0;
begin
   Put("Welcome to the Super Sieve");
   New_Line(2);
   Put("Type I for Integers, C for Complex, " &
       "P for Polynomials  ");
   Get(Char);
   Put("How many do you want to try?  ");
   Get(To_Try);
   New_Line;
   case Char is
      when 'I' | 'i' =>
         Integer_Sieve.Do_It(To_Try, Found);
      when 'C' | 'c' =>
         Complex_Sieve.Do_It(To_Try, Found);
      when 'P' | 'p' =>
         Poly_Sieve.Do_It(To_Try, Found);
      when others =>
         Put("No such sieve");
   end case;
   New_Line(2);
   Put("Number of primes found was ");
   Put(Found, 0);
   New_Line(2);
   Put_Line("Finished");
   Skip_Line(2);
end Super_Sieve;