Prime Numbers

JohnNichols · ‎08-05-2020

Dear Reader:

mecej4 very kindly turned a brute force and very inefficient code using OMP into a fast inefficient code.

I stumbled on this number theory problem whilst I was really bored with a awful proposal, still am.

It is really interesting -- and this is just a simple example

Late last night I started playing with the underlying numbers, as long as we remain in the set of integers we are stuck with a large set of inefficient numbers, but if we reduce to the set of primes it is a much smaller set.

each of the numbers in solution a^5+b^5 etc can be reduced to (the rhs can be represented as a five number - )

a = a1*a2*a3*a4 - where the ais are all small primes less than 19 for the first thousand integers.

It is easy to then create a table of these a and b numbers to the 5th power as a look up table -- by 100000 there density is 7% of the integers -- after that it is only lookup and addition -- to check the rhs -- turn the added result into a real -- take the fifth root and see if there is a ceiling remainder.

I had a look at the Lander 1967 method -- whilst it was quicker I thought there was a quicker way - it may not be

if the ceiling remainder is close to zero do an integer check --

So late last night I coded it -- runs in about 5 minutes, but it has a Fortran bug with the powers function -- I am making a mistake -- I will look at today

I added a new post as the old one was getting long

That is why the great mathematicians say -- a piece of paper is your best friend after all you lot of course

Here is the code -- it does not yest work -- I am slowly debugging it and then I will use mecej4 example and parallelize it -- I borrowed the prime code -- if it works it should be fast -- if it does not work then I will merely look like a failed mathematician as my year 10 teacher once said -- give up you and math do not fit.

PS I really appreciate the help

Here if you do not treat one as a prime you need to treat the a and bs as variable length and it is easy to just add one as prime -- here I am somewhat baffled by the pure math jocks who treat one as not prime, it has all of the characteristics of a prime for this type of problem and a lot of others -- but it is still prime like. If my grandmother can feed me - do I turn her away at the door because she is carrying food and more mother than the average mother --

Really the issue is the weird numbers 1, -1, i, e and pi. God in her infinite wisdom sent them as a challenge to mere mortals.

! ---------------------------------------------------------------
! This program finds all prime numbers in the range of 2 and an
! input integer.
! ---------------------------------------------------------------

PROGRAM  Primes
   IMPLICIT  NONE
   integer*8 a(100), b(7500), c(7500),f(4)

   INTEGER*8  :: Range, Number, Divisor, Count, n, i,j,k,l,m
   real*8  d,e

   a(1) = 1
    b = 0
    c = 0
   WRITE(*,1000) 
1000 format(   'What is the range ? ',\)
   DO                                   ! keep trying to read a good input
      READ(*,*)  Range                  ! ask for an input integer
      IF (Range >= 2)  EXIT             ! if it is GOOD, exit
      WRITE(*,*)  'The range value must be >= 2.  Your input = ', Range
      WRITE(*,*)  'Please try again:'   ! otherwise, bug the user
   END DO

   Count = 1                            ! input is correct. start counting
   WRITE(*,*)                           ! since 2 is a prime
   WRITE(*,*)  'Prime number #', Count, ': ', 2
   a(count+1) = 2 
   DO Number = 3, Range, 2              ! try all odd numbers 3, 5, 7, ...

      Divisor = 3                       ! divisor starts with 3
      DO
         IF (Divisor*Divisor > Number .OR. MOD(Number,Divisor) == 0)  EXIT
         Divisor = Divisor + 2          ! if does not evenly divide, next odd
      END DO

      IF (Divisor*Divisor > Number) THEN     ! are all divisor exhausted?
         Count = Count + 1              ! yes, this Number is a prime
         WRITE(*,*)  'Prime number #', Count, ': ', Number
         a(count+1) = number
      END IF
   END DO

   WRITE(*,*)
   WRITE(*,*)  'There are ', Count, ' primes in the range of 2 and ', Range
   
   n = count + 1
   
   do 100 i = 1,n
       do 200 j = 1,n
           do 300 k = 1,n
               do 400 l = 1,n
                count = a(i)*a(j)*a(k)*a(l)
                if(count .lt. 7501) then 
                b(count) = 1
                c(count) = count*count*count*count*count
                write(*,*)c(count)
                end if
400            end do
300        end do
200    end do
100 end do
    do 1100 i = 1,n
       do 1200 j = 1,n
           do 1300 k = 1,n
               do 1400 l = 1,n
                  count = c(i)+c(j)+c(k)+c(l)
                  f(1) = (27**5)  
                  f(2) =(84**5) !+(110**5)+(133**5)
                  d = real(count,8)
                  d = d**0.2d0
                  e = floor(d)
                  if((d - e) .lt. 0.00000001d0)then
                  write(*,*) d
                  end if
1400            end do
1300        end do
1200    end do
1100 end do
    
       write(*,*)'here'            

END PROGRAM  Primes

The final step is to do very large integers with Fortran code --

JohnNichols · ‎08-08-2020

It is much faster to create a list of possible numbers that can be a valid entry into Euler Conjecture equation using small primes, there are only 359 that will fit into the size of integer available on the Intel Fortran Compiler once raised to the 5th power, using integers to 19 and then not all combinations.

The next step is to loop through the list -- although we will find multiples so again it will be a reduced search -- I was looking at this stage at the brute force just to get the algorithm sorted out. The method I am using is to create the number for the LHS and take the 1/5th power, take the ceiling and then raise the number to the 5th power and compare the two integers. There are look up tables and the mecej4 assistance to speed it up more, plus the RHS has limited entries.

