1. Robert Smith
  2. computable-reals

Source

computable-reals / reals.lisp

Diff from to

File reals.lisp

  • Ignore whitespace
 ;;;; Computable real numbers
+
 ;;;; Michael Stoll
 ;;;; 1989-06-11, 1989-06-12, 1989-06-13, 1989-06-14, 1989-06-17, 1989-06-30
 
+;;;; Modified by Robert Smith
+;;;; 2011-12-07
 
 ;;;;   I N T E R N A L   S T R U C T U R E S   A N D   I N T E R F A C E
 ;;;;   -----------------------------------------------------------------
 
-;; Package for computable real numbers
-
-
 (in-package #:computable-reals)
 
 ;;; Computable reel numbers are rational numbers or structures:
 
 (defstruct (c-real (:copier nil)
                    (:print-function print-c-real))
-           (value     0      :type integer)
-           (precision -1     :type (integer -1 *))
-           (compute   nil    :type (function ((integer 0 *)) integer)
-                             :read-only t))
+  (value     0      :type integer)
+  (precision -1     :type (integer -1 *))
+  (compute   nil    :type (function ((integer 0 *)) integer)
+                    :read-only t))
 
 (deftype CREAL ()
   "type of the computable real numbers"
 
 ;;; A few shortcuts for signalling type errors.
 
+;;; TODO: these can be removed if conditions and asserts are used.
+
 (defun cr-error (fun x)
   (error "~S: ~S is not a computable real number" fun x))
 
                  (erg 0 (+ erg (round y n))))
                 ((< (abs y) n) (round-cr erg k0))))))))
 
-;;; Now log2:
-
-(defconstant LOG2-R (+r (ash-r (log-r2 1/7) 1) (log-r2 1/17))
-  "log(2) as CREAL")
-
-(get-approx log2-r 200)  ; precompute to ca. 60 decimal digits
-
 ;;; (log-r1 x) takes a creal x from [1,2] and returns log(x) as creal
 
 (defun log-r1 (x)
                  (erg 0 (+ erg (round y n))))
                 ((< (abs y) n) (round-cr erg k0))))))))
 
-;;; Now pi:
-
-(defconstant PI-R (-r (ash-r (atan-r1 1/10) 5)
-                      (ash-r (atan-r1 1/515) 4)
-                      (ash-r (atan-r1 1/239) 2))
-  "pi as CREAL")
-
-(defconstant 2PI-R (ash-r pi-r 1) "2*pi as CREAL")
-(defconstant PI/2-R (ash-r pi-r -1) "pi/2 as CREAL")
-(defconstant PI/4-R (ash-r pi-r -2) "pi/4 as CREAL")
-
-(get-approx 2pi-r 200) ; compute to ca. 60 decimal digits
-(get-approx pi-r 200)
-(get-approx pi/2-r 200)
-(get-approx pi/4-r 200)
-
 ;;; (atan-r0 x) takes a creal x and returns atan(x) as creal.
 
 (defun atan-r0 (x)
   (let ((a (get-approx x 3)))
     (cond ((<= -3 a 3) (atan-r1 x))
           ((< a -3) (minus-r (atan-r0 (minus-r x)))) ; atan(x) = -atan(-x)
-          ((< 3 a 17) (+r pi/4-r (atan-r1 (transf x))))
+          ((< 3 a 17) (+r +pi/4-r+ (atan-r1 (transf x))))
                     ; atan(x) = pi/4 + atan((x-1)/(x+1))
-          (t (-r pi/2-r (atan-r1 (invert-r x))))))) ; atan(x) = pi/2 - atan(1/x)
+          (t (-r +pi/2-r+ (atan-r1 (invert-r x))))))) ; atan(x) = pi/2 - atan(1/x)
 
 ;;; (atan-r x [y]) computes the arctangent of the creals x (and y if given)
 
                         (- (integer-length ay)) (- nx))))
           (cond ((and (plusp sx) (>= mx-my 0)) (atan-r0 (/r y x)))
                 ((and (plusp sy) (<= mx-my 0))
-                 (-r pi/2-r (atan-r0 (/r x y))))
+                 (-r +pi/2-r+ (atan-r0 (/r x y))))
                 ((and (minusp sy) (<= mx-my 0))
-                 (minus-r (+r (atan-r0 (/r x y)) pi/2-r)))
+                 (minus-r (+r (atan-r0 (/r x y)) +pi/2-r+)))
                 ((and (minusp sx) (minusp sy) (>= mx-my 0))
-                 (-r (atan-r0 (/r y x)) pi-r))
-                (t (+r (atan-r0 (/r y x)) pi-r))))))))
+                 (-r (atan-r0 (/r y x)) +pi-r+))
+                (t (+r (atan-r0 (/r y x)) +pi-r+))))))))
 
 ;;; (sin-r1 x) takes a creal x with |x|<4 and returns sin(x) as creal.
 
   ;; remember sin(k*2pi + y) = sin(y)
   (if (eql x 0)
     0
-    (multiple-value-bind (q r) (round-r x 2pi-r 10)
+    (multiple-value-bind (q r) (round-r x +2pi-r+ 10)
       (declare (ignore q))
       (sin-r1 r))))
 
   ;; remember cos(k*2pi + y) = cos(y)
   (if (eql x 0)
     1
-    (multiple-value-bind (q r) (round-r x 2pi-r 10)
+    (multiple-value-bind (q r) (round-r x +2pi-r+ 10)
       (declare (ignore q))
       (cos-r1 r))))