divert(-1)
############################################################
#
# Macro definitions for m4
#
############################################################

# Synopsis: def_sum(node, i, start, end, expr)
# node := sum(i=start..end, expr)
#
# Synopsis: def_product(node, i, start, end, expr)
# node := prod(i=start..end, expr)

define(`def_sum', `
	__$1[$2] = 0 if $2<$3;
	__$1[$2] = __$1[$2-1] + $5;
	$1 = __$1[$4];
')

define(`def_product', `
	__$1[$2] = 1 if $2<$3;
	__$1[$2] = __$1[$2-1].($5);
	$1 = __$1[$4];
')

# Synopsis: def_sumx(node, i, start, end, expr, x)
# This function is best described on an example.
# If you say
#	def_sumx(`foo', i, 1, 4, `i+10', `i+100')
# it defines
#	foo := 11 + 101.(12 + 102.(13 + 103.14))
#           == 11 + 101.12 + 101.102.13 + 101.102.103.14
# When x doesn't depend on i, it's equivalent to the evaluation of a
# polynomial on x using Horner's method.

define(`def_sumx', `
	__$1[$2] = 0 if $2>$4;
	__$1[$2] = $5 + ($6).__$1[$2+1];
	$1 = __$1[$3];
')

# Synopsis: peel_polynomial(node, i, expr, lit)
# Meaning: node[i] := coefficient of lit^i in expr

define(`peel_polynomial', `
	$1[$2] = 0 if $2<0;
	$1[$2] = (__$1[$2],$4->0);
	__$1[0] = $3;
	__$1[$2] = (__$1[$2-1]-$1[$2-1])/($4);
')

# Synopsis: inverse_series(node2, node1)
# Meaning: the series sum(i=0..infty, node1[i].X^i) is inverted, and the
# coefficients of the inverse series are put into the node2[i].

define(`inverse_series', `
	$1[0] = 1/$2[0];
	__$1[i,-1] = 0;
	__$1[i,j] = __$1[i,j-1] + $1[j].$2[i-j];
	$1[i] = -__$1[i,i-1]/$2[0];
')

# Synopsis: additive_moebius_inversion(node, i, expr)
# Synopsis: multiplicative_moebius_inversion(node, i, expr)
# See the examples below.

define(`additive_moebius_inversion', `
	$1[$2] = $3 - __$1[$2,$2/2];
	__$1[i,0] = 0;
	__$1[i,j] = __$1[i,j-1] if i%j;
	__$1[i,j] = __$1[i,j-1] + $1[j];
')

define(`multiplicative_moebius_inversion', `
	$1[$2] = ($3) / __$1[$2,$2/2];
	__$1[i,0] = 1;
	__$1[i,j] = __$1[i,j-1] if i%j;
	__$1[i,j] = __$1[i,j-1] * $1[j];
')

### Exponential of a series
#
#define(`exp_integral', `
#	$1[0] = 1;
#	$1[k] = $1[k,k-1] / k;
#	$1[k,-1] = 0;
#	$1[k,l] = $1[k,l-1] + $2[k-1-l].$1[l] ifelse($#,2,,`.$3');
#')
#
#define(`exp_series', `
#	$1[0] = 1;
#	$1[k] = $1[k,k-1] / k;
#	$1[k,-1] = 0;
#	$1[k,l] = $1[k,l-1] + (k-l).$2[k-l].$1[l] ifelse($#,2,,`.$3');
#')

#
# End of macro definitions
#
divert(0)dnl
############################################################
#
# Rules for induce
#
############################################################

# Factorials
fact[0] = 1;
fact[k] = k.fact[k-1];

# You can't define ff[k] = fact[fact[k]], so how would you compute ff[k] ?
# (problem posed by maltey@enst-bretagne.fr); easy! just define the
# auxiliary nodes ff[k,l] = fact[fact[k].l], then it is plain that you
# have the following relations. NOTE: even ff[8] requires plenty of
# CPU time. You cannot compute ff[13] or more, because fact[13] >= 2^31,
# and the index will wrap around.
#
# Note also that such cleverness is no longer needed since induce(1) supports
# external functions; you just have to write fact[fact(k)] or fact(fact[k])
# or fact(fact(k)).

ffact[k] = ffact[k,1];
ffact[k,l] = ffact[k-1,k.l] if k>0;
ffact[0,l] = Fact[l];

# The following definition is probably more efficient if you have only one
# huge factorial to compute, because it tries to multiply numbers which
# have roughly the same size. We define
# Fact[n,s] = s.(s+1)...(s+n-1). Then the following relations hold:

Fact[n] = Fact[n,1];
Fact[n,s] = Fact[n/2,s].Fact[n-n/2,s+n/2] if n>100;
Fact[n,s] = s.(s+1).(s+2).(s+3).(s+4).(s+5).(s+6).(s+7).(s+8)
	.(s+9).Fact[n-10,s+10] if n>10;
Fact[n,s] = s.Fact[n-1,s+1] if n>1;
Fact[1,s] = s;
Fact[0,s] = 1;

# It is a syntax error to put an index (or even a constant integral
# expression) in an exponent. It's probably a bad idea anyway, but
# here's a workaround for (-1)^k.
#
# alt[i] = (-1)^k

alt[i] = -1	if i % 2;
alt[i] =  1;	# Otherwise.

# Two ways to compute approximations of the constant `e'. Compare them
# to see if they are equally sensitive to rounding errors.
# In both cases we compute 1+...+1/n!, but with Horner's method in the first
# case. This block of code is skipped if induce(1) runs in "exact" mode,
# for obvious reasons.

ifdef(`__PRECISION__', `
	def_sumx(`approx_e1[n]', k, 0, n, `1', `1/(k+1)')
	approx_e2[0] = 1;
	approx_e2[n] = approx_e2[n-1] + 1/fact[n];
	e1_e2[n] = approx_e1[n] - approx_e2[n];
')

# Let p[n] be the number of partitions of the integer n. These numbers
# verify the following recurrence relation (for n>0):
#
# p[n] - p[n-1] - p[n-2] + p[n-5] + ...
#	+ (-1)^k.p[n-k(3k-1)/2] + (-1)^k.p[n-k(3k+1)/2] + ... = 0.
#
# The only problem here is to stop the sum when the indices are negative.
part[n] = part[n,1]	if n > 0;
part[0] = 1;
part[n] = 0;		# When n is negative.

part[n,k] = 0		if n < k.(3.k-1)/2;
part[n,k] = part[n-k.(3.k-1)/2] + part[n-k.(3.k+1)/2] - part[n,k+1];

# The (in)famous Fibonacci numbers. Induce is perfectly adapted to this
# kind of definitions.
fibo[n] = fibo[n-1] + fibo[n-2]		if n > 1;
fibo[0] = 0;
fibo[1] = 1;
fibo[n] = fibo[n+2] - fibo[n+1];	# when n is negative.

# This is mainly for fun. If Induce can't compute fun[n] for some positive n,
# publish your result! (Alas, 32-bit integers can easily overflow.)
fun[1] = 0;
fun[n] = 1 + fun[n/2]		if n%2 = 0;
fun[n] = 1 + fun[(3.n+1)/2];	# when n is odd.

# Euler's totient function phi[n] and the Moebius function mu[n] can both
# be obtained with the Moebius inversion formula.
additive_moebius_inversion(phi, n, `n')
additive_moebius_inversion(mu, n, `!(n-1)')

# Inverse a series twice, and compare to the original series
s1[k] = 10 / (k+1);
inverse_series(s2,s1)
inverse_series(s3,s2)
s3error[k] = s3[k] - s1[k];

# Binomial coefficients
binom[n,p] = 0		if (p<0)+(p>n);
binom[n,0] = 1;
binom[n,n] = 1;
binom[n,p] = binom[n-1,p-1]+binom[n-1,p];
Binom[n,p] = fact[n] / fact[p] / fact[n-p];

# Some calculations on the iterates of P(a,b): z->az+bz2+z3
# First define z_exp[k]=z^k
z_exp[0] = 1;
z_exp[k] = z.z_exp[k-1];

# We suppose here that `a' is a n-th primitive root of unity. Let Iter[k,n]
# be the k-th iterate of P(a,b), reduced modulo the cyclotomic polynomial,
# computed up to z^(n+1). We know that Iter[n,n] will have the form
# Iter[n,n] = z + Cab[n].z^(n+1), and we want to compute Cab[n].

Iter[0,n] = z;
Iter[i,n] = (Iter[i-1,n], z -> a.z+b.z^2+z^3, z_exp[n+2] -> 0,
	cmain[n] >> cmain[n]-cyclo[n]);
Iter[n] = Iter[n,n];
Cab[n] = (Iter[n] - z) / z_exp[n+1];

CabDisc[n] = (disc(Cab[n],b), cmain[n] >> cmain[n]-cyclo[n]);

CabDisc_e[n,x] = disc((Cab[n],a->x),b);
CabDisc_e[n,k,0] = CabDisc_e[n,k];
CabDisc_e[n,x,k] = (CabDisc_e[n,x+1,k-1] - CabDisc_e[n,x,k-1]) / k;
poch[0] = 1;
poch[k] = poch[k-1].(a-k+1);
CabDisc_i[n,-1] = 0;
CabDisc_i[n,k] = CabDisc_i[n,k-1] + CabDisc_e[n,0,k].poch[k];
CabDisc2[n,k] = (CabDisc_i[n,k], cmain[n] >> cmain[n]-cyclo[n]);
CabDisc2[n] = CabDisc2[n,(2.n-1).phi(n)];

# Cyclotomic polynomials (and their leading coefficient)
multiplicative_moebius_inversion(`cyclo', n, `(z_exp[n]-1,z->a)')
multiplicative_moebius_inversion(`cmain', n, `(z_exp[n],z->a)')

# Let phi(x) = x + O(1) such that phi(x^2) = phi(x)^2 + c.
# phi(x) = x + cphi[1]/x + cphi[2]/x^2 + ...
# In real mode we choose c=1/4.

ifdef(`__PRECISION__',`c=1/4;')
cphi[k] = 0	if k%2 = 0;
cphi[1] = -c/2;
cphi[k] = cphi[k/2]/2 - cphi[k-1,k-1]	if k%2;
cphi[k] = - cphi[k-1,k-1];		# Otherwise.
cphi[k,0] = 0;
cphi[k,l] = cphi[k,l-1] + cphi[l].cphi[k-l];

# Let arrg[m,n] be the number of colourings of a m.n rectangle with three
# colours, such that two adjacent cells haven't the same color.
# This number can be computed as follows.
arrg[n,0,j] = 1;
arrg[n,i,0] = arrg[n,i-1,n];
arrg[n,i,j] = (
	(arrg[n,i,j-1], c1[i-1,j]->0, c1[i,j-1]->0) . c1[i,j] +
	(arrg[n,i,j-1], c2[i-1,j]->0, c2[i,j-1]->0) . c2[i,j] +
	(arrg[n,i,j-1], c3[i-1,j]->0, c3[i,j-1]->0) . c3[i,j],
		c1[i-1,j] -> 1, c2[i-1,j] -> 1, c3[i-1,j] -> 1  );
arrg_[m,n,0] = arrg[n,m,n];
arrg_[m,n,i] = (arrg_[m,n,i-1], c1[m,i]->1, c2[m,i]->1, c3[m,i]->1);
arrg[m,n] = arrg_[m,n,n];
arrg[n] = arrg[n,n];

# You can handle matrices too, with the appropriate notations. Let's check
# the Cauchy-Hamilton theorem on small matrices, for instance.
# We're cheating slightly here, because the difficult calculations are done
# with the external functions disc() and mxprod(). The point here is that you
# can use induce as a "shell" to coordinate many small processes.

identity[0] = 0;
identity[s] = identity[s-1] + __e[s,s];
def_sum(`generic[s,i]', j, 1, s, `a[i,j].__e[i,j]')
def_sum(`generic[s]',   i, 1, s, `generic[s,i]')
# Generic matrices are too cpu-intensive? Define the symbol `HILBERT'
# if you simply want to check Cayley-Hamilton on Hilbert matrices.
ifdef(`HILBERT',`a[n,p] = 1/(n+p-1);')

# Let pcar[n,k] be the coefficient of X^k in pcar[n]
pcar[n] = det(X.identity[n]-generic[n],n);
peel_polynomial(`pcar[n]', k, `pcar[n]', X)
# Now we can prove the Cayley-Hamilton theorem
c_ham[n,n] = pcar[n,n].identity[n];
c_ham[n,k] = pcar[n,k].identity[n] + mxprod(generic[n],c_ham[n,k+1]);
c_ham[n] = c_ham[n,0];  # Should always be zero.

# End of file