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))