program WrgStTab;                 {Generate a TABle of STatistics concerning
                                   the number N(x) of integers n ó x that can
                                   be expressed as a sum of s k-th powers}
{$R+,N+,E+}
uses DOS, crt, printer;
{$I GetInput.i }
{$I Timer.i }

const
MaxN = 99999999999.0;             {10^11 - 1.  If k ò 3 then we store
                                   n^(1/k) k-th powers in an array.  Each
                                   entry occupies 8 bytes.  Thus we need
                                   8n^(1/3) < 64K}
MaxI0 = 1000;                     {maximum number of entries in array N}
MaxK = 10;
MaxS = 75;
PrinterLimit = 1000;              {This allows printing of the entire table}

var
k,                                {power}
s,                                {number of summands}
l,                                {length of string}
i,                                {index of array N}
i0,                               {N has been completed through i0}
i1,                               {value through which N is requested}
i2,                               {maximum allowable value of i1}
n0,                               {view of table begins at line n0}
n1                                {largest admissible value of n0}

  : integer;
h,                                {Display N(x) at multiples of h}
ll,                               {lower limit of table}
ul,                               {upper limit of table}
t                                 {temporary variable}
  : comp;
N                                 {values of N(x) at multiples of h}
  : array[0..MaxI0] of comp;
pwr                               {k-th powers, computed in advance}
  : array[0..4642] of comp;       {Powers are pre-computed only for k ò 3.
                                   Since ul < 10^11, we have x[i] ó 4641
                                   for all i.}
r                                 {r[i] = True if i+ll is represented}
  : array[0..15000] of Boolean;
OriginalMode                      {original video mode -- will be restored}
  : word;
Ch                                {character read from keyboard}
  : char;
St, St1,                          {prompts}
STi2                              {i2 expressed as a string}
  : string[80];
PrinterOn,                        {Is the printer responding?}
FuncKey,                          {Was the key pressed a function key?}
InputOk                           {Is the input ok?}
  : Boolean;

procedure MakeTable(ll, ul : comp);

var
m,                                {number represented}
y                                 {x[s]^k}
  : comp;
i,                                {an index}
i0,                               {marked index,
                                   used in lexicographic ordering}
MaxPwr                            {largest number whose k-th power is needed}
  : integer;
x                                 {variables---s of them}
  : array[0..MaxS] of longint;
t                                 {temporary variable}
  : extended;

begin
  MaxPwr := 1 + round(exp(ln(ul)/k));
  if k > 2
  then
    begin
      pwr[0] := 0;
      for i := 1 to MaxPwr do pwr[i] := exp(k*ln(i))
    end;
  for i := 0 to 15000 do r[i] := False;
  for i := 0 to s do x[i] := 0;
  if ll > 0
  then
    begin
      t := exp(ln(ll)/k);
      x[s] := round(t);
      if pwr[x[s]] < ll
      then x[s] := x[s]+1
    end;
  i0 := s;
  repeat
    m := 0;
    if k = 2                      {form the sum of s k-th powers}
    then for i := 1 to s do m := m + sqr(x[i])
    else for i := 1 to s do m := m + pwr[x[i]];
    if (m >= ll) and (m <= ul)
    then
      begin
        r[round(m-ll)] := True;
        x[s] := x[s] + 1
      end
    else
      begin
        i0 := s-1;
        while (x[i0] = x[i0+1]) and (i0 >= 1) do i0 := i0 - 1;
                   {find i0 so that x[i0] < x[i0+1] = ... = x[s]}
        if (i0 <= s - k) and (KeyPressed)    {allow termination}
        then
          begin
            TextMode(OriginalMode);
            Halt
          end;
        if i0 > 0
        then
          begin
            x[i0] := x[i0]+1;
            for i := i0 + 1 to s do x[i]  := x[i0]
          end;
        m := 0;
        if k = 2
        then
          begin
            for i := 1 to s - 1 do m := m + sqr(x[i]);
            y := sqr(x[s])
          end
        else
          begin
            for i := 1 to s - 1 do m := m + pwr[x[i]];
            y := pwr[x[s]]
          end;
        if m + y < ll
        then
          begin
            t := exp(ln(ll-m)/k);
            x[s] := round(t);
            if pwr[x[s]] < ll-m
            then x[s] := x[s]+1
          end
      end
  until i0 = 0
end;                              {of procedure MakeTable}

procedure Build(i1 : integer);    {build table of values of N through N(i1)}

var
  cl,                             {upper limit of calculation}
  dn,                             {change in N}
  t                               {temporary variable}
     : comp;

begin
  SetTimer;
  dn := 0;
  ul := h*i0;
  while i0 < i1 do
    begin
      ll := ul + 1;
      ul := ll + 15000;
      if ul > h*i1
      then ul := h*i1;
      MakeTable(ll, ul);
      for i := 0 to round(ul - ll) do
        begin
          if r[i] = True
          then dn := dn + 1;
          t := (i+ll)/h;
          if t*h = i + ll      {if h|(i+ll)}
          then
            begin
              N[i0+1] := N[i0] + dn;
              i0 := i0 + 1;
              dn := 0
            end
        end
    end;
  ReadTimer
end;                              {of procedure Build (array of values)}

procedure Display;

var
i,                                {index of rows in table}
j                                 {index of columns in table}
  : integer;

begin
  TextColor(15);
  TextBackground(0);
  ClrScr;
  if k = 2
  then
    begin
      Write('  NUMBER N(x) OF INTEGERS n ó x THAT CAN BE EXPRESSED AS A');
      WriteLn(' SUM OF ', s:1, ' SQUARES')
    end;
  if k = 3
  then
    begin
      Write('    NUMBER N(x) OF INTEGERS n ó x THAT CAN BE EXPRESSED AS A');
      WriteLn(' SUM OF ', s:1, ' CUBES')
    end;
  if k > 3
  then
    begin
      Write('    NUMBER N(x) OF INTEGERS n ó x THAT CAN BE EXPRESSED AS A');
      WriteLn(' SUM OF ', s:1, ' ',k:1, '-TH POWERS')
    end;
  WriteLn;
  Write('           x          N(x)       N(x)/x');
  WriteLn('        dN        dN/dx');
  i := n0;
  while (i <= i0) and (i- n0 <= 19) do
    begin
      t := h;
      t := i*t;
      Write(t:17:0, ' ');
      TextColor(14);
      TextBackground(1);
      Write(N[i]:11:0, '   ', N[i]/t:8:6);
      WriteLn('  ', (N[i]-N[i-1]):11:0,  '   ', (N[i]-N[i-1])/h:8:6, '  ');
      TextColor(15);
      TextBackground(0);
      i := i + 1
    end;
  if (not FuncKey) and (UpCase(Ch) IN ['E', 'F'])
  then
    begin
      Write('      It took ', TimerString, ' to ');
      if (not FuncKey) and (UpCase(Ch) = 'E')
      then Write('extend ')
      else Write('construct ');
      Write('this table.')
    end;
  GoToXY(1, 25);
  TextBackground(7);
  ClrEoL;
  if n0 > 1
  then TextColor(4)
  else TextColor(15);
  Write('     PgUp');
  if n0 + 19 < i0
  then TextColor(4)
  else TextColor(15);
  Write('        PgDn');
  TextColor(4);
  Write('            x');
  if i0 < 1000
  then
    begin
      TextColor(4);
      Write('         E');
      TextColor(0);
      Write('xtend')
    end
  else
    begin
      TextColor(15);
      Write('         Extend')
    end;
  if PrinterOn
  then
    begin
      TextColor(4);
      Write('        P');
      TextColor(0);
      Write('rint')
    end
  else
    begin
      TextColor(15);
      Write('        Print')
    end;
  TextColor(4);
  Write('        Esc');
  GoToXY(1, 25);                  {hide the cursor}
  TextColor(7);
  Write(' ');
  GoToXY(1, 25);
  TextColor(0)
end;                              {of procedure Display}

procedure PrinTab;                {print the table}

var
i, j                              {indices}
  : integer;

begin
  WriteLn(lst);
  if k = 2
  then
    begin
      Write(lst,'  NUMBER N(x) OF INTEGERS n <= x THAT CAN BE EXPRESSED ');
      WriteLn(lst, 'AS A SUM OF ', s:1, ' SQUARES')
    end;
  if k = 3
  then
    begin
      Write(lst, '    NUMBER N(x) OF INTEGERS n <= x THAT CAN BE EXPRESSED ');
      WriteLn(lst, 'AS A SUM OF ', s:1, ' CUBES')
    end;
  if k > 3
  then
    begin
      Write(lst, '    NUMBER N(x) OF INTEGERS n <= x THAT CAN BE EXPRESSED ');
      WriteLn(lst, 'AS A SUM OF ', s:1, ' ',k:1, '-TH POWERS')
    end;
  WriteLn(lst);
  Write(lst, '           x          N(x)       N(x)/x');
  WriteLn(lst, '        dN        dN/dx');
  for i := 1 to i0 do
    begin
      t := h;
      t := t*i;
      Write(lst, t:17:0, ' ');
      Write(lst, N[i]:11:0, '   ', N[i]/t:8:6);
      WriteLn(lst, '  ', (N[i]-N[i-1]):11:0, '   ', (N[i]-N[i-1])/h:8:6)
    end;
  WriteLn(lst)
end;                              {of procedure PrinTab}

