#
# Computation of the ridge and the directrix of a multihomogeneous ideal for
# Maple 16
#
# File: multihomogeneous_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:
#   mh_ridge ([(a+b)^3*x^2+c^3*(y+z)^2], 0, tdeg, [[a,b,c],[x,y,z]]);
#   mh_ridge ([(a+b)^3*x^2+c^3*(y+z)^2], 2, grlex, [[a,b,c],[x,y,z]]);
#   mh_ridge ([(a+b)^3*x^2+c^3*(y+z)^2], 3, "FGb", [[a,b,c],[x,y,z]]);
#   mh_directrix ([a^4*(x+y+z),(b+c)^2*z^3], 0, tdeg, [[a,b,c],[x,y,z]]);
#   mh_directrix ([a^4*(x+y+z),(b+c)^2*z^3], 2, grlex, [[a,b,c],[x,y,z]]);
#   mh_directrix ([a^4*(x+y+z),(b+c)^2*z^3], 3, "FGb", [[a,b,c],[x,y,z]]);
#   mh_ideal_retrieving ([a^4*(x+y+z),(b+c)^2*z^3], 0, tdeg,
#                        [[a,b,c],[x,y,z]]);
#   mh_ideal_retrieving ([a^4*(x+y+z),(b+c)^2*z^3], 2, grlex,
#                        [[a,b,c],[x,y,z]]);
#   mh_ideal_retrieving ([a^4*(x+y+z),(b+c)^2*z^3], 3, "FGb",
#                        [[a,b,c],[x,y,z]]);

read "homogeneous_ideal_ridge.mpl";

to_one_list := proc (J::list, m)
    return [seq (J[i][], i=1..m)];
end proc:

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

mh_Hasse_Schmidt := proc (J::list, vars::list, add_vars::list, m, n, p, ord)
local L, all_vars, i, repl, f, coeffs_f, g, d;
    L := [seq([], i=1..m)];
    all_vars := to_one_list (vars, m);
    for i from 1 to m do
        repl := {seq (vars[i][j] = vars[i][j] + add_vars[i][j], j=1..n[i])};
        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`), all_vars)];
            for g in coeffs_f do
                d := degree (g, add_vars[i]);
                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[i])))
                then L[i] := [L[i][], g];
                end if;
            end do;
        end do;
    end do;
    if (ord = "FGb")
    then return [seq (my_fgb_gbasis (L[i], p, add_vars[i],
                               max_degree (L[i], add_vars[i])[1]), i=1..m)];
    else return [seq (InterReduce (L[i], ord (add_vars[i][]),
                                   characteristic = p), i=1..m)];
    end if;
end proc:

# computation of the ridge of a multihomogeneous 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 list of variables
# output:
# list of list of additive polynomials
mh_ridge := proc (J::list, p, ord, vars)::list;
local J_x, m, n, v_x, v_y, G;
    J_x := expand (J);
    if (p > 0) then J_x := J_x mod p; end if;
    m    := nops (vars);
    n    := [seq (nops (vars [i]), i=1..m)];
    v_x  := [seq ([seq (x[i][j], j=1..n[i])], i=1..m)];
    v_y  := [seq ([seq (y[i][j], j=1..n[i])], i=1..m)];
    J_x  := subs ({seq (seq (vars[i][j] = x[i][j], j=1..n[i]), i=1..m)}, J_x);
    G    := mh_Hasse_Schmidt (mh_Giraud (J_x, v_x, m, p, ord),
                              v_x, v_y, m, n, p, ord);
    return subs ({seq (seq (y[i][j] = vars[i][j], j=1..n[i]), i=1..m)}, G);
end proc:


# computation of the directrix of a multihomogeneous 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 list of variables
# mh_rid the ridge of J
# output :
# list of list of linear forms
mh_directrix := proc (J::list, p, ord, vars,
                      mh_rid:=mh_ridge (J, p, ord, vars))::list;
local R, m, D;
    R := mh_rid;
    if (p = 0) then return R; end if;
    m := nops (vars);
    D := [seq ([seq (pol_to_linear (R[i][j], vars[i], nops (vars[i])),
                     j=1..nops(R[i]))], i=1..m)];
    if (ord = "FGb")
    then return [seq (my_fgb_gbasis (D[i], p, vars[i], 1), i=1..m)];
    else return [seq (InterReduce (D[i], ord(vars[i][]), characteristic = p),
                      i=1..m)];
    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 list of variables
# mh_dir the directrix of J
# output :
# list of polynomials in the linear forms of the directrix
mh_ideal_retrieving := proc (J::list, p, ord, vars,
                             mh_dir:=mh_directrix (J, p, ord, vars))::list;
local m, D, D_lt, c, i;
    m    := nops (vars);
    D    := to_one_list (mh_dir, m);
    c    := nops (D);
    D_lt := [seq (LeadingMonomial (D[i], grlex (to_one_list (vars, m)[])),
                  i=1..c)];
    return [seq (pol_retrieving (f, p, D, c, D_lt), f in J)];
end proc:


# f := (a+b+c+d)*(i+j+k)*(x+y):
# g := (a+b)^2*(i+j+k)^3*t^2*z + c^2*(i+j)^3*(x+y)^3:

# for p in [0, 2, 3] do
#     print ("hom ridge", ridge ([f], p, tdeg, [a,b,c,d,i,j,k,t,x,y,z]));
#     print ("m-h ridge", mh_ridge ([f], p, tdeg,
#                                   [[a,b,c,d],[i,j,k],[t,x,y,z]]));
#     print ("hom ridge", ridge ([g], p, grlex, [a,b,c,d,i,j,k,t,x,y,z]));
#     print ("m-h ridge", mh_ridge ([g], p, grlex,
#                                   [[a,b,c,d],[i,j,k],[t,x,y,z]]));
#     print ("hom ridge", ridge ([f,g], p, "FGb", [a,b,c,d,i,j,k,t, x,y,z]));
#     print ("m-h ridge", mh_ridge ([f,g], p, "FGb",
#                                   [[a,b,c,d],[i,j,k],[t,x,y,z]]));
# end do: