Pari/GP CAS (Computer Algebra System)

[HermannSW]

Just stumbled over Pari/GP Computer Algebra System:
https://en.wikipedia.org/wiki/PARI/GP
https://pari.math.u-bordeaux.fr/

Developed since the early eighties in the University of Bordeaux, PARI is free software, available under the GPL.

Pari is a C library, gp its command line interface.
I found this nice (4 pages) reference card:
https://math.mit.edu/~brubaker/PARI/PARIrefcard.pdf

You can play with Pari/GP without installation in your browser:
https://pari.math.u-bordeaux.fr/gp.html


I am interested in fast multiple precision integer arithmetic and integer factorization with my RSA_numbers_factored project in Python+gmpy2+sympy sofar.

Pari comes with all that, and while it has weaknesses (string handling, no memoization), this statement is killer argument for Pari/GP for me:
https://oeis.org/wiki/PARI/GP#Comparison_to_other_CAS

In the domain of number theory, PARI/GP is a strong rival to well established all-purpose CAS as Maple and Mathematica, mainly due to its computational speed, where it typically outperforms both of these well known commercial CAS, and its free availability.

After

Code: Select all

sudo apt install pari-gp

it can be used on Raspberry (64bit) PiOS, with its command line processor "gp".

Simple factorization example in command line interpreter:

Code: Select all

pi@pi400-64:~/Documents $ gp -q
? print(factorint(2^ 256-1))
[3, 1; 5, 1; 17, 1; 257, 1; 641, 1; 65537, 1; 274177, 1; 6700417, 1; 67280421310721, 1; 59649589127497217, 1; 5704689200685129054721, 1]
? 

gp can execute scripts as well, or like this:

Code: Select all

pi@pi400-64:~ $ gp -q <(echo "print(factorint(2^256-1));quit")
[3, 1; 5, 1; 17, 1; 257, 1; 641, 1; 65537, 1; 274177, 1; 6700417, 1; 67280421310721, 1; 59649589127497217, 1; 5704689200685129054721, 1]
pi@pi400-64:~ $ 

This is "factorint()" doc from page 193 of 675 pages manual:
https://pari.math.u-bordeaux.fr/pub/par ... f#page=193

Pari_factorint.png

Pari_factorint.png (159.44 KiB) Viewed 2699 times

Smallest RSA number from RSA_numbers_factored repo is 59 decimal digits number RSA-59 (biggest numbers are RSA-617 and RSA-2048):

Code: Select all

pi@pi400-64:~/RSA_numbers_factored/python $ python -q
>>> from RSA_numbers_factored import RSA, digits
>>> R = RSA()
>>> l,n,p,q = R.get(59)[0:4]
>>> print(l,n,p,q)
59 71641520761751435455133616475667090434063332228247871795429 200429218120815554269743635437 357440504101388365610785389017
>>> l == digits(n) and n == p * q
True
>>>

I started Python sympy factorint 65 minutes ago on Pi400, and it has not completed yet:

Code: Select all

pi@pi400-64:~/RSA_numbers_factored/python $ python <(echo "from sympy import factorint; print(factorint(71641520761751435455133616475667090434063332228247871795429))")

This is a wow for Pari/GP for me — only 12 seconds for factoring RSA-59 with Pari/GP on my (1.8GHz default) Raspberry Pi400:

Code: Select all

pi@pi400-64:~ $ time gp -q <(echo "print(factorint(71641520761751435455133616475667090434063332228247871795429));quit")
[200429218120815554269743635437, 1; 357440504101388365610785389017, 1]

real	0m12.048s
user	0m11.971s
sys	0m0.030s
pi@pi400-64:~ $ 

Just to be clear, Pari/GP factorint is nice, but surely is no help to factor any RSA number above biggest factored sofar RSA number RSA-250. But it will surely be a cool addition to Python+gmpy2+sympy in my tool chest!

[HermannSW]

