|
From: Peter S. <ps...@gm...> - 2008-07-25 09:50:54
|
Hello Attached is the code that I have written to be attached to long_extras along with the test code. Everything is passing the tests and I am confident that everything is working. I am fairly certain that some things can be optimized more which I will look into. Also could I ask for gcc to be upgraded on the server so that I can begin programing in OpenMP. Peter |
|
From: Bill H. <goo...@go...> - 2008-07-25 14:05:40
|
Excellent work!! I'll try to integrate these functions and test code this afternoon, and I'll try to put the right gcc on the server tomorrow. Bill. 2008/7/25 Peter Shrimpton <ps...@gm...>: > Hello > > Attached is the code that I have written to be attached to long_extras along > with the test code. Everything is passing the tests and I am confident that > everything is working. I am fairly certain that some things can be optimized > more which I will look into. Also could I ask for gcc to be upgraded on the > server so that I can begin programing in OpenMP. > > Peter > > ------------------------------------------------------------------------- > 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 > > |
|
From: Bill H. <goo...@go...> - 2008-08-08 15:12:24
|
Peter, I finally merged your Jacobi symbol code, but at the moment it hangs when trying to compute (282612/13). Can you take a look at it. I made a few changes. It seems that GMP returns 1 for (0/1), so I made your code return 1 when a is zero. I also noticed you were precomputing an inverse them using z_mod2_precomp. But this is slower than just doing a % b, since a precomputed inverse is only faster when you are computing more than one thing with the same b. It's possible some of the changes I made affected your code somehow. I also changed the test code to be more in line with the other tests in long_extras. Also, I noticed that you had something about GMP not handling negative numbers, but it actually does. If you use _si instead of _ui it will assume the number is signed. I didn't actually see any restriction on y being positive in the GMP documentation for mpz_jacobi. Does your code also handle negative y? I assumed for now that it doesn't. Bill. 2008/7/25 Peter Shrimpton <ps...@gm...>: > Hello > > Attached is the code that I have written to be attached to long_extras along > with the test code. Everything is passing the tests and I am confident that > everything is working. I am fairly certain that some things can be optimized > more which I will look into. Also could I ask for gcc to be upgraded on the > server so that I can begin programing in OpenMP. > > Peter > > ------------------------------------------------------------------------- > 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 > > |
|
From: Bill H. <goo...@go...> - 2008-08-08 15:25:42
|
Peter, I also noticed that in your code you multiply b*b*exp without regard for how large the result is, yet it still seems to work fine. Can you write some comments which explain why this still works, even if the product turns out to be longer than 64 bits. At first I thought the code would just fail, but it seems to actually work, even when x and y are very large. Bill. 2008/8/8 Bill Hart <goo...@go...>: > Peter, > > I finally merged your Jacobi symbol code, but at the moment it hangs > when trying to compute (282612/13). Can you take a look at it. > > I made a few changes. It seems that GMP returns 1 for (0/1), so I made > your code return 1 when a is zero. > > I also noticed you were precomputing an inverse them using > z_mod2_precomp. But this is slower than just doing a % b, since a > precomputed inverse is only faster when you are computing more than > one thing with the same b. > > It's possible some of the changes I made affected your code somehow. I > also changed the test code to be more in line with the other tests in > long_extras. > > Also, I noticed that you had something about GMP not handling negative > numbers, but it actually does. If you use _si instead of _ui it will > assume the number is signed. > > I didn't actually see any restriction on y being positive in the GMP > documentation for mpz_jacobi. Does your code also handle negative y? I > assumed for now that it doesn't. > > Bill. > > 2008/7/25 Peter Shrimpton <ps...@gm...>: >> Hello >> >> Attached is the code that I have written to be attached to long_extras along >> with the test code. Everything is passing the tests and I am confident that >> everything is working. I am fairly certain that some things can be optimized >> more which I will look into. Also could I ask for gcc to be upgraded on the >> server so that I can begin programing in OpenMP. >> >> Peter >> >> ------------------------------------------------------------------------- >> 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 >> >> > |
|
From: Peter S. <ps...@gm...> - 2008-08-08 16:33:38
|
Hey
I think there is a bug in the test code. I tried wrighting a seperate
program which just calculted (282612/13) and it worked fine. When I ran the
long_extras-test code it did hang but I think it was hanging in the loop:
do {
m = 2*z_randbits(bits2) + 1;
} while (z_gcd(a,m) > 1);
because I think a was 0. I replaced
a = z_randbits(bits1);
with
a = z_randbits(bits1)+1;
and that seamed to work for a bit. However it kept generating m's that were
negative. I looked at a few sources and all of them only define the Jacobi
symbol for m a positive odd integer. I think this is because the law of
quadratic reciprocity:
(a/m)(m/a)=(-1)^((m-1)(n-1)/4)
would break down or require complex numbers. So I replaced the line
ulong bits2 = z_randint(FLINT_BITS);
with
ulong bits2 = z_randint(FLINT_BITS-1);
and now it passes the test fine.
The reason why the code still works despite possible overflows is that we
only consider those statements mod 2. i.e. we only look at the last bit and
check to see if it is 1 or 0. Changing the lines:
if (((exp*(b*b-1))/8)%2 == 1)
and
if (((a-1)*(b-1)/4)%2 == 1)
with:
if (((exp*(b*b-1))/8)&1 == 1)
and
if (((a-1)*(b-1)/4)&1 == 1)
does exactly the same and is slightly clearer as to why it works.
In other news:
I have fixed the bug that caused the seg-fault in my code for doing multiple
Lucas tests at once (attached). However I ran some tests to see how quick it
was and it was between 2-10 times slower then just doing them individually.
The relative speed seamed to depend mainly on the size of N=n_1*n_2*n_3...
we can reduce this by decreasing the number of bits it preprocesses however
my code randomly stops working when the number of bits not preprocessed
reaches 16 and I have no idea why.
Some slightly better news is that I found a source today that says the $620
problem has been reduced slightly and we can replace the strong pseudoprime
test with an ordinary psudoprime test. i.e. we only need to check that
2^(n-1)=1 mod n and do a Lucas test. I wonder if it will be possible to
speed this up using the partial factorization we will have of n-1 from the
sieving? Or maybe do multiple tests at once?
Peter
2008/8/8 Bill Hart <goo...@go...>
> Peter,
>
> I also noticed that in your code you multiply b*b*exp without regard
> for how large the result is, yet it still seems to work fine. Can you
> write some comments which explain why this still works, even if the
> product turns out to be longer than 64 bits.
>
> At first I thought the code would just fail, but it seems to actually
> work, even when x and y are very large.
>
> Bill.
>
> 2008/8/8 Bill Hart <goo...@go...>:
> > Peter,
> >
> > I finally merged your Jacobi symbol code, but at the moment it hangs
> > when trying to compute (282612/13). Can you take a look at it.
> >
> > I made a few changes. It seems that GMP returns 1 for (0/1), so I made
> > your code return 1 when a is zero.
> >
> > I also noticed you were precomputing an inverse them using
> > z_mod2_precomp. But this is slower than just doing a % b, since a
> > precomputed inverse is only faster when you are computing more than
> > one thing with the same b.
> >
> > It's possible some of the changes I made affected your code somehow. I
> > also changed the test code to be more in line with the other tests in
> > long_extras.
> >
> > Also, I noticed that you had something about GMP not handling negative
> > numbers, but it actually does. If you use _si instead of _ui it will
> > assume the number is signed.
> >
> > I didn't actually see any restriction on y being positive in the GMP
> > documentation for mpz_jacobi. Does your code also handle negative y? I
> > assumed for now that it doesn't.
> >
> > Bill.
> >
> > 2008/7/25 Peter Shrimpton <ps...@gm...>:
> >> Hello
> >>
> >> Attached is the code that I have written to be attached to long_extras
> along
> >> with the test code. Everything is passing the tests and I am confident
> that
> >> everything is working. I am fairly certain that some things can be
> optimized
> >> more which I will look into. Also could I ask for gcc to be upgraded on
> the
> >> server so that I can begin programing in OpenMP.
> >>
> >> Peter
> >>
> >>
> -------------------------------------------------------------------------
> >> 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
>
|
|
From: Bill H. <goo...@go...> - 2008-08-08 16:46:14
|
OK thanks for finding that bug. Sorry it turned out not to be in your code.
That's good news about finding the segfault in your other code. I'm
sure we'll be able to speed it up, so don't be too concerned about
that.
Bill.
2008/8/8 Peter Shrimpton <ps...@gm...>:
> Hey
>
> I think there is a bug in the test code. I tried wrighting a seperate
> program which just calculted (282612/13) and it worked fine. When I ran the
> long_extras-test code it did hang but I think it was hanging in the loop:
>
> do {
> m = 2*z_randbits(bits2) + 1;
> } while (z_gcd(a,m) > 1);
>
> because I think a was 0. I replaced
>
> a = z_randbits(bits1);
>
> with
>
> a = z_randbits(bits1)+1;
>
> and that seamed to work for a bit. However it kept generating m's that were
> negative. I looked at a few sources and all of them only define the Jacobi
> symbol for m a positive odd integer. I think this is because the law of
> quadratic reciprocity:
>
> (a/m)(m/a)=(-1)^((m-1)(n-1)/4)
>
> would break down or require complex numbers. So I replaced the line
>
> ulong bits2 = z_randint(FLINT_BITS);
>
> with
>
> ulong bits2 = z_randint(FLINT_BITS-1);
>
> and now it passes the test fine.
>
> The reason why the code still works despite possible overflows is that we
> only consider those statements mod 2. i.e. we only look at the last bit and
> check to see if it is 1 or 0. Changing the lines:
>
> if (((exp*(b*b-1))/8)%2 == 1)
> and
> if (((a-1)*(b-1)/4)%2 == 1)
>
> with:
>
> if (((exp*(b*b-1))/8)&1 == 1)
> and
> if (((a-1)*(b-1)/4)&1 == 1)
>
> does exactly the same and is slightly clearer as to why it works.
>
> In other news:
>
> I have fixed the bug that caused the seg-fault in my code for doing multiple
> Lucas tests at once (attached). However I ran some tests to see how quick it
> was and it was between 2-10 times slower then just doing them individually.
> The relative speed seamed to depend mainly on the size of N=n_1*n_2*n_3...
> we can reduce this by decreasing the number of bits it preprocesses however
> my code randomly stops working when the number of bits not preprocessed
> reaches 16 and I have no idea why.
>
> Some slightly better news is that I found a source today that says the $620
> problem has been reduced slightly and we can replace the strong pseudoprime
> test with an ordinary psudoprime test. i.e. we only need to check that
> 2^(n-1)=1 mod n and do a Lucas test. I wonder if it will be possible to
> speed this up using the partial factorization we will have of n-1 from the
> sieving? Or maybe do multiple tests at once?
>
> Peter
>
> 2008/8/8 Bill Hart <goo...@go...>
>>
>> Peter,
>>
>> I also noticed that in your code you multiply b*b*exp without regard
>> for how large the result is, yet it still seems to work fine. Can you
>> write some comments which explain why this still works, even if the
>> product turns out to be longer than 64 bits.
>>
>> At first I thought the code would just fail, but it seems to actually
>> work, even when x and y are very large.
>>
>> Bill.
>>
>> 2008/8/8 Bill Hart <goo...@go...>:
>> > Peter,
>> >
>> > I finally merged your Jacobi symbol code, but at the moment it hangs
>> > when trying to compute (282612/13). Can you take a look at it.
>> >
>> > I made a few changes. It seems that GMP returns 1 for (0/1), so I made
>> > your code return 1 when a is zero.
>> >
>> > I also noticed you were precomputing an inverse them using
>> > z_mod2_precomp. But this is slower than just doing a % b, since a
>> > precomputed inverse is only faster when you are computing more than
>> > one thing with the same b.
>> >
>> > It's possible some of the changes I made affected your code somehow. I
>> > also changed the test code to be more in line with the other tests in
>> > long_extras.
>> >
>> > Also, I noticed that you had something about GMP not handling negative
>> > numbers, but it actually does. If you use _si instead of _ui it will
>> > assume the number is signed.
>> >
>> > I didn't actually see any restriction on y being positive in the GMP
>> > documentation for mpz_jacobi. Does your code also handle negative y? I
>> > assumed for now that it doesn't.
>> >
>> > Bill.
>> >
>> > 2008/7/25 Peter Shrimpton <ps...@gm...>:
>> >> Hello
>> >>
>> >> Attached is the code that I have written to be attached to long_extras
>> >> along
>> >> with the test code. Everything is passing the tests and I am confident
>> >> that
>> >> everything is working. I am fairly certain that some things can be
>> >> optimized
>> >> more which I will look into. Also could I ask for gcc to be upgraded on
>> >> the
>> >> server so that I can begin programing in OpenMP.
>> >>
>> >> Peter
>> >>
>> >>
>> >> -------------------------------------------------------------------------
>> >> 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
>
>
> -------------------------------------------------------------------------
> 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
>
>
|
|
From: Bill H. <goo...@go...> - 2008-08-23 13:32:53
|
Hi Peter, I merged you Lucas pseudoprime tests and wrote test code. However I don't know what values of a and b are allowed in the first version of the test. I removed some minor bugs, and it also seems that P was being set to 1 in the second test and therefore not being used, so I removed it. The second version works beautifully!! Bill. 2008/7/25 Peter Shrimpton <ps...@gm...>: > Hello > > Attached is the code that I have written to be attached to long_extras along > with the test code. Everything is passing the tests and I am confident that > everything is working. I am fairly certain that some things can be optimized > more which I will look into. Also could I ask for gcc to be upgraded on the > server so that I can begin programing in OpenMP. > > Peter > > ------------------------------------------------------------------------- > 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 > > |
|
From: Peter S. <ps...@gm...> - 2008-08-25 12:18:48
|
Hello. Sorry for the late response. The first version should work for everything except for a=1,b=1 and a=-1,b=1 although I am not sure how 'good' it is for large values of a and b. (i.e. how many composites get through). I noticed that the "bug" I thought was in the code was actually me misreading the $620 problem. There appear to be two different versions: 1) a composite congruent to 3 or 7 mod 10 which passes both a 2-Fermat-pesudoprime test and a Lucas test with parameters a=1,b=-1 2) any composite which passes a 2-Strong-pesudoprime test and the second version of the Lucas test. I think the first is significantly quicker especially when a and b are fixed because we can avoid doing the gcd and issquare tests. Also I did some more tests on the Montgomery powering that I wrote with mixed results. I did random 2-Fermat-pesudoprime tests to numbers of various lengths and found that it was quicker.......but only when the numbers reached 60 bits in length. We are only dealing with numbers between 53 and 57 bits. I am currently working on optimizing the sieve so that it only looks at odd numbers but still factorizes the even numbers. I need to calculate the smallest j such that: beginning+2*j = 0 mod p I have tried using z_invert and modular arrithmethic but it seams to be slower then a simple while (begining+2*j mod p !=0) j++; and there is no guarantee that they will give the smallest j which is required. Is there a quicker way of doing it? Peter 2008/8/23 Bill Hart <goo...@go...> > Hi Peter, > > I merged you Lucas pseudoprime tests and wrote test code. However I > don't know what values of a and b are allowed in the first version of > the test. > > I removed some minor bugs, and it also seems that P was being set to 1 > in the second test and therefore not being used, so I removed it. > > The second version works beautifully!! > > Bill. > > 2008/7/25 Peter Shrimpton <ps...@gm...>: > > Hello > > > > Attached is the code that I have written to be attached to long_extras > along > > with the test code. Everything is passing the tests and I am confident > that > > everything is working. I am fairly certain that some things can be > optimized > > more which I will look into. Also could I ask for gcc to be upgraded on > the > > server so that I can begin programing in OpenMP. > > > > Peter > > > > ------------------------------------------------------------------------- > > 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 > |
|
From: William H. <ha...@ya...> - 2008-08-25 15:50:32
|
Can you take a look at the current code and test code, because when I use anything other than a = 1 and b = 1, it doesn't declare all primes to be pseudoprime. I did try to fix some minor issues with the original code, but maybe I screwed something up. As for your question: If I want to do a/2 mod p, then if a is even I just do a/2, if a is odd I just do (a+p)/2. No inverse mod p is necessary. Therefore all I need to do is: 1) let beg = beginning mod p 2) let s = p - beg (let s = 0 if beg = 0) 3) compute s/2 mod p. Bill. Bill. --- Peter Shrimpton <ps...@gm...> wrote: > Hello. Sorry for the late response. > > The first version should work for everything except > for a=1,b=1 and a=-1,b=1 > although I am not sure how 'good' it is for large > values of a and b. (i.e. > how many composites get through). I noticed that the > "bug" I thought was in > the code was actually me misreading the $620 > problem. There appear to be two > different versions: > > 1) a composite congruent to 3 or 7 mod 10 which > passes both a > 2-Fermat-pesudoprime test and a Lucas test with > parameters a=1,b=-1 > > 2) any composite which passes a 2-Strong-pesudoprime > test and the second > version of the Lucas test. > > I think the first is significantly quicker > especially when a and b are fixed > because we can avoid doing the gcd and issquare > tests. > > Also I did some more tests on the Montgomery > powering that I wrote with > mixed results. I did random 2-Fermat-pesudoprime > tests to numbers of various > lengths and found that it was quicker.......but only > when the numbers > reached 60 bits in length. We are only dealing with > numbers between 53 and > 57 bits. > > I am currently working on optimizing the sieve so > that it only looks at odd > numbers but still factorizes the even numbers. I > need to calculate the > smallest j such that: > > beginning+2*j = 0 mod p > > I have tried using z_invert and modular arrithmethic > but it seams to be > slower then a simple while (begining+2*j mod p !=0) > j++; and there is no > guarantee that they will give the smallest j which > is required. Is there a > quicker way of doing it? > > Peter > > 2008/8/23 Bill Hart <goo...@go...> > > > Hi Peter, > > > > I merged you Lucas pseudoprime tests and wrote > test code. However I > > don't know what values of a and b are allowed in > the first version of > > the test. > > > > I removed some minor bugs, and it also seems that > P was being set to 1 > > in the second test and therefore not being used, > so I removed it. > > > > The second version works beautifully!! > > > > Bill. > > > > 2008/7/25 Peter Shrimpton <ps...@gm...>: > > > Hello > > > > > > Attached is the code that I have written to be > attached to long_extras > > along > > > with the test code. Everything is passing the > tests and I am confident > > that > > > everything is working. I am fairly certain that > some things can be > > optimized > > > more which I will look into. Also could I ask > for gcc to be upgraded on > > the > > > server so that I can begin programing in OpenMP. > > > > > > Peter > > > > > > > ------------------------------------------------------------------------- > > > 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 > > > > ------------------------------------------------------------------------- > 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 > |