User:Charles R Greathouse IV/Pari

Formatting

See User:Charles R Greathouse IV/Format for a script to format b-files.

Transforms

Christian G. Bower has a script file with code for the basic transformations used in the OEIS.

Recurrences

Here's some really ugly code that works out a recurrence relation, given the first few terms of the sequence. Usage: rec(1,1,2,3,5,8) or rec(v), where v is a vector.

rec(v[..])=if(#v==1&&type(v[1])=="t_VEC",v=v[1]);findrec(v)
addhelp(rec, "rec(v): Tries to find a homogeneous linear recurrence relation with constant coefficients to fit the vector v.");

findrec(v:vec, verbose:bool=1)={
	my(c,d = (#v - 1) >> 1, firstNonzero = 0);
	if (#v < 3,
		warning(Str("findrec: Not enough information to determine if the sequence is a recurrence relation: matrix is underdetermined. Give more terms and try again."));
		return
	);
	forstep(i=d,1,-1,
		if(v[i] != 0, firstNonzero = i)
	);
	if(firstNonzero == 0,
		if(verbose,
			print1("Recurrence relation is a(n) = 0.");
		);
		return([0]~);
	);
	for (i=firstNonzero,d,
		c = findrecd(v,i,verbose);
		if(c, return(c))
	);
	if(verbose,print("Cannot be described by a homogeneous linear recurrence relation with "d" or fewer coefficients."));
	0
};
addhelp(findrec, "findrec(v, {verbose=1}): Tries to find a homogeneous linear recurrence relation with constant coefficients to fit the vector v.");


findrecd(v:vec, d:int, verbose:bool=1)={
	my(M,c);
	if (#v < 2*d,
		warning(Str("findrec: Not enough information to determine if the sequence is a "d"-coefficient recurrence relation; matrix is underdetermined. Give more terms or reduce d and try again."));
		return
	);
	M = matrix(d,d,x,y,v[d+x-y]);
	\\print(M" * c = "vector(d,i,v[d+i])~);
	if(matdet(M) == 0, return(0));	\\ Non-invertible matrix, no solutions for this d
	c = matsolve(M,vector(d,i,v[d+i])~);
	for(n=2*d+1,#v,
		if(v[n] != sum(i=1,d,v[n-i] * c[i]),return(0))
	);
	if(verbose,
		my(init=1,s);
		print1("Recurrence relation is a(n) = ");
		for(i=1,#c,
			if(c[i] == 0, next);
			if(init,
				s = initial(c[i], Str("a(n-", i, ")"));
				init = 0
			,
				s = Str(s, medial(c[i], Str("a(n-", i, ")")))
			)
		);
		print(s".");
		if((vecmax(c) == 1 && vecmin(c) == 0 && vecsum(c) == 1) || c == [1]~,
			print("Sequence has period "#c".");		
		,
			my(g=0);
			for(i=1,#c,
				if(c[i] != 0, g = gcd(g, i))
			);
			if (g > 1,
				my(gvec = vector(#c/g, i, c[i*g]),s,init=1);
				for(i=1,#gvec,
					if(gvec[i] == 0, next);
					if(init,
						s = initial(gvec[i], Str("a(n - ", i, ")"));
						init = 0
					,
						s = Str(s, medial(gvec[i], Str("a(n - ", i, ")")))
					)
				);
				print("Can be though of as "g" interlocking sequences, each of the form a(n) = "s".")
			)
		);
		print1("<a href=\"/index/Rec#order_", if(#c<10,Str("0",#c),#c), "\">Index to sequences with linear recurrences with constant coefficients</a>, signature ("c[1]);
		for(i=2,#c,print1(","c[i]));
		print(").");
		print((#v-2*d)" d.f.")
	);
	c
};
addhelp(findrecd, "findrecd(v, d, {verbose=1}): Helper function for findrec. Tries to find a d-coefficient homogeneous linear recurrence relation with constant coefficients to fit the vector v.");

initial(n:int, s:str)={
	if (n == 0, return(""));
	if (n == -1, return(Str("-", s)));
	if (n == 1, return(s));
	Str(n, s)
};


medial(n:int, s:str)={
	if (n == 0, return(""));
	if (n == -1, return(Str(" - ", s)));
	if (n == 1, return(Str(" + ", s)));
	Str(if (n < 0, " - ", " + "), abs(n), s)
};


if(warning=='warning,warning(s:str)=print("Warning: "s));

See also

• Sequence Tools which includes other PARI/GP tools