Pari bootcamp

[CRGreathouse]

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 ultra-short "getting started" guide. Look up your operating system and follow the instructions.
1. Windows: Go to the download page and click on the self-installing 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 double-click the installer, and click "next" or the like until it's done. You may also wish to right-click 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 right-clicking.
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:

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.");

As you write your own programs, I strongly encourage you to save your work in a script file. You can load these functions from Pari with the read command (abbreviated "\r" if you really need to save those extra characters). In fact, you can even skip that step if you add the line

  read "foobar.gp"

(where "foobar.gp" is the name of your script file) to the end of your .gprc file, which is in your Pari directory (usually C:\Program Files\PARI\ in Windows).

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.

[CRGreathouse]

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 hard-to-read 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.

Faulhaber(e:small,a='x)={
	substpol(sum(i=0,e,binomial(e+1,i)*bernfrac(i)*'x^(e+1-i))/(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);
	s-B*(B-1)/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(k-8))
};
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 b-bit 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^(b-1);
			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(r-M))),q);
	q=2^(b-1)\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 b-bit 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.068528194e-1, 4.860838829e-2, 4.910925648e-3, 3.547247005e-4, 1.964969635e-5, 8.745669953e-7, 3.232069304e-8, 1.016248283e-9, 2.770171838e-11, 6.644809070e-13, 1.419713165e-14, 2.729189030e-16, 4.760639989e-18, 7.589908004e-20];
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.");

I have omitted most of the copious comments on the function rsp(), largely those which are just ideas for improving it. I think, for most purposes, this version is workable as written.

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...).

[CRGreathouse]

In response to a PM:

You can check the current primelimit, the size of the largest precomputed prime, with

default(primelimit)

So for example,

ff(x)={
  if(x>default(primelimit),
    error("primelimit too small, need "x)
  );
  my(s=0.);
  forprime(p=2,x,
    s += 1/p
  );
  s
};

[science_man_88]

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)

Though if you choose a n above default(primelimit) it slows down. One other thing I notice is since if you want the first n primes multiplied instead of all primes up to N it's not what you want. So I may try to add a variable in that if it is included it tells how many primes to multiply together rather than n but I bet this will make it tricky for me to code.

Note: CRG's post really helped me figure this out.

[R. Gerbicz]

Originally Posted by CRGreathouse

[code]

rp(b:small)={
b=2^b;
nextprime(b+random(b))
};
addhelp(rp, "rp(b): Returns a random b-bit prime.");

That's never true, depending on the random value it will be b+1 or b+2 bit prime number.

[CRGreathouse]

Originally Posted by R. Gerbicz

That's never true, depending on the random value it will be b+1 or b+2 bit prime number.

Good call. For some reason I left the line

  b--;

out; maybe I was copying from an older version. (My actual version, like for most of these functions, is written in C with the Pari library.)

I haven't entirely decided how to handle the fencepost issues. Both rp() and rsp() can give out-of-range 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.

>sum(i=1,1e5,isprime(rp(10)+2))
%1 = 16590
>sum(i=1,1e5,isprime(rp(10)-2))
%2 = 4788

The main purpose of the rp function was to provide random semiprimes for cranks with new factorizations methods, so precision wasn't important. I'd just do rp(100)*rp(100) and paste the result.

[CRGreathouse]

Originally Posted by science_man_88

One other thing I notice is since if you want the first n primes multiplied instead of all primes up to N it's not what you want.

I'm not going to say the other way is wrong, it's a A002110 vs. A034386 thing. But "product of the primes up to x" is much more common, and if you want to use the Dubner # you need to use that convention.

For what it's worth, here's my messy-but-fast 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.");

Here's a slightly-unfair comparison:

>primorial(10^6);
time = 220 ms.
>Primorial(10^6);
time = 1,830 ms.

(about 8x speedup)

Originally Posted by science_man_88

So I may try to add a variable in that if it is included it tells how many primes to multiply together rather than n but I bet this will make it tricky for me to code.

primorial2(n)=primorial(prime(n))

[science_man_88]

Originally Posted by CRGreathouse

primorial2(n)=primorial(prime(n))

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.

[R. Gerbicz]

Originally Posted by CRGreathouse

Here's a slightly-unfair comparison:

>primorial(10^6);
time = 220 ms.
>Primorial(10^6);
time = 1,830 ms.

(about 8x speedup)

Your code is still slow, a much faster (I don't know an asymptotically faster way to compute it):

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(len-mid,i,A[i+mid])))
Primorial2(n)=prodtree(primes(primepi(n)))

comparison:

? default(primelimit,10^7)
%117 = 10000000
? gettime();x=Primorial2(10^7);gettime()
%118 = 8580
? gettime();y=primorial(10^7);gettime()
%119 = 68219
? x-y
%120 = 0

[CRGreathouse]

Originally Posted by science_man_88

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.

That would be one approach. The other would be to warn the user that you don't have enough primes and to tell them how many you'd need to do the calculation.

primorialNaive(x)=my(pr=1,n=2);while(n<=x, pr *= n; n = nextprime(n+1)); pr

Time for primorialNaive(1e7): 320 seconds
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.

noop(x)=forprime(p=2,x,)
noopNaive(x)=my(n=2);while(n<=x, n = nextprime(n+1))

Going to 1e7, the first takes 20 milliseconds while the second takes 3.2 seconds (160 times faster).

[CRGreathouse]

Originally Posted by R. Gerbicz

Your code is still slow, a much faster (I don't know an asymptotically faster way to compute it):

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(len-mid,i,A[i+mid])))
Primorial2(n)=prodtree(primes(primepi(n)))

That version burst my stack! It's certainly a lot more memory-intensive. (It would be better with the library -- no need to duplicate the vector in that case.)

But yes, I get ~7 times faster with your version than with mine when I increase the stack appropriately.