Trailing-Edge
-
PDP-10 Archives
-
clisp
-
clisp/upsala/arith.clisp
There are no other files named arith.clisp in the archive.
;;; 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).
;;; **********************************************************************
;;; Arithmetic functions for Spice Lisp.
;;; These functions are part of the standard Spice Lisp environment.
;;; Written by Jim Large
;;; Modified by Scott Fahlman
;;; Updated by Walter van Roggen, November 1982
;;; *******************************************************************
(in-package 'lisp)
(export '(max min conjugate phase signum numerator denominator floor ceiling
round mod rem ffloor fceiling ftruncate fround complex realpart
imagpart cis logior logxor logand logeqv lognand lognor logandc1
logandc2 logorc1 logorc2 lognot boole-clr boole-set boole-1 boole-2
boole-c1 boole-c2 boole-and boole-ior boole-xor boole-eqv boole-nand
boole-nor boole-andc1 boole-andc2 boole-orc1 boole-orc2 logtest
logbitp logcount ldb-test deposit-field lcm float-digits
float-precision rationalize))
;;; Predicates, excluding type predicates
;20; lots of things in the kernel
(defun max (number &rest more-numbers)
"Returns the greatest of its arguments."
(do ((nlist more-numbers (cdr nlist))
(result number))
((null nlist) (return result))
(declare (list nlist))
(if (> (car nlist) result) (setq result (car nlist)))))
(defun min (number &rest more-numbers)
"Returns the least of its arguments."
(do ((nlist more-numbers (cdr nlist))
(result number))
((null nlist) (return result))
(declare (list nlist))
(if (< (car nlist) result) (setq result (car nlist)))))
;;; Arithmetic Operations
;;; Doesn't bug out automatically, as a straight recursive version would.
;;; (i.e. (+ 1.0 2.0) <=> (+ 0 1.0 2.0))
;20; all in the kernel
;;; Note: INCF and DECF are defined in the MACROS file.
(defun conjugate (number)
"Returns the complex conjugate of NUMBER. For non-complex numbers, this is
an identity."
(if (typep number 'complex)
(complex (realpart number) (- (imagpart number)))
number))
;;; Note: GCD and LCM are defined in the implementation dependent NUM file.
;;; Note: Transcendental and trig functions are in a file of their own.
(defun phase (number)
"Returns the angle part of the polar representation of a complex number.
For non-complex numbers, this is 0."
(if (complexp number)
(atan (realpart number) (imagpart number))
0))
(defun cis (theta)
"Returns cos(Theta) + i sin(Theta)."
(declare (ignore theta))
(error "Complex numbers are not implemented yet."))
(defun signum (number)
"If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER)).
Currently not implemented for complex numbers."
(if (zerop number)
number
(if (rationalp number)
(if (plusp number) 1 -1)
(/ number (abs number)))))
;;; Type conversions.
;20; float is in kernel
;;; Note: RATIONAL and RATIONALIZE are in the implementation dependent
;;; NUM files
;20; RATIONAL is in the kernel. Here is RATIONALIZE.
;;; Rationalize does a rational, but it assumes that floats
;;; are only accurate to their precision, and generates a good
;;; rational aproximation of them.
(defun rationalize (x)
"Converts any non-complex number to be a rational number. Assumes
that a floating-point number is accurate only to the precision of the
floating-point representation, and may return any rational number for
which the floating-point number is the best available approximation of
its format; in doing this it attempts to keep both numerator and
denominator small."
(typecase x
(rational x)
(short-float (rationalize-float x short-float-epsilon))
(long-float (rationalize-float x long-float-epsilon))
(otherwise (error "Argument not a non complex number, ~A."
x))))
;;; Thanks to Kim Fateman, who stole this function rationalize-float
;;; from macsyma's rational. Macsyma'a rationalize was written
;;; by the legendary Gosper (rwg). Gosper is now working for Symbolics.
;;; Guy Steele said about Gosper, "He has been called the
;;; only living 17th century mathematician and is also the best
;;; pdp-10 hacker I know." So, if you can understand or debug this
;;; code you win big.
(defun rationalize-float (x &optional (eps long-float-epsilon))
(cond ((minusp x) (- (rationalize (- x))))
((zerop x) 0)
(t (let ((y ())
(a ()))
(do ((xx x (setq y (/ 1.0 (- xx (float a x)))))
(num (setq a (truncate x))
(+ (* (setq a (truncate y)) num) onum))
(den 1 (+ (* a den) oden))
(onum 1 num)
(oden 0 den))
((and (not (zerop den))
(not (> (abs (/ (- x (/ (float num x)
(float den x)))
x))
eps)))
(/ num den)))))))
(defun numerator (rational)
"Extracts the numerator of an integer or a ratio."
(typecase rational
(integer rational)
(ratio (%sp-svref rational %ratio-numerator-slot ))
(t (error "~S is not a rational number." rational))))
(defun denominator (rational)
"Extracts the denominator of an integer or a ratio."
(typecase rational
(integer 1)
(ratio (%sp-svref rational %ratio-denominator-slot))
(t (error "~S is not a rational number." rational))))
(defun floor (number &optional (divisor 1))
"Returns the greatest integer not greater than number, or number/divisor.
The second returned value is (mod number divisor)."
(multiple-value-bind (tru rem) (truncate number divisor)
(if (and (not (zerop rem))
(if (minusp divisor)
(plusp number)
(minusp number)))
(values (1- tru) (+ rem divisor))
(values tru rem))))
(defun ceiling (number &optional (divisor 1))
"Returns the smallest integer not less than number, or number/divisor.
The second returned value is the remainder."
(multiple-value-bind (tru rem) (truncate number divisor)
(if (and (not (zerop rem))
(if (minusp divisor)
(minusp number)
(plusp number)))
(values (1+ tru) (- rem divisor))
(values tru rem))))
;20; truncate is in the kernel
(defun round (number &optional (divisor 1 divp) &aux thresh)
"Rounds number (or number/divisor) to nearest integer.
The second returned value is the remainder."
(multiple-value-bind (tru rem) (truncate number divisor)
(setq thresh (/ (abs divisor) 2))
(cond ((or (> rem thresh)
(and (= rem thresh) (oddp tru)))
(if (minusp divisor)
(values (- tru 1) (+ rem divisor))
(values (+ tru 1) (- rem divisor))))
((or (< rem (- thresh))
(and (= rem (- thresh)) (oddp tru)))
(if (minusp divisor)
(values (+ tru 1) (- rem divisor))
(values (- tru 1) (+ rem divisor))))
(t (values tru rem)))))
(defun mod (number divisor)
"Returns second result of FLOOR."
(let ((rem (rem number divisor)))
(if (and (not (zerop rem))
(if (minusp divisor)
(plusp number)
(minusp number)))
(+ rem divisor)
rem)))
(defun rem (number divisor)
"Returns second result of TRUNCATE."
(multiple-value-bind (tru rem) (truncate number divisor)
(declare (ignore tru))
rem))
(defun ffloor (number &optional (divisor 1))
"Same as FLOOR, but returns first value as a float."
(multiple-value-bind (flr rem) (floor number divisor)
(values (float flr) rem)))
(defun fceiling (number &optional (divisor 1))
"Same as CEILING, but returns first value as a float."
(multiple-value-bind (cei rem) (ceiling number divisor)
(values (float cei) rem)))
(defun ftruncate (number &optional (divisor 1))
"Same as TRUNCATE, but returns first value as a float."
(multiple-value-bind (tru rem) (truncate number divisor)
(values (float tru) rem)))
(defun fround (number &optional (divisor 1))
"Same as ROUND, but returns first value as a float."
(multiple-value-bind (rou rem) (round number divisor)
(values (float rou) rem)))
;;; Component extractions on numbers
;;; FLOAT-SIGNIFICAND, FLOAT-EXPONENT, SCALE-FLOAT, FLOAT-RADIX, and
;;; FLOAT-SIGN are defined in the implementation dependent NUM file.
(defun complex (realpart &optional (imagpart 0))
"Builds a complex number. Currently, complex numbers are not supported."
(error "Complex numbers are not yet implemented.")
(if (complexp realpart)
(error "~S is a complex number." realpart))
(if (complexp imagpart)
(error "~S is a complex number." imagpart))
(let ((x (%sp-alloc-complex)))
(%sp-svset x %complex-real-part-slot realpart)
(%sp-svset x %complex-imaginary-part-slot imagpart)
x))
(defun realpart (number)
"Extracts the real part of a number."
(if (complexp number)
(%sp-svref number %complex-real-part-slot)
number))
(defun imagpart (number)
"Extracts the imaginary part of a number."
(typecase number
(complex
(%sp-svref number %complex-imaginary-part-slot))
(number (coerce 0 (type-of number)))
(t (error "~S not a number." number))))
;;; Bit-wise logical functions map integers onto integers
(defun logior (&rest integers)
"Returns the bit-wise or of its arguments. Args must be integers."
(declare (list integers))
(if integers
(do* ((result (pop integers) (boole 7 result (pop integers))))
((null integers) result))
0))
(defun logxor (&rest integers)
"Returns the bit-wise exclusive or of its arguments. Args must be integers."
(declare (list integers))
(if integers
(do* ((result (pop integers) (boole 6 result (pop integers))))
((null integers) result))
0))
(defun logand (&rest integers)
"Returns the bit-wise and of its arguments. Args must be integers."
(declare (list integers))
(if integers
(do* ((result (pop integers) (boole 1 result (pop integers))))
((null integers) result))
-1))
(defun logeqv (&rest integers)
"Returns the bit-wise equivalence of its arguments. Args must be integers."
(declare (list integers))
(if integers
(do* ((result (pop integers) (boole #o11 result (pop integers))))
((null integers) result))
-1))
(defun lognand (integer1 integer2)
"Returns the complement of the logical AND of integer1 and integer2."
(boole #o16 integer1 integer2))
(defun lognor (integer1 integer2)
"Returns the complement of the logical OR of integer1 and integer2."
(boole #o10 integer1 integer2))
(defun logandc1 (integer1 integer2)
"Returns the logical AND of (LOGNOT integer1) and integer2."
(boole 2 integer1 integer2))
(defun logandc2 (integer1 integer2)
"Returns the logical AND of integer1 and (LOGNOT integer2)."
(boole 4 integer1 integer2))
(defun logorc1 (integer1 integer2)
"Returns the logical OR of (LOGNOT integer1) and integer2."
(boole #o13 integer1 integer2))
(defun logorc2 (integer1 integer2)
"Returns the logical OR of integer1 and (LOGNOT integer2)."
(boole #o15 integer1 integer2))
(defun lognot (number)
"Returns the bit-wise logical not of integer."
(boole #o12 number 0))
;;; The boole function dispaches to any of the above operations depending on
;;; the value of a variable. Presently, legal selector values are [0..15].
;;; boole is open coded for calls with a constant selector. or with calls
;;; using any of the constants declared below.
(defconstant boole-clr 0
"Boole function op, makes BOOLE return 0.")
(defconstant boole-set #o17
"Boole function op, makes BOOLE return -1.")
(defconstant boole-1 5
"Boole function op, makes BOOLE return integer1.")
(defconstant boole-2 3
"Boole function op, makes BOOLE return integer2.")
(defconstant boole-c1 #o12
"Boole function op, makes BOOLE return complement of integer1.")
(defconstant boole-c2 #o14
"Boole function op, makes BOOLE return complement of integer2.")
(defconstant boole-and 1
"Boole function op, makes BOOLE return logand of integer1 and integer2.")
(defconstant boole-ior 7
"Boole function op, makes BOOLE return logior of integer1 and integer2.")
(defconstant boole-xor 6
"Boole function op, makes BOOLE return logxor of integer1 and integer2.")
(defconstant boole-eqv #o11
"Boole function op, makes BOOLE return logeqv of integer1 and integer2.")
(defconstant boole-nand #o16
"Boole function op, makes BOOLE return log nand of integer1 and integer2.")
(defconstant boole-nor #o10
"Boole function op, makes BOOLE return lognor of integer1 and integer2.")
(defconstant boole-andc1 2
"Boole function op, makes BOOLE return logandc1 of integer1 and integer2.")
(defconstant boole-andc2 4
"Boole function op, makes BOOLE return logandc2 of integer1 and integer2.")
(defconstant boole-orc1 #o13
"Boole function op, makes BOOLE return logorc1 of integer1 and integer2.")
(defconstant boole-orc2 #o15
"Boole function op, makes BOOLE return logorc2 of integer1 and integer2.")
;;; Bit testing, Shifting, and assorted logcal operations
;20; not in kernel. for the moment write code for this
(defun logtest (integer1 integer2)
"Predicate which returns T if logand of integer1 and integer2 is not zero."
(not (zerop (boole 1 integer1 integer2))))
(defun logbitp (index integer)
"Predicate returns T if bit index of integer is a 1."
(ldb-test (byte 1 index) integer))
(defun logcount (integer)
"If INTEGER is negative, then # of 0 bits is returned,
else # of 1 bits is returned."
(let ((total 0))
(if (minusp integer)
(dotimes (i (integer-length integer))
(if (not (logbitp i integer)) (incf total)))
(dotimes (i (integer-length integer))
(if (logbitp i integer) (incf total))))
total))
;20; integer-length is in the kernel
;;; Byte manipulation functions operate on fields within integers.
;20; most of this is in the kernel
(defun ldb-test (bytespec integer)
"Returns T if any of the specified bits in integer are 1's."
(not (zerop (ldb bytespec integer))))
(defun deposit-field (newbyte bytespec integer)
"Returns new integer with newbyte in specified position, newbyte is not
right justified."
(dpb (ldb bytespec newbyte) bytespec integer))
;;; Lcm.
;;; Lcm2 is defined this way so that operations won't unnecessarily bignumify.
(defun lcm2 (n m)
"Least common multiple of two nonzero integers. No type checking."
(* (/ (max n m) (gcd n m)) (min n m)))
(defun lcm (&rest args)
"Returns the least common multiple of one or more integers."
(do ((args args (cdr args))
(lcm 1 (lcm2 lcm (car args))))
((null args) lcm)
(cond ((not (integerp (car args)))
(error "Lcm: argument not an integer, ~A"
(car args)))
((zerop (car args)) ;Result is zero.
(dolist (arg (cdr args))
(if (not (integerp arg))
(error "Lcm: argument not an integer, ~A" arg)))
(return 0)))))
;;; FLOAT-RADIX, FLOAT-SIGN, etc are in the kernel.
(defun float-digits (f)
"Returns, as a non-negative integer, the number of radix-b digits used in
the representation of its argument."
(typecase f
(short-float 23)
(long-float 62)
(t (error "Float-digits: ~A not a float" f))))
(defun float-precision (f)
"Returns, as a non-negative integer, the number of significant radix-b
digits present in the argument: if the argument is (a floating-point) zero,
then the result is (an integer) zero."
(if (zerop f)
0
(float-digits f)))