program euler
! Finds all the integer solutions of the Diophantine equation
!
!    p^5 + q^5 + r^5 + s^5 = n^5, p,q,r,s and n integers > 0
!
! for all values of n from 1 to 999
!
! Only one primitive solution has been found:
!
!   L.J. Lander, T.R. Parkin and J.L. Selfridge:
!      A Survey of Equal Sums of Like Powers
!           Mathematics of Computation
!     Vol. 21, No. 99 (Jul., 1967), pp. 446-459
!
implicit none
integer*8 s
integer nv, vals(4),g,gcf
character*10 tm
integer :: iu(7) = [999,955,900,840,760,640,1]  ! for 6 threads, assuming O(n^2)
integer i, j, nt, omp_get_thread_num, ithrd
logical found
!$OMP PARALLEL private (i,j,s,nv,vals,found,g)
!$OMP DO
do nt=1,6
do i=iu(nt),iu(nt+1)+1,-1
   ithrd = omp_get_thread_num()
   if(mod(i,10)==0)then
      call date_and_time(time=tm)
      print '(1x,i2,2x,A,i3,2x,3(A2,:,":"))',ithrd, &
         'Trying x = ',i,tm(1:2),tm(3:4),tm(5:6)
   end if
   s=i
   s=s*s*s*s*s
   nv=0
   call search(s,4,found,i-1,nv,vals)
   if(found)then
      g = vals(1)
      do j=2,4
        g = gcf(g,vals(j))
      end do
      if(g == 1)then
         print 10,ithrd,vals,i
      else
         print 10,ithrd,vals,i,g
      endif
   endif
10 format(1x,i2,3(2x,I3,'^5 + '),2x,I3,'^5 = ',I3,'^5',:,&
          2x,'NP, GCF = ',i2)
end do
end do
!$OMP END DO
!$OMP END PARALLEL
end program

recursive subroutine search(s,n,found,ilim,nv,vals)
implicit none
integer*8 s,j8,i8
integer n,j,im,ilim,nv,vals(4)
logical found

found=.false.
!search for n integers whose 5th power equals s
if(n == 1)then
   call search1(s,found,im,nv,vals)
   if(found)then
      nv=nv+1
      vals(nv)=im
   endif
   return
endif
do j8=ilim,1,-1
   i8=j8*j8*j8*j8*j8
   if(i8 > s)cycle
   j=j8
   call search(s-i8,n-1,found,j-1,nv,vals)
   if(found)then
      nv=nv+1
      vals(nv)=j8
      exit
   endif
end do

return
end subroutine

subroutine search1(s,found,im,nv,vals)
implicit none
integer*8 s
integer il,ir,im,i,nv,vals(4)
integer*8 i5,im5,vl,vr
logical found,start
data start/.true./
integer*8, save :: pow5tbl(1024)
if(start)then
   do i=1,1024
      i5=i
      pow5tbl(i)=i5*i5*i5*i5*i5
   end do
   start=.false.
endif
! Brute Force: il=1; ir=1024
il=nint(s**0.2); ir=il+1; il=il-1
vl = pow5tbl(il); vr = pow5tbl(ir)
if((vl < s .and. vr < s) .or. (vl > s .and. vr > s))then
   found=.false.
   return
endif
do while(il /= ir)
  im = (il+ir)/2
  im5=pow5tbl(im)
  if(s > im5)then
     il=im
  else if(s < im5)then
     ir=im
  else
     found=.true.
     return
  endif
  if(ir <= il+1)exit
end do
found = pow5tbl(im) == s
return
end subroutine

integer function gcf(l,m)
! GCF by Euclid's algorithm
integer l,m,i,j,t
i=l; j=m
if(i > j)then
   t = j
   j = i
   i = t
endif
do
   t = mod(j,i)
   if(t == 0)exit
   j = i
   i = t
end do
gcf =i
return
end function