# Enhanced ulong_extras.h/n_is_prime_pocklington #110

Merged
merged 11 commits into from Jan 27, 2015

## Conversation

Projects
None yet
4 participants

### kush789 commented Jan 25, 2015

 As mentioned in the todo, now n requires factorisation only up to n^1/3.

### kush789 added some commits Jan 18, 2015

``` Enhanced is_prime_pocklington. Now requires factorization of n only u… ```
`…pto n^1/3`
``` 731ded4 ```
``` Enhanced is_prime_pocklington ```
``` c907d4c ```
``` Enhanced ulong_extras/is_prime_pocklington ```
``` 9f7b6c5 ```
``` Cleaned up code. ```
``` 6d762cd ```
``` Cleaned up code. ```
``` 2d14b86 ```
``` Updated copyright header ```
``` d9f7852 ```
Owner

### wbhart commented Jan 25, 2015

Hi Kushagra,

thanks for this. Sorry for the delay in looking it over.

The test we're using is called the Brillhart-Lehmer-Selfridge (BLS) test.
It is only applicable if n^1/3 <= F < n^1/2.

• At the moment you are testing for F < floor(n^1/2) which is not quite the
same thing. You really want F <= n_sqrt(n), not F < n_sqrt(n), unless n is
a square. In the case n is a square we want to exclude F = sqrt(n) only.

You should change it to test if F <= n_sqrt(n). Then, since you compute F_F
anyway, it should be easy to add a test if n == F_F and return 0 if so to
rule out that one case that we want to exclude.

• You may as well reuse the computed value of F*F wherever you need it. You
can even declare it const inside the BLS block.
• The test also only works if F >= n^1/3. At present I don't see where we
are testing for that.
• It looks like the same code is duplicated twice: once if the cube root
test doesn't immediately return 0 and once if F > n^1/2. Both of these
sections of code could be combined into one, and use a goto to jump into it
from one of the locations where it is currently needed.

Once you make these changes and I've checked it, I'll be happy to sign off
on it and merge it.

It might also be worth updating the documentation in
ulong_extras/doc/ulong_extras.txt to say that n_is_prime_pocklington now
uses the BLS test.

Best Wishes,

Bill.

On 25 January 2015 at 04:56, Kushagra Singh notifications@github.com
wrote:

## As mentioned in the todo, now n requires factorisation only up to n^1/3.

You can view, comment on, or merge this pull request online at:

#110
Commit Summary

• Enhanced is_prime_pocklington. Now requires factorization of n only
upto n^1/3
• Enhanced is_prime_pocklington
• Enhanced ulong_extras/is_prime_pocklington
• Cleaned up code.
• Cleaned up code.

File Changes

Reply to this email directly or view it on GitHub
#110.

### danaj commented Jan 25, 2015

 A small request that you be more specific and point out exactly which theorem or corollary from the BLS75 paper is being used (there are over 30 in the paper). Ideally even quoting or paraphrasing it in the documentation. [Added:] Having the conditions clearly stated should help you as well as reviewers and later readers check that the code correctly tests all the conditions and hence results in a valid proof. I also wonder if the name should be changed, since theorem 5, that allows proof once factored past (N/2)^(1/3), isn't from Pocklington any more so seems misleading. Or leave is_prime_pocklington and add is_prime_bls75_t5 (or whatever it gets called). This is a bigger and longer term decision. Dana
Owner

### wbhart commented Jan 25, 2015

 The reference I have is the book by Pomerance and Crandall: "Primes: a Computational Perspective". It's Theorem 4.1.5 on pp.175. The proof is given in the book, so I think that is a pretty reasonable reference to give. The way the implementation is currently combined with the Pocklington test is better, I think, than having two separate functions. But it would be a good idea to document that our implementation of Pocklington makes use of this improvement. It actually only adds a very few lines of lines of code, and I think in the context it is used, it is better thought of as an improvement to that test, rather than a separate test. On 25 January 2015 at 22:09, Dana Jacobsen notifications@github.com wrote: A small request that you be more specific and point out exactly which theorem or corollary from the BLS75 paper is being used (there are over 30 in the paper). Ideally even quoting or paraphrasing it in the documentation. I also wonder if the name should be changed, since theorem 5, that allows proof once factored past (N/2)^(1/3), isn't from Pocklington any more so seems misleading. Or leave is_prime_pocklington and add is_prime_bls75_t5 (or whatever it gets called). This is a bigger and longer term decision. Dana — Reply to this email directly or view it on GitHub #110 (comment).

### danaj commented Jan 25, 2015

 That makes sense, and the book is certainly well known enough to use as a reference. It's a little less general than BLS75 T5 but not much different -- the thing that stands out most to me is requiring a single 'a' for eq 4.3, rather than a separate 'a' for each prime factor of F as done in BLS (I). That's tangential though. Moving from n^1/2 to n^1/3 while keeping most of the code intact is a very nice improvement.
