CycloSolve := proc(A::Matrix,b::Vector,M::polynom(rational),z::name,N::integer)
#
# Solve Ax=b mod M(z) where M(z) the cyclotomic polynomial of order N.
#
# From Cramer's rule we have  x[i] = det(A(i)) / det(A) mod M(z)
# where A(i) is A with the i'th column replaced by b.
# This code computes the solution vector x satisfying x[i] in Q[z] using rational
# number reconstruction, and, simulataneously, the vector y satisfying
# y[i] = det(A(i)) mod M(z) and determinant D = det(A) mod M(z).
# Notice that because M(z) is monic over Z, y[i] and D are over Z and hence
# can be recovered using Chinese remaindering only.
# So the code simultaneously reconstructs x[i], y[i] and D from images modulo primes.
# It stops when the first method succeeds.
# The advantage of the y[i], D representation is that D is not inverted mod M(z)
# which produces a blowup of a factor of d=deg(M(z))=phi(N), in general.
# That is the size of the integers in the vector y and D are d times smaller
# than the integers in x in general, thus requiring a factor of 2d fewer primes.
# However, it is possible, and it does occur in practice that the size of the
# integers in x are MUCH smaller than those in y and D, therefore, the desire
# to stop when the first succeeds.
#
# Author: Michael Monagan, June 2008.
#
local m,n,p,q,den,Ab,Ap,Aa,K,r,R,alpha,d,i,j,x,X,P,v,Y,inc,rnk,y,C,Z,upto,
      noMod,HA,Hm,Hb,Hy,T,rnkBnd,D,det,DET,a,adj,ADJ,deltaX,deltaD,HD;

  # On a 32 bit machine, 24 bit floating point primes are better than 15.5
  # bit integer primes both in size (we don't want to run out of primes)
  # and they are usually faster.  On a 64 bit machine 30 bit integer primes
  # may or may not be faster than floats.  They are on an AMD Opteron.
  if modp1(Prime(1)) < 2^16 then T := float[8]; else T := integer[8]; fi;

  m,n := op(1,A);
  if m<n then error "not enough equations" fi;
  if m<>op(1,b) then error "vector dimension is not compabible with matrix" fi;
  if not type(M,polynom(rational,z)) then 
     error sprintf("minimal polynomial must be in Q[%a]",z)
  fi;
  d := degree(M,z);
  userinfo(3,solve,sprintf("Solving A x = b mod M where A is %d by %d and M = %a",m,n,M));
  if T=float[8] then
     p := 2^24; # the largest floating point prime must be less than 2^25
                # but 2^24 gives a faster result by avoiding divisions mod p
     userinfo(3,solve,"modulo 24 bit primes using double precision floating point.");
  else # T=integer[8]
     p := 2^30; # modp1(Prime(1)) which is 31.5 bits on a 64 bit machine
                # 2^30 is faster as it allows 16 additions before division
                # mod p is necessary (assuming unsigned 64 bit arithmetic)
     userinfo(3,solve,"modulo 30 bit primes using 64 bit integer arithmetic.");
  fi;

  # Bounds.
  HA := max( seq( seq( maxnorm(A[i,j]), j=1..n ), i=1..m ) );
  Hb := max( seq( maxnorm(b[i]), i=1..m ) );
  Hm := maxnorm(M);

  # Theorem: M splits iff p == 1 mod N
  q := iquo(p,N); p := N*q+1+N;

  Ab := Matrix(1..m,1..n+1,A); for i to m do Ab[i,n+1] := b[i] od; # Ab = <A|b>
  Ap := Matrix(m,n+1,'datatype'=T,'order'='C_order');
  y := Array(1..d);
  C := 0; # count degeneracies
  P := 1; # product of moduli
  rnkBnd := 0; # lower bound on rank
  # because map(modp,Ab) is not very efficient I'm avoiding it
  noMod := true;
  for i to m while noMod do
      Z := indets({seq(Ab[i,j],j=1..n+1)},rational); 
      noMod := type(Z,set(integer)) and min(op(Z))>-p and max(op(Z))<p;
  od;
  upto := 1;
  use LinearAlgebra:-Modular in
  K := 0;
  do
     p := p-N; while irem(p,2)=0 or not isprime(p) do p := p-N od;
     R := Roots(M) mod p;
     #if nops(R)<d then error "this should not happen" fi;
     if nops(R)<d then
		# this should not happen
		error "unable to solve";
     fi;

     alpha := [seq(r[1], r=R)];
     userinfo(4,solve,sprintf("prime[%d]=%d",K,p));
     if noMod then Aa := Ab; else 
        if T=float[8] then Aa := map(modp,Ab,p);
        else #Aa := rtable(1..m,1..n+1,
             #     [seq([seq(`if`(Ab[i,j]=0,0,modp1(ConvertIn(Ab[i,j],z),p)), j=1..n+1)],i=1..m)],
             #     subtype=Matrix );
             Aa := Matrix(m,n+1);
             for i to m do for j to n+1 do 
                 if Ab[i,j]<>0 then Aa[i,j] := modp1(ConvertIn(Ab[i,j],z),p); fi;
             od od;
        fi;
     fi;
     for i to d do
         Mod(p,Aa,z=alpha[i],Ap);
         RowReduce(p,Ap,m,n+1,n,`if`(m=n,'det',0),0,'rnk',0,'inc',true);
         if inc>0 then return NULL; fi; # system is proven inconsistent
         if rnk<n then break; fi; # this could be just unlucky
         rnkBnd := n;
         y[i] := [seq(trunc(Ap[j,n+1]),j=1..n)];
         if m=n then a[i] := modp(det*y[i],p); D[i] := det; fi;
     od;
     if rnk<n then
        if rnkBnd=n then next; fi; # prime or alpha[i] is unlucky
        C := C+1;
        if C=2 then
           userinfo(3,solve,"system is probably singular");
           return FAIL;
        fi;
        next;
     fi;
     # Try to recover the solution vector using rational reconstruction.
     x := [seq( [seq( y[i][j], i=1..d )], j=1..n )]; # = transpose(y)
     x := [seq( Interp(alpha,x[j],z) mod p, j=1..n )];
     if assigned(X) # then Chinese remaindering
     then i := 1/P mod p; v := i*(x-X) mod p; X,P := X+P*v, p*P;
     else X,P := x,p
     fi;
     K := K+1;
     if K = upto then
        userinfo(3,solve,sprintf("attempting rational reconstruction after %d primes",K));
        Y := iratrecon(X,P,'maxquo','scaled');
        if Y<>FAIL then
           userinfo(3,solve,"rational reconstruction succeeded: checking the solution");
           Y := Vector(Y);
           den := ilcm( seq(denom(Y[i]),i=1..n) );
           Z := den*Y;
           Hy := max( seq( maxnorm(Z[i]), i=1..n ) );
           #Z := A.(den*Y)-den*b; # clear denominators from Y first
           #for i to n do if not divide(Z[i],M) then break fi; od;
           if P > 2*( n*d*HA*Hy*(1+Hm)^(d-1) + den*Hb ) then return Y; fi;
        else
           # Attempt rational reconstruction for K in {1,2,3,5,8,12,18,27,...}
           upto := iquo(3*upto+1,2);
        fi;
     fi;
     if m=n then # Try to recover y[i] and D using Chinese remaindering
        P := P/p;
        x := [seq( [seq( a[i][j], i=1..d )], j=1..n )]; # = transpose(a)
        x := [seq( modp(Interp(alpha,x[j],z),p), j=1..n )];
        det := modp(Interp(alpha,[seq(D[i],i=1..d)],z),p);
        adj := mods(x,p);
        det := mods(det,p);
        if assigned(ADJ) # then Chinese remaindering
        then i := modp(1/P,p);
             deltaX := mods(i*(x-mods(ADJ,p)),p); ADJ := ADJ+P*deltaX;
             deltaD := mods(i*(det-mods(DET,p)),p); DET := DET+P*deltaD;
             P := p*P;
        else DET,ADJ,P := det,adj,p
        fi;
        if convert(deltaX,set) = {0} and deltaD = 0 then
           userinfo(3,solve,"Chinese remaindering succeeded: checking the solution");
           Hy := max( seq(maxnorm(ADJ[i]),i=1..m) );
           HD := maxnorm(DET);
           if P > 2*(d*n*HA*Hy*(1+Hm)^(d-1) + d*HD*Hb*(1+Hm)^(d-1)) then
              Y := Vector([seq(normal(ADJ[i]/DET),i=1..n)]); # form ratios
              return Y;
           fi;
        fi;
     fi;
  od;
  end use;
end: