# last changed: 
#   20000711: default pn corrected from g+1 to 2g-1
#             deleted computation of Weierstrass places (silly idea, that)
#   20000712: individual pn for each place
#   20000719: corrected fast (=Newton) expansion routine mx

Read("prototypes");

calc_c := function(m,s)        #  m : size of matrices, s : no. of mat.

local pne;

b  := 2;                                  # base restricted to 2 for now

AlffInit(FF(b));                            # initialize func.flds. over F_b
ffe := FFElts(FF(b));                       #  need those two, when writing
cnv := l->List( l, e->Position(ffe,e)-1 );  #  the mats to file
                                            # beware: NOT eq `natural bij'

tp := function(l) return Sum(List(l,i->T^i)); end;
y4 := function(c4,c2,c1,c0) 
        return y^4*tp(c4)+y^2*tp(c2)+y*tp(c1)+tp(c0); 
      end;
yp := function(yl,tpl) 
        return Sum(List([1..Length(yl)], i-> y^yl[i]*tp(tpl[i]))); 
      end;

ffdata:=[                        # some of the collected defining eq.s :
                                 # most of those are optimal (N(F)=N_2(g(F)))
 y^2+y+(T^3+T),                             # N(F) = 5, g = 1
 #(T^3+T+1)*(y^2+y)+(T^2+T),                 # 6,2  
 y^2+(T^3+T+1)*y+tp([5,4,3,1]),              # 6,2  
 y^3+(T^2+T+1)*y^2+(T^3+T^2)*y+(T^2+T),      # 7,3
# y4([1,0],[0],[1],[4,3]),                   # 8,4
 y4([0],[1,0],[3,1],[7,3]),                  # 8,4 # the same, w. integral eq.
 y4([0],[2,1,0],[2,1],[7,3]),                # 9,5   
 y4([10,8,6,4,0],[10,9,6,5,4,2,0],[9,8,5,2],[12,11,8,5]),          # 10,6
# y4([0],[20,19,18,17,14,12,11,10,8,6,5,2,0], # the same, with integral eq.
#    [29,28,24,22,20,18,16,14,13,10,9,8,5,2],
#    [42,41,40,39,38,36,35,33,32,31,27,26,24,19,17,16,14,9,8,5]),   # 10,6
 y4([8,4,0],[9,6,5,1,0],[9,8,6,5,4,1],[8,7,4,3]),                  # 11,8
 y4([12,10,8,6,4,2,0],[14,13,11,9,7,5,3,1,0],[1..14],[11,9,7,5]),  # 12,9
 yp([8,4,2,1,0],[[0],[3,2,0],[6,4,3,2],[6,4],[11,7]]),             # 15,13
 yp([8,4,2,1,0],[[0],[3,2,0],[6,4,3,2],[6,4],[11,7]]),             # 15,13
 yp([8,4,2,1,0],[[0],[3,2,0],[6,4,3,2],[6,4],[11,7]]),             # 15,13
 yp([0,1,2,4,8],[[7,15],[4,8],[2,3,4,5,6,8],[0,2,3,5,6],[0]]),     # 17,17
 yp([0,1,2,4,8],[[7,15],[4,8],[2,3,4,5,6,8],[0,2,3,5,6],[0]]),     # 17,17
 yp([0..21],[ [10],[7,9],[6,8,9,10],[3,4,5,7,8,10],                # 21,21
 [2,3,4,5,6,8,9],[1,3,4,5,9],[0,1],[0,1,4,7,9],[3,6,8,9,10], 
 [1,6,9,10],[0,1,2,4,8,9],[0,1,2,4,5,6,8,9,10],[7,10],
 [2,6,9,10],[9],[10],[],[9,10],[8,9,10],[8,9,10],[],[10] ]),
 yp([0..21],[ [10],[7,9],[6,8,9,10],[3,4,5,7,8,10],             # 21,21
 [2,3,4,5,6,8,9],[1,3,4,5,9],[0,1],[0,1,4,7,9],[3,6,8,9,10], 
 [1,6,9,10],[0,1,2,4,8,9],[0,1,2,4,5,6,8,9,10],[7,10],
 [2,6,9,10],[9],[10],[],[9,10],[8,9,10],[8,9,10],[],[10] ]),
 yp([0..21],[ [10],[7,9],[6,8,9,10],[3,4,5,7,8,10],             # 21,21
 [2,3,4,5,6,8,9],[1,3,4,5,9],[0,1],[0,1,4,7,9],[3,6,8,9,10], 
 [1,6,9,10],[0,1,2,4,8,9],[0,1,2,4,5,6,8,9,10],[7,10],
 [2,6,9,10],[9],[10],[],[9,10],[8,9,10],[8,9,10],[],[10] ]),
 yp([0..21],[ [10],[7,9],[6,8,9,10],[3,4,5,7,8,10],             # 21,21
 [2,3,4,5,6,8,9],[1,3,4,5,9],[0,1],[0,1,4,7,9],[3,6,8,9,10], 
 [1,6,9,10],[0,1,2,4,8,9],[0,1,2,4,5,6,8,9,10],[7,10],
 [2,6,9,10],[9],[10],[],[9,10],[8,9,10],[8,9,10],[],[10] ]),
 ];

