Commits

Robert Smith committed 8b1cfe8

original computable reals library by Michael Stoll in 1989

  • Participants

Comments (0)

Files changed (2)

+;;;; The code below can be used freely under the following conditions:
+;;;; * Due credit is given to the original author.
+;;;; * You use it at your own risk. The code is experimental and meant
+;;;;   as a case study to show that "lazy reals" can be implemented.
+;;;;   In particular, no guarantee is given that the code conforms to
+;;;;   the specification given below.
+;;;;
+;;;; If you use it in publications (including software), please let me
+;;;; know.
+;;;;
+;;;; Note: The main focus here is on correct results (which are, however,
+;;;; not guaranteed, see above). It is certainly possible to improve the
+;;;; efficiency of the code quite a bit.
+;;;;
+;;;; Michael Stoll, December 11, 2009
+;;;; firstname.lastname@uni-bayreuth.de
+
+;;;; ======================================================================
+
+;;;; Computable real numbers
+;;;; Michael Stoll
+;;;; 1989-06-11, 1989-06-12, 1989-06-13, 1989-06-14, 1989-06-17, 1989-06-30
+
+
+;;;;                        I N T R O D U C T I O N
+;;;;                        =======================
+
+;;;; Computable real numbers x are interpreted as (potentially) infinite
+;;;; fractions in base 2 that are specified through a rule for computation
+;;;; of an integer a with |(2^k)*x - a| <= 1 for any k>=0.
+
+;;;; The internal data structure should not be accessed.
+;;;; The interface for the outside world is as follows:
+
+;;;; The type CREAL is a supertype of the type RATIONAL. (CREAL-P x), for
+;;;; an object x, returns T if x is of type CREAL, otherwise NIL.
+;;;;
+;;;; (APPROX-R x k), for a CREAL x and an integer k>=0, returns an integer
+;;;; a with |(2^k)*x - a| < 1.
+;;;;
+;;;; (MAKE-REAL fun) returns the real number given by fun. Here fun is a
+;;;; function taking an argument k, that computes a as above.
+
+;;;; CREALs are output by print etc. as a decimal fraction. The error hereby
+;;;; is at most one unit in the last digit that was output. The number of
+;;;; decimal digits after the decimal point is defined through the dynamic
+;;;; variable *PRINT-PREC*.
+;;;;
+;;;; For comparison operations etc. a precision threshold is used. It is
+;;;; defined through the dynamic variable *CREAL-TOLERANCE*. Its value should
+;;;; be a nonnegative integer n, meaning that numbers are considered equal
+;;;; if they differ by at most 2^(-n).
+
+
+;;;;      E X P O R T E D   F U N C T I O N S   A N D   C O N S T A N T S
+;;;;      ===============================================================
+
+;;;; The following functions, constants and variables are exported. (The
+;;;; package is named "REALS".)
+;;;; CREAL                  type        type of the computable real numbers
+;;;; CREAL-P object         function    tests for type CREAL
+;;;; *PRINT-PREC*           variable    specifies precision of output
+;;;; *CREAL-TOLERANCE*      variable    precision threshold for comparison
+;;;; APPROX-R x:creal k:int>=0
+;;;;                        function    returns approximation of x to k digits
+;;;; MAKE-REAL function     function    creates object of type CREAL
+;;;; RAW-APPROX-R x:creal   function    returns 3 values a,n,s with:
+;;;;                                    if a = 0: |x| <= 2^(-n), s = 0
+;;;;                                        and n >= *CREAL-TOLERANCE*
+;;;;                                    else: a0 integer > 4, n0 integer >=0,
+;;;;                                        s = +1 or -1, and sign(x) = s,
+;;;;                                        (a-1)*2^(-n) <= |x| <= (a+1)*2^(-n)
+;;;; PRINT-R x:creal k:int>=0 &optional (flag t)
+;;;;                        function    outputs x with k decimal digits.
+;;;;                                    If flag is true, first a newline.
+;;;; +R {creal}*            function    computes the sum of the arguments
+;;;; -R creal {creal}*      function    computes negative or difference
+;;;; *R {creal}*            function    computes the product of the arguments
+;;;; /R creal {creal}*      function    computes reciprocal or quotient
+;;;; SQRT-R creal           function    computes the square root
+;;;; LOG2-R                 constant    log(2)
+;;;; PI-R                   constant    pi
+;;;; 2PI-R                  constant    2*pi
+;;;; PI/2-R                 constant    pi/2
+;;;; PI/4-R                 constant    pi/4
+;;;; LOG-R x:creal &optional b:creal
+;;;;                        function    computes the logarithm of n in base b;
+;;;;                                    default is the natural logarithm
+;;;; EXP-R creal            function    computes the exponential function
+;;;; EXPT-R x:creal y:creal function    computes x^y
+;;;; SIN-R creal            function    computes the sine
+;;;; COS-R creal            function    computes the cosine
+;;;; TAN-R creal            function    computes the tangent
+;;;; ATAN-R x:creal &optional y:creal
+;;;;                        function    computes the arctangent of x or
+;;;;                                    the phase angle of (x,y)
+;;;; ASH-R x:creal n:int    function    computes x * 2^n
+;;;; ROUND-R x:creal &optional y:creal
+;;;;                        function    computes two values q (integer) and r
+;;;;                                    (creal) with x = q*y + r and |r|<=|y|/2
+;;;;                                    according to the precision specified by
+;;;;                                    *CREAL-TOLERANCE*
+;;;; FLOOR-R x:creal &optional y:creal
+;;;;                        function    like ROUND-R, corresponding to floor
+;;;; CEILING-R x:creal &optional y:creal
+;;;;                        function    like ROUND-R, corresponding to ceiling
+;;;; TRUNCATE-R x:creal &optional y:creal
+;;;;                        function    like ROUND-R, corresponding to truncate
+
+
+;;;; In the code, exported symbols at their location of definition are
+;;;; emphasized in UPPERCASE.
+
+;;;; ==========================================================================
+
+;;;;   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
+(defpackage "REALS" (:use "COMMON-LISP"))
+(in-package "REALS")
+
+(export '(creal approx-r make-real creal-p print-r +r -r *r /r sqrt-r
+          log2-r log-r exp-r pi-r 2pi-r sin-r cos-r *print-prec* round-r
+          *creal-tolerance* ash-r raw-approx-r floor-r ceiling-r truncate-r
+          pi/2-r pi/4-r atan-r expt-r tan-r
+)        )
+
+;;; 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
+)          )
+
+(deftype CREAL () "type of the computable real numbers"
+  '(or rational c-real)
+)
+
+(defun CREAL-P (x) (or (rationalp x) (c-real-p x)))
+
+;; If r is a c-real with (c-real-value r) = a and (c-real-precision r) = k,
+;; then a*2^(-k) is an approximation of the value of the number represented
+;; by r that deviates from the actual value by at most 2^(-k).
+;; (c-real-compute r) is a function taking an argument k, that returns an
+;; approximation of precision 2^(-k) and returns the corresponding value a.
+
+;;; make-real creates a c-real from a computation function.
+
+(defun MAKE-REAL (comp) "This function creates CREALs"
+  (declare (type (function ((integer 0 *)) integer) comp))
+  (make-c-real :compute comp)
+)
+
+;;; The following function takes an object of type creal and a number k,
+;;; and returns an integer a with |a*2^(-k) - x| <= 2^(-k), where x denotes
+;;; the corresponding real number.
+
+(defun APPROX-R (x k) "This function computes approximations for CREALs"
+  (unless (creal-p x) (cr-error 'approx-r x))
+  (unless (and (integerp k) (>= k 0)) (nat-error 'approx-r k))
+  (get-approx x k)
+)
+
+(defun get-approx (x k)
+  (declare (type creal x) (type (integer 0 *) k))
+  (cond ((integerp x) (ash x k))
+        ((rationalp x) (round (ash (numerator x) k) (denominator x)))
+        ((c-real-p x)
+         (if (>= (c-real-precision x) k)
+           (ash (c-real-value x) (- k (c-real-precision x)))
+           (let ((a (funcall (c-real-compute x) k)))
+             (setf (c-real-value x) a (c-real-precision x) k)
+             a
+        )) )
+        (t (cr-error 'get-approx x))
+) )
+
+;;; A few shortcuts for signalling type errors.
+
+(defun cr-error (fun x)
+  (error "~S: ~S is not a computable real number" fun x)
+)
+(defun int-error (fun x)
+  (error "~S: ~S is not an integer" fun x)
+)
+(defun nat-error (fun x)
+  (error "~S: ~S is not a nonnegative integer" fun x)
+)
+
+;;;; ==========================================================================
+
+;;;;                           V A R I A B L E S
+;;;;                           -----------------
+
+;;; *print-prec* specifies how many digits after the decimal point are output
+;;; (by print etc.)
+
+(defparameter *PRINT-PREC* 20
+  "number of decimal digits after the decimal point during output of CREALs"
+)
+
+;;; *creal-tolerance* specifies the precision of comparison operations
+
+(defparameter *CREAL-TOLERANCE* 100
+  "precision threshold for the comparison of CREALs,
+denoting the number of binary digits after the decimal point"
+)
+
+;;;;                  A U X I L I A R Y   F U N C T I O N S
+;;;;                  -------------------------------------
+
+;;; The following functions perform rounding, when less precision is needed.
+
+(defun round-cr (a k)
+  (declare (type integer a) (type (integer 0 *) k))
+  (if (eql k 0)
+    a
+    (if (logbitp (1- k) a) (1+ (ash a (- k))) (ash a (- k)))
+) )
+
+;;; Auxiliary function for approximating.
+
+(defun raw-approx-cr (x)
+  (declare (type creal x))
+  (do* ((k 0 (+ k 4))
+        (a (get-approx x 0) (get-approx x k))
+        (crt (+ 2 *creal-tolerance*))
+       )
+       ((> (abs a) 4) (values (abs a) k (signum a)))
+    (when (> k crt) (return (values 0 (- k 3) 0)))
+) )
+
+(defun RAW-APPROX-R (x) "This function returns an approximation for CREALs"
+  (unless (creal-p x) (cr-error 'raw-approx-r x))
+  (raw-approx-cr x)
+)
+
+;;;;                        P R I N T   F U N C T I O N
+;;;;                        ---------------------------
+
+;;; Small auxiliary function for avoiding repeated computation:
+
+(let* ((pp *print-prec*) (tenpowerpp (expt 10 pp)))
+  (declare (type (integer 0 *) pp tenpowerpp))
+  (defun tenpower (k)
+    (declare (type (integer 0 *) k))
+    (if (eql k pp)
+      tenpowerpp
+      (let ((zhk (expt 10 k)))
+        (when (eql k *print-prec*) (setq pp k tenpowerpp zhk))
+        zhk
+) ) ) )
+
+;;; The next function performs output to k digits after the decimal point,
+;;; ensuring an error of at most one unit on the last digit.
+
+(defun PRINT-R (x k &optional (flag t) (stream *standard-output*))
+  "output function for CREALs"
+  ;; flag /= NIL: the value is printed in a new line
+  ;; flag = NIL: no linefeed
+  (unless (creal-p x) (cr-error 'print-r x))
+  (unless (and (integerp k) (>= k 0)) (nat-error 'print-r k))
+  (unless (streamp stream) (error "~S: ~S is not a stream" 'print-r stream))
+  (creal-print x k flag stream)
+)
+
+(defun creal-print (x k flag stream)
+  (declare (type creal x) (type (integer 0 *) k) (type stream stream))
+  (let* ((k1 (tenpower k))
+         (n (1+ (integer-length k1)))
+         (x1 (get-approx x n))
+         (sign (signum x1))
+         (x2 (round-cr (* (abs x1) k1) n))
+         (*print-base* 10.)
+        )
+    (multiple-value-bind (vor nach) (floor x2 k1)
+      (when flag (terpri stream))
+      (write-char (if (minusp sign) #\- #\+) stream)
+      (prin1 vor stream)
+      (write-char #\. stream)
+      (let ((s (prin1-to-string nach)))
+        (write-string (make-string (- k (length s)) :initial-element #\0)
+                      stream
+        )
+        (write-string s stream)
+        (write-string "..." stream)
+        (values)
+) ) ) )
+
+(defun print-c-real (x stream d)
+  (declare (ignore d))
+  (creal-print x *print-prec* nil stream))
+
+;;;;                           A R I T H M E T I C
+;;;;                           -------------------
+
+;;; Now comes the addition.
+
+(defun +R (&rest args &aux (sn 0) (rl nil)) "addition of CREALs"
+  (declare (type rational sn) (type list #|(list creal)|# rl))
+  (dolist (x args)
+    (cond ((rationalp x) (setq sn (+ x sn)))
+          ((c-real-p x) (setq rl (cons x rl)))
+          (t (cr-error '+ x))
+  ) )
+  ;; sn = exact partial sum
+  ;; rl = list of the "real" real arguments
+  (let* ((n (length rl)) ; n = how many of them
+         (k1 (integer-length (if (integerp sn) n (1+ n))))
+           ; k1 = number of additional binary digits for the summands
+        )
+    (if (eql n 0)
+      sn        ; sum is exact
+      (make-real
+        #'(lambda (k &aux (k2 (+ k k1)))
+            (do ((sum (get-approx sn k2) (+ sum (get-approx (first l) k2)))
+                 (l rl (rest l))
+                )
+                ((null l) (round-cr sum k1))
+) ) ) )   ) )
+
+;;; Negation:
+
+(defun minus-r (x)
+  (cond ((rationalp x) (- x))
+        ((c-real-p x) (make-real #'(lambda (k) (- (get-approx x k)))))
+        (t (cr-error '- x))
+) )
+
+;;; Subtraction:
+
+(defun -R (x1 &rest args) "subtraction and negation of CREALs"
+  (if (null args)
+    (minus-r x1)
+    (+r x1 (minus-r (apply #'+r args)))
+) )
+
+;;; Now comes the multiplication.
+
+(defun *R (&rest args &aux (pn 1) (rl nil)) "Multiplication for CREALs"
+  (declare (type rational pn) (type list #|(list creal)|# rl))
+  (dolist (x args)
+    (cond ((rationalp x) (setq pn (* x pn)))
+          ((c-real-p x) (setq rl (cons x rl)))
+          (t (cr-error '* x))
+  ) )
+  ;; pn = product of the rational factors
+  ;; rl = list of the c-real factors
+  (when (or (eql pn 0) (null rl)) (return-from *r pn))
+  ;; If pn is a true fraction, handle it like a c-real.
+  (unless (integerp pn) (setq rl (cons pn rl) pn 1))
+  (let ((y (* (length rl) (abs pn))) (al nil) (nl nil) (ns 1) ll)
+    (dolist (x rl)
+      (multiple-value-bind (a0 n0) (raw-approx-cr x)
+        (setq al (cons (1+ a0) al)
+              nl (cons n0 nl)
+              y (* y (1+ a0))
+              ns (- ns n0)
+    ) ) )
+    (setq ll (mapcar #'(lambda (z m)
+                         (+ m ns (integer-length (1- (ceiling y z))))
+                       )
+                     al nl
+             )
+          rl (nreverse rl)
+    )
+    ;; rl = list of the factors (not including the integer pn)
+    ;; ll = list of the corresponding precision differences
+    ;; nl = list of the correspodning minimum precisions
+    (make-real
+      #'(lambda (k)
+          (let ((erg pn) (s (- k)) (rl rl) (ll ll) (nl nl) k1)
+            (loop (setq k1 (max (first nl) (+ k (first ll)))
+                        s (+ s k1)
+                        erg (* erg (get-approx (first rl) k1))
+                        rl (rest rl)
+                        ll (rest ll)
+                        nl (rest nl)
+                  )
+                  (when (null rl) 
+                    (return (if (minusp s)
+                              0
+                              (round-cr erg s)
+) ) )   ) ) )     ) )       )
+
+;;; Reciprocal:
+
+(defun invert-r (x)
+  (cond ((rationalp x) (/ x))
+        ((c-real-p x)
+         (multiple-value-bind (a0 n0) (raw-approx-cr x)
+           (when (eql a0 0) (error "division by 0"))
+           (let ((k1 (+ 4 (* 2 (- n0 (integer-length (1- a0))))))
+                 (k2 (1+ n0))
+                )
+             (make-real #'(lambda (k &aux (k0 (max k2 (+ k k1))))
+                            (round (ash 1 (+ k k0)) (get-approx x k0))
+        )) ) )            )
+        (t (cr-error '/ x))
+) )
+
+;;; Division:
+
+(defun /R (x1 &rest args) "division for CREALs"
+  (if (null args)
+    (invert-r x1)
+    (*r x1 (invert-r (apply #'*r args)))
+) )
+
+;;; Square root:
+
+(defun SQRT-R (x &aux s) "square root for CREALs"
+  (unless (creal-p x) (cr-error 'sqrt x))
+  (if (and (rationalp x) (>= x 0) (rationalp (setq s (sqrt x))))
+    s
+    (multiple-value-bind (a0 n0 s) (raw-approx-cr x)
+      (unless (plusp s)
+        (error "~S: attempting to compute the square root of a negative number"
+               'sqrt-r
+      ) )
+      (let ((k1 (1+ (ceiling (- n0 (integer-length (1- a0))) 2)))
+            (n1 (ceiling n0 2))
+           )
+        (make-real
+          #'(lambda (k &aux (k2 (max n1 (ceiling (+ k k1) 2)))
+                            (k3 (max 0 (- k -2 k1)))
+                    )
+              (round-cr (isqrt (ash (get-approx x (* 2 k2)) (* 2 k3)))
+                        (+ k3 k2 (- k))
+) ) ) ) )   ) )
+
+;;; Now comes a round function.
+;;; (round-r x y l) (x, y creal, l int>=0) returns two values q and r,
+;;; where q is an integer and r a creal, so that x = q*y + r and
+;;; |r| <= (1/2+2^(-l))*|y|. The default value of l is such that |r| exceeds
+;;; |y|/2 by at most 2^(- *CREAL-TOLERANCE*).
+;;; The third argument is specified only for internal purposes.
+
+(defun ROUND-R (x &optional (y 1) (l nil)) "round for CREALs"
+  (divide-r 'round #'round x y l)
+)
+
+(defun FLOOR-R (x &optional (y 1) (l nil)) "floor for CREALs"
+  (divide-r 'floor #'floor x y l)
+)
+
+(defun CEILING-R (x &optional (y 1) (l nil)) "ceiling for CREALs"
+  (divide-r 'ceiling #'ceiling x y l)
+)
+
+(defun TRUNCATE-R (x &optional (y 1) (l nil)) "truncate for CREALs"
+  (divide-r 'truncate #'truncate x y l)
+)
+
+(defun divide-r (name what x y l)
+  ; name = name of the calling function
+  ; what = #'round, #'floor, #'ceiling or #'truncate
+  (unless (creal-p x) (cr-error name x))
+  (unless (creal-p y) (cr-error name y))
+  (if (and (rationalp x) (rationalp y))
+    (funcall what x y)    ; for rational numbers use the common function
+    (multiple-value-bind (a0 n0) (raw-approx-cr y)
+      (when (eql a0 0) (error "~S: division by 0" name))
+      (when (null l)
+        (setq l (+ (integer-length a0) *creal-tolerance* (- n0)))
+      )
+      (let* ((x1 (abs (get-approx x n0)))
+             (m (max n0 (+ l 2 n0 (integer-length (+ x1 a0 -1))
+                           (* -2 (integer-length (1- a0)))
+             )  )       )
+             (q (funcall what (get-approx x m) (get-approx y m)))
+            )
+        (values q (rest-help-r x y (- q)))
+) ) ) )
+
+;; (rest-help-r x y q), with x,y creal, q integer, computes x + q*y.
+
+(defun rest-help-r (x y q)
+  (declare (type creal x y) (type integer q))
+  (if (eql q 0)
+    x
+    (let ((k1 (1+ (integer-length (1- (abs q))))))
+      (make-real
+        #'(lambda (k)
+            (round-cr (+ (ash (get-approx x (+ k 2)) (- k1 2))
+                         (* q (get-approx y (+ k k1)))
+                      )
+                      k1
+) ) ) )   ) )
+
+;;; Now comes the arithmetic shift function for infinite binary fractions:
+
+(defun ASH-R (x n) "shift function for CREALs"
+  (unless (creal-p x) (cr-error 'ash-r x))
+  (unless (integerp n) (int-error 'ash-r n))
+  (cond ((eql n 0) x)
+        ((integerp x)
+         (if (plusp n) (ash x n) (/ x (ash 1 (- n))))
+        )
+        ((rationalp x)
+         (if (plusp n)
+           (/ (ash (numerator x) n) (denominator x))
+           (/ (numerator x) (ash (denominator x) (- n)))
+        ))
+        ((plusp n) (make-real #'(lambda (k) (get-approx x (+ k n)))))
+        (t (make-real #'(lambda (k)
+                          (if (minusp (+ k n))
+                            (round-cr (get-approx x 0) (- (+ k n)))
+                            (get-approx x (+ k n))
+) )     )  )            ) )
+
+;;; Now we look at the most important transcendental functions.
+
+;;; (log-r2 x) takes a creal x |x|<=1/2 and returns log((1+x)/(1-x)) as creal.
+;;;   log((1+x)/(1-x)) = 2*(x + x^3/3 + x^5/5 + ... )
+
+(defun log-r2 (x)
+  (declare (type creal x))
+  (if (eql x 0)
+    0
+    (make-real
+      #'(lambda (k)
+          (let* ((k0 (integer-length (1- (integer-length k))))
+                     ; k0 = extra precision needed for partial sums
+                 (k1 (+ k k0 1)) ; k1 = total precision needed
+                                 ; (+1 because of factor 2)
+                 (ax (get-approx x (1+ k1)))
+                 (fx (round ax 2)) ; fx = k1-approximation of x
+                 (fx2 (round-cr (* ax ax) (+ k1 2))) ; fx2 = dito of x^2
+                )
+            (do ((n 1 (+ n 2))
+                 (y fx (round-cr (* y fx2) k1))
+                 (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)
+  (declare (type creal x))
+  (log-r2 (transf x))
+)
+
+;;; (transf x) takes a creal x from [1,2] and returns (x-1)/(x+1) as creal
+
+(defun transf (x)
+  (declare (type creal x))
+  (if (rationalp x)
+    (/ (1- x) (1+ x))
+    (make-real #'(lambda (k)
+                   (let ((a (get-approx x k)) (e (ash 1 k)))
+                     (round (ash (- a e) k) (+ a e))
+) ) )            ) )
+
+;;; Now the logarithm.
+
+(defun LOG-R (x &optional (b nil)) "logarithm for CREALs"
+  (unless (creal-p x) (cr-error 'log x))
+  (unless (or (null b) (creal-p b)) (cr-error 'log b))
+  (if b
+    (/r (log-r x) (log-r b))
+    ;; remember log(2^n * a) = n*log(2) + log(a)
+    (multiple-value-bind (a0 n0 s) (raw-approx-cr x)
+      (unless (plusp s)
+        (error "~S: attempt to compute the logarithm of a nonpositive number"
+               'log-r
+      ) )
+      (let ((shift (- (integer-length a0) 1 n0)))
+        (rest-help-r (log-r1 (ash-r x (- shift))) log2-r shift)
+) ) ) )
+
+;;; Now the exponential function.
+
+;;; (exp-r1 x) takes a creal x with |x| <= 1/2*log(2)
+;;; and returns exp(x) as creal
+
+(defun exp-r1 (x)
+  (declare (type creal x))
+  (make-real
+    #'(lambda (k)
+        (let ((m 3) (k2 (+ k 3)))
+          (loop (when (<= k2 (ash (- m 2) m)) (return))
+                (incf m)
+          )
+          (setq m (+ m 3) k2 (+ k m))
+          (do ((x1 (get-approx x k2))
+               (n 1 (1+ n))
+               (y (ash 1 k2) (round-cr (round (* y x1) n) k2))
+               (erg 0 (+ erg y))
+              )
+              ((eql y 0) (round-cr erg m))
+) )   ) ) )
+
+(defun EXP-R (x) "exponential function for CREALs"
+  (unless (creal-p x) (cr-error 'exp x))
+  ;; remember exp(a*log2 + b) = exp(b) * 2^a
+  (if (eql x 0)
+    1
+    (multiple-value-bind (q r) (round-r x log2-r 10)
+      (ash-r (exp-r1 r) q)
+) ) )
+
+;;; (expt-r x y) takes creals x,y and computes x^y
+
+(defun EXPT-R (x y &aux s) "exponentiation function for CREALs"
+  (unless (creal-p x) (cr-error 'expt x))
+  (unless (creal-p y) (cr-error 'expt y))
+  (cond ((eql y 0) 1)
+        ((integerp y)
+         (if (rationalp x) (expt x y) (expt-r1 x y))
+        )
+        ((and (rationalp y)
+              (eql 2 (denominator y))
+              (rationalp x)
+              (rationalp (setq s (sqrt x)))
+         )
+         (expt s (* 2 y))
+        )
+        (t (exp-r (*r y (log-r x))))
+) )
+
+(defun expt-r1 (x y)
+  (declare (type creal x) (integer y))
+  (cond ((minusp y) (expt-r1 (invert-r x) (- y)))
+        ((eql y 1) x)
+        ((evenp y) (expt-r1 (*r x x) (floor y 2)))
+        (t (*r x (expt-r1 (*r x x) (floor y 2))))
+) )
+
+;;; Now the trigonometric functions.
+
+;;; (atan-r1 x) takes a creal x with |x| <= 1/2 and returns atan(x) as creal
+
+(defun atan-r1 (x)
+  (declare (type creal x))
+  (if (eql x 0)
+    0
+    (make-real
+      #'(lambda (k)
+          (let* ((k0 (integer-length (1- (integer-length k))))
+                     ; k0 = extra precision needed for partial sums
+                 (k1 (+ k k0)) ; k1 = total precision needed
+                 (ax (get-approx x (1+ k1)))
+                 (fx (round ax 2)) ; fx = k1-approximation of x
+                 (fx2 (- (round-cr (* ax ax) (+ k1 2)))) ; fx2 = dito of -x^2
+                )
+            (do ((n 1 (+ n 2))
+                 (y fx (round-cr (* y fx2) k1))
+                 (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)
+  (declare (type creal 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))))
+                    ; 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)
+) ) )
+
+;;; (atan-r x [y]) computes the arctangent of the creals x (and y if given)
+
+(defun ATAN-R (x &optional (y nil)) "arctangent for CREALs"
+  (unless (creal-p x) (cr-error 'atan x))
+  (unless (or (null y) (creal-p y)) (cr-error 'atan y))
+  (if (null y)
+    (atan-r0 x)
+    (multiple-value-bind (ax nx sx) (raw-approx-cr x)
+      (multiple-value-bind (ay ny sy) (raw-approx-cr y)
+        (when (and (eql 0 sx) (eql 0 sy))
+          (error "~S: both arguments should not be zero"
+                 'atan
+        ) )
+        (let ((mx-my (+ (integer-length ax) ny 
+                        (- (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)))
+                )
+                ((and (minusp sy) (<= mx-my 0))
+                 (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))
+) ) ) ) ) )
+
+;;; (sin-r1 x) takes a creal x with |x|<4 and returns sin(x) as creal.
+
+(defun sin-r1 (x)
+  (declare (type creal x))
+  (make-real
+    #'(lambda (k)
+        (let ((m 3) (k2 (+ k 3)))
+          (loop (when (<= k2 (ash (- m 2) m)) (return))
+                (incf m)
+          )
+          (setq m (+ m 4) k2 (+ k m))
+          (let ((x0 (get-approx x k2)))
+            (do ((x1 (- (round-cr (* x0 x0) k2)))
+                 (n 2 (+ n 2))
+                 (y x0 (round-cr (round (* y x1) (* n (1+ n))) k2))
+                 (erg 0 (+ erg y))
+                )
+                ((eql y 0) (round-cr erg m))
+) )   ) ) ) )
+
+(defun SIN-R (x) "sine for CREALs"
+  (unless (creal-p x) (cr-error 'sin x))
+  ;; remember sin(k*2pi + y) = sin(y)
+  (if (eql x 0)
+    0
+    (multiple-value-bind (q r) (round-r x 2pi-r 10)
+      (declare (ignore q))
+      (sin-r1 r)
+) ) )
+
+;;; (cos-r1 x) takes a creal x with |x|<4 and returns cos(x) as creal.
+
+(defun cos-r1 (x)
+  (declare (type creal x))
+  (make-real
+    #'(lambda (k)
+        (let ((m 3) (k2 (+ k 3)))
+          (loop (when (<= k2 (ash (- m 2) m)) (return))
+                (incf m)
+          )
+          (setq m (+ m 4) k2 (+ k m))
+          (let ((x0 (get-approx x k2)))
+            (do ((x1 (- (round-cr (* x0 x0) k2)))
+                 (n 1 (+ n 2))
+                 (y (ash 1 k2) (round-cr (round (* y x1) (* n (1+ n))) k2))
+                 (erg 0 (+ erg y))
+                )
+                ((eql y 0) (round-cr erg m))
+) )   ) ) ) )
+
+(defun COS-R (x) "cosine for CREALs"
+  (unless (creal-p x) (cr-error 'cos x))
+  ;; remember cos(k*2pi + y) = cos(y)
+  (if (eql x 0)
+    1
+    (multiple-value-bind (q r) (round-r x 2pi-r 10)
+      (declare (ignore q))
+      (cos-r1 r)
+) ) )
+
+(defun TAN-R (x) "tangent for CREALs"
+  (unless (creal-p x) (cr-error 'tan x))
+  (/r (sin-r x) (cos-r x))
+)

File realstst.lisp

+(load "reals.lsp")
+(echo t)
+(use-package "REALS")
+(print-r pi-r 20)
+(print-r (sqrt-r 2) 20)
+(print-r pi-r 50)
+(print-r (setq e163 (exp-r (*r pi-r (sqrt-r 163)))) 20)
+(print-r (setq e58 (exp-r (*r pi-r (sqrt-r 58)))) 20)
+(defun get-koeffs (x n &aux (y x) q r)
+  (dotimes (i n)
+    (multiple-value-setq (q r) (round-r y))
+    (print q)
+    (setq y (*r x r))
+) )
+(get-koeffs e163 10)
+(get-koeffs e58 10)
+(print-r (sin-r pi-r) 20)
+(print-r (cos-r pi-r) 20)
+(print-r (sin-r pi/2-r) 20)
+(echo nil)
+