aboutsummaryrefslogtreecommitdiff
path: root/src/algebra/pgcd.spad.pamphlet
blob: f06e27af03ac18ea8222a7694ec7519ff03b1366 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra pgcd.spad}
\author{Michael Lucks, Patrizia Gianni, Barry Trager, Frederic Lehobey}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package PGCD PolynomialGcdPackage}
\subsection{failure case in gcdPrimitive(p1:SUPP,p2:SUPP) : SUPP}
Barry Trager has discovered and fixed a bug in pgcd.spad which occasionally
(depending on choices of random substitutions) could return the 
wrong gcd. The fix is simply to comment out one line in 
$gcdPrimitive$ which was causing the division test to be skipped.
The line used to read:
\begin{verbatim}
        degree result=d1 => result
\end{verbatim}
but has now been removed.
<<bmtfix>>=
@

<<package PGCD PolynomialGcdPackage>>=
)abbrev package PGCD PolynomialGcdPackage
++ Author: Michael Lucks, P. Gianni
++ Date Created:
++ Date Last Updated: 17 June 1996
++ Fix History: Moved unitCanonicals for performance (BMT);
++              Fixed a problem with gcd(x,0) (Frederic Lehobey)
++ Basic Functions: gcd, content
++ Related Constructors: Polynomial
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   This package computes multivariate polynomial gcd's using
++ a hensel lifting strategy. The contraint on the coefficient
++ domain is imposed by the lifting strategy. It is assumed that
++ the coefficient domain has the property that almost all specializations
++ preserve the degree of the gcd.
 
I        ==> Integer
NNI      ==> NonNegativeInteger
PI       ==> PositiveInteger
 
