# last changed: 
#   200204  : improved comments, "system independence"
#   200107  : added dimensions 17 to 33
#   20000711: default pn corrected from g+1 to 2g-1
#   20000712: individual pn for each place
#   20000719: corrected fast (=Newton) expansion routine mx (now called xn)

Read("prototypes");		# variable declarations
Read("ffdata");			# loading function field data

#### setbase 
## Arguments : 
##    b - a prime power
## Result    : 
##    initializes the rational function field over F_b
####

setbase := function(btemp)
  b  := btemp;                                  
  AlffInit(FF(b));
end;
setbase(2);

## just some convenient abbreviations

div := AlffDivisor;
lbas:= AlffDivisorLBasis;
ldim:= AlffDivisorLDim;
v   := AlffEltValuation;

####  sortedP  
## Arguments : 
##    F - a function field
## Result    : 
##    generates a "sorted" list of deg 1 places
## Comments  :
##    Difference to AlffPlaces : for each function field
##    generated by some fixed equation always the same
##    list is produced; with AlffPlaces it may vary.      
##    Implemented sorting here is lex according to screen output in KASH
####

sortedP:=function(F)
 local P,Psortf;
 P:=AlffPlaces(F,1);
 Psortf:=function(x,y) 
  if SPrint(AlffPlaceMin(x))<SPrint(AlffPlaceMin(y)) then
   return true;
 elif SPrint(AlffPlaceMin(x))=SPrint(AlffPlaceMin(y)) then
  return Reversed(SPrint(AlffPlaceIdeal(x))) < 
         Reversed(SPrint(AlffPlaceIdeal(y)));
 else return false;
 fi;
 end;
 Sort(P,Psortf);
 return P;
end;

#### divpol    ( division of polynomial truncations of series )
## Arguments :
##    pol? - dividend and divisor
## Result :
##    The series expansion of pol1/pol2
## Comments :
####

divpol:=function(pol1,pol2)
 local dg,pol,frac;

# first find the highest factor T^dg contained in pol2 and div by it

# dg   := Position(Reversed(PolyToList(pol2)),ffe[1]^0)-1; # for b=2
 dg := PolyDeg(pol2);
 dg := dg - PolyDeg(T^dg*Eval(pol2,1/T));
 frac := pol2/T^dg;

 pol  := PolyQuotRem(pol1,frac);

# for b = 2 : 
# emulate the inversion 1 / frac by the
# geometric series 1/frac ca. = 1+(frac+1)+(frac+1)^2+...+(frac+1)^255 = 
# = (1+(frac+1))(1+(frac+1)^2)...(1+(frac+1)^128) =
# = (frac)(1+(frac^2+1^2))...(1+(frac^128+1^128)) = frac^(1+2+...+128)
# for b > 2 : using (1+frac+...+frac^(b-1) )^(b^5-1) analogously

# if b=2 then
#  frac:= frac^255;
# else
#  frac := PolyQuotRem( frac^b-2, frac-2)[1]*frac^(b*(b^5-1)/(b-1));
# fi;

# or, simply (and safer ?) :

# frac := Sum([0..499],i->(1-frac)^i);

# faster, by geometric series formula:

 if PolyDeg(frac)=0
 then
  frac := 1/frac;
 else
  frac := (-frac)^(b-1)*Product([1..5],i->(2-frac^(b^i)));
 fi;

 frac := pol[2]*frac;
 pol  := pol[1]+frac;

# apply the remaining division by T^dg (disregarding fractional terms 
#  because there should be none in the applications below ! )

# attention : polyquotrem does not allow the second arg to be 1, so:

 if dg=0
 then
  return pol;
 else
  return PolyQuotRem(pol,T^dg)[1];
 fi;
end;

#### albas    ( "*a*scending" *L-bas*is )
## Arguments:
##   pl - a place
##   c0,c1 - range of coefficients for pl in a divisor
## Result:
##   a list of two things:
##   1) a basis of L(D+c1*pl), ordered such that the i-th vector 
##   is a basis vector of L(D+n_i*pl), where n_i >= c0+i is a 
##   pole number of pl at D (i.e. the dimension increases)
##   2) the n_i
####

