115 lines
4.2 KiB
Fortran
115 lines
4.2 KiB
Fortran
! Magic constant
|
|
! tested with Intel ifx (IFX) 2025.2.1 20250806 on Kubuntu 25.10
|
|
! GNU gfortran (Ubuntu 15.2.0-4ubuntu4) 15.2.0 on Kubuntu 25.10
|
|
! VSI Fortran x86-64 V8.7-001 on OpenVMS V9.2-3
|
|
! U.B., March 2026
|
|
!
|
|
program magicConstant
|
|
implicit none
|
|
integer,parameter :: DP=8, k_int=4
|
|
|
|
real(kind=DP), parameter :: realten=10._DP
|
|
|
|
integer(kind=k_int) :: order, ii
|
|
! Starting at order 3, show the first 20 magic constants.
|
|
|
|
write (*, '("Count Order Magic")')
|
|
order = 3
|
|
do ii=1,20
|
|
write (*,'(i5, 2x, i5,2x,i0)') ii, order, magicNumber (order)
|
|
order = order + 1
|
|
end do
|
|
|
|
! Show the 1000th magic constant. (Order 1003)
|
|
! 1st magic constant is for order 3, i.e. order = count+2
|
|
write (*,*)
|
|
order = 1002
|
|
ii = 1000
|
|
write (*,'(i5, 2x, i5,2x,i0,/)') ii, order, magicNumber (order)
|
|
!
|
|
! Find and show the order of the smallest N x N magic square whose constant is greater than 10^1 through 10^10.
|
|
! Stretch
|
|
! Find and show the order of the smallest N x N magic square whose constant is greater than 10^11 through 10^20.
|
|
!
|
|
write (*, '(" Limit Order")')
|
|
do i=1,20
|
|
order = unmagic (realten**i) ! Using REAL for all 10^i to avoid 132 bit integers
|
|
write (*,'("10**", i2, 2x, i0)') i, order
|
|
end do
|
|
|
|
contains
|
|
|
|
! =============================================
|
|
! Return the magic number of a sqare of order n
|
|
! =============================================
|
|
pure function magicNumber (n) result (retval)
|
|
implicit none
|
|
integer(kind=k_int), intent(in) :: n
|
|
integer(kind=k_int) :: retval
|
|
integer (kind=k_int) :: nn
|
|
nn = n
|
|
retval = (nn**3+nn) / 2
|
|
|
|
end function magicNumber
|
|
|
|
! =====================================================================================
|
|
! return the order of the smallest magic square with a magic number >= a given "magic"
|
|
! =====================================================================================
|
|
function unMagic (magic) result (order)
|
|
implicit none
|
|
real (kind=DP), intent(in) :: magic
|
|
|
|
real (kind=DP) :: guess
|
|
integer(kind=k_int) :: order
|
|
|
|
integer(kind=k_int) :: lo,hi,mid, mnc, mnf
|
|
|
|
! Following formula calculates an approximate value of the order we're looking for.
|
|
! This is because the term n^3 in the formula for the Magic Number is dominant,
|
|
! so the correct (real) value of the magic square's order is most likely between
|
|
! Floor and Ceiling of this first guess, so the integer result is CEILING (guess)
|
|
!
|
|
! It turned out this approximation and using ceiling(this guess) as function
|
|
! value is correct for all test values up to 10^20.
|
|
|
|
guess = (2._DP*magic)**(1._8/3._8)
|
|
|
|
! This check is only here because above reasoning is (maybe) plausible, but it is
|
|
! not a conclusive proof. This function allows us to check whether any results are incorrect.
|
|
if (checkGuess (magic, guess)) then
|
|
order = ceiling (guess, k_int)
|
|
else
|
|
order=-1 ! signal fault
|
|
print *, 'Problem: Guess for inverse magic fails for ', magic
|
|
endif
|
|
|
|
end function unMagic
|
|
|
|
! ========================================================================================================
|
|
! Check if the value of "guess" fulfills the condition
|
|
! magicNumber (floor(guess) <= magic < magicNumber(ceiling(guess)
|
|
! Special care has be taken to avoid integer overflows. Note that argument 'magic' can be as large as 10^20,
|
|
! and 'guess' is (2*guess)^(1/3)), which is .lt. 10^7. So guess , floor(guess) and ceiling(guess)
|
|
! can be represented by a 32-bin integer word, but 'magic' would cause an integer overflow.
|
|
! So we replicate the calculation of the integer function magicNumber(n) here, but we
|
|
! use double precision results instead of integers.
|
|
! ========================================================================================================
|
|
!
|
|
function checkGuess (magic, guess) result (guessIsOK)
|
|
real (kind=DP), intent(in) :: magic, guess
|
|
logical :: guessIsOK
|
|
|
|
real (kind=DP) :: FloorGuess, CeilingGuess
|
|
real (kind=dp) :: CeilingMagic,FloorMagic
|
|
|
|
FloorGuess = real (floor (guess, k_int) , DP)
|
|
FloorMagic = (Floorguess**3 + Floorguess) /2.0_DP
|
|
|
|
CeilingGuess = real (ceiling (guess, k_int) , DP)
|
|
CeilingMagic = (CeilingGuess**3 + Ceilingguess) /2.0_DP
|
|
|
|
guessIsOK = FloorMagic .le. magic .and. magic .lt. CeilingMagic
|
|
|
|
end function checkGuess
|
|
|
|
end program magicConstant
|