Python sympy "factorint()" completed RSA_59 factorization after hours of computation.
That does only demonstate that the algorithmic ideas of Pari "factorint()" are much better.

I currently think of transpiling RSA_numbers_factored.py to RSA_numbers_factored.gp manually, like I did to JavaScript/NodeJS:
https://github.com/Hermann-SW/RSA_numbe ... mplates.md

First I wanted syntax highlighting for Pari/GP, and my favorite editor vi (that I use since 1990) does provide that already:

Code: Select all

pi@pi400-64:~ $ head -5 /usr/share/vim/vim82/syntax/gp.vim 
" Vim syntax file
" Language:	gp (version 2.5)
" Maintainer:	Karim Belabas <[email protected]>
" Last change:	2012 Jan 08
" URL:		http://pari.math.u-bordeaux.fr
pi@pi400-64:~ $

I searched for Pari/GP samples on the web, and found some from Max Alekseyev:
https://home.gwu.edu/~maxal/gpscripts/

His invph.gp script shows that with "vi invphi.gp" syntax highlighting works:

pari_gp.syntax.png

pari_gp.syntax.png (27.18 KiB) Viewed 2623 times

I started with the first two simple RSA_numbers_factored.py helper functions:

Code: Select all

def bits(n: int) -> int:
...
def digits(n: int) -> int:

"digits()" seems to be a Pari/GP defined function:

Code: Select all

? ?digits
digits(x,{b=10}): gives the vector formed by the digits of x in base b (x and b 
integers).

? 

So I renamed "my" function to "digits_()".
This is initial script "r.gp":

Code: Select all

