aboutsummaryrefslogtreecommitdiff
path: root/src/algebra/gpgcd.spad.pamphlet
blob: c247b8cc905d4153e27f737914a0efba8dbf754c (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
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra gpgcd.spad}
\author{The Axiom Team}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package GENPGCD GeneralPolynomialGcdPackage}
<<package GENPGCD GeneralPolynomialGcdPackage>>=
)abbrev package GENPGCD GeneralPolynomialGcdPackage
++ Description:
++ This package provides operations for GCD computations
++ on polynomials 
GeneralPolynomialGcdPackage(E,OV,R,P):C == T where
    R     :  PolynomialFactorizationExplicit
    P     :  PolynomialCategory(R,E,OV)
    OV    :  OrderedSet
    E     :  OrderedAbelianMonoidSup
 
    SUPP      ==> SparseUnivariatePolynomial P
--JHD    ContPrim  ==> Record(cont:P,prim:P)
 
    C == with
           gcdPolynomial     :   (SUPP,SUPP) -> SUPP
		++ gcdPolynomial(p,q) returns the GCD of p and q
           randomR        : ()                ->R
		++ randomR() should be local but conditional
--JHD      gcd               :   (P,P)    -> P
--JHD      gcd               :   List P   -> P
--JHD      gcdprim           :   (P,P)    -> P
--JHD      gcdprim           :   List P   -> P
 
--JHD      gcdcofact         :   List P   -> List P
--JHD      gcdcofactprim     :   List P   -> List P
 
--JHD      primitate         :   (P,OV)   -> P
--JHD      primitate         :    SUPP    -> SUPP
 
--JHD      content           :     P      -> P
--JHD      content           :   List P   -> List P
--JHD      contprim          :   List P   -> List ContPrim
 
--JHD      monomContent      :   (P,OV)   -> P
--JHD      monomContent      :    SUPP    -> SUPP
 
 
    T == add
 
      SUPR     ==> SparseUnivariatePolynomial R
--JHD      SUPLGcd  ==> Record(locgcd:SUPP,goodint:List R)
--JHD      LGcd     ==> Record(locgcd:P,goodint:List R)
--JHD      UTerm    ==> Record(lpol:List SUPR,lint:List R,mpol:P)
--JHD--JHD      pmod:R   :=  (prevPrime(2**26)$IntegerPrimesPackage(Integer))::R
 
--JHD      import MultivariateLifting(E,OV,R,P,pmod)
      import UnivariatePolynomialCategoryFunctions2(R,SUPR,P,SUPP)
      import UnivariatePolynomialCategoryFunctions2(P,SUPP,R,SUPR)
                 --------  Local  Functions  --------
 
--JHD      abs             :      P       -> P
      better          :    (P,P)     -> Boolean
--JHD      failtest        :   (P,P,P)    -> Boolean
--JHD      gcdMonom        :  (P,P,OV)    -> P
--JHD      gcdTermList     :    (P,P)     -> P
--JHD      gcdPrim         :  (P,P,OV)    -> P
--JHD      gcdSameMainvar  :  (P,P,OV)    -> P
--JHD      internal        :  (P,P,OV)    -> P
--JHD      good            :  (P,List OV) -> Record(upol:SUPR,inval:List R)
--JHD      gcdPrs          : (P,P,NNI,OV) -> Union(P,"failed")
--JHD
--JHD      chooseVal       :     (P,P,List OV)          -> UTerm
--JHD      localgcd        :     (P,P,List OV)          -> LGcd
--JHD      notCoprime      :  (P,P, List NNI,List OV)   -> P
--JHD      imposelc        : (List SUPR,List OV,List R,List P)  -> List SUPR
 
--JHD      lift? :(P,P,UTerm,List NNI,List OV)  -> Union("failed",P)
--      lift  :(P,SUPR,SUPR,P,List OV,List NNI,List R) -> P
      lift  :  (SUPR,SUPP,SUPR,List OV,List R) -> Union(SUPP,"failed")
         -- lifts first and third arguments as factors of the second
         -- fourth is number of variables.
