# 8/3/18  modified to divide out A by X[i]-Y[i] before computing det(A) !!
# 17/6/19  modified to build matrix by saving coefficients in X

read detGJcount;

dixonmatrix := proc(F::list(polynom),X::list(name),Y::list(name))
local n,A,t,i,j,k,det,delta,multicoeff,monsX,mon,monsY,rows,cols,M,cof,B,cofsX,cofsY,numterms;
global Delta;

    n := nops(X);
    if nops({op(X)}) < n then error "X has duplicates" fi;

    if nargs=2 then return dixonmatrix(F,X,[seq(t[i],i=1..n)]) fi;

printf("Dixon: n=%d\n",n);

    if n <> nops(F)-1 then error "not implemented" fi;
    A := Matrix(n+1,n+1);
    for j to n+1 do A[1,j] := F[j] od;
    for i to n do
        for j to n+1 do
            A[i+1,j] := subs(X[i]=Y[i],A[i,j]);
        od;
    od;


    # Don't compute det(A) then divide by product( X[i]-Y[i], i=1..n )
    # Robert's trick : Row(A,i+1)-Row(A,i) is divisible by X[i]-Y[i]

    for i from n by -1 to 1 do
        for j from 1 to n+1 do
            if not divide( A[i+1,j]-A[i,j], X[i]-Y[i], evaln(A[i+1,j]) ) then print(i,j,bad) fi;
        od;
    od;

printf("Cancellation matrix done\n");

for j to n+1 do t := A[1,j]; A[1,j] := A[n+1,j]; A[n+1,j] := t; od;
#printf("Dixon: #minors=%a\n",determinant(A));

    numterms := proc(f) if f=0 then 0 elif type(f,`+`) then nops(f) else 1 fi; end;

    det := LinearAlgebra:-Determinant;
    delta := det(A,method=minor);
    delta := expand( delta ); # expand  because of coeffs

#if display then 
#   Delta := delta; print('Delta'=collect(delta,Y,distributed)); 
#fi;


printf("Dixon: #delta=%d\n",numterms(delta));

    cofsX := coeffs(delta,X,'monsX');
    cofsX := [cofsX];
    monsX := [monsX];
    rows := nops(monsX); 
    coeffs(delta,Y,'monsY');
    monsY := [monsY];
    cols := nops(monsY);
    for j to cols do M[monsY[j]] := j; od;

printf("Dixon: %d x %d\n",rows,cols);

    A := Matrix(rows,cols,storage=sparse); # free old A
    for i to rows do
        cof := cofsX[i];
        if cof=0 then next fi;
        cofsY := coeffs(cof,Y,'monsY');
        cofsY := [cofsY];
        monsY := [monsY];
        for k to nops(monsY) do j := M[monsY[k]]; A[i,j] := cofsY[k]; od;
    od;

    A; # the dixon matrix

end: