program IndDem;                   {Demonstrate Shanks'
                                   index-finding algorithm}
{$N+,E+}
uses CRT, nothy;
{$I GetInput.i }

type
Sequence = array[1..10000] of longint;

var
p,                                {the modulus; presumed to be prime}
g,                                {the base; presumed to be a primitive root}
h,                                {gh ð 1 (mod p)}
a,                                {number whose index is sought}
c,                                {g^s (mod p)}
d,                                {g^(i*s) (mod p)}
i                                 {exponent is  i*s + j}
  : comp;
z                                 {variable used for testing inputs}
  : extended;
B                                 {the sequence to be sorted}
  : Sequence;
Index                             {original index of member of sorted array}
  : array[1..10000] of integer;
j,                                {an index of terms in the sequence}
j0,                               {index returned by Search}
j1,                               {index of element hit, in sorted array
                                   if d < B[1] then j1 = -1;
                                   if d > B[s] then j1 = 1}
s,                                {the number of terms in the sequence}
i0,                               {lower endpoint of window, in search}
i1,                               {upper endpoint of window, in search}
l                                 {length of s, as a string}
  : integer;
t                                 {temporary variable used for swapping}
  : longint;
InputOk
  : Boolean;
code                              {error level when reading inputs}
  : word;
St                                {s, as a string}
  : string[20];

procedure Prompt;

var
x0, y0                    {coordinates of cursor location}
         : integer;
Ch       : Char;

begin
  x0 := WhereX;
  y0 := WhereY;
  if y0 = 25
  then
    begin
      WriteLn;
      y0 := 24
    end;
  GoToXY(1, 25);
  ClrEoL;
  Write('                       Press any key to continue . . . ');
  Ch := ReadKey;
  if Ch = #0
  then Ch := ReadKey;
  GoToXY(1, 25);
  ClrEoL;
  GoToXY(x0, y0)
end;                               {of procedure Prompt}

procedure HeapSort;

var
  i
  : integer;

procedure Drop(i, k : integer);

var
i0,
i1,
j,
t,
t1
   : longint;

begin
  i1 := i;
  t  := B[i];
  t1 := Index[i];
  repeat
    i0 := i1;
    j := 2*i1;
    if (j = k) and (B[k] > B[t])
    then
      begin
        B[i1] := B[k];
        Index[i1] := Index[k];
        i1 := k
      end;
    if j < k
    then
      begin
        if B[j+1] > B[j]
        then j := j + 1;
        if B[j] > t
        then
          begin
            B[i1] := B[j];
            Index[i1] := Index[j];
            i1 := j
          end
      end;
  until i0 = i1;
  B[i1] := t;
  Index[i1] := t1
end;                              {of Drop}

begin                             {main body of HeapSort}
  for i := s div 2 downto 1 do Drop(i, s);
  for i := s downto 2 do
    begin
      t := B[i];
      B[i] := B[1];
      B[1] := t;
      t := Index[i];
      Index[i] := Index[1];
      Index[1] := t;
      if i > 2
      then Drop(1, i-1)
    end
end;                              {of HeapSort}

function Search(d : comp) : integer;
                                  {Search for d in array B, return Index
                                   if found; else return -1}
var
im,                               {middle of window}
j                                 {the index sought}
  : integer;

begin
  i0 := 1;
  i1 := s;
  if d = B[1]
  then
    begin
      Search := Index[1];
      Exit
    end;
  if d = B[s]
  then
    begin
      Search := Index[s];
      Exit
    end;
  if d < B[1]
  then
    begin
      Search := -1;
      j1 := -1;
      Exit
    end;
  if d > B[s]
  then
    begin
      Search := -1;
      j1 := 1;
      Exit
    end;
  j := -1;
  repeat
    im := (i0+i1) div 2;
    if d = B[im]
    then
      begin
        j := Index[im];
        j1 := im
      end;
    if d < B[im]
    then i1 := im;
    if d > B[im]
    then i0 := im
  until (j >= 0) or (i1 - i0 = 1);
  Search := j
end;                              {of function Search}


