David Steuber wrote:
> Just for a lark, I implemented an infinite series for calculating Pi
> from David Blatner's "The Joy of Pi" (the title actually uses the
> Greek letter). I used one of Newton's series on page 42. This is the
> first time I've written a Lisp program that runs faster in CLisp than
> in SBCL.
CLISP does very well on programs that spend most of their time
in bignum arithmetic. Since most of that time is in a library, which
in CLISP is written in C and is highly optimized, the fact that CLISP
doesn't compile to machine code is not as important.
Paul
David Steuber <·····@david-steuber.com> wrote:
> Just for a lark, I implemented an infinite series for calculating Pi
> from David Blatner's "The Joy of Pi" (the title actually uses the
> Greek letter). I used one of Newton's series on page 42. This is the
> first time I've written a Lisp program that runs faster in CLisp than
> in SBCL. OpenMCL runs it very slowly.
First I've tested it with LispWorks on my Pentium 1.5 GHz on Windows,
which is VERY slow:
CL-USER > (time (pi-decimal 1025 2000))
Timing the evaluation of (PI-DECIMAL 1025 2000)
; Loading fasl file C:\Program Files\Xanalys\LispWorks\lib\4-3-0-0
\modules\util\callcoun.fsl
[...]
user time = 269.557
system time = 0.360
Elapsed time = 0:04:44
Allocation = 11562961928 bytes standard / 573507 bytes conses
0 Page faults
Calls to %EVAL 35
NIL
With CLISP it is nearly 20 times faster (compiled):
Real time: 15.902867 sec.
Run time: 15.302003 sec.
Space: 26241288 Bytes
GC: 50, GC time: 0.0801152 sec.
NIL
> It might be interesting to
> make this program fast. Mostly that would involve using a faster
> computing method.
you are right. Using a spigot algorithm (translated from tcl:
http://wiki.tcl.tk/8401 ) you can write it like this:
(defun pi-2000 ()
(let ((e 0)
(f (make-array 3501 :initial-element 2000)))
(loop for c from 3500 above 0 by 14 do
(let ((d 0))
(loop for b from c above 0 do
(let ((g (1- (* 2 b))))
(setf d (+ (* d b) (* (aref f b) 10000)))
(setf (aref f b) (mod d g))
(setf d (floor d g))))
(multiple-value-bind (div mod) (floor d 10000)
(format t "~4,'0D" (+ e div))
(setf e mod))))))
In LispWorks:
CL-USER > (time (pi-2000))
Timing the evaluation of (PI-2000)
[...]
user time = 0.881
system time = 0.000
Elapsed time = 0:00:01
Allocation = 12470624 bytes standard / 85030 bytes conses
0 Page faults
Calls to %EVAL 33
NIL
and in CLISP:
Real time: 4.987171 sec.
Run time: 4.897042 sec.
Space: 5802196 Bytes
GC: 11, GC time: 0.0100144 sec.
NIL
--
Frank Bu�, ··@frank-buss.de
http://www.frank-buss.de, http://www.it4-systems.de
Frank Buss <··@frank-buss.de> wrote:
> (defun pi-2000 ()
sorry, wrong name, it calculates only 1000 digits, like the other
algorithm, to compare the time. Use 7001 and 7000 to calculate 2000 digits,
which is fast, too.
--
Frank Bu�, ··@frank-buss.de
http://www.frank-buss.de, http://www.it4-systems.de
Frank Buss <··@frank-buss.de> writes:
> David Steuber <·····@david-steuber.com> wrote:
>
> > Just for a lark, I implemented an infinite series for calculating Pi
> > from David Blatner's "The Joy of Pi" (the title actually uses the
> > Greek letter). I used one of Newton's series on page 42. This is the
> > first time I've written a Lisp program that runs faster in CLisp than
> > in SBCL. OpenMCL runs it very slowly.
>
> First I've tested it with LispWorks on my Pentium 1.5 GHz on Windows,
> which is VERY slow:
Did you compile pi-decimal?
Gregm
Greg Menke <··········@toadmail.com> wrote:
> Did you compile pi-decimal?
yes. Perhaps the bignum support is slow in LispWorks or I did something
other wrong.
--
Frank Bu�, ··@frank-buss.de
http://www.frank-buss.de, http://www.it4-systems.de
Frank Buss <··@frank-buss.de> writes:
> Greg Menke <··········@toadmail.com> wrote:
>
> > Did you compile pi-decimal?
>
> yes. Perhaps the bignum support is slow in LispWorks or I did something
> other wrong.
>
Did adding (declare ... ) for the types have any effect? Interesting
you found such a big performance difference though. I've found
Lispworks vastly faster than clisp in other ways though. OTOH,
different implementations have different strengths...
Gregm
Here's a Series version of the "pi-decimal" program. It is just about
identical in performance to the non-series version.
Running under CMUCL on my machine (1.4GHz Athlon) it takes about 10
seconds. The CLISP version takes about 5 seconds. I couldn't get it
to compile under SBCL. Does the Series package compile under SBCL?
I also tried this with ACL 6.2 but it's about 2 orders of magnitude
slower. I guess bignum performance wasn't a priority in that
implementation.
Disclaimer: I'm not a Series expert. Possibly someone else could do
this in a clever way. Mine is more a brute-force translation of the
original. I do like the way the Series version cleanly unfolds the
pieces of the computation without compromising efficiency.
Anyway here it is:
;;; Calculate pi using one of Newton's series:
;;;
;;; pi/6 = 1/2 + 1/2(1/3x2^3) + ((1x3)/(2x4))(1/5x2^2) +
;;; ((1x3x5)/(2x4x6))(1/7x2^7) + ...
;;;
;;; Series version.
;;;
(eval-when (:compile-toplevel :load-toplevel)
;; This is if you have installed Series as a CMUCL subsystem.
#+:CMU
(require :series)
#-:CMU
(progn
(load "s-package")
(load "s-code")))
(eval-when (:compile-toplevel :load-toplevel)
(series::install))
(defun alternating-number-series (first-number)
(declare (optimizable-series-function 1))
(scan-range :from first-number :by 2 :type integer))
(defun numerator-series ()
(declare (optimizable-series-function 1))
(collecting-fn 'integer
#'(lambda () 1)
#'*
(alternating-number-series 1)))
(defun denominator-series ()
(declare (optimizable-series-function 1))
(collecting-fn 'integer
#'(lambda () 1)
#'*
(alternating-number-series 2)))
(defun coefficient-series ()
(declare (optimizable-series-function 1))
(alternating-number-series 3))
(defun term-series ()
(declare (optimizable-series-function 1))
(mapping ((n (numerator-series))
(d (denominator-series))
(c (coefficient-series)))
(* (/ n d) (/ 1 (* c (expt 2 c))))))
(defun newton-pi (terms)
(let ((term-series (term-series)))
(* 6 (+ 1/2 (collect-sum (subseries term-series 0 terms) 'rational)))))
(defun digit-series (number)
(declare (optimizable-series-function 1))
(scan-fn '(values fixnum rational)
#'(lambda () (floor number))
#'(lambda (number remainder)
(declare (ignore number))
(floor (* 10 remainder)))))
(defun pi-decimal (digits &optional (terms 100))
(let* ((pi-rational (newton-pi terms)))
(iterate ((digit (digit-series pi-rational))
(i (scan-range :from 0 :by 1 :below digits)))
(princ digit)
(when (and (> i 0) (= 0 (mod i 4)))
(princ #\Space))
(when (and (> i 0) (= 0 (mod i 64)))
(terpri))
(when (= i 0) (princ ".") (terpri)))))
--
Fred Gilham ······@csl.sri.com
"In the 20th century, more citizens were killed by their own
governments than by foreign enemies....totalitarianism first of all
regards its own people as the enemy." --- Arnold Beichman
I wrote:
> (defun newton-pi (terms)
> (let ((term-series (term-series)))
> (* 6 (+ 1/2 (collect-sum (subseries term-series 0 terms) 'rational)))))
which has a little left-over cruft from the development process.
(defun newton-pi (terms)
(* 6 (+ 1/2 (collect-sum (subseries (term-series) 0 terms) 'rational)))
Similarly,
> (defun pi-decimal (digits &optional (terms 100))
> (let* ((pi-rational (newton-pi terms)))
> (iterate ((digit (digit-series pi-rational))
> (i (scan-range :from 0 :by 1 :below digits)))
> (princ digit)
> (when (and (> i 0) (= 0 (mod i 4)))
> (princ #\Space))
> (when (and (> i 0) (= 0 (mod i 64)))
> (terpri))
> (when (= i 0) (princ ".") (terpri)))))
should be
(defun pi-decimal (digits &optional (terms 100))
(iterate ((digit (digit-series (newton-pi terms)))
(i (scan-range :from 0 :by 1 :below digits)))
(princ digit)
(when (and (> i 0) (= 0 (mod i 4)))
(princ #\Space))
(when (and (> i 0) (= 0 (mod i 64)))
(terpri))
(when (= i 0) (princ ".") (terpri)))))
--
Fred Gilham ······@csl.sri.com | If sophists such as Richard Rorty
are correct, and the West has moved from a post-religious age to a
post-metaphysical age, then there is literally nothing Western left
about the West to defend. --David Dieteman
I messed about with this some more today.
Frank Buss <··@frank-buss.de> writes:
> David Steuber <·····@david-steuber.com> wrote:
>
> > It might be interesting to make this program fast. Mostly that
> > would involve using a faster computing method.
>
> you are right. Using a spigot algorithm (translated from tcl:
> http://wiki.tcl.tk/8401 ) you can write it like this:
>
> (defun pi-2000 ()
> (let ((e 0)
> (f (make-array 3501 :initial-element 2000)))
> (loop for c from 3500 above 0 by 14 do
> (let ((d 0))
> (loop for b from c above 0 do
> (let ((g (1- (* 2 b))))
> (setf d (+ (* d b) (* (aref f b) 10000)))
> (setf (aref f b) (mod d g))
> (setf d (floor d g))))
> (multiple-value-bind (div mod) (floor d 10000)
> (format t "~4,'0D" (+ e div))
> (setf e mod))))))
This is a very interesting algorithm and I can't figure out how it
works. I may have to get the book, "Pi Unleashed" that is reviewed
here:
http://www.maa.org/reviews/piunleashed.html
I did some hunting around for the article in AMM that was cited on the
Wiki page as well as here:
http://home.att.net/~srschmitt/pi-spigot.html
But I've had no luck so far.
I also did some minor tweaking of the code:
(defun pi-spigot ()
(let ((e 0)
(f (make-array 3501 :initial-element 2000)))
(symbol-macrolet ((f[b] (aref f b)))
(loop for c from 3500 above 0 by 14 do
(let ((d 0))
(loop for b from c above 0 do
(let ((g (1- (* 2 b))))
(multiple-value-setq (d f[b])
(floor (+ (* d b) (* f[b] 10000)) g))))
(multiple-value-bind (div mod) (floor d 10000)
(format t "~4,'0D" (+ e div))
(setf e mod)))))))
The symbol-macrolet allowed me to dump the redundant call to mod and
use multiple-value-setq. I also dropped a setf to d. It runs a shade
faster, but the tweaks didn't help me gain any insight into the
algorithm. My hope was to find a way to do away with the large
array. Say that three times fast. The array holds state information
used for later digits. I guess the array doesn't count as part of the
outputs in, "A spigot algorithm produces its outputs incrementally,
and does not reuse them after producing them."
Bruno Haible is mentioned in the book review for writing a cln library
that is included with Pi Unleashed. I thought that was pretty neat.
--
An ideal world is left as an excercise to the reader.
--- Paul Graham, On Lisp 8.1
Following up on myself here, mostly to use Google as a back up for my
file :-)
I cleaned up some of the Pi code I have. Fred Gilham's use of series
is pretty interesting, but I don't have that installed. As I said
before, I also like Common Lisp's closures. I don't know what
performance hit I take using them, but they make these sums easier to
think about.
In the code below, I have three different series. I haven't done any
timings, but I suspect that the series that converge faster (ie,
produce more digits of Pi for n terms) are more efficient. The best
series seem to use roots which I have not figured out how to code.
Here is pi.lisp as it exists now:
;;; Calculate Pi using one of Newton's series:
;;;
;;; Pi/6 = 1/2 + 1/2(1/3x2^3) + ((1x3)/(2x4))(1/5x2^2) +
;;; ((1x3x5)/(2x4x6))(1/7x2^7) + ...
;;;
(defun make-alternating-number-generator (first-number)
(let ((num (- first-number 2)))
(lambda () (incf num 2))))
(defun make-term-generator ()
(let ((numerator-gen (make-alternating-number-generator 1))
(denominator-gen (make-alternating-number-generator 2))
(coeff-gen (make-alternating-number-generator 3)))
(let ((n 1) (d 1))
(lambda ()
(setf n (* n (funcall numerator-gen))
d (* d (funcall denominator-gen)))
(let ((c (funcall coeff-gen)))
(* (/ n d) (/ 1 (* c (expt 2 c)))))))))
(defun newton-pi (terms)
(let ((a 1/2)
(term (make-term-generator)))
(dotimes (i terms)
(incf a (funcall term)))
(values (* 6 a) (* 6 (+ a (funcall term))))))
;;; Calculate Pi using the formula on page 118 of Joy of Pi.
;;; It is also at this URL:
;;;
;;; http://mathworld.wolfram.com/BBPFormula.html
;;;
(defun make-bbp-term-generator (&optional (nth-term 0))
(let ((n (1- nth-term)))
(lambda ()
(incf n)
(let ((8n (* 8 n)))
(* (expt 1/16 n)
(- (/ 4 (+ 8n 1))
(/ 2 (+ 8n 4))
(/ 1 (+ 8n 5))
(/ 1 (+ 8n 6))))))))
(defun bbp-pi (terms)
(let ((a 0)
(termgen (make-bbp-term-generator)))
(dotimes (i terms)
(incf a (funcall termgen)))
(values a (+ a (funcall termgen)))))
;;; Calculate Pi using Rabinowitz and Wagon series:
;;;
;;; inf n!^2 x 2^(n+1)
;;; Pi = sigma --------------
;;; n=0 (2n + 1)!
;;;
(let ((factorials (make-array 100 :adjustable t :initial-element nil)))
(defun n! (n)
"Returns n!. Memoization via an array for speed is done."
(when (< (length factorials) n)
(adjust-array factorials (+ n 100) :initial-element nil))
(if (< n 1)
1
(let ((f (aref factorials (1- n))))
(if f
f
(setf (aref factorials (1- n))
(* n (n! (1- n)))))))))
(defun make-rabinowitz-wagon-term-generator (&optional (nth-term 0))
(let ((n (1- nth-term)))
(lambda ()
(incf n)
(/ (* (expt (n! n) 2)
(expt 2 (1+ n)))
(n! (1+ (* 2 n)))))))
(defun rabinowitz-wagon-pi (terms)
(let ((a 0)
(termgen (make-rabinowitz-wagon-term-generator)))
(dotimes (i terms)
(incf a (funcall termgen)))
(values a (+ a (funcall termgen)))))
;;; Decimal expansion of above series. Pi is returned as a rational
;;; approximation along with the next term so that
;;; print-pi-decimal-expansion can terminate when the digits diverge.
(defun make-digit-generator (r)
(lambda ()
(multiple-value-bind (digit rem) (floor r)
(setf r (* 10 rem))
digit)))
(defun print-pi-decimal-expansion (fn terms)
(multiple-value-bind (p1 p2) (funcall fn terms)
(let ((g1 (make-digit-generator p1))
(g2 (make-digit-generator p2)))
(loop for d1 = (funcall g1) then (funcall g1)
for d2 = (funcall g2) then (funcall g2)
for i from 0
while (equal d1 d2)
do (progn
(princ d1)
(when (and (> i 0) (= 0 (mod i 5)))
(princ #\Space))
(when (and (> i 0) (= 0 (mod i 60)))
(terpri))
(when (= i 0) (princ ".") (terpri)))
finally (return i)))))
;;; Spigot algorithm
(defun pi-1000 ()
(let ((e 0)
(f (make-array 3501 :initial-element 2000)))
(loop for c from 3500 above 0 by 14 do
(let ((d 0))
(loop for b from c above 0 do
(let ((g (1- (* 2 b))))
(setf d (+ (* d b) (* (aref f b) 10000)))
(setf (aref f b) (mod d g))
(setf d (floor d g))))
(multiple-value-bind (div mod) (floor d 10000)
(format t "~4,'0D" (+ e div))
(setf e mod))))))
(defun pi-spigot ()
(let ((e 0)
(f (make-array 3501 :initial-element 2000)))
(symbol-macrolet ((f[b] (aref f b)))
(loop for c from 3500 above 0 by 14 do
(let ((d 0))
(loop for b from c above 0 do
(let ((g (1- (* 2 b))))
(multiple-value-setq (d f[b])
(floor (+ (* d b) (* f[b] 10000)) g))))
(multiple-value-bind (div mod) (floor d 10000)
(format t "~4,'0D" (+ e div))
(setf e mod)))))))
--
An ideal world is left as an excercise to the reader.
--- Paul Graham, On Lisp 8.1
David Steuber wrote:
> I cleaned up some of the Pi code I have. Fred Gilham's use of
> series is pretty interesting, but I don't have that installed.
Someone took the Series package and complexicated it. So I used my
decomplexicator on it and put it on my ftp site.
ftp.csl.sri.com:/pub/users/gilham/series.lisp
This is the 2.2.7 version from Sourceforge after running it through
the decomplexicator. It runs under CMUCL and CLISP at least.
Just compile this file and load it. Then if you want the series
symbols to be in your current package type (series::install).
Doesn't include documentation --- you'll have to snag that off the web
yourself.
If you have CMUCL installed in, say, /usr/local, and you want to use
Series as a library package, you can copy series.x86f (or whatever the
fasl file is called for your CMUCL) to
/usr/local/lib/cmucl/lib/subsystems/series-library.x86f. Then you can
do
(require :series)
whenever you want to use the package. Very convenient. I do this for
a number of packages.
You can do this for pretty much any package by just cat-ing all the
fasl files together (in the right order, of course) and copying them
to the lib/subsystems directory, naming them <package>-library.x86f.
--
Fred Gilham ······@csl.sri.com
Top five reasons why the God of Jesus Christ makes a better God than
`Caesar':
5. God only wants 10%.
4. God has fewer laws.
3. God has a clear understanding of the difference between destroying
things and saving them.
2. God has a better retirement package.
1. Caesar wants you to send your sons to die for him. God sent his
Son to die for you.
>>>>> "Fred" == Fred Gilham <······@snapdragon.csl.sri.com> writes:
Fred> David Steuber wrote:
>> I cleaned up some of the Pi code I have. Fred Gilham's use of
>> series is pretty interesting, but I don't have that installed.
Fred> Someone took the Series package and complexicated it. So I used my
What does complexicated mean?
Ray
Ray Toy asked:
> What does complexicated mean?
Well, the Series package used to be a single file (not counting the
test code). Someone split it up so you have to load at least 2 files
to get it to run.
This has apparently proved to be an obstacle to people.
My decomplexicator, aka emacs, allowed me to put series back into a
single file that can be compiled, loaded and run in a simple fashion.
--
Fred Gilham ······@csl.sri.com
"Due to inclement weather and massive amounts of ice everywhere,
tonight's Healthy Environment Forum on Global Warming with Dr. Patz
has been postponed", wrote Sarah Doll of the Oregon Environmental
Council in an email sent to the press, "Sorry for any inconvenience
and hope you are staying warm." --- John Daly
>>>>> "Fred" == Fred Gilham <······@snapdragon.csl.sri.com> writes:
Fred> Ray Toy asked:
>> What does complexicated mean?
Fred> Well, the Series package used to be a single file (not counting the
Fred> test code). Someone split it up so you have to load at least 2 files
Fred> to get it to run.
Ah, I see. It's 3 files now. One for the package definitions, which
I think Tim Bradshaw said was a good idea; one for the series stuff;
and one to install it. I think it's a bad idea to install series
automatically, so we wanted the user to explicitly install it.
Ray