I am going to be a heretic -- the Euler Conjecture is a statistical statement -- the lhs of the equation provides an integer say T, we take the fifth root and we are left with a real number S, if the S is functionally equal to an integer we have a answer - first one is 144. But if we loop over all of the valid T's we get a set of S and these using the ceiling function gives us an array from 0 to 1. This will have a statistical distribution shape, the Conjecture is merely that we find the element with an exact value of 0. We know it is rare, but really what it is saying for the four integers raised to the fifth power is what it is the interesting probability that five integer set raised to the fifth power have an element that is zero.

The problem is the definition is positive integers only for the elements. If we include 0 it simply means is there a SHAPE with five sides where one side is zero, this is merely, for example, the fact that a triangle is four sided shape with a zero length.

Anyway the real problem is I need very large integers, playing up to 19 is a waste, and as the primes get more sparse it should easier but we need big integers in FORTRAN -- we can play with the RHS to limit its look up table and probably speed things up.

So consider that in the first

61917364224

numbers there is only one that fits into the equation that is a probability the sun will explode tomorrow level, not quite but you get the picture.

Here is the latest code -- it is far from perfect, but I do not want to have more fun until I can find a large Fortran int package.

    ! ---------------------------------------------------------------
    ! This program finds all prime numbers in the range of 2 and an
    ! input integer.
    ! ---------------------------------------------------------------

    PROGRAM  Primes

    IMPLICIT  NONE
    integer*8, parameter :: nn = 7500
    integer*8 a(100), b(7500), c(nn),f(4)

    INTEGER*8  :: Range, Number, Divisor, Count, n, i,j,k,l,m,counter
    real*8  d,e
    integer*2 jk

    a(1) = 1
    b = 0
    c = 0
    counter = 0
    jk = -1
    WRITE(*,1000)
1000 format(   'What is the range ? ',\)
    DO                                   ! keep trying to read a good input
        READ(*,*)  Range                  ! ask for an input integer
        IF (Range >= 2)  EXIT             ! if it is GOOD, exit
        WRITE(*,*)  'The range value must be >= 2.  Your input = ', Range
        WRITE(*,*)  'Please try again:'   ! otherwise, bug the user
    END DO

    Count = 1                            ! input is correct. start counting
    WRITE(*,*)                           ! since 2 is a prime
    WRITE(*,*)  'Prime number #', Count, ':           ', 2
    a(count+1) = 2
    DO Number = 3, Range, 2              ! try all odd numbers 3, 5, 7, ...

        Divisor = 3                       ! divisor starts with 3
        DO
            IF (Divisor*Divisor > Number .OR. MOD(Number,Divisor) == 0)  EXIT
            Divisor = Divisor + 2          ! if does not evenly divide, next odd and next
        END DO

        IF (Divisor*Divisor > Number) THEN     ! are all divisor exhausted?
            Count = Count + 1              ! yes, this Number is a prime
            WRITE(*,*)  'Prime number #', Count, ': ', Number
            a(count+1) = number
        END IF
    END DO

    WRITE(*,*)
    WRITE(*,50)Count,Range
50  Format('There are ', i4, ' primes in the range of 2 and ', i4, ' excluding one.')

    n = count + 1

    do 100 i = 1,n-3
        do 200 j = 1,n-3
            do 300 k = 1,n
                do 400 l = 1,n
                    count = a(i)*a(j)*a(k)*a(l)
                    if(count .lt. 6207) then
                        counter = counter + 1
                        b(count) = 1
                        c(count) = count*count*count*count*count
                        write(*,500)counter,a(i),a(j),a(k),a(l),c(count)
500                     Format(' Count :: ',i6,' : ',4(2x,i4),'   ',I23)
                    end if
400 end do
300 end do
200     end do
        write(*,250)huge(c(count))
250     Format(49x,i20)
100 end do

    call quicksortA(c,nn, counter)



    write(*,*)'here'

    END PROGRAM  Primes

    

    subroutine quicksortA(a,n, counter)
    integer*8 n, dA
    integer*8 a(n), counter
    integer*8 b(counter),dB,dC
    integer count, i,j,k,l
    real*8 d,e

    j = 0
    b = 0
    count = 0

    do 100 i = 1,n

        if(a(i) .gt. 0) then
            count = count + 1
            if(j .gt. 359) then
                write(*,*)count
            endif
            b(count) = a(i)

           ! write(*,*)i,a(i),count,b(count)
        endif
100 end do

    do 1100 i = 1,count
        do 1200 j = 1, count
            do 1300 k = 1, count
                do 1400 l = 1, count
                    dA = b(i) + b(j) + b(k) + b(l)
                    d = real(dA,8)
                    d = d**0.2d0
                    e = floor(d)
                    write(*,*)da,d,e
                    if((d - e) .lt. 0.00000000001d0)then
                        dB = int(d,8)
                        dB = dB*dB*dB*dB*dB
                        dC = dA - dB
                        
                        write(*,*) d,dC
                    end if
1400            end do
1300        end do
1200    end do
1100 end do

    return
    end subroutine quicksortA

GVautier · ‎08-09-2020

As I said before, years ago, I developed arithmetic functions operating in base 256 able to deal with unlimited integers. If you are interested, I can take a look and refresh this old code.

JohnNichols · ‎08-09-2020

I would really appreciate that please.

GVautier · ‎08-09-2020

OK, I will exhume that F77 code and take a look at it next week.

JohnNichols · ‎08-09-2020

I really appreciate that.

JohnNichols · ‎08-09-2020

This is the large number package from Smith at Loyola transferred to Intel Fortran - minor name changes in files and extensions

It ran first time -- took about 5 minutes to compile

How do I turn it into a library so it can be easily used?

25 years since I built a library

JohnNichols · ‎08-09-2020

1. Read the manual

2. Lib attached in zip file