# $Source: /u/maple/research/lib/DEtools/diffop/src/RCS/compute_modp,v $
# $Notify: hoeij@sci.kun.nl $

# Procedures for users
#    none
# Procedures for use only in the rest of the diffop package:
#    compute_modm, reconstruct, modulus2  (only used in try_factorization in DFactor)
# Procedures for use only in this file:
#    compute_modp, modp_table


# Purpose of this file: switch to modular arithmetic, and when done
# then repair the damage that was done to the tables that contain the
# values of the laurent series.


# macro's for g_ext:
macro(
	modulus=`DEtools/diffop/modulus`,
	x=`DEtools/diffop/x`,
	DF=`DEtools/diffop/DF`,
	xDF=`DEtools/diffop/xDF`
):

macro(
	g_conversion1=`algcurves/g_conversion1`,
	g_conversion2=`algcurves/g_conversion2`,
	rootof=`algcurves/rootof`,
	degree_ext=`algcurves/degree_ext`,
	g_evala=`algcurves/g_evala`,
	g_evala_rem=`algcurves/g_evala_rem`,
	g_zero_of=`algcurves/g_zero_of`,
	v_ext_m=`algcurves/v_ext_m`,
	ext_to_coeffs=`algcurves/e_to_coeff`,
	g_ext_r=`algcurves/g_ext_r`,
	g_gcdex=`algcurves/gcdex`,
	truncate=`algcurves/truncate`,
	g_factors=`algcurves/g_factors`
);
macro(
	new_laurent2=`DEtools/diffop/new_laurent2`,
	differentiate=`DEtools/diffop/differentiate`,
	g_ext=`DEtools/diffop/g_ext`,
	g_factors=`DEtools/diffop/g_factors`,
	g_normal=`DEtools/diffop/g_normal`,
	set_laurents=`DEtools/diffop/set_laurents`,
	description_laurent=`DEtools/diffop/description_laurent`,
	modm=`DEtools/diffop/modm`,
	g_solve=`DEtools/diffop/g_solve`,
	g_quotient=`DEtools/diffop/g_quotient`,
	g_rem=`DEtools/diffop/g_rem`,
	g_quo=`DEtools/diffop/g_quo`,
	truncate_x=`DEtools/diffop/truncate_x`,
	g_expand=`DEtools/diffop/g_expand`
):

# macro's for eval_laurent:
macro(	accuracy_laurent=`DEtools/diffop/accuracy_laurent`,
	value_laurent=`DEtools/diffop/value_laurent`,
	LL=_LL
):
macro(
	lowerbound_val=`DEtools/diffop/lowerbound_val`,
	max0_val=`DEtools/diffop/max0_val`,
	new_laurent=`DEtools/diffop/new_laurent`,
	nmterms_laurent=`DEtools/diffop/nmterms_laurent`,
	nterms_expression=`DEtools/diffop/nterms_expression`,
	nthterm_laurent=`DEtools/diffop/nthterm_laurent`,
	ramification_laur=`DEtools/diffop/ramification_laur`,
	lift_ramification_laur=`DEtools/diffop/lift_ramification_laur`,
	eval_laurent=`DEtools/diffop/eval_laurent`,
	subsvalueslaurents=`DEtools/diffop/subsvalueslaurents`,
	series_val=`DEtools/diffop/pseries_val`,
	valuation_laurent=`DEtools/diffop/valuation_laurent`,
	nt=`DEtools/diffop/nt`
):

# macro's for factor_op and substitute
macro(
	Newtonpolygon=`DEtools/diffop/Newtonpolygon`,
	coefs_operator=`DEtools/diffop/coefs_operator`,
	factor_newton=`DEtools/diffop/factor_newton`,
	factor_newton2=`DEtools/diffop/factor_newton2`,
	factor_op=`DEtools/diffop/factor_op`,
	factor_riccati=`DEtools/diffop/factor_riccati`,
	lift_newton=`DEtools/diffop/lift_newton`,
	lift_ramification=`DEtools/diffop/lift_ramification`,
	lift_substitute=`DEtools/diffop/lift_substitute`,
	nm_mult=`DEtools/diffop/nm_mult`,
	nm_block=`DEtools/diffop/nm_block`,
	op_with_slope=`DEtools/diffop/op_with_slope`,
	ram_laur=`DEtools/diffop/ram_laur`,
	ramification_of=`DEtools/diffop/ramification_of`,
	rem_lift=`DEtools/diffop/rem_lift`,
	skipped_factors=`DEtools/diffop/skipped_factors`,
	substitute=`DEtools/diffop/substitute`
):


