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))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 rnkFAIL 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: