;; -*- Mode:Common-Lisp;Package:mma; Base:10 -*- (in-package :mma) ;;; Common Lisp Complex Bigfloat Package, Part I: Basic arithmetic only. ;; copyright (c) 1990 Richard Fateman, UC Berkeley ;;; Complex Bigfloats are stored as structures containing two numbers, ;;; which are interpreted as real and imaginary parts of a quantity. ;;; Often, each part is a bigfloat, but could be some other number ;;; quantity (not a complex, please.) (require "mma") (require "bf") (provide 'cbf) (proclaim '(special bigfloatone bigfloatzero)) ;; declare the structure of complex bigfloats ;; we use a rectangular form, but could presumably use (r, theta) ;; if we wished to change a bunch of routines. (defstruct (cbigfloat (:constructor make-cbig (real imag)) (:print-function cbigfloatprintfunction) ) (real :type bigfloat :read-only t) ;;real (imag :type bigfloat :read-only t) ;; imag ) (defun cbigfloatprintfunction (x s pl) (declare (ignore pl)) (format s "(~s + ~s I)" (cbigfloat-real x)(cbigfloat-imag x))) (defun cbf-real(x)(typecase x (cbigfloat (cbigfloat-real x)) (bigfloat x) (number (realpart x)) ;CL number (t `(Real ,x)))) (defun cbf-imag(x)(typecase x (cbigfloat (cbigfloat-imag x)) (bigfloat 0) (number (imagpart x)) ;CL number (t `(Imag ,x)))) ;; test for equality of two cbigfloats (defun cbigfloat-eql(x y) (and (bigfloat-eql (cbf-real x)(cbf-real y)) (bigfloat-eql (cbf-imag x)(cbf-imag y)))) ;; takes whatever x is, and returns a complex bigfloat. (defun cbigfloat-convert(x) (cond ((and (typep x 'cbigfloat) (= bigfloat-bin-prec (bigfloat-precision (cbf-real x))) (= bigfloat-bin-prec (bigfloat-precision (cbf-imag x)))) x) (t (make-cbig (bigfloat-convert (cbf-real x)) (bigfloat-convert (cbf-imag x)))))) ;; add two bigfloats of any precision; converts each to global precision. ;; c = a+b. (defun cbigfloat-+ (a b) (make-cbig (bigfloat-+ (cbf-real a)(cbf-real b)) (bigfloat-+ (cbf-imag a)(cbf-imag b)))) ;; cbigfloat-* multiplies two cbigfloats of arbitrary (perhaps different) ;; precisions, and return a cbigfloat of global (bigfloat-bin-prec) ;; precision. ;; (note: (a+bi)*(c+di)= (a*c-b*d) +(a*d+b*c)i ) (defun cbigfloat-* (r s) (let ((a (cbf-real r)) (b (cbf-imag r)) (c (cbf-real s)) (d (cbf-imag s))) (make-cbig (bigfloat-- (bigfloat-* a c)(bigfloat-* b d)) (bigfloat-+ (bigfloat-* a d)(bigfloat-* b c))))) ;; cbigfloat-/ divides two cbigfloats of arbitrary (perhaps different) ;; precisions, and return a bigfloat of global (bigfloat-bin-prec) ;; precision. ;; (a+bi)/(c+di) = 1/(c^2+d^2) * ( (ac+bd) + (bc-ad) i). (defun cbigfloat-/ (r s) (let ((a (cbf-real r)) (b (cbf-imag r)) (c (cbf-real s)) (d (cbf-imag s)) denom) (setq denom (bigfloat-+ (bigfloat-* c c)(bigfloat-* d d))) (make-cbig (bigfloat-/(bigfloat-+ (bigfloat-* a c)(bigfloat-* b d))denom) (bigfloat-/(bigfloat-- (bigfloat-* b c)(bigfloat-* a d))denom)))) ;; coerce-cbigfloat is like the CL function coerce, but allows ;; the first argument to be of type bigfloat. It converts the ;; first argument to ratio, double-float, single-float (= float), ;; integer (= bignum or fixnum) (defun coerce-cbigfloat(x typ) (cond ;; convert bigfloat to bigfloat of global precision ((eq typ 'cbigfloat) (cbigfloat-convert x)) (complex (coerce-bigfloat (cbf-real x) typ) (coerce-bigfloat (cbf-imag x) typ))) ;; compute x-y or -x (if y is missing) (defun cbigfloat-- (x &optional (y nil)) (cond (y (cbigfloat-+ x (cbigfloat-- y)));; compute x-y ((cbigfloat-zerop x) x) (t (makecbig (bigfloat-- (cbf-real x)) (bigfloat-- (cbf-imag x)))))) ;; compute p^n for bigfloat p. nn is positive or negative integer. ;; Note that we do NOT allow bigfloat exponent here. (use log/exp for that) (defun cbigfloat-expt (p nn) (format t "cbigfloat-expt not implemented yet") ) (defun cbigfloat-abs (x) (let ((c (cbf-real x))(d (cbf-imag x))) (bigfloat-+ (bigfloat-* c c)(bigfloat-* d d)))) ;;how to print these guys?