############################################################
#
#  Linear feedback shift registers mod p
#
# series, derivatives of rational functions, ...
#10-2001
#
#  A. Hulpke code improvements added
#  10-3-2001 
############################################################

 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:
## 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:
## 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)
 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)
#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 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
## series(1/(1+x^3+x^4),0,45) mod 2;
#########  ran out of memory after several hours


################################################################
#Alexander Hulpke - modified code, using caching of intermediate results

    Dolfsr:=function(key,a,n,p)
    local k,cache,lfsr,lfsrcache,offset;
    k:=Size(key);
    cache:=[];
    offset:=1;

    # check whether value known
    lfsrcache:=function(n)
      if not IsBound(cache[n+offset]) then
        cache[n+offset]:=lfsr(n);
      fi;
      return cache[n+offset];
    end;

    lfsr:=function(n)
    local i;
    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]*lfsrcache(n-i)) mod p; 
    fi;
    if n>=k then 
          return Sum([1..k],i->a[i]*lfsrcache(n-i)); 
    fi;
    end;


    return lfsrcache(n);
    end;

    L:=[];
    key:=[1,0,0,1];a:=[0,0,1,1];p:=2;
    for i in [0..1000] do Add(L,Dolfsr(key,a,i,p)); od;
    return L;

---------------------------------------------------------------------
# modified polynomial code


# derivative is declared as operation in the library, not as attribute
# (this -- as well as the missing Value method -- will be added in the next
# release.)
MyDerivative:=NewAttribute("MyDerivative",IsUnivariateRationalFunction);

InstallMethod(MyDerivative,"for univariate rational functions",
  true,[IsUnivariateRationalFunction],0,
function(r)
#extends Derivative to rational fcns
 local f,g;
 # this is now done automatically by the method installation
 #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);

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 MyDerivative(r); fi;
if k>1 then return higher_derivative(MyDerivative(r),k-1); fi;
end;

# replacement for `evaluate'
InstallOtherMethod(Value,"univariate rational functions",true,
  [IsUnivariateRationalFunction,IsObject],0,
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:=Value(r,a);
for i in [1..n] do
 s:=s+Value(higher_derivative(r,i),a)*x^i/Factorial(i);
od;
return s;
end;