albas := function(pl,c0,c1)
  local nc,bas,i,vi,tbas;
  if pl=p[1] then        # since D=2*g*p[1], we need an offset if pl=p[1]
   vi := 2*g;
  else
   vi := 0;
  fi;
  nc  := [];

# compute the "ground" basis

  bas := lbas(D+c0*pl);

# the first basis is sorted by valuation at pl
# (although these vectors are never really needed, I think)

  Sort(bas, function(a,b) return v(pl,one*a) > v(pl,one*b); end );

# in the loop  we compute a fresh L-basis of the increased
# divisor D+i*pl ,  check wether or not dimension increases and if so
# add the vector and its valuation to the output lists

  for i in [c0+1..c1] do
   tbas := lbas(D+i*pl);
   if Length(tbas) > Length(bas) then
    Add( bas, First( tbas, e->v(pl,e) = -(i+vi) ));
    Add( nc, i);
   fi;
  od;
  return [bas, nc];
 end;

#### xp         (eXPansion)
## Arguments :
##    e - the function field element to be expanded
##    n - number of terms in the expansion
## Result : 
##    List of coefficients of expansion of e at p0 up to T^n
##    starting with 1 
## Comments :
##    A bit Lisp-ishly obfuscated, I'm afraid...
##    Nevertheless, the core is the use of the (imperfect and undocumented)
##    KASH-internal function EltExpansion. The result of this function is
##    a list of '[coefficient, exponent]' entries. This is converted to a 
##    "pure" coefficient list. (First converted to a
##    polynomial in 'T', then 'T^(n+1)' is added to gain the trailing zeroes,
##    then conversion back to a list and reversed.) 
####

xp  := function(e,n) 
     local el, eone;
     eone := AlffElt(eq,Concatenation([1],0*[2..fdeg]));
        return Reversed( PolyToList( T^(n+1)+
          Sum(List(
           AlffPlaceLocals.EltExpansion(p0,T,AlffEltMove(eone*e,o),n,false),
          i->i[1]*T^i[2])))){[1..n]};
       end;

#### Further auxiliary (debug) functions

showmats   := function() Exec("cat kashout"); end;
testmats   := function() Exec(SPrint("./tcalc kashout ",b)); end;

#######################################################################

Print("\nUsage: calc_c(m,s), produces (g(s),m,s-1)-seq gen. matrices\n"); 

############ calc_c   ( calculate the generating matrices )
## Arguments :
##    m - size of matrices (= log number of elements in generated point set)
##    s - number (+1) of matrices (= dimensions of generated point set)
## Result :
##    Generating matrices in GF(b)
############

calc_c_par := function(F,b,m,s)        #  m : size of matrices, s : no. of mat.

local pne,f8rev;


ffe := FFElts(FF(b));                       # list of fin.field elements
cnv := l->List( l, e->Position(ffe,e)-1 );  # biject. integers to fin.fld.elts
                                            # beware: NOT eq `natural bij'

o:=AlffOrderMaxFinite(F);

Print("\n Please hold on, while bases are calculated ...\n");

fdeg:= AlffDeg(F);
eq  := AlffOrderEqFinite(F);

## Y   : an element such that F = k(T,Y)
## one : multiplicative unit of F ( not of k !) 
## p   : a "sorted" list of places of degree one

Y   := AlffElt(eq,Concatenation([0,1],0*[3..fdeg])); 
one := AlffElt(o,Concatenation([1],0*[2..fdeg]));

p   := sortedP(F); 

##  p0 : a distinguished nonramified place above T (the P_\infty of the paper) 
##  p0ind : its index in p

p0ind    := PositionProperty(p, pl->
             (pl in AlffPlaceSplit(F,T)) and (AlffPlaceRam(pl) = 1) );
