close Warning: Can't synchronize with repository "(default)" (The repository directory has changed, you should resynchronize the repository with: trac-admin $ENV repository resync '(default)'). Look in the Trac log for more information.

Changeset 54


Ignore:
Timestamp:
2015-06-05T11:16:40-07:00 (10 years ago)
Author:
Marek Rychlik
Message:

* empty log message *

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/f4grobner/grobner.lisp

    r47 r54  
    3333($load "functs")
    3434
    35 ;; Macros for making lists with iterators - an exammple of GENSYM
    36 ;; MAKELIST-1 makes a list with one iterator, while MAKELIST accepts an
    37 ;; arbitrary number of iterators
    38 
    39 ;; Sample usage:
    40 ;; Without a step:
    41 ;; >(makelist-1 (* 2 i) i 0 10)
    42 ;; (0 2 4 6 8 10 12 14 16 18 20)
    43 ;; With a step of 3:
    44 ;; >(makelist-1 (* 2 i) i 0 10 3)
    45 ;; (0 6 12 18)
    46 
    47 ;; Generate sums of squares of numbers between 1 and 4:
    48 ;; >(makelist (+ (* i i) (* j j)) (i 1 4) (j 1 i))
    49 ;; (2 5 8 10 13 18 17 20 25 32)
    50 ;; >(makelist (list i j '---> (+ (* i i) (* j j))) (i 1 4) (j 1 i))
    51 ;; ((1 1 ---> 2) (2 1 ---> 5) (2 2 ---> 8) (3 1 ---> 10) (3 2 ---> 13)
    52 ;; (3 3 ---> 18) (4 1 ---> 17) (4 2 ---> 20) (4 3 ---> 25) (4 4 ---> 32))
    53 
    54 ;; Evaluate expression expr with variable set to lo, lo+1,... ,hi
    55 ;; and put the results in a list.
    56 (defmacro makelist-1 (expr var lo hi &optional (step 1))
    57   (let ((l (gensym)))
    58     `(do ((,var ,lo (+ ,var ,step))
    59           (,l nil (cons ,expr ,l)))
    60          ((> ,var ,hi) (reverse ,l))
    61        (declare (fixnum ,var)))))
    62 
    63 (defmacro makelist (expr (var lo hi &optional (step 1)) &rest more)
    64   (if (endp more)
    65       `(makelist-1 ,expr ,var ,lo ,hi ,step)
    66     (let* ((l (gensym)))
    67       `(do ((,var ,lo (+ ,var ,step))
    68             (,l nil (nconc ,l `,(makelist ,expr ,@more))))
    69            ((> ,var ,hi) ,l)
    70          (declare (fixnum ,var))))))
    71 
    72 ;;----------------------------------------------------------------
    73 ;; This package implements BASIC OPERATIONS ON MONOMIALS
    74 ;;----------------------------------------------------------------
    75 ;; DATA STRUCTURES: Monomials are represented as lists:
    76 ;;
    77 ;;      monom:  (n1 n2 ... nk) where ni are non-negative integers
    78 ;;
    79 ;; However, lists may be implemented as other sequence types,
    80 ;; so the flexibility to change the representation should be
    81 ;; maintained in the code to use general operations on sequences
    82 ;; whenever possible. The optimization for the actual representation
    83 ;; should be left to declarations and the compiler.
    84 ;;----------------------------------------------------------------
    85 ;; EXAMPLES: Suppose that variables are x and y. Then
    86 ;;
    87 ;;      Monom x*y^2 ---> (1 2)
    88 ;;
    89 ;;----------------------------------------------------------------
    90 
    91 (deftype exponent ()
    92   "Type of exponent in a monomial."
    93   'fixnum)
    94 
    95 (deftype monom (&optional dim)
    96   "Type of monomial."
    97   `(simple-array exponent (,dim)))
    98 
    99 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    100 ;;
    101 ;; Construction of monomials
    102 ;;
    103 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    104 
    105 (defmacro make-monom (dim &key (initial-contents nil initial-contents-supplied-p)
    106                                (initial-element 0 initial-element-supplied-p))
    107   "Make a monomial with DIM variables. Additional argument
    108 INITIAL-CONTENTS specifies the list of powers of the consecutive
    109 variables. The alternative additional argument INITIAL-ELEMENT
    110 specifies the common power for all variables."
    111   ;;(declare (fixnum dim))
    112   `(make-array ,dim
    113                :element-type 'exponent
    114                ,@(when initial-contents-supplied-p `(:initial-contents ,initial-contents))
    115                ,@(when initial-element-supplied-p `(:initial-element ,initial-element))))
    116 
    117 
    118 
    119 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    120 ;;
    121 ;; Operations on monomials
    122 ;;
    123 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    124 
    125 (defmacro monom-elt (m index)
    126   "Return the power in the monomial M of variable number INDEX."
    127   `(elt ,m ,index))
    128 
    129 (defun monom-dimension (m)
    130   "Return the number of variables in the monomial M."
    131   (length m))
    132 
    133 (defun monom-total-degree (m &optional (start 0) (end (length m)))
    134   "Return the todal degree of a monomoal M. Optinally, a range
    135 of variables may be specified with arguments START and END."
    136   (declare (type monom m) (fixnum start end))
    137   (reduce #'+ m :start start :end end))
    138 
    139 (defun monom-sugar (m &aux (start 0) (end (length m)))
    140   "Return the sugar of a monomial M. Optinally, a range
    141 of variables may be specified with arguments START and END."
    142   (declare (type monom m) (fixnum start end))
    143   (monom-total-degree m start end))
    144 
    145 (defun monom-div (m1 m2 &aux (result (copy-seq m1)))
    146   "Divide monomial M1 by monomial M2."
    147   (declare (type monom m1 m2 result))
    148   (map-into result #'- m1 m2))
    149 
    150 (defun monom-mul (m1 m2  &aux (result (copy-seq m1)))
    151   "Multiply monomial M1 by monomial M2."
    152   (declare (type monom m1 m2 result))
    153   (map-into result #'+ m1 m2))
    154 
    155 (defun monom-divides-p (m1 m2)
    156   "Returns T if monomial M1 divides monomial M2, NIL otherwise."
    157   (declare (type monom m1 m2))
    158   (every #'<= m1 m2))
    159 
    160 (defun monom-divides-monom-lcm-p (m1 m2 m3)
    161   "Returns T if monomial M1 divides MONOM-LCM(M2,M3), NIL otherwise."
    162   (declare (type monom m1 m2 m3))
    163   (every #'(lambda (x y z) (declare (type exponent x y z)) (<= x (max y z))) m1 m2 m3))
    164 
    165 (defun monom-lcm-divides-monom-lcm-p (m1 m2 m3 m4)
    166   "Returns T if monomial MONOM-LCM(M1,M2) divides MONOM-LCM(M3,M4), NIL otherwise."
    167   (declare (type monom m1 m2 m3 m4))
    168   (every #'(lambda (x y z w) (declare (type exponent x y z w)) (<= (max x y) (max z w))) m1 m2 m3 m4))
    169 
    170 (defun monom-lcm-equal-monom-lcm-p (m1 m2 m3 m4)
    171   "Returns T if monomial MONOM-LCM(M1,M2) equals MONOM-LCM(M3,M4), NIL otherwise."
    172   (declare (type monom m1 m2 m3 m4))
    173   (every #'(lambda (x y z w) (declare (type exponent x y z w)) (= (max x y) (max z w))) m1 m2 m3 m4))
    174 
    175 (defun monom-divisible-by-p (m1 m2)
    176   "Returns T if monomial M1 is divisible by monomial M2, NIL otherwise."
    177   (declare (type monom m1 m2))
    178    (every #'>= m1 m2))
    179 
    180 (defun monom-rel-prime-p (m1 m2)
    181   "Returns T if two monomials M1 and M2 are relatively prime (disjoint)."
    182   (declare (type monom m1 m2))
    183   (every #'(lambda (x y) (declare (type exponent x y)) (zerop (min x y))) m1 m2))
    184 
    185 (defun monom-equal-p (m1 m2)
    186   "Returns T if two monomials M1 and M2 are equal."
    187   (declare (type monom m1 m2))
    188   (every #'= m1 m2))
    189 
    190 (defun monom-lcm (m1 m2 &aux (result (copy-seq m1)))
    191   "Returns least common multiple of monomials M1 and M2."
    192   (declare (type monom m1 m2))
    193   (map-into result #'max m1 m2))
    194 
    195 (defun monom-gcd (m1 m2 &aux (result (copy-seq m1)))
    196   "Returns greatest common divisor of monomials M1 and M2."
    197   (declare (type monom m1 m2))
    198   (map-into result #'min m1 m2))
    199 
    200 (defun monom-depends-p (m k)
    201   "Return T if the monomial M depends on variable number K."
    202   (declare (type monom m) (fixnum k))
    203   (plusp (elt m k)))
    204 
    205 (defmacro monom-map (fun m &rest ml &aux (result `(copy-seq ,m)))
    206   `(map-into ,result ,fun ,m ,@ml))
    207 
    208 (defmacro monom-append (m1 m2)
    209   `(concatenate 'monom ,m1 ,m2))
    210 
    211 (defmacro monom-contract (k m)
    212   `(subseq ,m ,k))
    213 
    214 (defun monom-exponents (m)
    215   (declare (type monom m))
    216   (coerce m 'list))
    217 
    218 
    219 
    220 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    221 ;;
    222 ;; Implementations of various admissible monomial orders
    223 ;;
    224 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    225 
    226 ;; pure lexicographic
    227 (defun lex> (p q &optional (start 0) (end (monom-dimension  p)))
    228   "Return T if P>Q with respect to lexicographic order, otherwise NIL.
    229 The second returned value is T if P=Q, otherwise it is NIL."
    230   (declare (type monom p q) (type fixnum start end))
    231   (do ((i start (1+ i)))
    232       ((>= i end) (values nil t))
    233     (declare (type fixnum i))
    234     (cond
    235      ((> (monom-elt p i) (monom-elt q i))
    236       (return-from lex> (values t nil)))
    237      ((< (monom-elt p i) (monom-elt q i))
    238       (return-from lex> (values nil nil))))))
    239 
    240 ;; total degree order , ties broken by lexicographic
    241 (defun grlex> (p q &optional (start 0) (end (monom-dimension  p)))
    242   "Return T if P>Q with respect to graded lexicographic order, otherwise NIL.
    243 The second returned value is T if P=Q, otherwise it is NIL."
    244   (declare (type monom p q) (type fixnum start end))
    245   (let ((d1 (monom-total-degree p start end))
    246         (d2 (monom-total-degree q start end)))
    247     (cond
    248       ((> d1 d2) (values t nil))
    249       ((< d1 d2) (values nil nil))
    250       (t
    251         (lex> p q start end)))))
    252 
    253 
    254 ;; total degree, ties broken by reverse lexicographic
    255 (defun grevlex> (p q &optional (start 0) (end (monom-dimension  p)))
    256   "Return T if P>Q with respect to graded reverse lexicographic order,
    257 NIL otherwise. The second returned value is T if P=Q, otherwise it is NIL."
    258   (declare (type monom p q) (type fixnum start end))
    259   (let ((d1 (monom-total-degree p start end))
    260         (d2 (monom-total-degree q start end)))
    261     (cond
    262      ((> d1 d2) (values t nil))
    263      ((< d1 d2) (values nil nil))
    264      (t
    265       (revlex> p q start end)))))
    266 
    267 
    268 ;; reverse lexicographic
    269 (defun revlex> (p q &optional (start 0) (end (monom-dimension  p)))
    270   "Return T if P>Q with respect to reverse lexicographic order, NIL
    271 otherwise.  The second returned value is T if P=Q, otherwise it is
    272 NIL. This is not and admissible monomial order because some sets do
    273 not have a minimal element. This order is useful in constructing other
    274 orders."
    275   (declare (type monom p q) (type fixnum start end))
    276   (do ((i (1- end) (1- i)))
    277       ((< i start) (values nil t))
    278     (declare (type fixnum i))
    279     (cond
    280      ((< (monom-elt p i) (monom-elt q i))
    281       (return-from revlex> (values t nil)))
    282      ((> (monom-elt p i) (monom-elt q i))
    283       (return-from revlex> (values nil nil))))))
    284 
    285 
    286 (defun invlex> (p q &optional (start 0) (end (monom-dimension  p)))
    287   "Return T if P>Q with respect to inverse lexicographic order, NIL otherwise
    288 The second returned value is T if P=Q, otherwise it is NIL."
    289   (declare (type monom p q) (type fixnum start end))
    290   (do ((i (1- end) (1- i)))
    291       ((< i start) (values nil t))
    292     (declare (type fixnum i))
    293       (cond
    294          ((> (monom-elt p i) (monom-elt q i))
    295           (return-from invlex> (values t nil)))
    296          ((< (monom-elt p i) (monom-elt q i))
    297           (return-from invlex> (values nil nil))))))
    298 
    299 
    300 
    301 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    302 ;;
    303 ;; Order making functions
    304 ;;
    305 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    306 
    307 (defvar *monomial-order* #'lex>
    308   "Default order for monomial comparisons")
    309 
    310 (defmacro monomial-order (x y)
    311   `(funcall *monomial-order* ,x ,y))
    312 
    313 (defun reverse-monomial-order (x y)
    314   (monomial-order y x))
    315 
    316 (defvar *primary-elimination-order* #'lex>)
    317 
    318 (defvar *secondary-elimination-order* #'lex>)
    319 
    320 (defvar *elimination-order* nil
    321   "Default elimination order used in elimination-based functions.
    322 If not NIL, it is assumed to be a proper elimination order. If NIL,
    323 we will construct an elimination order using the values of
    324 *PRIMARY-ELIMINATION-ORDER* and *SECONDARY-ELIMINATION-ORDER*.")
    325 
    326 (defun elimination-order (k)
    327   "Return a predicate which compares monomials according to the
    328 K-th elimination order. Two variables *PRIMARY-ELIMINATION-ORDER*
    329 and *SECONDARY-ELIMINATION-ORDER* control the behavior on the first K
    330 and the remaining variables, respectively."
    331   (declare (type fixnum k))
    332   #'(lambda (p q &optional (start 0) (end (monom-dimension  p)))
    333       (declare (type monom p q) (type fixnum start end))
    334       (multiple-value-bind (primary equal)
    335            (funcall *primary-elimination-order* p q start k)
    336          (if equal
    337              (funcall *secondary-elimination-order* p q k end)
    338            (values primary nil)))))
    339 
    340 (defun elimination-order-1 (p q &optional (start 0) (end (monom-dimension  p)))
    341   "Equivalent to the function returned by the call to (ELIMINATION-ORDER 1)."
    342   (declare (type monom p q) (type fixnum start end))
    343   (cond
    344    ((> (monom-elt p start) (monom-elt q start)) (values t nil))
    345    ((< (monom-elt p start) (monom-elt q start)) (values nil nil))
    346    (t (funcall *secondary-elimination-order* p q (1+ start) end))))
    347 
    348 
    349 
    350 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    351 ;;
    352 ;; Priority queue stuff
    353 ;;
    354 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    355 
    356 (defparameter *priority-queue-allocation-size* 16)
    357 
    358 (defun priority-queue-make-heap (&key (element-type 'fixnum))
    359   (make-array *priority-queue-allocation-size* :element-type element-type :fill-pointer 1
    360               :adjustable t))
    361 
    362 (defstruct (priority-queue (:constructor priority-queue-construct))
    363   (heap (priority-queue-make-heap))
    364   test)
    365 
    366 (defun make-priority-queue (&key (element-type 'fixnum)
    367                             (test #'<=)
    368                             (element-key #'identity))
    369   (priority-queue-construct
    370    :heap (priority-queue-make-heap :element-type element-type)
    371    :test #'(lambda (x y) (funcall test (funcall element-key y) (funcall element-key x)))))
    372  
    373 (defun priority-queue-insert (pq item)
    374   (priority-queue-heap-insert (priority-queue-heap pq) item (priority-queue-test pq)))
    375 
    376 (defun priority-queue-remove (pq)
    377   (priority-queue-heap-remove (priority-queue-heap pq) (priority-queue-test pq)))
    378 
    379 (defun priority-queue-empty-p (pq)
    380   (priority-queue-heap-empty-p (priority-queue-heap pq)))
    381 
    382 (defun priority-queue-size (pq)
    383   (fill-pointer (priority-queue-heap pq)))
    384 
    385 (defun priority-queue-upheap (a k
    386                &optional
    387                (test #'<=)
    388                &aux  (v (aref a k)))
    389   (declare (fixnum k))
    390   (assert (< 0 k (fill-pointer a)))
    391   (loop
    392    (let ((parent (ash k -1)))
    393      (when (zerop parent) (return))
    394      (unless (funcall test (aref a parent) v)
    395        (return))
    396      (setf (aref a k) (aref a parent)
    397            k parent)))
    398   (setf (aref a k) v)
    399   a)
    400 
    401    
    402 (defun priority-queue-heap-insert (a item &optional (test #'<=))
    403   (vector-push-extend item a)
    404   (priority-queue-upheap a (1- (fill-pointer a)) test))
    405 
    406 (defun priority-queue-downheap (a k
    407                  &optional
    408                  (test #'<=)
    409                  &aux  (v (aref a k)) (j 0) (n (fill-pointer a)))
    410   (declare (fixnum k n j))
    411   (loop
    412    (unless (<= k (ash n -1))
    413      (return))
    414    (setf j (ash k 1))
    415    (if (and (< j n) (not (funcall test (aref a (1+ j)) (aref a j))))
    416        (incf j))
    417    (when (funcall test (aref a j) v)
    418      (return))
    419    (setf (aref a k) (aref a j)
    420          k j))
    421   (setf (aref a k) v)
    422   a)
    423 
    424 (defun priority-queue-heap-remove (a &optional (test #'<=) &aux (v (aref a 1)))
    425   (when (<= (fill-pointer a) 1) (error "Empty queue."))
    426   (setf (aref a 1) (vector-pop a))
    427   (priority-queue-downheap a 1 test)
    428   (values v a))
    429 
    430 (defun priority-queue-heap-empty-p (a)
    431   (<= (fill-pointer a) 1))
     35
     36
    43237
    43338
    … …  
    530135;;
    531136;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    532 
    533 (defun make-term-variable (ring nvars pos
    534                                 &optional
    535                                 (power 1)
    536                                 (coeff (funcall (ring-unit ring)))
    537                                 &aux
    538                                 (monom (make-monom nvars :initial-element 0)))
    539   (declare (fixnum nvars pos power))
    540   (incf (monom-elt monom pos) power)
    541   (make-term monom coeff))
    542 
    543 (defstruct (term
    544             (:constructor make-term (monom coeff))
    545             ;;(:constructor make-term-variable)
    546             ;;(:type list)
    547             )
    548   (monom (make-monom 0) :type monom)
    549   (coeff nil))
    550 
    551 (defun term-sugar (term)
    552   (monom-sugar (term-monom term)))
    553 
    554 (defun termlist-sugar (p &aux (sugar -1))
    555   (declare (fixnum sugar))
    556   (dolist (term p sugar)
    557     (setf sugar (max sugar (term-sugar term)))))
    558137
    559138
    … …  
    703282
    704283
    705 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    706 ;;
    707 ;; Additional structure operations on a list of terms
    708 ;;
    709 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    710 
    711 (defun termlist-contract (p &optional (k 1))
    712   "Eliminate first K variables from a polynomial P."
    713   (mapcar #'(lambda (term) (make-term (monom-contract k (term-monom term))
    714                                       (term-coeff term)))
    715           p))
    716 
    717 (defun termlist-extend (p &optional (m (make-monom 1 :initial-element 0)))
    718   "Extend every monomial in a polynomial P by inserting at the
    719 beginning of every monomial the list of powers M."
    720   (mapcar #'(lambda (term) (make-term (monom-append m (term-monom term))
    721                                       (term-coeff term)))
    722           p))
    723 
    724 (defun termlist-add-variables (p n)
    725   "Add N variables to a polynomial P by inserting zero powers
    726 at the beginning of each monomial."
    727   (declare (fixnum n))
    728   (mapcar #'(lambda (term)
    729               (make-term (monom-append (make-monom n :initial-element 0)
    730                                        (term-monom term))
    731                          (term-coeff term)))
    732           p))
    733 
    734 
    735 
    736 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    737 ;;
    738 ;; Arithmetic on polynomials
    739 ;;
    740 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    741 
    742 (defstruct (poly
    743             ;;BOA constructor, by default constructs zero polynomial
    744             (:constructor make-poly-from-termlist (termlist &optional (sugar (termlist-sugar termlist))))
    745             (:constructor make-poly-zero (&aux (termlist nil) (sugar -1)))
    746             ;;Constructor of polynomials representing a variable
    747             (:constructor make-variable (ring nvars pos &optional (power 1)
    748                                               &aux
    749                                               (termlist (list
    750                                                          (make-term-variable ring nvars pos power)))
    751                                               (sugar power)))
    752             (:constructor poly-unit (ring dimension
    753                                      &aux
    754                                      (termlist (termlist-unit ring dimension))
    755                                      (sugar 0))))
    756   (termlist nil :type list)
    757   (sugar -1 :type fixnum))
    758 
    759 ;; Leading term
    760 (defmacro poly-lt (p) `(car (poly-termlist ,p)))
    761 
    762 ;; Second term
    763 (defmacro poly-second-lt (p) `(cadar (poly-termlist ,p)))
    764 
    765 ;; Leading monomial
    766 (defun poly-lm (p) (term-monom (poly-lt p)))
    767 
    768 ;; Second monomial
    769 (defun poly-second-lm (p) (term-monom (poly-second-lt p)))
    770 
    771 ;; Leading coefficient
    772 (defun poly-lc (p) (term-coeff (poly-lt p)))
    773 
    774 ;; Second coefficient
    775 (defun poly-second-lc (p) (term-coeff (poly-second-lt p)))
    776 
    777 ;; Testing for a zero polynomial
    778 (defun poly-zerop (p) (null (poly-termlist p)))
    779 
    780 ;; The number of terms
    781 (defun poly-length (p) (length (poly-termlist p)))
    782 
    783 (defun scalar-times-poly (ring c p)
    784   (make-poly-from-termlist (scalar-times-termlist ring c (poly-termlist p)) (poly-sugar p)))
    785 
    786 ;; The scalar product omitting the head term
    787 (defun scalar-times-poly-1 (ring c p)
    788   (make-poly-from-termlist (scalar-times-termlist ring c (cdr (poly-termlist p))) (poly-sugar p)))
    789    
    790 (defun monom-times-poly (m p)
    791   (make-poly-from-termlist (monom-times-termlist m (poly-termlist p)) (+ (poly-sugar p) (monom-sugar m))))
    792 
    793 (defun term-times-poly (ring term p)
    794   (make-poly-from-termlist (term-times-termlist ring term (poly-termlist p)) (+ (poly-sugar p) (term-sugar term))))
    795 
    796 (defun poly-add (ring p q)
    797   (make-poly-from-termlist (termlist-add ring (poly-termlist p) (poly-termlist q)) (max (poly-sugar p) (poly-sugar q))))
    798 
    799 (defun poly-sub (ring p q)
    800   (make-poly-from-termlist (termlist-sub ring (poly-termlist p) (poly-termlist q)) (max (poly-sugar p) (poly-sugar q))))
    801 
    802 (defun poly-uminus (ring p)
    803   (make-poly-from-termlist (termlist-uminus ring (poly-termlist p)) (poly-sugar p)))
    804 
    805 (defun poly-mul (ring p q)
    806   (make-poly-from-termlist (termlist-mul ring (poly-termlist p) (poly-termlist q)) (+ (poly-sugar p) (poly-sugar q))))
    807 
    808 (defun poly-expt (ring p n)
    809   (make-poly-from-termlist (termlist-expt ring (poly-termlist p) n) (* n (poly-sugar p))))
    810 
    811 (defun poly-append (&rest plist)
    812   (make-poly-from-termlist (apply #'append (mapcar #'poly-termlist plist))
    813              (apply #'max (mapcar #'poly-sugar plist))))
    814 
    815 (defun poly-nreverse (p)
    816   (setf (poly-termlist p) (nreverse (poly-termlist p)))
    817   p)
    818 
    819 (defun poly-contract (p &optional (k 1))
    820   (make-poly-from-termlist (termlist-contract (poly-termlist p) k)
    821              (poly-sugar p)))
    822 
    823 (defun poly-extend (p &optional (m (make-monom 1 :initial-element 0)))
    824   (make-poly-from-termlist
    825    (termlist-extend (poly-termlist p) m)
    826    (+ (poly-sugar p) (monom-sugar m))))
    827 
    828 (defun poly-add-variables (p k)
    829   (setf (poly-termlist p) (termlist-add-variables (poly-termlist p) k))
    830   p)
    831 
    832 (defun poly-list-add-variables (plist k)
    833   (mapcar #'(lambda (p) (poly-add-variables p k)) plist))
    834 
    835 (defun poly-standard-extension (plist &aux (k (length plist)))
    836   "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]."
    837   (declare (list plist) (fixnum k))
    838   (labels ((incf-power (g i)
    839              (dolist (x (poly-termlist g))
    840                (incf (monom-elt (term-monom x) i)))
    841              (incf (poly-sugar g))))
    842     (setf plist (poly-list-add-variables plist k))
    843     (dotimes (i k plist)
    844       (incf-power (nth i plist) i))))
    845 
    846 (defun saturation-extension (ring f plist &aux (k (length plist)) (d (monom-dimension (poly-lm (car plist)))))
    847   "Calculate [F, U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK]."
    848   (setf f (poly-list-add-variables f k)
    849         plist (mapcar #'(lambda (x)
    850                           (setf (poly-termlist x) (nconc (poly-termlist x)
    851                                                          (list (make-term (make-monom d :initial-element 0)
    852                                                                           (funcall (ring-uminus ring) (funcall (ring-unit ring)))))))
    853                           x)
    854                       (poly-standard-extension plist)))
    855   (append f plist))
    856 
    857 
    858 (defun polysaturation-extension (ring f plist &aux (k (length plist))
    859                                                    (d (+ k (length (poly-lm (car plist))))))
    860   "Calculate [F, U1*P1+U2*P2+...+UK*PK-1], where PLIST=[P1,P2,...,PK]."
    861   (setf f (poly-list-add-variables f k)
    862         plist (apply #'poly-append (poly-standard-extension plist))
    863         (cdr (last (poly-termlist plist))) (list (make-term (make-monom d :initial-element 0)
    864                                                             (funcall (ring-uminus ring) (funcall (ring-unit ring))))))
    865   (append f (list plist)))
    866 
    867 (defun saturation-extension-1 (ring f p) (polysaturation-extension ring f (list p)))
    868 
    869 
    870 
    871 
    872 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    873 ;;
    874 ;; Evaluation of polynomial (prefix) expressions
    875 ;;
    876 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    877 
    878 (defun coerce-coeff (ring expr vars)
    879   "Coerce an element of the coefficient ring to a constant polynomial."
    880   ;; Modular arithmetic handler by rat
    881   (make-poly-from-termlist (list (make-term (make-monom (length vars) :initial-element 0)
    882                               (funcall (ring-parse ring) expr)))
    883              0))
    884 
    885 (defun poly-eval (ring expr vars &optional (list-marker '[))
    886   (labels ((p-eval (arg) (poly-eval ring arg vars))
    887            (p-eval-list (args) (mapcar #'p-eval args))
    888            (p-add (x y) (poly-add ring x y)))
    889     (cond
    890      ((eql expr 0) (make-poly-zero))
    891      ((member expr vars :test #'equalp)
    892       (let ((pos (position expr vars :test #'equalp)))
    893         (make-variable ring (length vars) pos)))
    894      ((atom expr)
    895       (coerce-coeff ring expr vars))
    896      ((eq (car expr) list-marker)
    897       (cons list-marker (p-eval-list (cdr expr))))
    898      (t
    899       (case (car expr)
    900         (+ (reduce #'p-add (p-eval-list (cdr expr))))
    901         (- (case (length expr)
    902              (1 (make-poly-zero))
    903              (2 (poly-uminus ring (p-eval (cadr expr))))
    904              (3 (poly-sub ring (p-eval (cadr expr)) (p-eval (caddr expr))))
    905              (otherwise (poly-sub ring (p-eval (cadr expr))
    906                                   (reduce #'p-add (p-eval-list (cddr expr)))))))
    907         (*
    908          (if (endp (cddr expr))         ;unary
    909              (p-eval (cdr expr))
    910            (reduce #'(lambda (p q) (poly-mul ring p q)) (p-eval-list (cdr expr)))))
    911         (expt
    912          (cond
    913           ((member (cadr expr) vars :test #'equalp)
    914            ;;Special handling of (expt var pow)
    915            (let ((pos (position (cadr expr) vars :test #'equalp)))
    916              (make-variable ring (length vars) pos (caddr expr))))
    917           ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
    918            ;; Negative power means division in coefficient ring
    919            ;; Non-integer power means non-polynomial coefficient
    920            (coerce-coeff ring expr vars))
    921           (t (poly-expt ring (p-eval (cadr expr)) (caddr expr)))))
    922         (otherwise
    923          (coerce-coeff ring expr vars)))))))
    924 
    925284
    926285
    … …  
    931290;;
    932291;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    933 
    934 
    935 
    936292(defmacro debug-cgb (&rest args)
    937293  `(when $poly_grobner_debug (format *terminal-io* ,@args)))
Note: See TracChangeset for help on using the changeset viewer.