1)abbrev package MULTSQFR MultivariateSquareFree
2++Author : P.Gianni
3++ This package provides the functions for the computation of the square
4++ free decomposition of a multivariate polynomial.
5++ It uses the package GenExEuclid for the resolution of
6++ the equation \spad{Af + Bg = h} and its generalization to n polynomials
7++ over an integral domain and the package \spad{MultivariateLifting}
8++ for the "multivariate" lifting.
9
10MultivariateSquareFree (E, OV, R, P) : C == T where
11 Z   ==> Integer
12 NNI ==> NonNegativeInteger
13 R   : EuclideanDomain
14 OV  : OrderedSet
15 E   : OrderedAbelianMonoidSup
16 P   : PolynomialCategory(R, E, OV)
17 SUP ==> SparseUnivariatePolynomial P
18
19 BP          ==> SparseUnivariatePolynomial(R)
20 fUnion      ==> Union("nil","sqfr","irred","prime")
21 ffSUP       ==> Record(flag : fUnion, factor : SUP, exponent : NNI)
22 ffP         ==> Record(flag : fUnion, factor : P, exponent : NNI)
23 FFE         ==> Record(factor : BP, exponent : NNI)
24 FFEP        ==> Record(factor : P, exponent : NNI)
25 FFES        ==> Record(factor : SUP, exponent : NNI)
26 Choice      ==> Record(upol : BP, Lval : List(R), Lfact : List FFE, ctpol : R)
27 squareForm  ==> Record(unitPart : P, suPart : List FFES)
28 Twopol      ==> Record(pol : SUP, polval : BP)
29 UPCF2       ==> UnivariatePolynomialCategoryFunctions2
30
31
32 C == with
33   squareFree      :      P     -> Factored P
34     ++ squareFree(p) computes the square free
35     ++ decomposition of a multivariate polynomial p.
36   squareFree      :      SUP     -> Factored SUP
37     ++ squareFree(p) computes the square free
38     ++ decomposition of a multivariate polynomial p presented as
39     ++ a univariate polynomial with multivariate coefficients.
40
41
42                    ----  local functions  ----
43   compdegd   :             List FFE                   -> Z
44        ++ compdegd should be local
45   univcase   :              (P, OV)                    -> Factored(P)
46        ++ univcase should be local
47   consnewpol :            (SUP, BP, Z)                  -> Twopol
48        ++ consnewpol should be local
49   nsqfree    :           (SUP, List(OV), List List R)  -> squareForm
50        ++ nsqfree should be local
51   intChoose  :           (SUP, List(OV), List List R)   -> Choice
52        ++ intChoose should be local
53   coefChoose : (Z, P, List FFEP) -> P
54        ++ coefChoose should be local
55   check      :     (List(FFE), List(FFE))              -> Boolean
56        ++ check should be local
57   lift : (SUP, BP, BP, P, List(OV), List(NNI), List(R), R
58          ) -> Union(List(SUP), "failed")
59        ++ lift should be local
60   myDegree   :       (SUP, List OV, NNI)                -> List NNI
61        ++ myDegree should be local
62   normDeriv2 :            (BP, Z)                      ->  BP
63        ++ normDeriv2 should be local
64
65
66
67 T == add
68
69   -- Make sure that pmod is prime over gaussian integers
70   pmod : Integer := 8388619
71
72
73   import from MultivariateLifting(E, OV, R, P)
74   import from PolynomialGcdPackage(E, OV, R, P)
75   import from FactoringUtilities(E, OV, R, P)
76   import from IntegerCombinatoricFunctions(Z)
77
78   next_mod(m : Integer) : Integer ==
79       repeat
80           m := nextPrime(m)$IntegerPrimesPackage(Integer)
81           if (m - 3) rem 4 = 0 then return m
82
83   ----  new square-free algorithm for primitive polynomial  ----
84   nsqfree1(oldf : SUP, lvar : List(OV), ltry : List List R,
85            npmod : Integer
86          ) : squareForm ==
87     f := oldf
88     univPol := intChoose(f, lvar, ltry)
89--     debug msg
90--     if not empty? ltry then output("ltry =", (ltry::OutputForm))$OutputPackage
91     f0 := univPol.upol
92     --the polynomial is square-free
93     f0 = 1$BP => [1$P, [[f, 1]$FFES]]$squareForm
94     lfact : List(FFE) := univPol.Lfact
95     lval := univPol.Lval
96     ctf := univPol.ctpol
97     leadpol : Boolean := false
98     sqdec : List FFES := empty()
99     exp0 : NNI := 0
100     unitsq : P := 1
101     lcf : P := leadingCoefficient f
102     if ctf ~= 1 then
103       f0 := ctf*f0
104       f := (ctf::P)*f
105       lcf := ctf*lcf
106     sqlead : List FFEP := empty()
107     sqlc : Factored P := 1
108     if lcf ~= 1$P then
109       leadpol := true
110       sqlc := squareFree lcf
111       unitsq := unitsq*(unit sqlc)
112       sqlead := factors sqlc
113     lfact := sort((z1 : FFE, z2 : FFE) : Boolean +-> z1.exponent > z2.exponent, lfact)
114     while lfact ~= [] repeat
115       pfact := lfact.first
116       (g0, exp0) := (pfact.factor, pfact.exponent)
117       lfact := lfact.rest
118       lfact=[] and exp0 =1 =>
119         f := (f exquo (ctf::P))::SUP
120         gg := unitNormal leadingCoefficient f
121         sqdec := cons([gg.associate*f, exp0], sqdec)
122         return  [gg.unit, sqdec]$squareForm
123       if ctf ~= 1 then g0 := ctf*g0
124       npol := consnewpol(f, f0, exp0)
125       (d, d0) := (npol.pol, npol.polval)
126       if leadpol then lcoef := coefChoose(exp0, unitsq, sqlead)
127       else lcoef := 1$P
128       ldeg := myDegree(f, lvar, exp0)
129       result := lift(d, g0, (d0 exquo g0)::BP, lcoef, lvar, ldeg, lval,
130                      npmod::R)
131       result case "failed" =>
132           return nsqfree1(oldf, lvar, ltry, next_mod(npmod))
133       result0 : SUP := (result::List SUP).1
134       r1 : SUP := result0^exp0
135       if (h := f exquo r1) case "failed" then
136           return nsqfree1(oldf, lvar, [], next_mod(npmod))
137       sqdec := cons([result0, exp0], sqdec)
138       f := h::SUP
139       f0 := completeEval(h, lvar, lval)
140       lcr : P := leadingCoefficient result0
141       nsqlead := sqlead
142       for lpfact in sqlead  while lcr ~= 1 repeat
143           ground? lcr =>
144             unitsq := (unitsq exquo lcr)::P
145             lcr := 1$P
146           repeat
147               lpfact.exponent < exp0 => break
148               (h1 := lcr exquo lpfact.factor) case P =>
149                   lcr := h1::P
150                   lpfact.exponent := (lpfact.exponent - exp0)::NNI
151               hh1 := gcd(lcr, lpfact.factor)
152               hh1 = 1$P => break
153               lcr := (lcr exquo hh1)::P
154               nff := (lpfact.factor exquo hh1)::P
155               nsqlead := cons([nff, lpfact.exponent]$FFEP, nsqlead)
156               lpfact.factor := hh1
157               lpfact.exponent := (lpfact.exponent - exp0)::NNI
158       sqlead := nsqlead
159     [((retract f) exquo ctf)::P, sqdec]$squareForm
160
161   nsqfree(oldf : SUP, lvar : List(OV), ltry : List List R) : squareForm ==
162       nsqfree1(oldf, lvar, ltry, pmod)
163
164
165   squareFree(f : SUP) : Factored SUP ==
166     degree f =0 =>
167       fu := squareFree retract f
168       makeFR((unit fu)::SUP, [["sqfr", ff.factor::SUP, ff.exponent]
169               for ff in factorList fu])
170     lvar := "setUnion"/[variables cf for cf in coefficients f]
171     empty? lvar =>  -- the polynomial is univariate
172       upol : BP := map(ground, f)$UPCF2(P, SUP, R, BP)
173       usqfr := squareFree upol
174       makeFR(map(coerce, unit usqfr)$UPCF2(R, BP, P, SUP),
175              [["sqfr", map(coerce, ff.factor)$UPCF2(R, BP, P, SUP), ff.exponent]
176                 for ff in factorList usqfr])
177
178     lcf := content f
179     f := (f exquo lcf) ::SUP
180     lcSq := squareFree lcf
181     lfs : List ffSUP := [["sqfr", ff.factor::SUP, ff.exponent]
182                          for ff in factorList lcSq]
183     partSq := nsqfree(f, lvar, empty())
184
185     lfs := append([["sqfr",fu.factor,fu.exponent]$ffSUP
186                    for fu in partSq.suPart], lfs)
187     makeFR((unit lcSq * partSq.unitPart) ::SUP, lfs)
188
189   squareFree(f : P) : Factored P ==
190     ground? f => makeFR(f, [])      ---   the polynomial is constant  ---
191
192     lvar : List(OV) := variables(f)
193     result1 : List ffP := empty()
194
195     lmdeg := minimumDegree(f, lvar)     ---       is the mindeg > 0 ? ---
196     p : P := 1$P
197     for im in 1..#lvar repeat
198       (n := lmdeg.im)=0 => "next im"
199       y := lvar.im
200       p := p*monomial(1$P, y, n)
201       result1 := cons(["sqfr",y::P,n],result1)
202     if p ~= 1$P then
203       f := (f exquo p)::P
204       if ground? f then return makeFR(f, result1)
205       lvar := variables(f)
206
207
208     #lvar = 1 =>                    ---  the polynomial is univariate ---
209       result := univcase(f, lvar.first)
210       makeFR(unit result, append(result1, factorList result))
211
212     ldeg := degree(f, lvar)          ---  general case ---
213     m := "min"/[j for j in ldeg|j ~= 0]
214     i : Z := 1
215     for j in ldeg while j>m repeat i := i+1
216     x := lvar.i
217     lvar := delete(lvar, i)
218     f0 := univariate (f, x)
219     lcont : P := content f0
220     nsqfftot := nsqfree((f0 exquo lcont)::SUP, lvar, empty())
221     nsqff:List ffP := [["sqfr",multivariate(fu.factor,x),fu.exponent]$ffP
222                          for fu in nsqfftot.suPart]
223     result1 := append(result1, nsqff)
224     ground? lcont => makeFR(lcont*nsqfftot.unitPart, result1)
225     sqlead := squareFree(lcont)
226     makeFR(unit sqlead*nsqfftot.unitPart, append(result1, factorList sqlead))
227
228  -- Choose the integer for the evaluation.                        --
229  -- If the polynomial is square-free the function returns upol=1. --
230
231   intChoose(f : SUP, lvar : List(OV), ltry : List List R) : Choice ==
232     degf := degree f
233     try_n : NNI := 0
234     nvr := #lvar
235     range : Z := 10
236     lfact1 : List(FFE) := []
237     lval1 : List R := []
238     lfact : List(FFE)
239     ctf1 : R := 1
240     f1 : BP := 1$BP
241     d1 : Z
242     while range<10000000000 repeat
243       range := 2*range
244       lval := [ran(range) for i in 1..nvr]
245       member?(lval,ltry) => "new integer"
246       ltry := cons(lval, ltry)
247       f0 := completeEval(f, lvar, lval)
248       degree f0 ~= degf  => "new integer"
249       ctf := content f0
250       lfact : List(FFE) := factors(squareFree((f0 exquo (ctf::R)::BP)::BP))
251
252            ----  the univariate polynomial is square-free  ----
253       if #lfact = 1 and (lfact.1).exponent = 1 then
254         return [1$BP, lval, lfact, 1$R]$Choice
255
256       d0 := compdegd lfact
257                 ----      inizialize lfact1      ----
258       try_n = 0 =>
259         f1 := f0
260         lfact1 := lfact
261         ctf1 := ctf
262         lval1 := lval
263         d1 := d0
264         try_n := 1
265       d0 = d1 =>
266         return [f1, lval1, lfact1, ctf1]$Choice
267       d0 < d1 =>
268         try_n := 1
269         f1 := f0
270         lfact1 := lfact
271         ctf1 := ctf
272         lval1 := lval
273         d1 := d0
274
275
276        ----  Choose the leading coefficient for the lifting  ----
277   coefChoose(exp : Z, unitsq : P, sqlead : List FFEP) : P ==
278     lcoef := unitsq
279     for term in sqlead repeat
280       texp := term.exponent
281       texp<exp => "next term"
282       texp = exp => lcoef := lcoef*term.factor
283       lcoef := lcoef*(term.factor)^((texp quo exp)::NNI)
284     lcoef
285
286        ----  Construction of the polynomials for the lifting  ----
287   consnewpol(g : SUP, g0 : BP, deg : Z) : Twopol ==
288     deg = 1 => [g, g0]$Twopol
289     deg := deg-1
290     [normalDeriv(g, deg), normDeriv2(g0, deg)]$Twopol
291
292         ----  lift the univariate square-free factor  ----
293   lift(ud : SUP, g0 : BP, g1 : BP, lcoef : P, lvar : List(OV),
294        ldeg : List(NNI), lval : List(R), npmod : R
295       ) : Union(List SUP,"failed") ==
296     leadpol : Boolean := false
297     lcd : P := leadingCoefficient ud
298     leadlist : List(P) := empty()
299
300     if not ground?(leadingCoefficient ud) then
301       leadpol := true
302       ud := lcoef*ud
303       lcg0 : R := leadingCoefficient g0
304       if ground? lcoef then g0 := (retract(lcoef) exquo lcg0)::R *g0
305       else g0 := (retract(eval(lcoef, lvar, lval)) exquo lcg0)::R * g0
306       g1 := lcg0*g1
307       leadlist := [lcoef, lcd]
308     plist := lifting(ud, lvar, [g0, g1], lval, leadlist, ldeg, npmod)
309     plist case "failed" => "failed"
310     (p0 : SUP, p1 : SUP) := ((plist::List SUP).1, (plist::List SUP).2)
311     if completeEval(p0, lvar, lval) ~= g0 then (p0, p1) := (p1, p0)
312     [primitivePart p0, primitivePart p1]
313
314                ----  the polynomial is univariate  ----
315   univcase(f : P, x : OV) : Factored(P) ==
316     uf := univariate f
317     cf := content uf
318     uf := (uf exquo cf)::BP
319     result : Factored BP := squareFree uf
320     makeFR(multivariate(cf*unit result, x),
321         [["sqfr",multivariate(term.factor,x),term.exponent]
322           for term in factorList result])
323
324   compdegd(lfact : List(FFE)) : Z ==
325     ris : Z := 0
326     for pfact in lfact repeat
327       ris := ris+(pfact.exponent -1)*degree pfact.factor
328     ris
329
330   normDeriv2(f : BP, m : Z) : BP ==
331     (n1 : Z := degree f) < m => 0$BP
332     n1 = m => (leadingCoefficient f)::BP
333     k := binomial(n1, m)
334     ris : BP := 0$BP
335     n : Z := n1
336     while n>= m repeat
337       while n1>n repeat
338         k := (k*(n1-m)) quo n1
339         n1 := n1-1
340       ris := ris+monomial(k*leadingCoefficient f, (n-m)::NNI)
341       f := reductum f
342       n := degree f
343     ris
344
345   myDegree(f : SUP, lvar : List OV, exp : NNI) : List NNI==
346     [n quo exp for n in degree(f, lvar)]
347
348--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
349--All rights reserved.
350--
351--Redistribution and use in source and binary forms, with or without
352--modification, are permitted provided that the following conditions are
353--met:
354--
355--    - Redistributions of source code must retain the above copyright
356--      notice, this list of conditions and the following disclaimer.
357--
358--    - Redistributions in binary form must reproduce the above copyright
359--      notice, this list of conditions and the following disclaimer in
360--      the documentation and/or other materials provided with the
361--      distribution.
362--
363--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
364--      names of its contributors may be used to endorse or promote products
365--      derived from this software without specific prior written permission.
366--
367--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
368--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
369--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
370--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
371--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
372--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
373--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
374--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
375--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
376--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
377--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
378