p0       := p[p0ind];  

## D : a prime divisor of degree twice the genus g, and with support {p[1]}
##     not containing p0. 

p[p0ind] := p[1];
p        := p{[2..s+1]};        
g   := AlffGenus(F);
D   := 2*g*p[1];

c   := [];                   # declaration of generating matrices 
cl  := []; 

##  the w and n are important data for the Niederreiter-Xing modified
##  local expansions that are used in this technique

wn := albas(p0, -(2*g+1), 0);
w  := Reversed(wn[1]);
n  := Reversed(-wn[2])+1;
 

######## xn : produces AlffEltGenY(F) by Newton iteration
#          (but we need AlffElt(eq,[0,1,0,...,0]), so xn isn't used )

xn:=function(prec,init)
 local te,tl1,tl2;
 te  := init;
 tl1 := init;
 repeat
  te  := te + Eval(ffdata[s-3],te) / Eval(PolyDeriv(ffdata[s-3]),te);
  tl2 := tl1;
  tl1 := divpol(QfeNum(te), QfeDen(te)) mod T^prec;
 until tl1=tl2;
 return tl1;
end;

## YP : list of expansions of powers of Y, = basis of finite equation order

# First attempt: (wrong, see xn descr., why)
#    YP:=xn(2*(m+g+1), Poly(kT,Reversed(xp(Y,1+v(p0,Y)))) ); 
# Second try (that worked!):

YP:=xp(Y,2*(m+g+2));                                 # Expansion of Y
YP:=Sum(List([1..2*(m+g+2)],i->YP[i]*T^(i-1)));      # Into poly
YP:=List([0..fdeg-1],i->YP^i mod T^(2*(m+g+2)) );    # Powers  

#### mx   ( "my" expansion routine )
## Arguments :
##    e - a function field element
##    n - degree of the resulting polynomial
## Results :
##    again, the local expansion
## Comments :
##    Lisp strikes again ... Similar processing steps as in the xp routine.
##    This routine circumvents the (imperfect?) internal expansion routine by
##    regarding the func.fld.elt.s as vectors over the finite equation order.
####

mx:=function(e,n) 
 local el, eone;
 eone := AlffElt(o,Concatenation([1],0*[2..fdeg]));
 el:=AlffEltToList(AlffEltMove(eone*e,eq));
 return Reversed(PolyToList( T^(n+1)+ 
   (divpol( Sum(List([1..fdeg], i->el[1][i]*YP[(i-1)+1])), el[2]) 
   mod T^(n+1)))){[1..n]};
end;

#### mxp   ( "my" expansion routine )
## Arguments :
##    e - a function field element
##    n - precision to expand to = number of terms
## Results :
##    again, an expansion routine
## Comments :
##    
####

mxp:=function(e,n) 
 local el, eone;
 eone := AlffElt(o,Concatenation([1],0*[2..fdeg]));
 el:=AlffEltToList(AlffEltMove(eone*e,eq));
 #local el;
 #el:=AlffEltToList(AlffEltMove(e,eq));
 return 
   (divpol( Sum(List([1..fdeg], i->el[1][i]*YP[(i-1)+1])), el[2]) 
   mod T^(n+1));
end;


#### nxp   ( Niederreiter-Xing expansion )
## Arguments :
##    e - a function field element
##    prec - precision to expand to = number of terms
## Result :
##    Niederreiter-Xing modified expansion
## Comments :
####

nxp := function(e,prec)
  local txp,te,i,xpp;    # txp - temp. expansion, te - temp. element
  xpp:=xp;               # xp is slower, mx "safer" (more likely to be correct)

# first we take the original  local expansion

  txp := xpp(e,prec);

# then we succesively replace the terms with exponent n[i] of the series 
# expansion by w[i] 

  te  := e - txp[n[1]]*w[1];
  for i in [2..g+1] do
   txp{[n[i-1]+1..prec]} := xpp(te,prec){[n[i-1]+1..prec]}; 
   te := te - txp[n[i]]*w[i];
  od;
  txp{[n[g+1]+1..prec]} := xpp(te,prec){[n[g+1]+1..prec]};
  return txp;
 end;


