Pi

[hsmyers]

Found this while wandering around the internet:

Code: Select all

; Copyright (c) 2008 the authors listed at the following URL, and/or
; the authors of referenced articles or incorporated external code:
; http://en.literateprograms.org/Pi_with_Machin's_formula_(Lisp)?action=history&offset=20060307031705
;
; Permission is hereby granted, free of charge, to any person obtaining
; a copy of this software and associated documentation files (the
; "Software"), to deal in the Software without restriction, including
; without limitation the rights to use, copy, modify, merge, publish,
; distribute, sublicense, and/or sell copies of the Software, and to
; permit persons to whom the Software is furnished to do so, subject to
; the following conditions:
;
; The above copyright notice and this permission notice shall be
; included in all copies or substantial portions of the Software.
;
; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
; MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
; IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
; CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
; TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
; SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
;
; Retrieved from: http://en.literateprograms.org/Pi_with_Machin's_formula_(Lisp)?oldid=2756

(define (arccot x unity)
  (let ((xpower (floor (/ unity x))))
    (arccot-plus-helper (* x x) 1 xpower)))

(define (arccot-plus-helper xsq n xpower)
  (let ((term (floor (/ xpower n))))
    (if (= term 0)
        0
      (+ (arccot-minus-helper xsq (+ n 2) (floor (/ xpower xsq)))
         term))))

(define (arccot-minus-helper xsq n xpower)
  (let ((term (floor (/ xpower n))))
    (if (= term 0)
        0
      (- (arccot-plus-helper xsq (+ n 2) (floor (/ xpower xsq)))
         term))))

(define (pidigits digits)
    (letn (
      (unity (pow 10 (+ digits 10)))
      (pi (* 4 (- (* 4 (arccot 5 unity)) (arccot 239 unity)))))
      (floor (/ pi (pow 10 10)))))

I converted it to newLISP by the obvious means (defun to (define. Didn't have a (expt function so I substituted (pow. On my 32bit machine (pidigits works up to 8 digits(yeah I know--- impressive!) then I get the following:

Code: Select all

> (pidigits 6)
3141592
> (pidigits 7)
31415926
> (pidigits 8)
314159265
> (pidigits 9)
-791741031

I'm curious if this is some sort of limit reached with 32bit math? Could someone with a 64bit system give this a run?

--hsm

[Lutz]

I'm curious if this is some sort of limit reached with 32bit math? Could someone with a 64bit system give this a run?

in newLISP all integer math using the +,-,*,/ is 64bit. But you are mixing it with floating point math when using the power and floor functions. Then, when dividing using / you truncate the values to integers again. Convert everthing to pure float using add,sub,mul and div and it will work.

Read here:

http://www.newlisp.org/downloads/newlis ... #int_float

and look here to see how to format the out-coming float number back to an integer:

http://www.newlisp.org/downloads/newlis ... tml#format

[hsmyers]

Lutz Good catch! Great improvement. Now the limit is understandable:

Code: Select all

C:\>newlisp H:\pi.lsp
newLISP v.9.3.6 on Win32 IPv4, execute 'newlisp -h' for more info.

> (pidigits 18)
"3.141592653589793280"
> (pidigits 19)
"9.223372036854775807"
>

Which I ascribe to the cap on 64 bit unsigned integers being very near 20 digits long at 18,446,744,073,709,551,615, and I presume the next step probably craps out accordingly.

--hsm