Re: [flint-devel] Comparisons between z_isprime and a Baillie-PSW test.

[Bill H. <goo...@go...>]

Sorry, I mean the numbers can exceed the size of a double, not an unsigned long.

Bill.


2008/8/20 William Hart <ha...@ya...>:
> Hi Peter,
>
> This is very interesting. There is a fairly dramatic
> jump at 1122004669633. This is the point where 5
> pseudoprime tests are used instead of 4, and not 10^16
> where Miller Rabin is used.
>
> The reason for the jump is that SPRP64 has to be used
> instead of SPRP since the numbers can exceed the size
> of an unsigned long.
>
> The Lucas test seems to be very good, and given that a
> counterexample to the conjecture that it plus a 2-SPRP
> test is a primality test is worth money, we may as
> well use it in FLINT for our probably prime testing up
> to 2^64.
>
> We'll use Pocklington-Lehmer for z_isprime after
> 10^16, I think. I'll be very interested to see what
> the time comparison is like for that, not that it
> matters, because we currently don't have any test in
> FLINT which is guaranteed to separate composites from
> primes after 10^16.
>
> Bill.
>
> --- Peter Shrimpton <ps...@gm...> wrote:
>
>> Hello
>>
>> Attached are spreadsheets (one excel the other
>> openoffice) showing the
>> relative speeds of z_isprime and a Baillie-PSW test.
>> In contradiction to my
>> earlier conversation with Bill (I was getting the
>> colours confused) the
>> Baillie-PSW test does perform a lot better then then
>> z_isprime as soon as
>> z_isprime has to do 5 pseudprime tests. However
>> whist testing this I found
>> that there is a bug in my lucas test code which
>> makes it return a handfull
>> of composites as prime. It only seems to affect
>> about 10 or so numbers which
>> are all 'quite small' (about 8 digits long).
>>
>> Also I have managed to get mongomery powering
>> working for unsigned longs
>> (see attached). I realised that as you only need the
>> first few bits for some
>> of the multiplications and so the whole opperation
>> can fit into 2 unsigned
>> longs. Currently it is slightly slower than z_powmod
>> in most places although
>> I can occasionally get improvements. I think It
>> depends on the relative
>> sizes of the parameters and/or how many 1's or 0's
>> they have but as there
>> are 3 of them I have no idea how to test it
>> completely to see which cases it
>> is faster.
>>
>> Peter
>> > #include <gmp.h>
>> #include "flint.h"
>> #include "long_extras.h"
>> #include <stdio.h>
>> #include "longlong.h"
>>
>> unsigned long mprod_precomp(unsigned long x,
>> unsigned long y, unsigned long n, unsigned long
>> nprime, unsigned long r1, unsigned long bits)
>> {
>>       unsigned long z,high,low,xyh,xyl;
>>
>>       umul_ppmm(xyh,xyl,x,y);
>>       z=(xyl)&r1;
>>       z=(z*(nprime))&r1;
>>       umul_ppmm(high,low,z,n);
>>       add_ssaaaa(high,low,xyh,xyl,high,low);
>>       z=(low>>bits)+(high<<(FLINT_BITS-bits));
>>
>>
>>       if (z>n){
>>               return z-n;
>>       }
>>       return z;
>> }
>>
>>
>> unsigned long mpow(unsigned long const x, unsigned
>> long const y, unsigned long const n)
>> {
>>       int bits=FLINT_BIT_COUNT(n);
>>       unsigned long r,xhat,phat,power,nprime;
>>       pre_inv2_t ninv;
>>       r=(1L<<bits);
>>       power=1L<<FLINT_BIT_COUNT(y);
>>       ninv=z_precompute_inverse(n);
>>       nprime=r-z_invert(n,r);
>>       phat=z_mod2_precomp(r,n,ninv);
>>       xhat=z_mulmod2_precomp(x,phat,n,ninv);
>>       r=r-1;
>>
>>       while (power)
>>       {
>>               phat=mprod_precomp(phat,phat,n,nprime,r,bits);
>>               if ((y&power)>0){
>>                       phat=mprod_precomp(phat,xhat,n,nprime,r,bits);
>>               }
>>               power=power>>1;
>>       }
>>       return mprod_precomp(phat,1L,n,nprime,r,bits);
>>
>> }
>>
>> int main()
>> {
>>       unsigned long x,y,n;
>>       int i,lengthx,lengthy,lengthn,minlen,maxlen;
>>       minlen=60;
>>       maxlen=63;
>>       //x=2L;
>>
>>       for (lengthx=minlen; lengthx<maxlen; lengthx++)
>>       {
>>               for (lengthy=minlen; lengthy<maxlen; lengthy++)
>>               {
>>                       for (lengthn=minlen; lengthn<maxlen; lengthn++)
>>                       {
>>                               //printf("%d^%d mod
>> %d\n",lengthx,lengthy,lengthn);
>>                               for (i=0; i<100000; i++)
>>                               {
>>                                       x=z_randbits(lengthx)+1;
>>                                       //y=z_randbits(lengthy)+1;
>>                                       n=2*z_randbits(lengthn-1)+1;
>>                                       y=n-1;
>>                                       //mpow(x,y,n);
>>                                       z_powmod2(x,y,n);
>>                               }
>>                       }
>>               }
>>       }
>> }
>>
>> >
> -------------------------------------------------------------------------
>> This SF.Net email is sponsored by the Moblin Your
>> Move Developer's challenge
>> Build the coolest Linux based applications with
>> Moblin SDK & win great prizes
>> Grand prize is a trip for two to an Open Source
>> event anywhere in the world
>>
> http://moblin-contest.org/redirect.php?banner_id=100&url=/>
> _______________________________________________
>> flint-devel mailing list
>> fli...@li...
>>
> https://lists.sourceforge.net/lists/listinfo/flint-devel
>>
>
>
>
>
>
> -------------------------------------------------------------------------
> This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
> Build the coolest Linux based applications with Moblin SDK & win great prizes
> Grand prize is a trip for two to an Open Source event anywhere in the world
> http://moblin-contest.org/redirect.php?banner_id=100&url=/
> _______________________________________________
> flint-devel mailing list
> fli...@li...
> https://lists.sourceforge.net/lists/listinfo/flint-devel
>