s  := s-1; 
if s-3>Length(ffdata) then
Print("\n Sorry, currently only dimensions up to ",Length(ffdata)+3,".\n");
 return;
fi;

F:=Alff( ffdata[s-3] );                 # initialize the func.fld.
o:=AlffOrderMaxFinite(F);

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

                                           # Some abbreviations:
div := AlffDivisor;
lbas:= AlffDivisorLBasis;
ldim:= AlffDivisorLDim;
v   := AlffEltValuation;
fdeg:= AlffDeg(F);
eq  := AlffOrderEqFinite(F);
Y   := AlffElt(eq,Concatenation([0,1],0*[3..fdeg]));
one := AlffElt(o,Concatenation([1],0*[2..fdeg]));

p   := AlffPlacesDegOne(F); 
p0ind    := PositionProperty(p, pl->
             (pl in AlffPlaceSplit(F,T)) and (AlffPlaceRam(pl) = 1) );
p0       := p[p0ind];   # p0 : distinct nonramified place above T ( = P_\infty) 
p[p0ind] := p[1];
p     := p{[2..s+1]};
g   := AlffGenus(F);
D   := (g+g)*p[1];
c   := [];                   # (declaration of matrices)
cl  := []; 

# albas : compute an ascending lbas with respect to Place pl
#         and coefficients ranging from c0 to c1 

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 := g+g;
  else
   vi := 0;
  fi;
  nc  := [];
  bas := lbas(D+c0*pl);
  Sort(bas, function(a,b) return v(pl,one*a) > v(pl,one*b); end );
  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;

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

# xp :  List of coefficients of expansion of e at p0 up to T^n
#       starting with 1

xp  := function(e,n) 
        return Reversed( PolyToList( T^(n+1)+
          Sum(List(
           AlffPlaceLocals.EltExpansion(p0,T,AlffEltMove(e,o),n,false),
          i->T^i[2])))){[1..n]};
       end;

divpol:=function(pol1,pol2)
 local dg,pol,frac;
 dg   := Position(Reversed(PolyToList(pol2)),ffe[1]^0)-1;
 frac := pol2/T^dg;
 pol  := PolyQuotRem(pol1,frac);
 frac := pol[2]*frac^255;
 pol  := pol[1]+frac;
 return PolyQuotRem(pol,T^dg)[1];
end;

######## 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 : powers of Y, = basis of finite equation order

#######First try: (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^(m+m+g+g+2+2));   # Powers  

mx:=function(e,n) 
 local el;
 el:=AlffEltToList(AlffEltMove(one*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:=function(e,n) 
 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 := function(e,prec)
  local txp,te,i,xpp;
  xpp:=mx;                          # xp is slower, but safer
  txp := xpp(e,prec);
  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;


# pt : the generators of L(D+i p[j])

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 := List( [1..s], pl->Concatenation( 
  albas(p[pl],0,pn[pl])[1]{[g+2..g+pn[pl]+1]},
  [ First(lbas(pn[pl]*p[pl]),e->v(p[pl],one*e)=-pn[pl]) ]) ); 
PT := pt;
pt := List(pt,e->List(e,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;

Print(" All preliminaries done, now for the expansions...\n");
 PrintTo("matrices",m," ",s+1,"\n",num_c(m),"\n");
 # sorry for the obfuscation here .... :
 Exec("sed \"s/\[\\[,\]/ /g\" matrices > matr.new");
 Exec("sed \"s/\]/ /g\" matr.new > matrices");
 Exec("rm matr.new");
 Exec("mv matrices kashout");
Print(" Matrices constructed and written to 'kashout'\n");
#Print(" Now verifying the t-value (should be less/eq ",g,")...\n");
# if b=2 then
#  return Exec("/home/pirsic/csrc/tcalc/tcalcgf2gauss kashout 2"); 
# else
#  return Exec(Concatenation(
#    "/home/pirsic/csrc/tcalc/tcalcgauss    kashout ", SPrint(b))); 
# fi;
end;