-- prob7a.adb   20,000 by 20,000 matrix
-- a(i,j)=i th prime if i=j
-- a(i,j)=1  if |i-j| is 1 or a power of 2
-- a(i,j) otherwise = 0
-- invert and get value of a(1,1) to 100 digits or more

-- note 1/2 the elements are zero and need no computation
-- note the matrix is symmetric and thus the inverse also
-- done with computing A(1,1) when k=N


with Ada.Text_IO ; use Ada.Text_IO ;

procedure Prob7a is
  subtype Real is Long_Float;
  package Real_IO is new float_IO(Real); use Real_Io;
  package Int_IO is new INTEGER_IO(Integer); use Int_Io;
  package Bool_Io is new Enumeration_Io(Boolean); use Bool_Io;

    A  : array(1..1000,1..1000) of real;
    P  : array(1..20000) of integer;
    P2 : array(0..20000) of boolean;
    T  : Integer;
    K  : Integer;
    I  : Integer;

    Pivot : Real;
begin

    Put_Line("prob7a.adb  \n");
    -- prime to compute primes in an array P(1)..P(N)
    P(1) := 2;
    P(2) := 3;
    P(3) := 5;
    P(4) := 7;
    K := 4;  -- number of primes found
    I := 11; -- next to check for prime
    loop
      for J in 1..K loop
        T := I/P(J);
        if T<P(J) then
          K := K+1;
          P(K) := I;
          exit;
        end if;
        if T*P(J)=I then
          exit; -- has a factor, not prime
        end if;
      end loop;
      if K=20000 then exit; end if;
      I := I+2;
    end loop;
    Put("K="); Put(K); Put("  P(K)="); Put(P(K)); New_line;

    -- P2  to compute powers of 2 P2(1)..P2(N)
    for J in 0..20000 loop P2(J) := False; end loop;
    I := 1;
    for J in 1..14 loop
      P2(I) := true;
      I := I+I;
      if I>20000 then exit; end if;
    end loop;
    Put("P2(1)=");  Put(P2(1)); New_line;
    Put("P2(2)=");  Put(P2(2)); New_line;
    Put("P2(3)=");  Put(P2(3)); New_line;
    Put("P2(4)=");  Put(P2(4)); New_line;


  for N in 1..20 loop
    -- set up A matrix
    for I in 1..N loop
      for J in 1 .. N loop
        A(i,j) := 0.0;
        if I=J then A(I,J) := Real(P(I)); end if;
        if P2(abs(I-J)) then A(I,J) := 1.0; end if;
      end loop;
    end loop;

    --for I in 1..N loop
    --  for J in 1..N loop
    --    Put("A("); Put(I); Put(","); Put(J); Put(")="); Put(A(I,J)); New_Line;
    --  end loop;
    --end loop;

    -- invert A matrix
    for K in 1..N loop
      Pivot := A(K,K);
      if Pivot = 0.0 then exit; end if; -- should not happen
      --  reduce about pivot
      A(k,k) := 1.0 / pivot ;
      for J in 1..N loop
        if J /= K and A(K,J) /= 0.0 then
           A(K,J) := A(K,J)*A(K,K);
        end if;
      end loop;
      --  inner reduction loops
      for I in 1..N loop
        if K /= I and A(I,K) /= 0.0 then  -- half the elements in this matrix
          for J in 1..N loop
            if K /=J and A(K,J) /= 0.0 then
               A(I,J) := A(I,J)-A(I,K)*A(K,J);
               if K=N and I=1 and J=1 then goto done; end if;
            end if;
          end loop;
          A(I,K) := -A(I,K)*A(K,K);
        end if;
      end loop;
    end loop;
<<Done>>
    -- matrix inverted
    --for I in 1..N loop
    --  for J in 1..N loop
    --    Put("AI("); Put(I); Put(","); Put(J); Put(")="); Put(A(I,J)); New_Line;
    --  end loop;
    --end loop;

    Put("N="); Put(N); Put("  A(1,1)="); Put(A(1,1)); New_Line;
  end loop;
end Prob7a;