# $Source: /u/maple/research/lib/mod/src/RCS/ReduceField,v $
# $Notify: mvanhoei@daisy.uwaterloo.ca
#
# Calling sequence:  ReduceField(A,vars) mod p
#
# Purpose: Reduction to a finite field.
#
#          If A is an element of the field K(vars) of rational functions in
#          vars with coefficients in K, where K is a field with characteristic
#          zero, then this procedure computes a modular image of A in F(vars)
#          where F is a finite field of characteristic p.
#
#          The input A may also be another expression over K, as long
#          as "A mod p" and indets(A,..) work. So lists and sets are OK, but
#          not tables or matrices.
#
#          The optional argument vars must be set of variables. If omitted
#          then vars={}. All indeterminates appearing in A, except the elements
#          of vars, will be replaced by random elements of Fp. Take p big in
#          order to avoid "division by zero" and such errors (20 digits should
#          be safe enough). If indets(A,'name') minus vars = {} then the only
#          randomness in this procedure is the choice of a factor in case a
#          RootOf factors in the reduced field. This will happen for example if
#          A has both RootOf(x^2-2) and RootOf(x^2-3), because the splitting
#          field of those two RootOf's over Fp has degree <= 2, so at least one
#          of them will factor over the other.
#
#          If there are several RootOf's, this procedure will write all of them
#          in terms of just 1 RootOf, so that the result is usable in other
#          Maple procedures (right now Maple only supports 1 RootOf in char p).
#          This can also be used to support multiple RootOf's mod p in Maple;
#          instead of "ERROR, not supported" if RootOf's [R1,R2,..] occur one
#          could also use ReduceField([R1,R2,...]) mod p instead of [R1,R2,..]
#          and continue the computation (although the reverse transformation is
#          not yet given in the output of this procedure).
#
# Author: Mark van Hoeij, May 1998.
#
#          Note: I wrote a procedure mod/Primfield that is the mod p
#          equivalent of evala/Primfield. To keep the syntax consistent with
#          evala, users should not use mod/ReduceField but use mod/Primfield
#          instead. Therefore only mod/Primfield will be documented by a
#          help page.


# Example application: compute the rank (not savelib'ed).
# Input: matrix. Nested RootOf's, radicals, and variables and
# such are allowed in the input, but not functions.
# Output: the rank
# Timings:
#   A := matrix(10,10,[seq(randpoly(x,degree=5),i=1..100)]):
#   time(); `linalg/rank_heu`(A); time(); linalg[rank](A); time();
`linalg/rank_heu`:=proc(A)
	local n,m,p,B,r,i,j;
	option `Copyright (c) 1998 Waterloo Maple Inc. All rights reserved.`;
	n:=linalg[rowdim](A);
	m:=linalg[coldim](A);
	p:=nextprime(rand(10^10..10^11)());
	# chance on wrong answer: approximately 1/p
	B:=[seq(seq(A[i,j],j=1..m),i=1..n)];
	Gausselim(matrix(n,m,ReduceField(B) mod p),'r') mod p;
	eval(r)
end:

# Reduce not necessarily irreducible RootOf's to one single
# irreducible RootOf by factorization.
#
# Remove parameters by random substitution in Fp. This is
# probabilistic. For this to work with high probability, the prime p
# must be big.
#
# vars=do not substitute values for these variables. They should
# not occur within the RootOfs.
#
`mod/ReduceField` := proc(A, vars, p)
local B, C, Bnew, sub, Br, F, f, x, i, roots, subC;
option `Copyright (c) 1998 Waterloo Maple Inc. All rights reserved.`;
    if nargs = 2 then return procname(A, {}, vars)
    elif not (type(p, posint) and vars=indets(vars,'name')) or
    has(indets(indets(A, 'RootOf'), 'name'), vars) then
        error "wrong input"
    end if;
    B := indets(A, 'name') minus vars;
    if B <> {} then
        return procname(subs(B[1] = rand(1 .. p)(), A), vars, p)
    end if;
    B := indets(A, {'RootOf', 'radical', 'nonreal'});
    if B = {} then A mod p
    elif member('radfielded', vars) or
    indets(A, {'radical', 'nonreal'}) = {} and
    mul(nops([op(i)]), i = B) = 1 then
		# B is set of RootOf's which are radfielded or which have
		# no indices.
        C := B minus `union`(seq(indets(op(1, i), 'RootOf'), i = B));
        C := C[1];
		# Find all RootOf's with the same minpoly (so the
		# only thing that differs is the index).
        C := [seq(`if`(op(1, i) = op(1, C), i, NULL), i = B)];
        B := `mod/ReduceField/sort_B`(B minus {op(C)});
        Bnew := procname(B, vars, p);
        sub := seq(B[i] = Bnew[i], i = 1 .. nops(B)), _Z = x;
        Br := op(indets(Bnew, 'RootOf')); # either NULL or just 1 RootOf.
        F := eval(subs(sub, op(1, C[1])));
        F := Factors(F, Br) mod p;
        subC := NULL;
        do
            F := {seq(i[1], i = F[2])};
		# pick factor of lowest degree.
            f := min(seq(degree(i, x), i = F));
            for i in F do if degree(i, x) = f then f := i; break end if
            end do;
            F := F minus {f};
            roots, sub :=
                `mod/ReduceField/roots_field`(f, {Br}, x, p, nops(C));
            Br := op(indets([roots, subs(sub, {Br})], 'RootOf'));
		# New field = Fp(Br)
            Bnew := subs(sub, Bnew); # This is B translated to the new field.
            f := min(nops(C), nops(roots));
            subC := seq(C[i] = roots[i], i = 1 .. f),
                seq(lhs(i) = subs(sub, rhs(i)), i = [subC]);
            C := C[f + 1 .. -1];
            if C = [] then return
                subs(subC, seq(B[i] = Bnew[i], i = 1 .. nops(B)), A) mod p
            elif F = {} then error "reduction failed"
		# Note that by trying a bit harder it would still be possible
		# to do the reduction.
            end if;
		# There are still some not-yet-assigned indexed RootOf's. Factor
		# the remaining factors of the polynomial in that RootOf. These
		# need not be irreducible anymore because the field Fp(Br) may
		# have changed.
            F := Factors(`*`(op(subs(sub, F))), Br) mod p
        end do
    else procname(subs(radfield([args])[1], A), vars union {'radfielded'},
        p)
    end if