``` cleaned up code, increased upper bound on F ```
``` d0890bd ```
Author

### kush789 commented Jan 26, 2015

 Hey, I have changed the inequality, so that it now checks for F <=sqrt(n). I have also added the related test for a perfect square you suggested. When we calculate F (or the co factor to be more precise) using n_factor_partial, the function makes sure that the product of the factors is greater than the limit (which I have provided as n^1/3). So i guess there is no need to check this constraint again. I have removed the duplicate code, the current version looks much cleaner. I have added the source of the theory behind this implementation of the pocklington test in the documentation. One thing which is troubling me is that the test fails for 5, 17, 19, 37, and 73 repeatedly when I try to do a make check. However when I ran a test loop till INT_MAX, it was giving me the exact same results, as is_prime, apart from a few undetermined verdicts. Thanks, Kushagra
Owner

### wbhart commented Jan 26, 2015

I think the cbrt function is dangerous to use here. If n is larger than 53
bits it will get truncated. The cube root will be computed well enough, but
only of the truncated value. Since the theorem is only guaranteed for F >
n^1/3 I think we need to add one to the computed cube root (if n >= 2^53)
for the theorem to be valid.

Also, I don't think the factor_partial function will ensure F is > n^1/3.
It will look for factors up to that limit. The value F is the product of
the factors it finds, and could easily be 1, even though limit was set to
n^1/3.

Basically, if F < n^1/3 you can't do anything. It should return
undetermined.

Best Wishes,

Bill.

On 26 January 2015 at 15:35, Kushagra Singh notifications@github.com
wrote:

Hey,

## I have removed the duplicate code, the current version looks much cleaner.

I have added the source of the theory behind this implementation of
the pocklington test in the documentation.

One thing which is troubling me is that the test fails for 5, 17, 19, 37,
and 73 repeatedly when I try to do a make check. However when I ran a test
loop till INT_MAX, it was giving me the exact same results, as is_prime,
apart from a few undetermined verdicts.

Thanks,
Kushagra

Reply to this email directly or view it on GitHub
#110 (comment).

``` Update is_prime_pocklington.c ```
``` fbd0641 ```
Author

### kush789 commented Jan 26, 2015

 Hey, I added the suggestion which you gave, regarding the limit if n exceeds 2^53. I also have added the test which returns undetermined if F < n^1/3. I was confused initially regarding this as the documentation says that n_factor_partial factors n till the product of the prime factors exceeds the limit. Thanks, Kushagra
Owner

### wbhart commented Jan 26, 2015

 Actually, I think you are right. n_factor_partial does factor until the product exceeds the limit. Sorry for misleading you there. It can fail, but produces an exception in that case (it's not been known to fail, though one day we ought to ensure it really can't). I guess what you had was correct in that case, so long as the limit is set correctly. Do you want to change that bit back? The add one thing was correct however. Bill. On 26 January 2015 at 18:39, Kushagra Singh notifications@github.com wrote: Hey, I added the suggestion which you gave, regarding the limit if n exceeds 2^53. I also have added the test which returns undetermined if F < n^1/3. I was confused initially regarding this as the documentation says that n_factor_partial factors n till the product of the prime factors exceeds the limit. Thanks, Kushagra — Reply to this email directly or view it on GitHub #110 (comment).
Author

### kush789 commented Jan 26, 2015

 Yeah, I'll do that and push it. Saves unnecessary comparisons. Thanks, Kushagra
``` Removed lower bound check on F ```
``` 6aa933c ```
Owner

### wbhart commented Jan 26, 2015

 It's actually fairly slow to force it to keep going until F > n^1/3 because it has to use some heavy duty factoring routines. It wouldn't be insane to allow it to only use light weight factoring routines and just fail if it doesn't make its n^1/3 quota. But let's not change it to do that. Most people will call the routine because they want primality proved, not because they want a really fast primality test (which we already have). On 26 January 2015 at 18:53, Kushagra Singh notifications@github.com wrote: Yeah, I'll do that and push it. Saves unnecessary comparisons. Thanks, Kushagra — Reply to this email directly or view it on GitHub #110 (comment).
Author

### kush789 commented Jan 26, 2015

 So should I make any further changes to this implementation? Kushagra
``` corrected function name ```
``` 28c2014 ```
Owner

### wbhart commented Jan 26, 2015

 I can't see any other problems with it. Could you perhaps add a note about why we add 1 when > 2^53. Just a mention of the fact that we are using doubles which can only store 53 bits of the mantissa when taking the square root. Do you still have test failures? If so, could you set n to one of the failing values in the test code, temporarily, and put some printf's in the code to try and figure out where the failure is occurring? On 26 January 2015 at 19:09, Kushagra Singh notifications@github.com wrote: So should I make any further changes to this implementation? Kushagra — Reply to this email directly or view it on GitHub #110 (comment).
