;;; numbers.lisp ;;; ;;; Copyright (C) 2003-2006 Peter Graves ;;; $Id$ ;;; ;;; This program is free software; you can redistribute it and/or ;;; modify it under the terms of the GNU General Public License ;;; as published by the Free Software Foundation; either version 2 ;;; of the License, or (at your option) any later version. ;;; ;;; This program is distributed in the hope that it will be useful, ;;; but WITHOUT ANY WARRANTY; without even the implied warranty of ;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ;;; GNU General Public License for more details. ;;; ;;; You should have received a copy of the GNU General Public License ;;; along with this program; if not, write to the Free Software ;;; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. ;;; ;;; As a special exception, the copyright holders of this library give you ;;; permission to link this library with independent modules to produce an ;;; executable, regardless of the license terms of these independent ;;; modules, and to copy and distribute the resulting executable under ;;; terms of your choice, provided that you also meet, for each linked ;;; independent module, the terms and conditions of the license of that ;;; module. An independent module is a module which is not derived from ;;; or based on this library. If you modify this library, you may extend ;;; this exception to your version of the library, but you are not ;;; obligated to do so. If you do not wish to do so, delete this ;;; exception statement from your version. ;;; Adapted from CMUCL/SBCL. (in-package "SYSTEM") (defun signum (number) "If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))." (if (zerop number) number (if (rationalp number) (if (plusp number) 1 -1) (/ number (abs number))))) (defun round (number &optional (divisor 1)) "Rounds number (or number/divisor) to nearest integer. The second returned value is the remainder." (multiple-value-bind (tru rem) (truncate number divisor) (if (zerop rem) (values tru rem) (let ((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)))) ((let ((-thresh (- thresh))) (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 ffloor (number &optional (divisor 1)) "Same as FLOOR, but returns first value as a float." (multiple-value-bind (tru rem) (ftruncate number divisor) (if (and (not (zerop rem)) (if (minusp divisor) (plusp number) (minusp number))) (values (1- tru) (+ rem divisor)) (values tru rem)))) (defun fceiling (number &optional (divisor 1)) "Same as CEILING, but returns first value as a float." (multiple-value-bind (tru rem) (ftruncate number divisor) (if (and (not (zerop rem)) (if (minusp divisor) (minusp number) (plusp number))) (values (+ tru 1) (- rem divisor)) (values tru rem)))) (defun fround (number &optional (divisor 1)) "Same as ROUND, but returns first value as a float." (multiple-value-bind (res rem) (round number divisor) (values (float res (if (floatp rem) rem 1.0)) rem))) ;;; FIXME (defun rationalize (number) (rational number)) (defun gcd (&rest integers) (cond ((null integers) 0) ((null (cdr integers)) (let ((n (car integers))) (if (integerp n) (abs n) (error 'type-error :datum n :expected-type 'integer)))) (t (do ((gcd (car integers) (gcd-2 gcd (car rest))) (rest (cdr integers) (cdr rest))) ((null rest) gcd))))) ;;; From discussion on comp.lang.lisp and Akira Kurihara. (defun isqrt (natural) "Returns the root of the nearest integer less than natural which is a perfect square." (unless (and (integerp natural) (not (minusp natural))) (error 'simple-type-error :format-control "The value ~A is not a non-negative real number." :format-arguments (list natural))) (if (and (fixnump natural) (<= natural 24)) (cond ((> natural 15) 4) ((> natural 8) 3) ((> natural 3) 2) ((> natural 0) 1) (t 0)) (let* ((n-len-quarter (ash (integer-length natural) -2)) (n-half (ash natural (- (ash n-len-quarter 1)))) (n-half-isqrt (isqrt n-half)) (init-value (ash (1+ n-half-isqrt) n-len-quarter))) (loop (let ((iterated-value (ash (+ init-value (truncate natural init-value)) -1))) (unless (< iterated-value init-value) (return init-value)) (setq init-value iterated-value)))))) ;; FIXME Need to add support for denormalized floats! ;; "FLOAT-PRECISION returns the number of significant radix b digits present in ;; FLOAT; if FLOAT is a float zero, then the result is an integer zero." ;; "For normalized floats, the results of FLOAT-DIGITS and FLOAT-PRECISION are ;; the same, but the precision is less than the number of representation digits ;; for a denormalized or zero number. (defun float-precision (float) (if (floatp float) (cond ((zerop float) 0) ((typep float 'single-float) 24) ((typep float 'double-float) 53) (t ;; Shouldn't get here! (aver nil))) (error 'simple-type-error :format-control "~S is not of type FLOAT." :format-arguments (list float)))) (defun decode-float (float) (cond ((typep float 'single-float) (decode-float-single float)) ((typep float 'double-float) (decode-float-double float)) (t (error 'simple-type-error :format-control "~S is neither SINGLE-FLOAT nor DOUBLE-FLOAT." :format-arguments (list float))))) ;;; From . Thanks Xophe! (defun sane-integer-decode-float (float) (multiple-value-bind (mantissa exp sign) (integer-decode-float float) (let ((fixup (- (integer-length mantissa) (float-precision float)))) (values (ash mantissa (- fixup)) (+ exp fixup) sign)))) (defun decode-float-single (float) ;; TODO memoize (let ((float-precision-single (float-precision 1f0))) (multiple-value-bind (significand exponent sign) (sane-integer-decode-float float) (values (coerce (/ significand (expt 2 float-precision-single)) 'single-float) (+ exponent float-precision-single) (if (minusp sign) -1f0 1f0))))) (defun decode-float-double (float) ;; TODO memoize (let ((float-precision-double (float-precision 1d0))) (multiple-value-bind (significand exponent sign) (sane-integer-decode-float float) (values (coerce (/ significand (expt 2 float-precision-double)) 'double-float) (+ exponent float-precision-double) (if (minusp sign) -1d0 1d0))))) (defun conjugate (number) (etypecase number (complex (complex (realpart number) (- (imagpart number)))) (number number))) (defun phase (number) "Returns the angle part of the polar representation of a complex number. For complex numbers, this is (atan (imagpart number) (realpart number)). For non-complex positive numbers, this is 0. For non-complex negative numbers this is PI." (etypecase number (rational (if (minusp number) (coerce pi 'single-float) 0.0f0)) (single-float (if (minusp (float-sign number)) (coerce pi 'single-float) 0.0f0)) (double-float (if (minusp (float-sign number)) (coerce pi 'double-float) 0.0d0)) (complex (if (zerop (realpart number)) (coerce (* (/ pi 2) (signum (imagpart number))) (if (typep (imagpart number) 'double-float) 'double-float 'single-float)) (atan (imagpart number) (realpart number))))))