#
# Computation of the ridge and the directrix of a homogeneous ideal for
# Maple 16
#
# File: homogeneous_ideal_ridge.mpl
# Author: Jérémy Berthomieu
# Address: LIP6, Université Pierre-et-Marie-Curie, Boîte Courrier 169,
# 	   4 place Jussieu, 75252 Paris Cedex 05, France
# Date: January 2013
# Licence: GPL v2
#
# Example:
#   ridge ([(a+b)^3*x^2+c^3*(y+z)^2], 0, tdeg);
#   ridge ([(a+b)^3*x^2+c^3*(y+z)^2], 2, grlex, [a,b,c,x,y,z]);
#   ridge ([(a+b)^3*x^2+c^3*(y+z)^2], 3, "FGb", [a,b,c,x,y,z]);
#   directrix ([a^4*(x+y+z),(b+c)^2*z^3], 0, tdeg);
#   directrix ([a^4*(x+y+z),(b+c)^2*z^3], 2, grlex, [a,b,c,x,y,z]);
#   directrix ([a^4*(x+y+z),(b+c)^2*z^3], 3, "FGb", [a,b,c,x,y,z]);
#   ideal_retrieving ([a^4*(x+y+z),(b+c)^2*z^3], 0, tdeg);
#   ideal_retrieving ([a^4*(x+y+z),(b+c)^2*z^3], 2, grlex, [a,b,c,x,y,z]);
#   ideal_retrieving ([a^4*(x+y+z),(b+c)^2*z^3], 3, "FGb", [a,b,c,x,y,z]);

with (Groebner):
with (FGb):

my_fgb_gbasis := proc (J::list, p, vars::list, d_max)
    return fgb_gbasis (J, p, vars, [], {"dlim"=d_max, "index"=2000000});
end proc:

max_degree := proc (J::list, vars::list)
local r, d_max, b, i, d;
    r     := nops (J);
    if (r = 0) then return [-1, 1]; end if;
    d_max := degree (J[1], vars);
    b     := 1;
    for i from 2 to r do
        d := degree (J[i], vars);
        if (d > d_max)
        then
            d_max := d;
            b     := 0;
        end if;
    end do;
    return [d_max, b];
end proc:

Giraud := proc (J::list, vars::list, p, ord)
local max_deg, Groeb, G, f;
    max_deg := max_degree (J, vars);
    if (max_deg[2] = 1)
    then return J;
    else
        if (ord = "FGb")
        then return my_fgb_gbasis (J, p, vars, max_deg[1]);
        else
            Groeb := Basis (J, ord (vars[]), characteristic = p);
            G     := [];
            for f in Groeb do;
                if (degree (f, vars) <= max_deg [1])
                then G := [G[], f];
                end if;
            end do;
            return G;
        end if;
    end if;
end proc:

Hasse_Schmidt := proc (J::list, vars::list, add_vars::list, n, p, ord)
local L, repl, f, coeffs_f, g, d;
    L    := [];
    repl := {seq (vars[i] = vars[i] + add_vars[i], i=1..n)};
    for f in J do
    	f := expand (subs (repl, f));
        if (p > 0) then f := f mod p; end if;
	coeffs_f := [coeffs (collect (f, `distributed`), vars)];
        for g in coeffs_f do
            d := degree (g, add_vars);
            if ((p = 0 and d = 1 and not (g in L)) or  # char 0
                (p > 0 and d <> 0 and                  # char p
                 is(log[p](d)::AndProp (integer, RealRange(0, infinity)))
                 and not (g in L)))
            then L := [L[], g];
            end if;
        end do;
    end do;
    if (ord = "FGb")
    then return my_fgb_gbasis (L, p, add_vars, max_degree (L, add_vars)[1]);
    else return InterReduce   (L, ord (add_vars[]), characteristic = p);
    end if;
end proc:

pol_to_linear := proc (F, vars::list, n)
local d, L, i;
    d := degree (F, vars);
    if (d = 1) then return F; end if;
    L := 0;
    for i from 1 to n do
        L := L + coeff (F, vars[i], d)*vars[i];
    end do;
    return L;
end proc:

pol_retrieving := proc (F, p, dir::list, c, dir_lt)
local f, i;
    f := collect (expand (F), `distributed`);
    for i from 1 to c do
        f := collect (expand (algsubs (dir[i]=dir_lt[i], f)), `distributed`);
    end do;
    if (p > 0) then f := f mod p; end if;
    return f;
end proc:

# computation of the ridge of a homogeneous ideal
# input:
# J     list of polynomials
# p     characteristic of the field of coefficients
# ord   monomial order or "FGb" to use FGb library
# vars  list of variables
# output:
# list of additive polynomials
ridge := proc (J::list, p, ord, vars := [op(indets (J))])::list;
local J_x, n, v_x, v_y, G;
    if (J = []) then return J; end if;
    J_x := expand (J);
    if (p > 0) then J_x := J_x mod p; end if;
    n    := nops (vars);
    v_x  := [seq (x[i], i=1..n)];
    v_y  := [seq (y[i], i=1..n)];
    J_x  := subs ({seq (vars[i] = x[i],  i=1..n)}, J_x);
    G    := Hasse_Schmidt (Giraud (J_x, v_x, p, ord), v_x, v_y, n, p, ord);
    return subs ({seq (y[i] = vars[i], i=1..n)}, G);
end proc:

# computation of the directrix of a homogeneous ideal
# input:
# J     list of polynomials
# p     characteristic of the field of coefficients
# ord   monomial order or "FGb" to use FGb library
# vars  list of variables
# rid   the ridge of J
# output :
# list of linear forms
directrix := proc (J::list, p, ord, vars := [op(indets (J))],
	           rid := ridge (J, p, ord, vars))::list;
local R, D;
    if (J = []) then return J; end if;
    R := rid;
    if (p = 0 or R = []) then return R; end if;
    D := [seq (pol_to_linear (R[i], vars, nops (vars)), i=1..nops(R))];
    if (ord = "FGb")
    then return my_fgb_gbasis (D, p, vars, 1);
    else return InterReduce (D, ord(vars[]), characteristic = p);
    end if;
end proc:

# computation of the original ideal
# input:
# J     list of polynomials
# p     characteristic of the field of coefficients
# ord   monomial order or "FGb" to use FGb library
# vars  list of variables
# dir   the directrix of J
# output :
# list of polynomials in the linear forms of the directrix
ideal_retrieving := proc (J::list, p, ord, vars := [op(indets (J))],
                          dir := directrix (J, p, ord, vars))::list;
local D, D_lt, c, i;
    if (J = []) then return []; end if;
    D    := dir;
    c    := nops (D);
    D_lt := [seq (LeadingMonomial (D[i], grlex (vars[])), i=1..c)];
    return [seq (pol_retrieving (f, p, D, c, D_lt), f in J)];
end proc:


# ridge ([(x+y)*z, (x+y)^2], 0, tdeg, [x,y,z]);
# ridge ([(x+y)^4+z^4, (x+y)^4*z^4], 2, grlex);
# ridge ([x^3*y^2*z, (x+y)^6], 0, tdeg);
# ridge ([x^3*y^2*z, (x+y)^6], 2, grlex, [x,y,z]);
# ridge ([x^3*y^2*z, (x+y)^6], 3, "FGb", [x,y,z]);