;;; make-real creates a c-real from a computation function.

+ "Create a c-real from a computation function COMP."

(declare (type (function ((integer 0 *)) integer) comp))

(make-c-real :compute comp))

"Computes approximations for CREALs"

- (unless (creal-p x) (cr-error 'approx-r x))

- (unless (and (integerp k) (>= k 0)) (nat-error 'approx-r k))

+ (assert (and (and integerp k)

- (declare (type creal x) (type (integer 0 *) k))

+ (declare (type creal x)

+ (type (integer 0 *) k))

(cond ((integerp x) (ash x k))

((rationalp x) (round (ash (numerator x) k) (denominator x)))

"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))

+ (assert (and (integerp k) (not (minusp k))))

+ (assert (streamp stream))

(creal-print x k flag stream))

(defun creal-print (x k flag stream)

- (~~unless~~ (creal-p x)~~ (cr-error 'sqrt x)~~)

(if (and (rationalp x) (>= x 0) (rationalp (setq s (sqrt x))))

(multiple-value-bind (a0 n0 s) (raw-approx-cr x)

(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)

"shift function for CREALs"

- (unless (creal-p x) (cr-error 'ash-r x))

- (unless (integerp n) (int-error 'ash-r n))

(if (plusp n) (ash x n) (/ x (ash 1 (- n)))))

(defun LOG-R (x &optional (b nil))

- (unless (creal-p x) (cr-error 'log x))

- (unless (or (null b) (creal-p b)) (cr-error 'log b))

+ (assert (or (null b) (creal-p 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)

- (error "~S: attempt to compute the logarithm of a nonpositive number"

- (let ((shift (- (integer-length a0) 1 n0)))

- (rest-help-r (log-r1 (ash-r x (- shift))) log2-r shift)))))

+ (/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)

+ (error "~S: attempt to compute the logarithm of a nonpositive number"

+ (let ((shift (- (integer-length a0) 1 n0)))

+ (rest-help-r (log-r1 (ash-r x (- shift))) log2-r shift)))))

;;; Now the exponential function.

(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))

(if (rationalp x) (expt x y) (expt-r1 x y)))

(defun ATAN-R (x &optional (y nil))

- (unless (creal-p x) (cr-error 'atan x))

- (unless (or (null y) (creal-p y)) (cr-error 'atan y))

+ (assert (or (null y) (creal-p y)))

(multiple-value-bind (ax nx sx) (raw-approx-cr x)

- (~~unless~~ (creal-p x)~~ (cr-error 'sin x)~~)

;; remember sin(k*2pi + y) = sin(y)

- (~~unless~~ (creal-p x)~~ (cr-error 'cos x)~~)

;; remember cos(k*2pi + y) = cos(y)

- (~~unless~~ (creal-p x)~~ (cr-error 'tan x)~~)

(/r (sin-r x) (cos-r x)))