Jacobi symbol computation

[nemothenoone]

Resolves #170 and NilFoundation#30. Depends on #161 (not really with the functionality, but it serves the same purpose as modular adaptor does) and #143.

[nemothenoone]

I haven't tried to investigate why, but running tommath through your jacobi code is slower still:

BM_jacobi<boost::multiprecision::cpp_int>/512         44.6 ms         43.5 ms           14
BM_jacobi<boost::multiprecision::cpp_int>/1024         108 ms          107 ms            7
BM_jacobi<boost::multiprecision::cpp_int>/2048         260 ms          260 ms            3
BM_jacobi<boost::multiprecision::cpp_int>/4096         543 ms          547 ms            1
BM_jacobi<boost::multiprecision::cpp_int>/8192        1680 ms         1641 ms            1
BM_jacobi<boost::multiprecision::cpp_int>/16384       4859 ms         4781 ms            1
BM_jacobi<boost::multiprecision::cpp_int>/32768      17869 ms        17688 ms            1
BM_jacobi<boost::multiprecision::cpp_int>_BigO       16.76 N^2       16.58 N^2
BM_jacobi<boost::multiprecision::cpp_int>_RMS            8 %             8 %
BM_jacobi<boost::multiprecision::mpz_int>/512         10.3 ms         10.2 ms           75
BM_jacobi<boost::multiprecision::mpz_int>/1024        19.0 ms         18.1 ms           45
BM_jacobi<boost::multiprecision::mpz_int>/2048        47.3 ms         46.9 ms           15
BM_jacobi<boost::multiprecision::mpz_int>/4096        82.6 ms         78.1 ms            7
BM_jacobi<boost::multiprecision::mpz_int>/8192         200 ms          198 ms            3
BM_jacobi<boost::multiprecision::mpz_int>/16384        687 ms          688 ms            1
BM_jacobi<boost::multiprecision::mpz_int>/32768       1452 ms         1422 ms            1
BM_jacobi<boost::multiprecision::mpz_int>_BigO     2910.28 NlgN    2860.64 NlgN
BM_jacobi<boost::multiprecision::mpz_int>_RMS           14 %            14 %
BM_jacobi<boost::multiprecision::tom_int>/512          508 ms          500 ms            2
BM_jacobi<boost::multiprecision::tom_int>/1024        1058 ms         1062 ms            1
BM_jacobi<boost::multiprecision::tom_int>/2048        2065 ms         2031 ms            1
BM_jacobi<boost::multiprecision::tom_int>/4096        3426 ms         3328 ms            1
BM_jacobi<boost::multiprecision::tom_int>/8192        5599 ms         5594 ms            1
BM_jacobi<boost::multiprecision::tom_int>/16384      14822 ms        14766 ms            1
BM_jacobi<boost::multiprecision::tom_int>/32768      45262 ms        45000 ms            1
BM_jacobi<boost::multiprecision::tom_int>_BigO       43.11 N^2       42.87 N^2
BM_jacobi<boost::multiprecision::tom_int>_RMS           20 %            20 %

Well, I've just introduced tommath-specific Jacobi symbol computation method usage, so the performance should improve.

[nemothenoone]

• Couldn't gmp_int use mpz_jacobi ?

Of course. Done.

• For personal preference I'd like to see some more tests, the big matrix of values at https://en.wikipedia.org/wiki/Jacobi_symbol would be a good start. Or use a "known good" (Wolfram Alpha?) to churn out a table of test values.

Agreed. The matrix-based testing is coming.

• As noted above, like it or not, Boost.Build is still our testing tool: the test will need to be added to the Jamfile (this also helps us check that CI passes OK).

My bad. Forgot to do that. Fixed. (And in modular adaptor PR as well).

• I believe Jacobi symbol computation is based on the Euclid GCD algorithm? Does the Lehmer GCD code help here (recently added via Better gcd: Lehmers algorithm #250)? It's way way faster than plain Euclid.

Nah, these are related algorithms, but they are not using each other. So we're stuck with this one.

• We also need documentation - if you really really don't want to mess with the doc quickbook sources, let me have some plain text or rst markup and I'll translate.

Agreed. The documentation should be here. I'ill do that myself in quickbook.

[jzmaddock]

I believe Jacobi symbol computation is based on the Euclid GCD algorithm? Does the Lehmer GCD code help here (recently added via #250)? It's way way faster than plain Euclid.

Nah, these are related algorithms, but they are not using each other. So we're stuck with this one.

See https://arxiv.org/abs/1907.07795

[nemothenoone]

See https://arxiv.org/abs/1907.07795

Oh. That makes things even better. I'ill reimplement the algorithm according to this paper.