Dear GAP Forum:

One example I've been trying to illustrate in my undergraduate
algebra class is the following:

 Let
     s(n+4) = s(n) + s(n+1) mod 2,
and s(0)=s(3)=1, s(1)=s(2)=0. This period 15 linear feedback 
shift register can be described using 
(a) the coefficients of x^i in the power series expansion
of 1/(1+x^3+x^4) (mod 2),
(b) the 1st entry of the vector A^i*k^t, where k=(1,0,0,1)
and A = [[0,1,0,0],[0,0,1,0],[0,0,0,1],[1,1,0,0]],
(c) s(i) = c1*r1^i+c2*r2^i+c3*r3^i+c4*r4^i, where 
r1, ..., r4 are the roots of x^4+x+1=0 (mod 2) in
GF(16) and c1, ..., c4 are constants (also in GF(16)).

GAP does all these, but (a) is the hardest. My code
for (a) is enclosed below. It is so slow that it is almost
impractical for classroom use (try computing 45 terms of the
power series for 1/(1+x^3+x^4) and GAP might crash from
lack of memory!). I had to extend several GAP functions
to rationals, perhaps I didn't do something right.

QUESTION: Does anyone have any suggestions for speeding up my code?

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Code for (a):

lfsr:=function(key,a,n,p)
local i,k;
k:=Size(key);
if (n<k and p>0) then 
 for i in [1..k] do
   if n=i-1 then return (key[i] mod p); fi;
 od;
fi;
if n<k then 
 for i in [1..k] do
   if n=i-1 then return key[i]; fi; 
 od;
fi;
if n>=k and p>0 then  
     return Sum([1..k],i->a[i]*lfsr(key,a,n-i,p)) mod p; 
fi;
if n>=k then 
     return Sum([1..k],i->a[i]*lfsr(key,a,n-i,p)); 
fi;
end;

## An example (Fibonacci):
## key:=[1,1];a:=[1,1];p:=0;
## L:=[];
## for i in [0..10] do Add(L,lfsr(key,a,i,p)); od;
## L; 

## Another example (the one mentioned above):
## L:=[];
## key:=[1,0,0,1];a:=[0,0,1,1];p:=2;
## for i in [0..10] do Add(L,lfsr(key,a,i,p)); od;
## L;

derivative:=function(r)
#extends Derivative to rational fcns
local f,g;
if not IsRationalFunction(r) then Print("wrong type\n"); return 0; fi;
f:=NumeratorOfRationalFunction(r);
g:=DenominatorOfRationalFunction(r);
return (Derivative(f)*g-f*Derivative(g))/g^2;
end;

## An example:
## x:=Indeterminate(Rationals,1);
## R:=PolynomialRing(Rationals,["x"]);
## f:=3+x^2+x^5;
## g:=1+x^2;
## derivative(f/g);
#######  Returns (-4*x+5*x^4+3*x^6)/(1+2*x^2+x^4)
#######  which is correct.

higher_derivative:=function(r,k)
## returns k-th der of rat fcn r
if k=0 then return r; fi;
if k=1 then return derivative(r); fi;
if k>1 then return derivative(higher_derivative(r,k-1)); fi;
end;

evaluate:=function(r,a)
#extends Value to rational fcns
#returns value of rat fcn r at a
local f,g;
f:=NumeratorOfRationalFunction(r);
g:=DenominatorOfRationalFunction(r);
return Value(f,a)/Value(g,a);
end;

series:=function(r,a,n)
##computes Taylor series of rat fcn r about x=a
local s,f,g,i;
s:=evaluate(r,a);
for i in [1..n] do
s:=s+evaluate(higher_derivative(r,i),a)*x^i/Factorial(i);
od;
return s;
end;

## An example:
## series(1/(1-x),0,10);  
######### returns 1+x+x^2+x^3+x^4+x^5+x^6+x^7+x^8+x^9+x^10
## series(1/(1+x^3+x^4),0,15);
######### returns 
######### 1-x^3-x^4+x^6+2*x^7+x^8-x^9-3*x^10-3*x^11+4*x^13+6*x^14+3*x^15
## series(1/(1+x^3+x^4),0,15) mod 2;
######### returns 1+x^3+x^4+x^6+x^8+x^9+x^10+x^11+x^15
#########   each taking about about 35 seconds on a 200 Mz celeron pc
under linux
## series(1/(1+x^3+x^4),0,45) mod 2;
#########  ran out of memory after several hours

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

--
Prof David Joyner, Mathematics Department
U. S. Naval Academy, Annapolis, MD 21402
phone: (410) 293-6738