{ bits(N) = #binary(N); }

\\ { digits_(N) = ceil(log(N+1)/log(10)); }
{ digits_(N) = length(digits(N)); }

Using it:

Code: Select all

pi@pi400-64:~/RSA_numbers_factored/pari $ gp -q r.gp 
? bits(127)
7
? bits(128)
8
? digits_(99999)
5
? digits_(100000)
6
? 

Can be used non-interactive as well ("quit" exits "gp"), here "d.gp" demo:

Code: Select all

print(bits(65535), " ", bits(65536))
print(digits_(99), " ", digits_(100))
quit

Sample execution of both (I found no Pari/GP include/import concept yet):

Code: Select all

pi@pi400-64:~/RSA_numbers_factored/pari $ gp -q r.gp d.gp 
16 17
2 3
pi@pi400-64:~/RSA_numbers_factored/pari $ 

P.S:
I tested factorization of RSA-59 with Mathematica:
viewtopic.php?t=352146

It does not take hours for factorization like Python with sympy "factorint()", but 13.5× slower than 12 seconds with Pari/GP:

Code: Select all

pi@pi400-64:~ $ time wolframscript -code "FactorInteger[71641520761751435455133616475667090434063332228247871795429]"
{{200429218120815554269743635437, 1}, {357440504101388365610785389017, 1}}

real	2m42.880s
user	0m1.478s
sys	0m1.959s
pi@pi400-64:~ $ 

P.P.S:
"##" is a nice command, found in left column of 1st page of 4page Pari/GP reference card:
https://math.mit.edu/~brubaker/PARI/PARIrefcard.pdf

You just execute a command (like factorization of RSA-59 below).
If you then, after command completed, want to know the command's execution time, "##" does that for you!

Code: Select all

? factorint(71641520761751435455133616475667090434063332228247871795429)
%7 = 
[200429218120815554269743635437 1]

[357440504101388365610785389017 1]

? ##
  ***   last result: cpu time 12,052 ms, real time 12,080 ms.
? 

[HermannSW]

I learned how to program Pari/GP from:

• https://skalatan.de/en/archive/pariguid ... er_GP.html
• https://pari.math.u-bordeaux.fr/pub/par ... torial.pdf
• https://home.gwu.edu/~maxal/gpscripts/binomod.gp

No idea whether code style is good, tried to mimic what I saw.
76 line gist with some Robert Chapman Python functions transpiled to Pari/GP:
https://gist.github.com/Hermann-SW/96f6 ... 0ffc736027
Code demonstrates multiple assignments, indices in vector are 1-based (not used here):

pari_gp,code.png

pari_gp,code.png (30.78 KiB) Viewed 2584 times

I found no "assert()" in Pari/GP, implemented it:
https://gist.github.com/Hermann-SW/96f6 ... -sq2-gp-L8

Code: Select all

{ assert(b, v, s) = if(!(b), error(Str(v) Str(s))); }

"sq2(p)" determines sum of squares for prime p = 1 (mod 4).
Calling with 97 works (97 == 9^2 + 4^2), with non-prime 85 asserts (with assert() from line 8):

pari_gp.demo.png

pari_gp.demo.png (25.23 KiB) Viewed 2584 times

P.S:
Pari/GP "sq2()" quite fast, here for both prime factors of RSA-768:

pari_gp.demo2.png

pari_gp.demo2.png (26.14 KiB) Viewed 2573 times

[HermannSW]

I did function definitions as I had seen elsewhere.
Today I read two short tutorials from Pari/GP website and corrected my gist:
https://gist.github.com/Hermann-SW/96f6 ... 0ffc736027

• https://pari.math.u-bordeaux.fr/Events/ ... s/init.pdf
• https://pari.math.u-bordeaux.fr/Events/ ... s/prog.pdf

"my(...)" is used to declare all local variables, but is not needed if there are none.
vi syntax highlighting looks much nicer now:

pari_gp,code2.png

pari_gp,code2.png (23.65 KiB) Viewed 2534 times

I started to transpile RSA:numbers_factored.py to RSA_numbers_factored.gp.

Python "factorint()" returns a factorization dictionary, so I used that in the "rsa" array entries last two components:

Code: Select all

rsa = [
    [
        59,
        71641520761751435455133616475667090434063332228247871795429,
        200429218120815554269743635437,
        357440504101388365610785389017,
        {2: 2, 3: 2, 946790500267: 1, 5880369817360553: 1},
        {2: 3, 41: 1, 149: 1, 1356913: 1, 2739881: 1, 1967251783951: 1},
    ],
...

In Pari/GP "factorint()" returns a matrix instead, and "#M" is the number of columns:

Code: Select all

? M = factorint(735)
%1 = 
[3 1]

[5 1]

[7 2]

? #M
%2 = 2
?

Leftmost column can be used with "#" to determine number of rows (Pari/GP really starts index with 1 ...).
The "~" after the vector says that it is a column vector (transposed), and not a row vector:

Code: Select all

? M[,1]
%4 = [3, 5, 7]~
? #M[,1]
%5 = 3
? 

So transpiled to Pari/GP RSA-59 entry now looks like below (Pari/GP needs "\" at end if something is not complete).
Matrix rows are separated by ";", column entries by ",":

Code: Select all

rsa = [\
    [\
        59,\
        71641520761751435455133616475667090434063332228247871795429,\
        200429218120815554269743635437,\
        357440504101388365610785389017,\
        [2, 2; 3, 2; 946790500267, 1; 5880369817360553, 1],\
        [2, 3; 41, 1; 149, 1; 1356913, 1; 2739881, 1; 1967251783951, 1]\
    ],\
...

Quite some work to be done.

[HermannSW]

I did more investigation on Pari/GP factorint() capabilities.
After having factored 59 decimal digit number RSA-59 in 12s on Pi400, I tried to factor RSA-79.
That did run out of stack, because the default stack size is too low:

Code: Select all

pi@pi400-64:~ $ gp
Reading GPRC: /etc/gprc
GPRC Done.
...
parisize = 8000000, primelimit = 500000, nbthreads = 4
? 

So I did set as instructed in /etc/gprc:

Code: Select all

pi@pi400-64:~ $ head -17 /etc/gprc | tail -3
\\ Set maximal stack size. A safe value is half of the system memory.
\\ parisizemax = 1G
parisizemax = 2G
pi@pi400-64:~ $ 

I was able to determine the runtimes on i7-11850H as well, after having downloaded stable version
https://pari.math.u-bordeaux.fr/download.html
and following instruction in its INSTALL file.

In addition I did factorize RSA-59 and RSA-79 with msieve which I wrote about before:
viewtopic.php?t=342544

I updated the table of CPUs and test results for RSA-59/RSA-79 and RSA-100 for i7-11850H:
https://github.com/Hermann-SW/msieve#readme


Summary for factoring RSA-59 and RSA-79:

• unsurprisingly i7-11850H runtimes are much better than Pi400 runtimes
• msieve outperforms Pari/GP factorint() [completely on Pi400 for RSA-79]

Nevertheless it is good to have Pari/GP factorint() compared to hours of Python sympy factorint():

factorint_msieve.png

factorint_msieve.png (8.3 KiB) Viewed 2429 times

[HermannSW]

Pari/GP is such a performant number theory language (mostly because of the many cool algorithms making the difference).

So much happened in the last days, it is only two days ago that I registered for Pari/GP pari-users mailing list, and posted my first question regarding a memory and definitely performance problem of "gcd()":
https://pari.math.u-bordeaux.fr/archive ... 00003.html

Bill Allombert (Pari/GP package maintainer) answered, accepted the bug and fixed it on master branch.
I asked how I can use master branch and he pointed to instructions:
https://pari.math.u-bordeaux.fr/archive ... 00013.html


First I did build on i7-11850H laptop Linux, later on Pi400 with 64bit PiOS as well.

If you are interested in Computer Algebra System that is highly performant in number theory and open source, start here:
https://pari.math.u-bordeaux.fr/

There is a big amount of tutorials:
https://pari.math.u-bordeaux.fr/tutorials.html

Huge documentation:
https://pari.math.u-bordeaux.fr/doc.html

If you want to play with the calculator, you can without installing in your browser!
https://pari.math.u-bordeaux.fr/gp.html


There is so much I learned from Bill the last days.
Especially, that when knowing sum of squares for prime p=x^2+y^2, then you can compute sqrt(-1) (mod p) by "Mod(x/y, p)" (no division, but multiplication with reciprocal that is guaranteed to exist for prime p).
Also that you can use blazing fast "halfgcd(sqrtm1, p)" to determine sum of squares.
Also how to determine sqrtm1 and sum of squares directly only having p.
I determined runtimes, here for 100355-digit prime (333,333 bits):
https://pari.math.u-bordeaux.fr/archive ... 00019.html
GraphvizFiddle share link:
https://stamm-wilbrandt.de/100355-digit_prime.html

100355-digits_prime.png

100355-digits_prime.png (57.76 KiB) Viewed 2348 times

Last I went to biggest known twin prime pair and used bigger of them (388342-digits or 1,290,000 bits).
I used precomputed sqrtm1 value to determine the very efficient runtimes on the right (on i7-11850H):
https://pari.math.u-bordeaux.fr/archive ... 00020.html

388342-digits_prime.png

388342-digits_prime.png (49 KiB) Viewed 2348 times

Only 123ms/61ms for such big prime number computations.
Of course I did run that on Pi400 as well.
And yes, the i7-11850H is 5× faster, but the Pi400 computations are still sub second !

sqrtm1_qfbcornacchia.biggest_twin_prime.gp.png

sqrtm1_qfbcornacchia.biggest_twin_prime.gp.png (78.41 KiB) Viewed 2348 times

The Pari/GP script is in my repo (it asserts for all required identities), other scripts are in same "pari" directory:
https://github.com/Hermann-SW/RSA_numbe ... n_prime.gp
I also commited first version of RSA_numbers_factored.gp transpiled from RSA_numbers_factored.py ...

[HermannSW]

I learned from Bill how to get rid of the 2 phases on average for determining sqrtm1:
https://pari.math.u-bordeaux.fr/archive ... 00046.html

sqrtm1.smallest_known_1million_digit_prime.gp.png

sqrtm1.smallest_known_1million_digit_prime.gp.png (47.23 KiB) Viewed 2110 times

(I really like the infinite loop syntax "forprime(t=2,oo, ...)" of PARI/GP, especially "oo" for infinity)

Later I learned how to speedup "lift(Mod(x/y, p))", and do the efficient transformations on the right of above diagram in C++ / Python / PARI/GP:
https://pari.math.u-bordeaux.fr/archive ... 00074.html

xy_sqrtm1__xy.transformation_times.png

xy_sqrtm1__xy.transformation_times.png (22.42 KiB) Viewed 2110 times

This time on "smallest known prime with 1,000,000 decimal digits".


P.S:
The list of 5000 largest known primes is here:
https://t5k.org/primes/lists/all.txt

There is quite some updating activity on that list (new primes get added, pseudo primes that were proven not to be prime removed).
To stay informed on changes I created "latest_new_primes" gist that does some "diff" woodoo:
https://gist.github.com/Hermann-SW/eb39 ... e64b0ff0c5

latest_new_primes.2.png

latest_new_primes.2.png (27.45 KiB) Viewed 2110 times

[HermannSW]

Enough playing with "small" numbers of having 1million digits only

I targeted biggest known prime p that is =1 (mod 4), the 9,383,761-digit number being rank 10 of largest known primes list:
https://t5k.org/primes/lists/all.txt

Code: Select all

...
    9  2^32582657-1                     9808358 G9    2006 Mersenne 44
   10  10223*2^31172165+1               9383761 SB12  2016 
   11  2^30402457-1                     9152052 G9    2005 Mersenne 43
   ...

Last Sunday I stopped computation of "sqrt(-1) (mod p)" for that huge prime after 9d 20.9h, because I found llr tool (based on gwnum library) that was able to compute that value in 13.2h instead of expected 75.4 days!
https://github.com/Hermann-SW/9383761-d ... ime#readme

The determined "sqrtm1" constant is in this new PARI/GP script:
https://github.com/Hermann-SW/RSA_numbe ... 4_prime.gp

Determining sum of squares from sqrtm1 takes 20.8s, determining sqrtm1 from x and y takes 31.3s:

9383761.gp.pi400.png

9383761.gp.pi400.png (54.41 KiB) Viewed 1459 times

(7600X CPU PC needs 2.9s/4.2s, but that is OK)

[HermannSW]

On 7/25/2023 new largest prime p =1 (mod 4) was proven and published.
I did compute sqrt(-1) (mod p) for that 11,887,192-digit prime with patched LLR version again.
And this time computation did take 6.7 days to complete!
Here are the details:
https://github.com/Hermann-SW/11887192- ... motivation

While the PARI/GP script converting sqrt(-1) (mod p) to unique sum of squares and back does take few 100 milliseconds on AMD 7600X CPU, Pi400 runtimes are not that bad either for computing on numbers with more than 11million decimal digits or 39million bits !

Code: Select all

pi@pi400-64:~/RSA_numbers_factored/pari $ gp -q < sqrtm1.11887192_digit.largest_known_1mod4_prime.gp 
11887192-digit prime p (39488395 bits)
[M,V] = halfgcd(sqrtm1, p)
  ***   last result: cpu time 1,524 ms, real time 1,564 ms.
[x,y] = [V[2], M[2,1]]
  ***   last result: cpu time 8 ms, real time 8 ms.
sqrtm1 = lift(Mod(x, p)/y)
  ***   last result: cpu time 3,190 ms, real time 3,268 ms.
done, all asserts OK
pi@pi400-64:~/RSA_numbers_factored/pari $