A magic square is a square grid containing consecutive integers from 1 to N², arranged so that every row, column and diagonal adds up to the same number. That number is a constant. There is no way to create a valid N x N magic square that does not sum to the associated constant.

Task
Magic constant
You are encouraged to solve this task according to the task description, using any language you may know.
EG

A 3 x 3 magic square always sums to 15.

    ┌───┬───┬───┐
    │ 2 │ 7 │ 6 │
    ├───┼───┼───┤
    │ 9 │ 5 │ 1 │
    ├───┼───┼───┤
    │ 4 │ 3 │ 8 │
    └───┴───┴───┘

A 4 x 4 magic square always sums to 34.

Traditionally, the sequence leaves off terms for n = 0 and n = 1 as the magic squares of order 0 and 1 are trivial; and a term for n = 2 because it is impossible to form a magic square of order 2.


Task
  • Starting at order 3, show the first 20 magic constants.
  • Show the 1000th magic constant. (Order 1003)
  • Find and show the order of the smallest N x N magic square whose constant is greater than 10¹ through 10¹⁰.


Stretch
  • Find and show the order of the smallest N x N magic square whose constant is greater than 10¹¹ through 10²⁰.


See also


11l

Translation of: Python
F a(=n)
   n += 2
   R n * (n ^ 2 + 1) / 2

F inv_a(x)
   V k = 0
   L k * (k ^ 2.0 + 1) / 2 + 2 < x
      k++
   R k

