1;;; -*- Package: MAXIMA; Base: 10.; Syntax: Common-lisp -*- 2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3;;; ;;;;; 4;;; Copyright (c) 1984 by William Schelter,University of Texas ;;;;; 5;;; All rights reserved ;;;;; 6;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 7 8(in-package :maxima) 9 10(eval-when 11 #+gcl (compile load eval) 12 #-gcl (:compile-toplevel :load-toplevel :execute) 13 14 (defmacro p-cof (f) ;leading coefficient of f 15 `(third ,f)) 16 17 (defmacro p-next-term (f) 18 `(cddr ,f)) 19 20 (defmacro p-deg (f) 21 `(second ,f)) 22 23 (defmacro term-cof (terms) 24 `(second ,terms)) 25 26 (defmacro term-deg (terms) 27 `(first ,terms))) 28 29;;(make-poly var x)==> (list x 1 1) 30 31(defun make-polynomial (&key var terms) 32 (psimp var terms)) 33 34(defun afc-quotient (f g) 35 (cquotient f g)) 36 37(defun fsignal (&rest l) 38 (error (car l))) 39 40(defmacro working-modulo (list-of-monic-polynomials &body body &aux (old-tellrats (make-symbol "old-tellrats"))) 41 "The computations in body are done modulo the list-of-monic-polynomials. The 42 results of squareing,multiplication, and exponentiating should be of lower degree in each of the 43 monic polynomials than the degree of the monic polynomial" 44 `(let ((,old-tellrats (list-previous-tellrats ,list-of-monic-polynomials))) 45 (unwind-protect 46 (progn 47 (set-tellrats ,list-of-monic-polynomials) 48 ,@body) 49 (undo-tellrats ,old-tellrats)))) 50 51;;sample usage 52;;(working-modulo (st-rat #$[x^2+x+1,y^4+y]$) (ptimes ...)) 53 54(defun list-previous-tellrats (new-tellrats) 55 (loop for v in new-tellrats 56 collecting (cons (car v) (get (car v) 'tellrat)))) 57 58(defun set-tellrats (new-tellrats) 59 (loop for v in new-tellrats 60 do (putprop (car v) (cdr v) 'tellrat))) 61 62(defun undo-tellrats (old-list) 63 (loop for v in old-list 64 when (null (cdr v)) 65 do (remprop (car v) 'tellrat) 66 else do (putprop (car v) (cdr v) 'tellrat))) 67 68;;version for possibly some zero terms resulting 69;;using nconc less space but slower since (list a b) is cdr coded and it has to fix up. 70(defmacro term-operation (f g operation-function &optional deg-shift-result) 71 "If f and g are polynomials this constructs polynomial whose main variable is the same as f ~ 72 with coefficients (operation-function f_i g)..[f_i means the i'th coefficient of f 73 degree (+ i deg-shift-result)" 74 (let ((cof (make-symbol "cof")) 75 (deg (make-symbol "deg")) 76 (tem (gensym))) 77 `(psimp (p-var ,f) 78 (loop for (,deg ,cof) on (cdr ,f) by #'cddr 79 with ,tem 80 do (setq ,tem (,operation-function ,cof ,g)) 81 when (not (pzerop ,tem)) 82 nconc (list , (cond (deg-shift-result `(+ ,deg ,deg-shift-result)) 83 (t deg)) 84 ,tem))))) 85 86(defmacro plain-term-operation (terms-f terms-g operation-function &optional deg-shift-result) 87 "constructs terms of a polynomial with coefficients (operation-function f_i terms-g) in 88 degree (+ i deg-shift-result)" 89 (let ((cof (make-symbol "cof")) 90 (deg (make-symbol "deg")) 91 (tem (gensym))) 92 `(loop for (,deg ,cof) on ,terms-f by #'cddr 93 with ,tem 94 do (setq ,tem (,operation-function ,cof ,terms-g)) 95 when (not (pzerop ,tem)) 96 nconc (list ,(cond (deg-shift-result `(+ ,deg ,deg-shift-result)) 97 (t deg)) 98 ,tem)))) 99 100(defmacro quotient-not-exact () 101 `(cond ((and (boundp '*testing*) 102 *testing*) 103 (throw 'quotient-not-exact t)) 104 (t (fsignal 'quotient-not-exact)))) 105 106 107(DEFUN afp-CQUOTIENT (A B) 108 (declare (special *testing*)) 109 (COND ((EQl A 0) 0) 110 ((NULL MODULUS) 111 (let ((quot (quotient a b))) 112 (cond ((eql (* quot b) a) 113 quot) 114 (t (quotient-not-exact))))) 115 (t (ctimes a (crecip b))))) 116 117 118;;The following works to eliminate zero terms when using modulus etc. 119;;It is also faster than ptimes by a slight margin. 120;;unlike ptimes which may introduce zero terms eg: 121;;(let ((modulus 4))(ptimes (x+y+z)^4 (x+y+z)^4) gives a bad result 122 123 124 125(defun afp-terms-times (terms-f terms-g &aux answ (g-orig terms-g) prod-exp tail prod-cof) 126 "Returns the terms of the polynomial which is the product of the two polynomials ~ 127 whose terms are terms-f and terms-g. If modulus is in effect the the result will have its ~ 128 coefficients reduced by modulus." 129 (cond 130 (terms-g (setq answ (plain-term-operation terms-g (term-cof terms-f) afp-times (term-deg terms-f))) 131 (loop for (f-exp f-cof ) on (cddr terms-f) by #'cddr 132 do 133 (setq terms-g g-orig) 134 (prog () 135 first-product 136 (cond ((null terms-g)(return nil))) 137 (setq prod-exp (+ f-exp (term-deg terms-g))) 138 (setq prod-cof (afp-times(term-cof terms-g) f-cof)) 139 (cond ((pzerop prod-cof) (setq terms-g (cddr terms-g)) (go first-product)) 140 ((or (null answ) (> prod-exp (term-deg answ))) 141 (setq answ (cons prod-exp (cons prod-cof answ))) 142 (setq tail (cdr answ)) (go next-product)) 143 ((eql prod-exp (term-deg answ)) 144 (setq prod-cof (pplus (term-cof answ) prod-cof)) 145 (cond ((pzerop prod-cof) (setq answ (cddr answ)) 146 (setq terms-g (cddr terms-g)) (go first-product)) 147 (t (setf (term-cof answ) prod-cof) 148 (setq tail (cdr answ)) (go next-product))))) 149 ;;below here assume answ not empty and (term-deg answ) 150 ;;greater any possible future prod-exp (until next 151 ;;f-exp,f-cof) Once below this point we stay below, 152 ;;until next f-exp,f-cof . Tail begins with the 153 ;;coefficient whose corresponding degree is definitely 154 ;;higher than any prod-exp to be encountered. 155 (setq tail (cdr answ)) 156 tail-certain 157 (cond ((and (cdr tail)(> (second tail) prod-exp)) 158 (setq tail (cddr tail)) (go tail-certain))) 159 (cond ((or (null (cdr tail))(< (second tail) prod-exp)) 160 (setf (cdr tail) (cons prod-exp (cons prod-cof (cdr tail)))) 161 (setq tail (cddr tail)) (go next-product))) 162 (cond ((pzerop (setq prod-cof (pplus (third tail) prod-cof))) 163 (setf (cdr tail) (cdddr tail))) 164 (t (setf (third tail) prod-cof) (setq tail (cddr tail)))) 165 next-product 166 (setq terms-g (cddr terms-g)) 167 (cond ((null terms-g) (return nil))) 168 (setq prod-exp (+ f-exp (car terms-g))) 169 (setq prod-cof (afp-times (second terms-g) f-cof)) 170 (cond ((pzerop prod-cof) (go next-product))) 171 (go tail-certain))))) 172 answ) 173 174(defmacro afp-main-plus-non-main (constant f-main) 175 "Adds a polynomial CONSTANT to a polynomial whose main variable is higher than 176 any in F-MAIN" 177 `(psimp (p-var ,f-main) (afp-constant-term-plus ,constant (cdr ,f-main)))) 178 179(defmacro afp-number-plus (number poly) 180 "Adds a NUMBER to a polynomial POLY, returning the result" 181 `(cond ((numberp ,poly)(cplus ,number ,poly)) 182 (t (afp-main-plus-non-main ,number ,poly)))) 183 184(defun afp-plus (f g) 185 "Returns the sum of the two polynomials f and g" 186 (cond ((numberp f) (afp-number-plus f g)) 187 ((numberp g) (afp-number-plus g f)) 188 ((eq (p-var f) (p-var g)) 189 (psimp (p-var f) 190 (afp-terms-plus (cdr f) (cdr g)))) 191 ((pointergp (p-var f) (p-var g)) (afp-main-plus-non-main g f)) 192 (t (afp-main-plus-non-main f g)))) 193 194(defun afp-constant-term-plus (constant terms) 195 "Adds a polynomial (CONSTANT) not involving the main variable of a polynomial whose 196 terms are TERMS. Naturally the main variable is assumed higher than any in the CONSTANT. ~ 197 The result is the terms of the sum polynomial" 198 (cond ((pzerop constant) terms) 199 ((null terms) 200 (list 0 constant)) 201 ((zerop (car terms))(setq constant (afp-plus constant (second terms))) 202 (cond ((pzerop constant)nil) 203 (t (list 0 constant)))) 204 (t (cons (car terms) (cons (second terms) (afp-constant-term-plus constant (cddr terms))))))) 205 206 207(defun afp-terms-plus (terms-f terms-g &aux e) 208 "Returns the terms of the polynomial which is the sum of f and g 209 if the terms of f and the terms of g are the two arguments." 210 (cond ((null terms-f) terms-g) 211 ((null terms-g) terms-f) 212 ((eql (car terms-f) (car terms-g)) 213 (setq e (afp-plus (second terms-f) (second terms-g))) 214 (cond ((pzerop e)(afp-terms-plus (cddr terms-f) (cddr terms-g))) 215 (t (cons (car terms-f) (cons e (afp-terms-plus (cddr terms-f) (cddr terms-g))))))) 216 ((> (car terms-f) (car terms-g ))(cons (car terms-f) (cons (second terms-f) 217 (afp-terms-plus terms-g (cddr terms-f))))) 218 (t(cons (car terms-g) (cons (second terms-g) 219 (afp-terms-plus terms-f (cddr terms-g))))))) 220 221 222(defmacro qfirstn (n l) 223 (case n 224 (1 `(list (car ,l))) 225 (2 `(list (car ,l) (second ,l))) 226 (t `(subseq ,l 0 ,n)))) 227 228(defun afp-minus (f) 229 "makes no check that keeping in the modulus range nor that the result is reduced, 230 but the result g will satisfy (afp-plus f g) ==> 0" 231 (cond ((numberp f) (cminus f)) 232 (t (cons (car f) (afp-terms-minus (cdr f)))))) 233 234(defun afp-terms-minus (terms-f) 235 (loop for (deg pol) on terms-f by #'cddr nconc (list deg (afp-minus pol)))) 236 237(defmacro add-one-term (deg cof terms) 238 (cond ((symbolp cof) 239 `(cond ((pzerop ,cof) ,terms) 240 (t (cons ,deg (cons ,cof ,terms))))) 241 (t 242 `(let ((.cof. ,cof)) 243 (cond ((pzerop .cof.) ,terms) 244 (t (cons ,deg (cons .cof. ,terms)))))))) 245 246(defun afp-difference (f g) 247 (cond ((numberp f) 248 (cond ((numberp g)(cdifference f g)) 249 (t (psimp (p-var g) (afp-terms-constant-main-differ f (cdr g)))))) 250 ((numberp g) 251 (psimp (p-var f) (afp-terms-main-constant-differ (cdr f) g))) 252 ((eq (p-var f) (p-var g)) 253 (psimp (p-var f) (afp-terms-differ (p-terms f) (p-terms g)))) 254 ((pointergp (p-var f) (p-var g)) 255 (psimp (p-var f) (afp-terms-main-constant-differ 256 (cdr f) g))) 257 (t(psimp (p-var g) (afp-terms-main-constant-differ 258 f (cdr g)))))) 259 260 261 262(defun afp-terms-differ (terms-f terms-g) 263 (cond ((null terms-f) (afp-terms-minus terms-g)) 264 ((null terms-g) terms-f) 265 ((eql (term-deg terms-f) (term-deg terms-g)) 266 (add-one-term (term-deg terms-f) (afp-difference (term-cof terms-f) (term-cof terms-g)) 267 (afp-terms-differ (cddr terms-f) (cddr terms-g)))) 268 ((> (term-deg terms-f) (term-deg terms-g)) 269 (cons (term-deg terms-f) (cons (term-cof terms-f) (afp-terms-differ (p-next-term terms-f) terms-g)))) 270 (t 271 (cons (term-deg terms-g) (cons (afp-minus (term-cof terms-g)) (afp-terms-differ terms-f (p-next-term terms-g))))))) 272 273(defun afp-terms-constant-main-differ (const main) 274 "main is terms const is a polynomial" 275 (cond ((null main)(cond ((pzerop const) nil) 276 (t (list 0 const)))) 277 ((zerop (car main)) 278 (add-one-term 0 (afp-difference (second main) const) nil)) 279 (t (cons (car main) (cons (afp-minus (second main)) (afp-terms-constant-main-differ const (cddr main))))))) 280 281(defun afp-terms-main-constant-differ (main const) 282 (cond ((null main)(cond ((pzerop const) nil) 283 (t (list 0 (afp-minus const))))) 284 ((zerop (car main)) 285 (add-one-term 0 (afp-difference (second main) const) nil)) 286 (t (cons (car main) (cons (second main) (afp-terms-main-constant-differ (cddr main) const)))))) 287 288;;assumes divides evenly, integer coefficients. 289;;signal error otherwise the flavor of the error should be specified. 290(defun afp-quotient (f g) 291 "Tries to divide the polynomial f by g. One has result*g=f, or else ~ 292 it signals an error." 293 (declare (special *testing*)) 294 (cond ((numberp f) 295 (cond ((numberp g) 296 (afp-cquotient f g)) 297 ((zerop f) 0) 298 (t (quotient-not-exact)))) 299 ((numberp g) 300 (term-operation f g afp-quotient)) 301 ((pointergp (p-var f) (p-var g)) 302 (term-operation f g afp-quotient)) 303 ((eq (p-var f) (p-var g)) 304 (loop 305 with quot with deg-dif with main-var = (p-var f) with minus-quot 306 do 307 (setq deg-dif (- (p-deg f) (p-deg g))) 308 while (>= deg-dif 0) 309 collecting deg-dif into q 310 collecting (setq quot (afp-quotient (p-cof f) (p-cof g))) 311 into q 312 do (setq minus-quot (pminus quot)) 313 (setq f (pplus 314 f 315 (term-operation g 316 ;;-fn/gm 317 minus-quot 318 ptimes 319 deg-dif))) 320 while (not (and (numberp f) (zerop f))) 321 when (or (numberp f) (not (eq main-var (p-var f))) 322 (< (p-deg f) (p-deg g))) 323 do (quotient-not-exact) 324 finally (return (make-polynomial :terms q :var main-var)))) 325 (t (quotient-not-exact)))) 326 327 328(defun afp-test-divide (f g &aux ( *testing* t) quot) 329 (declare (special *testing*)) 330 (catch 'quotient-not-exact (setq quot (afp-quotient f g))) 331 (cond (quot quot) 332 (t nil))) 333 334(defun afc-remainder (n divisor) 335 (multiple-value-bind (q r) (truncate n divisor) (values r q))) 336 337;;pseudo division as in Knuth's book. 338 339(defun afp-pseudo-quotient (f g &aux (creqd 1)) 340 "This function returns the values: quotient, remainder and creqd ~ 341 so that creqd*f=g*quotient+remainder. Creqd does not involve the 342 main variable of f. The remainder has degree lower than g with respect 343 to the main variable of f." 344 (cond ((numberp f) 345 (cond ((numberp g) 346 (apply #'values (append '(1) (multiple-value-list (afc-remainder f g))))) 347 (t (values 0 f 1)))) 348 ((numberp g) 349 (values f 0 g)) 350 ((pointergp (p-var f) (p-var g)) 351 (values f 0 g)) 352 ((eq (p-var f) (p-var g)) 353 (loop with quot with deg-dif with main-var = (p-var f) 354 with remainder = 355 (cond ((and (numberp (p-cof g)) modulus) f) 356 ((and (numberp (p-cof g)) 357 (eql (abs (p-cof g)) 1)) f) 358 (t (ptimes f (setq creqd (pexpt (p-cof g) 359 (+ 1 (- (p-deg f) (p-deg g)))))))) 360 do 361 (setq deg-dif (- (p-deg remainder) (p-deg g))) 362 while (>= deg-dif 0) 363 collecting deg-dif into q 364 collecting (setq quot (afp-quotient (p-cof remainder) (p-cof g))) 365 into q 366 do 367 (setq remainder (pplus 368 remainder 369 (term-operation g 370 ;;-fn/gm 371 (pminus quot) 372 ptimes deg-dif))) 373 while (and (not (numberp remainder)) 374 (eql (p-var remainder) main-var)) 375 finally 376 (setq quot (make-polynomial :terms q :var main-var)) 377 ; (iassert (eql 0 (pdifference (ptimes f creqd ) (Pplus remainder (ptimes quot g))))) 378 (return (values quot remainder creqd)))) 379 (t (values 0 f 1)))) 380 381(defmacro assume-pointerg-or-equal-p (f g &optional reverse-flag) 382 `(cond (( gen-pointergp ,f ,g) nil) 383 (t 384 ,@ (cond (reverse-flag (list `(setq ,reverse-flag (null ,reverse-flag))))) 385 (rotatef ,f ,g)))) 386 387(defun gen-pointergp (f g) 388 (cond ((numberp g) t) 389 ((numberp f) nil) 390 (t (pointergp (p-var f) (p-var g))))) 391 392(defun gen-degree (f) 393 (cond ((numberp f) 0) 394 (t (p-deg f)))) 395 396(defmacro assume-greater-equal-degree (f g &optional reverse-flag) 397 `(cond ((>= (gen-degree ,f) (gen-degree ,g)) nil) 398 (t 399 ,@ (cond (reverse-flag (list `(setq (,reverse-flag (null ,reverse-flag)))))) 400 (rotatef ,f ,g)))) 401 402 403(defun afp-content (f ) 404 "Returns the gcd of the coefficients of f (with respect to the main variable) ~ 405 if f is a polynomial" 406 (cond ((numberp f) f) 407 (t (loop for (deg cof) on (cdddr f) by #'cddr 408 with cont = (p-cof f) 409 do (setq cont (afp-gcd cont cof)) 410 finally (return cont))))) 411 412(defun principal-part (f) 413 (afp-quotient f (afp-content f))) 414 415(defvar *afp-gcd* 'subresultant) 416 417(defun afp-gcd (f g &aux answer) 418 "Returns the gcd of its two arguments which may be any polynomial ~ 419 with integer coefficients. " 420 (setq answer 421 (case *afp-gcd* 422 (euclidean (afp-euclidean-gcd f g)) 423 (subresultant (afp-subresultant-gcd f g)) 424 (t (merror "~%The value of the switch *afp-gcd* ~A is not legal." *afp-gcd*)))) 425 (cond (modulus (afp-try-make-monic answer)) 426 (t answer))) 427 428;;This is the Euclidean gcd algorithm, as in Knuth. 429(defun afp-euclidean-gcd (f g &aux u v r d contf contg) 430 "Returns the gcd of its two arguments which may be any polynomial ~ 431 with integer coefficients. Currently uses the Euclidean algorithm, but should 432 have a switch" 433 (assume-pointerg-or-equal-p f g) 434 (cond ((numberp f)(gcd f g)) 435 ((gen-pointergp f g) 436 (afp-gcd (afp-content f) g)) 437; (loop for (deg cof) on (p-terms f) by #'cddr 438; with gcd = g 439; do (setq gcd (afp-gcd cof gcd)) 440; finally (return gcd))) 441 ((eq (p-var f) (p-var g)) 442 (setq d (afp-gcd 443 (setq contf (afp-content f))(setq contg (afp-content g)))) 444 (setq u (afp-quotient f contf)) 445 (setq v (afp-quotient g contg)) 446 (assume-greater-equal-degree u v) 447 (loop with unused 448 do (multiple-value-setq (unused r) (afp-pseudo-quotient u v)) 449 when (pzerop r) 450 do (return (ptimes d v)) 451 when (gen-pointergp f r) 452 do (return d) 453 do (setq u v) 454 (setq v (afp-quotient r (afp-content r))))) 455 (t (fsignal 'should-not-get-here)))) 456 457 458;;This was about twice as fast as the regular maxima pgcd on 459;;some simple functions: 121 msec. as opposed to 250 msec. 460;;this was with the afp-content in terms of the old afp-gcd. 461;; (afp-subresultant-gcd-compare (st-rat #$8*3*(x^2+1)*(x+1)*(z^2+2)*(x+2)$)(st-rat #$15*(x^2-1)*(x^2+1)*(y^2+4)*(z^4+2*z+1)$)) 462;(user:compare-recursive-functions 463(defun afp-subresultant-gcd (f g &aux delta u v r d contf contg answ) 464 "Returns the gcd of its two arguments which may be any polynomial ~ 465 with integer coefficients. It uses the subresultant algorithm" 466 (assume-pointerg-or-equal-p f g) 467 (cond ((numberp f)(gcd f g)) 468 ((gen-pointergp f g) 469 (cond ((pzerop g) f) 470 (t 471 (afp-subresultant-gcd (afp-content f) g)))) 472 ((eq (p-var f) (p-var g)) 473 (setq d (afp-subresultant-gcd 474 (setq contf (afp-content f)) 475 (setq contg (afp-content g)))) 476 (setq u (afp-quotient f contf)) 477 (setq v (afp-quotient g contg)) 478 (assume-greater-equal-degree u v) 479 (setq answ (loop with gg = 1 with h = 1 with g^delta with unused 480 do 481 (setq delta (- (p-deg u) (p-deg v))) 482 (multiple-value-setq (unused r) (afp-pseudo-quotient u v)) 483 when (pzerop r) 484 do (return (ptimes d (principal-part v))) 485 when (gen-pointergp f r) 486 do (return d) 487 do (setq u v) 488 (setq v (afp-quotient r (ptimes gg (pexpt h delta)))) 489 (setq gg (p-cof u)) 490 (setq g^delta (pexpt gg delta)) 491 (setq h (cond ((eql delta 1) g^delta) 492 ((> delta 1) (afp-quotient g^delta 493 (pexpt h (- delta 1)))) 494 ;;here delta=0 495 (t (ptimes g^delta h)))))) 496 (afp-try-make-monic answ)))) 497 498 499(defun one-ptimes (f g) 500 (cond ((eql f 1) g) 501 ((eql g 1) f) 502 (t (ptimes f g)))) 503 504(defun exponent-product (&rest alternating-factor-exponent-list) 505 "Exponents may be positive or negative, but assumes result is poly" 506 (loop for (fact deg) on alternating-factor-exponent-list by #'cddr 507 with numer = 1 with denom = 1 508 when (= deg 1) 509 do (setq numer (one-ptimes numer fact)) 510 else 511 when (> deg 1) 512 do (setq numer (one-ptimes numer (pexpt fact deg))) 513 else 514 when (= deg -1) 515 do (setq denom (one-ptimes denom fact)) 516 else 517 when (< deg -1) 518 do (setq denom (one-ptimes denom (pexpt fact (- deg)))) 519 finally (return (afp-quotient numer denom)))) 520 521(defun same-main-and-degree (f g) 522 (cond ((numberp f)(numberp g)) 523 ((numberp g)(numberp f)) 524 (t 525 (and (eq (car f) (car g)) 526 (eq (second f) (second g)))))) 527;;unfinished. 528;(defun afp-sqfr (f &aux deriv d) 529; (cond ((numberp f) f) 530; (t 531; (loop 532; do 533; (setq deriv (pderivative f (p-var f))) 534; (setq d (afp-gcd f deriv)) 535; (cond ((numberp d ) f) 536; ((same-main-and-degree d f) 537; (make-polynomial :var (p-var f) 538; :teerms 539; (loop for (deg cof) on (cdr f) by #'cddr 540; collecting (quotient deg modulus) 541; collecting cof))) 542; (t (setq f (afp-quotient f d)))))))) 543 544(defun afp-big-gcd (f g &aux tem) 545 "The arguments may be polynomials with integer coefficients. ~ 546 Three values are returned: gcd , f/gcd, and g/gcd." 547 (values (setq tem (afp-subresultant-gcd f g)) (afp-quotient f tem) (afp-quotient g tem))) 548 549(defun afp-pgcdcofacts (f g) 550 (multiple-value-list (afp-big-gcd f g))) 551;;this used 3 times the space and was much slower for the (x+y+z)^10 552;;problem but if I put z=1 then it went much faster, and used less 553;;space, but still not as good as the subresultant algorithm. It was 554;;about as fast as my subresultant algorithm. 555;(defun afp-big-gcd (f g) 556; (declare (values gcd-of-f-g f-over-gcd g-over-gcd)) 557; (apply 'values (pgcdcofacts f g))) 558 559#+debug 560(defun test (f g) 561 (multiple-value-bind (d qf qg) 562 (afp-big-gcd f g) 563 (iassert (equal f (ptimes d qf))) (iassert (equal g (ptimes d qg))))) 564 565(defun afp-square-free-with-modulus (u &aux (vari (list-variables u)) root deriv quot agcd) 566 (check-arg vari (null (cdr vari)) "univariate when modulus") 567 (check-arg modulus (primep modulus) "prime") 568 (cond ((numberp u) u) 569 (t (setq agcd (afp-gcd u (setq deriv (pderivative u (p-var u))))) 570 (cond ((numberp agcd) (list u 1)) 571 ((> (p-deg u) (p-deg agcd)) 572 (setq quot (afp-quotient u agcd)) 573 (append (afp-square-free-with-modulus quot) (afp-square-free-with-modulus 574 agcd))) 575 ;;deriv is 0 576 (t (check-arg deriv (eql 0 deriv) "zero") 577 (setq root (psimp (p-var u) 578 (loop for (deg cof) on (cdr u) by #'cddr 579 collecting (quotient deg modulus) 580 collecting cof))) 581 (loop for (pol deg) on (afp-square-free-with-modulus root) by #'cddr 582 collecting pol collecting (* deg modulus))))))) 583 584;;timing on factoring the (x+y+z)^10 2.6 sec 10,047 words 585;;timing on factoring the (x+y+z)^20 20.2 sec 37,697 words 586 587(defun afp-square-free-factorization (u &aux d tx v1 w1 some-factors unit) 588 "returns an alternating list of factors and exponents. In the characteristic 0 case each factor is ~ 589 guaranteed square free and relatively prime to the other factors. Note that if 590 modulus is not zero then u should be univariate. Otherwise for eg mod 3, x^3-t 591 is not square free in the field with cube root of t adjoined, and it can't be factored 592 in Z/3Z[x,t]." 593 (cond (modulus 594 (afp-square-free-with-modulus u)) 595 (t 596 (cond ((numberp u) 597 u) 598 (t 599 (setq d (afp-content u)) 600 (setq u (term-operation u d afp-quotient)) 601 (multiple-value-setq (tx v1 w1) 602 (afp-big-gcd u (pderivative u (p-var u)))) 603 ; (show tx v1 w1) 604 (setq some-factors 605 (cond ((eql tx 1) 606 (list u 1)) 607 ((numberp tx) 608 (fsignal 'how-did-this-happen)) 609 (t 610 (loop for i from 1 611 with vi = v1 with wi = w1 with videriv with vi+1 with ui with wi+1 612 with main-var = (p-var u) 613 do ; (show i) 614 (setq videriv (pderivative vi main-var)) 615 ; (show factor-list) 616 when (equal wi videriv) 617 do 618 (return (append factor-list (list vi i))) 619 do 620 (multiple-value-setq (ui vi+1 wi+1) 621 (afp-big-gcd vi (pdifference wi 622 videriv))) 623 ; (show vi wi ui vi+1 wi+1) 624 when (not (eql ui 1)) 625 nconc (list ui i) into factor-list 626 do 627 (setq vi vi+1) 628 (setq wi wi+1) 629 )))) 630 ; (show some-factors) 631 ;;this is all to collect some numbers and fix the unit multiple. 632 (loop for (pol deg) on some-factors by #'cddr 633 with answ = d 634 do (setq answ (ptimes answ (pexpt (p-cof pol) deg))) 635 finally 636 (setq unit (afp-quotient (p-cof u) answ)) 637 (cond ((eql unit 1) nil) 638 ((eql unit -1) (loop for (pol1 deg1) on some-factors by #'cddr 639 for i from 0 by 2 640 when (oddp deg1) do (setf (nth i some-factors) 641 (pminus pol1)) 642 (return 'done) 643 finally (merror "no odd factors yet differs by minus "))) 644 (t (fsignal "not handled yet")))) 645 (cond ((eql d 1) some-factors) 646 (t (append (afp-square-free-factorization d) some-factors)))))))) 647 648 649 650;;tested on (x+y+z)^10 times itself and got about 651;;the same time as ptimes 3100 milliseconds. in temporary area. 652;;It used 10% more space. 205,000 for (x+y+z)^20*(x+y+z)^10 for ptimes. 653;;note on the st-rat form it only takes 510 msec. for (x+y+z)^10 and 654;;2.00 sec for (x+y+z)^20 as opposed to 3 sec and 15 sec resp with 655;;;50000 656; afp-square-free-factorization pfactor 657; poly (x+y+z)^10 cre .51 7,887 3.0 17000 658; poly (x+y+z)^10 general 1.8 28,00 4.5 68,000 659; poly (x+y+z)^20 cre 2.0 28,000 15. 252,000 660; poly (x+y+z)^20 general 7.2 105,000 20.8 325,000 661; 662; 663;(compare-functions 664;(defun af-fake-times (f g) 665; (user:tim (progn (afp-times f g) nil))) 666;(defun reg-fake-times (f g) 667; (user:tim (progn (ptimes f g) nil))) 668;) 669 670 671;;essentially same speed as regualar ptimes or maybe 5 %faster 672;;does (afp-times (x+y+z)^10 times itself in 2530 msec. and 2640 for ptimes. 673 674 675(defun afp-times (f g) 676 "The two arguments are polynomials and the result is the product polynomial" 677 (cond ((numberp f) 678 (cond ((zerop f) 0) 679 (t (afp-pctimes g f )))) 680 ((numberp g) 681 (cond ((zerop g) 0) 682 (t(afp-pctimes f g)))) 683 ((eq (p-var f) (p-var g)) 684 (palgsimp (p-var f) (afp-terms-times (cdr f)(cdr g)) (alg f))) 685 ((pointergp (p-var f) (p-var g)) 686 (psimp (p-var f) (plain-term-operation (cdr f) g afp-times))) 687 (t(psimp (p-var g) (plain-term-operation (cdr g) f afp-times))))) 688 689 690;;;the following shuffles the terms together and is about the same speed as the 691;;;corresponding maxima function ptimes1 692;(defun afp-terms-times (terms-f terms-g &aux prev to-add repl new-deg one-product) 693; "assumes same main variable" 694; (loop while terms-f 695; do 696; (setq one-product 697; (plain-term-operation terms-g (term-cof terms-f) afp-times (term-deg terms-f))) 698; until one-product 699; do (setq terms-f (cddr terms-f))) 700; (cond ((null one-product) nil) 701; (t 702; (loop for (deg-f cof-f) on (cddr terms-f) by 'cddr 703; do 704; (loop for (deg cof) on terms-g by #'cddr ;;sue 705; ;;computes the coeff of x^deg+deg-g and adds it in to 706; ;;the terms of one-product. Prev stands for the terms beginning 707; ;;with where the previous one was added on, or 708; initially 709; (setq prev one-product) 710; do (setq new-deg (+ deg deg-f)) 711; (setq to-add (afp-times cof-f cof)) 712; (cond ((pzerop to-add)) 713; ;;you should only get here when not in an integral domain 714; ((>= new-deg(car prev)) 715; (cond ((> new-deg (car prev)) 716; (setq one-product (setq prev (cons new-deg (cons to-add prev))))) 717; (t 718; ;;claim this can't happen unless prev = one-product 719; ;;Since otherwise have had a non-trivial to-add already in the sue 720; ;;loop and this would have added something in degree new-deg, but back 721; ;;one step, and prev would have that new-deg as its first element. 722; ;;That new-deg must be bigger than the current new-deg. 723; (iassert (eq prev one-product)) 724; (cond ((pzerop (setq repl (afp-plus (term-cof prev) to-add))) 725; (setq prev (setq one-product (cddr prev)))) 726; (t(setf (second prev) repl)))))) 727; (t ;;each to-add term has lower new-deg than the previous (or to-add=0) 728; (loop for vvv on (cddr prev) by #'cddr 729; do 730; (cond ((> (term-deg vvv) new-deg) (setq prev vvv)) 731; ((eql (term-deg vvv) new-deg) 732; (cond ((pzerop (setq repl (afp-plus (term-cof vvv) to-add))) 733; (setf (cddr prev)(cddr vvv))) 734; (t (setf (term-cof vvv) repl) 735; (setq prev (cddr prev)))) 736; (return nil)) 737; (t (setf (cddr prev) (cons new-deg (cons to-add vvv))) 738; (setq prev (cddr prev)) 739; (return nil))) 740; finally (setf (cddr prev)(list new-deg to-add)))))) 741; finally (return one-product))))) 742 743 744(defun test-decrease (terms) 745 (loop for (deg cof) on terms by #'cddr 746 with d0 = 1000 747 when (not (> d0 deg)) do (fsignal 'bad-order) 748 else do (setq d0 deg))) 749 750(defun afp-pctimes (poly number) 751 "Its first argument must be a polynomial and second argument a number,~ 752 and it returns the product" 753 (cond ((atom poly)(ctimes poly number)) 754 (t (term-operation poly number afp-pctimes)))) 755 756(defmacro butlastn (n list) 757 "knocks off the last n items of a list" 758 `(setf ,list (cond ((< ,n (length ,list)) 759 (setf (cdr (lastn (1+ ,n) ,list)) nil) ,list) 760 (t nil)))) 761 762(defun test-times (f g &key empty) 763 (cond (modulus (iassert (equal (tim (afp-times f g)) 764 (remove-zero-coefficients (tim (ptimes f g)))))) 765 (empty (tim (progn 766 (afp-times f g) 767 nil)) 768 (tim (progn 769 (ptimes f g) 770 nil))) 771 (t (iassert (equal (tim (afp-times f g)) 772 (tim (ptimes f g))))))) 773 774 775(defun remove-zero-coefficients (poly) 776 (cond ((numberp poly)poly) 777 (t (loop with v = poly 778 while (cdr v) 779 do (setf (third v) (remove-zero-coefficients (third v))) 780 when (pzerop (third v)) 781 do(setf (cdr v) (cdddr v)) 782 else 783 do (setq v (cddr v)) 784 finally (return (cond ((null (cdr poly)) 0) 785 ((zerop (second poly))(third poly)) 786 (t poly))))))) 787 788 789(defmacro with-area-used (&rest body) 790 `(progn 791 (prog1 792 (progn ,@body)))) 793 794(defun recursive-ideal-gcd1 (f g ) 795 "assumes that f and g are polynomials of one variable and that modulus is non trivial 796 and that deg f >= deg g gcd = a*f +b*g , and deg a < deg g, deg b < deg f" 797 (cond ((numberp g)(setq g (cmod g)) 798 (cond ((zerop g)(values f 1 0)) 799 (t (values 1 0 (crecip g))))) 800 ((not (eql (p-var g) (p-var f)))(merror 'not-function-of-one-variable)) 801 (t(multiple-value-bind (quot zl-rem creqd) 802 (afp-pseudo-quotient f g) 803 creqd ;;ignore 804 (multiple-value-bind (gcd a b) 805 (recursive-ideal-gcd1 g zl-rem) 806 (values gcd b (pdifference a (ptimes quot b)))))))) 807 808 809 810(defun recursive-ideal-gcd (f g &aux rev? fact) 811 "assumes that f and g are polynomials of one variable and that modulus is non trivial 812 It returns (gcd a b) such that gcd = a*f +b*g , and deg a < deg g, deg b < deg f where 813 gcd is the gcd of f and g in the polynomial ring modulo modulus." 814 (cond ((null modulus) (merror "polynomials over the integers are not a PID"))) 815 (assume-pointerg-or-equal-p f g rev?) 816 (cond ((numberp f)(values 1 (crecip f) 0)) 817 (t (multiple-value-bind (gcd a b) 818 (recursive-ideal-gcd1 f g) 819; (iassert (zerop (pdifference gcd (pplus (ptimes f a) (ptimes g b))))) 820; (iassert (numberp (afp-quotient gcd (pgcd f g)))) 821 (cond ((numberp gcd) 822 (cond ((not(equal gcd 1))(setq fact (crecip gcd))))) 823 (t (cond ((not(equal (p-cof gcd) 1))(setq fact (crecip (p-cof gcd))))))) 824 (cond (fact (setq gcd (ptimes fact gcd)) 825 (setq a (ptimes a fact)) 826 (setq b (ptimes b fact)))) 827 (cond (rev? (values gcd b a)) 828 (t (values gcd a b))))))) 829 830 831 832(defvar $e_poles nil) 833 834;;;do the following at macsyma level to specify which poles you want: 835;;; e_poles:[e^-2,e^-1]$ 836 837(defun $list_poles_of_sum (sum) 838 (cond ((numberp sum) '((mlist simp) 0 0)) 839 (t 840 (check-arg sum (and (consp sum)(eql (caar sum) 'mplus)) "macsyma sum") 841 (let ((pol $e_poles)) 842 (cons (car pol) (loop for u in (cdr pol) 843 collecting 844 (loop for v in (cdr sum) 845 when (equal u v) 846 collecting 1 into cof 847 else 848 when (and (listp v) (member u (cdr v) :test 'equal)) 849 collecting 850 (cond ((eql (length v) 3) 851 (loop for vv in (cdr v ) 852 when (not (equal vv u)) 853 do (return vv))) 854 (t (loop for vv in v 855 when (not (equal vv u)) 856 collecting vv))) 857 into cof 858 finally (return (cond ((null cof) 0) 859 ((eql (length cof) 1) (car cof)) 860 (t (apply 'add* cof))))))))))) 861 862(defun constant-psublis (alist polynomial) 863 (setq alist (sort alist #'pointergp :key #'car)) 864 (constant-psublis1 alist polynomial)) 865 866(defun constant-psublis1 (alist polynomial) 867 (cond ((numberp polynomial) polynomial) 868 ((null alist) polynomial) 869 (t 870 (block sue 871 (prog 872 (main-var tem) 873 (setq main-var (p-var polynomial)) 874 reduce-subs 875 (cond ((and alist (pointergp (caar alist) main-var)) 876 (setq alist (cdr alist)) (go reduce-subs))) 877 (cond ((null alist) (return polynomial)) 878 ((eql (caar alist) main-var) 879 (loop for (deg cof) on (p-next-term (p-terms polynomial)) 880 by #'p-next-term 881 with repl = (cdr (car alist)) 882 with answ = (afp-pctimes 883 (constant-psublis1 884 (cdr alist) (p-cof polynomial)) 885 (cexpt repl (p-deg polynomial))) 886 do (setq answ 887 (pplus answ 888 (cond ((zerop deg) 889 (constant-psublis1 (cdr alist) cof)) 890 (t (afp-pctimes (constant-psublis1 891 (cdr alist) cof) 892 (cexpt repl deg)))))) 893 finally (return-from sue answ))) 894 (t (return 895 (psimp main-var 896 (loop for (deg cof) 897 on (p-terms polynomial) by #'p-next-term 898 unless (pzerop (setq tem (constant-psublis1 alist cof))) 899 nconc (list deg tem))))))))))) 900 901;;afp-pcsubsty used 1/2 space and was twice as fast as pcsubsty on substituting for one variable y 902;;in (x+y+z)^10 (33.5 msec. and 70 msec respect). 903; 904;(compare-functions 905;(defun afp-pcsubsty (vals vars pol) 906; (constant-psublis (subs-for-psublis vars vals) pol)) 907;(defun reg-pcsubsty (vals vars pol) 908; (pcsubsty vals vars pol))) 909 910(defun my-evenp (n &aux nn) 911 (setq nn (ash n -1)) 912 (and (equal (ash nn 1) n) nn)) 913 914;;On symbolics the regular gcd is only 40% of the speed of fast-gcd. 915;;and I compared the values on several hundred random 916;;m n and it was correct. 917;;on (test 896745600000 7890012 1000) of 1000 repeats 918;;the fast-gcd was .725 sec and the regualar gcd was 5.6 sec. using 0 and 23000 words resp. 919;;On Explorer: Release 1.0)the fast-gcd was 1.9 sec and the regualar gcd was .102 sec. using 0 and 0 words resp. 920;(defun test (m n rep &aux a b) 921; (tim (loop for i below rep 922; do (setq a (fast-gcd m n)))) 923; (tim (loop for i below rep 924; do (setq b (\\\\ m n)))) 925; (assert (equal a b))) 926;;for testing two gcd's give same results. 927;(defun test ( rep &aux m n a b) 928; (loop for i below rep 929; do (setq m (* (random (expt2 32))(setq n (random (^ 2 32))))) 930; when (oddp i) do(setq n (- n)) 931; do 932; (setq n (* n (random (^ 2 32)))) 933; do (setq a (gcd m n)) 934; do (setq b (\\\\ m n)) 935; (show (list m n a)) 936; (assert (equal a b)))) 937 938 939(defun fast-gcd (m n) 940 (setq m (abs m) n (abs n)) 941 (cond ((< n m)nil) 942 (t (rotatef m n))) 943 (cond ((zerop n) m) 944 ((fixnump n) 945 (setq m (mod m n)) 946 (cond ((zerop m) n) 947 (t (bin-gcd m n)))) 948 (t (gcd n (mod m n))))) 949 950 951(defun bin-gcd (u v &aux (k 0)u2 v2 t2 tt) 952 (loop 953 do (setq u2 (ash u -1)) 954 when (not (eql (ash u2 1) u)) 955 do (return k) 956 do (setq v2 (ash v -1)) 957 when (not (eql (ash v2 1) v)) 958 do (return k) 959 do (setq u u2 v v2 k (1+ k))) 960 (prog () 961 b2 962 (cond ((oddp u) (setq tt (- v))) 963 (t (setq tt (ash u -1)))) 964 b3b4 965 (loop do (setq t2 (ash tt -1)) 966 when (eql (ash t2 1) tt) 967 do (setq tt t2) 968 else do (return nil)) 969 (cond ((> tt 0) (setq u tt)) 970 (t (setq v (- tt)))) 971 (setq tt (- u v)) 972 (cond ((zerop tt)(return (ash u k))) 973 (t (go b3b4))))) 974 975(defun poly-length (poly) 976 (cond ((numberp poly ) 0) 977 (t (length (cdr poly))))) 978 979(defun poly-to-row (poly &optional row &aux leng) 980 (cond (row 981 (cond ((< (array-total-size row) (length poly)) 982 (setq row (adjust-array row (+ 10 (poly-length poly)) :fill-pointer (fill-pointer row)))))) 983 (t 984 (setq row (make-array (+ 10 (setq leng (poly-length poly))) :fill-pointer 0 :adjustable t)))) 985 (cond ((numberp poly) 986 (vector-push 0 row) 987 (vector-push poly row)) 988 (t 989 (loop for u in (cdr poly) do 990 (vector-push u row)))) 991 row) 992 993(defun row-to-terms (row) 994 (listarray (sort-grouped-array row 2 '>))) 995 996;;seems afp-square is roughly twice as fast on univariate. 997;(user:compare-recursive-functions 998 999(defun afp-square (poly) 1000 (cond ((numberp poly)(ctimes poly poly)) 1001 (t (palgsimp (p-var poly) (afp-terms-square (p-terms poly)) (alg poly))))) 1002 1003;;comparison for squaring (x+y+z)^10 and (x+y+z)^4 (x+1)^10 1004;;;done in temp area: 1005; afp-square (time space) pexpt (time space) 1006;(x+y+z)^10 1450 ms. 25,305 1780ms 32,270 1007;(x+y+z)^4 80 ms. 1,443 111ms 2,099 1008;(x+1)^10 82 ms. 527 190ms 2,911 1009 1010;; (defun test-square (f &key empty) 1011;; (cond (modulus 1012;; (iassert (equal (tim (afp-square f)) (remove-zero-coefficients (tim (pexpt f 2 )))))) 1013;; (empty 1014;; (tim (progn 1015;; (afp-square f) 1016;; nil)) 1017;; (tim (progn 1018;; (pexpt f 2) 1019;; nil))) 1020;; (t 1021;; (iassert (equal (tim (afp-square f)) (tim (pexpt f 2))))))) 1022 1023;; (defun test-square (f &key empty) 1024;; (cond (modulus 1025;; (iassert (equal (tim (afp-square f)) (remove-zero-coefficients (tim (ptimes f f)))))) 1026;; (empty 1027;; (tim (progn 1028;; (afp-square f) 1029;; nil)) 1030;; (tim (progn 1031;; (pexpt f 2) 1032;; nil))) 1033;; (t 1034;; (iassert (equal (tim (afp-square f)) (tim (ptimes f f))))))) 1035 1036;;pexpt is not accurate for modulus=9 1037;;note example set al:(x+y+z)^4 in polynomial form. 1038;;then (let ((modulus 9)) (equal (ptimes al al) (pexpt al 2))) ==> nil 1039;;but afp-square does work. 1040 1041(defun afp-terms-square (p-terms) 1042 (prog (lead orig-p2-terms p2-terms prod-exp prod-cof tail answ) 1043 begin 1044 (cond (p-terms (setq answ (afp-square (term-cof p-terms))) 1045 (setq orig-p2-terms(setq p2-terms (plain-term-operation (cddr p-terms) 2 afp-pctimes))) 1046 (setq answ (plain-term-operation p2-terms (term-cof p-terms) afp-times (term-deg p-terms))) 1047 (setq lead (afp-square (term-cof p-terms))) 1048 (cond ((pzerop lead) nil) 1049 (t (setq answ (cons (* 2 (term-deg p-terms)) 1050 (cons lead 1051 answ))))) 1052 (cond ((null answ) 1053 (setq p-terms (cddr p-terms)) 1054 (setq p2-terms (cddr p2-terms)) 1055 (go begin)) 1056 (t (go second-leading-square)))) 1057 (t (return nil))) 1058 tail-certain ;;add in term to tail 1059 (cond ((and (cdr tail)(> (second tail) prod-exp)) 1060 (setq tail (cddr tail)) (go tail-certain))) 1061 (cond ((or (null (cdr tail))(< (second tail) prod-exp)) 1062 (setf (cdr tail) (cons prod-exp (cons prod-cof (cdr tail)))) 1063 (setq tail (cddr tail)) (go next-double-product))) 1064 (cond ((pzerop (setq prod-cof (pplus (third tail) prod-cof))) 1065 (setf (cdr tail) (cdddr tail))) 1066 (t (setf (third tail) prod-cof) (setq tail (cddr tail)))) 1067 next-double-product 1068 (setq p2-terms (cddr p2-terms)) 1069 (cond ((null p2-terms)(go next-leading-square))) 1070 (setq prod-cof (afp-times (term-cof p-terms) (term-cof p2-terms))) 1071 (cond ((pzerop prod-cof) (go next-double-product))) 1072 (setq prod-exp (+ (term-deg p-terms) (term-deg p2-terms))) 1073 (go tail-certain) 1074 next-leading-square 1075 (setq orig-p2-terms (cddr orig-p2-terms)) 1076 second-leading-square 1077 (setq p-terms (cddr p-terms)) 1078 (cond ((null p-terms)(return answ))) 1079 (setq prod-cof (afp-square (term-cof p-terms))) 1080 (setq p2-terms orig-p2-terms) 1081 (setq tail (cdr answ)) 1082 (cond ((pzerop prod-cof) (go next-double-product))) 1083 (setq prod-exp (* 2 (term-deg p-terms))) 1084 (go tail-certain))) 1085 1086 1087(defmacro def-test (f1 f2) 1088 `(defun ,(intern (format nil "~A-~A" '#:test f1)) (&rest rest-args) 1089 (let (empty (*print-level* 2)(*print-length* 3) ansa ansb) 1090 1091 1092 (cond ((member ':empty rest-args :test #'eq) 1093 (setq rest-args (subseq rest-args 0 (- (length rest-args) 1094 (length (member ':empty rest-args :test #'eq))))) 1095 (setq empty t))) 1096 (format t "~%For functions ~A and ~A respectively, ~%with argument list being ~A" ',f1 ',f2 rest-args) 1097 (cond (empty (format t "~%All computations done in a temporary area:"))) 1098 (cond ((and (null empty)modulus) 1099 (progn 1100 (iassert (equal (setq ansa(tim (apply ',f1 rest-args))) 1101 (setq ansb (remove-zero-coefficients (tim (apply ',f2 rest-args )))))))) 1102 (empty (tim (progn 1103 (apply ',f1 rest-args ) 1104 nil)) 1105 (tim (progn 1106 (apply ',f2 rest-args) 1107 nil))) 1108 (t (progn 1109 (iassert (equal (setq ansa (tim (apply ',f1 rest-args ))) 1110 (setq ansb(tim (apply ',f2 rest-args))))))))))) 1111 1112;;timings for afp-expt and pexpt respectively with 5 th power 1113;(x+y+z)^4 1114;For functions AFP-EXPT and PEXPT respectively, 1115;with argument list (POLY EXPONENT) being ((Z 4 1 ...) 5) 1116;All computations done in a temporary area: 1117;2223.207 msec. at priority 1 1118;Reclaiming 39,793 in polynomial space (total 4,176,070).14:20:53 1119;3029.429 msec. at priority 1 1120;Reclaiming 50,579 in polynomial space (total 4,226,649).14:20:56" 1121; 1122;(x+1)^10 1123;For functions AFP-EXPT and PEXPT respectively, 1124;with argument list (POLY EXPONENT) being ((X 10 1 ...) 10) 1125;All computations done in a temporary area: 1126;1972.551 msec. at priority 1 1127;Reclaiming 17,220 in polynomial space (total 3,825,741).14:19:00 1128;3312.806 msec. at priority 1 1129;Reclaiming 28,634 in polynomial space (total 3,854,375).14:19:03" 1130 1131;;note ( pexpt (x+y+z)^4 5) yields false result.. if modulus is 4. 1132;;I checked that afp-expt was ok, by comparing to the following 1133;;simple exponent function. 1134;(def-test afp-expt pexpt-simple) 1135;(defun pexpt-simple (u n) 1136; (loop for i from 1 to n 1137; with answ = 1 1138; do (setq answ (ptimes answ u)) 1139; finally (return answ))) 1140(defun afp-expt-modulo (poly exponent &rest work-modulo) 1141 (working-modulo work-modulo 1142 (afp-expt poly exponent))) 1143 1144(defun afp-expt (poly exponent) 1145 "Raises the polynomial POLY to EXPONENT. It never performs more 1146 than a squaring before simplifying, so that if modulus or tellrat are 1147 in effect, it will still be reasonable." 1148 (cond ((eql exponent 1) poly) 1149 ((eql exponent 0) poly) 1150 ((< exponent 0) (merror "Use positive exponents")) 1151 (t (cond ((numberp poly)(cexpt poly exponent)) 1152 ((or (cdddr poly) (alg poly)) 1153 ;;main case 1154 (loop for i from 0 1155 with 2^i-power-poly = poly 1156 with answer = 1 1157 do 1158 (cond ((oddp exponent) 1159 (cond ((eql answer 1) 1160 (setq answer 2^i-power-poly)) 1161 (t 1162 (setq answer (afp-times answer 2^i-power-poly)))))) 1163 1164 (setq exponent (ash exponent -1)) 1165 (cond ((zerop exponent)(return answer))) 1166 (setq 2^i-power-poly (afp-square 2^i-power-poly)))) 1167 (t ;;monomial of variable not tellrated 1168 (let ((pow (afp-expt (p-cof poly) exponent))) 1169 (cond ((pzerop pow) 0) 1170 (t 1171 (psimp (p-var poly) 1172 (list (* (p-deg poly) exponent) 1173 pow)))))))))) 1174(defmacro push-poly-number-deg-cof (rows poly-number deg cof &aux (row (make-symbol "row"))) 1175 `(let ((,row (aref ,rows ,deg))) 1176 (vector-push-extend ,poly-number ,row) 1177 (vector-push-extend ,cof ,row))) 1178 1179 ;;want to find sol'ns of v^p-v=0 (mod u) 1180(defun berlekamp-set-up-and-reduce-matrix ( u p &aux rows powers estimated-size sp) 1181 (setq rows (make-array (p-deg u) :fill-pointer (p-deg u))) 1182 (let ((modulus p)(tellratlist (list u))) 1183 (working-modulo (list u) 1184 (setq powers 1185 (loop for i from 0 1186 with pol = (afp-expt (list (p-var u) 1 1) p) 1187 with pow = 1 with belo = (1- (p-deg u)) 1188 collecting pow 1189 while (< i belo) 1190 when (eql pow 1) 1191 do (setq pow pol) 1192 else do (setq pow (afp-times pol pow))))) 1193 (loop for i below (max 5 (length powers)) 1194 for vv in powers 1195 summing (or (and (atom vv) 0) (length vv)) into count 1196 finally (setq estimated-size (+ 10 (quotient count 5)))) 1197 (loop for i below (fill-pointer rows) 1198 do (setf (aref rows i) (make-array estimated-size :fill-pointer 0 :adjustable t))) 1199 ;;putting the entries in the sparse matrix. Each polynomial is a column. 1200 (loop for vv in (cdr powers) 1201 for i from 1 1202 when (numberp vv) 1203 do 1204 (push-poly-number-deg-cof rows i 0 vv) 1205 (push-poly-number-deg-cof rows i i -1) 1206 else 1207 do 1208 (loop for (deg cof) on (p-terms vv) by #'p-next-term 1209 with subtracted-it 1210 when (eql i deg) 1211 do 1212 (setq cof (- cof 1)) 1213 (setq subtracted-it t) 1214 do 1215 (cond ((not (zerop cof)) 1216 (push-poly-number-deg-cof rows i deg cof))) 1217 finally (cond ((null subtracted-it) 1218 (push-poly-number-deg-cof rows i i -1)))) 1219 ) 1220 (setq sp (make-sparse-matrix :rows rows 1221 :type-of-entries p 1222 :columns-used-to-pivot (make-hash-table :test 'equal))) 1223 (sp-set-rows sp rows (loop for i from 1 below (p-deg u) collecting i)) 1224 (sp-reduce sp) 1225 sp)) 1226 1227(defun berlekamp-polynomial-solutions (sp polynomial &aux sp-sols) 1228 (sp-solve sp ) 1229 (setq sp-sols (sp-solutions sp)) 1230 (loop for i below (row-length (sp-rows sp-sols)) 1231 collecting (psimp (p-var polynomial)(listarray (sort-grouped-array (sp-row sp-sols i) 2 '>))))) 1232 1233 1234(defun berlekamp-get-factors-little-prime (u reduced-sparse-matrix prime &aux number-of-factors b-polys poly tem 1235 (factor-list (list u))) 1236 (setq number-of-factors (- (p-deg u) (sp-number-of-pivots reduced-sparse-matrix) )) 1237 (show number-of-factors) 1238 (sp-show-matrix reduced-sparse-matrix) 1239 (cond ((eql 1 number-of-factors) ;;irreducible 1240 factor-list) 1241 (t (setq b-polys (berlekamp-polynomial-solutions 1242 reduced-sparse-matrix u)) 1243 (loop named sue 1244 for v in b-polys 1245 with half-p = (ash prime -1) 1246 do 1247 (loop for j from (- half-p) to half-p 1248 do 1249 (loop for u-fact in factor-list 1250 do 1251 ;;make more efficient. 1252 (setq poly (pplus v j)) 1253 (setq tem (afp-gcd poly u-fact)) 1254 when (not (or (numberp tem) 1255 (eql (p-deg tem) (p-deg u-fact)))) 1256 do (setq tem (afp-make-monic tem)) 1257 (setq factor-list (cons tem 1258 (cons (afp-quotient u-fact 1259 tem) 1260 (delete u-fact factor-list :test #'equal)))) 1261 when (eql (length factor-list) number-of-factors) 1262 do (return-from sue factor-list))))))) 1263 1264(defun berlekamp-get-factors-big-prime (u reduced-sparse-matrix prime &aux number-of-factors b-polys 1265 (factor-list (list u))) 1266 (setq number-of-factors (- (p-deg u) (sp-number-of-pivots reduced-sparse-matrix) )) 1267 (show number-of-factors) 1268 1269 (cond ((eql 1 number-of-factors) ;;irreducible 1270 factor-list) 1271 (t (setq b-polys (berlekamp-polynomial-solutions 1272 reduced-sparse-matrix u)) 1273 (prog (half-p tem pow v) 1274 (setq half-p (ash prime -1)) 1275 added-factor 1276 (cond ((>= (length factor-list) number-of-factors)(return factor-list))) 1277 (setq v (loop for vv in (cdr b-polys) 1278 with answ = (afp-pctimes (car b-polys) (random (max 500 prime))) 1279 do (setq answ (pplus answ (afp-pctimes vv (random prime)))) 1280 finally (return answ))) 1281 (working-modulo 1282 (list u) 1283 (setq pow (pdifference (afp-expt v half-p) 1))) 1284 (setq factor-list 1285 (loop for w in factor-list 1286 do (setq tem (afp-make-monic (afp-gcd pow w))) 1287 when (not (or (numberp tem) 1288 (eql (p-deg tem) (p-deg w)))) 1289 collecting tem 1290 and 1291 collecting (afp-quotient w tem) 1292 else collecting w)) 1293 (go added-factor))))) 1294 1295(defun afp-make-monic (poly) 1296 (cond ((numberp poly) 1) 1297 ((not (or modulus (eql (abs (p-cof poly)) 1))) 1298 (merror "not finite modulus or unit coefficient")) 1299 (t (let ((inv (crecip (p-cof poly )))) 1300 (afp-pctimes poly inv))))) 1301 1302 1303 1304(defun afp-try-make-monic (poly) 1305 (cond ((numberp poly) 1) 1306 ((eql (p-cof poly) 1) poly) 1307 ((eql (p-cof poly) -1) 1308 (pminus poly)) 1309 (modulus 1310 (cond ((numberp (p-cof poly)) 1311 (afp-pctimes poly (crecip (p-cof poly)))) 1312 (t poly))) 1313 (t poly))) 1314 1315(defun afp-berlekamp-factor (pol) 1316 (afp-berlekamp-factor1 pol modulus :use-little (<= modulus 13))) 1317 1318(defun afp-berlekamp-factor1 (pol p &key use-little use-big &aux sp (modulus p) answ) 1319 (setq sp (berlekamp-set-up-and-reduce-matrix pol p) ) 1320 (cond ((and (null use-big)(or use-little (< p 13))) (setq answ (berlekamp-get-factors-little-prime pol sp p))) 1321 (t (setq answ (berlekamp-get-factors-big-prime pol sp p)))) 1322 (loop for v in answ 1323 with ans = 1 1324 do (setq ans (ptimes ans v)) 1325 finally (show ans) (iassert (equal (afp-mod pol) ans))) 1326 (show answ) 1327 (mapcar #'afp-mod answ)) 1328 1329(defvar *mod-p-factor* 'afp-berlekamp-factor) 1330;(defvar *mod-p-factor* 'afp-distinct-degree-factor) 1331 1332(defvar *sort-factors* t) 1333 1334(defun afp-factor (pol &key square-free-arg &aux factor-function facts vars lead answ) 1335 (cond (modulus (setq lead (p-cof pol)) 1336 (setq pol (afp-try-make-monic pol)))) 1337 (cond (square-free-arg (setq facts (list pol 1))) 1338 (t(setq facts (afp-square-free-factorization pol)))) 1339 (setq vars (list-variables pol)) 1340 (setq answ (cond (modulus 1341 (cond ((cdr vars) (fsignal 'not-univariate)) 1342 (t (setq factor-function *mod-p-factor* ) 1343 (setq answ 1344 (loop for (pol deg) on facts by #'cddr 1345 with const = 1 1346 when (consp pol) 1347 appending 1348 (loop for fa in (funcall factor-function pol) 1349 collecting fa 1350 collecting deg) into all 1351 else do (setq const (ctimes const pol)) 1352 (iassert (eql deg 1)) 1353 finally (setq const (* lead const)) 1354 (return(cond ((not (eql const 1)) 1355 (cons const (cons 1 all))) 1356 (t all)))))))) 1357 (t (fsignal 'not-yet-without-modulus)))) 1358 (cond (*sort-factors* (sort-grouped-list answ 2 'alphagreatp)) 1359 (t answ))) 1360 1361(defun afp-mod (pol) 1362 (afp-pctimes pol 1)) 1363 1364(defun get-factors (u use-for-gcd p &aux tem) 1365 (let ((modulus p)) 1366 (loop for w in use-for-gcd 1367 appending 1368 (loop for i below p 1369 when (not (numberp (setq tem (afp-gcd (pplus w i) u)))) 1370 do (show i) and 1371 collecting (monize tem))))) 1372 1373 1374(defun afp-distinct-degree-factor (u &optional ( prime modulus) &aux (v u )(ww (list (p-var u) 1 1)) (d 0) gd 1375 mon (modulus prime) answer) 1376 "U should be univariate and square free. Calculations are done modulo PRIME. It 1377 returns an alternating list of factors" 1378 (setq mon ww) 1379 (setq answer 1380 (loop do (cond ((numberp v) (return answ)) 1381 ((> (* 2 (1+ d)) (p-deg v))(return (cons v answ)))) 1382 (incf d) 1383 (working-modulo (list v) (setq ww(afp-expt ww prime))) 1384 (setq gd (afp-gcd (pdifference ww mon) v)) 1385 (cond ((and (consp gd) (> d (p-deg gd)))(fsignal 'big))) 1386 when (not (numberp gd)) 1387 do 1388 (setq v (afp-quotient v gd)) 1389 (and (consp v) (consp ww) (setq ww (palgsimp (p-var v) (cdr ww) (cdr v)))) 1390 when (not (numberp gd)) 1391 appending (one-degree-factors (afp-try-make-monic gd) d prime) into answ)) 1392 (iassert (eql (p-deg u) 1393 (loop for v in answer 1394 summing (p-deg v) ))) 1395 answer) 1396 1397#+debug 1398(defun tel (u n p &aux (modulus p)) 1399 (check-arg p (oddp p) "odd") 1400 (iassert (eql (p-deg u) 1401 (loop for v in answ 1402 summing (p-deg v) ))) 1403 1404 (eql (p-deg u) n)) 1405(defun one-degree-factors (u deg p &aux pow tt tem facts used-tt answ (modulus p)) 1406 "assumes u has its irreducible factors of degree DEG and works modulo P. U should 1407 be univariate and square free. The result is thus just a list of the factors (powers 1408 are unnecessary" 1409 (check-arg p (oddp p) "odd") 1410 (cond((eql (p-deg u) deg)(setq answ (list u ))) 1411 (t (setq pow (quotient (- (expt p deg) 1) 2)) 1412 (setq facts (list u)) 1413 (loop named sue for i from 1 1414 do 1415 (loop do (setq tt (generate-t-for-one-degree-factors u deg p i)) 1416 unless (or (numberp tt) (member tt used-tt :test #'equal)) 1417 do (push tt used-tt) (return tt)) 1418 (working-modulo (list u) 1419 (setq tt (afp-expt tt pow) )) 1420 (loop for v in facts 1421 do 1422 (setq tem (afp-gcd v (pplus tt 1))) 1423 (cond ((not (or (numberp tem) (eql (p-deg tem) (p-deg v)))) 1424 (cond ((eql (p-deg tem) deg)(push (afp-make-monic tem) answ)) 1425 (t (push tem facts))) 1426 (cond ((eql (p-deg (setq tem (afp-quotient v tem))) deg) 1427 (push (afp-make-monic tem) answ)) 1428; ((numberp tem)) 1429 (t (push tem facts))) 1430 (setq facts (delete v facts :test #'equal)))) 1431 when (null facts) do (return-from sue answ))))) 1432 (cond ((member 1 answ) (fsignal 'bad))) 1433 (iassert (eql (p-deg u) 1434 (loop for pol in answ 1435 when (and (consp pol)) 1436 summing (p-deg pol) ))) 1437 answ) 1438 1439(defun generate-t-for-one-degree-factors (u deg p i &aux tem) 1440 (cond ((< i 5) 1441 (list (p-var u) 1 1 0 (random p))) 1442 ((psimp (p-var u) (loop for j downfrom (setq tem (+ 1 (random (- (* 2 deg) 1)))) to 0 1443 ;;make semi sparse 1444 when (evenp (random 2)) 1445 collecting j and collecting (1+ (random (1- p)))))))) 1446 1447(defun ff (u &optional ( prime modulus) &aux fact1 facts) 1448 (let ((modulus prime)) 1449 (setq facts (afp-square-free-factorization u)) 1450 (sort-grouped-list 1451 (loop for (pol pow) on facts by #'cddr 1452 do (show pol pow) 1453 (setq fact1 (afp-distinct-degree-factor pol prime)) 1454 appending 1455 (loop for ww in fact1 1456 collecting ww collecting pow)) 1457 2 1458 'alphagreatp))) 1459 1460(defun nnpfactor (pol) 1461 (sort-grouped-list (npfactor pol) 2 'alphagreatp)) 1462 1463(def-test nnpfactor afp-factor) 1464 1465;(setq w2 (st-rat #$x^2+91*x+11$)) 1466;(setq v2 (st-rat #$x^2+19*x+11$)) 1467;(setq v1 (let ((modulus 7)) (afp-mod v2))) 1468;(setq w1 (let ((modulus 7)) (afp-mod w2))) 1469;(setq prod (ptimes v2 w2)) 1470 1471(defun hensel-lift (product v w prime up-to-size &aux a b gcd 1472 (facts (list v w))) 1473 "Lifts v and w which satisfy product=v*w mod(prime) to a list FACTS = (uu vv) 1474 satisfying product = uu*vv mod (up-to-size). Product, v, and w are assumed to 1475 have leading coefficient 1" 1476; (declare (values fact-pair power)) 1477 (let ((modulus prime)) (multiple-value-setq (gcd a b) 1478 (recursive-ideal-gcd v w)) 1479 (cond ((not (numberp gcd))(fsignal "must have gcd of factors a unit"))) 1480 (check-arg v (eql 1 (afp-mod (p-cof v))) "monic")) 1481 (loop with power = 1 1482 do (setq power (* power prime)) 1483 while (< power up-to-size) 1484 do 1485 (setq v (car facts) ) (setq w (second facts)) 1486 (setq facts (hensel-lift1 product v w power prime a b)) 1487 finally (return (values facts power)))) 1488 1489 1490(defun gen-afp-times (&rest lis) 1491 (setq lis (delete 1 (copy-list lis) :test #'equal)) 1492 (cond ((null lis) 1) 1493 ((null (cdr lis))(car lis)) 1494 (t (afp-times (car lis) (cond ((cddr lis) 1495 (apply 'gen-afp-times1 (cdr lis))) 1496 (t (second lis))))))) 1497 1498(defun gen-afp-times1 (&rest lis) 1499 (cond ((null lis) 1) 1500 ((null (cdr lis))(car lis)) 1501 (t (afp-times (car lis) (cond ((cddr lis) 1502 (apply 'gen-afp-times (cdr lis))) 1503 (t (second lis))))))) 1504 1505(defun smallest-power-bigger (p up-to) 1506 (loop with pow = p 1507 while (< pow up-to) 1508 do (setq pow (* p pow)) 1509 finally (return pow))) 1510 1511(defun hensel-lift-list (product factor-list prime up-to-size &aux lift ) 1512 (cond ((eql (length factor-list)1 ) (list product)) 1513 ((eql (length factor-list) 2) (hensel-lift product (first factor-list) 1514 (second factor-list) prime up-to-size)) 1515 (t (setq lift (hensel-lift product (first factor-list) (let((modulus prime))(apply 'gen-afp-times (cdr factor-list))) 1516 prime up-to-size)) 1517 (cons (car lift) (hensel-lift-list (second lift) (cdr factor-list) prime up-to-size))))) 1518 1519(defun hensel-lift1 (product ve we prev-modulus prime &optional a b &aux dif h kk quot zl-rem creqd new-modulus gcd) 1520 "lifts u=ve*we mod (p^e) to u=ve+1*we+1 mod (p^e+1) with ve=ve+1 and we=we+1 mod (p^e+1) 1521 and deg(ve+1)<=deg(ve) deg(we+1)<=deg(we) and returns the list of ve+1 and we+1" 1522; (declare (values (list ve+1 we+1))) 1523 (setq new-modulus (* prime prev-modulus)) 1524 (let ((modulus new-modulus)) (setq dif (pdifference product (ptimes ve we)))) 1525 (cond ((pzerop dif)(list ve we)) 1526 (t 1527 (let ((modulus prime) ) 1528 (cond ((null a) 1529 (multiple-value-setq (gcd a b) 1530 (recursive-ideal-gcd ve we))))) 1531 (let ((modulus new-modulus)) 1532 (setq h (ptimes b dif)) 1533 (setq kk (ptimes a dif)) 1534 (let ((modulus new-modulus)) 1535 (multiple-value-setq ( quot zl-rem creqd) 1536 (afp-pseudo-quotient h ve))) 1537 (setq h zl-rem) 1538 (setq kk (pplus kk (ptimes quot we))) 1539 (list (pplus ve h) (pplus we kk)))))) 1540 ;(let ((modulus (expt prime (1+ e)))) (values (pplus ve h) (pplus we kk)))) 1541 1542 1543(defun collect-number-factors (fact-list &aux answ) 1544 "Makes sure there are no factors with leading cof = -1 and collects all constants together 1545 and puts them first" 1546 (loop for (pol deg) on fact-list by #'cddr 1547 with constant = 1 1548 do 1549 (cond ((numberp pol) 1550 (setq constant (ctimes constant (cexpt pol deg)))) 1551 ((eql (p-cof pol) -1) 1552 (setq constant (ctimes constant (cexpt (p-cof pol) deg))) 1553 (setq pol (pminus pol)) 1554 (setq answ (nconc answ (list pol deg)))) 1555 (t(setq answ (nconc answ (list pol deg))))) 1556 finally (return 1557 (cond ((eql 1 constant) answ) 1558 (t (cons constant (cons 1 answ))))))) 1559 1560(defun integer-univariate-factor (poly &aux facts) 1561 "returns an alternating list of irreducible polynomials and their powers in the factorization over the integers 1562 there will be at most one factor which is a number and it will be to the first power. The other irreducible 1563 polynomials will be relatively prime" 1564 (setq facts (afp-square-free-factorization poly)) 1565 (show facts) 1566 (collect-number-factors (loop for (pol deg) on facts by #'cddr 1567 appending (loop for pol1 in (integer-univariate-factor1 pol) 1568 collecting pol1 collecting deg)))) 1569 1570;;(defvar *small-primes* '(3 5 7 11 13 17 19)) ; overrides *small-primes* in src/ifactor.lisp; not used 1571 1572(defun find-small-prime-so-square-free (poly &optional (variable (p-var poly)) (prime-start 13)&aux deriv cof the-gcd) 1573 "finds a small prime P so that the roots of poly with respect to variable VARIABLE are distinct and so 1574 that poly has the same degree reduced modulo P. It will continue for ever if the input is not square free over 1575 the integers." 1576 (setq deriv (pderivative poly variable)) 1577 (setq cof (pcoeff poly (list variable (pdegree poly variable) 1))) 1578 (cond ( (eql *mod-p-factor* 'afp-distinct-degree-factor) 1579 (setq prime-start (max prime-start 3)))) 1580 (loop for p from prime-start 1581 when (> (- p prime-start) 500) do (format t "~%I hope you have given me a square free polynomial!!") 1582 when (and (q-primep p) 1583 (let ((modulus p))(and (not (pzerop (afp-mod cof))) 1584 (or (numberp (setq the-gcd (afp-gcd (afp-mod poly) (afp-mod deriv)))) 1585 (not (member variable (list-variables the-gcd) :test #'eq)))))) 1586 do (return p))) 1587(defun integer-univariate-factor1 (pol &aux p ( up-to 1000) facts lifts mod) 1588 "Argument pol is square free" 1589 (setq p (find-small-prime-so-square-free pol )) 1590 (let ((modulus p)) (setq facts (afp-factor pol :square-free-arg t))) 1591 (setq mod (smallest-power-bigger p up-to)) 1592 (setq facts (loop for (pol deg) on facts by #'cddr 1593 when (not (eql deg 1)) do (fsignal "bad-degree-for-square-free") 1594 collecting pol)) 1595 (setq lifts (hensel-lift-list pol facts p up-to)) 1596 (correct-factors pol lifts p mod)) 1597 1598(defun sub-lists (leng list) 1599 "Returns list of ordered sublists of LIST of length LENG" 1600 (cond ((zerop leng) nil) 1601 ((eql 1 leng)(mapcar 'list list)) 1602 (t 1603 (loop for v on list 1604 appending (loop for w in (sub-lists (1- leng) (cdr v)) 1605 collecting (cons (car v) w)))))) 1606 1607(defun correct-factors (pol factors p mod &aux tried d answ quot prod) 1608 "Given the factors FACTORS of polynomial POL modulo MOD, we determine 1609 what the factors are over the integers. It is assumed that MOD is sufficiently 1610 large so that any real factors of POL lie in the range -MOD/2 to MOD/2. Also 1611 POL is assumed square free so that the factors are listed without their multiplicities" 1612 p 1613 (prog () 1614 (setq d 1) 1615 look-for-factors 1616 (cond ((> (* d 2) (length factors)) 1617 (push pol answ) (return answ))) 1618 (loop for v in (sub-lists d factors) 1619 when (member v tried :test #'equal) 1620 nconc nil 1621 else 1622 do (push v tried) 1623 (let ((modulus mod)) (setq prod (apply 'gen-ptimes v))) 1624; (setq prod (gen-afp-times (p-cof pol) prod)) 1625 (cond ((setq quot (afp-test-divide pol prod)) 1626 (setq pol quot) 1627 (push prod answ) 1628 (loop for vv in v 1629 do (setq factors (delete vv factors :count 1 :test #'equal))) 1630 (go look-for-factors))) 1631 finally (incf d) 1632 (go look-for-factors)))) 1633 1634 1635(defun q-primep (i &aux lis) 1636 (cond ((car (member i (setq lis '(2 3 5 7 11 13))))) 1637 ((evenp i) nil) 1638 ((loop for v in (cdr lis) 1639 when (zerop (mod i v)) 1640 do (return t)) 1641 nil) 1642 (t (primep i)))) 1643 1644 1645(defun test-integer-factor (pol &aux facts) 1646 (setq facts (integer-univariate-factor pol)) 1647 (iassert (equal pol (apply 'exponent-product facts))) 1648 facts) 1649 1650(defun subs-translate-sublis ( point &optional inv) 1651 (cond (inv (loop for v in point 1652 when (not (pzerop (cdr v))) 1653 collecting (cons (car v) 1654 (pplus (list (car v) 1 1) (cdr v))))) 1655 (t 1656 (loop for v in point 1657 when (not (pzerop (cdr v))) 1658 collecting (cons (car v) 1659 (pdifference (list (car v) 1 1) (cdr v))))))) 1660 1661;;change to a defremember. and then clear it when starting new problem. 1662;;or alternately don't really need to clear since the modulus is stored. 1663(defun express-x^i (f g k &optional (modulus modulus) &aux a b zl-REM quot gc mon case0 cre) 1664 "f and g should be univariate and the leading coefficient of g invertible modulo modulus. 1665 A list of two coefficients (a b) are returned such that x^i= a*f + b*g and deg(a)< deg (g)" 1666 (check-arg modulus (not (null modulus)) "non-trivial") 1667 (cond ((eql k 0) (multiple-value-setq (gc a b) 1668 (recursive-ideal-gcd f g)) 1669 (list a b)) 1670 (t (setq case0 (express-x^i f g 0)) 1671 (setq a (ptimes (car case0) (setq mon (list (car f) k 1)))) 1672 (setq b (ptimes (second case0) mon)) 1673 (multiple-value-setq (quot zl-rem cre) (afp-pseudo-quotient a g)) 1674 (setq a zl-rem) 1675 (setq b (pplus b 1676 (ptimes quot f))) 1677 ;;(iassert (equal mon (afp-plus (afp-times f a) (afp-times g b)))) 1678 (list a b)))) 1679 1680;(setq point (rerat '((z . 0) (y . 0) ))) 1681;(setq u (st-rat #$zzx^4+zzx^3*(3-z)+zzx^2*((y-3)*z-y^2-13) 1682; +zzx*((y^2+3*y+15)*z +6) 1683; +(-y^3-15*y)*z-2*y^2-30$)) 1684;(setq f (st-rat #$ (zzx^2+2)$)) 1685;(setq g (st-rat #$zzx^2+3*zzx-15$)) 1686;;the wr-lift worked for the above data 1687 1688;;speed hacks: The following has two areas which can be improved 1689;; 1. should not truncate the polynomial, should multiply with truncation. 1690;; 2. should make express-x^i a defremember. 1691(defun wr-lift (u fi gi up-to-k big-mod point 1692 &aux (modulus big-mod) (main-var (p-var u)) tem fi+1 gi+1 v f0 g0 varl w aw bw) 1693 (setq v (psublis (subs-translate-sublis point) 1 u)) 1694 (setq varl (delete (car u) (list-variables u) :test #'equal)) 1695 (setq f0 fi) 1696 (setq g0 gi) 1697 (loop for i from 1 to up-to-k 1698 do 1699 (setq w (pdifference (afp-times fi gi) u)) 1700 (setq w (afp-truncate-powers w varl i)) 1701 (mshow w) 1702 (cond ((pzerop w) 'nothing-to-do) 1703 (t 1704 (cond ((or (numberp w) (not (eql (p-var w) main-var))) 1705 (setq aw (afp-times w (first (setq tem (express-x^i f0 g0 0))))) 1706 (setq aw (afp-times w (second tem)))) 1707 (t 1708 (loop for (deg cof) on (cdr w) by #'cddr 1709 initially (setq aw 0 bw 0) 1710 do 1711 (setq aw (afp-plus aw (afp-times (first (setq tem (express-x^i f0 g0 deg))) cof))) 1712 (setq bw (afp-plus bw (afp-times (second tem) cof)))) 1713 1714 (setq fi+1 (pdifference fi bw)) 1715 (setq gi+1 (pdifference gi aw)) 1716 1717 (setq fi fi+1 gi gi+1) 1718 (mshow aw bw fi gi))))) 1719 finally 1720 (show (pdifference v (ptimes fi gi))) 1721 (return (list fi gi)))) 1722 1723(defun afp-truncate-powers (pol varl deg) 1724 (setq varl (sort varl 'pointergp)) 1725 ( afp-truncate-powers1 pol varl deg)) 1726 1727(defun afp-truncate-powers1 (pol varl above-deg &aux new-cof tem) 1728 (cond ((< above-deg 0) 0) 1729 ((numberp pol) pol) 1730 ((null varl) pol) 1731 ((member (p-var pol) varl :test #'eq) 1732 (psimp (p-var pol) 1733 (loop for (exp cof) on (cdr pol) by #'cddr 1734 do (setq tem (- above-deg exp)) 1735 when (< tem 0) 1736 nconc nil 1737 else 1738 do (setq new-cof (afp-truncate-powers1 cof (cdr varl) tem)) 1739 and when (not(pzerop new-cof)) 1740 nconc (list exp new-cof)))) 1741 (t(psimp (p-var pol) 1742 (loop for (exp cof) on (cdr pol) by #'cddr 1743 when(not (pzerop (setq new-cof (afp-truncate-powers1 cof varl above-deg)))) 1744 nconc (list exp new-cof)))))) 1745 1746 1747 1748(defun normalize-factor-list (factor-list &aux tot-deg) 1749 "checks factor-list and eliminates repeats. All integer factors are grouped at the beginning if any" 1750 (let (numerical-factor(facts (collect-number-factors factor-list))) 1751 (cond ((numberp (car facts)) 1752 (setq numerical-factor (subseq facts 0 2)) 1753 (setq facts (cddr facts)))) 1754 (nconc numerical-factor 1755 (loop for (pol deg) on facts by #'cddr 1756 for rest-facts on facts by #'cddr 1757 do (show facts deg) 1758 when deg 1759 do (setq tot-deg 1760 (+ deg 1761 (loop for v on (cddr rest-facts) by #'cddr 1762 when (and (second v) (equal pol (car v))) 1763 summing (second v) 1764 and do (setf (second v) nil)))) 1765 and 1766 collecting pol and collecting tot-deg)))) 1767 1768;; (defmacro while (test &body body) 1769;; `(loop (cond ((null ,test) (return))) 1770;; ,@ body)) 1771 1772(defun find-factor-list-given-irreducible-factors (pol fact-list &aux deg divisor answ final-answ) 1773 (while (setq divisor (car fact-list)) 1774 (setq deg 1) 1775 (while (setq answ (test-divide pol divisor)) 1776 (setq pol answ) 1777 (incf deg)) 1778 (setq final-answ (cons deg final-answ)) 1779 (setq final-answ (cons pol final-answ))) 1780 (unless (equal final-answ 1) 1781 (fsignal "bad factorization")) 1782 final-answ) 1783 1784;(defun afp-multi-factor (pol) 1785; (let ((content (afp-content pol))) 1786; (cond (content (normalize-factorlist 1787; (append (afp-multi-factor content) 1788; (afp-multi-factor (afp-quotient pol content ))))) 1789; (t (afp-))))) 1790 1791(defun try-multi-factor0 (pol &aux facs fac deg) 1792 "input may be non square free" 1793 (setq facs (afp-square-free-factorization pol)) 1794 (prog (result) 1795 next-factor 1796 (setq fac (car facs) deg (second facs) facs (cddr facs)) 1797 (prog 1798 ((newfacs (try-multi-factor1 fac))) 1799 next-factor 1800 (push (* (second newfacs) deg) result) 1801 (push (car newfacs) result) 1802 (setq newfacs (cddr newfacs)) 1803 (and newfacs (go next-factor))) 1804 (and facs (go next-factor)) 1805 (return result))) 1806 1807 1808(defun afp-degvector (pol vars &aux (result (make-list (length vars)))) 1809 (do ((v vars (cdr v)) 1810 (w result (cdr result))) 1811 ((null v) result) 1812 (setf (car w) (pdegree pol (car v))))) 1813 1814(defun afp-smallest-var (degvector list-variables &aux (best (expt 2 20)) best-variable ) 1815 (do ((v degvector (cdr v)) 1816 (w list-variables (cdr w))) 1817 ((null v) best-variable) 1818 (cond ((< (car v) best) 1819 (setq best-variable (car w)))))) 1820 1821;;description of eez algorithm with enhancements: 1822;; 1823;1) make U squarefree and content 1 1824;2) make main variable be smallest degree variable 1825;3) factor leading coefficient into F1,f2,..fk,fk+1 (fk+1 = content) 1826;4) find a point in x2,...,xn space such that 1827;there exist p1,..,pk+1, pi|fj iff i = j 1828;and 1829;such that we have the right number of factors of U mod (point, B) (B big bound). 1830;5) Match up pi, and the univariate factors. 1831;6) 1832 1833;(defun try-multi-factor1 (pol) 1834; "Pol is square free multivariate" 1835; (let* ((list-variables (list-variables pol)) 1836; (degvector (afp-degvector pol list-variables)) 1837; (best-variable (afp-smallest-var degvector list-variables)) 1838; (leading-cof (pcoeff pol best-variable)) 1839; (point (delq best-variable (copy-list list-variables)))) 1840; (do ((v point (cdr v))) 1841; (null v) 1842; (setf (car v (cons (car v) (nth (random 5) *small-primes*))))) 1843; (let* ((lead-cof-at-point (psublis point 1 leading-cof))) 1844; (cond ((eql 0 1845