20100815, 23:49  #1 
Aug 2006
1761_{16} Posts 
Pari bootcamp
This is a thread for discussion of Pari/GP, the Pari library, gp scripts, gp2c, and related topics.
For anyone new to Pari, here's my ultrashort "getting started" guide. Look up your operating system and follow the instructions. 1. Windows: Go to the download page and click on the selfinstalling binary distributions for Windows, in particular, the development version. If you download the stable distribution, just delete it and download the development version instead. Then doubleclick the installer, and click "next" or the like until it's done. You may also wish to rightclick your shortcut and check "Quick Paste mode" or the like, which will allow you to copy by selecting with the mouse and pressing Enter and paste by rightclicking. 2. Mac OS X: Go to the Darwin ports page and follow the instructions. 3. Ubuntu, Debian, Mint, etc.: Go to Add/Remove Applications, the Ubuntu Software Center, the Synaptic Package Manager, the Advanced Packaging Tool, etc. as appropriate, search for PARI/GP, and click Install. 4. Linux (other than the above, or for advanced users): Compile it yourself! I have instructions that may or may not be useful to you. Once you have Pari/GP installed, you can use it as a calculator (but with advanced number theory features). Even better, you can extend the program with your own scripts. For example, here's my program for calculating the radical of an integer: Code:
rad(n:int)={ if(n < 0, n = n); if (issquarefree(n), return(n)); vecprod(factor(n)[,1]) }; addhelp(rad, "rad(n): Radical of n, the largest squarefree number dividing n. Sloane's A007947."); Code:
read "foobar.gp" To make editing your scripts easier I have a small collection of syntax highlighting files for various text editors. Disclaimer: Pari isn't my program; I have no connection to the Pari group that maintains it. In particular, I'm neither Henri Cohen nor Karim Belabas. 
20100816, 00:57  #2 
Aug 2006
5985_{10} Posts 
I'll post some examples of my scripts. I encourage anyone following the thread to do the same! Simple scripts can help beginners, complicated scripts can help experts, and even hairy awful hardtoread scripts can be useful. (I have a section in my script file called "Verbose monstrosities" for just that sort of code...)
I strongly recommend getting in the habit of including addhelp() commands by every function, so you can be reminded of the purpose, the limitations of the algorithm(s) used (if any), and the input format. Code:
Faulhaber(e:small,a='x)={ substpol(sum(i=0,e,binomial(e+1,i)*bernfrac(i)*'x^(e+1i))/(e+1) + 'x^e, 'x, a) }; addhelp(Faulhaber, "Faulhaber(e,{a='x}): Returns the polynomial for the sum 1^e + 2^e + ... + x^e, evaluated at a."); countsemi(n:int)={ my(s=0,B=sqrtint(n)); forprime(p=2,B, s+=primepi(n\p) ); B=primepi(B); sB*(B1)/2 }; addhelp(countsemi, "countsemi(n): Number of semiprimes up to n. Sloane's A072000."); countPowerful(lim)={ sum(k=1,(lim+0.5)^(1/3), if(issquarefree(k), sqrtint(lim\k^3)) ) }; addhelp(countPowerful, "countPowerful(lim): Number of powerful numbers up to lim."); isTriangular(n:int)={ issquare(n<<3+1) }; addhelp(isTriangular, "isTriangular(n): Is n a triangular number? Sloane's A010054."); isFibonacci(n:int)={ my(k=n^2); k+=((k + 1) << 2); issquare(k)  (n > 0 && issquare(k8)) }; addhelp(isFibonacci, "isFibonacci(n): Is n a Fibonacci number? Sloane's A010056; characteristic function of Sloane's A000045."); rp(b:small)={ b=2^b; nextprime(b+random(b)) }; addhelp(rp, "rp(b): Returns a random bbit prime."); rsp(b:small,flag:small=0)={ if(b<26, my(n); return(if(b<8, if(b<3, error("No semiprimes with fewer than 3 bits")); if(b==3,n=[4,6]); if(b==4,n=[9,10,14,15]); if(b==5,n=[21,22,25,26]); if(b==6,n=[33,34,35,38,39,46,49,51,55,57,58,62]); if(b==7,n=[65,69,74,77,82,85,86,87,91,93,94,95,106,111,115,118,119,121,122,123]); n=n[random(#n)+1]; if(flag,factor(n),n) , b=2^(b1); while(bigomega(n=random(b)+b)!=2,); \\ Brute force n )) ); \\ Choose the small prime as 2 with probability 1/2, 3 with probability 1/3, etc. my(M=.261497212847642783755,weight=log(b*log(2)/2)+M,r=random(round(weight << 64))*1.>>64,p=precprime(exp(exp(rM))),q); q=2^(b1)\p; q=nextprime(q+random(q)); if(p > q, b = p; p = q; q = b; ); if(flag, if(p==q, matrix(1,2,i,j,if(j==1,p,2)) , [p,1;q,1] ) , p*q ) }; addhelp(rsp, "rsp(b, {flag=0}): Generates a random bbit semiprime. If flag is nonzero, returns the factorization of the number."); deBruijnXi(x)={ my(left, right, m); if (x < 1, error ("deBruijnXi: Can't find a xi given x < 1.")); if (x > 1, left = log(x), left = eps()); right = 1.35 * log(x) + 1; \\ Heuristic \\Bisection while (right  left > left * eps(), m = (left + right) / 2; if (exp(m)  1 > x * m, right = m, left = m) ); (left + right) / 2 }; addhelp(deBruijnXi, "deBruijnXi(x): Helper function for rhoest. Finds a xi such that e^xi  1 = x * xi."); rhoest(x)={ my(xi = deBruijnXi(x)); \\exp(Euler) / sqrt(2 * Pi * x) * exp(1  exp(xi) + intnum(s = eps(), xi, (exp(s)  1) / s)) exp(eint1(xi)  x * xi) / sqrt(2 * Pi * x) / xi }; addhelp(rhoest, "de Bruijn's asymptotic approximation for rho(x), rewritten as in van de Lune and Wattel 1969. Curiously, their paper shows values for this estimate that differ from those calculated by this function, often as soon as the second decimal place  but as the difference is in the direction of the true value, I have not looked further into this."); rhoTable = [1, 3.068528194e1, 4.860838829e2, 4.910925648e3, 3.547247005e4, 1.964969635e5, 8.745669953e7, 3.232069304e8, 1.016248283e9, 2.770171838e11, 6.644809070e13, 1.419713165e14, 2.729189030e16, 4.760639989e18, 7.589908004e20]; DickmanRho(x)={ local(left, right, scale); if (x <= 2, return (1  log(max(x, 1)))); if (x <= 3, return( 1  (1  log(x  1))*log(x) + real(dilog(1  x)) + Pi^2 / 12 )); \\ Asymptotic estimate (scaled for continuity) if (x > #rhoTable, scale = rhoTable[#rhoTable] / rhoest(#rhoTable); \\ Let the scale factor dwindle away, since the estimate is (presumably) \\ better in the long run than any scaled version of it. The exponent \\ of 0.25 has been chosen to give the best results for 10 < x < 100 \\ with a table size of 10. scale = (scale  1) * (#rhoTable / x)^.25 + 1; return (precision(rhoest(x) * scale, 9)) ); \\ Scaling factors: the factor by which the true value of rho differs from \\ the estimates at the endpoints. left = rhoTable[floor(x)] / rhoest(floor(x)); right = rhoTable[ceil(x)] / rhoest(ceil(x)); \\ Linear interpolation on the scale factors. scale = left + (right  left) * (x  floor(x)); \\ Return a result based on the scale factor and the asymptotic formula. precision(rhoest(x) * scale, 9) }; addhelp(DickmanRho, "Estimates the value of the Dickman rho function. For x <= 3 the exact values are used, up to rounding; up to "#rhoTable" the value is interpolated using known values and rhoest; after "#rhoTable" rhoest is used, along with a correction factor based on the last value in rhoTable."); For those curious: the type information is ignored by gp. I include it in my script files in case I want to compile it with gp2c later (or if I've already done so...). 
20100816, 20:56  #3 
Aug 2006
3^{2}·5·7·19 Posts 
In response to a PM:
You can check the current primelimit, the size of the largest precomputed prime, with Code:
default(primelimit) Code:
ff(x)={ if(x>default(primelimit), error("primelimit too small, need "x) ); my(s=0.); forprime(p=2,x, s += 1/p ); s }; 
20100816, 21:12  #4 
"Forget I exist"
Jul 2009
Dumbassville
20300_{8} Posts 
Code:
Primorial(n)=if(n>default(primelimit),a=1;forprime(x=1,default(primelimit),a=a*x);for(x=default(primelimit),n,if(isprime(x),a=a*x)),a=1;forprime(s=1,n,a=a*s));return(a) Note: CRG's post really helped me figure this out. Last fiddled with by science_man_88 on 20100816 at 21:12 
20100816, 21:26  #5 
"Robert Gerbicz"
Oct 2005
Hungary
2×733 Posts 

20100816, 21:37  #6  
Aug 2006
3^{2}·5·7·19 Posts 
Quote:
Code:
b; I haven't entirely decided how to handle the fencepost issues. Both rp() and rsp() can give outofrange values, but for large b this is increasingly rare. There are many other design decisions relating to quality vs. speed: for example, it would clearly be inappropriate to use rp() to determine the empirical density of the upper members of twin primes. Code:
>sum(i=1,1e5,isprime(rp(10)+2)) %1 = 16590 >sum(i=1,1e5,isprime(rp(10)2)) %2 = 4788 Last fiddled with by CRGreathouse on 20100816 at 21:51 

20100816, 21:44  #7  
Aug 2006
3^{2}·5·7·19 Posts 
Quote:
For what it's worth, here's my messybutfast code: Code:
primorial(n)={ my(t1=1,t2=1,t3=1,t4=1); forprime(p=2,n>>3,t1=t1*p); forprime(p=(n>>3)+1,n>>2,t2=t2*p); t1=t1*t2;t2=1; forprime(p=(n>>2)+1,(3*n)>>3,t2=t2*p); forprime(p=((3*n)>>3)+1,n>>1,t3=t3*p); t1=t1*(t2*t3);t2=1;t3=1; forprime(p=(n>>1)+1,(5*n)>>3,t2=t2*p); forprime(p=((5*n)>>3)+1,(3*n)>>2,t3=t3*p); t2=t2*t3;t3=1; forprime(p=((3*n)>>2)+1,(7*n)>>3,t3=t3*p); forprime(p=((7*n)>>3)+1,n,t4=t4*p); t1*(t2*(t3*t4)) }; addhelp(primorial, "Returns the product of each prime less than or equal to n. Sloane's A034386."); Code:
>primorial(10^6); time = 220 ms. >Primorial(10^6); time = 1,830 ms. Quote:
Code:
primorial2(n)=primorial(prime(n)) 

20100816, 21:52  #8 
"Forget I exist"
Jul 2009
Dumbassville
2^{6}·131 Posts 
I'll have to check the functions in your code just to refresh, the hard part using primorial2 Crazy Request God is if you want a prime(n) such that it is out of range of primelimit you still need to add code to check for such and or use a way to calculate prime(n) or use isprime to find it. but I think I can make this work lol.
Last fiddled with by science_man_88 on 20100816 at 21:56 
20100816, 22:05  #9  
"Robert Gerbicz"
Oct 2005
Hungary
2×733 Posts 
Quote:
Code:
prodtree(A)=local(i,len,mid);len=length(A);if(len<=8,return(prod(i=1,len,A[i])));\ mid=len\2;return(prodtree(vector(mid,i,A[i]))*prodtree(vector(lenmid,i,A[i+mid]))) Primorial2(n)=prodtree(primes(primepi(n))) Code:
? default(primelimit,10^7) %117 = 10000000 ? gettime();x=Primorial2(10^7);gettime() %118 = 8580 ? gettime();y=primorial(10^7);gettime() %119 = 68219 ? xy %120 = 0 

20100816, 22:23  #10  
Aug 2006
3^{2}·5·7·19 Posts 
Quote:
Code:
primorialNaive(x)=my(pr=1,n=2);while(n<=x, pr *= n; n = nextprime(n+1)); pr Time for default(primelimit, 10^7): 30 milliseconds Time for Primorial(1e7): 319 seconds OK, so in this calculation it doesn't much matter, since all the time is spent moving around huge chunks of memory with the partial product. But with a more efficient function it really starts to make a difference. And certainly you can see that the time needed for setting your primelimit to something reasonable isn't excessive, compared to crunching through directly. Code:
noop(x)=forprime(p=2,x,) noopNaive(x)=my(n=2);while(n<=x, n = nextprime(n+1)) 

20100816, 22:26  #11  
Aug 2006
3^{2}·5·7·19 Posts 
Quote:
But yes, I get ~7 times faster with your version than with mine when I increase the stack appropriately. Last fiddled with by CRGreathouse on 20100816 at 22:32 

Thread Tools  
Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Pari/GP to PFGW  carpetpool  Programming  6  20171221 06:04 
LLL in GP/Pari  paul0  Programming  2  20151117 13:04 
PARI vs GAP  skan  Miscellaneous Math  0  20121216 00:13 
pari  devarajkandadai  Programming  21  20120831 18:08 
64bit Pari?  CRGreathouse  Software  2  20090313 04:22 