print(‘The first 20 magic constants are:’)
L(n) 1..19
   print(Int(a(n)), end' ‘ ’)
print("\nThe 1,000th magic constant is: "Int(a(1000)))

L(e) 1..19
   print(‘10^’e‘: ’inv_a(10.0 ^ e))
Output:
The first 20 magic constants are:
15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 
The 1,000th magic constant is: 503006505
10^1: 3
10^2: 6
10^3: 13
10^4: 28
10^5: 59
10^6: 126
10^7: 272
10^8: 585
10^9: 1260
10^10: 2715
10^11: 5849
10^12: 12600
10^13: 27145
10^14: 58481
10^15: 125993
10^16: 271442
10^17: 584804
10^18: 1259922
10^19: 2714418

Ada

-- Magic constants
-- J. Carter     2024 May

with Ada.Text_IO;
with System;

procedure Magic_Constant is
   type Big_U is mod System.Max_Binary_Modulus;

   function Magic (Order : in Big_U) return Big_U with
      Pre => Order > 2;
   -- Returns the constant for the magic square of order Order

   function Order (Guess : in Big_U) return Big_U;
   -- Returns the order of the smallest magic square with constant > Guess

   function Image (N : in Big_U; Width : in Positive) return String;
   -- Returns a blank-filled image of N of at least width characters

   function Magic (Order : in Big_U) return Big_U is
      (Order * (Order ** 2 + 1) / 2);

   function Order (Guess : in Big_U) return Big_U is
      Min  : Big_U :=         3;
      Max  : Big_U := 5_850_000;
      Mid  : Big_U;
      Prev : Big_U := 0;
   begin -- Order
      Search : loop
         Mid := Min + (Max - Min) / 2;
         
         exit Search when Mid = Prev;
         
         Prev := Mid;
         
         if Magic (Mid) > Guess and (Mid = 3 or else Magic (Mid - 1) <= Guess) then
            return Mid;
         end if;
         
         if Magic (Mid) > Guess then
            Max := Mid;
         else
            Min := Mid;
         end if;
      end loop Search;

      raise Program_Error with "Order: no solution for" & Guess'Image;
   end Order;

   function Image (N : in Big_U; Width : in Positive) return String is
      Raw : String renames N'Image;
      Img : String renames Raw (2 .. Raw'Last);
   begin -- Image
      return (1 .. Width - Img'Length => ' ') & Img;
   end Image;
begin -- Magic_Constant
   First_20 : for N in Big_U range 3 .. 22 loop
      Ada.Text_IO.Put_Line (Item => Image (N, 2) & Image (Magic (N), 5) );
   end loop First_20;

   Ada.Text_IO.New_Line;
   Ada.Text_IO.Put_Line (Item => "1000" & Magic (1000)'Image);
   Ada.Text_IO.New_Line;

   Ten_To : for P in 1 .. (if Big_U'Size < 64 then 9 elsif Big_U'Size < 128 then 18 else 20) loop
      Ada.Text_IO.Put_Line (Item => "10 ** " & Image (Big_U (P), 2) & Image (Order (10 ** P), 8) );
   end loop Ten_To;
end Magic_Constant;
Output:
 3   15
 4   34
 5   65
 6  111
 7  175
 8  260
 9  369
10  505
11  671
12  870
13 1105
14 1379
15 1695
16 2056
17 2465
18 2925
19 3439
20 4010
21 4641
22 5335

1000 500000500

10 **  1       3
10 **  2       6
10 **  3      13
10 **  4      28
10 **  5      59
10 **  6     126
10 **  7     272
10 **  8     585
10 **  9    1260
10 ** 10    2715
10 ** 11    5849
10 ** 12   12600
10 ** 13   27145
10 ** 14   58481
10 ** 15  125993
10 ** 16  271442
10 ** 17  584804
10 ** 18 1259922
10 ** 19 2714418
10 ** 20 5848036

ALGOL 68

Translation of: FreeBasic

... with the inverse magic constant routine as in the Julia/Wren samples.

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32

Uses ALGOL 68G's LONG LONG INT whose default precision is large enough to cope with 10^20.

BEGIN # find some magic constants - the row, column and diagonal sums of a magin square #
      # translation of the Free Basic sample with the Julia/Wren inverse function #
    # returns the magic constant of a magic square of order n + 2 #
    PROC a     = ( INT n )LONG LONG INT:
         BEGIN
            LONG LONG INT n2 = n + 2;
            ( n2 * ( ( n2 * n2 ) + 1 ) ) OVER 2
         END # a # ;
    # returns the order of the magic square whose magic constant is at least x #
    PROC inv a = ( LONG LONG INT x )LONG LONG INT:
         ENTIER long long exp( long long ln( x * 2 ) / 3 ) + 1;

    print( ( "The first 20 magic constants are " ) );
    FOR n TO 20 DO
        print( ( whole( a( n ), 0 ), " " ) )
    OD;
    print( ( newline ) );
    print( ( "The 1,000th magic constant is ", whole( a( 1000 ), 0 ), newline ) );
    LONG LONG INT e := 1;
    FOR n TO 20 DO
        e *:= 10;
        print( ( "10^", whole( n, -2 ), ": ", whole( inv a( e ), -9 ), newline ) )
    OD
END
Output:
The first 20 magic constants are 15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335 
The 1,000th magic constant is 503006505
10^ 1:         3
10^ 2:         6
10^ 3:        13
10^ 4:        28
10^ 5:        59
10^ 6:       126
10^ 7:       272
10^ 8:       585
10^ 9:      1260
10^10:      2715
10^11:      5849
10^12:     12600
10^13:     27145
10^14:     58481
10^15:    125993
10^16:    271442
10^17:    584804
10^18:   1259922
10^19:   2714418
10^20:   5848036

Arturo

a: function [n][
    n: n+2
    return (n*(1 + n^2))/2
]

aInv: function [x][
    k: new 0
    while [x > 2 + k*(1+k^2)/2]
        -> inc 'k
    return k
]
print "The first 20 magic constants are:"
print map 1..19 => a

print ""
print "The 1,000th magic constant is:"
print a 1000

print ""
loop 1..19 'z ->
    print ["10 ^" z "=>" aInv 10^z]
Output:
The first 20 magic constants are:
15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 

The 1,000th magic constant is:
503006505

10 ^ 1 => 3 
10 ^ 2 => 6 
10 ^ 3 => 13 
10 ^ 4 => 28 
10 ^ 5 => 59 
10 ^ 6 => 126 
10 ^ 7 => 272 
10 ^ 8 => 585 
10 ^ 9 => 1260 
10 ^ 10 => 2715 
10 ^ 11 => 5849 
10 ^ 12 => 12600 
10 ^ 13 => 27145 
10 ^ 14 => 58481 
10 ^ 15 => 125993 
10 ^ 16 => 271442 
10 ^ 17 => 584804 
10 ^ 18 => 1259922 
10 ^ 19 => 2714418

AWK

# syntax: GAWK -f MAGIC_CONSTANT.AWK
# converted from FreeBASIC
BEGIN {
    printf("The first 20 magic constants are:")
    for (i=1; i<=20; i++) {
      printf(" %d",a(i))
    }
    printf("\n")
    printf("The 1,000th magic constant is: %d\n",a(1000))
    for (i=1; i<=20; i++) {
      printf("10^%02d: %8d\n",i,inv_a(10^i))
    }
    exit(0)
}
function a(n) {
    n += 2
    return(n*(n^2+1)/2)
}
function inv_a(x,  k) {
    while (k*(k^2+1)/2+2 < x) {
      k++
    }
    return(k)
}
Output:
The first 20 magic constants are: 15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335
The 1,000th magic constant is: 503006505
10^01:        3
10^02:        6
10^03:       13
10^04:       28
10^05:       59
10^06:      126
10^07:      272
10^08:      585
10^09:     1260
10^10:     2715
10^11:     5849
10^12:    12600
10^13:    27145
10^14:    58481
10^15:   125993
10^16:   271442
10^17:   584804
10^18:  1259922
10^19:  2714418
10^20:  5848036

Basic

FreeBASIC

function a(byval n as uinteger) as ulongint
    n+=2
    return n*(n^2 + 1)/2
end function

function inv_a(x as double) as ulongint
    dim as ulongint k = 0
    while k*(k^2+1)/2+2 < x
        k+=1
    wend
    return k
end function

dim as ulongint n
print "The first 20 magic constants are ":
for n = 1 to 20
    print a(n);" ";
next n
print
print "The 1,000th magic constant is ";a(1000)

for e as uinteger = 1 to 20
    print using "10^##: #########";e;inv_a(10^cast(double,e))
next e
Output:

The first 20 magic constants are 15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335 The 1,000th magic constant is 503006505 10^ 1: 3 10^ 2: 6 10^ 3: 13 10^ 4: 28 10^ 5: 59 10^ 6: 126 10^ 7: 272 10^ 8: 585 10^ 9: 1260 10^10: 2715 10^11: 5849 10^12: 12600 10^13: 27145 10^14: 58481 10^15: 125993 10^16: 271442 10^17: 584804 10^18: 1259922 10^19: 2714418 10^20: 5848036

QB64

$NOPREFIX

DIM order AS INTEGER
DIM target AS INTEGER64

PRINT "First 20 magic constants:"
FOR i = 3 TO 22
    PRINT USING "####,  "; MagicSum(i);
    IF i MOD 5 = 2 THEN PRINT
NEXT i
PRINT
PRINT USING "1000th magic constant: ##########,"; MagicSum(1002)
PRINT
PRINT "Smallest order magic square with a constant greater than:"
FOR i = 1 TO 13 ' 64-bit integers can take us no further, unsigned or not
    target = 10 ^ i
    DO
        order = order + 1
    LOOP UNTIL MagicSum(order) > target
    PRINT USING "10^**: #####,"; i; order
    order = order * 2 - 1
NEXT i

FUNCTION MagicSum&& (n AS INTEGER)
    MagicSum&& = (n * n + 1) / 2 * n
END FUNCTION
Output:
First 20 magic constants:
   15     34     65    111    175
  260    369    505    671    870
1,105  1,379  1,695  2,056  2,465
2,925  3,439  4,010  4,641  5,335

1000th magic constant: 503,006,505

Smallest order magic square with a constant greater than:
10^*1:      3
10^*2:      6
10^*3:     13
10^*4:     28
10^*5:     59
10^*6:    126
10^*7:    272
10^*8:    585
10^*9:  1,260
10^10:  2,715
10^11:  5,849
10^12: 12,600
10^13: 27,145

BASIC256

function a(n)
	n = n + 2
	return n*(n^2 + 1)/2
end function

function inv_a(x)
	k = 0
	while k*(k^2+1)/2+2 < x
		k += 1
	end while
	return k
end function

print "The first 20 magic constants are:"
for n = 1 to 20
	print int(a(n));" ";
next n
print : print
print "The 1,000th magic constant is "; int(a(1000)); chr(10)

for e = 1 to 20
	print "10^"; e; ": "; chr(9); inv_a(10^e)
next e
end

PureBasic

Procedure.i a(n.i)
  n + 2
  ProcedureReturn n*(Pow(n,2) + 1)/2
EndProcedure

Procedure.i inv_a(x.i)
  k.i = 0
  While k*(Pow(k,2)+1)/2+2 < x
    k + 1
  Wend
  ProcedureReturn k
EndProcedure

OpenConsole()
PrintN("The first 20 magic constants are:")
For n.i = 1 To 20
  Print(Str(a(n)) + " ")
Next n
PrintN("") : PrintN("")
PrintN("The 1,000th magic constant is " + Str(a(1000)))

For e.i = 1 To 20
  PrintN("10^" + Str(e) + ": " + #TAB$ + Str(inv_a(Pow(10,e))))
Next e
CloseConsole()

QBasic

Works with: QBasic version 1.1
Works with: QuickBasic version 4.5
FUNCTION a (n)
    n = n + 2
    a = n * (n ^ 2 + 1) / 2
END FUNCTION

FUNCTION inva (x)
    k = 0
    WHILE k * (k ^ 2 + 1) / 2 + 2 < x
        k = k + 1
    WEND
    inva = k
END FUNCTION

PRINT "The first 20 magic constants are: ";
FOR n = 1 TO 20
    PRINT a(n); " ";
NEXT n
PRINT
PRINT "The 1,000th magic constant is "; a(1000)
PRINT
FOR e = 1 TO 20
    PRINT USING "10^##: #########"; e; inva(10 ^ e)
NEXT e
END

True BASIC

FUNCTION a(n)
    LET n = n + 2
    LET a = n*(n^2 + 1)/2
END FUNCTION

FUNCTION inv_a(x)
    LET k = 0
    DO WHILE k*(k^2+1)/2+2 < x
       LET k = k + 1
    LOOP
    LET inv_a = k
END FUNCTION

PRINT "The first 20 magic constants are: ";
FOR n = 1 TO 20
    PRINT a(n);" ";
NEXT n
PRINT
PRINT "The 1,000th magic constant is "; a(1000)

FOR e = 1 TO 20
    PRINT USING "10^##": e;
    PRINT USING": #########": inv_a(10^e)
NEXT e
END

Yabasic

sub a(n)
    n = n + 2
    return n*(n^2 + 1)/2
end sub
 
sub inv_a(x)
    k = 0
    while k*(k^2+1)/2+2 < x
        k = k + 1
    wend
    return k
end sub
 
print "The first 20 magic constants are: "
for n = 1 to 20
    print a(n), " ";
next n
print "\nThe 1,000th magic constant is ", a(1000), "\n"
 
for e = 1 to 20
    print "10^", e using"##", ": ", inv_a(10^e) using "#########"
next e
end


C#

Translation of: Java
using System;

public class MagicConstant {

    private const int OrderFirstMagicSquare = 3;

    public static void Main(string[] args) {
        Console.WriteLine("The first 20 magic constants:");
        for (int i = 1; i <= 20; i++) {
            Console.Write(" " + MagicConstantValue(Order(i)));
        }
        Console.WriteLine("\n");

        Console.WriteLine("The 1,000th magic constant: " + MagicConstantValue(Order(1_000)) + "\n");

        Console.WriteLine("Order of the smallest magic square whose constant is greater than:");
        for (int i = 1; i <= 20; i++) {
            string powerOf10 = "10^" + i + ":";
            Console.WriteLine($"{powerOf10,6}{MinimumOrder(i),8}");
        }
    }

    // Return the magic constant for a magic square of the given order 
    private static int MagicConstantValue(int n) {
        return n * (n * n + 1) / 2;
    }

    // Return the smallest order of a magic square such that its magic constant is greater than 10 to the given power
    private static int MinimumOrder(int n) {
        return (int)Math.Exp((Math.Log(2.0)+ n * Math.Log(10.0)) / 3) + 1;
    }

    // Return the order of the magic square at the given index
    private static int Order(int index) {
        return OrderFirstMagicSquare + index - 1;
    }
}
Output:
The first 20 magic constants:
 15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335

The 1,000th magic constant: 503006505

Order of the smallest magic square whose constant is greater than:
 10^1:       3
 10^2:       6
 10^3:      13
 10^4:      28
 10^5:      59
 10^6:     126
 10^7:     272
 10^8:     585
 10^9:    1260
10^10:    2715
10^11:    5849
10^12:   12600
10^13:   27145
10^14:   58481
10^15:  125993
10^16:  271442
10^17:  584804
10^18: 1259922
10^19: 2714418
10^20: 5848036

C++

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <string>

constexpr int32_t ORDER_FIRST_MAGIC_SQUARE = 3;
constexpr double LN2 = log(2.0);
constexpr double LN10 = log(10.0);

// Return the magic constant for a magic square of the given order
int32_t magicConstant(int32_t n) {
	return n * ( n * n + 1 ) / 2;
}

// Return the smallest order of a magic square such that its magic constant is greater than 10 to the given power
int32_t minimumOrder(int32_t n) {
	return (int) exp( ( LN2 + n * LN10 ) / 3 ) + 1;
}

// Return the order of the magic square at the given index
int32_t order(int32_t index) {
	return ORDER_FIRST_MAGIC_SQUARE + index - 1;
}

int main() {
	std::cout << "The first 20 magic constants:" << std::endl;
	for ( int32_t i = 1; i <= 20; ++i ) {
		std::cout << " " << magicConstant(order(i));
	}
	std::cout << std::endl << std::endl;

	std::cout << "The 1,000th magic constant: " << magicConstant(order(1'000)) << std::endl << std::endl;

	std::cout << "Order of the smallest magic square whose constant is greater than:" << std::endl;
	for ( int32_t i = 1; i <= 20; ++i ) {
		std::string power_of_10 = "10^" + std::to_string(i) + ":";
		std::cout << std::setw(6) << power_of_10 << std::setw(8) << minimumOrder(i) << std::endl;
	}
}
Output:
The first 20 magic constants:
 15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335

The 1,000th magic constant: 503006505

Order of the smallest magic square whose constant is greater than:
 10^1:       3
 10^2:       6
 10^3:      13
 10^4:      28
 10^5:      59
 10^6:     126
 10^7:     272
 10^8:     585
 10^9:    1260
10^10:    2715
10^11:    5849
10^12:   12600
10^13:   27145
10^14:   58481
10^15:  125993
10^16:  271442
10^17:  584804
10^18: 1259922
10^19: 2714418
10^20: 5848036

Delphi

Works with: Delphi version 6.0


function GetMagicNumber(N: double): double;
begin
Result:=N * (((N * N) + 1) / 2);
end;

function GetNumberLess(N: double): integer;
var M: double;
begin
for Result:=1 to High(Integer) do
	begin
	M:=GetMagicNumber(Result);
	if M>N then break;
	end;
end;

procedure ShowMagicNumber(Memo: TMemo);
var I,J: integer;
var N,M: double;
var S: string;
begin
S:='';
for I:=3 to 23 do
	begin
	S:=S+Format('%8.0n',[GetMagicNumber(I)]);
	if (I mod 5)=2 then S:=S+#$0D#$0A;
	end;
Memo.Lines.Add(S);
Memo.Lines.Add('');
Memo.Lines.Add('1000th: '+Format('%8.0n',[GetMagicNumber(1002)]));
Memo.Lines.Add('');
N:=10;
for I:=1 to 20 do
	begin
	J:=GetNumberLess(N);
	Memo.Lines.Add('M^'+Format('%d%8d',[I,J]));
	N:=N * 10;
	end;
end;
Output:
      15      34      65     111     175
     260     369     505     671     870
   1,105   1,379   1,695   2,056   2,465
   2,925   3,439   4,010   4,641   5,335
   6,095

1000th: 503,006,505

M^1       3
M^2       6
M^3      13
M^4      28
M^5      59
M^6     126
M^7     272
M^8     585
M^9    1260
M^10    2715
M^11    5849
M^12   12600
M^13   27145
M^14   58481
M^15  125993
M^16  271442
M^17  584804
M^18 1259922
M^19 2714418
M^20 5848036


EasyLang

func a n .
   n += 2
   return n * (n * n + 1) / 2
.
func inva x .
   while k * (k * k + 1) / 2 + 2 < x
      k += 1
   .
   return k
.
write "The first 20 magic constants: "
for n to 20
   write a n & " "
.
print ""
print ""
print "The 1,000th magic constant: " & a 1000
print ""
print "Smallest magic square with constant greater than:"
for e to 10
   print "10^" & e & ": " & inva pow 10 e
.

Factor

Works with: Factor version 0.99 2021-06-02
USING: formatting io kernel math math.functions.integer-logs
math.ranges prettyprint sequences ;

: magic ( m -- n ) dup sq 1 + 2 / * ;

"First 20 magic constants:" print
3 22 [a,b] [ bl ] [ magic pprint ] interleave nl
nl
"1000th magic constant: " write 1002 magic .
nl
"Smallest order magic square with a constant greater than:" print
1 0 20 [
    [ 10 * ] dip
    [ dup magic pick < ] [ 1 + ] while
    over integer-log10 over "10^%02d: %d\n" printf
    dup + 1 -
] times 2drop
Output:
First 20 magic constants:
15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335

1000th magic constant: 503006505

Smallest order magic square with a constant greater than:
10^01: 3
10^02: 6
10^03: 13
10^04: 28
10^05: 59
10^06: 126
10^07: 272
10^08: 585
10^09: 1260
10^10: 2715
10^11: 5849
10^12: 12600
10^13: 27145
10^14: 58481
10^15: 125993
10^16: 271442
10^17: 584804
10^18: 1259922
10^19: 2714418
10^20: 5848036

FutureBasic

Translation of: FreeBASIC
long def fn MagicConstant( n as long ) = n*(n^2+1)/2

UInt64 local fn InvA( x as double )
  UInt64 k = 0
  while ( k*(k^2+1)/2+2 < x )
    k++
  wend
end fn = k

void local fn DoIt
  window 1, @"Magic constant", (0,0,640,400)
  print @"First 20 magic constants:"
  long n = 3
  while n <= 1002
    if ( n <= 22 )
      print fn MagicConstant( n );@" ";
    end if
    if ( n == 1002 ) then print @"\n\n1000th magic constant: ";fn MagicConstant(1002)
    n++
  wend
  
  print
  
  for long e = 1 to 20
    printf @"10^%ld: %ld",e,fn InvA(10^e)
  next
end fn

fn DoIt

HandleEvents
Output:
First 20 magic constants:
15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335 

1000th magic constant: 503006505

10^1: 3
10^2: 6
10^3: 13
10^4: 28
10^5: 59
10^6: 126
10^7: 272
10^8: 585
10^9: 1260
10^10: 2715
10^11: 5849
10^12: 12600
10^13: 27145
10^14: 58481
10^15: 125993
10^16: 271442
10^17: 584804
10^18: 1259922
10^19: 2714418
10^20: 5848036

Go

Translation of: Wren
Library: Go-rcu
package main

import (
    "fmt"
    "math"
    "rcu"
)

func magicConstant(n int) int {
    return (n*n + 1) * n / 2
}

var ss = []string{
    "\u2070", "\u00b9", "\u00b2", "\u00b3", "\u2074",
    "\u2075", "\u2076", "\u2077", "\u2078", "\u2079",
}

func superscript(n int) string {
    if n < 10 {
        return ss[n]
    }
    if n < 20 {
        return ss[1] + ss[n-10]
    }
    return ss[2] + ss[0]
}

func main() {
    fmt.Println("First 20 magic constants:")
    for n := 3; n <= 22; n++ {
        fmt.Printf("%5d ", magicConstant(n))
        if (n-2)%10 == 0 {
            fmt.Println()
        }
    }

    fmt.Println("\n1,000th magic constant:", rcu.Commatize(magicConstant(1002)))

    fmt.Println("\nSmallest order magic square with a constant greater than:")
    for i := 1; i <= 20; i++ {
        goal := math.Pow(10, float64(i))
        order := int(math.Cbrt(goal*2)) + 1
        fmt.Printf("10%-2s : %9s\n", superscript(i), rcu.Commatize(order))
    }
}
Output:
Same as Wren example.

J

Implementation:

mgc=: 0 0.5 0 0.5&p.

In other words, the magic constant for a magic square of order x is the result of the polynomial (0.5*x)+(0.5*x^3)

Task examples:

   mgc 3+i.20
15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335
   mgc 1003x
504514015
   (#\,.],.mgc) x:(mgc i.3000) I.10^1+i.10 
 1    3          15
 2    6         111
 3   13        1105
 4   28       10990
 5   59      102719
 6  126     1000251
 7  272    10061960
 8  585   100101105
 9 1260  1000188630
10 2715 10006439295

stretch example:

   ((10+#\),.],.mgc) x:(mgc i.6e6) I.10^11+i.10 
11    5849          100049490449
12   12600         1000188006300
13   27145        10000910550385
14   58481       100003310078561
15  125993      1000021311323825
16  271442     10000026341777165
17  584804    100000232056567634
18 1259922   1000002262299152685
19 2714418  10000004237431278525
20 5848036 100000026858987459346

Java

public final class MagicConstant {

	public static void main(String[] aArgs) {
		System.out.println("The first 20 magic constants:");
		for ( int i = 1; i <= 20; i++ ) {
			System.out.print(" " + magicConstant(order(i)));
		}
		System.out.println(System.lineSeparator());
		
		System.out.println("The 1,000th magic constant: " + magicConstant(order(1_000)) + System.lineSeparator());
		
		System.out.println("Order of the smallest magic square whose constant is greater than:");
		for ( int i = 1; i <= 20; i++ ) {
			String powerOf10 = "10^" + i + ":";
			System.out.println(String.format("%6s%8s", powerOf10, minimumOrder(i)));
		}
	}
	
	// Return the magic constant for a magic square of the given order 
	private static int magicConstant(int aN) {	  
		return aN * ( aN * aN + 1 ) / 2;
	}

	// Return the smallest order of a magic square such that its magic constant is greater than 10 to the given power
	private static int minimumOrder(int aN) {		
		return (int) Math.exp( ( LN2 + aN * LN10 ) / 3 ) + 1;
	}
	
	// Return the order of the magic square at the given index
	private static int order(int aIndex) {
		return ORDER_FIRST_MAGIC_SQUARE + aIndex - 1;
	}
	
    private static final int ORDER_FIRST_MAGIC_SQUARE = 3;
	private static final double LN2 = Math.log(2.0);
	private static final double LN10 = Math.log(10.0);

}
Output:
The first 20 magic constants:
 15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335

The 1,000th magic constant: 503006505

Order of the smallest magic square whose constant is greater than:
 10^1:       3
 10^2:       6
 10^3:      13
 10^4:      28
 10^5:      59
 10^6:     126
 10^7:     272
 10^8:     585
 10^9:    1260
10^10:    2715
10^11:    5849
10^12:   12600
10^13:   27145
10^14:   58481
10^15:  125993
10^16:  271442
10^17:  584804
10^18: 1259922
10^19: 2714418
10^20: 5848036

jq

Adapted from Wren
Works with jq (*)
Works with gojq, the Go implementation of jq

(*) The arithmetic precision of the C implementation of jq is insufficient for the extended task.

Preliminaries

# To take advantage of gojq's arbitrary-precision integer arithmetic:
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);

# nth-root
def iroot($n):
  . as $in
  | if $n == 1 then .
    else (. < 0) as $neg
    | if $neg and (n % 2) == 0
      then "Cannot take the \($n)th root of a negative number." | error
      else ($n-1) as $n
      | {t: (if $neg then -. else . end)}
      | .s = .t + 1
      | .u = .t
      | until (.u >= .s;
          .s = .u
          | .u = ((.u * $n) + (.t / (.u|power($n)))) / ($n + 1) )
      | if $neg then - .s else .s end
      end
    end;

# input: an array
# output: a stream of arrays of size size except possibly for the last array
def group(size):
  recurse( .[size:]; length>0) | .[0:size];

def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;

def ss : ["\u2070", "\u00b9", "\u00b2", "\u00b3", "\u2074", 
          "\u2075", "\u2076", "\u2077", "\u2078", "\u2079"];
 
def superscript:
  if . < 10 then ss[.]
  elif . < 20 then ss[1] + ss[. - 10]
  else ss[2] + ss[0]
  end;

The Task

def magicConstant: (.*. + 1) * . / 2;
 
"First 20 magic constants:",
 ([range(3;23)
  | magicConstant]
  | group(10) | map(lpad(5)) | join(" ")),
 "",
 "1,000th magic constant: \( 1002| magicConstant)",
 "",
 "Smallest order magic square with a constant greater than:",
 (range(1; 21) as $i
  | (10 | power($i)) as $goal
  | ((($goal * 2)|iroot(3) + 1) | floor) as $order
  | ("10\($i|superscript)" | lpad(5)) + ": \($order|lpad(9))"  )
Output:

As for Wren except for commatization.

Julia

Uses the inverse of the magic constant function for the last part of the task.

using Lazy

magic(x) = (1 + x^2) * x ÷ 2
magics = @>> Lazy.range() map(magic) filter(x -> x > 10) # first 2 values are filtered out
println("First 20 magic constants: ", Int.(take(20, magics)))
println("Thousandth magic constant is: ", collect(take(1000, magics))[end])

println("Smallest magic square with constant greater than:")
for expo in 1:20
    goal = big"10"^expo
    ordr = Int(floor((2 * goal)^(1/3))) + 1
    println("10^", string(expo, pad=2), "    ", ordr)
end
Output:
First 20 magic constants: [15, 34, 65, 111, 175, 260, 369, 505, 671, 870, 1105, 1379, 1695, 2056, 2465, 2925, 3439, 4010, 4641, 5335]
Thousandth magic constant is: 503006505
Smallest magic square with constant greater than:
10^01    3
10^02    6
10^03    13
10^04    28
10^05    59
10^06    126
10^07    272
10^08    585
10^09    1260
10^10    2715
10^11    5849
10^12    12600
10^13    27145
10^14    58481
10^15    125993
10^16    271442
10^17    584804
10^18    1259922
10^19    2714418
10^20    5848036

Lua

function magic (x)
    return x * (1 + x^2) / 2
end

print("Magic constants of orders 3 to 22:")
for i = 3, 22 do
    io.write(magic(i) .. " ")
end

print("\n\nMagic constant 1003: " .. magic(1003) .. "\n")

print("Orders of smallest magic constant greater than...")
print("-----\t-----\nValue\tOrder\n-----\t-----")
local order = 1
for i = 1, 20 do
    repeat 
        order = order + 1
    until magic(order) > 10 ^ i
    print("10^" .. i, order)
end
Output:
Magic constants of orders 3 to 22:
15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335

Magic constant 1003: 504514015

Orders of smallest magic constant greater than...
-----   -----
Value   Order
-----   -----
10^1    3
10^2    6
10^3    13
10^4    28
10^5    59
10^6    126
10^7    272
10^8    585
10^9    1260
10^10   2715
10^11   5849
10^12   12600
10^13   27145
10^14   58481
10^15   125993
10^16   271442
10^17   584804
10^18   1259922
10^19   2714418
10^20   5848036

Mathematica /Wolfram Language

ClearAll[i, n, MagicSumHelper, MagicSum, InverseMagicSum]
MagicSumHelper[n_] = Sum[i, {i, n^2}]/n;
MagicSum[n_] := MagicSumHelper[n + 2]
InverseMagicSum[lim_] := Ceiling[-(1/(3^(1/3) (9 lim + Sqrt[3] Sqrt[1 + 27 lim^2])^(1/3))) + (9 lim + Sqrt[3] Sqrt[1 + 27 lim^2])^(1/3)/3^(2/3)]

MagicSum /@ Range[20]
MagicSum[1000]

exps = Range[1, 50];
nums = 10^exps;
Transpose[{Superscript[10, #] & /@ exps, InverseMagicSum[nums]}] // Grid
Output:
{15, 34, 65, 111, 175, 260, 369, 505, 671, 870, 1105, 1379, 1695, 2056, 2465, 2925, 3439, 4010, 4641, 5335}
503006505

10^1	3
10^2	6
10^3	13
10^4	28
10^5	59
10^6	126
10^7	272
10^8	585
10^9	1260
10^10	2715
10^11	5849
10^12	12600
10^13	27145
10^14	58481
10^15	125993
10^16	271442
10^17	584804
10^18	1259922
10^19	2714418
10^20	5848036
10^21	12599211
10^22	27144177
10^23	58480355
10^24	125992105
10^25	271441762
10^26	584803548
10^27	1259921050
10^28	2714417617
10^29	5848035477
10^30	12599210499
10^31	27144176166
10^32	58480354765
10^33	125992104990
10^34	271441761660
10^35	584803547643
10^36	1259921049895
10^37	2714417616595
10^38	5848035476426
10^39	12599210498949
10^40	27144176165950
10^41	58480354764258
10^42	125992104989488
10^43	271441761659491
10^44	584803547642574
10^45	1259921049894874
10^46	2714417616594907
10^47	5848035476425733
10^48	12599210498948732
10^49	27144176165949066
10^50	58480354764257322

Nim

import std/[math, unicode]

func magicConstant(n: int): int =
  ## Return the magic constant for a magic square of order "n".
  n * (n * n + 1) div 2

func minOrder(n: int): int =
  ## Return the smallest order such as the magic constant is greater than "10^n".
  const Ln2 = ln(2.0)
  const Ln10 = ln(10.0)
  result = int(exp((Ln2 + n.toFloat * Ln10) / 3)) + 1

const First = 3
const Superscripts: array['0'..'9', string] = ["⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹"]

template order(idx: Positive): int =
  ## Compute the order of the magic square at index "idx".
  idx + (First - 1)

func superscript(n: Natural): string =
  ## Return the Unicode string to use to represent an exponent.
  for d in $n:
    result.add(Superscripts[d])

echo "First 20 magic constants:"
for idx in 1..20:
  stdout.write ' ', order(idx).magicConstant
echo()

echo "\n1000th magic constant: ", order(1000).magicConstant

echo "\nOrder of the smallest magic square whose constant is greater than:"
for n in 1..20:
  let left = "10" & n.superscript & ':'
  echo left.alignLeft(6), ($minOrder(n)).align(7)
Output:
First 20 magic constants:
 15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335

1000th magic constant: 503006505

Order of the smallest magic square whose constant is greater than:
10¹:        3
10²:        6
10³:       13
10⁴:       28
10⁵:       59
10⁶:      126
10⁷:      272
10⁸:      585
10⁹:     1260
10¹⁰:    2715
10¹¹:    5849
10¹²:   12600
10¹³:   27145
10¹⁴:   58481
10¹⁵:  125993
10¹⁶:  271442
10¹⁷:  584804
10¹⁸: 1259922
10¹⁹: 2714418
10²⁰: 5848036

Pascal

Free Pascal

program MagicConst;
{$IFDEF FPC}{$MODE DELPHI}{$OPTIMIZATION ON,ALL}{$ENDIF}
{$IFDEF WINDOWS}{$APPTYPE CONSOLE}{$ENDIF}

function MagicSum(n :Uint32):Uint64; inline;
var
  k : Uint64;
begin
  k := n*Uint64(n);
  result := (k*k+k) DIV 2;
end;

function MagSumPerRow(n:Uint32):Uint32;
begin
//result := MagicSum(n) DIV n;
  //(n^3 + n) /2
  result := ((Uint64(n)*n+1)*n) DIV 2;
end;
var
  s : String[31];
  i : Uint32;
  lmt,rowcnt : extended;
Begin
  writeln('First Magic constants 3..20');
  For i := 3 to 20 do
     write(MagSumPerRow(i),' ');
  writeln;
  writeln('1000.th ',MagSumPerRow(1002));

  writeln('First Magic constants > 10^xx');
  //lmt = (rowcnt^3 + rowcnt) /2 -> rowcnt > (lmt*2 )^(1/3)
  lmt := 2.0 * 10.0;
  For i :=  1 to 50 do
  begin
    rowcnt := Int(exp(ln(lmt)/3))+1.0;//+1 suffices 
    str(trunc(rowcnt),s);
    writeln('10^',i:2,#9,s:18);
    f := 10.0*lmt;
  end;
end.
Output:
First Magic constants 3..20
15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010
1000.th 503006505
First Magic constants > 10^xx
10^ 1	                 3
10^ 2	                 6
10^ 3	                13
10^ 4	                28
10^ 5	                59
10^ 6	               126
10^ 7	               272
10^ 8	               585
10^ 9	              1260
10^10	              2715
10^11	              5849
10^12	             12600
10^13	             27145
10^14	             58481
10^15	            125993
10^16	            271442
10^17	            584804
10^18	           1259922
10^19	           2714418
10^20	           5848036
... same as Mathematica
10^46	  2714417616594907
10^47	  5848035476425733
10^48	 12599210498948732
10^49	 27144176165949066
10^50	 58480354764257322

Perl

#!/usr/bin/perl

use strict; # https://rosettacode.org/wiki/Magic_constant
use warnings;

my @twenty = map $_ * ( $_ ** 2 + 1 ) / 2, 3 .. 22;
print "first twenty: @twenty\n\n" =~ s/.{50}\K /\n/gr;

my $thousandth = 1002 * ( 1002 ** 2 + 1 ) / 2;
print "thousandth: $thousandth\n\n";

print "10**N   order\n";
for my $i ( 1 .. 20 )
  {
  printf "%3d %9d\n", $i, (10 ** $i * 2) ** ( 1 / 3 ) + 1;
  }
Output:
first twenty: 15 34 65 111 175 260 369 505 671 870
1105 1379 1695 2056 2465 2925 3439 4010 4641 5335

thousandth: 503006505

10**N   order
  1         3
  2         6
  3        13
  4        28
  5        59
  6       126
  7       272
  8       585
  9      1260
 10      2715
 11      5849
 12     12600
 13     27145
 14     58481
 15    125993
 16    271442
 17    584804
 18   1259922
 19   2714418
 20   5848036

Phix

with javascript_semantics

function magic(integer nth)
    integer order = nth+2
    return (order*order+1)/2 * order
end function
printf(1,"First 20 magic constants: %V\n",{apply(tagset(20),magic)})
printf(1,"1000th magic constant: %,d\n",{magic(1000)})

include mpfr.e

mpz {goal, order} = mpz_inits(2)
for i=1 to 20 do
    mpz_ui_pow_ui(goal,10,i)
    mpz_mul_si(order,goal,2)
    mpz_nthroot(order,order,3)
    mpz_add_si(order,order,1)
    printf(1,"1e%d: %s\n",{i,mpz_get_str(order,10,true)})
end for
Output:
First 20 magic constants: {15,34,65,111,175,260,369,505,671,870,1105,1379,1695,2056,2465,2925,3439,4010,4641,5335}
1000th magic constant: 503,006,505
1e1: 3
1e2: 6
1e3: 13
1e4: 28
1e5: 59
1e6: 126
1e7: 272
1e8: 585
1e9: 1,260
1e10: 2,715
1e11: 5,849
1e12: 12,600
1e13: 27,145
1e14: 58,481
1e15: 125,993
1e16: 271,442
1e17: 584,804
1e18: 1,259,922
1e19: 2,714,418
1e20: 5,848,036

Prolog

Minimalistic, efficient approach

m(X,Y):- Y is X*(X*X+1)/2.

l(L,R,T,X):- L > R -> X is L; M is div(L+R,2), m(M,F),
    (T < F -> R_ is M-1, l(L,R_,T,X); L_ is M+1, l(L_,R,T,X)).
l(B,X):- l(1,B,B,X).

task:-
    write("First 20 magic constants are:"), forall(between(3,22,N), (m(N,X), format(" ~d",X))), nl,
    write("The 1000th magic constant is:"), forall(m(1002,X), format(" ~d",X)), nl,
    forall(between(1,20,N), (l(10**N,X), format("10^~d:\t~d\n",[N,X]))).
Output:
?- task. 
First 20 magic constants are: 15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335
The 1000th magic constant is: 503006505
10^1:   3
10^2:   6
10^3:   13
10^4:   28
10^5:   59
10^6:   126
10^7:   272
10^8:   585
10^9:   1260
10^10:  2715
10^11:  5849
10^12:  12600
10^13:  27145
10^14:  58481
10^15:  125993
10^16:  271442
10^17:  584804
10^18:  1259922
10^19:  2714418
10^20:  5848036
true.


Python

#!/usr/bin/python

def a(n):
    n += 2
    return n*(n**2 + 1)/2
 
def inv_a(x):
    k = 0
    while k*(k**2+1)/2+2 < x:
        k+=1
    return k

 
if __name__ == '__main__':
    print("The first 20 magic constants are:");
    for n in range(1, 20):
        print(int(a(n)), end = " ");
    print("\nThe 1,000th magic constant is:",int(a(1000)));
     
    for e in range(1, 20):
        print(f'10^{e}: {inv_a(10**e)}');
Output:
The first 20 magic constants are:
15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 
The 1,000th magic constant is: 503006505
10^1: 3
10^2: 6
10^3: 13
10^4: 28
10^5: 59
10^6: 126
10^7: 272
10^8: 585
10^9: 1260
10^10: 2715
10^11: 5849
10^12: 12600
10^13: 27145
10^14: 58481
10^15: 125993
10^16: 271442
10^17: 584804
10^18: 1259922
10^19: 2714418

Quackery

  [ 3 + dup 3 ** + 2 / ] is magicconstant ( n --> n )
 
  20 times [ i^ magicconstant echo sp ] cr cr
 
  1000 magicconstant echo cr cr
 
  0 1 
  [ over magicconstant over > if
      [ over 3 + echo cr
        10 * ] 
    dip 1+
    [ 10 21 ** ] constant 
    over = until ] 
  2drop
Output:
15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335 

504514015

3
4
6
13
28
59
126
272
585
1260
2715
5849
12600
27145
58481
125993
271442
584804
1259922
2714418
5848036

Raku

use Lingua::EN::Numbers:ver<2.8+>;

my @magic-constants = lazy (3..∞).hyper.map: { (1 + .²) * $_ / 2 };

put "First 20 magic constants: ", @magic-constants[^20]».&comma;
say "1000th magic constant: ", @magic-constants[999].&comma;
say "\nSmallest order magic square with a constant greater than:";

(1..20).map: -> $p {printf "10%-2s: %s\n", $p.&super, comma 3 + @magic-constants.first( * > exp($p, 10), :k ) }
Output:
First 20 magic constants: 15 34 65 111 175 260 369 505 671 870 1,105 1,379 1,695 2,056 2,465 2,925 3,439 4,010 4,641 5,335
1000th magic constant: 503,006,505

Smallest order magic square with a constant greater than:
10¹ : 3
10² : 6
10³ : 13
10⁴ : 28
10⁵ : 59
10⁶ : 126
10⁷ : 272
10⁸ : 585
10⁹ : 1,260
10¹⁰: 2,715
10¹¹: 5,849
10¹²: 12,600
10¹³: 27,145
10¹⁴: 58,481
10¹⁵: 125,993
10¹⁶: 271,442
10¹⁷: 584,804
10¹⁸: 1,259,922
10¹⁹: 2,714,418
10²⁰: 5,848,036

RPL

This task can be solved through a few one-liners, by using algebraic and equation-solving features.

First, let's store the equation of a(n):

'X*(SQ(X)+1)/2' 'A6003' STO

Then, evaluate it for n=3 to 22 to get the first 20 magic constants:

≪ {} 3 22 FOR j j 'X' STO A6003 EVAL + NEXT ≫ EVAL

We need now to define the inverse function of a(n):

A6003 OVER - 'X' ROT ROOT CEIL ≫ 'A6003→' STO

And finally, look for a-1( 10 ^ j ), for j=1 to 10:

≪ {} 1 10 FOR j 10 j ^ A6003→ + NEXT ≫ EVAL

which is a one-minute job for a basic HP-28S.

Output:
2: { 15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335 }
1: { 3  6  13  28 59 126 272 585 1260 2715 }

Sidef

func f(n) {
    (n+2) * ((n+2)**2 + 1) / 2
}

func order(n) {
    iroot(2*n, 3) + 1
}

say ("First 20 terms: ", f.map(1..20).join(' '))
say ("1000th term: ", f(1000), " with order ", order(f(1000)))

for n in (1 .. 20) {
    printf("order(10^%-2s) = %s\n", n, order(10**n))
}
Output:
First 20 terms: 15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335
1000th term: 503006505 with order 1003
order(10^1 ) = 3
order(10^2 ) = 6
order(10^3 ) = 13
order(10^4 ) = 28
order(10^5 ) = 59
order(10^6 ) = 126
order(10^7 ) = 272
order(10^8 ) = 585
order(10^9 ) = 1260
order(10^10) = 2715
order(10^11) = 5849
order(10^12) = 12600
order(10^13) = 27145
order(10^14) = 58481
order(10^15) = 125993
order(10^16) = 271442
order(10^17) = 584804
order(10^18) = 1259922
order(10^19) = 2714418
order(10^20) = 5848036

Wren

Library: Wren-seq
Library: Wren-fmt

This uses Julia's approach for the final parts.

import "./seq" for Lst
import "./fmt" for Fmt
 
var magicConstant = Fn.new { |n| (n*n + 1) * n / 2 }
 
var ss = ["\u2070", "\u00b9", "\u00b2", "\u00b3", "\u2074", 
          "\u2075", "\u2076", "\u2077", "\u2078", "\u2079"]
 
var superscript = Fn.new { |n| (n < 10) ? ss[n] : (n < 20) ? ss[1] + ss[n - 10] : ss[2] + ss[0] }
 
System.print("First 20 magic constants:")
var mc20 = (3..22).map { |n| magicConstant.call(n) }.toList
for (chunk in Lst.chunks(mc20, 10)) Fmt.print("$5d", chunk)
 
Fmt.print("\n1,000th magic constant: $,d", magicConstant.call(1002))
 
System.print("\nSmallest order magic square with a constant greater than:")
for (i in 1..20) {
    var goal = 10.pow(i)
    var order = (goal * 2).cbrt.floor + 1
    Fmt.print("10$-2s : $,9d", superscript.call(i), order)
}
Output:
First 20 magic constants:
   15    34    65   111   175   260   369   505   671   870
 1105  1379  1695  2056  2465  2925  3439  4010  4641  5335

1,000th magic constant: 503,006,505

Smallest order magic square with a constant greater than:
10¹  :         3
10²  :         6
10³  :        13
10⁴  :        28
10⁵  :        59
10⁶  :       126
10⁷  :       272
10⁸  :       585
10⁹  :     1,260
10¹⁰ :     2,715
10¹¹ :     5,849
10¹² :    12,600
10¹³ :    27,145
10¹⁴ :    58,481
10¹⁵ :   125,993
10¹⁶ :   271,442
10¹⁷ :   584,804
10¹⁸ : 1,259,922
10¹⁹ : 2,714,418
10²⁰ : 5,848,036

XPL0

A magic square of side N contains N^2 items. The sum of a sequence 1..N^2 is given by: Sum = (N^2+1) * N^2 / 2. A grid row adds to the magic constant, and N rows add to the Sum. Thus the magic constant = Sum/N = (N^3+N)/2.

int  N, X;
real M, Thresh, MC;
[Text(0, "First 20 magic constants:^M^J");
for N:= 3 to 20+3-1 do
    [IntOut(0, (N*N*N+N)/2);  ChOut(0, ^ )];
CrLf(0);
Text(0, "1000th magic constant: ");
N:= 1000+3-1;
IntOut(0, (N*N*N+N)/2);
CrLf(0);
Text(0, "Smallest order magic square with a constant greater than:^M^J");
Thresh:= 10.;
M:= 3.;
Format(1, 0);
for X:= 1 to 10 do
    [repeat MC:= (M*M*M+M)/2.;
            M:= M+1.;
    until   MC > Thresh;
    Text(0, "10^^");
    if X < 10 then ChOut(0, ^0);
    IntOut(0, X);
    Text(0, ": ");
    RlOut(0, M-1.);
    CrLf(0);
    Thresh:= Thresh*10.;
    ];
]
Output:
First 20 magic constants:
15 34 65 111 175 260 369 505 671 870 1105 1379 1695 2056 2465 2925 3439 4010 4641 5335 
1000th magic constant: 503006505
Smallest order magic square with a constant greater than:
10^01: 3
10^02: 6
10^03: 13
10^04: 28
10^05: 59
10^06: 126
10^07: 272
10^08: 585
10^09: 1260
10^10: 2715