[Maxima-commits] [git] Maxima CAS branch, master, updated. branch-5_31-base-80-g9f2bf3d
Subject: [Maxima-commits] [git] Maxima CAS branch, master, updated. branch-5_31-base-80-g9f2bf3d
From: Raymond Toy
Date: Thu, 24 Oct 2013 21:33:27 -0700
Unfortunately, this doesn't help. Cmucl is still slow. :-(
I'll try to take a look this weekend and see what's causing this problem.
On Thu, Oct 24, 2013 at 11:59 AM, Rupert Swarbrick <
rswarbrick at users.sourceforge.net> wrote:
> This is an automated email from the git hooks/post-receive script. It was
> generated because a ref change was pushed to the repository containing
> the project "Maxima CAS".
>
> The branch, master has been updated
> via 9f2bf3dc73250c88152d346918053b5b16aab866 (commit)
> from f5e6d072e3a41bf217e58abfb2eb2dea8cffa3c4 (commit)
>
> Those revisions listed above that are new to this repository have
> not appeared on any other notification email; so we list those
> revisions in full, below.
>
> - Log -----------------------------------------------------------------
> commit 9f2bf3dc73250c88152d346918053b5b16aab866
> Author: Rupert Swarbrick <rswarbrick at gmail.com>
> Date: Thu Oct 24 19:47:49 2013 +0100
>
> Revert cleanup changes to rat package
>
> Starting with the one that nuked errjfflag, since it caused a massive
> slowdown on CMUCL.
>
> diff --git a/share/affine/sheafa.lisp b/share/affine/sheafa.lisp
> index 466918d..705ae3d 100644
> --- a/share/affine/sheafa.lisp
> +++ b/share/affine/sheafa.lisp
> @@ -2003,7 +2003,8 @@
> ;;where zi=qi(x1,x2, xn)/gp and xi=pi(z1,..,zn)/gq
>
> (defmfun my-testdivide (x y)
> - (ignore-rat-err (pquotient x y)))
> + (let ((errrjfflag t))
> + (catch 'raterr (pquotient x y))))
>
> (defun new-testdivide (f g &aux quot)
> (setq quot (ratreduce f g))
> diff --git a/src/maxima.system b/src/maxima.system
> index 0e554ce..a3e46ef 100644
> --- a/src/maxima.system
> +++ b/src/maxima.system
> @@ -514,8 +514,7 @@
> (:file "rat3c")
> (:file "rat3e")
> (:file "nrat4")
> - (:file "ratout"))
> - :depends-on (rat-macros other-macros))
> + (:file "ratout")))
> (:module maxima-language-compiler :source-pathname ""
> :components ((:file "transl")
> (:file "transs")
> diff --git a/src/nrat4.lisp b/src/nrat4.lisp
> index f35546d..3f5c40c 100644
> --- a/src/nrat4.lisp
> +++ b/src/nrat4.lisp
> @@ -15,6 +15,8 @@
> (declare-top (special $ratsimpexpons *exp *exp2 *radsubst *loglist
> $radsubstflag
> $radexpand $logsimp *v *var fr-factor radcanp
> ratsubvl))
>
> +(load-macsyma-macros rzmac ratmac)
> +
> (defmvar $radsubstflag nil
> "`radsubstflag' `t' makes `ratsubs' call `radcan' when it appears
> useful")
>
> diff --git a/src/rat3a.lisp b/src/rat3a.lisp
> index a8cb706..ef84c4f 100644
> --- a/src/rat3a.lisp
> +++ b/src/rat3a.lisp
> @@ -25,169 +25,94 @@
> ;;slow it down on lispm. We also eliminated the special
> ;;from ptimes2--wfs
>
> +(load-macsyma-macros ratmac)
> +
> ;; Global variables referenced throughout the rational function package.
>
> (defmvar modulus nil "Global switch for doing modular arithmetic")
>
> -;; CQUOTIENT
> -;;
> -;; Calculate the quotient of two coefficients, which should be numbers. If
> -;; MODULUS is non-nil, we try to take the reciprocal of A with respect to
> the
> -;; modulus (using CRECIP) and then multiply by B. Note that this fails if
> B
> -;; divides A as an integer, but B is not a unit in the ring of integers
> modulo
> -;; MODULUS. For example,
> -;;
> -;; (let ((modulus 20)) (cquotient 10 5)) => ERROR
> -;;
> -;; If MODULUS is nil, then we work over the ring of integers when A and B
> are
> -;; integers, and raise a RAT-ERROR if A is not divisible by B. If either
> A or B
> -;; is a float then the division is done in floating point. Floats can get
> as far
> -;; as the rational function code if $KEEPFLOAT is true.
> +(defmacro bctimes (&rest l)
> + `(rem (* , at l) modulus))
> +
> +;; coefficient quotient a / b
> +;; a and b may be integers (possibly with modulus) or floats if
> keepfloat=true
> (defun cquotient (a b)
> (cond ((equal a 0) 0)
> ((null modulus)
> - (if (or (floatp a) (floatp b)
> - (zerop (rem a b)))
> - (/ a b)
> - (rat-error "quotient is not exact")))
> + (cond ((equal 0 (cremainder a b)) (/ a b))
> + (t (rat-error "quotient is not exact"))))
> (t (ctimes a (crecip b)))))
>
> -;; ALG
> -;;
> -;; Get any value stored on the tellrat property of (car l). Returns NIL
> if L
> -;; turns out not to be a list or if $ALGEBRAIC is false.
> (defun alg (l)
> - (unless (atom l) (algv (car l))))
> + (and $algebraic (not (atom l)) (get (car l) 'tellrat)))
>
> -;; PACOEFP
> -;;
> -;; Return T if either X is a bare coefficient or X is a polynomial whose
> main
> -;; variable has a declared value as an algebraic integer. Otherwise
> return NIL.
> (defun pacoefp (x)
> - (and (or (pcoefp x) (alg x))
> - T))
> + (or (pcoefp x) (alg x)))
>
> -;; LEADTERM
> -;;
> -;; Return the leading term of POLY as a polynomial itself.
> (defun leadterm (poly)
> - (if (pcoefp poly)
> - poly
> - (make-poly (p-var poly) (p-le poly) (p-lc poly))))
> -
> -;; CBEXPT
> -;;
> -;; Raise an number to a positive integral power. P should be a number (and
> -;; really only makes sense if it is an integer). N should be a
> non-negative
> -;; integer.
> + (cond ((pcoefp poly) poly)
> + (t (make-poly (p-var poly) (p-le poly) (p-lc poly)))))
> +
> +(defun cremainder (a b)
> + (cond ((or modulus (floatp a) (floatp b)) 0)
> + ((rem a b))))
> +
> (defun cbexpt (p n)
> (do ((n (ash n -1) (ash n -1))
> (s (if (oddp n) p 1)))
> ((zerop n) s)
> - (setq p (rem (* p p) modulus))
> - (when (oddp n) (setq s (rem (* s p) modulus)))))
> + (setq p (bctimes p p))
> + (and (oddp n) (setq s (bctimes s p)))))
> +
>
> ;; Coefficient Arithmetic -- coefficients are assumed to be something
> ;; that is NUMBERP in lisp. If MODULUS is non-NIL, then all coefficients
> ;; are assumed to be less than its value. Some functions use LOGAND
> ;; when MODULUS = 2. This will not work with bignum coefficients.
>
> -;; CRECIP
> -;;
> -;; Takes the inverse of an integer N mod MODULUS. If there is a modulus
> then the
> -;; result is constrained to lie in (-modulus/2, modulus/2]
> -;;
> -;; This just uses the extended Euclidean algorithm: once you have found
> a,b such
> -;; that a*n + b*modulus = 1 then a must be the reciprocal you're after.
> -;;
> -;; When MODULUS is greater than 2^15, we use exactly the same algorithm in
> -;; CRECIP-GENERAL, but it can't use fixnum arithmetic. Note: There's no
> -;; particular reason to target 32 bits except that trying to work out the
> right
> -;; types on the fly looks complicated and this lisp compiler, at least,
> uses 32
> -;; bit words. Since we have to take a product, we constrain the types to
> 16 bit
> -;; numbers.
> -(defun crecip (n)
> - ;; Punt on anything complicated
> - (unless (and modulus (typep modulus '(unsigned-byte 15)))
> - (return-from crecip (crecip-general n)))
> -
> - ;; And make sure that -MODULUS < N < MODULUS
> - (unless (<= (- modulus) n modulus)
> - (merror "N is out of range [-MODULUS, MODULUS] in crecip."))
> -
> - ;; N in [0, MODULUS]
> - (when (minusp n) (setf n (+ n modulus)))
> -
> - ;; The mod-copy parameter stores a copy of MODULUS on the stack, which
> is
> - ;; useful because the lisp implementation doesn't know that the special
> - ;; variable MODULUS is still an (unsigned-byte 15) when we get to the
> end
> - ;; (since it can't tell that our function calls don't change it behind
> our
> - ;; backs, I guess)
> - (let ((mod modulus) (remainder n) (a 1) (b 0)
> - (mod-copy modulus))
> - ;; On SBCL in 2013 at least, the compiler doesn't spot that MOD and
> - ;; REMAINDER are unsigned and bounded above by MODULUS, a 16-bit
> integer. So
> - ;; we have to tell it. Also, the lisp implementation can't really be
> - ;; expected to know that Bezout coefficients are bounded by the
> modulus and
> - ;; remainder, so we have to tell it that too.
> - (declare (type (unsigned-byte 15) mod mod-copy remainder)
> - (type (signed-byte 16) a b))
> -
> - (loop
> - until (= remainder 1)
> -
> - when (zerop remainder) do
> - (merror (intl:gettext "CRECIP: attempted inverse of zero (mod
> ~M)")
> - mod)
> - doing
> - (multiple-value-bind (quot rem)
> - (truncate mod remainder)
> - (setf mod remainder
> - remainder rem)
> - (psetf a (- b (* a quot))
> - b a))
> -
> - finally
> - ;; Since this isn't some general purpose Euclidean algorithm, but
> - ;; instead is trying to find a modulo inverse, we need to ensure
> that
> - ;; the Bezout coefficient we found (called A) is actually in [0,
> - ;; MODULUS).
> - ;;
> - ;; The general code calls CMOD here, but that doesn't know about
> the
> - ;; types of A and MODULUS, so we do it by hand, special-casing
> the easy
> - ;; case of modulus=2.
> - (return
> - (if (= mod-copy 2)
> - (logand a 1)
> - (let ((nn (mod a mod-copy)))
> - ;; nn here is in [0, modulus)
> - (if (<= (* 2 nn) mod-copy)
> - nn
> - (- nn mod-copy))))))))
> -
> -;; CRECIP-GENERAL
> -;;
> -;; The general algorithm for CRECIP, valid when the modulus is any
> integer. See
> -;; CRECIP for more details.
> -(defun crecip-general (n)
> - ;; We assume that |n| < modulus, so n+modulus is always positive
> - (let ((mod modulus)
> - (remainder (if (minusp n) (+ n modulus) n))
> - (a 1) (b 0))
> - (loop
> - until (= remainder 1)
> -
> - when (zerop remainder) do
> - (merror (intl:gettext "CRECIP: attempted inverse of zero (mod
> ~M)")
> - mod)
> - doing
> - (let ((quotient (truncate mod remainder)))
> - (psetf mod remainder
> - remainder (- mod (* quotient remainder)))
> - (psetf a (- b (* a quotient))
> - b a))
> -
> - finally (return (cmod a)))))
> +;; Takes the inverse of an integer N mod P. Solve N*X + P*Y = 1
> +;; I suspect that N is guaranteed to be less than P, since in the case
> +;; where P is a fixnum, N is also assumed to be one.
> +
> +(defmfun crecip (n)
> + (cond ((bignump modulus) ;; Have to use bignum arithmetic if
> modulus is a bignum
> + (prog (a1 a2 y1 y2 q (big-n n))
> + (if (minusp big-n) (setq big-n (+ big-n modulus)))
> + (setq a1 modulus a2 big-n)
> + (setq y1 0 y2 1)
> + (go step3)
> + step2 (setq q (truncate a1 a2))
> + (psetq a1 a2 a2 (- a1 (* a2 q)))
> + (psetq y1 y2 y2 (- y1 (* y2 q)))
> + step3 (cond ((zerop a2) (merror (intl:gettext "CRECIP:
> attempted inverse of zero (mod ~M)") modulus))
> + ((not (equal a2 1)) (go step2)))
> + (return (cmod y2))))
> + ;; Here we can use fixnum arithmetic
> + (t (prog ((a1 0) (a2 0) (y1 0) (y2 0) (q 0) (nn 0) (pp 0))
> + (declare (fixnum a1 a2 y1 y2 q nn pp))
> + (setq nn n pp modulus)
> + (cond ((minusp nn) (setq nn (+ nn pp))))
> + (setq a1 pp a2 nn)
> + (setq y1 0 y2 1)
> + (go step3)
> + step2 (setq q (truncate a1 a2))
> + (psetq a1 a2 a2 (rem a1 a2))
> + (psetq y1 y2 y2 (- y1 (* y2 q)))
> + step3 (cond ((= a2 0) (merror (intl:gettext "CRECIP:
> attempted inverse of zero (mod ~M)") modulus))
> + ((not (= a2 1)) (go step2)))
> + ;; Is there any reason why this can't be (RETURN (CMOD Y2))
> ? -cwh
> + (return (cmod y2)
> + ;; (COND ((= PP 2) (LOGAND 1 Y2))
> + ;; (T (LET ((NN (rem Y2
> PP)))
> + ;; (DECLARE (FIXNUM
> NN))
> + ;; (COND ((MINUSP NN)
> + ;; (AND (< NN
> (- (ASH PP -1)))
> + ;; (SETQ
> NN (+ NN PP))))
> + ;; ((> NN (ASH
> PP -1))
> + ;; (SETQ NN (-
> NN PP))))
> + ;; NN)))
> + )
> + ))))
>
> (defun cexpt (n e)
> (cond ((null modulus) (expt n e))
> diff --git a/src/rat3b.lisp b/src/rat3b.lisp
> index ca7c5e0..178542d 100644
> --- a/src/rat3b.lisp
> +++ b/src/rat3b.lisp
> @@ -17,6 +17,8 @@
>
> (declare-top (special $algebraic $ratfac $keepfloat $float))
>
> +(load-macsyma-macros ratmac)
> +
> (defmvar $ratwtlvl nil)
> (defmvar $ratalgdenom t) ;If T then denominator is rationalized.
>
> diff --git a/src/rat3c.lisp b/src/rat3c.lisp
> index 2d48c44..a92e2d5 100644
> --- a/src/rat3c.lisp
> +++ b/src/rat3c.lisp
> @@ -15,6 +15,8 @@
> ;; THIS IS THE NEW RATIONAL FUNCTION PACKAGE PART 3.
> ;; IT INCLUDES THE GCD ROUTINES AND THEIR SUPPORTING FUNCTIONS
>
> +(load-macsyma-macros ratmac)
> +
> (declare-top (special $float $keepfloat $algebraic $ratfac genvar))
>
> ;; List of GCD algorithms. Default one is first.
> @@ -283,12 +285,7 @@
> (go a)))
>
> (defun prem (p q)
> - (cond ((pcoefp p)
> - (if (pcoefp q)
> - (if (or modulus (floatp p) (floatp q))
> - 0
> - (rem p q))
> - p))
> + (cond ((pcoefp p) (if (pcoefp q) (cremainder p q) p))
> ((pcoefp q) (pzero))
> (t (psimp (p-var p) (pgcd1 (p-terms p) (p-terms q))))))
>
> diff --git a/src/rat3d.lisp b/src/rat3d.lisp
> index d19151f..58cf5c2 100644
> --- a/src/rat3d.lisp
> +++ b/src/rat3d.lisp
> @@ -12,6 +12,9 @@
>
> (macsyma-module rat3d)
>
> +(load-macsyma-macros ratmac)
> +
> +
> ;; THIS IS THE NEW RATIONAL FUNCTION PACKAGE PART 4.
> ;; IT INCLUDES THE POLYNOMIAL FACTORING ROUTINES.
>
> diff --git a/src/rat3e.lisp b/src/rat3e.lisp
> index 39c5439..a226c85 100644
> --- a/src/rat3e.lisp
> +++ b/src/rat3e.lisp
> @@ -20,6 +20,8 @@
> vlist scanmapp radlist expsumsplit *ratsimp* mplc*
> $ratsimpexpons $expop $expon $negdistrib $gcd))
>
> +(load-macsyma-macros rzmac ratmac)
> +
> (defmvar genvar nil
> "List of gensyms used to point to kernels from within polynomials.
> The values cell and property lists of these symbols are used to
> diff --git a/src/ratout.lisp b/src/ratout.lisp
> index 1d87f62..84c9ab9 100644
> --- a/src/ratout.lisp
> +++ b/src/ratout.lisp
> @@ -18,6 +18,8 @@
> genvar *a* *alpha *var* *x* *p *max *var *res *chk
> *l $intfaclim
> $ratfac u* $ratwtlvl *ratweights $ratweights
> $keepfloat))
>
> +(load-macsyma-macros ratmac)
> +
> (declare-top (special $gcd xv bigf1 bigf2 nonlindeg $linhack
> $intfaclim bigf1tilde bigf2tilde
> gcd $factorflag *gcdl* last-good-prime))
>
> -----------------------------------------------------------------------
>
> Summary of changes:
> share/affine/sheafa.lisp | 3 +-
> src/maxima.system | 3 +-
> src/nrat4.lisp | 2 +
> src/rat3a.lisp | 203
> +++++++++++++++-------------------------------
> src/rat3b.lisp | 2 +
> src/rat3c.lisp | 9 +--
> src/rat3d.lisp | 3 +
> src/rat3e.lisp | 2 +
> src/ratout.lisp | 2 +
> 9 files changed, 81 insertions(+), 148 deletions(-)
>
>
> hooks/post-receive
> --
> Maxima CAS
>
>
> ------------------------------------------------------------------------------
> October Webinars: Code for Performance
> Free Intel webinars can help you accelerate application performance.
> Explore tips for MPI, OpenMP, advanced profiling, and more. Get the most
> from
> the latest Intel processors and coprocessors. See abstracts and register >
> http://pubads.g.doubleclick.net/gampad/clk?id=60135991&iu=/4140/ostg.clktrk
> _______________________________________________
> Maxima-commits mailing list
> Maxima-commits at lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/maxima-commits
>