macro(DUP=DenseUnivariatePolynomial);
SquareFreeNorm := proc(f,D,so)

     local Qxa, Qx, K, Q, Qa, alpha, g, m, r, rp, s, x, Qax, sub, X, Alpha;

# f is the polynomial to compute the square free norm of
# D is the domain (the algebraic extension)
#
# Returns a list of two things.  The first is an polynomial in the base
# domain (below the algebraic extension) which is the square free norm,
# the second is an integer which says what the transformation is that
# was used.
#
# D is univariate polynomials over the extension
# K is the algebraic extension
#
# x is the variable for polynomials in D
# m is the minimal polynomial
# alpha is the variable in Qa
#
# Q is the base field
# Qa is the polynomial domain for the minimal polynomial
# Qx is polynomials in x over the base field
# Qxa Q[x][alpha]
# Qax Q[alpha][x]

    K := D['CoefficientRing'];
    x := D['Variable'];
    m := K['Extension'];
    Qa := K['Representation'];
    alpha := Qa['Variable'];
    Q := K['BaseField'];
    Qx := DUP(Q,x);
    Qxa := DUP(Qx,alpha);
    Qax := DUP(Qa,x);
    X := Qax['Monomial']();
    Alpha := Qax['Constant'](Qa['Monomial']());

    m := Qax['Constant'](m);
    m := swapvar(m,Qax);
    s := 0;
    sub := Qax[`-`](X, Qax[`*`](s, Alpha));
    g := Qax['PolyEval'](f,sub);
    g := swapvar(g,Qax);
    r := Qxa['Resultant'](g,m);
    rp := Qx['Diff'](r);
    while not Qx['RelativelyPrime'](r,rp) do
        s := s + 1;
        sub := Qax[`-`](X, Qax[`*`](s, Alpha));
        g := Qax['PolyEval'](f,sub);
        g := swapvar(g,Qax);
        r := Qxa['Resultant'](g,m);
        rp := Qx['Diff'](r);
    od;
    so := s;
    r;
end:

swapvar := proc(p,Dxy)

    local Dx, Dy, Dyx, c, dy, i, j, p_out, yi_coeff, D;

    Dx := Dxy['CoefficientRing'];
    D := Dx['CoefficientRing'];
    Dy := DenseUnivariatePolynomial(D,Dxy['Variable']);
    Dyx := DenseUnivariatePolynomial(Dy,Dx['Variable']);
    
    dy := Dxy['Degree'](p);

    p_out := Dyx[0];
    for i from 0 to dy do
	yi_coeff := Dxy['Coeff'](p,i);
	for j from 0 to Dx['Degree'](yi_coeff) do
	    c := Dx['Coeff'](yi_coeff,j);
	    c := Dy['Monomial'](c,i);
	    c := Dyx['Monomial'](c,j);
	    p_out := Dyx[`+`](p_out,c)
	od
    od;
    p_out;
end:

save `SFN.m`;
quit