begin                             {main body}
  OriginalMode := LastMode;
  TextBackground(1);
  TextColor(14);
  ClrScr;
  GoToXY(6, 10);
  Write('               STATISTICS ON WARING''S PROBLEM');
  GoToXY(6,12);
  WriteLn('Will construct a table of statistics concerning the number ');
  GoToXY(6, 13);
  WriteLn('N(x) of numbers n, 1 ó n ó x, that can be expressed as a sum');
  GoToXY(6, 14);
  WriteLn('of s k-th powers.  Values are reported when x = h, 2h, 3h, ...');
  WriteLn;
  str(MaxS:1, St);
  St := '                                 (1 ó s ó ' + St + ')';
  s := round(GetInput(8, WhereY, 'Choose the number of summands, s = ',
             St, 1, MaxS));
  GoToXY(1, WhereY-1);
  ClrEoL;
  WriteLn('                                    s = ', s:1);
  str(MaxK, St);
  St := '                                   (2 ó k ó ' + St + ')';
  k := round(GetInput(16, WhereY, 'Select the exponent, k = ',
                      St, 2, MaxK));
  GoToXY(1, WhereY - 1);
  ClrEoL;
  GoToXY(37, WhereY);
  WriteLn('k = ', k:1);
  h := GetInput(9, WhereY, 'Select frequency of report, h = ',
      '                               (1 ó h ó 1000000000)', 1, 1000000000);
  GoToXY(1, WhereY-1);
  ClrEoL;
  GoToXY(37, WhereY);
  WriteLn('h = ', h:1:0);
  i2 := 1000;
  if h > 10000000
  then i2 := round(1E11/h);
  if i2*h > 1E11
  then i2 := i2-1;
  str(h:1:0, St1);
  St1 := 'Compute up to x = ' + St1 + 'm, where m = ';
  while length(St1) < 40 do St1 := ' ' + St1;
  l := 40;
  St := '(1 ó m ó ';
  while length(St) < l do St := ' ' + St;
  str(i2:1, STi2);
  St := St + STi2 + ')';
  i1 := round(GetInput(1, WhereY, St1, St, 1, i2));
  GoToXY(1, WhereY-1);
  ClrEoL;
  GoToXY(37, WhereY);
  WriteLn('m = ', i1:1);
  WriteLn;
  Write('           Computing valuesÄÄÄpress any key to terminate . . . ');
  N[0] := 0;
  i0 := 0;
  Build(i1);
  Ch := 'F';
  FuncKey := False;
  n0 := 1;
  Repeat
    {$I-}
    Write(lst,#0);
    if IOResult = 0
    then PrinterOn := True
    else PrinterOn := False;
    {$I+}
    if (not FuncKey) and (UpCase(Ch) = 'X')
    then
      begin
        n1 := i0 - 19;
        str(n1:1, St);
        St := '                  (1 ó m ó ' + St + ')';
        n0 := round(GetInput(1, 24, 'Display x = mh, where  m = ',
                    St, 1, n1))
      end;
    if (not FuncKey) and (UpCase(Ch) = 'P')
    then
      begin
        ClrEoL;
        Write('              Printing table . . . ');
        Printab
      end;
    if (FuncKey) and (Ch = #73)     {PgUp}
    then
      begin
        n0 := n0 - 20;
        if n0 < 1 then n0 := 1
      end;
    if (FuncKey) and (Ch = #81)     {PgDn}
    then
      begin
        n0 := n0 + 20;
        if n0 + 19 > i0
        then n0 := i0 - 19
      end;
    if (not FuncKey) and (UpCase(Ch) = 'E')
    then                        {Extend the table}
      begin
        str(h:1:0, St1);
        St1 := '      Extend the table to x = ' + St1 + 'm, where m = ';
        l := length(St1);
        str(i0:1, St);
        St := '(' + St + ' < m ó ';
        while length(St) < l do St := ' ' + St;
        St := St + STi2 + ')';
        i1 := round(GetInput(1, 24, St1, St, i0+1, i2));
        ClrEoL;
        Write('                Extending the table . . . ');
        Build(i1)
      end;
    Display;
    repeat
      Ch := ReadKey;
      if Ch <> #0
      then FuncKey := False
      else
        begin
          FuncKey := True;
          Ch := ReadKey
        end;
      InputOk := False;
      if FuncKey and (Ch = #73) and (n0 > 1)
      then InputOk := True;
      if FuncKey and (Ch = #81) and (n0 + 19 < i0)
      then InputOk := True;
      if (not FuncKey) and (UpCase(Ch) = 'X') and (i0 > 20)
      then InputOk := True;
      if (not FuncKey) and (UpCase(Ch) = 'E') and (i0 < 1000)
      then InputOk := True;
      if (not FuncKey) and (Ch = #27)
      then InputOk := true;
      if PrinterOn and (not FuncKey) and (UpCase(Ch) = 'P')
      then InputOk := True
    until InputOk
  until (not FuncKey) and (Ch = #27);       {Esc was pressed}
  TextMode(OriginalMode)
end.