#######################
#  mod p computation  #-----------------------------------------------------
#######################				diffop_modp

macro(
	compute_modm=`DEtools/diffop/compute_modm`,
	reconstruct=`DEtools/diffop/reconstruct`,
	modulus2=`DEtools/diffop/modulus2`
):

macro(
	compute_modp=`DEtools/diffop/compute_modp`,
	iratrecon=readlib('iratrecon'),
	modp_table=`DEtools/diffop/modp_table`,
	reconstruct=`DEtools/diffop/reconstruct`
):

compute_modp:=proc()
	global  accuracy_laurent, description_laurent, set_laurents,
	 value_laurent, rem_lift, modp_table;
	local i,procedure,ac_old,de_old,va_old;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if nargs=0 then # no modular computation
		accuracy_laurent:=table(modp_table[0][1]);
		description_laurent:=table(modp_table[0][2]);
		set_laurents:=modp_table[0][3];
		value_laurent:=table(modp_table[0][4]);
		rem_lift:=table(modp_table[0][5]);
		RETURN()
	fi;
	userinfo(6,'diffop',`computing modulo`,modulus);
	modp_table[0]:=[op(op(accuracy_laurent)),op(op(
	 description_laurent)),set_laurents,op(op(value_laurent))
	 ,op(op(rem_lift))];
	if modp_table[modulus]<>evaln(modp_table[modulus]) then
		accuracy_laurent:=table(modp_table[modulus][1]);
		description_laurent:=table(modp_table[modulus][2]);
		value_laurent:=table(modp_table[modulus][4]);
		rem_lift:=table(modp_table[modulus][5]);
		ac_old:=table(modp_table[0][1]);
		de_old:=table(modp_table[0][2]);
		va_old:=table(modp_table[0][4]);
		for i in set_laurents minus modp_table[modulus][3] do
			accuracy_laurent[i]:=ac_old[i];
			description_laurent[i]:=de_old[i];
			value_laurent[i]:=va_old[i]
		od;
		for i in modp_table[0][5] do if rem_lift[op(i)[1]]
		 =evaln(rem_lift[op(i)[1]]) then
			rem_lift[op(i)[1]]:=op(i)[2]
		fi od
	else
		compute_modp()
	fi;
	procedure:=args[1];
	i:=traperror(procedure(args[2..nargs]));
	if i=lasterror then
		userinfo(2,'diffop',modulus,`gives error`,i);
		compute_modp();
		RETURN(`new prime needed`)
	else
		modp_table[modulus]:=[op(op(accuracy_laurent)),op(op(
		 description_laurent)),set_laurents,op(op(value_laurent))
		 ,op(op(rem_lift))]
	fi;
	compute_modp();
	i
end:


# Compute modulo prime power.
compute_modm:=proc(exponent)
	global  modulus, modulus2;
	local v;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	modulus:=modulus2^exponent;
	v:=compute_modp(args[2..nargs]);
	while v=`new prime needed` do
		modulus2:=nextprime(modulus2*5);
		modulus:=modulus2^exponent;
		v:=compute_modp(args[2..nargs])
	od;
	modulus:=0;
	v
end:

# Input: v is a list of polynomials, with coefficients in Z/(modl Z).
# Output: a corresponding list with rational number coefficients.
reconstruct:=proc(v,modl,N)
	local w,i,t,r;
option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
if nargs=2 then
	reconstruct(args,floor(evalf(sqrt(modl/5))))
	# 5 instead of 2, because that way we can be almost sure that if
	# an answer is returned, it will be the right answer. That way we
	# don't compute a rightdivision in vain.
elif type(v,list) then
	w:=NULL;
	for i in v do
		w:=w,reconstruct(i,modl,N);
		if w[nops([w])]='failed' then RETURN('failed') fi
	od;
	[w]
elif indets(v,'name')<>{} then
	t:=[op(indets(v,'name'))];
	t:=t[1];
	w:=0;
	for i from ldegree(v,t) to degree(v,t) do
		r:=reconstruct(coeff(v,t,i),modl,N);
		if r='failed' then RETURN(r) fi;
		w:=w+r*t^i
	od;
	w
elif iratrecon(v,modl,N,N,'r','t') then
	eval(r/t)
else
	'failed'
fi
end:

#savelib('compute_modp','compute_modm','reconstruct',\
  '`DEtools/diffop/compute_modp.m`'):