--JHD      monomContent      :   (P,OV)   -> P
      monomContentSup   :    SUPP    -> SUPP
--
--JHD  gcdcofact         :   List P   -> List P
 
      gcdTrivial      :  (SUPP,SUPP)   -> SUPP
      gcdSameVariables:  (SUPP,SUPP,List OV)    -> SUPP
      recursivelyGCDCoefficients: (SUPP,List OV,SUPP,List OV) -> SUPP
      flatten         : (SUPP,List OV) -> SUPP
                        -- evaluates out all variables in the second
                        -- argument, leaving a polynomial of the same
                        -- degree
--    eval            : (SUPP,List OV,List R) -> SUPP
      variables       :  SUPP          -> List OV
                     ---- JHD's exported functions ---
      gcdPolynomial(p1:SUPP,p2:SUPP) ==
        zero? p1 => p2
        zero? p2 => p1
        0=degree p1  => gcdTrivial(p1,p2)
        0=degree p2  => gcdTrivial(p2,p1)
        if degree p1 < degree p2 then (p1,p2):=(p2,p1)
        p1 exquo p2 case SUPP  => (unitNormal p2).canonical
        c1:= monomContentSup(p1)
        c2:= monomContentSup(p2)
        p1:= (p1 exquo c1)::SUPP
        p2:= (p2 exquo c2)::SUPP
        (p1 exquo p2) case SUPP => (unitNormal p2).canonical * gcd(c1,c2)
        vp1:=variables p1
        vp2:=variables p2
        v1:=setDifference(vp1,vp2)
        v2:=setDifference(vp2,vp1)
        #v1 = 0 and #v2 = 0 => gcdSameVariables(p1,p2,vp1)*gcd(c1,c2)
                 -- all variables are in common
        v:=setDifference(vp1,v1)
        pp1:=flatten(p1,v1)
        pp2:=flatten(p2,v2)
        g:=gcdSameVariables(pp1,pp2,v)
        one? g => gcd(c1,c2)::SUPP
        (#v1 = 0 or not (p1 exquo g) case "failed") and
                     -- if #vi = 0 then pp1 = p1, so we know g divides
              (#v2 = 0 or not (p2 exquo g) case "failed")
            => g*gcd(c1,c2)  -- divdes them both, so is the gcd
        -- OK, so it's not the gcd: try again
        v:=variables g -- there can be at most these variables in answer
        v1:=setDifference(vp1,v)
        v2:=setDifference(vp2,v)
        if (#v1 = 0) then g:= gcdSameVariables(g,flatten(p2,v2),v)
        else if (#v2=0) then g:=gcdSameVariables(g,flatten(p1,v1),v)
        else g:=gcdSameVariables(g,flatten(p1,v1)-flatten(p2,v2),v)
        one? g => gcd(c1,c2)::SUPP
        (#v1 = 0 or not (p1 exquo g) case "failed") and
              (#v2 = 0 or not (p2 exquo g) case "failed")
            => g*gcd(c1,c2)::SUPP  -- divdes them both, so is the gcd
        v:=variables g -- there can be at most these variables in answer
        v1:=setDifference(vp1,v)
        if #v1 ~= 0 then
           g:=recursivelyGCDCoefficients(g,v,p1,v1)
           one? g => return gcd(c1,c2)::SUPP
           v:=variables g -- there can be at most these variables in answer
        v2:=setDifference(vp2,v)
        recursivelyGCDCoefficients(g,v,p2,v2)*gcd(c1,c2)
      if R has StepThrough then
         randomCount:R := init()
         randomR() ==
            (v:=nextItem(randomCount)) case R =>
                randomCount:=v
                v
            SAY("Taking next stepthrough range in GeneralPolynomialGcdPackage")$Lisp
            randomCount:=init()
            randomCount
      else
            randomR() == (random()$Integer rem 100)::R
                     ---- JHD's local functions ---
      gcdSameVariables(p1:SUPP,p2:SUPP,lv:List OV) ==
            -- two non-trivial primitive (or, at least, we don't care
            -- about content)
            -- polynomials with precisely the same degree
          #lv = 0 => map(#1::P,gcdPolynomial(map(ground,p1),
                                             map(ground,p2)))
          degree p2 = 1 =>
            p1 exquo p2 case SUPP => p2
            1
          gcdLC:=gcd(leadingCoefficient p1,leadingCoefficient p2)
          lr:=[randomR() for vv in lv]
          count:NonNegativeInteger:=0
          while count<10 repeat
	    while zero? eval(gcdLC,lv,lr) and count<10 repeat
		  lr:=[randomR() for vv in lv]
		  count:=count+1
	    count = 10 => error "too many evaluations in GCD code"
	    up1:SUPR:=map(ground eval(#1,lv,lr),p1)
	    up2:SUPR:=map(ground eval(#1,lv,lr),p2)
	    u:=gcdPolynomial(up1,up2)
	    degree u = 0 => return 1
            -- let's pick a second one, just to check
            lrr:=[randomR() for vv in lv]
	    while zero? eval(gcdLC,lv,lrr) and count<10 repeat
		  lrr:=[randomR() for vv in lv]
		  count:=count+1
	    count = 10 => error "too many evaluations in GCD code"
	    vp1:SUPR:=map(ground eval(#1,lv,lrr),p1)
	    vp2:SUPR:=map(ground eval(#1,lv,lrr),p2)
	    v:=gcdPolynomial(vp1,vp2)
	    degree v = 0 => return 1
            if degree v < degree u then
               u:=v
               up1:=vp1
               up2:=vp2
               lr:=lrr
	    up1:=(up1 exquo u)::SUPR
	    degree gcd(u,up1) = 0 =>
                ans:=lift(u,p1,up1,lv,lr)
                ans case SUPP => return ans
                "next"
	    up2:=(up2 exquo u)::SUPR
	    degree gcd(u,up2) = 0 =>
                ans:=lift(u,p2,up2,lv,lr)
                ans case SUPP => return ans
                "next"
	    -- so neither cofactor is relatively prime
	    count:=0
	    while count < 10 repeat
	       r:=randomR()
	       uu:=up1+r*up2
	       degree gcd(u,uu)=0 =>
                 ans:= lift(u,p1+r::P *p2,uu,lv,lr)
                 ans case SUPP => return ans
                 "next"
	    error "too many evaluations in GCD code"
          count >= 10 => error "too many evaluations in GCD code"
      lift(gR:SUPR,p:SUPP,cfR:SUPR,lv:List OV,lr:List R) ==
        -- lift the coprime factorisation gR*cfR = (univariate of p)
        -- where the variables lv have been evaluated at lr
        lcp:=leadingCoefficient p
        g:=monomial(lcp,degree gR)+map(#1::P,reductum gR)
        cf:=monomial(lcp,degree cfR)+map(#1::P,reductum cfR)
        p:=lcp*p       -- impose leaidng coefficient of p on each factor
        while lv ~= [] repeat
           v:=first lv
           r:=first lr
           lv:=rest lv
           lr:=rest lr
           thisp:=map(eval(#1,lv,lr),p)
           d:="max"/[degree(c,v) for c in coefficients p]
           prime:=v::P - r::P
           pn:=prime
           origFactors:=[g,cf]::List SUPP
           for n in 1..d repeat
              Ecart:=(thisp- g*cf) exquo pn
              Ecart case "failed" =>
                 error "failed lifting in hensel in Complex Polynomial GCD"
              zero? Ecart => leave
              step:=solveLinearPolynomialEquation(origFactors,
                                                  map(eval(#1,v,r),Ecart::SUPP))
              step case "failed" => return "failed"
              g:=g+pn*first step
              cf:=cf+pn*second step
              pn:=pn*prime
           thisp ~= g*cf => return "failed"
        g
      recursivelyGCDCoefficients(g:SUPP,v:List OV,p:SUPP,pv:List OV) ==
         mv:=first pv   -- take each coefficient w.r.t. mv
         pv:=rest pv    -- and recurse on pv as necessary
         d:="max"/[degree(u,mv) for u in coefficients p]
         for i in 0..d repeat
             p1:=map(coefficient(#1,mv,i),p)
             oldg:=g
             if pv = [] then g:=gcdSameVariables(g,p1,v)
             else g:=recursivelyGCDCoefficients(p,v,p1,pv)
             one? g => return 1
             g~=oldg =>
                oldv:=v
                v:=variables g
                pv:=setUnion(pv,setDifference(v,oldv))
         g
      flatten(p1:SUPP,lv:List OV) ==
         #lv = 0 => p1
         lr:=[ randomR() for vv in lv]
         dg:=degree p1
         ans : SUPP
         while dg ~= degree (ans:= map(eval(#1,lv,lr),p1)) repeat
           lr:=[ randomR() for vv in lv]
         ans
--      eval(p1:SUPP,lv:List OV,lr:List R) == map(eval(#1,lv,lr),p1)
      variables(p1:SUPP) ==
        removeDuplicates ("concat"/[variables u for u in coefficients p1])
      gcdTrivial(p1:SUPP,p2:SUPP) ==
        -- p1 is non-zero, but has degree zero
        -- p2 is non-zero
        cp1:=leadingCoefficient p1
        one? cp1 => 1
        degree p2 = 0 => gcd(cp1,leadingCoefficient p2)::SUPP
        un?:=unit? cp1
        while not zero? p2 and not un? repeat
           cp1:=gcd(leadingCoefficient p2,cp1)
           un?:=unit? cp1
           p2:=reductum p2
        un? => 1
        cp1::SUPP
 
                     ---- Local  functions ----
--JHD    -- test if something wrong happened in the gcd
--JHD      failtest(f:P,p1:P,p2:P) : Boolean ==
--JHD        (p1 exquo f) case "failed" or (p2 exquo f) case "failed"
--JHD
--JHD    -- Choose the integers
--JHD      chooseVal(p1:P,p2:P,lvar:List OV):UTerm ==
--JHD        x:OV:=lvar.first
--JHD        lvr:=lvar.rest
--JHD        d1:=degree(p1,x)
--JHD        d2:=degree(p2,x)
--JHD        dd:NNI:=0$NNI
--JHD        nvr:NNI:=#lvr
--JHD        lval:List R :=[]
--JHD        range:I:=8
--JHD        for i in 1..  repeat
--JHD          range:=2*range
--JHD          lval:=[(random()$I rem (2*range) - range)::R  for i in 1..nvr]
--JHD          uf1:SUPR:=univariate eval(p1,lvr,lval)
--JHD          degree uf1 ~= d1 => "new point"
--JHD          uf2:SUPR:=univariate eval(p2,lvr,lval)
--JHD          degree uf2 ~= d2 => "new point"
--JHD          u:=gcd(uf1,uf2)
--JHD          du:=degree u
--JHD         --the univariate gcd is 1
--JHD          if du=0 then return [[1$SUPR],lval,0$P]$UTerm
--JHD
--JHD          ugcd:List SUPR:=[u,(uf1 exquo u)::SUPR,(uf2 exquo u)::SUPR]
--JHD          uterm:=[ugcd,lval,0$P]$UTerm
--JHD          dd=0 => dd:=du
--JHD
--JHD        --the degree is not changed
--JHD          du=dd =>
--JHD
--JHD           --test if one of the polynomials is the gcd
--JHD            dd=d1 =>
--JHD              if ^((f:=p2 exquo p1) case "failed") then
--JHD                return [[u],lval,p1]$UTerm
--JHD              if dd~=d2 then dd:=(dd-1)::NNI
--JHD
--JHD            dd=d2 =>
--JHD              if ^((f:=p1 exquo p2) case "failed") then
--JHD                return [[u],lval,p2]$UTerm
--JHD              dd:=(dd-1)::NNI
--JHD            return uterm
--JHD
--JHD         --the new gcd has degree less
--JHD          du<dd => dd:=du
--JHD
--JHD      good(f:P,lvr:List OV):Record(upol:SUPR,inval:List R) ==
--JHD        nvr:NNI:=#lvr
--JHD        range:I:=1
--JHD        ltry:List List R:=[]
--JHD        while true repeat
--JHD          range:=2*range
--JHD          lval:=[(random()$I rem (2*range) -range)::R  for i in 1..nvr]
--JHD          member?(lval,ltry) => "new point"
--JHD          ltry:=cons(lval,ltry)
--JHD          uf:=univariate eval(f,lvr,lval)
--JHD          if degree gcd(uf,differentiate uf)=0 then return [uf,lval]
--JHD
--JHD      -- impose the right lc
--JHD      imposelc(lipol:List SUPR,
--JHD               lvar:List OV,lval:List R,leadc:List P):List SUPR ==
--JHD        result:List SUPR :=[]
--JHD        lvar:=lvar.rest
--JHD        for pol in lipol for leadpol in leadc repeat
--JHD           p1:= univariate eval(leadpol,lvar,lval) * pol
--JHD           result:= cons((p1 exquo leadingCoefficient pol)::SUPR,result)
--JHD        reverse result
--JHD
--JHD    --Compute the gcd between not coprime polynomials
--JHD      notCoprime(g:P,p2:P,ldeg:List NNI,lvar:List OV) : P ==
--JHD        x:OV:=lvar.first
--JHD        lvar1:List OV:=lvar.rest
--JHD        lg1:=gcdcofact([g,differentiate(g,x)])
--JHD        g1:=lg1.1
--JHD        lg:LGcd:=localgcd(g1,p2,lvar)
--JHD        (l,lval):=(lg.locgcd,lg.goodint)
--JHD        p2:=(p2 exquo l)::P
--JHD        (gd1,gd2):=(l,l)
--JHD        ul:=univariate(eval(l,lvar1,lval))
--JHD        dl:=degree ul
--JHD        if degree gcd(ul,differentiate ul) ~=0 then
--JHD          newchoice:=good(l,lvar.rest)
--JHD          ul:=newchoice.upol
--JHD          lval:=newchoice.inval
--JHD        ug1:=univariate(eval(g1,lvar1,lval))
--JHD        ulist:=[ug1,univariate eval(p2,lvar1,lval)]
--JHD        lcpol:=[leadingCoefficient univariate(g1,x),
--JHD                leadingCoefficient univariate(p2,x)]
--JHD        while true repeat
--JHD          d:SUPR:=gcd(cons(ul,ulist))
--JHD          if degree d =0 then return gd1
--JHD          lquo:=(ul exquo d)::SUPR
--JHD          if degree lquo ~=0 then
--JHD            lgcd:=gcd(cons(leadingCoefficient univariate(l,x),lcpol))
--JHD            gd2:=lift(l,d,lquo,lgcd,lvar,ldeg,lval)
--JHD            l:=gd2
--JHD            ul:=univariate(eval(l,lvar1,lval))
--JHD            dl:=degree ul
--JHD          gd1:=gd1*gd2
--JHD          ulist:=[(uf exquo d)::SUPR for uf in ulist]
--JHD
--JHD -- we suppose that the poly have the same mainvar, deg p1<deg p2 and the
--JHD -- polys primitive
--JHD      internal(p1:P,p2:P,x:OV) : P ==
--JHD        lvar:List OV:=sort(#1>#2,setUnion(variables p1,variables p2))
--JHD        d1:=degree(p1,x)
--JHD        d2:=degree(p2,x)
--JHD        result: P:=localgcd(p1,p2,lvar).locgcd
--JHD        -- special cases
--JHD        result=1 => 1$P
--JHD        (dr:=degree(result,x))=d1 or dr=d2  => result
--JHD        while failtest(result,p1,p2) repeat
--JHD          SAY$Lisp  "retrying gcd"
--JHD          result:=localgcd(p1,p2,lvar).locgcd
--JHD        result
--JHD
--JHD    --local function for the gcd : it returns the evaluation point too
--JHD      localgcd(p1:P,p2:P,lvar:List(OV)) : LGcd ==
--JHD        x:OV:=lvar.first
--JHD        uterm:=chooseVal(p1,p2,lvar)
--JHD        listpol:= uterm.lpol
--JHD        ud:=listpol.first
--JHD        dd:= degree ud
--JHD
--JHD        --the univariate gcd is 1
--JHD        dd=0 => [1$P,uterm.lint]$LGcd
--JHD
--JHD        --one of the polynomials is the gcd
--JHD        dd=degree(p1,x) or dd=degree(p2,x) =>
--JHD                                           [uterm.mpol,uterm.lint]$LGcd
--JHD        ldeg:List NNI:=map(min,degree(p1,lvar),degree(p2,lvar))
--JHD
--JHD       -- if there is a polynomial g s.t. g/gcd and gcd are coprime ...
--JHD        -- I can lift
--JHD        (h:=lift?(p1,p2,uterm,ldeg,lvar)) case "failed" =>
--JHD          [notCoprime(p1,p2,ldeg,lvar),uterm.lint]$LGcd
--JHD        [h::P,uterm.lint]$LGcd
--JHD
--JHD
--JHD  -- content, internal functions return the poly if it is a monomial
--JHD      monomContent(p:P,var:OV):P ==
--JHD        ground? p => 1$P
--JHD        md:= minimumDegree(p,var)
--JHD        ((var::P)**md)*(gcd sort(better,coefficients univariate(p,var)))
 
      monomContentSup(u:SUPP):SUPP ==
        degree(u) = 0$NonNegativeInteger => 1$SUPP
        md:= minimumDegree u
        gcd(sort(better,coefficients u)) * monomial(1$P,md)$SUPP
 
--JHD  -- change the polynomials to have positive lc
--JHD      abs(p:P): P == unitNormal(p).canonical
 
  -- 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)
 
   -- PRS algorithm
   -- gcdPrs(p1:P,p2:P,d:NNI,var:OV):Union(P,"failed") ==
   --   u1:= univariate(p1,var)
   --   u2:= univariate(p2,var)
   --   finished:Boolean:= false
   --   until finished repeat
   --     dd:NNI:=(degree u1 - degree u2)::NNI
   --     lc1:SUPP:=leadingCoefficient u2 * reductum u1
   --     lc2:SUPP:=leadingCoefficient u1 * reductum u2
   --     u3:SUPP:= primitate((lc1-lc2)*monomial(1$P,dd))$%
   --     (d3:=degree(u3)) <= d => finished:= true
   --     u1:= u2
   --     u2:= u3
   --     if d3 > degree(u1) then (u1,u2):= (u2,u1)
   --   g:= (u2 exquo u3)
   --   g case SUPP => abs multivariate(u3,var)
   --   "failed"
 
  -- Gcd between polynomial p1 and p2 with
  -- mainVariable p1 < x=mainVariable p2
--JHD      gcdTermList(p1:P,p2:P) : P ==
--JHD        termList:=sort(better,
--JHD           cons(p1,coefficients univariate(p2,(mainVariable p2)::OV)))
--JHD        q:P:=termList.first
--JHD        for term in termList.rest until q = 1$P repeat q:= gcd(q,term)
--JHD        q
--JHD
--JHD  -- Gcd between polynomials with the same mainVariable
--JHD      gcdSameMainvar(p1:P,p2:P,mvar:OV): P ==
--JHD        if degree(p1,mvar) < degree(p2,mvar) then (p1,p2):= (p2,p1)
--JHD        (p1 exquo p2) case P => abs p2
--JHD        c1:= monomContent(p1,mvar)$%
--JHD        c1 = p1 => gcdMonom(p1,p2,mvar)
--JHD        c2:= monomContent(p2,mvar)$%
--JHD        c2 = p2 => gcdMonom(p2,p1,mvar)
--JHD        p1:= (p1 exquo c1)::P
--JHD        p2:= (p2 exquo c2)::P
--JHD        if degree(p1,mvar) < degree(p2,mvar) then (p1,p2):= (p2,p1)
--JHD        (p1 exquo p2) case P => abs(p2) * gcd(c1,c2)
--JHD        abs(gcdPrim(p1,p2,mvar)) * gcd(c1,c2)
--JHD
--JHD   --  make the polynomial primitive with respect to var
--JHD      primitate(p:P,var:OV):P == (p exquo monomContent(p,var))::P
--JHD
--JHD      primitate(u:SUPP):SUPP == (u exquo monomContentSup u)::SUPP
--JHD
--JHD   -- gcd between primitive polynomials with the same mainVariable
--JHD      gcdPrim(p1:P,p2:P,mvar:OV):P ==
--JHD        vars:= removeDuplicates append(variables p1,variables p2)
--JHD        #vars=1 => multivariate(gcd(univariate p1,univariate p2),mvar)
--JHD        vars:=delete(vars,position(mvar,vars))
--JHD        --d:= degModGcd(p1,p2,mvar,vars)
--JHD        --d case "failed" => internal(p2,p1,mvar)
--JHD        --deg:= d:NNI
--JHD        --deg = 0$NNI => 1$P
--JHD        --deg = degree(p1,mvar) =>
--JHD        --  (p2 exquo p1) case P => abs(p1)  -- already know that
--JHD                                           -- ^(p1 exquo p2)
--JHD        --  internal(p2,p1,mvar)
--JHD        --cheapPrs?(p1,p2,deg,mvar) =>
--JHD        --  g:= gcdPrs(p1,p2,deg,mvar)
--JHD        --  g case P => g::P
--JHD        --  internal(p2,p1,mvar)
--JHD        internal(p2,p1,mvar)
--JHD
--JHD   -- gcd between a monomial and a polynomial
--JHD      gcdMonom(m:P,p:P,var:OV):P ==
--JHD        ((var::P) ** min(minimumDegree(m,var),minimumDegree(p,var))) *
--JHD          gcdTermList(leadingCoefficient(univariate(m,var)),p)
--JHD
--JHD    --If there is a pol s.t. pol/gcd and gcd are coprime I can lift
--JHD      lift?(p1:P,p2:P,uterm:UTerm,ldeg:List NNI,
--JHD                     lvar:List OV) : Union("failed",P) ==
--JHD        x:OV:=lvar.first
--JHD        leadpol:Boolean:=false
--JHD        (listpol,lval):=(uterm.lpol,uterm.lint)
--JHD        d:=listpol.first
--JHD        listpol:=listpol.rest
--JHD        nolift:Boolean:=true
--JHD        for uf in listpol repeat
--JHD              --note uf and d not necessarily primitive
--JHD          degree gcd(uf,d) =0 => nolift:=false
--JHD        nolift => "failed"
--JHD        f:P:=([p1,p2]$List(P)).(position(uf,listpol))
--JHD        lgcd:=gcd(leadingCoefficient univariate(p1,x),
--JHD                  leadingCoefficient univariate(p2,x))
--JHD        lift(f,d,uf,lgcd,lvar,ldeg,lval)
--JHD
--JHD   -- interface with the general "lifting" function
--JHD      lift(f:P,d:SUPR,uf:SUPR,lgcd:P,lvar:List OV,
--JHD                        ldeg:List NNI,lval:List R):P ==
--JHD        x:OV:=lvar.first
--JHD        leadpol:Boolean:=false
--JHD        lcf:P
--JHD        lcf:=leadingCoefficient univariate(f,x)
--JHD        df:=degree(f,x)
--JHD        leadlist:List(P):=[]
--JHD
--JHD        if lgcd~=1$P then
--JHD          leadpol:=true
--JHD          f:=lgcd*f
--JHD          ldeg:=[n0+n1 for n0 in ldeg for n1 in degree(lgcd,lvar)]
--JHD          lcd:R:=leadingCoefficient d
--JHD          if ground? lgcd then d:=((retract lgcd) *d exquo lcd)::SUPR
--JHD          else d:=(retract(eval(lgcd,lvar.rest,lval)) * d exquo lcd)::SUPR
--JHD          uf:=lcd*uf
--JHD        leadlist:=[lgcd,lcf]
--JHD        lg:=imposelc([d,uf],lvar,lval,leadlist)
--JHD        plist:=lifting(univariate(f,x),lvar,lg,lval,leadlist,ldeg)::List P
--JHD        (p0:P,p1:P):=(plist.first,plist.2)
--JHD        if univariate eval(p0,rest lvar,lval) ~= lg.first then
--JHD           (p0,p1):=(p1,p0)
--JHD        ^leadpol => p0
--JHD        cprim:=contprim([p0])
--JHD        cprim.first.prim
--JHD
--JHD  -- Gcd for two multivariate polynomials
--JHD      gcd(p1:P,p2:P) : P ==
--JHD        (p1:= abs(p1)) = (p2:= abs(p2)) => p1
--JHD        ground? p1 =>
--JHD          p1 = 1$P => p1
--JHD          p1 = 0$P => p2
--JHD          ground? p2 => gcd((retract p1)@R,(retract p2)@R)::P
--JHD          gcdTermList(p1,p2)
--JHD        ground? p2 =>
--JHD          p2 = 1$P => p2
--JHD          p2 = 0$P => p1
--JHD          gcdTermList(p2,p1)
--JHD        mv1:= mainVariable(p1)::OV
--JHD        mv2:= mainVariable(p2)::OV
--JHD        mv1 = mv2 => gcdSameMainvar(p1,p2,mv1)
--JHD        mv1 < mv2 => gcdTermList(p1,p2)
--JHD        gcdTermList(p2,p1)
--JHD
--JHD  -- Gcd for a list of multivariate polynomials
--JHD      gcd(listp:List P) : P ==
--JHD        lf:=sort(better,listp)
--JHD        f:=lf.first
--JHD        for g in lf.rest repeat
--JHD          f:=gcd(f,g)
--JHD          if f=1$P then return f
--JHD        f
--JHD   -- Gcd and cofactors for a list of polynomials
--JHD      gcdcofact(listp : List P) : List P ==
--JHD         h:=gcd listp
--JHD         cons(h,[(f exquo h) :: P  for f in listp])
--JHD
--JHD   -- Gcd for primitive polynomials
--JHD      gcdprim(p1:P,p2:P):P ==
--JHD        (p1:= abs(p1)) = (p2:= abs(p2)) => p1
--JHD        ground? p1 =>
--JHD          ground? p2 => gcd((retract p1)@R,(retract p2)@R)::P
--JHD          p1 = 0$P => p2
--JHD          1$P
--JHD        ground? p2 =>
--JHD          p2 = 0$P => p1
--JHD          1$P
--JHD        mv1:= mainVariable(p1)::OV
--JHD        mv2:= mainVariable(p2)::OV
--JHD        mv1 = mv2 =>
--JHD          md:=min(minimumDegree(p1,mv1),minimumDegree(p2,mv1))
--JHD          mp:=1$P
--JHD          if md>1 then
--JHD            mp:=(mv1::P)**md
--JHD            p1:=(p1 exquo mp)::P
--JHD            p2:=(p2 exquo mp)::P
--JHD          mp*gcdPrim(p1,p2,mv1)
--JHD        1$P
--JHD
--JHD  -- Gcd for a list of primitive multivariate polynomials
--JHD      gcdprim(listp:List P) : P ==
--JHD        lf:=sort(better,listp)
--JHD        f:=lf.first
--JHD        for g in lf.rest repeat
--JHD          f:=gcdprim(f,g)
--JHD          if f=1$P then return f
--JHD        f
--JHD   -- Gcd and cofactors for a list of primitive polynomials
--JHD      gcdcofactprim(listp : List P) : List P ==
--JHD         h:=gcdprim listp
--JHD         cons(h,[(f exquo h) :: P  for f in listp])
--JHD
--JHD   -- content of a polynomial (with respect to its main var)
--JHD      content(f:P):P ==
--JHD        ground? f => f
--JHD        x:OV:=(mainVariable f)::OV
--JHD        gcd sort(better,coefficients univariate(f,x))
--JHD
--JHD   -- contents of a list of polynomials
--JHD      content(listf:List P) : List P == [content f for f in listf]
--JHD
--JHD   -- contents and primitive parts of a list  of polynomials
--JHD      contprim(listf:List P) : List ContPrim ==
--JHD        prelim :List P := content listf
--JHD        [[q,(f exquo q)::P]$ContPrim for q in prelim for f in listf]
--JHD

@
\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 GENPGCD GeneralPolynomialGcdPackage>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}