MODULE MagicConstant;

FROM STextIO IMPORT WriteString, WriteLn;
FROM SWholeIO IMPORT WriteInt;
FROM RealMath IMPORT exp, ln;

VAR
  N: CARDINAL;
  E: REAL;

(* Returns the magic constant of a magic square of order N + 2 *)
PROCEDURE A(N: CARDINAL): CARDINAL;
VAR
  N2: CARDINAL;
BEGIN
  N2 := N + 2;
  RETURN (N2 * ((N2 * N2) + 1)) DIV 2
END A;

(* Returns the order of the magic square whose magic constant is at least X *)
PROCEDURE InvA(X: REAL): CARDINAL;
BEGIN
  RETURN VAL(INTEGER, exp(ln((2. * X)) / 3.) + 1.)
  (* Use of power(2. * X, 1. / 3.) + 1. does not give enough precision due to rounded exponent *)
END InvA;

BEGIN
  WriteString("The first 20 magic constants are ");
  FOR N := 1 TO 20 DO
    WriteInt(A(N), 1);
    WriteString(" ")
  END;
  WriteLn;
  WriteString("The 1,000th magic constant is ");
  WriteInt(A(1000), 1);
  WriteLn;
  E := 1.;
  FOR N := 1 TO 20 DO (* The results are equal to these from Mathematica for N <= 23 *)
    E := E * 10.;
    WriteString("10^");
    WriteInt(N, 2);
    WriteString(": ");
    WriteInt(InvA(E), 9);
    WriteLn
  END
END MagicConstant.