program Hensel;                   {Lift solutions of a congruence (mod p) to
                                   higher powers of p.}
{$N+,E+}
uses CRT, nothy;
{$I GetInput.i }

const
MaxTerms = 20;
MaxZeros = 2000;
MaxMod = 2000;

var
a,                                {a residue class}
c,                                {a residue class (mod p)}
p,                                {the prime modulus}
m,                                {a modulus}
r,                                {a root mod p^j}
r1                                {a root mod p^k lying below r}
    : comp;
x   : extended;
h,                                {index of residues (mod p), used in
                                   singular case}
is,                               {index of singular roots}
ins,                              {index of non-singular roots}
j,                                {the power of p in top row of display
                                   in the case of a singular root;
                                   the power of p at the bottom of the
                                   display for non-singular roots}
k,                                {an index of powers of p}
MaxJ,                             {largest j s.t. p^j <= 1E18}
l,                                {length of a string}
t,                                {number of terms in the polynomial}
x0, y0,                           {cursor coordinates}
ind,                              {amount to indent by}
NoNSR,                            {number of non-singular roots mod p}
NoSR                              {number of singular roots mod p}
    : longint;
Coefficient,                      {coefficient of monomial of polynomial}
Exponent                          {exponent of monomial of polynomial}
    : array[1..MaxTerms] of comp;
NSR,                              {the non-singular roots}
SR,                               {the singular roots}
u                                 {factor used in lifting}
    : array[1..MaxZeros] of comp;
pj                                {pj[j] = p^j}
    : array[0..64] of comp;
ModSet,                           {prime modulus successfully defined}
FuncKey,                          {Was a function key pressed?}
InputOk,                          {Is the key pressed admissible?}
Singular                          {Display singular roots?}
    : Boolean;
OriginalMode
    : word;
Ch                              {input for keyboard}
    : char;
St1, St2, St3,                  {prompts}
Stfc,                           {coefficients of f(x) expressed as a string}
Stfe                            {exponents of f(x) expressed as a string}
    : string[255];

procedure DefPoly;

var
i
    : integer;
c
    : comp;

begin
  for i := 21 to 25 do
    begin
      GoToXY(1, i);
      ClrEoL
    end;
  x0 := 7;
  y0 := 22;
  GoToXY(7, 23);
  Write('Enter monomial terms ax^k ');
  Write(' with exponents k strictly decreasing.');
  t := 0;
  c := 1;
  while ((t < MaxTerms) and (c <> 0)) or (t = 0) do
    begin
      str((t+1):1, St1);
      St2 := 'Enter coefficient a of term #' + St1 + ', a = ';
      c := GetInput(7, 24, St2,
           '                                     (a = 0 to terminate)',
           -1E18, 1E18);
      if c <> 0
      then
        begin
          t := t + 1;
          coefficient[t] := c;
          St2 := 'Enter exponent k of term #' + St1 + ', k = ';
          if t = 1
          then c := GetInput(7, 24, St2, '', 0, 1E18)
          else
            begin
              str((Exponent[t-1]-1):1:0, St3);
              St3 := '                              (0 ó k ó ' + St3 + ')';
              c := GetInput(7, 24, St2, St3, 0, Exponent[t-1] - 1)
            end;
          exponent[t] := c;
          str(Abs(coefficient[t]):1:0, St1);
          l := length(St1);
          str(exponent[t]:1:0, St1);
          l := l + length(St1);
          if l + x0 > 75
          then
            begin
              GoToXY(1, 23);
              ClrEoL;
              GoToXY(1, 24);
              ClrEoL;
              GoToXY(1, 25);
              ClrEoL;
              WriteLn;
              WriteLn;
              WriteLn;
              x0 := 9;
              y0 := 22
            end;
          GoToXY(x0, y0);
          if (coefficient[t] > 0) and (t > 1)
          then Write(' + ');
          if coefficient[t] < 0
          then
            if t > 1
            then Write(' - ')
            else Write(' -');
          if (Abs(coefficient[t]) > 1) or (exponent[t] = 0)
          then Write(Abs(coefficient[t]):1:0);
          if exponent[t] > 0
          then Write('x');
          if exponent[t] > 1
          then
            begin
              GoToXY(WhereX, WhereY - 1);
              Write(exponent[t]:1:0);
              GoToXY(WhereX, WhereY + 1)
            end;
          x0 := WhereX;
          y0 := WhereY;
        end;
    end;