``` Fixed code to remove build error ```
``` 027cc5c ```
Author

### kush789 commented Jan 27, 2015

 Hey, I have fixed the code. There was a minor bug, n_factor_partial was returning 1 as a factor also, which could not be used as it is not prime. It is not giving any test errors now. Kushagra
Owner

### wbhart commented Jan 27, 2015

 That's great! I'll try to merge this into my repository today. Thanks for your hard work! On 27 January 2015 at 06:00, Kushagra Singh notifications@github.com wrote: Hey, I have fixed the code. There was a minor bug, n_factor_partial was returning 1 as a factor also, which could not be used as it is not prime. It is not giving any test errors now. Kushagra — Reply to this email directly or view it on GitHub #110 (comment).

Owner

### wbhart commented Jan 27, 2015

 I've now merged this. But we should try to fix the bug in the factorisation code. It shouldn't be putting the factor 1 in the factor struct. If you have some time, you could track this down for us. Either way I will open a ticket for this. On 27 January 2015 at 10:47, Bill Hart goodwillhart@googlemail.com wrote: That's great! I'll try to merge this into my repository today. Thanks for your hard work! On 27 January 2015 at 06:00, Kushagra Singh notifications@github.com wrote: Hey, I have fixed the code. There was a minor bug, n_factor_partial was returning 1 as a factor also, which could not be used as it is not prime. It is not giving any test errors now. Kushagra — Reply to this email directly or view it on GitHub #110 (comment).
Owner

### wbhart commented Feb 16, 2015

 I just discovered a minor problem with this implementation. The cbrt function is only available in C99, but flint only uses C90. This leads to a compilation warning. It would be good if the call to cbrt could be replaced with pow(x, 1.0/3), which should be supported on all systems. On 27 January 2015 at 15:34, Bill Hart goodwillhart@googlemail.com wrote: I've now merged this. But we should try to fix the bug in the factorisation code. It shouldn't be putting the factor 1 in the factor struct. If you have some time, you could track this down for us. Either way I will open a ticket for this. On 27 January 2015 at 10:47, Bill Hart goodwillhart@googlemail.com wrote: That's great! I'll try to merge this into my repository today. Thanks for your hard work! On 27 January 2015 at 06:00, Kushagra Singh notifications@github.com wrote: Hey, I have fixed the code. There was a minor bug, n_factor_partial was returning 1 as a factor also, which could not be used as it is not prime. It is not giving any test errors now. Kushagra — Reply to this email directly or view it on GitHub #110 (comment).
Owner

### wbhart commented Feb 20, 2015

 I think there is still one more issue we need to address. 1.0/3 is represented in binary as 0.01010101010101010101... This means, I think, that when represented with 53 bits, the value will be slightly less than the actual value of 1/3. For that reason I think we need to unconditionally add 1 to limit, rather than just for > 2^52. That will accommodate for the fact that we aren't quite taking the cube root. You could of course go through every number x from 1 to 10^11 or something, and see if (pow(x, 1.0/3) + 1)^3 > x. If so then you would only need to add 1 to limit from 10^11 onwards. If there is an exception, limit only needs to have 1 added to it from that point onwards. Another concern is the pow function however. This might be implemented poorly in some places. If so, then perhaps it is best to add 1 unconditionally to limit.

### fredrik-johansson commented Feb 20, 2015

 pow() can legally return whatever it wants. I think we should have an n_nthroot function that guarantees computing floor(x^(1/n)), for example by computing a floating-point approximation and then incrementing or decrementing until it's correct. Is there a reason that the condition cannot be written as x^3 < y instead of x < y^(1/3), though?
Owner

### wbhart commented Feb 20, 2015

 The condition can be written x^3 < y. That seems equivalent to me. So i's not really important how good an approximation pow gives, only that the integer part is not off by more than one. On 20 February 2015 at 18:31, Fredrik Johansson notifications@github.com wrote: pow() can legally return whatever it wants. I think we should have an n_nthroot function that guarantees computing floor(x^(1/n)), for example by computing a floating-point approximation and then incrementing or decrementing until it's correct. Is there a reason that the condition cannot be written as x^3 < y instead of x < y^(1/3), though? — Reply to this email directly or view it on GitHub #110 (comment).
Author

### kush789 commented Feb 21, 2015

 Hi, Would it be wrong if we just tested whether limit^3 < n ? If it is less, then we can increment limit till (limit)^3 becomes greater then n. This would also take care of the condition on n-1 that we were checking earlier (whether n-1 fits in a double or not ), as we would be using just mp_limb_t types to check this, not doubles. A corner case would be for a value of n = 2^64. To ensure there is no overflow, in the loop where I am increasing limit if necessary, I am also checking that limit doesn't cross 2642246, as 2642246^3 is just greater than 2^64, which is the upper bound for n in this case. I am opening a new pull request regarding the same, kindly review the code and let me know if something looks off. Thanks Kushagra