BeginPackage["Polynom`"] (* Copyright: Erich Kaltofen, 4/27/1990. Permission to use provided the copyright notice is not removed. *) PTrace::usage = "Type \"Polynom`PTrace[s__]:=Print[s]\" to trace this package." MonicLC::usage= "MonicLC[f_, x_] returns {f/lc[f, x], lc[f, x]}, the latter being the the leading coefficient of f in x. For f===0 it returns {0, 1}." MakeMonic::usage= "MakeMonic[f_, x_] returns the monic associate of f with respect to the variable x." MakeMonicMod::usage= "MakeMonic[f_, x_, p_] returns the monic associate of f with respect to the variable x modulo the prime p." RandomPoly::usage= "RandomPoly[d_, s_, x_] returns a polynomial in x of degree d, such that each coefficient is a random integer between -s to s." PlotRoots::usage= "PlotRoots[f_, x_] plots the complex roots of the polynomial f in x." RootBound::usage= "RootBound[f_, x_] returns an integer bound b such that all roots r of the real or complex polynomial f in x satisfy Abs[r] < b." FactorCoefficientBound::usage= "FactorCoefficientBound[f_, x_] gives a bound for the coefficient magnitude of all but one integer polynomial factors of the integer polynomial f in x." ModBalanced::usage= "ModBalanced[f_, x_, p_] balances the coefficients of f in x w.r.t. the modulus p." Begin["`private`"] MakeMonic[f_, x_]:= If[f === 0, f, Expand[f/Coefficient[f, x, Exponent[f, x] ] ] ] MakeMonicMod[f_, x_, p_]:= Block[{fmodp, lc, lcinv}, fmodp = PolynomialMod[f, p]; lc = Part[MonicLC[fmodp, x], 2]; lcinv = PowerMod[lc, -1, p]; Return[ PolynomialMod[ Expand[ lcinv * f ], p ] ] ] MonicLC[f_, x_]:= Block[{lc}, If[f === 0, List[f, 1], lc = Coefficient[f, x, Exponent[f, x] ]; List[Expand[f/lc], lc] ] ] ListRoots[E_]:= Block[{i, L = {}}, For[i=1,i<=Length[E],i++, L=Append[L,{Re[ E[[i,2]] ],Im[ E[[i,2]] ]}] ]; Return[L] ] RandomPoly[d_,s_,x_]:= Block[{G,i}, G = 0; For[i=0, i<=d, i++, G = G + Random[Integer,{-s,s}] * x^i ]; Return[G] ] PlotRoots[f_,x_,opts___]:=ListPlot[ListRoots[ NRoots[f==0, x] ], opts, PlotRange->{{-1.5, 1.5}, {-1.5, 1.5}}, AspectRatio->1, PlotStyle->{PointSize[0.03]} ] RootBound[f_, x_]:= Block[{coeffs, lc, B, trav, n, newB}, coeffs = Reverse[ CoefficientList[f, x] ]; lc = Abs[ First[coeffs] ]; coeffs = Rest[coeffs]; coeffs = Map[Abs, coeffs]; coeffs = Map[Ceiling, coeffs]; (* speeds up Horner eval *) (* go from B = 1 until bound is reached *) B = 1; rhs = 1; n = Length[coeffs]; While[True, (* Horner evaluate the reductum of f *) sum = 0; trav = coeffs; While[trav != {}, sum = sum * B + First[trav]; trav = Rest[ trav ] ]; lhs = sum / lc; (* evaluation of lead term *) If[rhs <= lhs, (* double B and try again *) B = 2 B; rhs = 2^n rhs, (* else *) Break[] ] ]; Return[B] ] (* end RootBound *) FactorCoefficientBound[f_, x_]:= Block[{coeflist, deg, norm2, (* || f || ^ 2 *) prec, (* precision for numerical results *) lc, (* absolute value of leading coefficient of f *) tc, (* absolute value of trailing coefficient of f *) BMlc (* bound for measure(f) * lc *) }, coeflist = CoefficientList[f, x]; deg = Length[coeflist] - 1; lc = Abs[Last[coeflist] ]; tc = Abs[First[coeflist] ]; PTrace[StringForm["deg = ``, coeflist = ``, tc = ``, lc = ``", deg, coeflist, tc, lc] ]; coeflist = Map[#^2&, coeflist]; norm2 = Apply[Plus, coeflist]; prec = Ceiling[ N[Log[10, norm2] ] ]; (* See Knuth 1981, answer to Exercise 20 of 4.6.2 *) BMlc = Ceiling[ N[ Sqrt[ 1/2 * (norm2 + N[ Sqrt[norm2^2 - 4 * lc^2 * tc^2], prec]) ], prec ] ]; (* This bound is based on the fact that if f = g[[1]]...g[[r]] then height[g[[1]] ] ... height[g[[r]] ] <= Binomial[ deg[f], deg[f]/2 ] * BMlc *) bound = Ceiling[ N[ Sqrt[ Binomial[deg, Floor[deg/2] ] * BMlc ], prec] ]; PTrace[StringForm["norm2 = ``, prec = ``, BMlc = ``, bound = ``", norm2, prec, BMlc, bound ] ]; Return[bound] ] (* end FactorCoefficientBound *) ModBalanced[f_, x_, p_]:= Block[{fmodp, coeflist, newlist, newpoly, i}, fmodp = Mod[f, p]; (* just to be save *) coeflist = CoefficientList[f, x]; coeflist = Map[ Function[ Mod[#, p] ], coeflist]; newlist = Map[ Function[If[#>(p-1)/2, #-p, #] ], coeflist]; newpoly = 0; i = 0; While[newlist != {}, newpoly = newpoly + First[newlist] * x^i; newlist = Rest[newlist]; i = i+1 ]; Return[newpoly] ] (* end ModBalanced *) End[] (* `private` *) SwinnertonDyer::usage= "SwinnertonDyer[x_, r_, p_, opts___] returns f[r; p[[1]],...,p[[n]] ] or f^*[r; p[[1]],...,p[[n]] ], the generalized Swinnerton-Dyer polynomials [Kaltofen et el, SIAM J. Comput. 12, pp. 473-485 (1983)], as a polynomial in x. opts are: Starred->True and Modulus->p" TwoGroupPoly::usage= "TwoGroupPoly[x_, k_] returns a defining equation for Sqrt[1+Sqrt[1+...+Sqrt[1+Sqrt[2] ] ] ] with k-1 1s." SpecialWrite::usage= "Austin knows!" Begin["`special`"] Options[SwinnertonDyer] = {Global`Starred->False, Global`Modulus->0} SwinnertonDyer[x_, r_, p_, opts___]:= Block[{flag, SD1, SD1subs, y, parity, g, SD}, flag = Global`Starred/.{opts}/.Options[SwinnertonDyer]; If[p == {}, If[flag, Return[Cyclotomic[r, x] ], Return[x] ] ]; SD1 = SwinnertonDyer[x, r, Drop[p, -1], opts]; parity = (-1)^(r * Exponent[SD1, x]); SD1subs = Expand[ SD1/.x->(x-y) ]; g = y^r - Last[p]; SD = Resultant[SD1subs, g, y, System`Modulus->(Global`Modulus/.{opts}/.Options[SwinnertonDyer]) ]; Return[ Expand[parity * SD] ] ] (* end SwinnertonDyer *) SpecialWrite[plist_,mod_]:= Block[{stream,poly,coeffs,i}, stream = OpenWrite["p.txt"]; poly = SwinnertonDyer[x, 2, plist, Global`Modulus->mod]; Write[stream, mod]; WriteString[stream, "["]; coeffs = CoefficientList[poly,x]; For[i=1,i<=Length[coeffs],i++, WriteString[stream, StringForm["`` ", coeffs[[i]] ] ]; ]; WriteString[stream, "]"]; Close[stream] ] (* end SpecialWrite *) TwoGroupPoly[x_, k_]:= Block[{i, f}, f = x; For[i=1, i<=k-1, i++, f = Expand[f^2 - 1] ]; Return[Expand[f^2 - 2] ] ] (* end TwoGroupPoly *) End[] (* `special` *) EndPackage[] (* Polynom` *) Print["Package Polynom` loaded"]