end;                               {of procedure DefPoly}

function poly(x, m : comp) : comp;

var
a,                                {a residue class (mod m)}
d,                                {difference between successive exponents}
y                                 {a power of x}
             : comp;
i                                 {an index}
             : integer;

begin
  a := coefficient[1];
  for i := 1 to t - 1 do
    begin
      d := exponent[i] - exponent[i+1];
      y := power(x, d, m);
      a := mult(a, y, m);
      a := condition(a + coefficient[i+1], m);
    end;
  y := power(x, exponent[t], m);
  poly := mult(a, y, m)
end;                                          {of function poly}

function deriv(x, m : comp) : comp;

var
a, b,                             {residue classes (mod m)}
d,                                {difference between successive exponents}
y                                 {a power of x}
             : comp;
i                                 {an index}
             : integer;

begin
  a := coefficient[1]*exponent[1];
  for i := 1 to t - 1 do
    begin
      d := exponent[i] - exponent[i+1];
      if exponent[i+1] = 0
      then d := d - 1;
      y := power(x, d, m);
      a := mult(a, y, m);
      b := mult(coefficient[i+1], exponent[i+1], m);
      a := condition(a + b, m);
    end;
  if exponent[t] > 1
  then
    begin
      y := power(x, exponent[t] - 1, m);
      deriv := mult(a, y, m)
    end
  else
    deriv := a;
end;                                          {of function deriv}

procedure Title;

var
i, j,                             {indices}
l                                 {length of string}
     : integer;

begin
  Stfc := 'ROOTS OF ';
  Stfe := '         ';
  for i := 1 to t do
    begin
      if coefficient[i] > 0
      then
        begin
          if i > 1
          then
            begin
              Stfc := Stfc + ' + ';
              Stfe := Stfe + '   '
            end
        end
      else
        begin
          if i > 1
          then
            begin
              Stfc := Stfc + ' - ';
              Stfe := Stfe + '   '
            end
          else
            begin
              Stfc := Stfc + ' -';
              Stfe := Stfe + '  '
            end
        end;
      if (Abs(coefficient[i]) > 1) or (exponent[i] = 0)
      then
        begin
          str(Abs(coefficient[i]):1:0, St1);
          Stfc := Stfc + St1;
          for j := 1 to length(St1) do Stfe := Stfe + ' '
        end;
      if exponent[i] > 0
      then
        begin
          Stfc := Stfc + 'x';
          Stfe := Stfe + ' '
        end;
      if exponent[i] > 1
      then
        begin
          str(exponent[i]:1:0, St1);
          Stfe := Stfe + St1;
          for j := 1 to length(St1) do Stfc := Stfc + ' '
        end;
    end;
  Stfc := Stfc + ' ð 0 (mod ';
  Stfe := Stfe + '          ';
  str(p:1:0, St1);
  Stfc := Stfc + St1 + ' )';
  for j := 1 to length(St1) do Stfe := Stfe + ' ';
  Stfe := Stfe + 'j';
  l := length(Stfc);
  if l > 78
  then
    begin
      Stfc := '                           ROOTS OF f(x) ð 0 (mod p )';
      Stfe := '                                                   j'
    end
  else
    begin
      for j := 1 to (78 - l) div 2 do
        begin
          Stfc := ' ' + Stfc;
          Stfe := ' ' + Stfe
        end
    end
end;                              {of Title}

procedure Roots;                  {find roots (mod p), and list them
                                   according to whether they are singular
                                   or non-singular}
var
 i                                {index of residue classes (mod p)}
    : longint;
