182 lines
4.8 KiB
Common Lisp
182 lines
4.8 KiB
Common Lisp
;;; math.lisp -- various numerical operations
|
|
;;
|
|
;; Time-stamp: <2004-01-05 emarsden>
|
|
;;
|
|
;; some basic mathematical benchmarks
|
|
|
|
|
|
(in-package :cl-bench.math)
|
|
|
|
|
|
(defun factorial (n)
|
|
(declare (type integer n))
|
|
(if (zerop n) 1
|
|
(* n (factorial (1- n)))))
|
|
|
|
(defun run-factorial ()
|
|
(declare (inline factorial))
|
|
(factorial 500))
|
|
|
|
(defun fib (n)
|
|
(declare (type integer n))
|
|
(if (< n 2) 1 (+ (fib (- n 1)) (fib (- n 2)))))
|
|
|
|
(defun run-fib ()
|
|
(declare (inline fib))
|
|
(fib 25))
|
|
|
|
(defun fib-ratio (n)
|
|
(declare (type integer n))
|
|
(labels ((fr (n)
|
|
(if (= n 1) 1
|
|
(1+ (/ (fr (- n 1)))))))
|
|
(numerator (fr n))))
|
|
|
|
(defun run-fib-ratio ()
|
|
(declare (inline fib-ratio))
|
|
(fib-ratio 150))
|
|
|
|
|
|
;; The Ackermann function is the simplest example of a well-defined total
|
|
;; function which is computable but not primitive recursive, providing a
|
|
;; counterexample to the belief in the early 1900s that every computable
|
|
;; function was also primitive recursive (Dtzel 1991). It grows faster
|
|
;; than an exponential function, or even a multiple exponential function.
|
|
(defun ackermann (m n)
|
|
(declare (type integer m n))
|
|
(cond
|
|
((zerop m) (1+ n))
|
|
((zerop n) (ackermann (1- m) 1))
|
|
(t (ackermann (1- m) (ackermann m (1- n))))))
|
|
|
|
(defun run-ackermann ()
|
|
(ackermann 3 11))
|
|
|
|
|
|
;; calculate the "level" of a point in the Mandebrot Set, which is the
|
|
;; number of iterations taken to escape to "infinity" (points that
|
|
;; don't escape are included in the Mandelbrot Set). This version is
|
|
;; intended to test performance when programming in naïve math-style.
|
|
(defun mset-level/complex (c)
|
|
(declare (type complex c))
|
|
(loop :for z = #c(0 0) :then (+ (* z z) c)
|
|
:for iter :from 1 :to 300
|
|
:until (> (abs z) 4.0)
|
|
:finally (return iter)))
|
|
|
|
;; this version is intended to test lower-level performance-oriented
|
|
;; coding of the same function; hence the extra declarations and the
|
|
;; decoding of the operations on complex numbers.
|
|
(defun mset-level/dfloat (c1 c2)
|
|
(declare (type double-float c1 c2))
|
|
(let ((z1 0.0d0)
|
|
(z2 0.0d0)
|
|
(aux 0.0d0))
|
|
(declare (double-float z1 z2 aux))
|
|
(do ((iter 0 (1+ iter)))
|
|
((or (> (abs (+ (* z1 z1) (* z2 z2))) 4.0)
|
|
(> iter 300))
|
|
iter)
|
|
(setq aux z1
|
|
z1 (+ (* z1 z1) (- (* z2 z2)) c1)
|
|
z2 (+ (* 2.0d0 z2 aux) c2)))))
|
|
|
|
(defun run-mandelbrot/complex ()
|
|
(let ((n 100)
|
|
(sum 0))
|
|
(dotimes (i n)
|
|
(incf sum (mset-level/complex (complex 0.0001d0 (/ i n 0.25d0)))))))
|
|
|
|
(defun run-mandelbrot/dfloat ()
|
|
(let ((n 100)
|
|
(sum 0))
|
|
(dotimes (i n)
|
|
(incf sum (mset-level/dfloat 0.0001d0 (/ i n 0.25d0))))))
|
|
|
|
|
|
|
|
;; Common Lisp implementation of the multiple recursive random number
|
|
;; generator (MRG) of l'Ecuyer. Written by Raymond Toy.
|
|
|
|
(eval-when (:compile-toplevel :load-toplevel)
|
|
(defconstant +m1+ 4294967087d0)
|
|
(defconstant +m2+ 4294944443d0)
|
|
(defconstant +a12+ 1403580d0)
|
|
(defconstant +a13n+ 810728d0)
|
|
(defconstant +a21+ 527612d0)
|
|
(defconstant +a23n+ 1370589d0)
|
|
(defconstant +norm+ (/ (1+ +m1+))))
|
|
|
|
(declaim (inline mrg32k3a-comp-1 mrg32k3a-comp-2))
|
|
|
|
(defun mrg32k3a-comp-1 (state)
|
|
(declare (type (simple-array double-float (6)) state)
|
|
(optimize (speed 3) (safety 0)))
|
|
(let ((s10 (aref state 0))
|
|
(s11 (aref state 1)))
|
|
(declare (type (double-float 0d0 4294967086d0) s10 s11))
|
|
(let* ((p1 (- (* s11 +a12+) (* s10 +a13n+)))
|
|
(k (ftruncate (/ p1 +m1+)))
|
|
(p1b (- p1 (* k +m1+)))
|
|
(p1c (if (< p1b 0)
|
|
(+ p1b +m1+)
|
|
p1b)))
|
|
(shiftf (aref state 0)
|
|
(aref state 1)
|
|
(aref state 2)
|
|
p1c)
|
|
p1c)))
|
|
|
|
(defun mrg32k3a-comp-2 (state)
|
|
(declare (type (simple-array double-float (6)) state)
|
|
(optimize (speed 3) (safety 0)))
|
|
(let ((s20 (aref state 3))
|
|
(s22 (aref state 5)))
|
|
(declare (type (double-float 0d0 4294944442d0) s20 s22))
|
|
(let* ((p2 (- (* s22 +a21+) (* s20 +a23n+)))
|
|
(k (ftruncate (/ p2 +m2+)))
|
|
(p2b (- p2 (* k +m2+)))
|
|
(p2c (if (< p2b 0)
|
|
(+ p2b +m2+)
|
|
p2b)))
|
|
(shiftf (aref state 3)
|
|
(aref state 4)
|
|
(aref state 5)
|
|
p2c)
|
|
p2c)))
|
|
|
|
(declaim (inline mrg32k3a))
|
|
|
|
(defun mrg32k3a (state)
|
|
(declare (type (simple-array double-float (6)) state)
|
|
(optimize (speed 3) (safety 0)))
|
|
(let ((p1 (mrg32k3a-comp-1 state))
|
|
(p2 (mrg32k3a-comp-2 state)))
|
|
(if (<= p1 p2)
|
|
(* (+ (- p1 p2) +m1+) +norm+)
|
|
(* (- p1 p2) +norm+))))
|
|
|
|
(defun gen-mrg32 (n state)
|
|
(declare (fixnum n)
|
|
(optimize (speed 3) (safety 0)))
|
|
(let ((y 0d0))
|
|
(dotimes (k n)
|
|
(declare (fixnum k))
|
|
(setf y (mrg32k3a state)))
|
|
(* 3 y)))
|
|
|
|
(defun gen-ran (n)
|
|
(declare (fixnum n)
|
|
(optimize (speed 3) (safety 0)))
|
|
(let ((y 0d0))
|
|
(dotimes (k n)
|
|
(declare (fixnum k))
|
|
(setf y (random 1d0)))
|
|
(* 3 y)))
|
|
|
|
(defun run-mrg32k3a ()
|
|
(declare (inline gen-ran))
|
|
(gen-ran 1000000))
|
|
|
|
;; EOF
|