begin                             {main}
  InputOk := False;
  if ParamCount = 3
  then
    begin
      InputOk := True;
      Val(ParamStr(1), z, code);
      if (Abs(z) <= 999999999) and (frac(z) = 0) and (code = 0)
      then g := z
      else InputOk := False;
      Val(ParamStr(2), z, code);
      if (Abs(z) <= 999999999) and (frac(z) = 0) and (code = 0)
      then a := z
      else InputOk := False;
      Val(ParamStr(3), z, code);
      if (z > 1) and (z <= 999999999) and (frac(z) = 0) and (code = 0)
      then p := z
      else InputOk := True;
    end;
  if not InputOk
  then
    begin
      ClrScr;
      WriteLn;
      Write('Will demonstrate method of solving g');
      GoToXY(WhereX, WhereY-1);
      Write('x');
      GoToXY(WhereX, WhereY+1);
      WriteLn(' ð a (mod p), so that x = ind a.');
      g := GetInput(WhereX, WhereY, 'Enter g = ', '   (|g| ó 999999999)',
                    -999999999, 999999999);
      a := GetInput(WhereX, WhereY, 'Enter a = ', '   (|a| ó 999999999)',
                    -999999999, 999999999);
      p := GetInput(WhereX, WhereY, 'Enter p = ', ' (2 ó p ó 999999999)',
                    2, 999999999)
    end
  else
    begin
      ClrScr;
      WriteLn;
      Write('Will demonstrate method of solving g');
      GoToXY(WhereX, WhereY-1);
      Write('x');
      GoToXY(WhereX, WhereY+1);
      WriteLn(' ð a (mod p), so that x = ind a.');
      WriteLn('     g = ', g:1:0);
      WriteLn('     a = ', a:1:0);
      WriteLn('     p = ', p:1:0)
    end;
  WriteLn;
  Prompt;
  z := gcd(g, p);
  if z > 1
  then
    begin
      WriteLn('ERROR!  This program works only when (g, p) = 1, but');
      WriteLn('(', g:1:0, ', ', p:1:0, ') = ', z:1:0, ' > 1.');
      Halt
    end;
  WriteLn('Having checked that (g, p) = 1, we use LinCon to find h so that');
  WriteLn('g*h ð 1 (mod p)');
  lincon(g, 1, p, h, p);      {gh ð 1 (mod p)}
  Writeln('      h = ', h:1:0);
  Prompt;
  if p > 100000000
  then s := 10000
  else s := round(sqrt(p));
  if s < 3 then s := 3;
  WriteLn('Next we determine the value of s.');
  WriteLn('      s = ', s:1);
  str(s:1, St);
  l := length(St);
  Prompt;
  Write('Next we construct an array of numbers a*h^j (mod p),');
  WriteLn(' for j = 0 .. s-1.');
  WriteLn('As we proceed, we check for an occurance of a*h^j ð 1 (mod p),');
  WriteLn('in which case j is the desired exponent.  We also check for ');
  WriteLn('an occurance of a*h^j ð a (mod p), in which case j is the order');
  WriteLn('of g (mod p), and a is not a power of g (mod p).');
  B[1] := round(condition(a, p));
  Write('B[1] = ', B[1]:1);
  Index[1] := 0;
  for j := 2 to s do
    begin
      Index[j] := j - 1;
      B[j] := round(mult(B[j-1], h, p));
      GoToXY(1, WhereY);
      ClrEoL;
      Write('B[', j:1, '] = ', B[j]:1);
      if B[j] = 1
      then
        begin
          WriteLn;
          WriteLn('While constructing array, discovered that');
          WriteLn;
          Write(g:1:0);
          GoToXY(WhereX, WhereY-1);
          Write((j-1):1);
          GoToXY(WhereX, WhereY+1);
          WriteLn(' ð ', a:1:0, ' (mod ', p:1:0, ')');
          Halt
        end;
      if B[j] = B[1]
      then
        begin
          WriteLn;
          WriteLn('While constructing array, discovered that');
          Write(a:1:0, ' is not a power of ', g:1:0);
          WriteLn(' modulo ', p:1:0, ', because ', g:1:0, ' has order ');
          Write((j-1):1, ' (mod ', p:1:0, '), and ', a:1:0, ' is not');
          WriteLn(' among these powers.');
          Halt
        end
    end;
  GoToXY(1, WhereY);
  ClrEoL;
  WriteLn('The array has been constructed.  The numbers in the array are: ');
  for j := 1 to s do
    begin
      Write(B[j]:1);
      if j < s
      then Write(', ');
      if WhereX + l > 78
      then WriteLn
    end;
  WriteLn;
  Prompt;
  Write('Next will use HeapSort to ');
  WriteLn('sort these numbers.');
  Prompt;
  Write('Sorting . . . ');
  HeapSort;
  WriteLn('done.  The numbers are now:');
  for j := 1 to s do
    begin
      Write(B[j]:1);
      if j < s
      then Write(', ');
      if WhereX + l > 78
      then WriteLn
    end;
  WriteLn;
  Prompt;
  i := 1;
  c := power(g, s, p);
  d := c;
  WriteLn('Next we use binary subdivisions to search for g^(i*s) (mod p)');
  WriteLn('in the array, for i = 1, 2, ..., until i*s > p.');
  Prompt;
  while i*s <= p do
    begin
      Write('Searching for g^(', i:1:0);
      WriteLn('*s) = ', d:1:0, ' in the array,');
      j0 := Search(d);
      if j0 >= 0
      then
        begin
          Write('and we found it! ', d:1:0, ' = B[', j1:1, '] ð a*h^');
          WriteLn(j0:1, ' (mod p).  Hence');
          Write('g^(', i:1:0, '*s) ð a*h^', j0:1);
          WriteLn(' (mod p),');
          WriteLn('so that g^(', i:1:0, '*s + ', j0:1, ') ð a (mod p).');
          WriteLn('That is, ');
          WriteLn;
          Write(g:1:0);
          GoToXY(WhereX, WhereY-1);
          Write((i*s + j0):1:0);
          GoToXY(WhereX, WhereY+1);
          WriteLn(' ð ', a:1:0, ' (mod ', p:1:0, ')');
          Halt
        end
      else
        begin
          if i1 = i0 + 1
          then
            begin
              Write('but we find that ', d:1:0, ' lies between consecutive');
              WriteLn(' members,');
              Write('B[', i0:1, '] = ', B[i0]:1, ' < ', d:1:0);
              WriteLn(' < ', B[i1]:1 , ' = B[', i1:1, ']');
              Prompt
            end;
          if (i1 > i0 + 1) and (j1 = 1)
          then
            begin
              WriteLn('but we find that ', d:1:0, ' > ', B[s]:1, ' = B[s]');
              Prompt
            end;
          if (i1 > i0 + 1) and (j1 = -1)
          then
            begin
              WriteLn('but we find that ', d:1:0, ' < ', B[1]:1, ' = B[1]');
              Prompt
            end
        end;
      d := mult(c, d, p);
      i := i + 1
    end;
  WriteLn('Since no number g^(i*s) (mod p) was found in the array,');
  WriteLn('for 0 ó i*s < p, we conclude that');
  Write(a:1:0, ' is not a power of ', g:1:0);
  WriteLn(' modulo ', p:1:0, '.')
end.