end proc:


# Sort the RootOf's. The nested RootOf's come first,
# so that they will be the first in the substitution.
`mod/ReduceField/sort_B` := proc(B)
local C, i;
option `Copyright (c) 1998 Waterloo Maple Inc. All rights reserved.`;
    if B = {} then []
    else
        C := `union`(seq(indets(op(1, i), 'RootOf'), i = B));
        [op(B minus C), op(procname(C))]
    end if
end proc:


# Determine all roots of an irreducible polynomial f in the splitting
# field over the field given by F.  F={} or F={one RootOf}
`mod/ReduceField/roots_field` := proc(f, F, x, p, m)
local i, R, k, j;
option `Copyright (c) 1998 Waterloo Maple Inc. All rights reserved.`;
    if degree(f, x) = 1 or F = {} then
        R := RootOf(f, x) mod p;
        return `if`(degree(f, x) = 1 or m = 1, [R],
            [seq(i[1], i = Roots(f, R) mod p)]), 0 = 0
    elif m = 1 then
        i := `mod/ReduceField/generator`(f, x, F[1], p);
        if i <> FAIL then return i end if
    end if;
    do
        k := Randprime(degree(f, x)*degree(op(1, op(F)), _Z), x) mod p;
        R := RootOf(k, x) mod p;
        k := Resultant(subs(_Z = x + _Z, op(1, op(F))), op(R), _Z) mod p;
        k := Factors(k) mod p;
        if 1 < k[2][1][2] then next end if;
        k := Gcd(subs(x = x - R, k[2][1][1]), subs(_Z = x, op(1, op(F))))
            mod p;
        if degree(k, x) <> 1 then error  end if;
        k := op(F) = (-coeff(k, x, 0)/coeff(k, x, 1)) mod p;
        i := Expand(subs(x = x + R, subs(k, f))) mod p;
        i := Resultant(subs(R = _Z, i), op(R), _Z) mod p;
        i := Factors(i) mod p;
        if i[2][1][2] = 1 then return [seq(subs(`solve/linear`(
            {Gcd(subs(x = x - R, i[2][j][1]), subs(k, f)) mod p}, {x}), x)
            , j = 1 .. min(m, nops(i[2])))], k
        end if
    end do
end proc:


`mod/ReduceField/generator` := proc(f, x, R, p)
local n, r, g, y, A, q, AA, S, i;
    n := degree(f, x);
    if irem(n, p) = 0 then return FAIL end if;
    r := rand(1 .. p - 1)();
    g := subs(_Z = y, op(R));
    A := Resultant(subs(x = x - r*y, R = y, f), g, y) mod p;
    if 0 = Discrim(A, x) mod p then return FAIL end if;
    A := RootOf(A) mod p;
    q := p^degree(g, y);
    AA := A;
    S := A;
    for i to n - 1 do AA := Power(AA, q) mod p; S := S + AA end do;
    i := Expand((S + coeff(f, x, n - 1)/coeff(f, x, n))/(n*r)) mod p;
    [(A - r*i) mod p], R = i
end proc:

#savelib('`mod/ReduceField`',\
	 '`mod/ReduceField/sort_B`',\
	 '`mod/ReduceField/roots_field`',\
	 '`mod/ReduceField/generator`'):