\\ This pari program implements multiple precision arithmetic
\\ using pari! This may seem silly but is a good way of introducing 
\\ pari programming *and* showing that it works.

\\ Multiple precision numbers are implemented as polynomials
\\ since pari has no other linked list type! The coefficient of 
\\ X^n is of course M^n in this representation.

\\ The base 
M=10000;

\\ A timer to keep track of number of operations
time=0;

\\ The pari function polcoeff implements taking the the ith digit
digit(u,i)=
{
	polcoeff(u,i);
}

\\ The pari function poldegree gives the size
size(u)=
{
	poldegree(u);
}

\\ We need to have a nice input and output
pinput(n)=
{
	local(t,w);
	w=0;t=0;
	while(n,
		w=w+(n%M)*X^t;
		t=t+1;
		n=n\M;
	);
	w;
}

poutput(u)=
{
	subst(u,X,M);
}

\\ Single precision operations
\\ You can ignore this part since it is the "hidden part"
\\ inside a real computer.
\r assembly.gp

\\ Now the fun begins.

r=0; \\ This will denote a global carry/borrow! We can always check here for an
     \\ overflow

addition(u,v)=
{
	local(m,i,x,w);
	m=max(size(u),size(v));

	r=0;w=0;
	for(i=0, m,
	  x = padd( digit(u,i), digit(v,i), r); \\ x contains both parts
	  w = w + x[1]*X^i; \\ The first part is the extra "digit"
	  r = x[2]; \\ The rest is the carry over
	);
	w; \\ Print the "answer" note that r may contain a carry!
}

subtraction(u,v)=
{
	local(m,i,x,w);
	m=max(size(u),size(v));

	r=0;w=0;
	for(i=0, m,
	   x = psub( digit(u,i), digit(v,i), r);
	   w = w + x[1]*X^i;
	   r = x[2]; \\ r denotes the borrow
	);
	w; \\ Note that the algorithm is the same as addition with
	   \\ padd replaced by psub. Is that surprising or not?!
}

\\ Easy so far. Now comes the first complication
\\ Compute a.v+t when a is a digit and t has more digits than v
mulbydigit(a,v,t)=
{
	local(m,i,x,y,z,w);
	m=max(size(v)+1,size(t))+1;

	w=0;z=[digit(t,0),0];y=[0,0];
	for(i=0, m,
	   x = pmul( a, digit(v,i) );
	   y = padd( z[1], x[1], y[2]); \\ y[2] is one carry over
	   w = w + y[1]*X^i;
	   z = padd( digit(t,i+1),x[2],z[2]); \\ z[2] is another
	);
	w;
}

\\ Now the "actual" multiplication is easy

multiplication(u,v)=
{
	local(p,i,z,w);
	p=size(u);

	w=0;z=0;
	for(i=0, p,
		z = mulbydigit( digit(u,i), v, z);
		w = w + digit(z,0)*X^i; 
		z = (z - digit(z,0))/X;
	);
	w + z*X^(p+1);
}

\\ Things start to get hairy. First of all let us dispense with
\\ short division
shortdiv(a,u)=
{
	local(p,i,x,w);
	p=size(u);

	x=[0,0];w=0;
	for(i=0, p,
		x = pdiv( a, x[2], digit(u,p-i));
		w = w + x[1]*X^(p-i);
	);
	[ w, x[2] ];
}

\\ Long division by assuming the leading digit is "good"
\\ to improve the guesses. here v may have one digit more than
\\ u bu that digit should be smaller than the leading
\\ digit of u.
guessdiv(u,v)=
{
	local(p,q,x,w);
	p=size(u);q=max(size(v),p+1);

	r = 0;
	w = pdiv( digit(u,p), digit(v,q), digit(v,q-1))[1];
	x = mulbydigit(w,u,0);
	x = subtraction(v,x*X^(q-p-1));
	if( r == 0,
	    [ w*X^(q-p-1), x ],
	    x = addition(u*X^(q-p-1),x);
	    if ( r == 0,
		[ (w-2)*X^(q-p-1), addition(u*X^(q-p-1),x) ],
		[ (w-1)*X^(q-p-1), x]
	    );
	);
}	 

\\ Long division at last
division(u,v)=
{
	local(p,m,x,d,w);
	m=size(v)-size(u);p=size(u);

	x = padd( digit(u,p), 0 , 1);
	if( x[2] == 0,
		d = pdiv(x[1], 1, 0)[1],
		d = 1);
	u = mulbydigit(d,u,0);
	v = mulbydigit(d,v,0);
	w = [ 0, v ];
	while(m >= 0,
		w = [ w[1], 0 ] + guessdiv(u,w[2]);
		m--;
	);
	w = [w[1], shortdiv(d,w[2])[1]]; \\ Since d divides w[2]
			\\ there is no remainder.
}