615 lines
22 KiB
Fortran
615 lines
22 KiB
Fortran
!
|
|
! Cyclops Numbers
|
|
! tested with Intel ifx (IFX) 2025.2.1 20250806 on Kubuntu 25.10
|
|
! GNU Fortran (Ubuntu 15.2.0-4ubuntu4) 15.2.0 on Kubuntu 25.10
|
|
! VSI Fortran x86-64 V8.6-001 on OpenVMS x86_64 V9.2-3
|
|
! No Non-standard features used, should compile on any fairly recent Fortran.
|
|
! Original Pascal code converted to Fortran
|
|
! U.B., October 2025
|
|
!==============================================================================
|
|
|
|
module CyclopsNumbers
|
|
implicit none
|
|
|
|
! Constants
|
|
integer, parameter :: BIGLIMIT = 10000000 ! Upper limit for number search
|
|
integer, parameter :: MAX_DIGITS = 10 ! Maximum digits in base-9 representation
|
|
|
|
! Type for numbers in base-9 representation
|
|
type Num9
|
|
integer :: nmdgts(0:MAX_DIGITS) ! Digits in base-9 (0-8)
|
|
integer :: nmMaxDgtIdx ! Highest digit index used
|
|
integer(kind=8) :: nmNum ! Actual numerical value
|
|
end type Num9
|
|
|
|
! Type for cyclops numbers (left + zero + right parts)
|
|
type CyclopsNum
|
|
type(Num9) :: cnRight, cnLeft ! Right and left parts around zero
|
|
integer(kind=8) :: cnNum ! Complete cyclops number
|
|
integer :: cndigits ! Number of digits in each side
|
|
integer :: cnIdx ! Index of cyclops number
|
|
end type CyclopsNum
|
|
|
|
! Global arrays for precomputed values
|
|
integer(kind=8) :: cnMin(0:15) ! Minimum values for digit counts
|
|
integer(kind=8) :: cnPow10Shift(0:15) ! Powers of 10 for shifting
|
|
integer(kind=8) :: cnPow9(0:15) ! Powers of 9 for indexing
|
|
|
|
! Dynamic array to store found cyclops numbers
|
|
integer(kind=8), allocatable :: Cyclops(:)
|
|
|
|
contains
|
|
|
|
!============================================================================
|
|
! Initialize precomputed arrays for minimum values and powers
|
|
! cnMin: (0, 1, 11, 111, 1111, ...) - minimum values for n-digit numbers
|
|
! cnPow10Shift: (100, 1000, 10000, ...) - powers of 10 for position shifting
|
|
! cnPow9: (1, 9, 81, 729, ...) - powers of 9 for base conversion
|
|
!============================================================================
|
|
subroutine InitCnMinPow()
|
|
integer :: i
|
|
integer(kind=8) :: min_val, pow, pow9
|
|
|
|
min_val = 0
|
|
pow = 100
|
|
pow9 = 1
|
|
|
|
do i = 0, 15
|
|
cnMin(i) = min_val
|
|
min_val = 10 * min_val + 1 ! Build numbers like 1, 11, 111, ...
|
|
cnPow10Shift(i) = pow
|
|
pow = pow * 10 ! Next power of 10
|
|
cnPow9(i) = pow9
|
|
pow9 = pow9 * 9 ! Next power of 9
|
|
end do
|
|
end subroutine InitCnMinPow
|
|
|
|
!============================================================================
|
|
! Reset a Num9 type to minimum value for given digit count
|
|
! @param tn: Num9 object to clear
|
|
! @param idx: Digit count index
|
|
!============================================================================
|
|
subroutine ClearNum9 (tn, idx)
|
|
type(Num9), intent(out) :: tn
|
|
integer, intent(in) :: idx
|
|
integer :: i
|
|
|
|
tn%nmdgts = 0 ! Clear all digits
|
|
tn%nmMaxDgtIdx = 0 ! Reset max digit index
|
|
tn%nmNum = cnMin(idx + 1) ! Set to minimum value
|
|
end subroutine ClearNum9
|
|
|
|
!============================================================================
|
|
! Initialize a cyclops number to starting state
|
|
! @param cn: CyclopsNum object to initialize
|
|
!============================================================================
|
|
subroutine InitCycNum(cn)
|
|
type(CyclopsNum), intent(out) :: cn
|
|
|
|
cn%cndigits = 0 ! Start with 0 digits each side
|
|
call ClearNum9(cn%cnLeft, 0) ! Clear left part
|
|
call ClearNum9(cn%cnRight, 0) ! Clear right part
|
|
cn%cnNum = 0 ! Zero value
|
|
cn%cnIdx = 0 ! Start index
|
|
end subroutine InitCycNum
|
|
|
|
!============================================================================
|
|
! Increment a base-9 number, handling carry-over
|
|
! @param tn: Num9 object to increment
|
|
!============================================================================
|
|
subroutine IncNum9(tn)
|
|
type(Num9), intent(inout) :: tn
|
|
integer :: idx, fac, n, i
|
|
|
|
idx = 0
|
|
fac = 1
|
|
n = tn%nmdgts(0) + 1 ! Increment least significant digit
|
|
tn%nmNum = tn%nmNum + 1 ! Increment numerical value
|
|
|
|
! Handle carry-over propagation
|
|
do
|
|
if (n < 9) exit ! No carry-over needed
|
|
|
|
! Carry-over occurs
|
|
tn%nmNum = tn%nmNum + fac ! Adjust numerical value
|
|
tn%nmdgts(idx) = 0 ! Reset current digit
|
|
idx = idx + 1 ! Move to next digit
|
|
fac = fac * 10 ! Increase factor for numerical adjustment
|
|
|
|
if (idx > MAX_DIGITS) exit ! Check bounds
|
|
n = tn%nmdgts(idx) + 1 ! Get next digit + 1
|
|
end do
|
|
|
|
! Store final digit value
|
|
if (idx <= MAX_DIGITS) then
|
|
tn%nmdgts(idx) = n
|
|
if (tn%nmMaxDgtIdx < idx) then
|
|
tn%nmMaxDgtIdx = idx ! Update max digit index if needed
|
|
end if
|
|
end if
|
|
end subroutine IncNum9
|
|
|
|
!============================================================================
|
|
! Generate the next cyclops number in sequence
|
|
! @param cycnum: Current cyclops number to advance
|
|
!============================================================================
|
|
subroutine NextCycNum(cycnum)
|
|
type(CyclopsNum), intent(inout) :: cycnum
|
|
|
|
if (cycnum%cnIdx /= 0) then
|
|
! Normal case: increment from previous number
|
|
call IncNum9(cycnum%cnRight) ! Increment right part
|
|
|
|
! Check if right part overflowed (needs more digits)
|
|
if (cycnum%cnRight%nmMaxDgtIdx > cycnum%cndigits) then
|
|
call ClearNum9(cycnum%cnRight, cycnum%cndigits) ! Reset right to minimum
|
|
call IncNum9(cycnum%cnLeft) ! Increment left part
|
|
|
|
! Check if left part overflowed (needs more digits)
|
|
if (cycnum%cnLeft%nmMaxDgtIdx > cycnum%cndigits) then
|
|
cycnum%cndigits = cycnum%cndigits + 1 ! Increase digit count
|
|
call ClearNum9(cycnum%cnLeft, cycnum%cndigits) ! Reset left
|
|
call ClearNum9(cycnum%cnRight, cycnum%cndigits) ! Reset right
|
|
if (cycnum%cndigits > MAX_DIGITS) then
|
|
cycnum%cndigits = MAX_DIGITS ! Enforce maximum
|
|
end if
|
|
end if
|
|
end if
|
|
|
|
! Reconstruct complete cyclops number: left * 10^(digits+1) + right
|
|
cycnum%cnNum = cycnum%cnLeft%nmNum * cnPow10Shift(cycnum%cndigits) + cycnum%cnRight%nmNum
|
|
cycnum%cnIdx = cycnum%cnIdx + 1
|
|
else
|
|
! Special case: first cyclops number
|
|
cycnum%cnNum = 101
|
|
cycnum%cnIdx = 1
|
|
end if
|
|
end subroutine NextCycNum
|
|
|
|
|
|
!============================================================================
|
|
! Make a complete Palindromic Cyclops number when left half is already setup.
|
|
! Argument is already initialized to have its left half setup as needed.
|
|
! Fill digits of right half in reversed order.
|
|
! @param cycnum: half-filled Cyclops number to be completed as palindrome.
|
|
!============================================================================
|
|
|
|
subroutine MakePalinCycNum (cycnum)
|
|
type(CyclopsNum), intent(inout) :: cycnum
|
|
|
|
! make right to have reversed digits of left
|
|
integer :: n,dgt
|
|
integer :: i,j
|
|
|
|
n = 0
|
|
i = 0
|
|
do j = cycnum%cnDigits, 0 , -1
|
|
dgt = cycnum%cnLeft%nmdgts(i)
|
|
cycnum%cnRight%nmdgts(j) = dgt
|
|
n = 10*n+(dgt+1)
|
|
i = i + 1
|
|
end do
|
|
cycnum%cnRight%nmNum = n
|
|
cycnum%cnNum = cycnum%cnLeft%nmNum*cnPow10Shift(cycnum%cndigits)+n
|
|
end subroutine MakePalinCycNum
|
|
|
|
|
|
|
|
!============================================================================
|
|
! Increment left half of a Cyclops number, handling carry-over
|
|
!============================================================================
|
|
|
|
subroutine IncLeftCn( cn )
|
|
|
|
type(CyclopsNum), intent(inout) :: cn
|
|
|
|
! set right digits to minimum
|
|
call ClearNum9 (cn%cnRight,cn%cndigits)
|
|
! increment left digits
|
|
call IncNum9 (cn%cnLeft)
|
|
! One more digit in left half?
|
|
if (cn%cnLeft%nmMaxDgtIdx > cn%cndigits) then ! left has more digits than before?
|
|
cn%cndigits = cn%cndigits + 1 ! Increment for entire Cyclops number
|
|
call ClearNum9(cn%cnLeft,cn%cndigits) ! Set left to minimum for new number of digits
|
|
call ClearNum9(cn%cnRight,cn%cndigits) ! Set right to minimum for new number of digits
|
|
if (cn%cndigits>MAX_DIGITS) then
|
|
cn%cndigits = MAX_DIGITS
|
|
endif
|
|
endif
|
|
cn%cnNum = cn%cnLeft%nmNum*cnPow10Shift(cn%cndigits)+cn%cnRight%nmNUm
|
|
end subroutine IncLeftCn
|
|
|
|
!============================================================================
|
|
! Convert a zero-based index to the corresponding cyclops number
|
|
! Uses combinatorial counting based on digit positions in base-9
|
|
! @param n: Zero-based index
|
|
! @return: Corresponding cyclops number
|
|
!============================================================================
|
|
function IndexToCyclops(n) result(result_cn)
|
|
integer(kind=8), intent(in) :: n
|
|
type(CyclopsNum) :: result_cn
|
|
integer :: dgtCnt, i
|
|
integer(kind=8) :: p9, q, num, local_n
|
|
|
|
call InitCycNum(result_cn)
|
|
if (n == 0) return ! Index 0 returns initialized (0)
|
|
|
|
result_cn%cnIdx = n
|
|
dgtCnt = 0
|
|
local_n = n
|
|
|
|
! Determine number of digits needed
|
|
do
|
|
p9 = cnPow9(dgtCnt) * cnPow9(dgtCnt) ! Combinations for dgtCnt digits
|
|
if (local_n < p9) exit ! Found correct digit count
|
|
local_n = local_n - p9 ! Subtract counted combinations
|
|
dgtCnt = dgtCnt + 1 ! Try next digit count
|
|
if (dgtCnt > 10) exit ! Safety limit
|
|
end do
|
|
|
|
dgtCnt = dgtCnt - 1 ! Adjust to actual digit count
|
|
|
|
! Process right digits (convert from base-9 to base-10)
|
|
result_cn%cnRight%nmMaxDgtIdx = dgtCnt
|
|
do i = 0, dgtCnt
|
|
q = local_n / 9 ! Integer division
|
|
result_cn%cnRight%nmdgts(i) = local_n - 9 * q ! Remainder gives base-9 digit
|
|
local_n = q ! Continue with quotient
|
|
end do
|
|
|
|
! Convert base-9 digits to base-10 number (adding 1 to each digit)
|
|
num = 0
|
|
do i = dgtCnt, 0, -1 ! Most significant digit first
|
|
num = num * 10 + result_cn%cnRight%nmdgts(i) + 1
|
|
end do
|
|
result_cn%cnRight%nmNum = num
|
|
result_cn%cnNum = num ! Temporary: only right part
|
|
|
|
! Process left digits (same algorithm)
|
|
result_cn%cnLeft%nmMaxDgtIdx = dgtCnt
|
|
do i = 0, dgtCnt
|
|
q = local_n / 9
|
|
result_cn%cnLeft%nmdgts(i) = local_n - 9 * q
|
|
local_n = q
|
|
end do
|
|
|
|
! Convert left digits to base-10 number
|
|
num = 0
|
|
do i = dgtCnt, 0, -1
|
|
num = num * 10 + result_cn%cnLeft%nmdgts(i) + 1
|
|
end do
|
|
result_cn%cnLeft%nmNum = num
|
|
|
|
! Combine left and right with zero in middle: left * 10^(dgtCnt+1) + right
|
|
result_cn%cnNum = result_cn%cnNum + num * cnPow10Shift(dgtCnt)
|
|
result_cn%cndigits = dgtCnt
|
|
end function IndexToCyclops
|
|
|
|
!============================================================================
|
|
! Display cyclops numbers in formatted columns
|
|
! @param cl: Array of cyclops numbers
|
|
! @param colw: Column width for formatting
|
|
! @param colc: Number of columns per row
|
|
!============================================================================
|
|
subroutine Out_Cyclops(cl, colw, colc)
|
|
integer(kind=8), intent(in) :: cl(:)
|
|
integer, intent(in) :: colw, colc
|
|
integer :: i, n
|
|
character (len=50) :: colFormat
|
|
|
|
n = size(cl)
|
|
if (n > 100) n = 100
|
|
|
|
! Create suitable format string to nicely printout columns using colw as width
|
|
write (colformat, '("(I", I0, " )" )') colw
|
|
|
|
! Print all "n=size(cl)" Numbers of array cl
|
|
do i = 1, n
|
|
write(*, colformat, advance='no') cl(i) ! Print number using calculated format string
|
|
if (mod(i, colc) == 0) then
|
|
write(*, *) ! New line after colc numbers
|
|
else
|
|
write(*, '(A)', advance='no') ' ' ! Space between numbers
|
|
end if
|
|
end do
|
|
|
|
! Terminate line if incomplete
|
|
if (mod(i, colc) /= 0) write(*, *)
|
|
end subroutine Out_Cyclops
|
|
|
|
!============================================================================
|
|
! Check if a number is prime (simple trial division)
|
|
! @param n: Number to check
|
|
! @return: .true. if prime, .false. otherwise
|
|
!============================================================================
|
|
logical function isPrime(n)
|
|
integer(kind=8), intent(in) :: n
|
|
integer(kind=8) :: p, q
|
|
|
|
! Handle small numbers and simple cases
|
|
if (n <= 1) then
|
|
isPrime = .false.
|
|
return
|
|
end if
|
|
|
|
if (n == 2 .or. n == 3 .or. n == 5) then
|
|
isPrime = .true.
|
|
return
|
|
end if
|
|
|
|
! Even numbers and multiples of 3
|
|
if (mod(n, 2) == 0 .or. mod(n, 3) == 0) then
|
|
isPrime = .false.
|
|
return
|
|
end if
|
|
|
|
! Trial division for remaining numbers
|
|
p = 5
|
|
do
|
|
q = n / p
|
|
if (n - q * p == 0) then ! Divisible by p
|
|
isPrime = .false.
|
|
return
|
|
end if
|
|
p = p + 2 ! Next odd number
|
|
|
|
q = n / p
|
|
if (n - q * p == 0) then ! Divisible by p+2
|
|
isPrime = .false.
|
|
return
|
|
end if
|
|
|
|
if (q < p) exit ! No need to check further
|
|
p = p + 4 ! Skip to next candidate
|
|
end do
|
|
|
|
isPrime = .true. ! No divisors found
|
|
end function isPrime
|
|
|
|
!============================================================================
|
|
! Generate the first 'cnt' basic cyclops numbers
|
|
! @param cnt: Number of cyclops numbers to generate
|
|
! @return: Last cyclops number generated (for continuation)
|
|
!============================================================================
|
|
function FirstCyclops(cnt) result(result_cn)
|
|
integer, intent(in) :: cnt
|
|
type(CyclopsNum) :: result_cn
|
|
integer :: i
|
|
|
|
! Allocate storage for results
|
|
if (allocated(Cyclops)) deallocate(Cyclops)
|
|
allocate(Cyclops(cnt))
|
|
|
|
i = 1
|
|
call InitCycNum(result_cn)
|
|
|
|
! Generate requested number of cyclops numbers
|
|
do while (i <= cnt)
|
|
Cyclops(i) = result_cn%cnNum
|
|
i = i + 1
|
|
call NextCycNum(result_cn)
|
|
end do
|
|
|
|
! Continue until beyond BIGLIMIT for reporting
|
|
do while (result_cn%cnNum <= BIGLIMIT)
|
|
call NextCycNum(result_cn)
|
|
end do
|
|
end function FirstCyclops
|
|
|
|
!============================================================================
|
|
! Generate the first 'cnt' prime cyclops numbers
|
|
! @param cnt: Number of prime cyclops numbers to generate
|
|
! @return: Last prime cyclops number with updated index
|
|
! Algorithm as in function FirstCyclops, but check for each generated
|
|
! Cyclops number if its prime.
|
|
!============================================================================
|
|
function FirstPrimeCyclops(cnt) result(result_cn)
|
|
integer, intent(in) :: cnt
|
|
type(CyclopsNum) :: result_cn
|
|
integer :: i
|
|
|
|
if (allocated(Cyclops)) deallocate(Cyclops)
|
|
allocate(Cyclops(cnt))
|
|
|
|
i = 1
|
|
call InitCycNum(result_cn)
|
|
|
|
do while (i <= cnt)
|
|
if (isPrime(result_cn%cnNum)) then ! only consider prime Cyclops numbers
|
|
Cyclops(i) = result_cn%cnNum
|
|
i = i + 1
|
|
end if
|
|
call NextCycNum(result_cn)
|
|
end do
|
|
|
|
|
|
! Continue until prime Cyclops number is beyond BIGLIMIT for reporting
|
|
do
|
|
if (isPrime(result_cn%cnNum)) then
|
|
result_cn%cnIdx = i ! result_cn%cnIdx + 1
|
|
if (result_cn%cnNum .gt. BIGLIMIT) exit
|
|
i = i + 1
|
|
end if
|
|
call NextCycNum(result_cn)
|
|
end do
|
|
end function FirstPrimeCyclops
|
|
|
|
|
|
|
|
!============================================================================
|
|
! Generate the first 'cnt' prime blind cyclops numbers
|
|
! @param cnt: Number of prime cyclops numbers to generate
|
|
! @return: Last prime cyclops number with updated index
|
|
! Algorithm as in function FirstCyclops, but check for each generated
|
|
! Cyclops number if its prime, then construct blind number from this, and
|
|
! check again if result is prime as well.
|
|
!============================================================================
|
|
function FirstPrimeBlindCyclops(cnt) result(result_cn)
|
|
integer, intent(in) :: cnt
|
|
type(CyclopsNum) :: result_cn
|
|
integer :: i
|
|
integer (kind=8) :: n
|
|
|
|
if (allocated(Cyclops)) deallocate(Cyclops)
|
|
allocate(Cyclops(cnt))
|
|
|
|
i = 1
|
|
call InitCycNum(result_cn)
|
|
|
|
do while (i <= cnt)
|
|
if (isPrime(result_cn%cnNum)) then
|
|
|
|
n = result_cn%cnRight%nmnum
|
|
if (result_cn%cndigits .gt. 0) then
|
|
n = n + result_cn%cnLeft%nmnum*cnPow10Shift(result_cn%cndigits-1)
|
|
else
|
|
n = n + result_cn%cnLeft%nmnum*10
|
|
endif
|
|
if (isPrime (n)) then
|
|
Cyclops(i) =result_cn%cnNum
|
|
i = i + 1
|
|
endif
|
|
end if
|
|
call NextCycNum(result_cn)
|
|
end do
|
|
|
|
do
|
|
if (isPrime(result_cn%cnNum)) then
|
|
n = result_cn%cnRight%nmnum
|
|
if (result_cn%cndigits .gt. 0) then
|
|
n = n + result_cn%cnLeft%nmnum*cnPow10Shift(result_cn%cndigits-1)
|
|
else
|
|
n = n + result_cn%cnLeft%nmnum*10
|
|
endif
|
|
if (isPrime (n)) then
|
|
result_cn%cnIdx = i
|
|
if (result_cn%cnNum .gt. BIGLIMIT) exit
|
|
i = i + 1
|
|
end if
|
|
end if
|
|
call NextCycNum(result_cn)
|
|
end do
|
|
end function FirstPrimeBlindCyclops
|
|
|
|
!============================================================================
|
|
! Generate the first 'cnt' prime palindromic cyclops numbers
|
|
! @param cnt: Number of prime palindromic cyclops numbers to generate
|
|
! @return: Last prime cyclops number with updated index
|
|
! Algorithm: directly create palindromic cyclops numbers and look for primes.
|
|
!============================================================================
|
|
function FirstPrimePalindromicCyclops(cnt) result(result_cn)
|
|
integer, intent(in) :: cnt
|
|
type(CyclopsNum) :: result_cn
|
|
integer :: i
|
|
integer (kind=8) :: n
|
|
|
|
if (allocated(Cyclops)) deallocate(Cyclops)
|
|
allocate(Cyclops(cnt))
|
|
|
|
i = 1
|
|
call InitCycNum(result_cn)
|
|
|
|
do while (i <= cnt)
|
|
call MakePalinCycNum (result_cn) ! Complete prepare Cyclops number to a palindrome
|
|
if (isPrime(result_cn%cnNum)) then ! Consider only if its prime.
|
|
Cyclops(i) =result_cn%cnNum
|
|
i = i + 1
|
|
end if
|
|
call IncLeftCn (result_cn) ! Next one.
|
|
|
|
! Last digit cannot be even because we need only prime numbers
|
|
do while (result_cn%cnLeft%nmdgts(result_cn%cnDigits)+1 .ne. 1 .AND. &
|
|
result_cn%cnLeft%nmdgts(result_cn%cnDigits)+1 .ne. 3 .and. &
|
|
result_cn%cnLeft%nmdgts(result_cn%cnDigits)+1 .ne. 5 .and. &
|
|
result_cn%cnLeft%nmdgts(result_cn%cnDigits)+1 .ne. 7 .and. &
|
|
result_cn%cnLeft%nmdgts(result_cn%cnDigits)+1 .ne. 9)
|
|
call IncLeftCn(result_cn)
|
|
end do
|
|
end do
|
|
|
|
do
|
|
call MakePalinCycNum (result_cn)
|
|
if (isPrime(result_cn%cnNum)) then
|
|
result_cn%cnIdx = i
|
|
if (result_cn%cnNum .gt. BIGLIMIT) exit
|
|
i = i + 1
|
|
end if
|
|
call IncLeftCn (result_cn)
|
|
|
|
do while (result_cn%cnLeft%nmdgts(result_cn%cnDigits)+1 .ne. 1 .AND. &
|
|
result_cn%cnLeft%nmdgts(result_cn%cnDigits)+1 .ne. 3 .and. &
|
|
result_cn%cnLeft%nmdgts(result_cn%cnDigits)+1 .ne. 5 .and. &
|
|
result_cn%cnLeft%nmdgts(result_cn%cnDigits)+1 .ne. 7 .and. &
|
|
result_cn%cnLeft%nmdgts(result_cn%cnDigits)+1 .ne. 9)
|
|
call IncLeftCn(result_cn)
|
|
end do
|
|
end do
|
|
|
|
end function FirstPrimePalindromicCyclops
|
|
|
|
end module CyclopsNumbers
|
|
|
|
program main
|
|
use CyclopsNumbers
|
|
implicit none
|
|
|
|
type(CyclopsNum) :: cycnum
|
|
integer :: cnt
|
|
integer(kind=8) :: i, large_cnt
|
|
|
|
call InitCnMinPow()
|
|
|
|
! First 50 cyclops numbers
|
|
cnt = 50
|
|
write(*, "( 'The first ', I0, ' cyclops numbers are:')") cnt
|
|
cycnum = FirstCyclops(cnt)
|
|
call Out_Cyclops(Cyclops, 5, 10)
|
|
write(*, "('First cyclops number > ', I0, ' is ', I0, ' at 1-based index ', I0, / )") BIGLIMIT, cycnum%cnNum, cycnum%cnIdx
|
|
|
|
! First 50 prime cyclops numbers
|
|
cnt = 50
|
|
write(*, "('The first ', I0, ' prime cyclops numbers are:')") cnt
|
|
cycnum = FirstPrimeCyclops(cnt)
|
|
call Out_Cyclops(Cyclops, 7, 10)
|
|
write(*, "('First prime cyclops number > ', I0, ' is ', I0, ' at 1-based index ', I0, / )") BIGLIMIT, cycnum%cnNum, cycnum%cnIdx
|
|
write(*, *)
|
|
|
|
|
|
! First 50 prime blind cyclops numbers
|
|
cnt = 50
|
|
write(*, "('The first ', I0, ' prime blind cyclops numbers are:')") cnt
|
|
cycnum = FirstPrimeBlindCyclops(cnt)
|
|
call Out_Cyclops(Cyclops, 7, 10)
|
|
write(*, "('First prime blind cyclops number > ', I0, ' is ', I0, ' at 1-based index ', I0, / )") BIGLIMIT, cycnum%cnNum, cycnum%cnIdx
|
|
|
|
! First 50 prime palindromic cyclops numbers
|
|
cnt = 50
|
|
write(*, "( 'The first ', I0, ' prime palindromic cyclops numbers are:')") cnt
|
|
cycnum = FirstPrimePalindromicCyclops(cnt)
|
|
call Out_Cyclops(Cyclops, 7, 10)
|
|
write(*, "('First palindromic prime cyclops number > ', I0, ' is ', I0, ' at 1-based index ', I0, / )") BIGLIMIT, cycnum%cnNum, cycnum%cnIdx
|
|
|
|
|
|
|
|
|
|
! Additional test cases for large indices
|
|
write(*, '(A/)') 'Demonstrating Cyclops calculation for large indices:'
|
|
write(*, '(A)') ' Index Cyclop Number Calculated Cyclop Number '
|
|
large_cnt = 100
|
|
do while (large_cnt <= 1000000000000000_8)
|
|
write(*, '(I17, A)', advance='no') large_cnt, ' '
|
|
if (large_cnt <= 10000) then ! For small index values evaluate Cyclops number sequentially
|
|
call InitCycNum(cycnum)
|
|
do i = 1, large_cnt-1 ! at index 1 we have 0, is already there by initialisation
|
|
call NextCycNum(cycnum) ! hence we need only (large_cnt-1) next sequence elements
|
|
end do
|
|
write(*, '(I7)', advance='no') cycnum%cnNum
|
|
else
|
|
write (*, '(" ")', advance='no')
|
|
end if
|
|
cycnum = IndexToCyclops(large_cnt-1) ! direct calculation, for comparison when index is small
|
|
write(*, *) ' ', cycnum%cnNum
|
|
large_cnt = large_cnt * 10 ! next order of magnitude
|
|
end do
|
|
|
|
end program main
|