begin
  NoSR := 0;
  NoNSR := 0;
  GoToXY(1, 25);
  ClrEoL;
  Write('                    Locating roots (mod ', p:1:0, ') . . . ');
  for i := 0 to round(p) - 1 do
    if poly(i, p) = 0
    then
      begin
        a := deriv(i, p);
        if a = 0
        then
          begin
            NoSR := NoSR + 1;
            SR[NoSR] := i
          end
        else
          begin
            NoNSR := NoNSR + 1;
            NSR[NoNSR] := i;
            lincon(-a, 1, p, a, m);
            u[NoNSR] := a
          end
      end;
  is := 1;
  ins := 1;
  j := 1
end;                              {of procedure Roots}

procedure Display;

begin
  TextColor(15);
  TextBackground(0);
  ClrScr;
  Title;
  WriteLn(Stfe);
  WriteLn(Stfc);
  GoToXY(19, 3);
  WriteLn('j             x                   c  ');
  if Singular
  then
    begin
      k := 1;
      if j - k > 9
      then k := j - 9;
      GoToXY(1, 23);
      while k <= j do
        begin
          r1 := condition(r, pj[k]);
          str(r1:1:0, St1);
          l := length(St1);
          for ind := 1 to (21 - l) div 2 do St1 := ' ' + St1;
          Write(k:20, '  ');
          TextColor(14);
          TextBackground(1);
          Write(St1);
          while WhereX < 47 do Write(' ');
          c := r1/pj[k-1];
          if c*pj[k-1] > r1
          then c := c-1;
          Write(c:8:0, '  ');
          if k < j
          then
            begin
              GoToXY(23, WhereY-1);
              if p = 2
              then Write('        \   /                     ');
              if p = 3
              then Write('        \ ³ /                     ');
              if p > 3
              then Write('        \\³//                     ');
              GoToXY(1, WhereY-1)
            end;
          if (k = j) and (j < MaxJ)
          then
            begin
              GoToXY(23, WhereY-1);
              if (poly(r, pj[j+1]) = 0)
              then
                begin
                  if p = 2
                  then Write('        \   /                     ');
                  if p = 3
                  then Write('        \ ³ /                     ');
                  if p > 3
                  then Write('        \\³//                     ')
                end
              else Write('                                  ')
            end;
          TextColor(15);
          TextBackground(0);
          k := k + 1
        end
    end
  else
    begin
      k := j;
      GoToXY(1, 23);
      while (k <= MaxJ) and (k <= j + 9) do
      begin
        str(r:1:0, St1);
        l := length(St1);
        for ind := 1 to (21 - l) div 2 do St1 := ' ' + St1;
        Write(k:20, '  ');
        TextColor(14);
        TextBackground(1);
        Write(St1);
        While WhereX < 47 do Write(' ');
        c := r/pj[k-1];
        if c*pj[k-1] > r
        then c := c - 1;
        Write(c:8:0, '  ');
        GoToXY(23, WhereY-1);
        Write('          ³                       ');
        GoToXY(1, WhereY-1);
        TextColor(15);
        TextBackground(0);
        if k < MaxJ
        then
          begin
            c := mult(u[ins], poly(r, pj[k+1])/pj[k], p);
            r := r + c*pj[k]
          end;
        k := k + 1
      end
    end;
  GoToXY(1, 24);
  Write('       f(x) ð 0 (mod ', p:1:0, ') has ');
  if NoNSR > 0
  then Write(NoNSR:1)
  else Write('no');
  Write(' non-singular root');
  if NoNSR > 1
  then Write('s');
  Write(' and ');
  if NoSR > 0
  then Write(NoSR:1)
  else Write('no');
  Write(' singular root');
  if NoSR > 1
  then Write('s');
  TextBackground(7);
  GoToXY(1, 25);
  ClrEoL;
  if (Singular and (j < MaxJ) and (poly(r, pj[j+1]) = 0))
      or ((not Singular) and (j + 10 < MaxJ))
  then TextColor(4)
  else TextColor(15);
  Write('   ', #24);              {UpArrow}
  if j > 1
  then TextColor(4)
  else TextColor(15);
  Write('      ', #25);             {DownArrow}
  if ((ins > 1) and (not Singular)) or ((is > 1) and (j = 1) and Singular)
      or ((h > 0) and Singular and (j > 1))
  then TextColor(4)
  else TextColor(15);
  Write('      ', #27);
  if (Singular and (j = 1) and (is < NoSR))
         or (Singular and (j > 1) and (h < p - 1))
         or ((not Singular) and (ins < NoNSR))
  then TextColor(4)
  else TextColor(15);
  Write('      ', #26);
  if Singular
  then
    if NoNSR > 0
    then
      begin
        TextColor(4);
        Write('      N');
        TextColor(0);
        Write('on-singular')
      end
    else
      begin
        TextColor(15);
        Write('      non-singular')
      end
  else
    if NoSR > 0
    then
      begin
        TextColor(4);
        Write('          S');
        TextColor(0);
        Write('ingular')
      end
    else
      begin
        TextColor(15);
        Write('          Singular')
      end;
  TextColor(4);
  Write('      D');
  TextColor(0);
  Write('efine poly');
  TextColor(4);
  Write('      p');
  TextColor(4);
  Write('      Esc');
  GoToXY(1, 25);                  {hide the cursor}
  TextColor(7);
  Write(' ');
  GoToXY(1, 25);
  TextColor(0)
end;                              {of procedure Display}

begin                             {main body}
  OriginalMode := LastMode;
  Ch := 'F';
  TextColor(14);
  TextBackground(1);
  ClrScr;
  GoToXY(23, 3);
  WriteLn('LIFT ROOTS VIA HENSEL''S LEMMA');
  GoToXY(1, 6);
  Write('      Given a polynomial f(x),');
  WriteLn(' and a prime number p, will locate all');
  Write('      roots of the congruence ');
  WriteLn('f(x) ð 0 (mod p).  Singular roots');
  Write('      (those for which f''(x) ð 0 (mod p))');
  WriteLn(' are distinguished from ');
  Writeln('      non-singular roots.  ');
  Write('           If c(1) is a non-singular root (mod p), then there');
  WriteLn(' is a ');
  Write('      unique c(2), 0 ó c(2) < p, such that c(1) + c(2)p ');
  WriteLn('is a root (mod pý).');
  WriteLn('      Roots are lifted in this way to yield roots (mod p^j);');
  Write('      they are expressed both as residue ');
  WriteLn('classes (mod p^j) and by the');
  Write('      associated sequence of coefficients c(j).');
  GoToXY(1, 15);
  DefPoly;
  repeat                          {next define an initial prime modulus}
    str(MaxMod:1, St1);
    St1 := '                           (0 < p ó ' + St1 + ')';
    p := round(GetInput(7, 23, '      Enter prime modulus p = ',
                         St1, 2, MaxMod));
    ModSet := True;
    if (p > 2) and (not spsp(2, p)) then ModSet := False;
    if (p > 3) and (not spsp(3, p)) then ModSet := False;
                                  {Thus p is proved to be prime,
                                   provided that p < 1373653}
    if ModSet
    then
      begin
        Roots;
        if (NoSR = 0) and (NoNSR = 0)
        then
          begin
            ModSet := False;
            GoToXY(6, 25);
            ClrEoL;
            TextColor(12);
            Write('       The congruence f(x) ð 0 (mod ');
            Write(p:1:0, ') has no root.');
            TextColor(14)
          end
      end
    else
      begin
        GoToXY(7, 25);
        ClrEoL;
        TextColor(12);
        Write('      The number ', p:1:0, ' is composite.');
        TextColor(14)
      end
  until ModSet;
  MaxJ := -1;                     {determine the highest admissible    }
  x := 1;                         {power of p                          }
  while x < MaxAllow do
    begin
      MaxJ := MaxJ + 1;
      pj[MaxJ] := x;
      x := x*p
    end;
  if NoNSR > 0
  then
    begin
      Singular := False;
      r := NSR[1]
    end
  else
    begin
      Singular := True;
      r := SR[1]
    end;
  repeat
    if (not FuncKey) and (UpCase(Ch) = 'D')     {Define Poly}
    then
      begin
        DefPoly;
        Roots;
        if (NoNSR = 0) and (NoSR = 0)
        then
          begin
            GoToXY(1, 25);
            ClrEoL;
            TextColor(4);
            Write('             The congruence f(x) ð 0 (mod ');
            Write(p:1:0, ') has no root.');
            TextColor(0);
            repeat                {choose p so that f(x)ð0(p) has a root}
              str(MaxMod:1, St1);
              St1 := '                            (0 < p ó ' + St1 + ')';
              p := round(GetInput(7, 23, '       Enter prime modulus p = ',
                                   St1, 2, MaxMod));
              ModSet := True;
              if (p > 2) and (not spsp(2, p)) then ModSet := False;
              if (p > 3) and (not spsp(3, p)) then ModSet := False;
              if ModSet
              then
                begin
                  Roots;
                  if (NoSR = 0) and (NoNSR = 0)
                  then
                    begin
                      ModSet := False;
                      GoToXY(7, 25);
                      ClrEoL;
                      TextColor(4);
                      Write('      The congruence f(x) ð 0 (mod ');
                      Write(p:1:0, ') has no root.');
                      TextColor(0)
                    end
                end
              else
                begin
                  GoToXY(1, 25);
                  ClrEoL;
                  TextColor(4);
                  Write('             The number ');
                  Write(p:1:0, ' is composite.');
                  TextColor(0)
                end
            until ModSet;
          end;
        MaxJ := -1;           {determine the highest admissible    }
        x := 1;               {power of p                          }
        while x < MaxAllow do
          begin
            MaxJ := MaxJ + 1;
            pj[MaxJ] := x;
            x := x*p
          end;
        if NoNSR > 0
        then
          begin
            Singular := False;
            r := NSR[1]
          end
        else
          begin
            Singular := True;
            r := SR[1]
          end
      end;                        {of define poly}
    if (FuncKey) and (Ch = #72)                          {lift}
    then
      begin
        if Singular
        then
          begin
            j := j + 1;
            h := 0
          end
        else
          begin
            j := j + 10;
            if j + 9 > MaxJ
            then
              begin
                j := MaxJ - 9;
                r := condition(r, pj[j])
              end
          end
      end;                        {of lifting}
    if (FuncKey) and (Ch = #80)   {drop}
    then
      begin
        if Singular
        then
          begin
            j := j - 1;
            r := condition(r, pj[j]);
            if j > 1
            then
              begin
                h := round(r/pj[j-1]);
                if h*pj[j-1] > r
                then h := h - 1
              end
          end
        else
          begin
            j := j - 10;
            if j < 1
            then j := 1;
            r := condition(r, pj[j])
          end
      end;                        {of dropping}
    if (FuncKey) and (Ch = #75)   {left}
    then
      begin
        if not Singular
        then
          begin
            ins := ins - 1;
            j := 1;
            r := NSR[ins]
          end;
        if Singular and (j > 1)
        then
          begin
            h := h - 1;
            r := r - pj[j-1]
          end;
        if Singular and (j = 1)
        then
          begin
            is := is - 1;
            r := SR[is]
          end
      end;                        {of left}
    if (FuncKey) and (Ch = #77)   {move right}
    then
      begin
        if not Singular
        then
          begin
            ins := ins + 1;
            j := 1;
            r := NSR[ins]
          end;
        if Singular and (j = 1)
        then
          begin
            is := is + 1;
            r := SR[is]
          end;
        if Singular and (j > 1)
        then
          begin
            h := h + 1;
            r := r + pj[j-1]
          end
      end;                        {of right}
    if (not FuncKey) and (UpCase(Ch) = 'P')              {Set modulus p}
    then
      begin
        ClrEoL;
        GoToXY(1, 23);
        ClrEoL;
        str(MaxMod:1, St1);
        St1 := '                            (0 < p ó ' + St1 + ')';
        repeat
          p := round(GetInput(1, 23,
                              '             Enter prime modulus p = ',
                              St1, 2, MaxMod));
          ModSet := True;
          if (p > 2) and (not spsp(2, p)) then ModSet := False;
          if (p > 3) and (not spsp(3, p)) then ModSet := False;
          if ModSet
          then
            begin
              Roots;
              if (NoNSR = 0) and (NoSR = 0)
              then
                begin
                  ModSet := False;
                  GoToXY(1, 25);
                  ClrEoL;
                  TextColor(4);
                  Write('             The congruence f(x) ð 0 (mod ');
                  Write(p:1:0, ') has no root.');
                  TextColor(0)
                end
            end
          else
            begin
              GoToXY(1, 25);
              ClrEoL;
              TextColor(4);
              Write('             The number ', p:1:0, ' is composite.');
              TextColor(0)
            end
        until ModSet;
        MaxJ := -1;
        x := 1;
        while x < MaxAllow do
          begin
            MaxJ := MaxJ + 1;
            pj[MaxJ] := x;
            x := x*p
          end;
        if NoNSR > 0
        then
          begin
            Singular := False;
            r := NSR[1]
          end
        else
          begin
            Singular := True;
            r := SR[1]
          end
      end;                                  {modulus p set}
    if (not FuncKey) and (UpCase(Ch) = 'S') {switch to singular roots}
    then
      begin
        Singular := True;
        j := 1;
        r := SR[is]
      end;                                  {of switch to singular}
    if (not FuncKey) and (UpCase(Ch) = 'N') {switch to non-singular roots}
    then
      begin
        Singular := False;
        j := 1;
        r := NSR[ins]
      end;                                  {of switch to non-singular}
    Display;
    repeat                        {evaluate input from keyboard}
      Ch := ReadKey;
      if Ch <> #0
      then FuncKey := false
      else
        begin
          FuncKey := true;
          Ch := ReadKey
        end;
      InputOk := false;
      if (not FuncKey) and (Ch = #27)
      then InputOk := True;                      {Esc}
      if (not FuncKey) and (UpCase(Ch) IN ['D', 'P'])
      then InputOk := true;                      {Define poly or choose p}
      if (not Singular) and (j + 10 < MaxJ) and (FuncKey) and (Ch = #72)
      then InputOk := True;                      {UpArrow}
      if Singular and (j < MaxJ) and (FuncKey)
          and (Ch = #72) and (poly(r, pj[j+1])=0)
      then InputOk := True;                      {UpArrow}
      if (j > 1) and (FuncKey) and (Ch = #80)
      then InputOk := True;                      {DownArrow}
      if (not Singular) and (ins > 1) and (FuncKey) and (Ch = #75)
      then InputOk := True;                      {LeftArrow}
      if Singular and (h > 0) and (j > 1) and (FuncKey) and (Ch = #75)
      then InputOk := True;                      {LeftArrow}
      if Singular and (is > 1) and (j = 1) and (FuncKey) and (Ch = #75)
      then InputOk := True;                      {LeftArrow}
      if (not Singular) and (ins < NoNSR) and (FuncKey) and (Ch = #77)
      then InputOk := True;                      {RightArrow}
      if Singular and (is < NoSR) and (j = 1) and (FuncKey) and (Ch = #77)
      then InputOk := True;                      {RightArrow}
      if Singular and (h < p - 1) and (j > 1) and (FuncKey) and (Ch = #77)
      then InputOk := true;                      {RightArrow}
      if Singular and (NoNSR > 0)
            and (not FuncKey) and (UpCase(Ch) = 'N')
      then InputOk := True;                      {switch to non-singular}
      if (not Singular) and (NoSR > 0)
            and (not FuncKey) and (UpCase(Ch) = 'S')
      then InputOk := True                       {switch to singular}
    until InputOk;
  until (not funcKey) and (Ch = #27);            {until Esc pressed}
  TextMode(OriginalMode)
end.