#### pn, pne ( place "number" entry )
## Arguments :
##    pne : pl - a place
## Result :
##    The second pole number of pl, i.e. the first integer i such
##    that dim (i*pl) > 1. 
##    pn returns a list of pne for all elts. of p.
## Comment :
##    we are interested in L( pne(pl)*pl ), because the nontrivial
##    basis element of it can be used to construct arbitrarily
##    large L( D+i*pl ), if enough bases for fixed i are already known
####

pne:=function(pl) 
 local i; i:=1; 
 while ldim(i*pl)=1 do i:=i+1; od; 
 return i; end;
pn:=List(p, pne);

#### pt  ( (Place) L-space generators table )
## Arguments :
##
## Result :
##    a list of generators of L( D+i*p[j] ) 
####

pt := List( [1..s], pl->Concatenation( 
  albas(p[pl],0,pn[pl])[1]{[g+2..g+pn[pl]+1]},
  [ First(lbas(pn[pl]*p[pl]),ee->v(p[pl],one*ee)=-pn[pl]) ]) ); 
PT := pt;
pt := List(pt,eee->List(eee,ee->mxp(ee,m+g+1+5)));

ind:=function(pn,i)
  if i mod (pn)=0 then
   return pn;
  else
   return i mod pn;
  fi;
 end;

make_c  := 
 function(m)
  local xind;
  xind:=[1..m+g+1];
  SetSubtract(xind,n);
  c:=List([1..s], ptn->
      List([1..m],          it->
   nxp( pt[ptn][pn[ptn]+1]^Floor((it-1)/(pn[ptn]))*
        pt[ptn][ind(pn[ptn],it)], m+g+1){xind} ));
  return c;
 end;

num_c := function(m)
  local i,j,bpow;
  make_c(m);
  bpow:=Reversed(List([1..m],i->2^(i-1)));
  for i in [1..s] do
   cl[i]:=[];
   for j in [1..m] do
    cl[i][j] := cnv(c[i][j])*bpow;
   od;
  od;
  cl[s+1] := Reversed(bpow); 
  return cl;
 end;

############# After all the above definitions, this is 
############# where the main routine ( num_c ) is called

Print(" All preliminaries done, now for the expansions...\n");

 cl := num_c(m);
 PrintTo("kashout",m," ",s);
 for i in [1..s] do
  AppendTo("kashout","\n");
  for j in [1..m] do
   AppendTo("kashout", cl[i][j]," ");
  od;
 od;
 AppendTo("kashout", "\n");
 PrintTo("kashout2", c);

Print(" Matrices constructed and written to 'kashout'\n");

# These lines are unix-specific, modify according to taste (or replace x@y 
# with your e-mailaddress )

#Exec("mail -s 'Matrizen fertig!!!' x@y < kashout");
#Exec("./t2calc kashout && mail -s 't2calc fertig!!!' x@y < tresults &");

# If you are aware how it might fail, use this line for batch processing 
# to kill the kash process

#Exec("kill `/sbin/pidof kash`");

end;

calc_c := function(m,s)
 AlffInit(FF(b));                            # initialize func.flds. over F_b
 initffdata();
 s  := s-1; 
 if s-3>Length(ffdata) then
  Print("\n Sorry, currently only dimensions up to ",Length(ffdata)+4,".\n");
 return;
 fi;
 F:=Alff( ffdata[s-3] );       # initialize the func.fld. with the equation
 calc_c_par(F,b,m,s);
 showmats();
 testmats();
end;

b3testcase := function() 
 b:=3;
 AlffInit(FF(b));
 F:=Alff(y^2-(T^3-T+1));
 calc_c_par(F,b,13,5);
 showmats();
 testmats();
end;

#calc_c(13,5); Exec("cat kashout");Exec("kill `/sbin/pidof kash`");