It was suggested a long time ago that, instead of scanning a table of fifth powers of integers, one way to answer the question of whether an integer is a fifth power is to evaluate modulo(integer,11).
I think this is stated somewhat incorrectly. It is true that mod(i**5,11) is an element of {0,1,10} for all i , but there are integers such that mod(N,11) is an element of {0,1,10} that are not in the form i**5 . N ={11,12,21} are three examples.
I modified the earlier code based on this mod relation to cycle the outer loops early. Here is the code I added.
...
T = mod( n5(i), 11 ) ! target modulus has values 0, 1, and 10.
...
m11 = mod(sk,11)
if ( exclude_k( m11) ) cycle
...
m11 = mod(sl,11)
if ( exclude_l( m11 ) ) cycle
...
logical function exclude_k( m11 )
! determine if the index k should be excluded.
integer(ik), intent(in) :: m11
select case (T)
case (0); exclude_k = 3 <= m11 .and. m11 <= 8
case (1); exclude_k = 3 <= m11 .and. m11 <= 9
case (10); exclude_k = 2 <= m11 .and. m11 <= 8
end select
return
end function exclude_k
logical function exclude_l( m11 )
! determine if the index l should be excluded.
integer(ik), intent(in) :: m11
select case (T)
case (0); exclude_l = 2 <= m11 .and. m11 <= 9
case (1); exclude_l = 3 <= m11 .and. m11 <= 10
case (10); exclude_l = 1 <= m11 .and. m11 <= 8
end select
return
end function exclude_l
I computed those ranges by hand, so they might not be exactly correct, but I think they are at least close enough to see if they make a timing difference.
Here is how I derived them. Take the target modulus T=0 as an example. The sl partial sum modulus can then only be {0,1,10}. No other values can possibly lead to the target T=0. Backing up one loop level, the sk partial sum modulus can only have the values {0,1,2,9,10}. The inner loops can be skipped for all other sk values. There are no restrictions on the j values, all of them can lead to T=0. I did the same thing for target values of T=1 and T=10.
With these loop restrictions, along with the other checks discussed already above, the total number of innermost loop iterations is reduced from 333794145 to 198143894.
Here are the timings for an Apple M1 cpu on MacOS.
$gfortran -O3 euler.F90
$ time a.out
i^5 = 61917364224
133 110 84 27 144
real 0m0.133s
user 0m0.122s
sys 0m0.002s
That is some 30x faster than the orginal code without any loop exits. I can post the whole code if necessary. Here I just wanted to focus on the modulus relation of @mecej4 to show that effect.
1 Like
Please do post the latest code. As it is, someone who wants to try your code has to merge several code fragments to get a working program. As of now, your program, with the mod 11 additions, is the fastest, so it should be available for future readers.
Similarly, if (j,k,l,m) is not a solution, then x*(j,k,l,m) is also not a solution. Perhaps this could be used to eliminate some tests, although this proportion of future tests may be small ?
I squeezed a bit of performance out of Ron’s latest posting to get Finish 7.00049996E-02 with the following
module timer
integer*8 :: start_clock, clock_rate, end_clock
contains
subroutine start_time
call system_clock ( start_clock, clock_rate )
write (*,*) 'Start ',real(start_clock) / real(clock_rate )
end subroutine start_time
subroutine end_time
call system_clock ( end_clock )
write (*,*) 'Finish',real(end_clock-start_clock) / real(clock_rate )
end subroutine end_time
end module timer
program sum_conjecture
! sum_conjecture.f90
!
! 27^5 + 84^5 + 110^5 + 133^5 = 144^5
!
use, intrinsic :: iso_fortran_env, only: int64
use timer
implicit none
integer, parameter :: nmax = 6000
integer :: j, k, l, n
integer(int64) :: i, n5(nmax), sk, sl, sn
character(*), parameter :: cfmt = '(*(g0,1x))'
call start_time
outer: do i = 1, nmax
n5(i) = i**5
do j = 1, i
do k = 1, j
sk = n5(j) + n5(k)
if ( sk >= n5(i) ) exit
do l = 1, k
sl = sk + n5(l)
if ( sl >= n5(i) ) exit
do n = 1, l
sn = sl + n5(n)
if ( sn < n5(i) ) cycle
if ( sn == n5(i) ) then
print cfmt, "i^5 = ", n5(i)
print cfmt, j, k, l, n, i
exit outer
end if
exit
end do
end do
end do
end do
end do outer
call end_time
end program sum_conjecture
!$ gfortran -O3 sum_conjecture.f90
!$ time a.out
!i^5 = 61917364224
!133 110 84 27 144
!
!real 0m0.188s
!user 0m0.174s
!sys 0m0.007s
PROCESSOR_DESCRIPTION=AMD Ryzen 9 5900X 12-Core Processor 64GB 3200MHz memory
options=-fimplicit-none -O3 -march=native
gcc_path=C:\Program Files (x86)\gcc_eq\gcc_15.1.0\bin;C:\Program Files (x86)\gcc_eq\gcc_15.1.0\libexec\gcc\x86_64-w64-mingw32\15.1.0
Start 14077.1953
i^5 = 61917364224
133 110 84 27 144
Finish 7.00049996E-02
perhaps loop “do n = 1,l” could be limited with ?
n1 = max ( int( real(n5(i)-sl)**0.2 ), 1 )
do n = n1, l
Although this is not effective ?
Another thought on why they used a CDC in 1966, rather than, say, an IBM 7090. The IBMs had 36-bit words, so they could not handle integers such as 144^5. I don’t know if IBMs had multiple precision software support.
No doubt. There was a huge IBM vs CDC rivalry back then - culminating in the famous lawsuit where CDC won a massive settlement against IBM.
Probably the biggest mistake Mr Cray made on the 6600 was not making it a full 64 bit machine. Given the time frame of when it was designed, it was understandable. The CDC 8600 prototype, and then the Cray Research machines corrected that flaw.
OTOH, a 60 bit machine has the interesting quality that a word is divisible by 2, 3, 4, 5, and 6.
I need to get the latest version of Desktop Cyber running on the computer I’m typing this on. It’d be fun to see how fast the code runs. I do have an older version of dtcyber running on a much older and slower machine out in my workshop. (Also have a CDC emulator which I wrote about 25 years ago. It can deadstart and run older CDC OSes. But not newer ones. It also doesn’t 100% implement the floating point instructions. I abandoned it when Tom Hunter wrote dtcyber.)
Ok, here it is.
program sum_conjecture
! sum_conjecture.f90
!
! 27^5 + 84^5 + 110^5 + 133^5 = 144^5
!
use, intrinsic :: iso_fortran_env, only: ik => int64
implicit none
integer(ik), parameter :: nmax = 6208_ik ! largest value that avoids overflow.
integer(ik) :: i, j, k, l, m, sk, sl, sm
integer(ik) :: n5(nmax)
integer :: T
character(*), parameter :: cfmt = '(*(g0,1x))'
outer: do i = 1_ik, nmax
n5(i) = i**5
T = mod( n5(i), 11_ik ) ! target modulus; T = 0, 1, or 10.
do j = 1_ik, i
do k = 1_ik, j
sk = n5(j) + n5(k)
if ( sk > n5(i) ) exit
if ( exclude_k(mod(sk,11_ik)) ) cycle
do l = 1_ik, k
sl = sk + n5(l)
if ( sl > n5(i) ) exit
if ( exclude_l(mod(sl,11_ik)) ) cycle
do m = 1_ik, l
sm = sl + n5(m)
if ( sm < n5(i) ) then
cycle
elseif ( sm > n5(i) ) then
exit
else ! sm == n5(i)
print cfmt, "i^5 = ", n5(i)
print cfmt, j, k, l, m, i
exit outer
endif
enddo
enddo
enddo
enddo
enddo outer
contains
pure logical function exclude_k( m11 )
! determine if the index k should be excluded.
integer(ik), intent(in) :: m11
select case (T)
case (0); exclude_k = 3 <= m11 .and. m11 <= 8
case (1); exclude_k = 3 <= m11 .and. m11 <= 9
case (10); exclude_k = 2 <= m11 .and. m11 <= 8
case default
error stop 'illegal T in exclude_k()'
end select
return
end function exclude_k
pure logical function exclude_l( m11 )
! determine if the index l should be excluded.
integer(ik), intent(in) :: m11
select case (T)
case (0); exclude_l = 2 <= m11 .and. m11 <= 9
case (1); exclude_l = 3 <= m11 .and. m11 <= 10
case (10); exclude_l = 1 <= m11 .and. m11 <= 8
case default
error stop 'illegal T in exclude_l()'
end select
return
end function exclude_l
end program sum_conjecture
$ gfortran -O3 sum_conjecture.f90 && time a.out
i^5 = 61917364224
133 110 84 27 144
real 0m0.072s
user 0m0.070s
sys 0m0.001s
I tried to clean up the code a little before posting, hopefully I didn’t introduce any new errors in that process. These times are for an Apple M2 cpu with MacOS, so they are a little faster than the previous M1 times. Obviously, these do index tests make the code a little more complicated than with the simple loops. Also, it might be good if someone else can verify that the two exclude_*() functions are correct. If they are too restrictive, then the timings will be too small, even if the single correct solution is still found. And if they not restrictive enough, then the timings could be further reduced legitimately.
I can see how this might be useful, like a Sieve of Eratosthenes type algorithm to eliminate possibilities. Then whatever is left over is a solution.
I had also tried this, but like you I found that it was not effective. I think the problem is that it only can be applied to the innermost loop, which has the smallest range of the nested loops. So the floating point operations to compute the starting point end up costing more than the integer additions that it replaces. However, a user function to compute the integer fifth root might still work here. It does not need to be accurate, but it needs to err on the small side rather than the large side of the exact starting value.
The following change does simplify the inner loop with -O3
do n = 1, l
if ( n5(n) < n5(i) - sl ) cycle
if ( n5(n) == n5(i) - sl ) then
Here is a Fortran-66 version for you to use in your CDC emulator or DTCyber. You may need to remove the “*8” in the variable declarations line.
C*NOLIST
INTEGER*8 I,P5(145)
INTEGER J,K,L,M
DO 10 I=1,145
10 P5(I)=I*I*I*I*I
DO 100 I = 1, 145
DO 100 J = 1, I
DO 100 K = 1, J
DO 100 L = 1, K
DO 60 M = 1, L
IF(P5(J) + P5(K) + P5(L) + P5(M)-P5(I))60,20,60
20 PRINT 777, J,K,L,M,I
GOTO 200
60 CONTINUE
70 CONTINUE
100 CONTINUE
200 STOP
777 FORMAT(27HEULER CONJECTURE DISPROVED.,5I10)
END
1 Like
There are a number of known results of the same type. Lander and his colleagues also found:
See this Wolfram page.
2 Likes
Nice, good luck trying to fit 85359**5 into a 64 bit integer:
>>> 85359**5
4531548087264753520490799
>>> 2**64
18446744073709551616
So that means they had to write arbitrary-precision integer implementation in 1966 in Fortran. Unless there are some tricks that one can do to avoid large integers.
1 Like
The CDC 6600 was a 60 bit (1-s complement) computer, 4 bits less than is common in today’s 64 bit machines. All the more reason for my wanting to find out what software was used in the calculations in 1966-1967. About 100 CDC 6600s were built. There was one at the University of Minnesota, where they also had a CDC 6400.
Ones-complement (which meant there was a negative zero as well as a positive zero, and they were sometimes not equal to each other) is probably unfamiliar to most readers today.
“There was a strange instruction on Cybers called the Population Count instruction. It returned the number of ‘1’ bits in the argument. It was rumored that this instruction was added to Cybers by the designer, Seymour Cray, at the request of the National Security Agency (NSA)”. Today’s X64 CPUs have a POPCNT instruction!
1 Like
According to the tables in Fredric Stuart’s book, “FORTRAN PROGRAMMING” from 1970 (my first Fortran programming book), there were two compilers used on the CDC 6600 at that time. These were called “Chippewa Fortran” and “Extended Fortran”. The integer range for both compilers was +/- 5.76460752303423487e17 (or 18 digits) which I think is 2**59-1.
I have vague memories of reading that the number of effective bits fell to 48 if integers were multiplied, since even integer multiplication was done using floating point hardware. In the present application, where integers have to be raised to the fifth power, this may imply that without additional MP support in software, the calculations require a range that was not available on the CDC.
In that time period, CDC offered the RUN compiler - which was probably what the benchmark was run with. Then they also had FTN (aka Fortran Extended 4), a newer compiler that did a lot more global optimization. At the end of the 1970s, there was FTN5 - which was a full Fortran 77 compiler. Steve Jasik re-implemented a lot of the optimization in it. FTN5 was probably the second Fortran 77 compiler to hit the streets - following the AT&T f77 compiler on unix.
There were a few Fortran compilers developed by universities. Notably MNF and later M77 from the U of Minnesota. There was also RUNT, a variant of RUN, from the U of Washington.
I can claim to have used all of the above…
Yes, the absolute value of operands and results of integer multiplication and division had to be less than 2^48 because of the use of the floating point multiply functional unit. Addition using the Long Add functional unit could use the full 2^59. (Later Cray Research systems had similar restrictions. However they had a compiler option to generate a longer sequence of instructions to compute using the full 64 bits.)
Looking at this code, the multiplications are only used when initially setting up table P5. So they don’t really matter - timing wise.
The inner loop is a simple search loop - sequentially loading values from P5 and comparing them to a temporary consisting of P5(J) + P5(K) + P5(L). So just a load, subtraction (or maybe even a boolean comparison) and possibly branch out of the loop. Then bump the index and do a backwards branch to do the next iteration. This loop would easily fit in the 6600’s 8-word instruction stack (basically a small i-cache).
While I’m not sure RUN could do it, FTN with optimization turned on would probably unroll the loop at least once. And no doubt do an optimization called “bottom loading” - where it would start the loads of a couple of the operands for the next iteration of the loop before branching backwards to hide the load times under the branch. Branching on the 6600 was relatively slow. So these optimizations could result in nice improvements.
I don’t have time to get dtcyber running today. Maybe tomorrow or sometime next week. I’m kind of interested in what the code that FTN and FTN5 will generate.
2 Likes
Stuart’s table says that mixed mode expressions were supported. If so I would think that even if the floating point hardware was used the calculations would be done to the required precision. One problem with the floating point calculations was only a 15 digit mantissa (significand) was supported.
I worked some more on this little program, and one thing I did was to change the declarations to
integer, parameter :: nmax = 6208
integer(ik) :: sk, sl, sm
integer(ik) :: n5(nmax)
integer :: i, j, k, l, m, T
This is because it is only the partial sums and the i**5 computations that need the extra precision, all the other integers need only range up to that 6208 value. This reduces the computational effort and the memory bandwidth a little, and that change reduces the Apple M2 run times down to
i^5 = 61917364224
133 110 84 27 144
real 0m0.059s
user 0m0.057s
sys 0m0.002s
which is something like a 15% speedup.
Then it occurred to me that there are many cases where the do m = 1, l loop cannot possibly find a solution because l is too small. There are two ways to avoid the innermost loop in that case, one is with a straightforward test beforehand, and the other is to switch the loop order to do m = l, 1, -1, and then change the inner tests accordingly. This allows early exit when the value is too small, which occurs on the first m loop cycle in the best case. That change then gets the timings down to
i^5 = 61917364224
133 110 84 27 144
real 0m0.031s
user 0m0.029s
sys 0m0.002s
So that looks like good progress.
There are still cases where l is much too large, so I thought back to @JohnCampbell’s comment about using an approximate fifth root to find a better starting point for the loop. Now that the loop is reversed, we need an upper bound rather than a lower bound. The floating point expression using **0.2 is too expensive. This is because the inner loop now only contains an addition, a couple of array lookups, and the comparisons. The comparisons are already ordered correctly, from most frequent to least frequent, so there’s nothing left to squeeze out there. But that means that the upper bound calculation needs to involve only a few integer operations, with no floating point operations or other expensive steps. So I tried the following approach. The expression
nbits = bit_size(1_ik) - leadz(n5(i)-sl)
gives the number of bits in the target value m**5. The expression 2**((nbits+4)/5)*5 rounds that value up to a number that is a fifth power of 2. Taking the fifth root then gives 2**((nbits+4)/5) as a possible starting point for the reversed loop. That number might be larger than l, in which case you might as well use l, but sometimes it is smaller, so it can save some loop iteratiions. The only problem left is that the exponention itself might be done in an expensive way, with recursive multiplications. But since it is a power of 2, it can also be computed with the ibset() intrinsic. So all of that gives finally the inner loop expression
do m = min(l,ibset(0,((nbits+4)/5))), 1, -1
Unfortunately, this change has no overall effect on the timings. I assume it is the leadz() step and/or the integer division, /5, makes even this integer expression too expensive. I left in this discussion here in case someone else wants to try it with their hardware+compiler to see if it makes any improvement in their environment.
So ignoring that last part about integer fifth roots, here is the best code I have so far.
program sum_conjecture
! sum_conjecture.f90
!
! 27^5 + 84^5 + 110^5 + 133^5 = 144^5
!
use, intrinsic :: iso_fortran_env, only: ik => int64
implicit none
integer, parameter :: nmax = 6208
integer(ik) :: sk, sl, sm
integer(ik) :: n5(nmax)
integer :: i, j, k, l, m, T, nbits
character(*), parameter :: cfmt = '(*(g0,1x))'
outer: do i = 1, nmax
n5(i) = int(i,ik)**5
T = mod(n5(i),11_ik) ! target modulus; T = 0, 1, or 10.
do j = 1, i
do k = 1, j
sk = n5(j) + n5(k)
if ( sk > n5(i) ) exit
if ( exclude_k(mod(sk,11_ik)) ) cycle
do l = 1, k
sl = sk + n5(l)
if ( sl > n5(i) ) exit
if ( exclude_l(mod(sl,11_ik)) ) cycle
!nbits = bit_size(1_ik) - leadz(n5(i)-sl)
!do m = min(l,ibset(0,((nbits+4)/5))), 1, -1
do m = l, 1, -1
sm = sl + n5(m)
if ( sm > n5(i) ) then
cycle
elseif ( sm < n5(i) ) then
exit
else ! sm == n5(i)
print cfmt, "i^5 = ", n5(i)
print cfmt, j, k, l, m, i
exit outer
endif
enddo
enddo
enddo
enddo
enddo outer
contains
pure logical function exclude_k( m11 )
! determine if the index k should be excluded.
integer(ik), intent(in) :: m11
select case (T)
case (0); exclude_k = 3 <= m11 .and. m11 <= 8
case (1); exclude_k = 3 <= m11 .and. m11 <= 9
case (10); exclude_k = 2 <= m11 .and. m11 <= 8
case default
error stop 'illegal T in exclude_k()'
end select
return
end function exclude_k
pure logical function exclude_l( m11 )
! determine if the index l should be excluded.
integer(ik), intent(in) :: m11
select case (T)
case (0); exclude_l = 2 <= m11 .and. m11 <= 9
case (1); exclude_l = 3 <= m11 .and. m11 <= 10
case (10); exclude_l = 1 <= m11 .and. m11 <= 8
case default
error stop 'illegal T in exclude_l()'
end select
return
end function exclude_l
end program sum_conjecture
2 Likes
I experimented a little more with the code. Here is a version where the inner do loop search over m is replaced with a binary search. I’m a little surprised that this one is competitive, but it has good timings, so here it is. I suspect that for some hardware+compiler combinations the loop search version I gave previously will be best, but for my Apple M2 machine, this one is the best so far, about a 25% improvement over the loop search version.
program sum_conjecture
! sum_conjecture.f90
!
! 27^5 + 84^5 + 110^5 + 133^5 = 144^5
!
use, intrinsic :: iso_fortran_env, only: ik => int64
implicit none
integer, parameter :: nmax = 6208
integer(ik) :: i5, sk, sl, sm, diff
integer(ik) :: n5(nmax)
integer :: i, j, k, l, m, T, nbits, mmin, mmax
character(*), parameter :: cfmt = '(*(g0,1x))'
intrinsic :: int, mod, bit_size, leadz, ibset
outer: do i = 1, nmax
i5 = int(i,ik)**5
n5(i) = i5
T = mod(n5(i),11_ik) ! target modulus; T = 0, 1, or 10.
do j = 1, i
do k = 1, j
sk = n5(j) + n5(k)
if ( sk > i5 ) exit
if ( exclude_k(mod(sk,11_ik)) ) cycle
do l = 1, k
sl = sk + n5(l)
diff = i5 - sl ! target m**5 value.
if ( diff < 0 ) exit
if ( exclude_l(mod(sl,11_ik)) ) cycle
if ( n5(l) < diff ) cycle ! no possible solution in 1...l.
! set up for binary search for m**5 = diff.
nbits = bit_size(diff) - leadz(diff)
mmin = ibset(0,((nbits-1)/5)) ! lower bound to m.
mmax = min( l, ibset(0,((nbits+4)/5))) ! upper bound to m.
if ( n5(mmin) == diff ) then ! check for shortcuts.
mmax = mmin
elseif ( n5(mmax) == diff ) then
mmin = mmax
endif
do
m = (mmin + mmax) / 2
if ( n5(m) > diff ) then
mmax = m
elseif ( n5(m) < diff ) then
if ( m == mmin ) exit ! no solution.
mmin = m
else ! found.
print cfmt, "i^5 = ", i5
print cfmt, j, k, l, m, i
exit outer
endif
enddo
enddo
enddo
enddo
enddo outer
contains
pure logical function exclude_k( m11 )
! determine if the index k should be excluded.
integer(ik), intent(in) :: m11
select case (T)
case (0); exclude_k = 3 <= m11 .and. m11 <= 8
case (1); exclude_k = 3 <= m11 .and. m11 <= 9
case (10); exclude_k = 2 <= m11 .and. m11 <= 8
case default
error stop 'illegal T in exclude_k()'
end select
return
end function exclude_k
pure logical function exclude_l( m11 )
! determine if the index l should be excluded.
integer(ik), intent(in) :: m11
select case (T)
case (0); exclude_l = 2 <= m11 .and. m11 <= 9
case (1); exclude_l = 3 <= m11 .and. m11 <= 10
case (10); exclude_l = 1 <= m11 .and. m11 <= 8
case default
error stop 'illegal T in exclude_l()'
end select
return
end function exclude_l
end program sum_conjecture
gfortran -O3 sum_conjecture.f90 && time a.out
i^5 = 61917364224
133 110 84 27 144
real 0m0.026s
user 0m0.024s
sys 0m0.002s
I also replaced the m=(mmin+mmax)/2 statement with m=ishft(mmin+mmax,-1), but the timings were the same. That is valid if you know that the argument is nonnegative, so it is not something that the compiler will typically do.
Instead of a straight binary search, I also experimented a little with a Newton search. However, that requires some integer multiplications and divides, so it costs a little too much to be competitive with timings. However, there might be some future to that approach too. I have some code somewhere where I combined a Newton search with a binary search, so I might dig that out and try that if I see a way to avoid he integer divisions.
edit: I wanted to try a 128-bit version with flang, but it does not support the i**5 calculation or the ibset() and leadz() intrinsics. So I made the following changes to the code:
use, intrinsic :: iso_fortran_env, only: ik => int128
...
i5 = i * i ! i**2
i5 = i5 * i5 ! i**4
i5 = i * i5 ! i**5
...
mmin = 1
mmax = l
...
With those changes, I get:
$ flang -O3 sum_conjecture.f90 && time a.out
i^5 = 61917364224
133 110 84 27 144
real 0m0.079s
user 0m0.076s
sys 0m0.003s
With int128, huge=170141183460469231731687303715884105727 which would allow in principle searches up to nmax=44275350. Even with what is now close to an N^4 algorithm, that would still take a long time to find that other (85282, 28969, 3183, 55, 85359) solution with this algorithm.
I should have mention this before, but setting mmin=1 and mmax=l gives very similar timings to the more restricted search domain. I think it is again a tradeoff between paying for the bounds calculations and just accepting the larger domain and letting the binary search go longer. If each step of the search required more effort, then the narrow domain setup would be more desirable.
1 Like
As an undergraduate engineer, I got to use a CDC 6600 at University of Sydney in 1974 when they built a new computer centre. (In 1973 we only had access to IBM 7040 or an ? Algol machine in the Physics dept.)
In about 1976, it was later upgraded to a CDC 7600, or we were allowed to use the 7600. Not sure witch. I think the Fortran version changed between the machines.
However, most of us moved to a Pr1me 300, which didn’t have a cost code and did have a better Fortran. With a default integer*2, not a lot of interest in the Euler Conjecture there !
I remember a frightening experience when, in the 1960s, I put my card deck in the input queue and went back a few hours later to pick up the deck and the printer output as usual, but found only a slip of paper with “your program executed data” scribbled on it. I did not understand that such an act was possible and how I had committed it!
1 Like