- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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 --
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
OK, I will exhume that F77 code and take a look at it next week.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
1. Read the manual
2. Lib attached in zip file
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page