Enter your search terms Submit search form
	Web pdp-10.trailing-edge.com

;;; This is a -*-Lisp-*- file.
;;;
;;; **********************************************************************
;;; This code was written as part of the Spice Lisp project at
;;; Carnegie-Mellon University, and has been placed in the public domain.
;;; If you want to use this code or any part of Spice Lisp, please contact
;;; Scott Fahlman (FAHLMAN@CMUC). 
;;; **********************************************************************
;;;
;;; Functions to implement irrational functions for Spice Lisp 
;;; Written by David Adam.
;;; Maintained by Steven Handerson.
;;;
;;; The irrational functions are part of the standard Spicelisp environment.
;;;
;;; **********************************************************************
;;;

(in-package 'lisp)

(export '(expt isqrt asinh acosh atanh pi))

;;; Integer square root - (<= (expt result 2) input).
;;; Performs in (log input) time.

;;; Successively approximates the result using two bounds and their average,
;;; repeated until the bounds differ by at most 1 (for input=1, both start
;;; equal).

(defun isqrt (x)
  "Returns the integer square root; ie. (<= (expt result 2) input)."
  (if (and (integerp x) (not (minusp x)))
      (if (zerop  x) 0
	  (do* ((p () (<= (* m m) x))
		(b 1 (if p m b))
		(h x (if p h m))
		(m (ash x -1) (ash (+ b h) -1)))
	       ((<= h (1+ b)) b)))
      (error "Isqrt: ~S argument must be a nonnegative integer" x)))

    
;;; Function calculates the value of x raised to the nth power.
;;; This function calculates the successive squares of base,
;;; storing them in newbase, halving n at the same time.  If
;;; n is odd total is multiplied by base, at any given time (fix later)

;(proclaim '(inline intexp))

(defun intexp (base power)
  (cond ((minusp power)
	 (/ (intexp base (- power))))
	((and (rationalp base)
	      (not (integerp base)))  ; Could be a make-rational.
	 (/ (intexp (numerator base) power)
	    (intexp (denominator base) power)))
	((and (integerp base) (= base 2))
	 (ash 1 power))
	(t (do ((nextn (ash power -1) (ash power -1))
		(total (if (oddp power) base 1)
		       (if (oddp power) (* base total) total)))
	       ((zerop nextn) total)
	     (setq base (* base base))
	     (setq power nextn)))))

;;; This function calculates x raised to the nth power.  It separates
;;; the  cases by the type of n, for efficiency reasons, as powers can
;;; be calculated more efficiently if n is a positive integer,  Therefore,
;;; All integer values of n are calculated as positive integers, and
;;; inverted if negative.

(defun expt (x n)
  "Returns X raised to the Nth power."
  (cond ((and (rationalp x) (integerp n)) (intexp x n))
	((zerop x) (if (plusp n) x
		       (error "~A to a non-positive power ~A." x n)))
	((and (not (integerp n)) (minusp x))
	 (error "Negative number ~A to non-integral power ~A." x n))
	(t (%sp-expt x n))))

;;; The inverse hyperbolic functions will be coded by definition like
;;; the hyperbolic functions.

(defun asinh (x)
  "Returns the hyperbolic arc sine of the argument."
  (log (+ x (sqrt (+ (* x x) 1)))))

(defun acosh (x)
  "Returns the hyperbolic arc cosine of the argument."
  (if (plusp x)
      (log (+ x (sqrt (- (* x x) 1))))
      (error "~S argument out of range." x)))

(defun atanh (x)
  "Returns the hyperbolic arc tangent of the argument."
  (if (< -1 x 1)
      (* 0.5 (log (/ (1+ x) (- 1 x))))
      (error "~S argument out of range." x)))

;;this file gets loaded twice.  We unbind PI so that we won't get
;; a complaint the second time that it is already bound.
(makunbound 'pi)
(defconstant pi 3.1415926535897932384L0 "pi, as a long real.")