11
11
;;; ** (C) COPYRIGHT 1982 MASSACHUSETTS INSTITUTE OF TECHNOLOGY *****
12
12
;;; ****** THIS IS A READ-ONLY FILE! (ALL WRITES RESERVED) **********
13
13
;;; *****************************************************************
16
16
(macsyma-module spgcd)
18
18
(declare-top (special modulus temp genvar *alpha *which-factor*
19
$algebraic algfac* $GCD)
19
$algebraic algfac* $gcd)
23
23
(load-macsyma-macros ratmac)
25
25
(defvar *alpha nil)
27
(defmvar $POINTBOUND *alpha)
27
(defmvar $pointbound *alpha)
29
29
(defmacro 0? (x) `(= ,x 0))
31
31
(defmacro melt (a . indices) `(arraycall fixnum ,a . ,indices))
33
(defmacro ARRAYTYPE (m) `(array-type ,m))
35
(defmacro NCOLS (m) `(array-dimension-n ,m 1))
37
(defmacro NROWS (m) `(array-dimension-n ,m 2))
39
(defmacro LEN (lobj) `(cadr ,lobj))
41
(defmacro SKEL (lobj) `(caddr ,lobj))
43
(defmacro MATRIX (lobj) `(cadddr ,lobj))
45
(defmacro CURROW (lobj) `(cadr (cdddr ,lobj)))
47
(defmacro BADROWS (lobj) `(cadr (cddddr ,lobj)))
33
(defmacro arraytype (m) `(array-type ,m))
35
(defmacro ncols (m) `(array-dimension-n ,m 1))
37
(defmacro nrows (m) `(array-dimension-n ,m 2))
39
(defmacro len (lobj) `(cadr ,lobj))
41
(defmacro skel (lobj) `(caddr ,lobj))
43
(defmacro matrix (lobj) `(cadddr ,lobj))
45
(defmacro currow (lobj) `(cadr (cdddr ,lobj)))
47
(defmacro badrows (lobj) `(cadr (cddddr ,lobj)))
50
(defun PINTERP (x pts vals)
53
(pctimes (crecip (pcsubstz (car xk) qk))
56
(pcsubstz (car xk) u)))))
57
(uk (cdr vals) (cdr uk))
58
(qk (list x 1 1 0 (cminus (car pts)))
59
(ptimes qk (list x 1 1 0 (cminus (car xk)))))
60
(xk (cdr pts) (cdr xk)))
50
(defun pinterp (x pts vals)
53
(pctimes (crecip (pcsubstz (car xk) qk))
56
(pcsubstz (car xk) u)))))
57
(uk (cdr vals) (cdr uk))
58
(qk (list x 1 1 0 (cminus (car pts)))
59
(ptimes qk (list x 1 1 0 (cminus (car xk)))))
60
(xk (cdr pts) (cdr xk)))
63
(defun PCSUBSTZ (val p)
65
(do ((l (p-terms p) (pt-red l))
69
(cexpt val (f- (pt-le l) (pt-le (pt-red l)))))))
63
(defun pcsubstz (val p)
65
(do ((l (p-terms p) (pt-red l))
72
68
(cplus ans (pt-lc l))
75
(cexpt val (pt-le l))))))))
77
;skeletons consist of list of exponent vectors
78
;lifting structures are as follows:
80
;(matrix <n> ;length of <skel>
81
; <skel> ;skeleton of poly being lifted
82
; <matrix> ; partially diagonalized matrix used to obtain sol.
83
; <row> ;we have diagonalized this far
84
; <bad rows>) ; after we get
86
(defun EVAL-MON (mon pt)
91
(cexpt (car ll) (car l)))
69
(cexpt val (f- (pt-le l) (pt-le (pt-red l)))))))
75
(cexpt val (pt-le l))))))))
77
;;skeletons consist of list of exponent vectors
78
;;lifting structures are as follows:
80
;;(matrix <n> ;length of <skel>
81
;; <skel> ;skeleton of poly being lifted
82
;; <matrix> ; partially diagonalized matrix used to obtain sol.
83
;; <row> ;we have diagonalized this far
84
;; <bad rows>) ; after we get
86
(defun eval-mon (mon pt)
91
(cexpt (car ll) (car l)))
95
; ONE-STEP causes the (row,col) element of mat to be made to be zero. It is
96
; assumed that row > col, and that the matrix is diagonal above the row.
97
; n indicates how wide the rows are suppoded to be.
95
;; ONE-STEP causes the (row,col) element of mat to be made to be zero. It is
96
;; assumed that row > col, and that the matrix is diagonal above the row.
97
;; n indicates how wide the rows are suppoded to be.
99
(defun ONE-STEP (mat row col n)
101
(c (arraycall fixnum mat row col))) ;Got to save this away before it is
99
(defun one-step (mat row col n)
101
(c (arraycall fixnum mat row col))) ;Got to save this away before it is
104
(store (arraycall fixnum mat row i)
105
(cdifference (arraycall fixnum mat row i)
106
(ctimes c (arraycall fixnum mat col i))))))
108
; MONICIZE-ROW assumes the (row,row) element is non-zero and that this element
109
; is to be made equal to one. It merely divides the entire row by the
110
; (row,row) element. It is assumed that the (row,n) elements with n < row are
111
; already zero. Again n is the width of the rows.
113
(defun MONICIZE-ROW (mat row n)
115
(c (CRECIP (arraycall fixnum mat row row))))
117
(store (arraycall fixnum mat row row) 1)) ;Don't bother doing the
104
(store (arraycall fixnum mat row i)
105
(cdifference (arraycall fixnum mat row i)
106
(ctimes c (arraycall fixnum mat col i))))))
108
;; MONICIZE-ROW assumes the (row,row) element is non-zero and that this element
109
;; is to be made equal to one. It merely divides the entire row by the
110
;; (row,row) element. It is assumed that the (row,n) elements with n < row are
111
;; already zero. Again n is the width of the rows.
113
(defun monicize-row (mat row n)
115
(c (crecip (arraycall fixnum mat row row))))
117
(store (arraycall fixnum mat row row) 1)) ;Don't bother doing the
118
118
;division on the diagonal
119
(store (arraycall fixnum mat row i) (ctimes (arraycall fixnum mat row i) c))))
121
; FILL-ROW is given a point and the value of the skeleton at the point and
122
; fills the apropriate row in the matrix. with the values of the monomials
123
; n is the length of the skeleton
125
(defun FILL-ROW (skel mat n row pt val)
128
((= i n) ;l should be nil now,
129
(if (not (null l)) ;but lets check just to make sure
130
(merror "Skeleton too long in FILL-ROW: ~S"
132
(store (arraycall fixnum mat row n) val)) ;The value is stored in the
119
(store (arraycall fixnum mat row i) (ctimes (arraycall fixnum mat row i) c))))
121
;; FILL-ROW is given a point and the value of the skeleton at the point and
122
;; fills the apropriate row in the matrix. with the values of the monomials
123
;; n is the length of the skeleton
125
(defun fill-row (skel mat n row pt val)
128
((= i n) ;l should be nil now,
129
(if (not (null l)) ;but lets check just to make sure
130
(merror "Skeleton too long in `fill-row': ~S"
132
(store (arraycall fixnum mat row n) val)) ;The value is stored in the
134
(store (arraycall fixnum mat row i)
135
(eval-mon (car l) pt)))) ;Should think about a better
134
(store (arraycall fixnum mat row i)
135
(eval-mon (car l) pt)))) ;Should think about a better
136
136
;way to do this evaluation.
139
(defun SWAP-ROWS (mat m n) ;Interchange row m and n
143
(store (arraycall fixnum mat m k)
144
(prog2 0 (melt mat n k)
145
(store (melt mat n k) (melt mat m k))))))
139
(defun swap-rows (mat m n) ;Interchange row m and n
143
(store (arraycall fixnum mat m k)
144
(prog2 0 (melt mat n k)
145
(store (melt mat n k) (melt mat m k))))))
147
; PARTIAL-DIAG fills in one more row of the matrix and tries to diagonalize
148
; what it has so far. If the row which has just been introduced has a
149
; zero element on the diagonal then it is scanned for its first non-zero
150
; element and it is placed in the matrix so that the non-zero element is
151
; on the diagonal. The rows which are missing are added to the list in
152
; badrows. The current row pointer may skip a couple of numbers here, so
153
; when it is equal to n, the only empty rows to add thing to are on the
147
;; PARTIAL-DIAG fills in one more row of the matrix and tries to diagonalize
148
;; what it has so far. If the row which has just been introduced has a
149
;; zero element on the diagonal then it is scanned for its first non-zero
150
;; element and it is placed in the matrix so that the non-zero element is
151
;; on the diagonal. The rows which are missing are added to the list in
152
;; badrows. The current row pointer may skip a couple of numbers here, so
153
;; when it is equal to n, the only empty rows to add thing to are on the
156
(defun PARTIAL-DIAG (lobj pt val) ; Does one step of diagonalization of
157
(let ((n (len lobj)) ;the matrix in lobj
159
(mat (matrix lobj)) ;The matrix, obviously
160
(row (currow lobj)) ;This is the row which is to be
156
(defun partial-diag (lobj pt val) ; Does one step of diagonalization of
157
(let ((n (len lobj)) ;the matrix in lobj
159
(mat (matrix lobj)) ;The matrix, obviously
160
(row (currow lobj)) ;This is the row which is to be
162
(badrows (badrows lobj)) ;Rows which could not be used when
162
(badrows (badrows lobj)) ;Rows which could not be used when
163
163
;their time came, and will be used later
165
(cond ((= row n) ;All rows already done, must start
165
(cond ((= row n) ;All rows already done, must start
166
166
;using the rows in badrows.
167
(fill-row skel mat n (setq crow (car badrows)) pt val))
168
((fill-row skel mat n row pt val) ;Fill up the data
171
;; This loop kills all the elements in the row up to the diagonal.
173
(do ((i 0 (f1+ i)) ;runs over the columns of mat
174
(l (setq badrows (cons nil badrows)))
177
(setq badrows (cdr badrows))) ;update badrows
178
(if (cdr l) ;Any bad rows around?
179
(if (= (cadr l) i) ;No one for this row,
180
(if (and (null flag) ;So if this guy can fill in
181
(not (0? (melt mat crow i))))
183
(swap-rows mat crow i) ;Put this guy in the right spot
185
(setq flag t crow i)) ; and make him a winner
186
(setq l (cdr l))) ;At any rate this guy isn't
167
(fill-row skel mat n (setq crow (car badrows)) pt val))
168
((fill-row skel mat n row pt val) ;Fill up the data
171
;; This loop kills all the elements in the row up to the diagonal.
173
(do ((i 0 (f1+ i)) ;runs over the columns of mat
174
(l (setq badrows (cons nil badrows)))
177
(setq badrows (cdr badrows))) ;update badrows
178
(if (cdr l) ;Any bad rows around?
179
(if (= (cadr l) i) ;No one for this row,
180
(if (and (null flag) ;So if this guy can fill in
181
(not (0? (melt mat crow i))))
183
(swap-rows mat crow i) ;Put this guy in the right spot
185
(setq flag t crow i)) ; and make him a winner
186
(setq l (cdr l))) ;At any rate this guy isn't
188
(one-step mat crow i n))
189
(one-step mat crow i n)))
188
(one-step mat crow i n))
189
(one-step mat crow i n)))
191
(if (0? (melt mat crow crow)) ;diagonal is zero?
192
(setq badrows (cons crow badrows))
194
(monicize-row mat crow n) ;Monicize the diagonal element on this
191
(if (0? (melt mat crow crow)) ;diagonal is zero?
192
(setq badrows (cons crow badrows))
194
(monicize-row mat crow n) ;Monicize the diagonal element on this
196
(do ((j 0 (f1+ j))) ;For each element in the rows above
196
(do ((j 0 (f1+ j))) ;For each element in the rows above
197
197
;this one zero the the crow column
198
((= j crow)) ;Don't zero the diagonal element
199
(one-step mat j crow n))))
200
(cond ((and (= (f1+ row) n)
201
(progn (setq row (f1- row)) t) ;Decrement it just in case
202
(null (cdr badrows)))
203
(do ((l nil (cons (melt mat i n) l))
206
(list 'SOLUTION n skel mat l))))
207
(t (list 'MATRIX n skel mat (f1+ row) badrows)))))
198
((= j crow)) ;Don't zero the diagonal element
199
(one-step mat j crow n))))
200
(cond ((and (= (f1+ row) n)
201
(progn (setq row (f1- row)) t) ;Decrement it just in case
202
(null (cdr badrows)))
203
(do ((l nil (cons (melt mat i n) l))
206
(list 'solution n skel mat l))))
207
(t (list 'matrix n skel mat (f1+ row) badrows)))))
209
(defun GEN-POINT (vars)
210
(do ((l vars (cdr l))
211
(ans nil (cons (cmod (random $POINTBOUND)) ans)))
209
(defun gen-point (vars)
210
(do ((l vars (cdr l))
211
(ans nil (cons (cmod (random $pointbound)) ans)))
214
; PDIAG-ALL introduces a new row in each matrix in the list of l-lobjs.
215
; The RHS of the equations is stripped off of poly.
214
;; PDIAG-ALL introduces a new row in each matrix in the list of l-lobjs.
215
;; The RHS of the equations is stripped off of poly.
217
(defun PDIAG-ALL (l-lobjs poly pt)
218
(do ((l l-lobjs (cdr l))
222
(if solved (cons 'SOLVED l-lobjs)
224
(if (and lp (= (caar l) (pt-le lp))) ;This term corresponds to the
217
(defun pdiag-all (l-lobjs poly pt)
218
(do ((l l-lobjs (cdr l))
222
(if solved (cons 'solved l-lobjs)
224
(if (and lp (= (caar l) (pt-le lp))) ;This term corresponds to the
225
225
;current lobj, set c to the coefficient
226
(setq c (pt-lc lp) lp (pt-red lp))
228
;FIXTHIS ;Should put it some check for extra
226
(setq c (pt-lc lp) lp (pt-red lp))
228
;;FIXTHIS ;Should put it some check for extra
229
229
;terms in the polynomial
230
(if (not (eq (cadar l) 'SOLUTION))
231
(progn (rplacd (car l)
232
(partial-diag (cdar l) pt c))
233
(and solved (null (eq (cadar l) 'SOLUTION))
234
(setq solved nil))))))
230
(if (not (eq (cadar l) 'solution))
231
(progn (rplacd (car l)
232
(partial-diag (cdar l) pt c))
233
(and solved (null (eq (cadar l) 'solution))
234
(setq solved nil))))))
236
236
;; not currently called
237
237
;; (defun CREATE-INTVECT (h)
253
253
;; (t (error '|Bad Evaluation point - MERGE-INTVECT|)))))
256
(defun MERGE-SKEL (mon poly)
259
((do ((l (cdr poly) (cddr l))
261
(cons (cons (car l) mon) ans)))
264
(defun NEW-SKEL (skel polys)
266
(mapcan (fn (mon poly) (merge-skel mon poly))
269
(cond ((pcoefp q) (list q))
270
((do ((l (cdr q) (cddr l))
271
(ans nil (cons (cadr l) ans)))
275
(defun CREATE-LOBJS (prev-lift)
256
(defun merge-skel (mon poly)
259
((do ((l (cdr poly) (cddr l))
261
(cons (cons (car l) mon) ans)))
264
(defun new-skel (skel polys)
266
(mapcan (fn (mon poly) (merge-skel mon poly))
269
(cond ((pcoefp q) (list q))
270
((do ((l (cdr q) (cddr l))
271
(ans nil (cons (cadr l) ans)))
275
(defun create-lobjs (prev-lift)
276
276
(mapcar #'(lambda (q)
277
277
(let ((n (length (cadr q))))
279
(list 'MATRIX n (cadr q)
279
(list 'matrix n (cadr q)
280
280
(*array nil 'fixnum n (f1+ n))
284
(defun CLEAR-LOBJS (lobjs)
284
(defun clear-lobjs (lobjs)
285
285
(mapcar #'(lambda (q)
287
(list 'MATRIX (caddr q) (cadddr q)
287
(list 'matrix (caddr q) (cadddr q)
288
288
(caddr (cddr q)) 0 nil)))
291
(defun SPARSE-LIFT (c f g l-lobjs vars)
292
(do ((pt (gen-point vars) (gen-point vars))
294
((eq (car l-lobjs) 'SOLVED)
296
(setq gcd (lifting-factors-image
297
(pcsub c pt vars) (pcsub f pt vars) (pcsub g pt vars)))
299
(not (= (pt-le (p-terms gcd)) (caar l-lobjs))))
300
(throw 'Bad-Point nil)
301
(setq l-lobjs (pdiag-all l-lobjs gcd pt)))))
291
(defun sparse-lift (c f g l-lobjs vars)
292
(do ((pt (gen-point vars) (gen-point vars))
294
((eq (car l-lobjs) 'solved)
296
(setq gcd (lifting-factors-image
297
(pcsub c pt vars) (pcsub f pt vars) (pcsub g pt vars)))
299
(not (= (pt-le (p-terms gcd)) (caar l-lobjs))))
300
(throw 'bad-point nil)
301
(setq l-lobjs (pdiag-all l-lobjs gcd pt)))))
303
(defun LIFTING-FACTORS-IMAGE (c f g)
304
(let ((gcd (pgcdu f g)))
307
(2 (pquotient f gcd))
308
(3 (pquotient g gcd)))))
303
(defun lifting-factors-image (c f g)
304
(let ((gcd (pgcdu f g)))
307
(2 (pquotient f gcd))
308
(3 (pquotient g gcd)))))
310
(defun ZGCD-LIFT* (c f g vars degb)
311
(do ((vals (gen-point vars) (gen-point vars))
317
(ZGCD-LIFT c f g vars vals degb)))))
319
; ZGCD-LIFT returns objects called lifts. These have the the following
321
; ((n <skel> <poly>) ... )
322
; n corresponds to the degree in the main variable to which this guy
325
(defun ZGCD-LIFT (c f g vars vals degb)
326
(cond ((null vars) ;No variables left, just the main one
327
(let ((p (lifting-factors-image c f g))) ;Compute factor and
310
(defun zgcd-lift* (c f g vars degb)
311
(do ((vals (gen-point vars) (gen-point vars))
317
(zgcd-lift c f g vars vals degb)))))
319
;; ZGCD-LIFT returns objects called lifts. These have the the following
321
;; ((n <skel> <poly>) ... )
322
;; n corresponds to the degree in the main variable to which this guy
325
(defun zgcd-lift (c f g vars vals degb)
326
(cond ((null vars) ;No variables left, just the main one
327
(let ((p (lifting-factors-image c f g))) ;Compute factor and
328
328
;enforce leading coefficient
329
(if (pcoefp p) (throw 'relprime 1) ;if the GCD is 1 quit
330
(do ((l (p-terms p) (pt-red l)) ;otherwise march
329
(if (pcoefp p) (throw 'relprime 1) ;if the GCD is 1 quit
330
(do ((l (p-terms p) (pt-red l)) ;otherwise march
331
331
;though the polynomial
332
(ans nil ;constructing a lift for each term.
333
(cons (list (pt-le l) '(nil) (list (pt-lc l)))
337
((let ((prev-lift ;Recurse if possible
338
(zgcd-lift (pcsubsty (car vals) (car vars) c)
339
(pcsubsty (car vals) (car vars) f)
340
(pcsubsty (car vals) (car vars) g)
341
(cdr vars) (cdr vals) (cdr degb))))
342
(do ((i 0 (f1+ i)) ;counts to the degree bound
343
(lobjs (create-lobjs prev-lift) ;need to create
332
(ans nil ;constructing a lift for each term.
333
(cons (list (pt-le l) '(nil) (list (pt-lc l)))
337
((let ((prev-lift ;Recurse if possible
338
(zgcd-lift (pcsubsty (car vals) (car vars) c)
339
(pcsubsty (car vals) (car vars) f)
340
(pcsubsty (car vals) (car vars) g)
341
(cdr vars) (cdr vals) (cdr degb))))
342
(do ((i 0 (f1+ i)) ;counts to the degree bound
343
(lobjs (create-lobjs prev-lift) ;need to create
344
344
;the appropriate matrices
345
(clear-lobjs lobjs)) ;but reuse them at each
345
(clear-lobjs lobjs)) ;but reuse them at each
347
(pts (add-point (list (car vals))) ;List of random
348
(add-point pts)) ;points
349
(linsols (mapcar 'make-linsols prev-lift)
353
(pcsubsty (car pts) (car vars) c)
354
(pcsubsty (car pts) (car vars) f)
355
(pcsubsty (car pts) (car vars) g)
358
(interp-polys linsols (cdr pts) (car vars))))))))
347
(pts (add-point (list (car vals))) ;List of random
348
(add-point pts)) ;points
349
(linsols (mapcar 'make-linsols prev-lift)
353
(pcsubsty (car pts) (car vars) c)
354
(pcsubsty (car pts) (car vars) f)
355
(pcsubsty (car pts) (car vars) g)
358
(interp-polys linsols (cdr pts) (car vars))))))))
360
(defun MAKE-LINSOLS (prev-lift)
360
(defun make-linsols (prev-lift)
361
361
(list (car prev-lift)
363
363
(mapcan #'(lambda (q)
393
(defun ZGCD (f g &aux $algebraic algfac*)
394
(let ((f (oldcontent f)) ;This is a good spot to
395
(g (oldcontent g)) ;initialize random
398
;; *WHICH-FACTOR* is interpreted as follows. It is set fairly deep in the
399
;; algorithm, inside ZGCD1.
401
;; 2 -> Lift the cofactor of F
402
;; 3 -> Lift the cofactor of G
405
(setq mon (pgcd (car f) (car g)) ;compute GCD of content
406
f (cadr f) g (cadr g)) ;f and g are now primitive
407
(if (or (pcoefp f) (pcoefp g)
408
(not (eq (car f) (car g))))
409
(merror "Bad args to ZGCD"))
413
(setq gcd (catch 'relprime (zgcd1 f g)))
416
(1 (testdivide f gcd))
417
(2 (testdivide f gcd))
418
(3 (testdivide g gcd))))
419
(cond ((not (null test))
421
(cond ((equal *which-factor* 1)
424
((not (null modulus))
425
(return (let (($GCD '$RED))
429
(let* ((modulus modulus)
432
(vars (sort (union1 (listovars f) (listovars g))
434
(GENVAR (REVERSE VARS))
436
;; the elements of the following degree vectors all correspond to the
437
;; contents of var. Thus there may be variables missing that are in
439
;; (f-degv (zpdegreevector f vars)) ;;WHY NOT JUST PUT GENVAR THE REVERSE
440
;; (g-degv (zpdegreevector g vars)) ;;THE REVERSE OF VARS AND JUST USE PDEGREEVECTOR--wfs
441
(F-DEGV (PDEGREEVECTOR F))
442
(G-DEGV (PDEGREEVECTOR G))
443
(gcd-degv (gcd-degree-vector f g vars)))
445
;; First we try to decide which of the gcd and the cofactors of f and g
446
;; is smallest. The result of this decision is indicated by *which-factor*.
447
;; Then the leading coefficient that is to be enforced is changed if a
448
;; cofactor has been chosen.
449
(case (setq *which-factor*
450
(determine-lifting-factor f-degv g-degv gcd-degv))
451
(1 (setq c (pgcd (p-lc f) (p-lc g))))
452
(2 (setq c (p-lc f)))
453
(3 (setq c (p-lc g))))
455
;; Next a degree bound must be chosen.
461
(2 (mapcar #'- f-degv gcd-degv))
462
(3 (mapcar #'- g-degv gcd-degv)))
463
(zpdegreevector c vars))))
465
(cond ((not (null modulus))
466
(lobj->poly (car vars) (cdr vars)
467
(zgcd-lift* c f g (cdr vars) (cdr degb))))
469
(setq h (times (maxcoefficient f)
471
(setq modulus *alpha) ;Really should randomize
473
(zgcd-lift* (pmod c) (pmod f) (pmod g)
474
(cdr vars) (cdr degb)))
475
(do ((linsols (mapcar #'(lambda (q)
477
(new-skel (cadr q) (caddr q))))
479
(merge-sol-lin-z linsols
480
(sparse-lift cm fm gm lobjs (cdr vars))
482
(crecip (cmod coef-bound)))
483
(times modulus coef-bound)))
484
(lobjs (create-lobjs first-lift)
486
(coef-bound *alpha (times modulus coef-bound))
488
((greaterp coef-bound H)
490
(lobj->poly (car vars) (cdr vars) linsols))
491
(setq modulus (newprime modulus))
496
(defun LOBJ->POLY (var vars lobj)
502
(do ((x (cadr q) (cdr x))
503
(y (caddr q) (cdr y))
506
(disrep-monom (cdar x) (car y)
511
(defun DISREP-MONOM (monom c vars)
512
(cond ((null monom) c)
513
((equal 0 (car monom))
514
(disrep-monom (cdr monom) c (cdr vars)))
515
((list (car vars) (car monom)
516
(disrep-monom (cdr monom) c (cdr vars))))))
518
(defun MERGE-SOL-LIN-Z (l1 l2 c new-coef-bound)
521
(modulus new-coef-bound)
524
(cond ((= (caar l) (caar l2))
393
(defun zgcd (f g &aux $algebraic algfac*)
394
(let ((f (oldcontent f)) ;This is a good spot to
395
(g (oldcontent g)) ;initialize random
398
;; *WHICH-FACTOR* is interpreted as follows. It is set fairly deep in the
399
;; algorithm, inside ZGCD1.
401
;; 2 -> Lift the cofactor of F
402
;; 3 -> Lift the cofactor of G
405
(setq mon (pgcd (car f) (car g)) ;compute GCD of content
406
f (cadr f) g (cadr g)) ;f and g are now primitive
407
(if (or (pcoefp f) (pcoefp g)
408
(not (eq (car f) (car g))))
409
(merror "Bad args to `zgcd'"))
413
(setq gcd (catch 'relprime (zgcd1 f g)))
416
(1 (testdivide f gcd))
417
(2 (testdivide f gcd))
418
(3 (testdivide g gcd))))
419
(cond ((not (null test))
421
(cond ((equal *which-factor* 1)
424
((not (null modulus))
425
(return (let (($gcd '$red))
429
(let* ((modulus modulus)
432
(vars (sort (union1 (listovars f) (listovars g))
434
(genvar (reverse vars))
436
;; the elements of the following degree vectors all correspond to the
437
;; contents of var. Thus there may be variables missing that are in
439
;; (f-degv (zpdegreevector f vars)) ;;WHY NOT JUST PUT GENVAR THE REVERSE
440
;; (g-degv (zpdegreevector g vars)) ;;THE REVERSE OF VARS AND JUST USE PDEGREEVECTOR--wfs
441
(f-degv (pdegreevector f))
442
(g-degv (pdegreevector g))
443
(gcd-degv (gcd-degree-vector f g vars)))
445
;; First we try to decide which of the gcd and the cofactors of f and g
446
;; is smallest. The result of this decision is indicated by *which-factor*.
447
;; Then the leading coefficient that is to be enforced is changed if a
448
;; cofactor has been chosen.
449
(case (setq *which-factor*
450
(determine-lifting-factor f-degv g-degv gcd-degv))
451
(1 (setq c (pgcd (p-lc f) (p-lc g))))
452
(2 (setq c (p-lc f)))
453
(3 (setq c (p-lc g))))
455
;; Next a degree bound must be chosen.
461
(2 (mapcar #'- f-degv gcd-degv))
462
(3 (mapcar #'- g-degv gcd-degv)))
463
(zpdegreevector c vars))))
465
(cond ((not (null modulus))
466
(lobj->poly (car vars) (cdr vars)
467
(zgcd-lift* c f g (cdr vars) (cdr degb))))
469
(setq h (times (maxcoefficient f)
471
(setq modulus *alpha) ;Really should randomize
473
(zgcd-lift* (pmod c) (pmod f) (pmod g)
474
(cdr vars) (cdr degb)))
475
(do ((linsols (mapcar #'(lambda (q)
477
(new-skel (cadr q) (caddr q))))
479
(merge-sol-lin-z linsols
480
(sparse-lift cm fm gm lobjs (cdr vars))
482
(crecip (cmod coef-bound)))
483
(times modulus coef-bound)))
484
(lobjs (create-lobjs first-lift)
486
(coef-bound *alpha (times modulus coef-bound))
488
((greaterp coef-bound h)
490
(lobj->poly (car vars) (cdr vars) linsols))
491
(setq modulus (newprime modulus))
496
(defun lobj->poly (var vars lobj)
502
(do ((x (cadr q) (cdr x))
503
(y (caddr q) (cdr y))
506
(disrep-monom (cdar x) (car y)
511
(defun disrep-monom (monom c vars)
512
(cond ((null monom) c)
513
((equal 0 (car monom))
514
(disrep-monom (cdr monom) c (cdr vars)))
515
((list (car vars) (car monom)
516
(disrep-monom (cdr monom) c (cdr vars))))))
518
(defun merge-sol-lin-z (l1 l2 c new-coef-bound)
521
(modulus new-coef-bound)
524
(cond ((= (caar l) (caar l2))
531
(cplus b (ctimes c (cdifference a b)))))
535
(cadddr (cddar l2)) (caddar l)))))))
531
(cplus b (ctimes c (cdifference a b)))))
535
(cadddr (cddar l2)) (caddar l)))))))
537
537
;; The following function tries to determine the degree of gcd(f, g) in each
538
538
;; variable. This is done in the following manner: All but one of the
543
543
;; The univariate gcd's are computed with modulus=*alpha.
545
(defun GCD-DEGREE-VECTOR (f g vars)
546
(let ((modulus *alpha))
547
(setq f (pmod f) g (pmod g))
548
(do ((vn (cdr vars) (cdr vn))
549
(vs (zl-DELETE (car vars) (copy1 vars))
550
(zl-DELETE (car vn) (copy1 vars)))
551
(l) (f*) (g*) (gcd*) (rand))
553
(setq rand (gen-point vs))
554
(setq f* (pcsub f rand vs)
555
g* (pcsub g rand vs))
556
(cond ((or (pcoefp f*) (pcoefp g*)
557
(pcoefp (setq gcd* (pgcdu f* g*))))
559
(t (push (pt-le (p-terms gcd*)) l)))
561
(return l)))))) ;No reverse needed here
545
(defun gcd-degree-vector (f g vars)
546
(let ((modulus *alpha))
547
(setq f (pmod f) g (pmod g))
548
(do ((vn (cdr vars) (cdr vn))
549
(vs (zl-delete (car vars) (copy1 vars))
550
(zl-delete (car vn) (copy1 vars)))
551
(l) (f*) (g*) (gcd*) (rand))
553
(setq rand (gen-point vs))
554
(setq f* (pcsub f rand vs)
555
g* (pcsub g rand vs))
556
(cond ((or (pcoefp f*) (pcoefp g*)
557
(pcoefp (setq gcd* (pgcdu f* g*))))
559
(t (push (pt-le (p-terms gcd*)) l)))
561
(return l)))))) ;No reverse needed here
563
; DETERMINE-LIFTING-FACTOR returns a number indicating which factor of f or g
563
;; DETERMINE-LIFTING-FACTOR returns a number indicating which factor of f or g
566
566
(defun dlf-mumblify (a b)
567
(sloop for x in a for y in b sum (difference x y)))
567
(loop for x in a for y in b sum (difference x y)))
569
(defun DETERMINE-LIFTING-FACTOR (f-degv g-degv gcd-degv)
570
(let* ((fv ;(apply '+ (mapcar '- f-degv gcd-degv))
571
(dlf-mumblify f-degv gcd-degv))
572
(gv ;(apply '+ (mapcar '- g-degv gcd-degv))
573
(dlf-mumblify g-degv gcd-degv))
569
(defun determine-lifting-factor (f-degv g-degv gcd-degv)
570
(let* ((fv ;(apply '+ (mapcar '- f-degv gcd-degv))
571
(dlf-mumblify f-degv gcd-degv))
572
(gv ;(apply '+ (mapcar '- g-degv gcd-degv))
573
(dlf-mumblify g-degv gcd-degv))
574
574
(gcdv (apply '+ gcd-degv)))
575
575
(if (lessp fv gcdv)
576
576
(if (lessp fv gv) 2 3)
577
577
(if (lessp gv gcdv) 3 1))))
578
;(defun DETERMINE-LIFTING-FACTOR (f-degv g-degv gcd-degv)
579
; (let* ((fv (apply '+ (mapcar '- f-degv gcd-degv)))
580
; (gv (apply '+ (mapcar '- g-degv gcd-degv)))
581
; (gcdv (apply '+ gcd-degv)))
582
; (if (lessp fv gcdv)
583
; (if (lessp fv gv) 2 3)
584
; (if (lessp gv gcdv) 3 1))))
578
;;(defun DETERMINE-LIFTING-FACTOR (f-degv g-degv gcd-degv)
579
;; (let* ((fv (apply '+ (mapcar '- f-degv gcd-degv)))
580
;; (gv (apply '+ (mapcar '- g-degv gcd-degv)))
581
;; (gcdv (apply '+ gcd-degv)))
582
;; (if (lessp fv gcdv)
583
;; (if (lessp fv gv) 2 3)
584
;; (if (lessp gv gcdv) 3 1))))
587
(defun EXCISE-EXTRA-VARIABLES (degv vars)
588
(do ((l (reverse degv) (cdr l))
589
(lv (reverse genvar) (cdr lv))
593
(cond ((eq (car lv) (car vars))
595
(setq vars (cdr vars))))))
587
(defun excise-extra-variables (degv vars)
588
(do ((l (reverse degv) (cdr l))
589
(lv (reverse genvar) (cdr lv))
593
(cond ((eq (car lv) (car vars))
595
(setq vars (cdr vars))))))
597
(defun ZPDEGREEVECTOR (p vars)
598
(excise-extra-variables (pdegreevector p) vars))
597
(defun zpdegreevector (p vars)
598
(excise-extra-variables (pdegreevector p) vars))
added
removed
src/spgcd.lisp
