/*This makes the Fleming designs from a plane, and constructs
the p-ary code and the automorphism group in each case */

//K.Mackenzie-Fleming "A recursive construction for 2-designs",
//Designs, Codes and Cryptography, 13 (1998), 159--164

/* Input p and s to get q=p^s and a symmetric design on
q^3+q^2+q+1 points and with block size q^2+q+1,
and then a symmetric design on q^4+q^3+q^2+q+1 points and 
with block size q^3+q^2+q+1,*/

  d:=2;t:=1;

//Some functions
/*to construct a t-dimensional space x_1=x_2=...x_d-t=0
in AG(d,q)*/
  afspace:=function(pts,vv,d,t);
  ttsp:=sub<vv| {Basis(vv)[i]: i in [(d-t+1)..d]}>;
  asp:={ Position(pts,x): x in ttsp };
  return asp;
  end function;

/*to construct a projective t-dimensional space 
x_1=x_2=...x_d-t=0 in PG(d,q) */
  pgspace:=func<pts,vv,d,t|
  { Position(pts,x): x in sub<vv|
  {Basis(vv)[i]: i in [(d-t+1)..d+1]}> } diff { 0 }>;

  SQ:=SetToSequence;
  QS:=SequenceToSet;
  blocksp:=func<blox,pt|[blox[i]:i in [1..#blox]| pt in blox[i]]>;
  qwt:=func<u,q|  &+[Intseq(u,q)[j]: j in [1..#Intseq(u,q)]]>;
  grm:=func<d,t,p|
  &+[&+[(-1)^k *Binomial(d,k)*Binomial(i-k*p+d-1,i-k*p):
  k in [0..(i div p)]]:i in [0..(d-t)*(p-1)]]>;

  pgrm:=func<d,t,p|&+[ grm(j,t,p) : j in [t..d]]>;

/*Want Rank_p(PG_{m,r}(F_q)) where q=p^t*/
  PGrank:=function(m,r,p,t);
  if t eq 1 then
   return pgrm(m,r,p);
    else
     m:=m+1;rr:=r+1;
      q := p^t;
       n := (q^m)-1;
        np :=  n div (q-1);
         wts := [qwt(u,q) : u in [1..n]];
          k := ( m-rr)*(q-1);
           we:=&meet[ {i : i in [1..q^m -2]
           |wts[(p^j)*i mod  n] le k }: j in [0..t-1]];
           wee:={ i: i in we | (i mod (q-1)) eq 0};
   end if;
  return #wee+1;
  end function;
//now the real thing

/* feed in s,d,t and p where q=p^s for AG_{d,t}(F_q)*/
q:=p^s;
"geometry",p,"-rank=",PGrank(3,2,p,s);
  "q=",q,"p=",p;
  f := GaloisField(q);
  w:=PrimitiveElement(f);
  vv := VectorSpace(f, d);
  gg, pts := AffineGeneralLinearGroup(vv);
  tsp:=afspace(pts,vv,d,t);
  blox:=SetToSequence(tsp^gg);

  afdes:=AffinePlane<q^2|blox>;
  afblox:=LineSet(afdes);
  afb:=[Support(afblox.i):i in [1..q^2+q]];
  n:=Position(afb,{1..q});
  pc:=ParallelClass(afblox.n);
  spc:={Support(x):x in pc};

  pa:=ParallelClasses(afdes);
  spa:=SQ(pa);
  sspa:=[SQ(spa[i]):i in [1..#spa]];
  mspa:=[[[{x+(q^2)*k:x in Support(sspa[i][j])}:j in [1..#sspa[i]]]:
       i in [1..#sspa]]:k in [0..q]];
  ablox:=blox;

/*making the design of points and t-spaces in PG(d,q)..
pts and blox, as usual*/
  vv := VectorSpace(f, d+1);
  gg, pts := ProjectiveGeneralLinearGroup(vv);

  tsp:=pgspace(pts,vv,d,t);
  gblox:=SetToSequence(tsp^gg);

  pblox:=[ {x+q^3:x in gblox[i]}:i in [1..#gblox]];
  /*pencil ordering from line at infinity*/
  npblox:=&cat[Remove(blocksp(pblox,x),q+1):x in pblox[q^2+q+1]];
  npblox:=Append(npblox,pblox[q^2+q+1]);
  pblox:=npblox;
  for l:=1 to q^2+q do
  pblox:=Rotate(Remove(pblox,q^2+q+1),l) cat [npblox[q^2+q+1]];
  "Design No.",l;
  nblox:= [{i*q^2+1..i*q^2+q^2} join pblox[q^2+q+1]:i in [0..q-1]]
  cat [{q^3+1..q^3+q^2+q+1}] ;
  qm:=[QS(mspa[1][i]):i in [1..#mspa[1]]];
  cspa:=Remove(mspa[1],Position(qm,spc));

  zs:=&cat[[<i,j>:j in [1..q]]:i in [1..q]];
  csp:=cspa[1];
  ccsp:=[Sort(SQ(csp[i])):i in [1..q]];
  nseq:=&cat[[ccsp[i][j]:i in [1..q]]:j in [1..q]];

  tblox:=[];
  for m:=1 to q+1 do
    for j:=1 to q do
      for k:=1 to q do
       cl:=&join[mspa[zs[Position(nseq,x)][1]][m][zs[Position(nseq,x)][2]]:
       x in cspa[j][k]] join pblox[(m-1)*q+j];
       tblox:=Append(tblox,cl);
      end for;
    end for;
  end for;

  fblox:=tblox cat nblox;
  #fblox,"= no. of blocks constructed";
  des:=Design<2,q^3+q^2+q+1|fblox>;
  "Constructed a symmetric",des;
  dc:=LinearCode(des,GF(p));
  "    ",Dimension(dc),"=",p,"-rank";
  au:=PointGroup(des);
  "       automorphism group size=",Order(au);
"Now construct the bigger symmetric design...";
"Finite geometry design",p,"-rank=",PGrank(4,3,p,s);
B:=BlockSet(des);

rdes:=Residual(des,B.(q^3+q^2+q+1));
t,resol:=IsResolvable(rdes);
pres:=SQ(resol);
pres:=[SQ(pres[i]): i in [1..#pres]];

sspa:=pres;
  nmspa:=[[[{x+(q^3)*k:x in Support(sspa[i][j])}:j in [1..#sspa[i]]]:
       i in [1..#sspa]]:k in [0..q]];
Blox:=[Support(B.i):i in [1..#B]];
np:=Position(Blox,{q^3+1..q^3+q^2+q+1});

  Pblox:=[ {x+q^4:x in Blox[i]}:i in [1..#Blox]];
  Pblox:=Append(Remove(Pblox,np),{q^4+q^3+1..q^4+q^3+q^2+q+1});

  nblox:= [{i*q^3+1..i*q^3+q^3} join Pblox[q^3+q^2+q+1]:i in [0..q-1]]
  cat [{q^4+1..q^4+q^3+q^2+q+1}] ;

  tblox:=[];
  for m:=1 to q^2+q+1 do
    for j:=1 to q do
      for k:=1 to q do
       cl:=&join[nmspa[zs[Position(nseq,x)][1]][m][zs[Position(nseq,x)][2]]:
       x in cspa[j][k]] join Pblox[(m-1)*q+j];
       tblox:=Append(tblox,cl);
      end for;
    end for;
  end for;

  fblox:=tblox cat nblox;
  #fblox,"= no. of blocks constructed";
  ndes:=Design<2,q^4+q^3+q^2+q+1|fblox>;
  "Constructed a symmetric",ndes;
  dc:=LinearCode(ndes,GF(p));
  "    ",Dimension(dc),"=",p,"-rank";
  au:=PointGroup(ndes);
  "       automorphism group size=",Order(au);
   "-o-o-o-o-o-o-o-o-o-o-";
end for;