## Pari/GP CAS (Computer Algebra System)

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Pari/GP CAS (Computer Algebra System)

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 9577 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!
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

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 9501 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.
? ``````
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

I learned how to program Pari/GP from: 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 9462 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 9462 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 9451 times
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

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 "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 9412 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.
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

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
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
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:

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 9307 times
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

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
https://stamm-wilbrandt.de/100355-digit_prime.html
100355-digits_prime.png
100355-digits_prime.png (57.76 KiB) Viewed 9226 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 9226 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 9226 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 ...
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

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 8988 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 8988 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 8988 times
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

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!

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 8337 times
(7600X CPU PC needs 2.9s/4.2s, but that is OK)
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

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 \$ ``````
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

PARI/GP has factorial operator:

Code: Select all

``````? 7!
%5 = 5040
? 7*6*5*4*3*2
%6 = 5040
? ``````

I requested addition of primorial operator "#".
Initial response was not agreeing.
But finally the operator got added to development branch:
https://pari.math.u-bordeaux.fr/archive ... 00010.html

Works as it should:

Code: Select all

``````? (3*4+1)#
30030
? 13#
30030
? 2*3*5*7*11*13
30030
?``````

And allows to play with complex primorial expressions of 3753-digit prime quadruplet:
(from list of 5000 largest known primes: https://t5k.org/primes/lists/all.txt)

Code: Select all

``````\$ ./gp -q
? \v
GP/PARI CALCULATOR Version 2.16.1 (development 28704-e9753897f)
amd64 running linux (x86-64/GMP-6.1.2 kernel) 64-bit version
...
? p=(1049713153083*2917#*(567*2917#+1)+2310)*(567*2917#-1)/210+9;
? #digits(p)
3753
? [ispseudoprime(p-d)|d<-[0,2,6,8]]
[1, 1, 1, 1]
?``````

I had built PARI/GP dev from git according tutorial, with:

Code: Select all

``./Configure --mt=pthread``
https://pari.math.u-bordeaux.fr/Events/ ... ources.pdf

Code: Select all

``````? \v
GP/PARI CALCULATOR Version 2.16.1 (development 28707-376901920)
arm64 running linux (aarch64/GMP-6.2.1 kernel) 64-bit version
compiled: Oct  2 2023, gcc version 10.2.1 20210110 (Debian 10.2.1-6)
(readline v8.1 enabled, extended help enabled)
?
``````

Today I first time played with "parallel for" — parallel programming can be so easy!

Code: Select all

``````? ?parfor
parfor(i=a,{b},expr1,{r},{expr2}): evaluates the expression expr1 in parallel
for all i between a and b (if b is set to +oo, the loop will not stop),
resulting in as many values; if the formal variables r and expr2 are present,
evaluate sequentially expr2, in which r has been replaced by the different
results obtained for expr1 and i with the corresponding arguments.

?``````

As baseline I printed the 25 primes below 100 sequentially first:

Code: Select all

``````? for(i=2,100,if(isprime(i),print1(i,",")))
2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
? ``````

That was easy, now same with parfor, and each execution shows a different output order:

Code: Select all

``````pi@raspberrypi400:~ \$ gp -q
? parfor(i=2,100,isprime(i)*i,b,if(b,print1(b,",")))
2,3,7,11,13,19,23,17,29,31,37,43,41,47,59,61,53,67,71,73,79,5,89,97,83,
? parfor(i=2,100,isprime(i)*i,b,if(b,print1(b,",")))
2,3,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,5,83,89,97,
? parfor(i=2,100,isprime(i)*i,b,if(b,print1(b,",")))
2,3,7,11,13,17,19,23,29,31,5,43,37,47,41,53,59,61,67,71,73,79,83,89,97,
? parfor(i=2,100,isprime(i)*i,b,if(b,print1(b,",")))
2,3,7,5,11,13,17,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,19,
? ``````
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

I asked on pari-users email list wrt drawing in PARI/GP and got help.

This is plt.3.gp.
The foreach assert in line 5 verifies all hyperellratpoints property.

Code: Select all

``````assert(b,s)=if(!(b), error(Str(s)));
h=eval(getenv("h"));
e=eval(getenv("e"));
H=hyperellratpoints(e,h);
foreach(H,R,assert(subst(e,z,R[1])==R[2]^2));

print([R|R<-H,type(R[2])=="t_INT"]);
X=[R[1]|R<-H];
Y=[R[2]|R<-H];
Xi=[R[1]|R<-H,type(R[2])=="t_INT"];
Yi=[R[2]|R<-H,type(R[2])=="t_INT"];``````

Code: Select all

``````print(vecmin(X)," ",vecmax(X)," ",vecmin(Y)," ",vecmax(Y));
dx=(vecmax(X)-vecmin(X))/400;
dy=(vecmax(Y)-vecmin(Y))/240;

plotinit(1);
plotscale(1,vecmin(X)-2*dx,vecmax(X)+2*dx,vecmin(Y)-2*dy,vecmax(Y)+2*dy);
for(i=1,#X,plotmove(1,X[i]-dx/2,Y[i]-dy/2);plotrbox(1,dx,dy,1));
for(i=1,#Xi,plotmove(1,Xi[i]-3*dx/2,Yi[i]-3*dy/2);plotrbox(1,3*dx,3*dy,1));
plotdraw(1);``````

You execute it like this:

Code: Select all

``````pi@raspberrypi400:~ \$ e="353*z^2 + 188*z + 36" h="40" gp -q < plt.3.gp
[[-5, 89], [-5, -89], [0, 6], [0, -6], [5, 99], [5, -99]]
-5 5 -99 99
pi@raspberrypi400:~ \$``````

This is screenshot of the drawing window.
The plotrbox() does draw points of different size.
I am interested in the big points, [z,w] values in hyperellratpoints() result with integer w:
hyperellratpoints.jpg (28.3 KiB) Viewed 3166 times

There exists definitely better software for drawing.
But the commands are part of PARI/GP and easily allow to display stuff computed with PARI/GP.

(you can get short help with eg. ?plotmove, longer help with ??plotmove, and ???plot lists all commands related to plot)

P.S:
This thread tells you which font package to install so that needed 7x13 font is available:
viewtopic.php?t=360459
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

For my planar graph playground work I did output OpenSCAD files, as well as later JSCAD files:
viewtopic.php?t=333342&start=75#p2027793

Recently I started new repo for providing 3D output for PARI/GP:
https://github.com/Hermann-SW/GP3D

JSCAD is cool by its own, just open JSCAD menu top left and click on the different examples:

The website allowed to pass a model in the URL (https://jscad.app/#http...).
And it did allow to pass base64 data url as well.
I opened a feature request issue on allowing for "...#data,application/gzip;base64,..." as well.
Today the pull request was pushed to server:

I updated my tool to make use of new mechanism:
https://github.com/Hermann-SW/GP3D/blob ... view_gzb64

Code: Select all

``````#!/bin/bash
browser=firefox
browser=chromium-browser

\$browser "https://jscad.app/#data:application/gzip;base64,`gzip -c \$1 | base64 -w 0`"``````

I did run with 27303-bytes model I created with GP3D:

Code: Select all

``````pi@raspberrypi5:~/GP3D \$ tools/view_gzb64 gp.jscad
Opening in existing browser session.
pi@raspberrypi5:~/GP3D \$``````
The resulting base64(gzip(gp.jscad)) URL is a 6552 bytes data URL, just click and you can see my model in your browser.

If you select the matching options in dialog left bottom, it should look similar to this:
JSCAD.half_vertex.png (119.22 KiB) Viewed 1267 times

P.S:
I did performance measurements for that model today, on jscad.app and on openjscad.xyz.
It turned out that my Intel Celeron j4125 miniPC running Ubunutu is 50% slower than my Pi5 (default setup, no overclocking).
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

I worked a lot on GP3D, and learned in the meantime how to output "same" vertex fast in 3D by translate and using copy that was rendered once. Since recently jscad.app has "smooth render" option that makes view so much nicer.
Adding parameters is so easy, and allows to see what I really wanted.
With disabling "mod" for colors from palette I was able to see the sloped same level (same color) planes I searched for.
This animated .gif was recorded with peek screenrecorder (under X11 and not Wayfire) of chromium browser on Raspberry Pi5, generated with "n=61 make tqf_3D.2" of current GP3D repo:

P.S:
Just checked, and I am running Pi5 overclocked currently:

Code: Select all

``````pi@raspberrypi5:~/GP3D \$ freq
min=cur=3000000=max
pi@raspberrypi5:~/GP3D \$
``````
There are currently quite some demos in the repo, but README.md does not reflect that at present.
Makefile has these demo targets:

Code: Select all

``````pi@raspberrypi5:~/GP3D \$ grep : Makefile
demo:
demo2:
demo3:
sphere_embedding:
tqf_3D.1:
tqf_3D.2:
tqf_3D.3:
clean:
pi@raspberrypi5:~/GP3D \$
``````
"tqf*" demos require "n" environment variable to be set, and error out if not.

P.P.S:
"make demo2" demonstrates the spherical embedding (slower) vertices, great circle and straight edges, spherical triangle fill and vertex label features of GP3D:
demo2.gp.jpg
demo2.gp.jpg (31.27 KiB) Viewed 625 times

Only 18 lines of PARI/GP are needed for this demo:
demo2.gp.png
demo2.gp.png (68.07 KiB) Viewed 628 times
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

pidd
Posts: 5577
Joined: Fri May 29, 2020 8:29 pm
Location: Wirral, UK

### Re: Pari/GP CAS (Computer Algebra System)

Have you compared the speed of pari against gmp/BigNum for factorisation or prime test? I haven't tried gmp with the Pi5 yet, it wasn't very fast on the Pi4, it only uses a single thread.

I achieved a T5K a couple of weeks ago through SRBase, unfortunately I didn't notice until after the five days registration window so it was then registered as anonymous.

210060*91^331939-1

210060*91^374955-1 was also found by someone else on the same day, I've not seen similar primes occur on the same day before.

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

pidd wrote:
Tue Feb 06, 2024 3:25 am
Have you compared the speed of pari against gmp/BigNum for factorisation or prime test? I haven't tried gmp with the Pi5 yet, it wasn't very fast on the Pi4, it only uses a single thread.
PARI/GP is based on GMP by default (aarch64/GMP-6.2.1 kernel):

Code: Select all

``````? \v
GP/PARI CALCULATOR Version 2.15.4 (released)
arm64 running linux (aarch64/GMP-6.2.1 kernel) 64-bit version
compiled: Dec  9 2023, gcc version 12.2.0 (Debian 12.2.0-14)
(readline v8.2 enabled, extended help enabled)
?
``````
But for factorization GMP is not good, I once determined "sqrt(-1) (mod p)" 9.1million decimal digit prime p (largest prime =1 (mod 4) at that time) on 6C/12T AMD 7600X CPU, but as you said GMP is single threaded. I stopped after 10 days (knowing expected runtime would be 75.4 days). I did computation with patched version of Jean Penne's LLR software instead, in only 10.8h(!):
https://github.com/Hermann-SW/9383761-d ... -p-1-mod-4

For factoring primes I use multithreaded cado-nfs, that was able to factor RSA-100 on Pi400 in 12.1h first, then overclocked and with a ARM bug I reported fixed, took only 1.8h(!) on Pi400 :
Later I did measure factoring times on Pi400 and only 13:30min on overclocked Pi5(!):

Code: Select all

``````|     |  Pi400  |   Pi5   |
| RSA |  2.2GHz |   3GHz  |
|-----+---------+---------+
| 100 | 1:49:54 | 0:13:30 |
| 110 | 4:41:35 | 1:27:27 |
| 120 |    —    | 5:05:17 |``````
Of course factoring is much much faster on x86 CPUs:
https://github.com/Hermann-SW/RSA_numbe ... ring-rsa-x

I achieved a T5K a couple of weeks ago through SRBase, unfortunately I didn't notice until after the five days registration window so it was then registered as anonymous.

210060*91^331939-1

210060*91^374955-1 was also found by someone else on the same day, I've not seen similar primes occur on the same day before.
Cool -- how long did the computation take for you on what CPU?

Code: Select all

``````? #digits(210060*91^331939-1)
%1 = 650288
?
``````
I just used Jean's LLR 4.0.5:
http://jpenne.free.fr/
Runs on x86 architecture only because of gwnum lib, which is able to parallelize computation of single multiplication(!) of big numbers.
I learned that I have to force LLR to run on chiplet0 of my 16C/32T 7950X CPU, just did that and proved your number as prime in only 21min ...

Code: Select all

``````hermann@7950x:~\$ taskset -c 0-7,16-23 ./sllr64 -t16 -d -q"210060*91^331939-1"
Base factorized as : 7*13
Base prime factor(s) taken : 13
Resuming N+1 prime test of 210060*91^331939-1 at bit 1005663 [46.55%]
Using zero-padded AVX-512 FFT length 288K, Pass1=768, Pass2=384, clm=2, 16 threads, a = 3
210060*91^331939-1 may be prime. Starting Lucas sequence...
Using zero-padded AVX-512 FFT length 288K, Pass1=768, Pass2=384, clm=2, 16 threads, P = 6
210060*91^331939-1 may be prime, trying to compute gcd's
U((N+1)/13) is coprime to N!
210060*91^331939-1 is prime! (650288 decimal digits, P = 6)  Time : 1260.674 sec.
hermann@7950x:~\$ s
``````
This is snapshot in 2nd Lucas phase:

Code: Select all

``````hermann@7950x:~\$ ps -o etime= -p 2591
03:20
hermann@7950x:~\$
top - 10:17:39 up 13 min,  2 users,  load average: 10.34, 5.59, 2.39
Tasks: 464 total,   2 running, 462 sleeping,   0 stopped,   0 zombie
%Cpu(s):  0.0 us,  3.7 sy, 25.2 ni, 71.0 id,  0.0 wa,  0.0 hi,  0.0 si,  0.0 st
MiB Mem :  31190.9 total,  29557.6 free,    798.1 used,    835.2 buff/cache
MiB Swap:   2048.0 total,   2048.0 free,      0.0 used.  30006.2 avail Mem

PID USER      PR  NI    VIRT    RES    SHR S  %CPU  %MEM     TIME+ COMMAND
2591 hermann   39  19 1995820  84544   2816 R  1068   0.3  36:57.92 sllr64  ``````
"only" 1068% CPU, not enough work for 16 threads.
I recently figured out that >950,000 decimal digits are needed to see 16 threads being used:

So what do I use PARI/GP for?
Not factoring, and not prime proving.
Although I learned that its precomputed prime tables allow to really fast use of "prime(n)" to compute the n-th prime :
https://pari.math.u-bordeaux.fr/archive ... 00030.html

Code: Select all

``````pi@raspberrypi5:~ \$ freq
min=cur=3000000=max
pi@raspberrypi5:~ \$ gp -q -p4M
? for(n=1,2^18, prime(n));
? ##
***   last result: cpu time 1,023 ms, real time 1,028 ms.
? prime(2^18)
3681131
?
``````
I use it for complex mathematical number theoretic functions.
This small code does compute all xyz coordinates for n=61:
https://github.com/Hermann-SW/GP3D/blob ... gp#L32-L41

Code: Select all

``````    assert(getenv("n")!=0);
n=eval(getenv("n"));
print("n=",n);
assert(n%4!=0&&n%8!=7);

Q=qflllgram(get_tqf(n))^-1;
M=Q~*Q;
S=[x|x<-Vec(qfminim(M,n)[3]),qfeval(M,x)<=n];

``````
(the get_tqf() function is 25 lines long computation based on work from Dirichlet from 1850
https://github.com/Hermann-SW/GP3D/blob ... .gp#L5-L29 )

If I would replace "qfeval(M,x)<=n" by "qfeval(M,x)==n" I would get only those 3-tuples [x,y,z] with x^2+y^2+z^2==n.
I can see those by enabling white sphere surface with radius sqrt(61) for above animation:
tqf_3D.2.gp.jpg
tqf_3D.2.gp.jpg (38.26 KiB) Viewed 567 times

P.S:
One last thing, why did I buy Pi5?
Not because of PARI/GP, hat runs on PiOS as well as faster on x86_64 CPUs.
But wolframscript is free to use on any Raspberry computer, and Pi5 overclocked with 3GHz (I had luck in chip lottery) is 2.8x faster than my Pi400 for long wolframscript computations like:
viewtopic.php?p=2166024#p2166024

Code: Select all

``Solve[14*a^2+(130*b+2*c)*a+(302*b^2+101*c^2)==101,{a,b,c},Integers]``
(the "Integers" restriction is the hard part)
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

pidd
Posts: 5577
Joined: Fri May 29, 2020 8:29 pm
Location: Wirral, UK

### Re: Pari/GP CAS (Computer Algebra System)

Wolfram is amazing, my main use has been for its symbolic computation which avoids my mistakes when simplifying or re-arranging formula.

I love the stuff you do but I'm no longer capable of understanding all of it. These days I mostly donate computer time but also play with lossless binary compression which has been a lifetime project/interest. Programming is getting harder as my short-term memory is not what it was.

The T5K was either on a Ryzen-7 5800H or I7-4790S, I didn't catch it in time to see the online report, I will try and look through my local logs for it, I had five computers running 52 tasks at a time on that project. I have Pi4's and Pi5's running as well which do some tasks very well because of their incredible integer speed, unfortunately few projects are optimised for them and memory access tends to be a bottleneck as well.

Whatever time it was, it was a lot longer than 14 minutes, that is absolutely incredible for a prime of that size, it makes you wonder what Governments can actually do these days.

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

I did determine maximal number of vertices in a single plane (same level, same color), then determined normal vector for that plane (I learned on PARI/GP email list how to do that with this nice short command) ...

Code: Select all

``v=matker(Mat(apply(x->concat(x,1),P~)))[,1];``
...and the point in plane with minimal distance to origin (the small black sphere). Finally I did draw the plane itself transparent. There was no primitive for that in JSCAD, but using a very thin cuboid did the job! I really like the alpha slider for immediately changing transparency. This model was generated with:
https://github.com/Hermann-SW/GP3D/blob ... qf_3D.2.gp

Code: Select all

``````pi@raspberrypi5:~/GP3D \$ n=37 make tqf_3D.2
gp -q < tqf_3D.2.gp
n=37
pi@raspberrypi5:~/GP3D \$
``````
tqf_3D.2.gp.n_37.jpg
tqf_3D.2.gp.n_37.jpg (65.29 KiB) Viewed 293 times
The idea of generating JSCAD model from PARI/GP was that you can view, translate, rotate and zoom in your browser. Works nicely. I used my personal website as URL shortener for the 7543 bytes long application/gzip jscad.app data URL. Click, and you can play with the model in your browser:
https://stamm-wilbrandt.de/tqf_3D.2.gp.n_37.html

HermannSW wrote:
Tue Feb 06, 2024 12:46 am
This animated .gif was recorded with peek screenrecorder (under X11 and not Wayfire) of chromium browser on Raspberry Pi5, generated with "n=61 make tqf_3D.2" of current GP3D repo:
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/

HermannSW
Posts: 6247
Joined: Fri Jul 22, 2016 9:09 pm
Location: Eberbach, Germany

### Re: Pari/GP CAS (Computer Algebra System)

I like prototypical jscadui parameter development (that small model was generated with "n=3 gp -q < tqf_3D.2.gp").
First I had "plane" choice on whether to show plane (through max #vertices) or not.
Now I wanted to show plane through origin with mapped point from single red vertex (as very dark gray) as well.
I changed "plane" to choice and adapted logic, was only few lines (PARI/GP) code change that generates the JSCAD model ...
https://github.com/Hermann-SW/GP3D/comm ... 55cedc32a0
... I like it

And the generated jscad model in browser ui editor shows me the computed values (small black/very dark gray vertices are PO/DO):

Code: Select all

``````...
function main(params) {
for(i=0;i< 7 ;++i)
palette[4+i]=hexToRgb(params['p_'+i])
PO = [-1/9, 2/9, 2/9]
DO = [8/9, 11/9, -7/9]
N =  [116.56505117707798935157219372045329467, 48.189685104221401934142083269421740469]
...``````
https://github.com/Hermann-SW/RSA_numbers_factored
https://stamm-wilbrandt.de/GS_cam_1152x192@304fps
https://hermann-sw.github.io/planar_graph_playground
https://github.com/Hermann-SW/Raspberry_v1_camera_global_external_shutter
https://stamm-wilbrandt.de/