Update Mlucas.c, mi64.c for correct modular division (small scalars) #23
Conversation
Small fix, for PRP bases other than 3, to ensure the modular division by the square of the base at the end of a PRP test is not subject to a rounding error (bases 7, 27, 28, and 29 are susceptible among the possible PRP bases less than 32).
Make a tiny addition to `fquo` in the second limb of the modular division, in the case where a rounding error from the division results in the quotient lying between 0.999...999 (14 nines) and 1.0.
Tweaking the previous assertion for bases other than 3 to apply only to Suyama tests.
src/mi64.c
Outdated
There was a problem hiding this comment.
More discussion on this fix (whether this is an appropriate way to resolve the float to int conversion issue) is over at issue #18.
A few tweaks to the comments
Added warning to stderr if round-off error close to an integer encountered.
Extensive documentation on testing Fermat numbers with Mlucas (this pulls together a number of documented and undocumented information about Mlucas from various sources).
|
I’ve added in the Fermat documentation (and some small tweaks to the config script) from issue #1 as this is now ready to be added to the repository, whatever final disposition we reach on the status of whatever modification is desirable to |
|
|
||
| The `worktodo` format for the Suyama test is: | ||
| ``` | ||
| PRP=N/A,1,2,<2^m>,+1,99,0,3,5,"known_factor_1[,known_factor_2,...]" |
There was a problem hiding this comment.
It might help to explain more of these parameters, especially if we now support other PRP bases:
PRP=N/A,1,2,<2^m>,1[,how_far_factored,tests_saved[,base,residue_type[,"known_factors"]]]
The plus in the +1 is optional and assumed.
There was a problem hiding this comment.
I’m a bit on the fence about this; I don’t want to raise the issue of different bases in this context since the worktodo format for the Pépin test does not support bases other than 3; and choosing a different base for composite numbers always involves obtaining a completely different set of final residues.
(Also I like the plus in +1 to clearly indicate what it’s doing following the exponent.)
There was a problem hiding this comment.
(Yes, I like the +1 as well, I just thought I should note that it is not strictly required by the format.)
There was a problem hiding this comment.
Also: known_factors are mandatory for running a Suyama test, so the optional brackets disappear, and presently only residue_type 5 happens to be supported. Since you typically don’t want Mlucas to spontaneously decide to run P–1 prior to the PRP test, this sort of determines all of the variables preceding the factors:
PRP=N/A,1,2,<2^m>,1,how_far_factored,tests_saved,base,residue_type,"known_factors"
There was a problem hiding this comment.
Some users may want to run a PRP test on one of the (three) supported Fermat numbers without known factors, in which case maybe some of the parameters could potenchally be useful.
There was a problem hiding this comment.
Oh, interesting, thanks for the clarification.
Just a random aside from reading #4, since you write a lot of complex math formulas, GitHub supports LaTeX syntax, if you have used that before. See preda/gpuowl#278 for some examples.
There was a problem hiding this comment.
The last time I seriously used
I’m investigating how the Pépin test could be adapted for other bases (not all possible bases give a determinative result), so have been running PRP tests with the cofactor test following on; the main problem is the lack of specification for the base in the worktodo format. I coded other bases into my customised version of Gary’s cofact utility and it does appear the Mlucas code works correctly for Fermat modular squaring and cofactor testing in other bases, as I obtained identical Suyama results for base 5 tests of
There was a problem hiding this comment.
Interesting, thanks again for testing all of this. Regarding support for other bases, probably the easiest solution would be to update the existing PRP=/PRPDC= formats to fully support Fermat numbers, but as you noted above, this could still require a lot of changes.
I had not heard of Gary’s cofact utility, but we could always create a repository in this GitHub organization to store it, along with your customized version, if that would be helpful.
There was a problem hiding this comment.
I am amenable to adding a repository for GitHub but I would like Gary’s approval before doing so – and he was responsible for all versions before 0.6, which is the earliest I have a copy of. It’s a relatively simple Makefile and single C source file, but it would be nice to have all of the previous versions for populating blame with useful stuff; my fork splits from 0.6 and is now at 0.74c (the c stands for “Cathy”) while Gary is somewhere beyond 0.8.1.
There was a problem hiding this comment.
I can send Gary a message and ask if he would be OK moving the development of cofact to GitHub. We could then host your fork either in a separate repository or in a separate branch in the same repository, which ever you preferred. Unlike with Mlucas which is currently community maintained, you and Gary would have complete control over the cofact repository.
Even if Gary does not want to use GitHub, we could always host a mirror, as I do for the other Prime Search related programs. However, in this case I do not believe it could automatically be updated, as the downloads are only on the forum and not a website that could be scrapped, so it would be better to move the development to GitHub if possible.
Link within document as per suggestion Co-authored-by: Teal Dulcet <[email protected]>
A few tweaks as suggested by Teal.
One small amendment mprime actually does support Pépin tests (though it is handled somewhat awkwardly as it does not generate a certificate proof of the computation; but at least you do get a Res64). Added a link to docs/Fermat-testing.md
Explicitly replacing K with kilobytes (FFT length)
Update to Gerbicz’s nsquares integrity check, which assumes the PRP base is always 3.
Fix an ASSERT statement I mangled in the previous commit (that I mistakenly thought behaved like printf)
Successfully tested F13 with the 2K FFT, 8 8 16 radset
Eliminated 1K FFT, as I am now convinced 2K is the smallest possible workable FFT, which can support F13, F14, and F15. However, 2K requires fiddly tweaking of the Mlucas -f command that makes it a “one-off”.
Updated to add 2K FFT
if statement for F15 / 2K FFT
Consistency with assert around line 4100
Error in 2K FFT documentation
Tweak to last paragraph
Another 2K FFT tweak
|
I made a simplification to the For future reference, when you are just updating the documentation or code comments and do not need to run the CI to test the changes, you can add |
Substitute printing PRP_base in ASSERT with actual value by using cbuf.
[skip ci] Tweak to help.txt info
The if condition needed to be inverted, we want to possibly invoke an ASSERT if we do not have equality.
[skip ci] Adding Teal's suggestions.
Length of the residue in the arrtmp array is j, not i, so the SH 2^35-1 and 2^36-1 residues are not being calculated from the full residue.
| @@ -2432,7 +2433,7 @@ with the default #threads = 1 and affinity set to logical core 0, unless user ov | |||
| // by base and check that remainder 0. Note that we do want the quotient now, as that is our reside/base: | |||
| mi64_div(c_uint64_ptr, &itmp64, j+1,1, arrtmp,&rmodb); ASSERT(HERE, rmodb == 0ull,"After short-div, R != 0 (mod B)"); | |||
| // And recompute the S-H residues: | |||
| res_SH(arrtmp,i,&Res64,&Res35m1,&Res36m1); | |||
| res_SH(arrtmp,j,&Res64,&Res35m1,&Res36m1); | |||
There was a problem hiding this comment.
I think this value should be j, judging by the data I’ve looked at j usually is equal to i+1, so the calculation of the mod 2^35-1 and mod 2^36-1 residues is not encompassing the entire length of the arrtmp array. This partially resolves #27.
We refresh the 2^35-1 and 2^36-1 residues for Mersenne exponents (these values were not passed in by the Suyama_CF_PRP call).
| @@ -3310,6 +3311,8 @@ uint32 Suyama_CF_PRP(uint64 p, uint64*Res64, uint32 nfac, double a[], double b[] | |||
| /*A*/ ierr = func_mod_square(a, (int*)ci, n, ilo,ihi, 0ull, p, scrnFlag, tdiff, TRUE, 0x0); | |||
| convert_res_FP_bytewise(a, (uint8*)ci, n, p, Res64, &Res35m1, &Res36m1); // Overwrite passed-in Pepin-Res64 with Fermat-PRP one | |||
| snprintf_nowarn(cbuf,STR_MAX_LEN,"MaxErr = %10.9f\n",MME); mlucas_fprint(cbuf,1); | |||
| } else if (MODULUS_TYPE == MODULUS_TYPE_MERSENNE) { // Mersenne PRP-CF doesn't have the Res35m1 or Res36m1 values passed in, | |||
| res_SH(ci,n,&itmp64,&Res35m1,&Res36m1); // so we refresh these; see https://github.com/primesearch/Mlucas/issues/27 | |||
There was a problem hiding this comment.
This resolves the remaining part of #27, as the Suyama (A) value should print with correct Res35m1 and Res36m1 values. (Res64 is already correct so we don’t “re_update” it by referencing itmp64 instead.)
Assert condition covered by if loop.
[skip ci] UTC not needed in JSON string; see https://www.mersenneforum.org/showthread.php?p=653672#post653672
Just note that you would need to check the "Approve" radio button when writing the review to allow @xanthe-cat to merge this. Each PR needs at least one approving review. |
I would if I could see it :-) Sorry, I'm a newbie. EDIT: finally found it, sorry again... |
|
This is in a section of the code I am not familiar with. So I will trust Cathy's changes! |
Small fix, for PRP bases other than 3, to ensure the modular division by the square of the base at the end of a PRP test is not subject to a rounding error (bases 7, 14, 27, 28, and 29 are susceptible among the possible PRP bases less than 32). The rounding error for these bases appears in the
mi64_div_by_scalar64function ofmi64.c. We also modify theMlucas.ccode to permit PRP bases other than 3 by disabling theASSERTat line 698 (which would apply to both Mersenne and Fermat exponents) and instituting a similar check later, but only for Fermat numbers (since Mlucas always uses base 3 for the Pépin test).This fixes #18, fixes #27, as well as committing the Fermat testing documentation to the repository which fixes #1.