Annotation of GNUtools/emacs/lisp/float.el, revision 1.1

1.1     ! root        1: ;; Copyright (C) 1986 Free Software Foundation, Inc.
        !             2: ;; Author Bill Rosenblatt
        !             3: 
        !             4: ;; This file is part of GNU Emacs.
        !             5: 
        !             6: ;; GNU Emacs is free software; you can redistribute it and/or modify
        !             7: ;; it under the terms of the GNU General Public License as published by
        !             8: ;; the Free Software Foundation; either version 1, or (at your option)
        !             9: ;; any later version.
        !            10: 
        !            11: ;; GNU Emacs is distributed in the hope that it will be useful,
        !            12: ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
        !            13: ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
        !            14: ;; GNU General Public License for more details.
        !            15: 
        !            16: ;; You should have received a copy of the GNU General Public License
        !            17: ;; along with GNU Emacs; see the file COPYING.  If not, write to
        !            18: ;; the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
        !            19: 
        !            20: ;; Floating point arithmetic package.
        !            21: ;;
        !            22: ;; Floating point numbers are represented by dot-pairs (mant . exp)
        !            23: ;; where mant is the 24-bit signed integral mantissa and exp is the
        !            24: ;; base 2 exponent.
        !            25: ;;
        !            26: ;; Emacs LISP supports a 24-bit signed integer data type, which has a
        !            27: ;; range of -(2**23) to +(2**23)-1, or -8388608 to 8388607 decimal.
        !            28: ;; This gives six significant decimal digit accuracy.  Exponents can
        !            29: ;; be anything in the range -(2**23) to +(2**23)-1.
        !            30: ;;
        !            31: ;; User interface:
        !            32: ;; function f converts from integer to floating point
        !            33: ;; function string-to-float converts from string to floating point
        !            34: ;; function fint converts a floating point to integer (with truncation)
        !            35: ;; function float-to-string converts from floating point to string
        !            36: ;;                   
        !            37: ;; Caveats:
        !            38: ;; -  Exponents outside of the range of +/-100 or so will cause certain 
        !            39: ;;    functions (especially conversion routines) to take forever.
        !            40: ;; -  Very little checking is done for fixed point overflow/underflow.
        !            41: ;; -  No checking is done for over/underflow of the exponent
        !            42: ;;    (hardly necessary when exponent can be 2**23).
        !            43: ;; 
        !            44: ;;
        !            45: ;; Bill Rosenblatt
        !            46: ;; June 20, 1986
        !            47: ;;
        !            48: 
        !            49: (provide 'float)
        !            50: 
        !            51: ;; fundamental implementation constants
        !            52: (defconst exp-base 2
        !            53:   "Base of exponent in this floating point representation.")
        !            54: 
        !            55: (defconst mantissa-bits 24
        !            56:   "Number of significant bits in this floating point representation.")
        !            57: 
        !            58: (defconst decimal-digits 6
        !            59:   "Number of decimal digits expected to be accurate.")
        !            60: 
        !            61: (defconst expt-digits 2
        !            62:   "Maximum permitted digits in a scientific notation exponent.")
        !            63: 
        !            64: ;; other constants
        !            65: (defconst maxbit (1- mantissa-bits)
        !            66:   "Number of highest bit")
        !            67: 
        !            68: (defconst mantissa-maxval (1- (ash 1 maxbit))
        !            69:   "Maximum permissable value of mantissa")
        !            70: 
        !            71: ;;; Note that this value can't be plain (ash 1 maxbit), since
        !            72: ;;; (- (ash 1 maxbit)) = (ash 1 maxbit) - it overflows.
        !            73: (defconst mantissa-minval (1- (ash 1 maxbit))
        !            74:   "Minimum permissable value of mantissa")
        !            75: 
        !            76: ;;; This is used when normalizing negative numbers; if the number is
        !            77: ;;; less than this, multiplying it by 2 will overflow past
        !            78: ;;; mantissa-minval.
        !            79: (defconst mantissa-half-minval (ash (ash 1 maxbit) -1))
        !            80: 
        !            81: (defconst floating-point-regexp
        !            82:   "^[ \t]*\\(-?\\)\\([0-9]*\\)\
        !            83: \\(\\.\\([0-9]*\\)\\|\\)\
        !            84: \\(\\(\\([Ee]\\)\\(-?\\)\\([0-9][0-9]*\\)\\)\\|\\)[ \t]*$"
        !            85:   "Regular expression to match floating point numbers.  Extract matches:
        !            86: 1 - minus sign
        !            87: 2 - integer part
        !            88: 4 - fractional part
        !            89: 8 - minus sign for power of ten
        !            90: 9 - power of ten
        !            91: ")
        !            92: 
        !            93: (defconst high-bit-mask (ash 1 maxbit)
        !            94:   "Masks all bits except the high-order (sign) bit.")
        !            95: 
        !            96: (defconst second-bit-mask (ash 1 (1- maxbit))
        !            97:   "Masks all bits except the highest-order magnitude bit")
        !            98: 
        !            99: ;; various useful floating point constants
        !           100: (setq _f0 '(0 . 1))
        !           101: 
        !           102: (setq _f1/2 '(4194304 . -23))
        !           103: 
        !           104: (setq _f1 '(4194304 . -22))
        !           105: 
        !           106: (setq _f10 '(5242880 . -19))
        !           107: 
        !           108: ;; support for decimal conversion routines
        !           109: (setq powers-of-10 (make-vector (1+ decimal-digits) _f1))
        !           110: (aset powers-of-10 1 _f10)
        !           111: (aset powers-of-10 2 '(6553600 . -16))
        !           112: (aset powers-of-10 3 '(8192000 . -13))
        !           113: (aset powers-of-10 4 '(5120000 . -9))
        !           114: (aset powers-of-10 5 '(6400000 . -6))
        !           115: (aset powers-of-10 6 '(8000000 . -3))
        !           116: 
        !           117: (setq all-decimal-digs-minval (aref powers-of-10 (1- decimal-digits))
        !           118:       highest-power-of-10 (aref powers-of-10 decimal-digits))
        !           119: 
        !           120: (defun fashl (fnum)                    ; floating-point arithmetic shift left
        !           121:   (cons (ash (car fnum) 1) (1- (cdr fnum))))
        !           122: 
        !           123: (defun fashr (fnum)                    ; floating point arithmetic shift right
        !           124:   (cons (ash (car fnum) -1) (1+ (cdr fnum))))
        !           125: 
        !           126: (defun normalize (fnum)
        !           127:   (if (> (car fnum) 0)                 ; make sure next-to-highest bit is set
        !           128:       (while (zerop (logand (car fnum) second-bit-mask))
        !           129:        (setq fnum (fashl fnum)))
        !           130:     (if (< (car fnum) 0)               ; make sure next-to-highest bit is
        !           131:                                        ; zero, but fnum /= mantissa-minval.
        !           132:        (while (> (car fnum) mantissa-half-minval)
        !           133:          (setq fnum (fashl fnum)))
        !           134:       (setq fnum _f0)))                        ; "standard 0"
        !           135:   fnum)
        !           136:       
        !           137: (defun abs (n)                         ; integer absolute value
        !           138:   (if (natnump n) n (- n)))
        !           139: 
        !           140: (defun fabs (fnum)                     ; re-normalize after taking abs value
        !           141:   (normalize (cons (abs (car fnum)) (cdr fnum))))
        !           142: 
        !           143: (defun xor (a b)                       ; logical exclusive or
        !           144:   (and (or a b) (not (and a b))))
        !           145: 
        !           146: (defun same-sign (a b)                 ; two f-p numbers have same sign?
        !           147:   (not (xor (natnump (car a)) (natnump (car b)))))
        !           148: 
        !           149: (defun extract-match (str i)           ; used after string-match
        !           150:   (condition-case ()
        !           151:       (substring str (match-beginning i) (match-end i))
        !           152:     (error "")))
        !           153: 
        !           154: ;; support for the multiplication function
        !           155: (setq halfword-bits (/ mantissa-bits 2)        ; bits in a halfword
        !           156:       masklo (1- (ash 1 halfword-bits)) ; isolate the lower halfword
        !           157:       maskhi (lognot masklo)           ; isolate the upper halfword
        !           158:       round-limit (ash 1 (/ halfword-bits 2)))
        !           159: 
        !           160: (defun hihalf (n)                      ; return high halfword, shifted down
        !           161:   (ash (logand n maskhi) (- halfword-bits)))
        !           162: 
        !           163: (defun lohalf (n)                      ; return low halfword
        !           164:   (logand n masklo))
        !           165: 
        !           166: ;; Visible functions
        !           167: 
        !           168: ;; Arithmetic functions
        !           169: (defun f+ (a1 a2)
        !           170:   "Returns the sum of two floating point numbers."
        !           171:   (let ((f1 (if (> (cdr a1) (cdr a2)) a1 a2))
        !           172:        (f2 (if (> (cdr a1) (cdr a2)) a2 a1)))
        !           173:     (if (same-sign a1 a2)
        !           174:        (setq f1 (fashr f1)             ; shift right to avoid overflow
        !           175:              f2 (fashr f2)))
        !           176:     (normalize
        !           177:      (cons (+ (car f1) (ash (car f2) (- (cdr f2) (cdr f1))))
        !           178:           (cdr f1)))))
        !           179: 
        !           180: (defun f- (a1 &optional a2)            ; unary or binary minus
        !           181:   "Returns the difference of two floating point numbers."
        !           182:   (if a2
        !           183:       (f+ a1 (f- a2))
        !           184:     (normalize (cons (- (car a1)) (cdr a1)))))
        !           185: 
        !           186: (defun f* (a1 a2)                      ; multiply in halfword chunks
        !           187:   "Returns the product of two floating point numbers."
        !           188:   (let* ((i1 (car (fabs a1)))
        !           189:         (i2 (car (fabs a2)))
        !           190:         (sign (not (same-sign a1 a2)))
        !           191:         (prodlo (+ (hihalf (* (lohalf i1) (lohalf i2)))
        !           192:                    (lohalf (* (hihalf i1) (lohalf i2)))
        !           193:                    (lohalf (* (lohalf i1) (hihalf i2)))))
        !           194:         (prodhi (+ (* (hihalf i1) (hihalf i2))
        !           195:                    (hihalf (* (hihalf i1) (lohalf i2)))
        !           196:                    (hihalf (* (lohalf i1) (hihalf i2)))
        !           197:                    (hihalf prodlo))))
        !           198:     (if (> (lohalf prodlo) round-limit)
        !           199:        (setq prodhi (1+ prodhi)))      ; round off truncated bits
        !           200:     (normalize
        !           201:      (cons (if sign (- prodhi) prodhi)
        !           202:           (+ (cdr (fabs a1)) (cdr (fabs a2)) mantissa-bits)))))
        !           203: 
        !           204: (defun f/ (a1 a2)                      ; SLOW subtract-and-shift algorithm
        !           205:   "Returns the quotient of two floating point numbers."
        !           206:   (if (zerop (car a2))                 ; if divide by 0
        !           207:       (signal 'arith-error (list "attempt to divide by zero" a1 a2))
        !           208:     (let ((bits (1- maxbit))
        !           209:          (quotient 0) 
        !           210:          (dividend (car (fabs a1)))
        !           211:          (divisor (car (fabs a2)))
        !           212:          (sign (not (same-sign a1 a2))))
        !           213:       (while (natnump bits)
        !           214:        (if (< (- dividend divisor) 0)
        !           215:            (setq quotient (ash quotient 1))
        !           216:          (setq quotient (1+ (ash quotient 1))
        !           217:                dividend (- dividend divisor)))
        !           218:        (setq dividend (ash dividend 1)
        !           219:              bits (1- bits)))
        !           220:       (normalize
        !           221:        (cons (if sign (- quotient) quotient)
        !           222:             (- (cdr (fabs a1)) (cdr (fabs a2)) (1- maxbit)))))))
        !           223:   
        !           224: (defun f% (a1 a2)
        !           225:   "Returns the remainder of first floating point number divided by second."
        !           226:   (f- a1 (f* (ftrunc (f/ a1 a2)) a2)))
        !           227:          
        !           228: 
        !           229: ;; Comparison functions
        !           230: (defun f= (a1 a2)
        !           231:   "Returns t if two floating point numbers are equal, nil otherwise."
        !           232:   (equal a1 a2))
        !           233: 
        !           234: (defun f> (a1 a2)
        !           235:   "Returns t if first floating point number is greater than second,
        !           236: nil otherwise."
        !           237:   (cond ((and (natnump (car a1)) (< (car a2) 0)) 
        !           238:         t)                             ; a1 nonnegative, a2 negative
        !           239:        ((and (> (car a1) 0) (<= (car a2) 0))
        !           240:         t)                             ; a1 positive, a2 nonpositive
        !           241:        ((and (<= (car a1) 0) (natnump (car a2)))
        !           242:         nil)                           ; a1 nonpos, a2 nonneg
        !           243:        ((/= (cdr a1) (cdr a2))         ; same signs.  exponents differ
        !           244:         (> (cdr a1) (cdr a2)))         ; compare the mantissas.
        !           245:        (t
        !           246:         (> (car a1) (car a2)))))       ; same exponents.
        !           247: 
        !           248: (defun f>= (a1 a2)
        !           249:   "Returns t if first floating point number is greater than or equal to 
        !           250: second, nil otherwise."
        !           251:   (or (f> a1 a2) (f= a1 a2)))
        !           252: 
        !           253: (defun f< (a1 a2)
        !           254:   "Returns t if first floating point number is less than second,
        !           255: nil otherwise."
        !           256:   (not (f>= a1 a2)))
        !           257: 
        !           258: (defun f<= (a1 a2)
        !           259:   "Returns t if first floating point number is less than or equal to
        !           260: second, nil otherwise."
        !           261:   (not (f> a1 a2)))
        !           262: 
        !           263: (defun f/= (a1 a2)
        !           264:   "Returns t if first floating point number is not equal to second,
        !           265: nil otherwise."
        !           266:   (not (f= a1 a2)))
        !           267: 
        !           268: (defun fmin (a1 a2)
        !           269:   "Returns the minimum of two floating point numbers."
        !           270:   (if (f< a1 a2) a1 a2))
        !           271: 
        !           272: (defun fmax (a1 a2)
        !           273:   "Returns the maximum of two floating point numbers."
        !           274:   (if (f> a1 a2) a1 a2))
        !           275:       
        !           276: (defun fzerop (fnum)
        !           277:   "Returns t if the floating point number is zero, nil otherwise."
        !           278:   (= (car fnum) 0))
        !           279: 
        !           280: (defun floatp (fnum)
        !           281:   "Returns t if the arg is a floating point number, nil otherwise."
        !           282:   (and (consp fnum) (integerp (car fnum)) (integerp (cdr fnum))))
        !           283: 
        !           284: ;; Conversion routines
        !           285: (defun f (int)
        !           286:   "Convert the integer argument to floating point, like a C cast operator."
        !           287:   (normalize (cons int '0)))
        !           288: 
        !           289: (defun int-to-hex-string (int)
        !           290:   "Convert the integer argument to a C-style hexadecimal string."
        !           291:   (let ((shiftval -20)
        !           292:        (str "0x")
        !           293:        (hex-chars "0123456789ABCDEF"))
        !           294:     (while (<= shiftval 0)
        !           295:       (setq str (concat str (char-to-string 
        !           296:                        (aref hex-chars
        !           297:                              (logand (lsh int shiftval) 15))))
        !           298:            shiftval (+ shiftval 4)))
        !           299:     str))
        !           300: 
        !           301: (defun ftrunc (fnum)                   ; truncate fractional part
        !           302:   "Truncate the fractional part of a floating point number."
        !           303:   (cond ((natnump (cdr fnum))          ; it's all integer, return number as is
        !           304:         fnum)
        !           305:        ((<= (cdr fnum) (- maxbit))     ; it's all fractional, return 0
        !           306:         '(0 . 1))
        !           307:        (t                              ; otherwise mask out fractional bits
        !           308:         (let ((mant (car fnum)) (exp (cdr fnum)))
        !           309:           (normalize 
        !           310:            (cons (if (natnump mant)    ; if negative, use absolute value
        !           311:                      (ash (ash mant exp) (- exp))
        !           312:                    (- (ash (ash (- mant) exp) (- exp))))
        !           313:                  exp))))))
        !           314: 
        !           315: (defun fint (fnum)                     ; truncate and convert to integer
        !           316:   "Convert the floating point number to integer, with truncation, 
        !           317: like a C cast operator."
        !           318:   (let* ((tf (ftrunc fnum)) (tint (car tf)) (texp (cdr tf)))
        !           319:     (cond ((>= texp mantissa-bits)     ; too high, return "maxint"
        !           320:           mantissa-maxval)
        !           321:          ((<= texp (- mantissa-bits))  ; too low, return "minint"
        !           322:           mantissa-minval)
        !           323:          (t                            ; in range
        !           324:           (ash tint texp)))))          ; shift so that exponent is 0
        !           325: 
        !           326: (defun float-to-string (fnum &optional sci)
        !           327:   "Convert the floating point number to a decimal string.
        !           328: Optional second argument non-nil means use scientific notation."
        !           329:   (let* ((value (fabs fnum)) (sign (< (car fnum) 0))
        !           330:         (power 0) (result 0) (str "") 
        !           331:         (temp 0) (pow10 _f1))
        !           332: 
        !           333:     (if (f= fnum _f0)
        !           334:        "0"
        !           335:       (if (f>= value _f1)                      ; find largest power of 10 <= value
        !           336:          (progn                                ; value >= 1, power is positive
        !           337:            (while (f<= (setq temp (f* pow10 highest-power-of-10)) value)
        !           338:              (setq pow10 temp
        !           339:                    power (+ power decimal-digits)))
        !           340:            (while (f<= (setq temp (f* pow10 _f10)) value)
        !           341:              (setq pow10 temp
        !           342:                    power (1+ power))))
        !           343:        (progn                          ; value < 1, power is negative
        !           344:          (while (f> (setq temp (f/ pow10 highest-power-of-10)) value)
        !           345:            (setq pow10 temp
        !           346:                  power (- power decimal-digits)))
        !           347:          (while (f> pow10 value)
        !           348:            (setq pow10 (f/ pow10 _f10)
        !           349:                  power (1- power)))))
        !           350:                                          ; get value in range 100000 to 999999
        !           351:       (setq value (f* (f/ value pow10) all-decimal-digs-minval)
        !           352:            result (ftrunc value))
        !           353:       (let (int)
        !           354:        (if (f> (f- value result) _f1/2)        ; round up if remainder > 0.5
        !           355:            (setq int (1+ (fint result)))
        !           356:          (setq int (fint result)))
        !           357:        (setq str (int-to-string int))
        !           358:        (if (>= int 1000000)
        !           359:            (setq power (1+ power))))
        !           360: 
        !           361:       (if sci                          ; scientific notation
        !           362:          (setq str (concat (substring str 0 1) "." (substring str 1)
        !           363:                            "E" (int-to-string power)))
        !           364: 
        !           365:                                          ; regular decimal string
        !           366:        (cond ((>= power (1- decimal-digits))
        !           367:                                          ; large power, append zeroes
        !           368:               (let ((zeroes (- power decimal-digits)))
        !           369:                 (while (natnump zeroes)
        !           370:                   (setq str (concat str "0")
        !           371:                         zeroes (1- zeroes)))))
        !           372: 
        !           373:                                          ; negative power, prepend decimal
        !           374:              ((< power 0)              ; point and zeroes
        !           375:               (let ((zeroes (- (- power) 2)))
        !           376:                 (while (natnump zeroes)
        !           377:                   (setq str (concat "0" str)
        !           378:                         zeroes (1- zeroes)))
        !           379:                 (setq str (concat "0." str))))
        !           380: 
        !           381:              (t                                ; in range, insert decimal point
        !           382:               (setq str (concat
        !           383:                          (substring str 0 (1+ power))
        !           384:                          "."
        !           385:                          (substring str (1+ power)))))))
        !           386: 
        !           387:       (if sign                         ; if negative, prepend minus sign
        !           388:          (concat "-" str)
        !           389:        str))))
        !           390: 
        !           391:     
        !           392: ;; string to float conversion.
        !           393: ;; accepts scientific notation, but ignores anything after the first two
        !           394: ;; digits of the exponent.
        !           395: (defun string-to-float (str)
        !           396:   "Convert the string to a floating point number.
        !           397: Accepts a decimal string in scientific notation, 
        !           398: with exponent preceded by either E or e.
        !           399: Only the 6 most significant digits of the integer and fractional parts
        !           400: are used; only the first two digits of the exponent are used.
        !           401: Negative signs preceding both the decimal number and the exponent
        !           402: are recognized."
        !           403: 
        !           404:   (if (string-match floating-point-regexp str 0)
        !           405:       (let (power)
        !           406:        (f*
        !           407:         ; calculate the mantissa
        !           408:         (let* ((int-subst (extract-match str 2))
        !           409:                (fract-subst (extract-match str 4))
        !           410:                (digit-string (concat int-subst fract-subst))
        !           411:                (mant-sign (equal (extract-match str 1) "-"))
        !           412:                (leading-0s 0) (round-up nil))
        !           413: 
        !           414:           ; get rid of leading 0's
        !           415:           (setq power (- (length int-subst) decimal-digits))
        !           416:           (while (and (< leading-0s (length digit-string))
        !           417:                       (= (aref digit-string leading-0s) ?0))
        !           418:             (setq leading-0s (1+ leading-0s)))
        !           419:           (setq power (- power leading-0s)
        !           420:                 digit-string (substring digit-string leading-0s))
        !           421:           
        !           422:           ; if more than 6 digits, round off
        !           423:           (if (> (length digit-string) decimal-digits)
        !           424:               (setq round-up (>= (aref digit-string decimal-digits) ?5)
        !           425:                     digit-string (substring digit-string 0 decimal-digits))
        !           426:             (setq power (+ power (- decimal-digits (length digit-string)))))
        !           427: 
        !           428:           ; round up and add minus sign, if necessary
        !           429:           (f (* (+ (string-to-int digit-string)
        !           430:                    (if round-up 1 0))
        !           431:                 (if mant-sign -1 1))))
        !           432:           
        !           433:         ; calculate the exponent (power of ten)
        !           434:         (let* ((expt-subst (extract-match str 9))
        !           435:                (expt-sign (equal (extract-match str 8) "-"))
        !           436:                (expt 0) (chunks 0) (tens 0) (exponent _f1)
        !           437:                (func 'f*))
        !           438:  
        !           439:           (setq expt (+ (* (string-to-int
        !           440:                             (substring expt-subst 0
        !           441:                                        (min expt-digits (length expt-subst))))
        !           442:                            (if expt-sign -1 1))
        !           443:                         power))
        !           444:           (if (< expt 0)               ; if power of 10 negative
        !           445:               (setq expt (- expt)      ; take abs val of exponent
        !           446:                     func 'f/))         ; and set up to divide, not multiply
        !           447: 
        !           448:           (setq chunks (/ expt decimal-digits)
        !           449:                 tens (% expt decimal-digits))
        !           450:           ; divide or multiply by "chunks" of 10**6
        !           451:           (while (> chunks 0)  
        !           452:             (setq exponent (funcall func exponent highest-power-of-10)
        !           453:                   chunks (1- chunks)))
        !           454:           ; divide or multiply by remaining power of ten
        !           455:           (funcall func exponent (aref powers-of-10 tens)))))
        !           456:                  
        !           457:     _f0))                              ; if invalid, return 0
        !           458: 
        !           459: 

unix.superglobalmegacorp.com

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.