PolynomialGcdPackage(E,OV,R,P):C == T where
    R     :  EuclideanDomain
    P     :  PolynomialCategory(R,E,OV)
    OV    :  OrderedSet
    E     :  OrderedAbelianMonoidSup
 
    SUPP      ==> SparseUnivariatePolynomial P
 
    C == with
      gcd               :   (P,P)    -> P
        ++ gcd(p,q) computes the gcd of the two polynomials p and q.
      gcd               :   List P   -> P
        ++ gcd(lp) computes the gcd of the list of polynomials lp.
      gcd               :   (SUPP,SUPP)    -> SUPP
        ++ gcd(p,q) computes the gcd of the two polynomials p and q.
      gcd               :   List SUPP   -> SUPP
        ++ gcd(lp) computes the gcd of the list of polynomials lp.
      gcdPrimitive      :   (P,P)    -> P
        ++ gcdPrimitive(p,q) computes the gcd of the primitive polynomials
        ++ p and q.
      gcdPrimitive      :   (SUPP,SUPP)    -> SUPP
        ++ gcdPrimitive(p,q) computes the gcd of the primitive polynomials
        ++ p and q.
      gcdPrimitive      :   List P   -> P
        ++ gcdPrimitive lp computes the gcd of the list of primitive
        ++ polynomials lp.

    T == add
 
      SUP      ==> SparseUnivariatePolynomial R
 
      LGcd     ==> Record(locgcd:SUPP,goodint:List List R)
      UTerm    ==> Record(lpol:List SUP,lint:List List R,mpol:SUPP)
      pmod:R   :=  (prevPrime(2**26)$IntegerPrimesPackage(Integer))::R
 
      import MultivariateLifting(E,OV,R,P)
      import FactoringUtilities(E,OV,R,P)
 
                 --------  Local  Functions  --------
 
      myran           :    Integer   -> Union(R,"failed")
      better          :    (P,P)     -> Boolean
      failtest        :   (SUPP,SUPP,SUPP)    -> Boolean
      monomContent    :   (SUPP)   -> SUPP
      gcdMonom        :  (SUPP,SUPP)    -> SUPP
      gcdTermList     :    (P,P)     -> P
      good            :  (SUPP,List OV,List List R) -> Record(upol:SUP,inval:List List R)
 
      chooseVal       :  (SUPP,SUPP,List OV,List List R) -> Union(UTerm,"failed")
      localgcd        :  (SUPP,SUPP,List OV,List List R) -> LGcd
      notCoprime      :  (SUPP,SUPP, List NNI,List OV,List List R)   -> SUPP
      imposelc        : (List SUP,List OV,List R,List P) -> List SUP
 
      lift? :(SUPP,SUPP,UTerm,List NNI,List OV) -> Union(s:SUPP,failed:"failed",notCoprime:"notCoprime")
      lift  :(SUPP,SUP,SUP,P,List OV,List NNI,List R) -> Union(SUPP,"failed")
 
                     ---- Local  functions ----
    -- test if something wrong happened in the gcd
      failtest(f:SUPP,p1:SUPP,p2:SUPP) : Boolean ==
        (p1 exquo f) case "failed" or (p2 exquo f) case "failed"
 
    -- Choose the integers
      chooseVal(p1:SUPP,p2:SUPP,lvr:List OV,ltry:List List R):Union(UTerm,"failed") ==
        d1:=degree(p1)
        d2:=degree(p2)
        dd:NNI:=0$NNI
        nvr:NNI:=#lvr
        lval:List R :=[]
        range:I:=8
        repeat
          range:=2*range
          lval:=[ran(range) for i in 1..nvr]
          member?(lval,ltry) => "new point"
          ltry:=cons(lval,ltry)
          uf1:SUP:=completeEval(p1,lvr,lval)
          degree uf1 ~= d1 => "new point"
          uf2:SUP:= completeEval(p2,lvr,lval)
          degree uf2 ~= d2 => "new point"
          u:=gcd(uf1,uf2)
          du:=degree u
         --the univariate gcd is 1
          if du=0 then return [[1$SUP],ltry,0$SUPP]$UTerm
 
          ugcd:List SUP:=[u,(uf1 exquo u)::SUP,(uf2 exquo u)::SUP]
          uterm:=[ugcd,ltry,0$SUPP]$UTerm
          dd=0 => dd:=du
 
        --the degree is not changed
          du=dd =>
 
           --test if one of the polynomials is the gcd
            dd=d1 =>
              if not ((f:=p2 exquo p1) case "failed") then
                return [[u],ltry,p1]$UTerm
              if dd~=d2 then dd:=(dd-1)::NNI
 
            dd=d2 =>
              if not ((f:=p1 exquo p2) case "failed") then
                return [[u],ltry,p2]$UTerm
              dd:=(dd-1)::NNI
            return uterm
 
         --the new gcd has degree less
          du<dd => dd:=du
 
      good(f:SUPP,lvr:List OV,ltry:List List R):Record(upol:SUP,inval:List List R) ==
        nvr:NNI:=#lvr
        range:I:=1
        while true repeat
          range:=2*range
          lval:=[ran(range) for i in 1..nvr]
          member?(lval,ltry) => "new point"
          ltry:=cons(lval,ltry)
          uf:=completeEval(f,lvr,lval)
          if degree gcd(uf,differentiate uf)=0 then return [uf,ltry]
 
      -- impose the right lc
      imposelc(lipol:List SUP,
               lvar:List OV,lval:List R,leadc:List P):List SUP ==
        result:List SUP :=[]
        for pol in lipol for leadpol in leadc repeat
           p1:= univariate eval(leadpol,lvar,lval) * pol
           result:= cons((p1 exquo leadingCoefficient pol)::SUP,result)
        reverse result
 
    --Compute the gcd between not coprime polynomials
      notCoprime(g:SUPP,p2:SUPP,ldeg:List NNI,lvar1:List OV,ltry:List List R) : SUPP ==
        g1:=gcd(g,differentiate g)
        l1 := (g exquo g1)::SUPP
        lg:LGcd:=localgcd(l1,p2,lvar1,ltry)
        (l,ltry):=(lg.locgcd,lg.goodint)
        lval:=ltry.first
        p2l:=(p2 exquo l)::SUPP
        (gd1,gd2):=(l,l)
        ul:=completeEval(l,lvar1,lval)
        dl:=degree ul
        if degree gcd(ul,differentiate ul) ~=0 then
          newchoice:=good(l,lvar1,ltry)
          ul:=newchoice.upol
          ltry:=newchoice.inval
          lval:=ltry.first
        ug1:=completeEval(g1,lvar1,lval)
        ulist:=[ug1,completeEval(p2l,lvar1,lval)]
        lcpol:List P:=[leadingCoefficient g1, leadingCoefficient p2]
        while true repeat
          d:SUP:=gcd(cons(ul,ulist))
          if degree d =0 then return gd1
          lquo:=(ul exquo d)::SUP
          if degree lquo ~=0 then
            lgcd:=gcd(cons(leadingCoefficient l,lcpol))
            (gdl:=lift(l,d,lquo,lgcd,lvar1,ldeg,lval)) case "failed" =>
               return notCoprime(g,p2,ldeg,lvar1,ltry)
            l:=gd2:=gdl::SUPP
            ul:=completeEval(l,lvar1,lval)
            dl:=degree ul
          gd1:=gd1*gd2
          ulist:=[(uf exquo d)::SUP for uf in ulist]
 
      gcdPrimitive(p1:SUPP,p2:SUPP) : SUPP ==
        if (d1:=degree(p1)) > (d2:=degree(p2)) then
            (p1,p2):= (p2,p1)
            (d1,d2):= (d2,d1)
        degree p1 = 0 =>
            p1 = 0 => unitCanonical p2
            unitCanonical p1
        lvar:List OV:=sort(#1>#2,setUnion(variables p1,variables p2))
        empty? lvar =>
           raisePolynomial(gcd(lowerPolynomial p1,lowerPolynomial p2))
        (p2 exquo p1) case SUPP => unitCanonical p1
        ltry:List List R:=empty()
        totResult:=localgcd(p1,p2,lvar,ltry)
        result: SUPP:=totResult.locgcd
        -- special cases
        result=1 => 1$SUPP
<<bmtfix>>
        while failtest(result,p1,p2) repeat
--        SAY$Lisp  "retrying gcd"
          ltry:=totResult.goodint
          totResult:=localgcd(p1,p2,lvar,ltry)
          result:=totResult.locgcd
        result
 
    --local function for the gcd : it returns the evaluation point too
      localgcd(p1:SUPP,p2:SUPP,lvar:List(OV),ltry:List List R) : LGcd ==
        uterm:=chooseVal(p1,p2,lvar,ltry)::UTerm
        ltry:=uterm.lint
        listpol:= uterm.lpol
        ud:=listpol.first
        dd:= degree ud
 
        --the univariate gcd is 1
        dd=0 => [1$SUPP,ltry]$LGcd
 
        --one of the polynomials is the gcd
        dd=degree(p1) or dd=degree(p2) =>
                         [uterm.mpol,ltry]$LGcd
        ldeg:List NNI:=map(min,degree(p1,lvar),degree(p2,lvar))
 
       -- if there is a polynomial g s.t. g/gcd and gcd are coprime ...
        -- I can lift
        (h:=lift?(p1,p2,uterm,ldeg,lvar)) case notCoprime =>
          [notCoprime(p1,p2,ldeg,lvar,ltry),ltry]$LGcd
        h case failed => localgcd(p1,p2,lvar,ltry) -- skip bad values?
        [h.s,ltry]$LGcd
 
 
  -- content, internal functions return the poly if it is a monomial
      monomContent(p:SUPP):SUPP ==
        degree(p)=0 => 1
        md:= minimumDegree(p)
        monomial(gcd sort(better,coefficients p),md)
 
  -- Ordering for gcd purposes
      better(p1:P,p2:P):Boolean ==
        ground? p1 => true
        ground? p2 => false
        degree(p1,mainVariable(p1)::OV) < degree(p2,mainVariable(p2)::OV)
 
  -- Gcd between polynomial p1 and p2 with
  -- mainVariable p1 < x=mainVariable p2
      gcdTermList(p1:P,p2:P) : P ==
        termList:=sort(better,
           cons(p1,coefficients univariate(p2,(mainVariable p2)::OV)))
        q:P:=termList.first
        for term in termList.rest until q = 1$P repeat q:= gcd(q,term)
        q
 
  -- Gcd between polynomials with the same mainVariable
      gcd(p1:SUPP,p2:SUPP): SUPP ==
        if degree(p1) > degree(p2) then (p1,p2):= (p2,p1)
        degree p1 = 0 =>
           p1 = 0 => unitCanonical p2
           p1 = 1 => unitCanonical p1
           gcd(leadingCoefficient p1, content p2)::SUPP
        reductum(p1)=0 => gcdMonom(p1,monomContent p2)
        c1:= monomContent(p1)
        reductum(p2)=0 => gcdMonom(c1,p2)
        c2:= monomContent(p2)
        p1:= (p1 exquo c1)::SUPP
        p2:= (p2 exquo c2)::SUPP
        gcdPrimitive(p1,p2) * gcdMonom(c1,c2)
 
   -- gcd between 2 monomials
      gcdMonom(m1:SUPP,m2:SUPP):SUPP ==
        monomial(gcd(leadingCoefficient(m1),leadingCoefficient(m2)),
                 min(degree(m1),degree(m2)))

 
    --If there is a pol s.t. pol/gcd and gcd are coprime I can lift
      lift?(p1:SUPP,p2:SUPP,uterm:UTerm,ldeg:List NNI,
                     lvar:List OV) : Union(s:SUPP,failed:"failed",notCoprime:"notCoprime") ==
        leadpol:Boolean:=false
        (listpol,lval):=(uterm.lpol,uterm.lint.first)
        d:=listpol.first
        listpol:=listpol.rest
        nolift:Boolean:=true
        uf: SUP
        for fs in tails listpol while nolift repeat
          uf := first fs
          --note uf and d not necessarily primitive
          degree gcd(uf,d) =0 => nolift:=false
        nolift => ["notCoprime"]
        f:SUPP:=([p1,p2]$List(SUPP)).(position(uf,listpol))
        lgcd:=gcd(leadingCoefficient p1, leadingCoefficient p2)
        (l:=lift(f,d,uf,lgcd,lvar,ldeg,lval)) case "failed" => ["failed"]
        [l :: SUPP]
 
   -- interface with the general "lifting" function
      lift(f:SUPP,d:SUP,uf:SUP,lgcd:P,lvar:List OV,
                        ldeg:List NNI,lval:List R):Union(SUPP,"failed") ==
        leadpol:Boolean:=false
        lcf:P
        lcf:=leadingCoefficient f
        df:=degree f
        leadlist:List(P):=[]
 
        if lgcd~=1 then
          leadpol:=true
          f:=lgcd*f
          ldeg:=[n0+n1 for n0 in ldeg for n1 in degree(lgcd,lvar)]
          lcd:R:=leadingCoefficient d
          if degree(lgcd)=0 then d:=((retract lgcd) *d exquo lcd)::SUP
          else d:=(retract(eval(lgcd,lvar,lval)) * d exquo lcd)::SUP
          uf:=lcd*uf
        leadlist:=[lgcd,lcf]
        lg:=imposelc([d,uf],lvar,lval,leadlist)
        (pl:=lifting(f,lvar,lg,lval,leadlist,ldeg,pmod)) case "failed" =>
               "failed"
        plist := pl :: List SUPP
        (p0:SUPP,p1:SUPP):=(plist.first,plist.2)
        if completeEval(p0,lvar,lval) ~= lg.first then
           (p0,p1):=(p1,p0)
        not leadpol => p0
        p0 exquo content(p0)
 
  -- Gcd for two multivariate polynomials
      gcd(p1:P,p2:P) : P ==
        ground? p1 =>
          p1 := unitCanonical p1
          p1 = 1$P => p1
          p1 = 0$P => unitCanonical p2
          ground? p2 => gcd((retract p1)@R,(retract p2)@R)::P
          gcdTermList(p1,p2)
        ground? p2 =>
          p2 := unitCanonical p2
          p2 = 1$P => p2
          p2 = 0$P => unitCanonical p1
          gcdTermList(p2,p1)
        (p1:= unitCanonical(p1)) = (p2:= unitCanonical(p2)) => p1
        mv1:= mainVariable(p1)::OV
        mv2:= mainVariable(p2)::OV
        mv1 = mv2 => multivariate(gcd(univariate(p1,mv1),
                                      univariate(p2,mv1)),mv1)
        mv1 < mv2 => gcdTermList(p1,p2)
        gcdTermList(p2,p1)
 
  -- Gcd for a list of multivariate polynomials
      gcd(listp:List P) : P ==
        lf:=sort(better,listp)
        f:=lf.first
        for g in lf.rest repeat
          f:=gcd(f,g)
          if f=1$P then return f
        f

      gcd(listp:List SUPP) : SUPP ==
        lf:=sort(degree(#1)<degree(#2),listp)
        f:=lf.first
        for g in lf.rest repeat
          f:=gcd(f,g)
          if f=1 then return f
        f

 
   -- Gcd for primitive polynomials
      gcdPrimitive(p1:P,p2:P):P ==
        (p1:= unitCanonical(p1)) = (p2:= unitCanonical(p2)) => p1
        ground? p1 =>
          ground? p2 => gcd((retract p1)@R,(retract p2)@R)::P
          p1 = 0$P => p2
          1$P
        ground? p2 =>
          p2 = 0$P => p1
          1$P
        mv1:= mainVariable(p1)::OV
        mv2:= mainVariable(p2)::OV
        mv1 = mv2 =>
          md:=min(minimumDegree(p1,mv1),minimumDegree(p2,mv2))
          mp:=1$P
          if md>1 then
            mp:=(mv1::P)**md
            p1:=(p1 exquo mp)::P
            p2:=(p2 exquo mp)::P
          up1 := univariate(p1,mv1)
          up2 := univariate(p2,mv2)
          mp*multivariate(gcdPrimitive(up1,up2),mv1)
        1$P
 
  -- Gcd for a list of primitive multivariate polynomials
      gcdPrimitive(listp:List P) : P ==
        lf:=sort(better,listp)
        f:=lf.first
        for g in lf.rest repeat
          f:=gcdPrimitive(f,g)
          if f=1$P then return f
        f

@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>
 
<<package PGCD PolynomialGcdPackage>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}