Home > Archive > Fortran > August 2005 > little isprime challenge
You are viewing an archived Text-only version of the thread.
To view this thread in it's original format and/or if you want to reply to
this thread please [click here]
| Author |
little isprime challenge
|
|
| Bart Vandewoestyne 2005-08-23, 7:00 pm |
| I know there has been a kind of `string-challenge' in this group
which I have not been following... but it's been taking a while
and from what I've read I have the feeling that some people got
bored of it and now are definitely in need for something else...
so I have another little `isprime-challlenge' :-)
It's just 100% pure for the fun of it, and maybe it will also
teach some more novices and `less experienced' young fortran
programmers as me good optimization skills. It might also be usefull for
some people in need for some prime routines.
The `assignment' would be (I will try to formulate it as formal
as possible in order to avoid confusion...):
"Write a standard Fortran 95 function isprime(p) that also adheres to
the F-language and that checks if the number p is a prime. The routine
should be as fast as possible (to be measured with the CPU_TIME intrinsic.
2 is considered to be a prime. The speed of your routine will be measured
by checking the numbers 1,n where n is a very big number, whatever
`very big' may be ;-)"
I've put up a skeleton file that is available online at
http://www.cs.kuleuven.ac.be/~bartv...ime_isprime.f95
This is the file where you should add your version to.
My first result on an Intel(R) Pentium(R) 4 CPU 2.40GHz running
a stable Debian GNU/Linux and compiled with the latest
F-compiler:
bartv@vonneumann:~/fortran$ make time_isprime
F -o time_isprime -O4 time_isprime.f95
bartv@vonneumann:~/fortran$ ./time_isprime
Enter the number of points you want to time for: 10000000
Bart's method found 664579 primes in 27.120001 seconds.
Good luck!
Bart
:-)
--
"Share what you know. Learn what you don't."
| |
| Steven G. Kargl 2005-08-23, 7:00 pm |
| In article <1124835304.440418@seven.kulnet.kuleuven.ac.be>,
Bart Vandewoestyne <MyFirstName.MyLastName@telenet.be> writes:
>
> 2 is considered to be a prime.
>
This is an odd statement. Has 2 ever been considered
a non-prime number?
--
Steve
http://troutmask.apl.washington.edu/~kargl/
| |
| e p chandler 2005-08-23, 9:57 pm |
|
Bart Vandewoestyne wrote:
> I know there has been a kind of `string-challenge' in this group
> which I have not been following... but it's been taking a while
> and from what I've read I have the feeling that some people got
> bored of it and now are definitely in need for something else...
> so I have another little `isprime-challlenge' :-)
>
> It's just 100% pure for the fun of it, and maybe it will also
> teach some more novices and `less experienced' young fortran
> programmers as me good optimization skills. It might also be usefull for
> some people in need for some prime routines.
>
> The `assignment' would be (I will try to formulate it as formal
> as possible in order to avoid confusion...):
>
> "Write a standard Fortran 95 function isprime(p) that also adheres to
> the F-language and that checks if the number p is a prime. The routine
> should be as fast as possible (to be measured with the CPU_TIME intrinsic.
> 2 is considered to be a prime. The speed of your routine will be measured
> by checking the numbers 1,n where n is a very big number, whatever
> `very big' may be ;-)"
>
> I've put up a skeleton file that is available online at
>
> http://www.cs.kuleuven.ac.be/~bartv...ime_isprime.f95
>
> This is the file where you should add your version to.
>
> My first result on an Intel(R) Pentium(R) 4 CPU 2.40GHz running
> a stable Debian GNU/Linux and compiled with the latest
> F-compiler:
>
> bartv@vonneumann:~/fortran$ make time_isprime
> F -o time_isprime -O4 time_isprime.f95
> bartv@vonneumann:~/fortran$ ./time_isprime
> Enter the number of points you want to time for: 10000000
> Bart's method found 664579 primes in 27.120001 seconds.
>
> Good luck!
> Bart
> :-)
>
> --
> "Share what you know. Learn what you don't."
I'm going to go for the "fairly easy" and post the result of a simple
modification to your function. Instead of testing the candidate against
all numbers from 2 and up, why not test for 2, then test for
divisibility by 2, then test against odd numbers from 3 and up?
The core of my function is:
res = .false.
if (p < 2) then
res = .false.
return
else if (p == 2) then
res = .true.
return
else if (modulo(p,2) == 0) then
res = .false.
return
else
do i=3,p,2
if (i*i > p) then
res = .true.
return
else if (modulo(p,i) == 0) then
res = .false.
return
end if
end do
end if
Here's a run on a 2.8 GHz P4-HT under WinXP(SP2):
C:\temp>F -O4 -o test.exe test.f95
C:\temp>test
Enter the number of points you want to time for: 10000000
Bart's method found 664579 primes in 23.156250 seconds.
Elliot's method found 664579 primes in 11.656250
seconds.
This took only a few minutes to write and test but it cut the run time
nearly in half.
If I was doing this program in ye olden days, I would have put in a
COMMON block which contains an array holding the primes found so far.
Then I would have divided by the list of primes already found. Of
course F does not allow COMMON blocks. :-).
| |
| James Van Buskirk 2005-08-23, 9:57 pm |
| "Bart Vandewoestyne" <MyFirstName.MyLastName@telenet.be> wrote in message
news:1124835304.440418@seven.kulnet.kuleuven.ac.be...
> bartv@vonneumann:~/fortran$ make time_isprime
> F -o time_isprime -O4 time_isprime.f95
> bartv@vonneumann:~/fortran$ ./time_isprime
> Enter the number of points you want to time for: 10000000
> Bart's method found 664579 primes in 27.120001 seconds.
Well, your testing methodology leads one to write a Sieve of
Eratosthenes rather that a primality tester. See
http://home.comcast.net/~kmbtib/Sieve30h/ for example code of
the former type.
And I was hoping to get a chance to reference my newly
uploaded http://home.comcast.net/~kmbtib/real_ft.ZIP
but your challenge doesn't seem to be directly related...
--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end
| |
| Arjen Markus 2005-08-24, 3:57 am |
| I guess Bart is confusing 2 with the number 1. Given the definition, it
_is_ a prime :)
Regards,
Arjen
| |
| Bart Vandewoestyne 2005-08-24, 3:57 am |
| On 2005-08-24, e p chandler <epc8@juno.com> wrote:
>
> I'm going to go for the "fairly easy" and post the result of a simple
> modification to your function. Instead of testing the candidate against
> all numbers from 2 and up, why not test for 2, then test for
> divisibility by 2, then test against odd numbers from 3 and up?
>
> The core of my function is:
>
> <snip core>
>
> Here's a run on a 2.8 GHz P4-HT under WinXP(SP2):
>
> C:\temp>F -O4 -o test.exe test.f95
>
> C:\temp>test
> Enter the number of points you want to time for: 10000000
> Bart's method found 664579 primes in 23.156250 seconds.
> Elliot's method found 664579 primes in 11.656250
> seconds.
>
> This took only a few minutes to write and test but it cut the run time
> nearly in half.
Nice, very nice! I'm Looking forward to the next faster one :-)
Regards,
Bart
--
"Share what you know. Learn what you don't."
| |
| Bart Vandewoestyne 2005-08-24, 3:57 am |
| On 2005-08-24, Arjen Markus <arjen.markus@wldelft.nl> wrote:
> I guess Bart is confusing 2 with the number 1. Given the definition, it
> _is_ a prime :)
I guess i was indeed confusing... i changed it in 'The Mission'
statement on top of my testing program. I have also updated the
skeleton program a bit. The new version is where the previous
version was:
http://www.cs.kuleuven.ac.be/~bartv...ime_isprime.f95
Regards,
Bart
--
"Share what you know. Learn what you don't."
| |
| Bart Vandewoestyne 2005-08-24, 3:57 am |
| On 2005-08-24, James Van Buskirk <not_valid@comcast.net> wrote:
>
> Well, your testing methodology leads one to write a Sieve of
> Eratosthenes rather that a primality tester. See
> http://home.comcast.net/~kmbtib/Sieve30h/ for example code of
> the former type.
I have some other prime-related functions in mind that I might
also give as a challenge here, but for now, the procedure to get as
fast as possible (for the ones interested of course ;-) is the
function isprime(p). I thought the simplest way to test for it's
speed is to let it check a bunch of numbers from 1 to n. I am
always happy to hear of a better method of measuring it's speed.
> And I was hoping to get a chance to reference my newly
> uploaded http://home.comcast.net/~kmbtib/real_ft.ZIP
> but your challenge doesn't seem to be directly related...
I have other prime-number related challenges in mind (maybe one a
w or month or so ;-). But for now, i would first like to
stick to the isprime(p) function. I believe my prime-challenges
will be interesting for a lot of the people here because prime
stuff is something we all need once in a while, isn't it? :-)
Regards,
Bart
--
"Share what you know. Learn what you don't."
| |
| James Parsly 2005-08-24, 6:59 pm |
| Would it be faster to compute the square root of p, and use it as the loop
limit?
This would eliminate the need for the 'i*i' test.
If you want to keep the I*I test, there's a way to turn the multiplication
into
addition. Initial a new variable ISQUARE to 9 before the beginning of the
loop
and then update at the end with ISQUARE = ISQUARE + I + I + 1
or ISQUARE = ISQUARE + 2*I+1 (which is better depends on how
clever your compiler is, since multiplying by 2 can often be optimized into
a bit shift operation). Whether this will actually be faster is hard to
say, try it
and see.
"e p chandler" <epc8@juno.com> wrote in message
news:1124848022.896820.218980@f14g2000cwb.googlegroups.com...
>
> Bart Vandewoestyne wrote:
>
> I'm going to go for the "fairly easy" and post the result of a simple
> modification to your function. Instead of testing the candidate against
> all numbers from 2 and up, why not test for 2, then test for
> divisibility by 2, then test against odd numbers from 3 and up?
>
> The core of my function is:
>
> res = .false.
>
> if (p < 2) then
> res = .false.
> return
> else if (p == 2) then
> res = .true.
> return
> else if (modulo(p,2) == 0) then
> res = .false.
> return
> else
> do i=3,p,2
> if (i*i > p) then
> res = .true.
> return
> else if (modulo(p,i) == 0) then
> res = .false.
> return
> end if
> end do
> end if
>
> Here's a run on a 2.8 GHz P4-HT under WinXP(SP2):
>
> C:\temp>F -O4 -o test.exe test.f95
>
> C:\temp>test
> Enter the number of points you want to time for: 10000000
> Bart's method found 664579 primes in 23.156250 seconds.
> Elliot's method found 664579 primes in 11.656250
> seconds.
>
> This took only a few minutes to write and test but it cut the run time
> nearly in half.
>
> If I was doing this program in ye olden days, I would have put in a
> COMMON block which contains an array holding the primes found so far.
> Then I would have divided by the list of primes already found. Of
> course F does not allow COMMON blocks. :-).
>
| |
| Rich Townsend 2005-08-24, 6:59 pm |
| James Parsly wrote:
> Would it be faster to compute the square root of p, and use it as the loop
> limit?
> This would eliminate the need for the 'i*i' test.
Yep, this gives a modest speed up (modest, because SQRT() is slow, and
requires floating-point arithmetic). Here's the code:
function isprime_rich(p) result (res)
integer(kind=i4b), intent(in) :: p
integer(kind=i4b) :: i
logical :: res
res = .true.
if (p < 2) then
res = .false.
return
else if (p == 2) then
return
else if (modulo(p,2) == 0) then
res = .false.
return
else
do i=3,CEILING(SQRT(REAL(p))),2
if (modulo(p,i) == 0) then
res = .false.
return
end if
end do
end if
end function isprime_rich
Results (1.6GHz Pentium M Dothan):
Enter the number of points you want to time for: 10000000
Bart's method found 664579 primes in 29.443523 seconds.
Elliot's method found 664579 primes in 14.835743 seconds.
Rich's method found 664579 primes in 13.664925 seconds.
cheers,
Rich
| |
|
| And mod is faster than modulo:
Enter the number of points you want to time for: 3000000
Elliot's method found 216816 primes in 2.51 seconds.
Bart's method found 216816 primes in 4.99 seconds.
Rich's method found 216816 primes in 3.56 seconds.
Joost's method found 216816 primes in 1.92 seconds.
function isprime_joost(p) result (res)
integer(kind=i4b), intent(in) :: p
logical :: res
integer(kind=i4b) :: i,top
res = .true.
if (p < 2) then
res = .false.
return
else if (p == 2) then
res = .true.
return
else if (mod(p,2) == 0) then
res = .false.
return
else
do i=3,p,2
if (i*i>p) exit
if (mod(p,i) == 0) then
res = .false.
return
end if
end do
end if
end function isprime_joost
Joost
| |
| Rich Townsend 2005-08-24, 6:59 pm |
| Joost wrote:
> And mod is faster than modulo:
>
> Enter the number of points you want to time for: 3000000
> Elliot's method found 216816 primes in 2.51 seconds.
> Bart's method found 216816 primes in 4.99 seconds.
> Rich's method found 216816 primes in 3.56 seconds.
> Joost's method found 216816 primes in 1.92 seconds.
I'm not finding the same speedup from switching to MOD() on Intel ifort.
Also, I'm suprised my method (well, actually James Parsly's method,
implemented by me) using SQRT() is *slower* than Elliot's method, when I
found the reverse.
Which compiler are you using, on what platform?
cheers,
Rich
| |
|
| I was also surpised to see that you trick was slower, I don't really
understand why. That mod is faster than modulo I've observed before.
for these timing I have been using
Intel(R) Fortran Compiler for 32-bit applications, Version 9.0 Build
20050624Z Package ID: l_fc_c_9.0.024
on
vendor_id : GenuineIntel
cpu family : 15
model : 2
model name : Intel(R) Xeon(TM) CPU 3.06GHz
stepping : 9
I guess these micro-optimisations are CPU specific anyway, somebody
might have a clean explanation.
Cheers,
Joost
| |
| Bart Vandewoestyne 2005-08-24, 6:59 pm |
| On 2005-08-24, Joost <jv244@cam.ac.uk> wrote:
> And mod is faster than modulo:
>
> Joost
Sorry Joost, but MOD is not an F-intrinsic and the solution must be
F-compatible...
bartv@vonneumann:~/fortran$ grep mod\( example_modulo.f95
print *, "mod(7, 2) = ", mod(7, 2)
print *, "mod(13, 5) = ", mod(13, 5)
bartv@vonneumann:~/fortran$ make example_modulo
F -o example_modulo -O4 example_modulo.f95
Error: example_modulo.f95, line 28: MOD is not an F intrinsic procedure
detected at EXAMPLE_MODULO@<end-of-statement>
[F terminated - errors found by pass 1]
make: *** [example_modulo] Error 2
Regards,
Bart
--
"Share what you know. Learn what you don't."
| |
| Bart Vandewoestyne 2005-08-24, 6:59 pm |
| On 2005-08-24, Rich Townsend <rhdt@barVOIDtol.udel.edu> wrote:
>
> Also, I'm suprised my method (well, actually James Parsly's method,
> implemented by me) using SQRT() is *slower* than Elliot's method, when I
> found the reverse.
>
> Which compiler are you using, on what platform?
I also found your method to be slower than Elliot's:
bartv@vonneumann:~/fortran$ ./time_isprime
Enter the number of points you want to time for: 10000000
Elliot's method found 664579 primes in 13.76 seconds.
Bart's method found 664579 primes in 27.22 seconds.
Rich's method found 664579 primes in 14.90 seconds.
This was compiled with Fortran Company/NAG F compiler Release 20031017 using
F -o time_isprime -O4 time_isprime.f95
on an Intel(R) Pentium(R) 4 CPU 2.40GHz on Debian GNU/Linux (stable).
I have a feeling that if we're still looking for a *real* speedup, it will
have to be something algorithmic rather than these kind of micro-changes...
Greetzzz,
Bart
--
"Share what you know. Learn what you don't."
| |
| Greg Lindahl 2005-08-24, 6:59 pm |
| In article <1124905105.559996.155750@f14g2000cwb.googlegroups.com>,
Joost <jv244@cam.ac.uk> wrote:
>I guess these micro-optimisations are CPU specific anyway, somebody
>might have a clean explanation.
I'd guess that they're compiler-specific. A good compiler could turn
i*i into an induction variable, which is what the "+ 2*i + 1"
"optimization" is doing.
-- greg
| |
| Rich Townsend 2005-08-24, 6:59 pm |
| Greg Lindahl wrote:
> In article <1124905105.559996.155750@f14g2000cwb.googlegroups.com>,
> Joost <jv244@cam.ac.uk> wrote:
>
>
>
>
> I'd guess that they're compiler-specific. A good compiler could turn
> i*i into an induction variable, which is what the "+ 2*i + 1"
> "optimization" is doing.
Why does that help? You substitute one multiply by a multiply and two
adds, n'est ce pas?
cheers,
Rich
| |
| Brooks Moses 2005-08-24, 6:59 pm |
| NNTP-Posting-Host: dnab42a4a5.stanford.edu
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit
X-Trace: news.Stanford.EDU 1124912640 21733 171.66.164.165 (24 Aug 2005 19:44:00 GMT)
X-Complaints-To: news@news.stanford.edu
User-Agent: Mozilla Thunderbird 1.0.2 (Windows/20050317)
X-Accept-Language: en-us, en
In-Reply-To: <1124906128.671957@seven.kulnet.kuleuven.ac.be>
Xref: number1.nntp.dca.giganews.com comp.lang.fortran:136324
Bart Vandewoestyne wrote:
> I have a feeling that if we're still looking for a *real* speedup, it will
> have to be something algorithmic rather than these kind of micro-changes...
I tried a bit of one -- unrolling the loop in Elliot's so it has
hardwired checks against all the primes below 100, rather than simply
against 2. That produces the following result (on a 1.1GHz P3):
Enter the number of points you want to time for: 10000000
Elliot's method found 664579 primes in 30.39 seconds.
Bart's method found 664579 primes in 60.50 seconds.
Brooks's method found 664579 primes in 27.73 seconds.
A bit of checking finds that almost 90% of the numbers are trapped
somewhere in the hardwired checks. However, short-circuiting the check
for the remainder cuts my time down to 0.70 seconds, so this is clearly
where the time goes.
So, I hardwire in checks for all primes below 400 (yes, I _am_
autogenerating the code; why do you ask? :) and get this result:
Enter the number of points you want to time for: 10000000
Elliot's method found 664579 primes in 30.38 seconds.
Bart's method found 664579 primes in 60.14 seconds.
Brooks's method found 664579 primes in 23.06 seconds.
That's promising, so I hardwire in checks for all below 1000. (I also
dropped Bart's method from the test, as I was tired of waiting on it.)
Enter the number of points you want to time for: 10000000
Elliot's method found 664579 primes in 30.36 seconds.
Brooks's method found 664579 primes in 16.11 seconds.
Finally, hardwiring in everything below 3300, which should take care of
essentially anything in range. That seems to dramatically increase the
compilation time, but the computed result is substantially lower.
Enter the number of points you want to time for: 10000000
Elliot's method found 664579 primes in 30.27 seconds.
Brooks's method found 664579 primes in 8.08 seconds.
And there you have it; I'll submit that as my entry. Since I figure
nobody will want to retype all of my code (and, at 1600 lines, it's
_far_ too long to put in a Usenet message), I've put it online here:
http://dpdx.net/temp/time_isprime_brooks.f90
Of course, Matlab has us all beat, no doubt with careful assembly code
(and probably with optimizations for the fact that it's calculating
multiple primes). One the same machine:
[color=darkred]
ans = 1 664579
elapsed_time = 2.7030
- Brooks
--
The "bmoses-nospam" address is valid; no unmunging needed.
| |
| Glenn Forney 2005-08-24, 6:59 pm |
| >> I'd guess that they're compiler-specific. A good compiler could turn
>
>
> Why does that help? You substitute one multiply by a multiply and two
> adds, n'est ce pas?
>
> cheers,
>
> Rich
maybe 2*i can be done more cheaply than a general multiply - as either
i+i or a left shift of i
glenn
--
Glenn Forney
National Institute of Standards and Technology
100 Bureau Drive, Stop 8663
Gaithersburg MD 20899-8663
Telephone: (301) 975 2313
FAX: (301) 975 4052
----== Posted via Newsfeeds.Com - Unlimited-Uncensored-Secure Usenet News==----
http://www.newsfeeds.com The #1 Newsgroup Service in the World! 120,000+ Newsgroups
----= East and West-Coast Server Farms - Total Privacy via Encryption =----
| |
| Rich Townsend 2005-08-24, 6:59 pm |
| Glenn Forney wrote:
>
>
> maybe 2*i can be done more cheaply than a general multiply - as either
> i+i or a left shift of i
Oh, of course. D'oh!
cheers,
Rich
| |
| Brooks Moses 2005-08-24, 6:59 pm |
| Brooks Moses wrote:
> Enter the number of points you want to time for: 10000000
> Elliot's method found 664579 primes in 30.27 seconds.
> Brooks's method found 664579 primes in 8.08 seconds.
>
> And there you have it; I'll submit that as my entry. Since I figure
> nobody will want to retype all of my code (and, at 1600 lines, it's
> _far_ too long to put in a Usenet message), I've put it online here:
> http://dpdx.net/temp/time_isprime_brooks.f90
I forgot something; the previously-mentioned version unrolls the
mod(i,p) loop, but it doesn't include any checks if i*i>p, and so it
checks every p against the whole unrolled loop. This takes a good bit
of extra time, and so I put in a few checks that stop the checking and
declare it a prime if p is less than the relevant i*i. That produces
this result:
Enter the number of points you want to time for: 10000000
Elliot's method found 664579 primes in 30.17 seconds.
Bart's method found 664579 primes in 60.00 seconds.
Brooks's method found 664579 primes in 5.22 seconds.
I suspect this can be optimized even further by adding more of the
aforementioned checks, though certainly at some point one has too many
of them.
Again, the code (now revised) is at:
http://dpdx.net/temp/time_isprime_brooks.f90
As a sampling of what it looks like:
! Brooks Moses's function
!
function isprime_moses(p) result (res)
[...]
if (p < 3299) then
[...]
[here, we essentially have isprime_elliot, encapsulated]
[...]
else if (modulo(p,2) == 0) then
res = .false.
return
else if (modulo(p,3) == 0) then
res = .false.
return
[...]
else if (modulo(p,997) == 0) then
res = .false.
return
else if (p < 1000000) then
res = .true.
return
else if (modulo(p,1009) == 0) then
res = .false.
return
[...]
else if (modulo(p,3271) == 0) then
res = .false.
return
else
do i=3299,p,2
if (i*i > p) then
res = .true.
return
else if (modulo(p,i) == 0) then
res = .false.
return
end if
end do
end if
end function isprime_moses
- Brooks
--
The "bmoses-nospam" address is valid; no unmunging needed.
| |
| Bart Vandewoestyne 2005-08-24, 6:59 pm |
| On 2005-08-24, Brooks Moses <bmoses-nospam@cits1.stanford.edu> wrote:
>
> Enter the number of points you want to time for: 10000000
> Elliot's method found 664579 primes in 30.17 seconds.
> Bart's method found 664579 primes in 60.00 seconds.
> Brooks's method found 664579 primes in 5.22 seconds.
Hmm... what was the compiler/environment for this? Using F on my
Linux system, i get no speedup for your method:
bartv@vonneumann:~/fortran$ make time_isprime
F -o time_isprime -O4 time_isprime.f95
bartv@vonneumann:~/fortran$ ./time_isprime
Enter the number of points you want to time for: 10000000
Moses's method found 664579 primes in 13.75 seconds.
Elliot's method found 664579 primes in 13.75 seconds.
Rich's method found 664579 primes in 14.90 seconds.
Bart's method found 664579 primes in 27.18 seconds.
g95 also gives a negligible speedup:
bartv@vonneumann:~/fortran$ make time_isprime
g95 -o time_isprime -O4 time_isprime.f95
bartv@vonneumann:~/fortran$ ./time_isprime
Enter the number of points you want to time for: 10000000
Moses's method found 664579 primes in 21.49 seconds.
Elliot's method found 664579 primes in 21.50 seconds.
Rich's method found 664579 primes in 22.56 seconds.
Bart's method found 664579 primes in 43.44 seconds.
It's also  to see how the *Intel* compiler version 7.1 on my *Intel*
processor is slower than both F and g95 when using -O2:
bartv@vonneumann:~/fortran$ make time_isprime
ifc -o time_isprime -O2 -Tf time_isprime.f95
module MOD_NUMERIC_KINDS
module MOD_ISPRIME
module function ISPRIME_BART
module function ISPRIME_ELLIOT
module function ISPRIME_RICH
module function ISPRIME_MOSES
module subroutine TIME_FUNCTION
program TIME_ISPRIME
1645 Lines Compiled
bartv@vonneumann:~/fortran$ ./time_isprime
Enter the number of points you want to time for: 10000000
Moses's method found 664579 primes in 33.55 seconds.
Elliot's method found 664579 primes in 33.55 seconds.
Rich's method found 664579 primes in 33.30 seconds.
Bart's method found 664579 primes in 66.45 seconds.
Conclusion so far: Eliott and Moses are the `overall' fastest ones, with
Eliott winning for cleaner coding style and lesser number of lines ;-)
Any other suggestions? Let's say the competition ends each Sunday at 23:59
local time :-)
Regards,
Bart
--
"Share what you know. Learn what you don't."
| |
| Bart Vandewoestyne 2005-08-24, 6:59 pm |
| On 2005-08-24, Brooks Moses <bmoses-nospam@cits1.stanford.edu> wrote:
>
> Of course, Matlab has us all beat, no doubt with careful assembly code
> (and probably with optimizations for the fact that it's calculating
> multiple primes). One the same machine:
>
> ans = 1 664579
> elapsed_time = 2.7030
Hmm... i was keeping the primes(...) function for a later
challenge... so let's focus on isprimes(n) for now...
I'm not quite sure if it's fair to compare the matlab code as you
wrote it to what we have written. What we are doing for our
speed test is running over all numbers 1..n and checking if they
are prime.
We could also implement a primes(...) routine (subject of a later
challenge by the way ;-) using a specific prime-searching
algorithm like the Sieve of Eratothenes and we might come up with
a faster routine...
doing
from 1 to n
isprime(n)
end
is not the same as
primes(n)
We are really checking all n, the primes(n) routine could come up
with a cleverer sieve algorithm excluding n-values to be checked... right?
.... testing things in Matlab right now...
Hmm... take this testscript in matlab:
function time_isprime(n)
nbfound = 0;
t = cputime;
for i=1:n,
if isprime(i)
nbfound = nbfound + 1;
end
end;
cputime - t
Using my fortran code, i have:
bartv@vonneumann:~/fortran$ ./time_isprime
Enter the number of points you want to time for: 100000
Moses's method found 9592 primes in 0.05 seconds.
Elliot's method found 9592 primes in 0.06 seconds.
Rich's method found 9592 primes in 0.06 seconds.
Bart's method found 9592 primes in 0.10 seconds.
Using matlab on the same machine, i have:
[color=darkred]
ans =
14.8800
Just in case the cputime intrinsics are not the same: wallclock time is
also considerably higher for my matlab script... so are you sure we're not
beating the Matlab routines?
Regards,
Bart
--
"Share what you know. Learn what you don't."
| |
| Brooks Moses 2005-08-24, 6:59 pm |
| Bart Vandewoestyne wrote:
> On 2005-08-24, Brooks Moses <bmoses-nospam@cits1.stanford.edu> wrote:
>
>
> Hmm... what was the compiler/environment for this? Using F on my
> Linux system, i get no speedup for your method:
>
> bartv@vonneumann:~/fortran$ make time_isprime
> F -o time_isprime -O4 time_isprime.f95
> bartv@vonneumann:~/fortran$ ./time_isprime
> Enter the number of points you want to time for: 10000000
> Moses's method found 664579 primes in 13.75 seconds.
> Elliot's method found 664579 primes in 13.75 seconds.
> Rich's method found 664579 primes in 14.90 seconds.
> Bart's method found 664579 primes in 27.18 seconds.
How very odd. The fact that you're consistently getting the same answer
for my method and Elliot's method makes me wonder if there's something
strange going on in your version of the code. (Particularly since mine
uses the code from Elliot's as a fallback -- if for some reason it were
always accessing the fallback, this is the result I'd expect.) I find
it particularly odd that my absolute times are so much lower than the
ones you get, given that your computer is more than twice as fast as yours.
In any case, I'm using Intel's compiler, version 8.1, with the standard
"Release" compilation configuration. (I'm not sure what optimization
level that corresponds to; it just calls it "Optimize Speed", as opposed
to "Optimize Speed with Additional Optimizations".) The computer is a
1.1MHz P3, running Windows 2003.
- Brooks
--
The "bmoses-nospam" address is valid; no unmunging needed.
| |
| Bart Vandewoestyne 2005-08-24, 6:59 pm |
| On 2005-08-24, Brooks Moses <bmoses-nospam@cits1.stanford.edu> wrote:
>
> How very odd. The fact that you're consistently getting the same answer
> for my method and Elliot's method makes me wonder if there's something
> strange going on in your version of the code. (Particularly since mine
> uses the code from Elliot's as a fallback -- if for some reason it were
> always accessing the fallback, this is the result I'd expect.) I find
> it particularly odd that my absolute times are so much lower than the
> ones you get, given that your computer is more than twice as fast as yours.
*whoops*... I'm sorry... i made a copy-paste error here... these
are my new timings with the F compiler and -O4:
bartv@vonneumann:~/fortran$ make time_isprime
F -o time_isprime -O4 time_isprime.f95
bartv@vonneumann:~/fortran$ ./time_isprime
Enter the number of points you want to time for: 10000000
Moses's method found 664579 primes in 2.54 seconds.
Elliot's method found 664579 primes in 13.75 seconds.
Rich's method found 664579 primes in 14.90 seconds.
Bart's method found 664579 primes in 27.21 seconds.
Moses wins for now... there's time up till Sunday evening for improvement :-)
Btw, Moses... can you explain a bit more to me why your version has this
speedup? I'm not that familiar with loop-unrolling techniques... a simple
example might help me in understanding... or a good reference or something...
Regards,
Bart
--
"Share what you know. Learn what you don't."
| |
| Brooks Moses 2005-08-24, 6:59 pm |
| Bart Vandewoestyne wrote:
> *whoops*... I'm sorry... i made a copy-paste error here... these
> are my new timings with the F compiler and -O4:
>
> bartv@vonneumann:~/fortran$ make time_isprime
> F -o time_isprime -O4 time_isprime.f95
> bartv@vonneumann:~/fortran$ ./time_isprime
> Enter the number of points you want to time for: 10000000
> Moses's method found 664579 primes in 2.54 seconds.
> Elliot's method found 664579 primes in 13.75 seconds.
> Rich's method found 664579 primes in 14.90 seconds.
> Bart's method found 664579 primes in 27.21 seconds.
>
> Moses wins for now... there's time up till Sunday evening for improvement :-)
(Incidentally, my first name is actually Brooks -- though I admit to
having  things by using my last name in the function name, where
you and Elliot were using your first names!)
> Btw, Moses... can you explain a bit more to me why your version has this
> speedup? I'm not that familiar with loop-unrolling techniques... a simple
> example might help me in understanding... or a good reference or something...
This isn't exactly a straight loop-unrolling, but the fundamental
concept is pretty simple. Suppose you have the following loop:
do i=1,500
if modulo(i,p) then
res = .true.
return
end if
end do
Now, for each loop through, there's one addition (to the loop counter),
one comparison (of the loop counter to the loop end), and one set of
loop internals. If you adjust this slightly, to this:
do i=1,2,499
if modulo(i,p) then
res = .true.
return
end if
if modulo(i+1,p) then
res = .true.
return
end if
end do
Doing that saves you one comparison for every two passes through the
original loop. That may or may not be especially significant. What it
also may do is allow the two passes to be pipelined better, though in
this particular case with the embedded return statements, I'm not sure
how the compiler handles it.
You can further completely unroll and hardwire the loop, and write
if modulo(1,p) then
res = .true.
return
end if
if modulo(2,p) then
res = .true.
return
end if
[...]
if modulo(500,p) then
res = .true.
return
end if
This removes all of the additions, and all of the loop-variable
comparisons, at the expense of adding rather a lot of code, and making
it difficult to change things. There are probably fairly few cases
where something like this is useful.
Finally, though, there's an improvement here that is only incidentally
related to the loop-unrolling. To check a number for primeness, we only
need do modulo(i,p)==0 checks for the values of i that are prime.
Looping through every odd number is wasteful, except that in the simple
codes there's no easy way to avoid it. But what we'd really like is to
just loop i through the prime numbers.
Someone had mentioned that it would be an improvement to have the
function save a list of the already-found low prime numbers to use for
this purpose, but that the F-compatible rules prevented using globals to
accomplish this. However, the rules don't prohibit simply finding these
primes ahead of time and hardcoding them into the program!
That's were my code produces its greatest speedup, I suspect. For a
number p around 10000000, the maximum i is about 3300, so the simple
code that checks every odd i does about 1650 modulos and comparisons.
My code only checks the prime numbers, and thus does only 463 modulos
and comparisons.
- Brooks
--
The "bmoses-nospam" address is valid; no unmunging needed.
| |
| glen herrmannsfeldt 2005-08-24, 6:59 pm |
| Brooks Moses wrote:
(snip)
> This isn't exactly a straight loop-unrolling, but the fundamental
> concept is pretty simple. Suppose you have the following loop:
>
> do i=1,500
> if modulo(i,p) then
> res = .true.
> return
> end if
> end do
> Now, for each loop through, there's one addition (to the loop counter),
> one comparison (of the loop counter to the loop end), and one set of
> loop internals. If you adjust this slightly, to this:
>
> do i=1,2,499
> if modulo(i,p) then
> res = .true.
> return
> end if
> if modulo(i+1,p) then
> res = .true.
> return
> end if
> end do
I have seen loop unrolling done with the PL/I preprocessor.
I still wonder sometimes about the decrease in the power of
preprocessors from PL/I to C to Java. It seems that they are
going out of style. In any case, one can do:
%DO I=1 to 500;
IF MOD(I,P)=0 THEN DO;
RES='1'B;
RETURN;
END;
%END;
which is, to the compiler, 500 IF statements. The %DO loop
is done at compile time and the expanded result sent to be
compiled.
One could do something pretty complicated including testing
against some number of primes in a table generated at compile
time by the preprocessor.
One should even be able to generate freeform Fortran from the
PL/I preprocessor, maybe even a large table of flags indicating
which numbers are prime....
-- glen
| |
| e p chandler 2005-08-25, 3:57 am |
| Brooks Moses wrote:
[snip]
> That's were my code produces its greatest speedup, I suspect. For a
> number p around 10000000, the maximum i is about 3300, so the simple
> code that checks every odd i does about 1650 modulos and comparisons.
> My code only checks the prime numbers, and thus does only 463 modulos
> and comparisons.
While trying to better your time and still write a legal F program, I
looked for other ways to determine primality.
Remembering a bit of number theory, I realized that while interesting,
Wilson's Theorem would be of no practical use.
see http://primes.utm.edu/notes/proofs/Wilsons.html
However, while browsing one of the other fine pages at this web site
http://primes.utm.edu/
I found an application of the Little Fermat Therorem
http://primes.utm.edu/notes/proofs/...tleTheorem.html
to testing primality
http://primes.utm.edu/prove/prove2_3.html
This involves testing whether a**(p-1) = 1 (mod p) but only for
selected values of a depending on an upper bound for p.
What amazed me is that if p is small enough, only {2,3} need be tested.
For p in the range of the original program (1e7) only {2,3,5} need be
tested.
This changes the situation from p being a probable prime to p certainly
being a prime since the test is conducted for numbers in a limited
range.
Here are my timing results - 2.8 GHz P4-HT, WinXP(SP2):
C:\temp>make
C:\temp>F -O4 test.f95 -o test.exe
C:\temp>test
Enter the number of points you want to time for: 10000000
Elliot's method found 664579 primes in 11.69 seconds.
Bart's method found 664579 primes in 23.16 seconds.
Brooks's method found 664579 primes in 3.80 seconds.
Pierre's method found 664579 primes in 9.06 seconds.
C:\temp>
Although it was slower than your program, I was not really
disappointed, since the program actually worked without using extensive
divisibility testing. There was one gotcha - it requires double length
integer arithmetic (8 byte integers).
The following must be added to the declarations at the top of the
program:
public :: isprime_fermat
public :: calc_power
Here's the rest of the program listing:
! replace later divisibility tests by strong probable prime test
! based on little Fermat theorem
! p is .le. 1e7 so testing a={2,3,5} suffices to prove primality
function isprime_fermat(p) result (res)
integer(kind=i4b), intent(in) :: p
logical :: res
integer(kind=i4b) :: a, s, d, t, r
res = .false.
if (p < 2) then
res = .false.
return
else if (p == 2) then
res = .true.
return
else if (modulo(p,2) == 0) then
res = .false.
return
else if (p == 3) then
res = .true.
return
else if (modulo(p,3) == 0) then
res = .false.
return
else if (p == 5) then
res = .true.
return
else if (modulo(p,5) == 0) then
res = .false.
return
else if (p == 7) then
res = .true.
return
else if (modulo(p,7) == 0) then
res = .false.
return
else
! find d * (2 ** s) = p - 1 where d is odd
s = 0
d = p - 1
do
if (modulo(d,2) == 1) exit
s = s + 1
d = d / 2
end do
end if
res = .false.
a = 2
t = calc_power(a, d, p) ! a ** d
if (t == 1) res = .true.
if (.not. res) then
do r = 0, s - 1
if (t == p - 1) then
res = .true. ! (a ** d) ** ( 2 ** r)
exit
end if
t = calc_power(t, 2, p)
end do
end if
if (.not. res) return
res = .false.
a = 3
t = calc_power(a, d, p) ! a ** d
if (t == 1) res = .true.
if (.not. res) then
do r = 0, s - 1
if (t == p - 1) then
res = .true. ! (a ** d) ** ( 2 ** r)
exit
end if
t = calc_power(t, 2, p)
end do
end if
if (.not. res) return
res = .false.
a = 5
t = calc_power(a, d, p) ! a ** d
if (t == 1) res = .true.
if (.not. res) then
do r = 0, s - 1
if (t == p - 1) then
res = .true. ! (a ** d) ** ( 2 ** r)
exit
end if
t = calc_power(t, 2, p)
end do
end if
return
end function isprime_fermat
function calc_power(a, d, p) result (r)
integer, parameter :: i8b = selected_int_kind(14)
integer(kind=i4b), intent(in) :: a, d, p
integer(kind=i4b) :: r
integer(kind=i4b) :: b
integer(kind=i8b) :: m, t, e
! calculate a**d mod p
! generate successive squares of a
! multiplying the total by squares corresponding
! to 1 bits in the binary representation of d
b = d
e = a
t = 1
m = p
do
if (modulo(b,2) == 1) then
t = t * e
t = modulo(t,m)
end if
e = e * e
e = modulo(e,m)
b = b / 2
if (b == 0) exit
end do
r = t
return
end function calc_power
--------
Elliot
e-mail: epc8@juno.com
| |
| Brooks Moses 2005-08-25, 3:57 am |
| e p chandler wrote:
> I found an application of the Little Fermat Therorem
>
> http://primes.utm.edu/notes/proofs/...tleTheorem.html
>
> to testing primality
>
> http://primes.utm.edu/prove/prove2_3.html
>
> This involves testing whether a**(p-1) = 1 (mod p) but only for
> selected values of a depending on an upper bound for p.
>
> What amazed me is that if p is small enough, only {2,3} need be tested.
> For p in the range of the original program (1e7) only {2,3,5} need be
> tested.
[...]
> Brooks's method found 664579 primes in 3.80 seconds.
> Pierre's method found 664579 primes in 9.06 seconds.
That definitely seems promising. One thing that I noticed, in
programing my version, is that checking for divisibility by the prime
numbers below 100 would eliminate 90% of the numbers, and took only 0.70
seconds (on my computer, which would probably translate to about 0.5
seconds on yours) if I skipped everything after that. For the "brute
force" algorithms, the vast majority of the time is spent on the primes
and the numbers with very large divisors.
On the other hand, your method looks like it takes within a factor of
three of the same amount of time for all numbers (except ones divisible
by 2, 3, 5, or 7), regardless of primality. So it may well be faster
for those numbers -- particularly if there's a bit of optimization to be
done.
Thus, I'm wondering if a hybrid approach, which does hardwired modulo
checks against all the primes below some N (for N around 100 or 200,
say) and then uses your method to test primality for whatever gets
through that, would be fast enough to beat my method. I wouldn't be at
all surprised if it did, but I don't want to stay up late enough tonight
to write it!
- Brooks
P.S. How small does the number have to be for one to only need to check
against (2,3)? If that's any significant fraction of 10000000, it may
be worth putting in a check that calls p prime after only those checks
(rather than checking against 5 as well) if p is small enough.
--
The "bmoses-nospam" address is valid; no unmunging needed.
| |
| Bart Vandewoestyne 2005-08-25, 3:57 am |
| On 2005-08-24, Brooks Moses <bmoses-nospam@cits1.stanford.edu> wrote:
>
> [...]
> That's were my code produces its greatest speedup, I suspect. For a
> number p around 10000000, the maximum i is about 3300, so the simple
> code that checks every odd i does about 1650 modulos and comparisons.
> My code only checks the prime numbers, and thus does only 463 modulos
> and comparisons.
Maybe I should add this: in the ideal case, the range
n=1..huge(n) would be the range to be checked (*). I just picked
10000000 because it didn't run that long...
I don't know if this influences the solutions much...
Regards,
Bart
(*) Actually, the really really ideal case would be 1..\infty,
but then we would have found all possible prime numbers :-)
--
"Share what you know. Learn what you don't."
| |
| Michel OLAGNON 2005-08-25, 3:57 am |
|
Bart Vandewoestyne wrote:
> On 2005-08-24, Brooks Moses <bmoses-nospam@cits1.stanford.edu> wrote:
>
>
>
> Maybe I should add this: in the ideal case, the range
> n=1..huge(n) would be the range to be checked (*). I just picked
> 10000000 because it didn't run that long...
Then, why not n=huge(n)-10000000, huge(n)-1 ?
| |
| Bart Vandewoestyne 2005-08-25, 3:57 am |
| On 2005-08-25, e p chandler <epc8@juno.com> wrote:
>
> [...]
> There was one gotcha - it requires double length integer arithmetic (8 byte
> integers).
>
> [...]
>
> function calc_power(a, d, p) result (r)
> integer, parameter :: i8b = selected_int_kind(14)
This is going to be off-topic, but I do not understand why you
only put in 14...
To setup my integer kinds, i have always reasoned as follows:
! 1 byte => max 2^8 = 2.56e+02 so 10^2 is assured
integer, parameter, public :: i1b = selected_int_kind(2)
! 2 bytes => max 2^16 = 6.5536e+04, so 10^4 is assured
integer, parameter, public :: i2b = selected_int_kind(4)
! 4 bytes => max 2^32 = 4.294967296e+09, so 10^9 is assured
integer, parameter, public :: i4b = selected_int_kind(9)
where the 'max' values represent the maximum value to be represented in base 2
if you have that many digits available.
So, my reasoning for 8-byte integers would be:
! 8 bytes => max 2^64 = 1.8446744073709551616e+19, so 10^19 is assured
integer, parameter, public :: i8b = selected_int_kind(19)
but that appears not to be available on my Intel(R) Pentium(R) 4 CPU 2.40GHz
and compiling with my Fortran Company/NAG F compiler Release 20031017.
integer, parameter, public :: i8b = selected_int_kind(14)
however appears to be available.
Can somebody explain me why the 14 is for 8-byte integers instead of the 19?
I think it is because I should not reason in the `school model' of 8 bytes
being really 64 bits that I can set. I guess in the `real' model, the 64 bits
do not only contain the bits of the integer value, but also a sign bit and
maybe some other extra bits are used for other purposes?
=> ... doing some checking and experimenting now... <=
Now that I think of it... there's probably one sign bit, so that
would mean that there are still 63 bits left for the values.
This would make my reasoning as follows:
! 1 byte => max 2^7 = 1.28e+02 so 10^2 is assured
! 2 bytes => max 2^15 = 3.2768e+04, so 10^4 is assured
! 4 bytes => max 2^31 = 2.147483648e+09, so 10^9 is assured
! 8 bytes => max 2^63 = 9.223372036854775808e+18, so 10^18 is assured
And this seems to be ok, because selected_int_kind(18) is indeed available
on my setup.
Is this reasoning OK and is it indeed because of the sign bit that i have
10^18 assured instead of 10^19?
By the way, what would be the document to look for to see how the integers
are represented on my architecture? Simply an Intel-document or something?
For floating point there's `IEEE Standard 754', but is there a similar
standard for integers or does one simply have to look in the
manuals for the architecture one's working on?
Regards,
Bart
--
"Share what you know. Learn what you don't."
| |
| Janne Blomqvist 2005-08-25, 8:00 am |
| Bart Vandewoestyne wrote:
> By the way, what would be the document to look for to see how the integers
> are represented on my architecture? Simply an Intel-document or something?
> For floating point there's `IEEE Standard 754', but is there a similar
> standard for integers or does one simply have to look in the
> manuals for the architecture one's working on?
No official standard AFAIK, but the so-called "two's complement"
representation is pretty much the defacto standard these days
(e.g. x86 uses it). Google and ye shall find. I'm sure architecture
manuals also mention it.
--
Janne Blomqvist
| |
| Bart Vandewoestyne 2005-08-25, 8:00 am |
| On 2005-08-25, Michel OLAGNON <molagnon@ifremer-a-oter.fr> wrote:
>
>
> Then, why not n=huge(n)-10000000, huge(n)-1 ?
Hmm... this is indeed an interesting remark... i have changed the
skeleton a bit so you now can select whether you want to check n
numbers starting from one, or n numbres up to huge(n)-1.
The code is at
http://www.cs.kuleuven.ac.be/~bartv...ime_isprime.f95
Interesting observations appear... checking 10000000 numbers from 1 and up
favor's Brooks method:
bartv@vonneumann:~/fortran$ ./time_isprime
How many numbers should be tested? 10000000
What method (from_1|from_huge) ? from_1
from_1: Brooks's method found 664579 primes in 2.55 seconds.
from_1: Fermat's method found 664579 primes in 8.24 seconds.
from_1: Elliot's method found 664579 primes in 13.81 seconds.
from_1: Rich's method found 664579 primes in 14.87 seconds.
from_1: Bart's method found 664579 primes in 27.19 seconds.
However, if we do the same test for the last 25 numbers of our
integer range, I get:
bartv@vonneumann:~/fortran$ ./time_isprime
How many numbers should be tested? 25
What method (from_1|from_huge) ? from_huge
from_huge: Brooks's method found 1 primes in 4.52 seconds.
from_huge: Fermat's method found 1 primes in 0.00 seconds.
from_huge: Elliot's method found 1 primes in 4.57 seconds.
from_huge: Rich's method found 1 primes in 0.00 seconds.
from_huge: Bart's method found 1 primes in 9.20 seconds.
.... Hmm... although now that I think of it... although Fermat and
Rich's method seem extremely fast now, the 0.00 is probably due to overflow...
or not?
So far, I would go for Elliot's method since it's quite readable
(Sorry Brooks ;-), it doesn't require 8 byte integers, and it doesn't suffer
from the overflow problems that Rich and Fermat's method probably suffer
from...
But there's ofcourse still time for improvement until Sunday evening... ;-)
Regards,
Bart
--
"Share what you know. Learn what you don't."
| |
| Bart Vandewoestyne 2005-08-25, 8:00 am |
| On 2005-08-25, Janne Blomqvist <foo@bar.invalid> wrote:
>
> No official standard AFAIK, but the so-called "two's complement"
> representation is pretty much the defacto standard these days
> (e.g. x86 uses it). Google and ye shall find. I'm sure architecture
> manuals also mention it.
OK. This makes sense. I'm glad some cs-principles I've learned long
time ago finally come into my practice ;-)
Let me phrase it in my own words:
Because Intel uses two's complement notation, the integers come
from the closed interval
[-2^{n-1}...2^{n-1}-1]
where n is the number of bits used. Therefore, the maximum absolute
value that can occur is 2^{n-1}, which leads to the following table:
1 byte => max 2^7 = 1.28e+02 so 10^2 is assured
2 bytes => max 2^15 = 3.2768e+04, so 10^4 is assured
4 bytes => max 2^31 = 2.147483648e+09, so 10^9 is assured
8 bytes => max 2^63 = 9.223372036854775808e+18, so 10^18 is assured
and this matches with the selected_int_kind values I am expecting.
Regards,
Bart
--
"Share what you know. Learn what you don't."
| |
| e p chandler 2005-08-25, 8:00 am |
| Brooks Moses wrote:
> e p chandler wrote:
> [...]
>
> That definitely seems promising. One thing that I noticed, in
> programing my version, is that checking for divisibility by the prime
> numbers below 100 would eliminate 90% of the numbers, and took only 0.70
> seconds (on my computer, which would probably translate to about 0.5
> seconds on yours) if I skipped everything after that. For the "brute
> force" algorithms, the vast majority of the time is spent on the primes
> and the numbers with very large divisors.
>
> On the other hand, your method looks like it takes within a factor of
> three of the same amount of time for all numbers (except ones divisible
> by 2, 3, 5, or 7), regardless of primality. So it may well be faster
> for those numbers -- particularly if there's a bit of optimization to be
> done.
>
> Thus, I'm wondering if a hybrid approach, which does hardwired modulo
> checks against all the primes below some N (for N around 100 or 200,
> say) and then uses your method to test primality for whatever gets
> through that, would be fast enough to beat my method. I wouldn't be at
> all surprised if it did, but I don't want to stay up late enough tonight
> to write it!
>
The authors at one of the web pages I referenced recommended N of 257.
> - Brooks
>
> P.S. How small does the number have to be for one to only need to check
> against (2,3)? If that's any significant fraction of 10000000, it may
> be worth putting in a check that calls p prime after only those checks
> (rather than checking against 5 as well) if p is small enough.
>
Yes, it might speed the program up a bit. I'm going to try putting in
the first check listed below.
Let me extract one of the tables from the Strong PRP page:
if n < 1,373,653 then a={2,3}
if n < 25,326,001 then a={2,3,5}
if n < 118,670,087,467) and a={2,3,5,7} then either n=3,215,031,751 or
n is prime
......
some other conditions for n < 1e12
While these are large numbers, IIRC, they are still considered small by
cryptographic standards. :-).
| |
| e p chandler 2005-08-25, 8:00 am |
| Brooks Moses wrote:
> P.S. How small does the number have to be for one to only need to check
> against (2,3)? If that's any significant fraction of 10000000, it may
> be worth putting in a check that calls p prime after only those checks
> (rather than checking against 5 as well) if p is small enough.
Making this one modification changed the run time from 9.06 to 8.84
seconds.
I think I'll try adding more trial divisors, but that's going to be it.
-------
Elliot
| |
| e p chandler 2005-08-25, 8:00 am |
|
e p chandler wrote:
> Brooks Moses wrote:
>
>
> Making this one modification changed the run time from 9.06 to 8.84
> seconds.
>
> I think I'll try adding more trial divisors, but that's going to be it.
>
[excuse reply to self]
Adding in trial divisors from 11 to 251 brought the run time down to
6.36 seconds.
Note that I am using each trial divisor *twice*. Once to test if p is
that prime and again to test if p is divisible by that prime.
It's not obvious to me how you avoided using each prime twice in your
program. [ violating dictum "Never program when sleepy." :-) ]
-----
Elliot
| |
| glen herrmannsfeldt 2005-08-25, 6:59 pm |
| Janne Blomqvist wrote:
> Bart Vandewoestyne wrote:
[color=darkred]
> No official standard AFAIK, but the so-called "two's complement"
> representation is pretty much the defacto standard these days
> (e.g. x86 uses it). Google and ye shall find. I'm sure architecture
> manuals also mention it.
As I understand it, Unisys is still building some ones complement
machines. CDC CYBER machines were also ones complement, but those
may all be gone by now.
Otherwise, there is endianness. The bytes of the number may be
stored in memory either with the least significant first
(little endian) or most significant first (big endian).
Floating point data, while still IEEE 754, may also be stored
little or big endian, and possibly different than integers.
-- glen
| |
| Brooks Moses 2005-08-25, 6:59 pm |
| e p chandler wrote:
> Adding in trial divisors from 11 to 251 brought the run time down to
> 6.36 seconds.
>
> Note that I am using each trial divisor *twice*. Once to test if p is
> that prime and again to test if p is divisible by that prime.
>
> It's not obvious to me how you avoided using each prime twice in your
> program. [ violating dictum "Never program when sleepy." :-) ]
At the beginning of the whole thing:
if p < largest trial divisior
do simple brute-force test for primality
[else ...]
Given that the square root of the largest trial divisor is around 50 or
so even with my huge list, the brute-force test takes a trivial amount
of time, and there's no point in optimizing further.
- Brooks
--
The "bmoses-nospam" address is valid; no unmunging needed.
| |
| Brooks Moses 2005-08-25, 6:59 pm |
| Bart Vandewoestyne wrote:
> Hmm... this is indeed an interesting remark... i have changed the
> skeleton a bit so you now can select whether you want to check n
> numbers starting from one, or n numbres up to huge(n)-1.
[...]
> So far, I would go for Elliot's method since it's quite readable
> (Sorry Brooks ;-), it doesn't require 8 byte integers, and it doesn't suffer
> from the overflow problems that Rich and Fermat's method probably suffer
> from...
I think that requiring that the method work for everything up to huge(n)
without overflow is a fairly significant constraint on the methods,
since it severely limits the mathematical operations that one can do
with the number one is testing.
Do the non-brute-force methods work if one uses something like
huge(n)/32 rather than huge(n)?
(I also must -- in good humor and no real seriousness, mind you! --
protest the changing of the judging criteria after the contest was
announced, from "as fast as possible" to "fastest without being
'excessively long'". And I hereby enact that protest by including by
reference all versions of my program for highest hardwired divisors less
than 3300. The ones with slightly more than the first four primes will
be faster than Elliot's method and still readable, and thus there must
be an optimum score somewhere between there and the one of mine that you
already have. :)
Actually, if you were going to test for all the numbers between 1 and
huge(n), I'd be strongly tempted to write a code autogenerator, and let
you include the compile time in with the score. The autogenerator would
be quite readable, and the compile-time negligible by comparison.
- Brooks
--
The "bmoses-nospam" address is valid; no unmunging needed.
| |
| e p chandler 2005-08-25, 6:59 pm |
|
Brooks Moses wrote:
> e p chandler wrote:
>
> At the beginning of the whole thing:
>
> if p < largest trial divisior
> do simple brute-force test for primality
> [else ...]
>
> Given that the square root of the largest trial divisor is around 50 or
> so even with my huge list, the brute-force test takes a trivial amount
> of time, and there's no point in optimizing further.
>
> - Brooks
Ouch. I forgot about the enclosing conditional.
I think it's time to put this program to bed.
| |
| Ian Chivers 2005-08-25, 6:59 pm |
| the following might be of interest
program primes
implicit none
integer :: n = 10000000
integer :: i
integer :: k
integer , dimension(8) :: timing
logical, allocatable, dimension(:) :: prime
write(*,*) ' value of n (' , n , ')'
read (*,*) n
call date_and_time(values=timing)
print *,' initial
',timing(6),timing(7),timing(8)
allocate (prime(n))
prime(:) = .true.
do i = 2, n/2
if (prime(i)) then
prime(2*i:n:i) = .false.
endif
enddo
k = count(prime)
call date_and_time(values=timing)
print *,' after
',timing(6),timing(7),timing(8)
print *,k
end
could this be included in the test harness?
"e p chandler" <epc8@juno.com> wrote in message
news:1125003402.610068.128590@g44g2000cwa.googlegroups.com...
>
> Brooks Moses wrote:
>
> Ouch. I forgot about the enclosing conditional.
>
> I think it's time to put this program to bed.
>
| |
| Mark B Stucky 2005-08-25, 6:59 pm |
| I realize that this is a Fortran challenge and not a Matlab challenge, but
just to be "fair" to Matlab, I reworked your time_isprime routine in a
quick attempt to vectorize it :
function [nbfound] = time_isprime_v(n)
nbfound = 0;
t = cputime;
irng = 1:n;
result = isprime(irng);
nbfound = sum(result);
cputime - t
and got the following :
ans = 14.6250
n1 = 9592
[color=darkred]
ans = 3.4370
n1 = 9592
Now back to the Fortran challenge.... :)
--Mark
"Bart Vandewoestyne" <MyFirstName.MyLastName@telenet.be> wrote in message
news:1124917662.124456@seven.kulnet.kuleuven.ac.be...[color=darkred]
> On 2005-08-24, Brooks Moses <bmoses-nospam@cits1.stanford.edu> wrote:
>
> Hmm... i was keeping the primes(...) function for a later
> challenge... so let's focus on isprimes(n) for now...
>
> I'm not quite sure if it's fair to compare the matlab code as you
> wrote it to what we have written. What we are doing for our
> speed test is running over all numbers 1..n and checking if they
> are prime.
>
> We could also implement a primes(...) routine (subject of a later
> challenge by the way ;-) using a specific prime-searching
> algorithm like the Sieve of Eratothenes and we might come up with
> a faster routine...
>
> doing
>
> from 1 to n
> isprime(n)
> end
>
> is not the same as
>
> primes(n)
>
> We are really checking all n, the primes(n) routine could come up
> with a cleverer sieve algorithm excluding n-values to be checked... right?
>
> ... testing things in Matlab right now...
>
> Hmm... take this testscript in matlab:
>
> function time_isprime(n)
>
> nbfound = 0;
> t = cputime;
> for i=1:n,
> if isprime(i)
> nbfound = nbfound + 1;
> end
> end;
> cputime - t
>
>
> Using my fortran code, i have:
>
> bartv@vonneumann:~/fortran$ ./time_isprime
> Enter the number of points you want to time for: 100000
> Moses's method found 9592 primes in 0.05 seconds.
> Elliot's method found 9592 primes in 0.06 seconds.
> Rich's method found 9592 primes in 0.06 seconds.
> Bart's method found 9592 primes in 0.10 seconds.
>
> Using matlab on the same machine, i have:
>
>
> ans =
>
> 14.8800
>
>
> Just in case the cputime intrinsics are not the same: wallclock time is
> also considerably higher for my matlab script... so are you sure we're not
> beating the Matlab routines?
>
> Regards,
> Bart
>
> --
> "Share what you know. Learn what you don't."
| |
| Bart Vandewoestyne 2005-08-26, 3:56 am |
| On 2005-08-25, Mark B Stucky <mstucky5@cox.net> wrote:
> I realize that this is a Fortran challenge and not a Matlab challenge, but
> just to be "fair" to Matlab, I reworked your time_isprime routine in a
> quick attempt to vectorize it :
>
> function [nbfound] = time_isprime_v(n)
>
> nbfound = 0;
> t = cputime;
> irng = 1:n;
> result = isprime(irng);
> nbfound = sum(result);
> cputime - t
>
> and got the following :
>
>
> ans = 14.6250
>
> n1 = 9592
>
>
> ans = 3.4370
>
> n1 = 9592
Is that completely `fair'? I haven't checked how hard it would
be to make our is_prime an elementalfunction so it can act on arrays also,
but our current is_prime only accepts scalars so i thought it
would be fair to also just use a scalar version in my Matlab
script... what you do is kind of eleminating the loop around the
isprime function... i just thought it was fair to also keep that
loop in the Matlab script...
I could be wrong, and if I am then I'm glad to hear why... i always like to
learn from others and their comments on what I write ;-)
Regards,
Bart
--
"Share what you know. Learn what you don't."
| |
| Mark B Stucky 2005-08-26, 6:59 pm |
|
"Bart Vandewoestyne" <MyFirstName.MyLastName@telenet.be> wrote in message
news:1125041535.254177@seven.kulnet.kuleuven.ac.be...
> On 2005-08-25, Mark B Stucky <mstucky5@cox.net> wrote:
>
> Is that completely `fair'? I haven't checked how hard it would
> be to make our is_prime an elementalfunction so it can act on arrays also,
> but our current is_prime only accepts scalars so i thought it
> would be fair to also just use a scalar version in my Matlab
> script... what you do is kind of eleminating the loop around the
> isprime function... i just thought it was fair to also keep that
> loop in the Matlab script...
>
> I could be wrong, and if I am then I'm glad to hear why... i always like
> to
> learn from others and their comments on what I write ;-)
>
> Regards,
> Bart
My only reason for removing the loop was because Matlab "for" loops
are slow and I that you were trying to compare the Matlab "isprime"
function with the Fortran codes. I see the point of your reason for
using a scalar version, but it seemed to me that the version with the
"for" loop was more of a comparison between the interpreted Matlab
"for" loop with the compiled Fortran "do" loop.
Anyhow, it seems that we're getting OT for this Fortran isprime
challenge.
--Mark
| |
| Bart Vandewoestyne 2005-08-26, 6:59 pm |
| On 2005-08-26, Mark B Stucky <mstucky5@cox.net> wrote:
>
> My only reason for removing the loop was because Matlab "for" loops
> are slow and I that you were trying to compare the Matlab "isprime"
> function with the Fortran codes. I see the point of your reason for
> using a scalar version, but it seemed to me that the version with the
> "for" loop was more of a comparison between the interpreted Matlab
> "for" loop with the compiled Fortran "do" loop.
OK.
> Anyhow, it seems that we're getting OT for this Fortran isprime
> challenge.
True, but anyhow... there's still about 46 hours 15 minutes left
(in Belgian time that is, the Americans have 6 or 7 more hours ;-)
Regards,
Bart
--
"Share what you know. Learn what you don't."
| |
| Bart Vandewoestyne 2005-08-29, 7:56 am |
| On 2005-08-23, Bart Vandewoestyne <MyFirstName.MyLastName@telenet.be> wrote:
> I know there has been a kind of `string-challenge' in this group
> which I have not been following... but it's been taking a while
> and from what I've read I have the feeling that some people got
> bored of it and now are definitely in need for something else...
> so I have another little `isprime-challlenge' :-)
>
> It's just 100% pure for the fun of it, and maybe it will also
> teach some more novices and `less experienced' young fortran
> programmers as me good optimization skills. It might also be usefull for
> some people in need for some prime routines.
OK. The little isprime challenge has come to an end. My final testprogram is
at http://www.cs.kuleuven.ac.be/~bartv...ime_isprime.f95
These are the timing results on my computer:
bartv@archimedes:~/fortran$ make time_isprime
F -o time_isprime -O4 time_isprime.f95
bartv@archimedes:~/fortran$ ./time_isprime
How many numbers should be tested? 1000000
What method (from_1|from_huge) ? from_1
from_1: Brooks's method found 78498 primes in 0.15 seconds.
from_1: Fermat's method found 78498 primes in 0.74 seconds.
from_1: Elliot's method found 78498 primes in 0.56 seconds.
from_1: Rich's method found 78498 primes in 0.61 seconds.
from_1: Bart's method found 78498 primes in 1.08 seconds.
How many numbers should be tested? 1000000
What method (from_1|from_huge) ? from_huge
from_huge: Brooks's method found 46602 primes in 66.60 seconds.
from_huge: Fermat's method found 46602 primes in 1.08 seconds.
from_huge: Elliot's method found 46602 primes in 69.20 seconds.
from_huge: Rich's method found 46602 primes in 21.62 seconds.
from_huge: Bart's method found 46602 primes in 98.68 seconds.
I'll leave it to the reader to decide which version is preferable. What I
have learned from this challenge is the following:
* First of all, look for a fast theoretical algorithm :-)
* The loop-unrolling-alike technique as it is used by Brooks can speedup
dramatically
* I wonder if Fermat's method suffers from overflow problems since we are
testing up until huge(n)-1 and this algorithm requires 8-byte integers...
that would maybe explain the speed of the method, but apparently, it finds
as much primes as the other ones so it seems correct? I haven't looked at
the value of the primes yet, but this could be something for further
investigation... :-)
If I had to choose, I think I would go for Elliot's method because it doesn't
require 8-byte integers and it is definitely more readable than Brooks's method ;-)
See you in the next challenge (coming soon... :-)!
Bart
--
"Share what you know. Learn what you don't."
| |
| billyboyer@fastmail.fm 2005-08-29, 7:01 pm |
| READ THIS! I'm Glad I did! Best way for Quick easy cash, for REAL!
LOTS OF CASH, FAST AND COMPLETELY LEGAL, THIS REALLY WORKS!!
THIS REALLY CAN MAKE YOU EASY MONEY!! IT WORKS!!! BUT YOU
HAVE TO FOLLOW IT TO THE LETTER FOR IT TO WORK!!!!
A little while back, I was browsing through newsgroups, just like you
are now, and came across an article similar to this that said you could
make thousands of dollars within ws with only an initial investment
of $6.00! So I thought," Yeah, right, this must be a scam", but like
most of us, I was curious, so I kept reading. Anyway, it said that you
send $1.00 to each of the 6 names and
address stated in the article. You then place your own name and address
in the bottom of
the list at #6, and post the article in at least 200 newsgroups. (There
are thousands) No
catch, that was it. So after thinking it over, and talking to a few
people first, I thought about trying it. I figured what have I got to
lose except 6 stamps and $6.00, right? Like most
of us I was a little skeptical and a little worried about the legal
aspects of it all. So I checked it out with the U.S. Post Office
(1-800-725-2161)
and they confirmed that it is indeed legal! Then I invested the measly
$6.00. Well GUESS WHAT!!... within 7 days, I started getting money in
the mail! I was shocked! I figured
it would end soon, but the money just kept coming in. In my first w,
I made about $25.00.
By the end of the second w I had made a total of over $1,000.00! In
the third w I
had over $10,000.00 and it's still growing. This is now my fourth w
and I
have made a total of just over $42,000.00 and it's still coming in
rapidly. It's certainly worth $6.00, and 6 stamps, I have spent more
than that on
the lottery!! Let me tell you how this works and most importantly, why
it works....also, make sure you print a copy of this article NOW, so
you can get the information off
of it as you need it. A cousin of mine did this a few months back, and
he is now enjoying life in his brand new Tahoe Sport!! ($30,000+). So
it does work!!
STEP 1: Get 6 separate pieces of paper and write the following on each
piece of paper "PLEASE PUT ME ON YOUR MAILING LIST." Now get 6 US $1.00
bills and place ONE inside EACH of the 6 pieces of paper so the bill
will not be seen
through the envelope to prevent thievery. Next, place one paper in each
of the 6
envelopes and seal them. You should now have 6 sealed envelopes, each
with a piece of paper stating the above phrase, your name and address,
and a
$1.00 bill. What you are doing is creating a service by this. THIS IS
ABSOLUTELY LEGAL! Mail the 6 envelopes to the following addresses:
1- A.L.
: Rimon 631/24
: Natzareth illit, T.D. 13123
: ISRAEL (17000)
:
: 2- Jared Franklin
: PO Box 11
: McCune, KS 66753
: USA
:
: 3- Wallden Enterprises Inc.
: 14329 - 68th Avenue, Surrey, B.C. V3W 2H5
: Canada
:
: 4- Mark Gomez
: 11108 Shoreline Dr.
: El Paso, TX 79936
: USA
:
: 5- Robin Hanson
: PO Box 596
: Virden MB.
: ROM 2CO CANADA
6- Bill Boyer
: 8776 Fern Glen Dr.
: St. Louis, MO 63126
: USA
:
STEP 2: Now take the #1 name off the list that you see above, move the
other names up (6 becomes 5, 5 becomes 4, etc...) and add YOUR Name as
number 6 on
the list. STEP 3: Change anything you need to, but try to keep this
article as close
to original as possible. Now, post your amended article to at least 200
newsgroups. (I think there are close to 24,000 groups) All you need is
200, but remember, the more you post,
the more money you make! ---DIRECTIONS -----HOW TO POST TO
NEWSGROUPS------------ Step 1) You do not need to re-type this entire
letter to do your own
posting. Simply put your cursor at the beginning of this letter and
drag your cursor to the bottom of this document, and select 'copy' from
the
edit menu. This will copy the entire letter into the computers memory.
Step 2) Open a blank
'notepad' file and place your cursor at the top of the blank page. From
the 'edit' menu select
'paste'. This will paste a copy of the letter into notepad so that you
can add your name to
the list. Step 3) Save your new notepad file as a .txt file. If you
want to do your postings
in different sittings, you'll always have this file to go back to. Step
4) Use Netscape or Internet explorer and try searching for
various newsgroups (on-line forums, message boards, chat sites,
discussions.) Step 5) Visit these
message boards and post this article as a new message by highlighting
the text of this
letter and selecting paste from the edit menu. Fill in the Subject,
this will be the header that
everyone sees as they scroll through the list of postings in a
particular group, click the
post message button. You're done with your first one!
Congratulations...THAT'S IT! All you have to do is jump to different
newsgroupes and post away, after you get the hang of it, it will take
about 30 seconds for each newsgroup!
**REMEMBER, THE MORE NEWSGROUPS YOU POST IN, THE MORE
MONEY YOU WILL MAKE!! BUT YOU HAVE TO POST A MINIMUM OF
200** That's it! You will begin reciving money from around the world
within days! You may eventually wany to rent a P.O.Box due to the large
amount of mail you will receive. If
you wish to stay anonymous, you can invent a name to use, as
long as the postman will deliver it. **JUST MAKE SURE ALL THE
ADDRESSES ARE CORRECT.** Now the WHY part: Out of 200 postings, say I
receive only 5 replies (a very low example). So then I made $5.00 with
my name at #6 on the letter. Now,
each of the 5 persons who just sent me $1.00 make the MINIMUM 200
postings, each with my name at #5 and only 5 persons respond to each of
the original 5, that is
another $25.00 for me, now those 25 each make 200 MINIMUM posts with my
name at #4 and only
5 replies each, I will bring in an additional $125.00! Now, those 125
persons turn around and post the MINIMUM 200
with my name at #3 and only receive 5 replies each, I will make an
additional $626.00! OK, now here is the fun part, each of those 625
persons post a MINIMUM 200
letters with my name at #2 and they each only receive 5 replies, that
just made me $3,125.00!!! Those 3,125 persons will all deliver this
message to 200 newsgroups with my
name at #1 and if still 5 persons per 200 newsgroups react I will
receive $15,625,00! With a original investment of only $6.00! AMAZING!
When your name is no longer on the
list, you just take the latest posting in the newsgroups, and send out
another $6.00
to names on the list, putting your name at number 6 again. And start
posting again. The thing to remember is, do you realize that thousands
of people all over the world are joining the internet and reading these
articles everyday, JUST LIKE YOU are now!! So can you afford $6.00 and
see if it really works?? I think so... People have said, "what if the
plan is played out and no one sends you the money? So what! What are
the chances of that happening when there
are tons of new honest users and new honest people who are joining the
internet and newsgroups everyday and are willing to give it a try?
Estimates are at 20,000 to 50,000
new users, every day, with thousands of those joining the actual
internet. Remember, play FAIRLY and HONESTLY and this will work. Also,
for the skeptics, you're looking at this and thinking about doing it
right? imagine how many millions are looking and thinking the
EXACT same thing!!
: Enjoy the money,
: Rick
| |
| Brooks Moses 2005-08-29, 7:01 pm |
| Bart Vandewoestyne wrote:
> On 2005-08-23, Bart Vandewoestyne <MyFirstName.MyLastName@telenet.be> wrote:
> I'll leave it to the reader to decide which version is preferable. What I
> have learned from this challenge is the following:
>
> * First of all, look for a fast theoretical algorithm :-)
> * The loop-unrolling-alike technique as it is used by Brooks can speedup
> dramatically
> * I wonder if Fermat's method suffers from overflow problems since we are
> testing up until huge(n)-1 and this algorithm requires 8-byte integers...
> that would maybe explain the speed of the method, but apparently, it finds
> as much primes as the other ones so it seems correct? I haven't looked at
> the value of the primes yet, but this could be something for further
> investigation... :-)
>
> If I had to choose, I think I would go for Elliot's method because it doesn't
> require 8-byte integers and it is definitely more readable than Brooks's method ;-)
Another thing worth noting here: my loop-unrolling technique was
specifically aimed at the primes-below-10000000 version of the problem.
It performs well there, taking about 30% of the time of the next
comparable method. However, when the contest parameters are changed
slightly, to use numbers just below huge(n), it does only a small
fraction better than the method it's based on -- whereas Fermat's
method, which _is_ optimized for very large numbers, does startlingly
well, taking only a little more time than it required for small numbers,
whereas everything else takes more than an order of magnitude longer.
I think the lesson _I'd_ draw from this is that, when one's looking at
the fine details of optimizing routines like this, a problem
specification such as "calculate individual primes quickly" is not
really adequate. The "right answer" can vary dramatically depending on
the finer details.
In particular, consider this: I added more than a thousand lines of code
in order to optimize for what I saw as the detailed goal. Then the
details of the goal shifted (or, more accurately, my perception of the
goal was wrong!), and all of that optimization and complication really
didn't do any good. Here, it was just fun and games. In a real
programming situation, it would be a rather more costly mistake.
....
Meanwhile, were I to choose a general-purpose algorithm for raw speed
over the whole range, I'd probably use about the first part of mine, and
then tack Fermat's method on the end rather than Elliot's method (which
it currently uses) -- thereby having some measure of optimization at
both ends. But there are few cases where I'd be choosing a method for
raw speed; most times I'd want readability and quickness of programming,
and I agree that Elliot's method is very good for that.
- Brooks
--
The "bmoses-nospam" address is valid; no unmunging needed.
| |
| e p chandler 2005-08-29, 7:01 pm |
| Bart Vandewoestyne wrote:
> I'll leave it to the reader to decide which version is preferable. What I
> have learned from this challenge is the following:
>
> * First of all, look for a fast theoretical algorithm :-)
> * The loop-unrolling-alike technique as it is used by Brooks can speedup
> dramatically
> * I wonder if Fermat's method suffers from overflow problems since we are
> testing up until huge(n)-1 and this algorithm requires 8-byte integers...
> that would maybe explain the speed of the method, but apparently, it finds
> as much primes as the other ones so it seems correct? I haven't looked at
> the value of the primes yet, but this could be something for further
> investigation... :-)
>
Note that the value of huge depends on the *kind* of its argument. All
of Bart's numbers are of selected_int_kind(9), which I expect to be 4
byte integers. Hence huge(n) would be 2**31 - 1. The only place where I
think there is a danger of overflow in the Fermat routine is in
calc_power, where numbers of selected_int_kind(14) are used. Since n is
not an 8 byte integer, I do believe that overflow will not take place.
Thus I believe that the speed up seen was in fact genuine for numbers
close to huge(n). In addition, the numbers less than or equal to 2**31
- 1 are within the range where {2,3,5,7} (or subsets thereof) are the
proper test. (The single exceptional value cited is > 2**31 - 1.)
Overall I am not quite comfortable with the Fermat method. It clearly
only works for limited values of n. I do not know whether the tests for
primality cited are empirical or the result of mathematical proof.
------
Elliot
|
|
|
|
|