home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / float.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  38.2 KB  |  1,099 lines
  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8. ;;;     (c) Copyright 1982 Massachusetts Institute of Technology         ;;;
  9. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10.  
  11. (in-package "MAXIMA")
  12. (macsyma-module float)
  13.  
  14. ;; EXPERIMENTAL BIGFLOAT PACKAGE VERSION 2- USING BINARY MANTISSA 
  15. ;; AND POWER-OF-2 EXPONENT.  EXPONENTS MAY BE BIG NUMBERS NOW (AUG. 1975 --RJF)
  16. ;; Modified:    July 1979 by CWH to run on the Lisp Machine and to comment
  17. ;;              the code.
  18. ;;        August 1980 by CWH to run on Multics and to install
  19. ;;        new FIXFLOAT.
  20. ;;        December 1980 by JIM to fix BIGLSH not to pass LSH a second
  21. ;;        argument with magnitude greater than MACHINE-FIXNUM-PRECISION.
  22.  
  23. ;; Number of bits of precision in a fixnum and in the fields of a flonum for
  24. ;; a particular machine.  These variables should only be around at eval
  25. ;; and compile time.  These variables should probably be set up in a prelude
  26. ;; file so they can be accessible to all Macsyma files.
  27.  
  28. #.(SETQ MACHINE-FIXNUM-PRECISION
  29.     #+(OR PDP10 H6180)   36.
  30.     #+cl (integer-length most-positive-fixnum)
  31.     ;#+LISPM             24.
  32.     #+NIL             30.
  33.     #+Franz             32.
  34.  
  35. #|
  36.     MACHINE-MANTISSA-PRECISION
  37.     #+(OR PDP10 H6180)   27.
  38.     #+cl(integer-length (integer-decode-float most-positive-double-float))
  39.     ;#+LISPM             32.
  40.     #+(OR NIL Franz)     56.    ;double-float.  Long would be 113.
  41. |#
  42.     ;; Not used anymore, but keep it around anyway in case
  43.     ;; we need it later.
  44.  
  45.     MACHINE-EXPONENT-PRECISION
  46.     #+(OR PDP10 H6180)    8.
  47.     #+cl
  48.     (integer-length
  49.      (multiple-value-bind (a b)
  50.          (integer-decode-float most-positive-double-float)
  51.        b))
  52.     ;#+LISPM             11.
  53.     #+(OR NIL Franz)      8.    ;Double float.  Long would be 15.
  54.     )
  55.  
  56. ;; External variables
  57.  
  58. (DEFMVAR $FLOAT2BF NIL
  59.   "If TRUE, no MAXIMA-ERROR message is printed when a floating point number is
  60. converted to a bigfloat number.")
  61.  
  62. (DEFMVAR $BFTORAT NIL
  63.   "Controls the conversion of bigfloat numbers to rational numbers.  If
  64. FALSE, RATEPSILON will be used to control the conversion (this results in
  65. relatively small rational numbers).  If TRUE, the rational number generated
  66. will accurately represent the bigfloat.")
  67.  
  68. (DEFMVAR $BFTRUNC T
  69.   "Needs to be documented")
  70.  
  71. (DEFMVAR $FPPRINTPREC 0
  72.   "Needs to be documented"
  73.   FIXNUM)
  74.  
  75. (DEFMVAR $FPPREC 16.
  76.   "Number of decimal digits of precision to use when creating new bigfloats.
  77. One extra decimal digit in actual representation for rounding purposes.")
  78.  
  79. (DEFMVAR BIGFLOATZERO '((BIGFLOAT SIMP 56.) 0 0)
  80.   "Bigfloat representation of 0" IN-CORE)
  81. (DEFMVAR BIGFLOATONE  '((BIGFLOAT SIMP 56.) #.(EXPT 2 55.) 1)
  82.   "Bigfloat representation of 1" IN-CORE)
  83. (DEFMVAR BFHALF          '((BIGFLOAT SIMP 56.) #.(EXPT 2 55.) 0)
  84.   "Bigfloat representation of 1/2")
  85. (DEFMVAR BFMHALF      '((BIGFLOAT SIMP 56.) #.(MINUS (EXPT 2 55.)) 0)
  86.   "Bigfloat representation of -1/2")
  87. (DEFMVAR BIGFLOAT%E   '((BIGFLOAT SIMP 56.) 48968212118944587. 2)
  88.   "Bigfloat representation of %E")
  89. (DEFMVAR BIGFLOAT%PI  '((BIGFLOAT SIMP 56.) 56593902016227522. 2)
  90.   "Bigfloat representation of %PI")
  91.  
  92. ;; Internal specials
  93.  
  94. ;; Number of bits of precision in the mantissa of newly created bigfloats. 
  95. ;; FPPREC = ($FPPREC+1)*(Log base 2 of 10)
  96.  
  97. (DEFVAR FPPREC)
  98. (declare-top (FIXNUM FPPREC))
  99.  
  100. ;; FPROUND uses this to return a second value, i.e. it sets it before
  101. ;; returning.  This number represents the number of binary digits its input
  102. ;; bignum had to be shifted right to be aligned into the mantissa.  For
  103. ;; example, aligning 1 would mean shifting it FPPREC-1 places left, and
  104. ;; aligning 7 would mean shifting FPPREC-3 places left.
  105.  
  106. (DEFVAR *M)
  107. (declare-top (FIXNUM *M))
  108.  
  109. ;; *DECFP = T if the computation is being done in decimal radix.  NIL implies
  110. ;; base 2.  Decimal radix is used only during output.
  111.  
  112. (DEFVAR *DECFP NIL)
  113.  
  114. (DEFVAR MAX-BFLOAT-%PI BIGFLOAT%PI)
  115. (DEFVAR MAX-BFLOAT-%E  BIGFLOAT%E)
  116.  
  117. (declare-top (SPECIAL *CANCELLED $FLOAT $BFLOAT $RATPRINT $RATEPSILON
  118.           $DOMAIN $M1PBRANCH ADJUST)
  119.      ;; *** Local fixnum declarations ***
  120.      ;; *** Be careful of this brain-damage ***
  121.      (FIXNUM I N EXTRADIGS)
  122.      (*EXPR $BFLOAT $FLOAT)
  123.      (MUZZLED T)) 
  124.  
  125. ;; Representation of a Bigfloat:  ((BIGFLOAT SIMP precision) mantissa exponent)
  126. ;; precision -- number of bits of precision in the mantissa.  
  127. ;;         precision = (haulong mantissa)
  128. ;; mantissa -- a signed integer representing a fractional portion computed by
  129. ;;            fraction = (// mantissa (^ 2 precision)).
  130. ;; exponent -- a signed integer representing the scale of the number.
  131. ;;            The actual number represented is (f* fraction (^ 2 exponent)).
  132.  
  133. (DEFUN HIPART (X NN) (COND ((BIGP NN) (ABS X)) (T (HAIPART X NN))))
  134.  
  135. (DEFUN FPPREC1 (ASSIGN-VAR Q) 
  136.       ASSIGN-VAR ; ignored
  137.       (IF (OR (NOT (FIXNUMP Q)) (< Q 1))
  138.       (MERROR "Improper value for FPPREC:~%~M" Q))
  139.       (SETQ FPPREC (f+ 2 (HAULONG (EXPT 10. Q)))
  140.         BIGFLOATONE ($BFLOAT 1) BIGFLOATZERO ($BFLOAT 0)
  141.         BFHALF (LIST (CAR BIGFLOATONE) (CADR BIGFLOATONE) 0)
  142.         BFMHALF (LIST (CAR BIGFLOATONE) (MINUS (CADR BIGFLOATONE)) 0))
  143.       Q) 
  144.  
  145. ;; FPSCAN is called by lexical scan when a
  146. ;; bigfloat is encountered.  For example, 12.01B-3
  147. ;; would be the result of (FPSCAN '(/1 /2) '(/0 /1) '(/- /3))
  148. ;; Arguments to FPSCAN are a list of characters to the left of the
  149. ;; decimal point, to the right of the decimal point, and in the exponent.
  150.  
  151. (DEFUN FPSCAN (LFT RT EXP &AUX (*read-base* 10.) (*M 1) (*CANCELLED 0))
  152.        (SETQ EXP (READLIST EXP))
  153.        ;; Log[2](10) is 3.3219 ...
  154.        ;; This should be computed at compile time.
  155.        (BIGFLOATP
  156.     (LET ((FPPREC (PLUS 4 FPPREC (HAULONG EXP)
  157.                 (FIX (ADD1 (TIMES 3.322 (LENGTH LFT))))))
  158.           $FLOAT TEMP)
  159.          (SETQ TEMP (ADD (READLIST LFT)
  160.                  (DIV (READLIST RT) (EXPT 10. (LENGTH RT)))))
  161.          ($BFLOAT (COND ((GREATERP (ABS EXP) 1000.)
  162.                  (CONS '(MTIMES) (LIST TEMP (LIST '(MEXPT) 10. EXP))))
  163.                 (T (MUL2 TEMP (POWER 10. EXP))))))))
  164.  
  165. (DEFUN DIM-BIGFLOAT (FORM RESULT) (DIMENSION-ATOM (MAKNAM (FPFORMAT FORM)) RESULT))
  166.  
  167. (DEFUN FPFORMAT (L)
  168.  (IF (NOT (MEMQ 'SIMP (CDAR L)))
  169.      (SETQ L (CONS (CONS (CAAR L) (CONS 'SIMP (CDAR L))) (CDR L))))
  170.  (COND ((EQUAL (CADR L) 0)
  171.     (IF (NOT (EQUAL (CADDR L) 0))
  172.         (MTELL "Warning - an incorrect form for 0.0B0 has been generated."))
  173.         (LIST '|0| '|.| '|0| 'B '|0|))
  174.  (T   ;; L IS ALWAYS POSITIVE FP NUMBER
  175.   (LET ((EXTRADIGS (FIX (ADD1 (QUOTIENT (HAULONG (CADDR L)) 3.32))))
  176.     (*M 1) (*CANCELLED 0)) 
  177.        (SETQ L
  178.          ((LAMBDA (*DECFP FPPREC OF L EXPON)
  179.               (SETQ EXPON (DIFFERENCE (CADR L) OF))
  180.               (SETQ L
  181.                 (COND ((MINUSP EXPON)
  182.                    (FPQUOTIENT (INTOFP (CAR L))
  183.                            (FPINTEXPT 2 (MINUS expon) OF)))
  184.                   (T (FPTIMES* (INTOFP (CAR L))
  185.                            (FPINTEXPT 2 expon OF)))))
  186.               (SETQ FPPREC (PLUS (MINUS EXTRADIGS) FPPREC))
  187.               (LIST (FPROUND (CAR L))
  188.                 (PLUS (MINUS EXTRADIGS) *M (CADR L))))
  189.            T
  190.            (PLUS EXTRADIGS (DECIMALSIN (DIFFERENCE (CADDAR L) 2)))
  191.            (CADDAR L)
  192.            (CDR L)
  193.            NIL)))
  194.   (LET (#-CL (*print-base* 10.)   #+CL (*PRINT-BASE* 10.)
  195.     #-CL (*NOPOINT T) #+CL *PRINT-RADIX*
  196.     (L1 NIL))
  197.     (SETQ L1 (COND ((NOT $BFTRUNC) (EXPLODEC (CAR L)))
  198.            (T (DO ((L (NREVERSE (EXPLODEC (CAR L))) (CDR L)))
  199.               ((NOT (EQ '|0| (CAR L))) (NREVERSE L))))))
  200.     (NCONC (NCONS (CAR L1)) (NCONS '|.|)
  201.        (OR (AND (CDR L1)
  202.             (COND ((OR (ZEROP $FPPRINTPREC)
  203.                    (NOT (< $FPPRINTPREC $FPPREC))
  204.                    (NULL (CDDR L1)))
  205.                (CDR L1))
  206.               (T (SETQ L1 (CDR L1))
  207.                  (DO ((I $FPPRINTPREC (f1- I)) (L2))
  208.                  ((OR (< I 2) (NULL (CDR L1)))
  209.                   (COND ((NOT $BFTRUNC) (NREVERSE L2))
  210.                     (T (DO ((L3 L2 (CDR L3)))
  211.                            ((NOT (EQ '|0| (CAR L3)))
  212.                         (NREVERSE L3))))))
  213.                    (SETQ L2 (CONS (CAR L1) L2) L1 (CDR L1))))))
  214.            (NCONS '|0|))
  215.        (NCONS 'B)
  216.        (EXPLODEC (SUB1 (CADR L))))))))
  217.  
  218.  
  219. (DEFUN BIGFLOATP (X) 
  220.  (PROG NIL
  221.        (COND ((NOT ($BFLOATP X)) (RETURN NIL))
  222.          ((= FPPREC (CADDAR X)) (RETURN X))
  223.          ((> FPPREC (CADDAR X))
  224.           (SETQ X (BCONS (LIST (FPSHIFT (CADR X) (DIFFERENCE FPPREC (CADDAR X)))
  225.                    (CADDR X)))))
  226.          (T (SETQ X (BCONS (LIST (FPROUND (CADR X))
  227.                      (PLUS (CADDR X) *M FPPREC (MINUS (CADDAR X))))))))
  228.        (RETURN (COND ((EQUAL (CADR X) 0) (BCONS (LIST 0 0))) (T X))))) 
  229.  
  230. (DEFUN BIGFLOAT2RAT (X)
  231.  (SETQ X (BIGFLOATP X))
  232.  ((LAMBDA ($FLOAT2BF EXP Y SIGN) 
  233.    (SETQ EXP (COND ((MINUSP (CADR X))
  234.             (SETQ SIGN T Y (FPRATION1 (CONS (CAR X) (FPABS (CDR X)))))
  235.             (RPLACA Y (TIMES -1 (CAR Y))))
  236.            (T (FPRATION1 X))))
  237.    (COND ($RATPRINT (PRINC "RAT replaced ")
  238.             (COND (SIGN (PRINC "-")))
  239.             (PRINC (MAKNAM (FPFORMAT (CONS (CAR X) (FPABS (CDR X))))))
  240.             (PRINC " by ") (PRINC (CAR EXP)) (TYO #. forward-slash-char) (PRINC (CDR EXP))
  241.             (PRINC " = ") (SETQ X ($BFLOAT (LIST '(RAT SIMP) (CAR EXP) (CDR EXP))))
  242.             (COND (SIGN (PRINC "-")))
  243.             (PRINC (MAKNAM (FPFORMAT (CONS (CAR X) (FPABS (CDR X))))))
  244.             (TERPRI)))
  245.    EXP) 
  246.   T NIL NIL NIL))
  247.  
  248. (DEFUN FPRATION1 (X)
  249.  ((LAMBDA (FPRATEPS)
  250.        (OR (AND (EQUAL X BIGFLOATZERO) (CONS 0 1))
  251.        (PROG (Y A)
  252.          (RETURN (DO ((XX X (SETQ Y (INVERTBIGFLOAT
  253.                          (BCONS (FPDIFFERENCE (CDR XX) (CDR ($BFLOAT A)))))))
  254.                   (NUM (SETQ A (FPENTIER X))
  255.                    (PLUS (TIMES (SETQ A (FPENTIER Y)) NUM) ONUM))
  256.                   (DEN 1 (PLUS (TIMES A DEN) ODEN))
  257.                   (ONUM 1 NUM)
  258.                   (ODEN 0 DEN))
  259.                  ((AND (NOT (ZEROP DEN))
  260.                    (NOT (FPGREATERP
  261.                      (FPABS (FPQUOTIENT
  262.                             (FPDIFFERENCE (CDR X)
  263.                                   (FPQUOTIENT (CDR ($BFLOAT NUM))
  264.                                           (CDR ($BFLOAT DEN))))
  265.                             (CDR X)))
  266.                         FPRATEPS)))
  267.                   (CONS NUM DEN)))))))
  268.   (CDR ($BFLOAT (COND ($BFTORAT (LIST '(RAT SIMP) 1 (EXPTRL 2 (f1- FPPREC))))
  269.               (T $RATEPSILON))))))
  270.  
  271. ;; Convert a floating point number into a bigfloat.
  272. (DEFUN FLOATTOFP (X) 
  273.        (UNLESS $FLOAT2BF
  274.            (MTELL "Warning:  Float to bigfloat conversion of ~S~%" X))
  275.        (SETQ X (FIXFLOAT X))
  276.        (FPQUOTIENT (INTOFP (CAR X)) (INTOFP (CDR X))))
  277.  
  278. ;; Convert a bigfloat into a floating point number.
  279. (DEFMFUN FP2FLO (L)
  280.   (LET ((PRECISION (CADDAR L))
  281.     (MANTISSA (CADR L))
  282.     (EXPONENT (CADDR L))
  283.     (FPPREC #.MACHINE-MANTISSA-PRECISION)
  284.     (*M 0))
  285.     ;;Round the mantissa to the number of bits of precision of the machine,
  286.     ;;and then convert it to a floating point fraction.
  287.     (SETQ MANTISSA (QUOTIENT (FPROUND MANTISSA)
  288.                  #.(EXPT 2.0 MACHINE-MANTISSA-PRECISION)))
  289.     ;;Multiply the mantissa by the exponent portion.  I'm not sure
  290.     ;;why the exponent computation is so complicated.
  291.     (SETQ PRECISION
  292.       (ERRSET (TIMES MANTISSA (EXPT 2.0 (f+ EXPONENT (MINUS PRECISION) *M
  293.                            #.MACHINE-MANTISSA-PRECISION)))
  294.           NIL))
  295.     (IF PRECISION
  296.     (CAR PRECISION)
  297.     (MERROR "Floating point overflow in converting ~:M to flonum" L))))
  298.  
  299. ;; New machine-independent version of FIXFLOAT.  This may be buggy. - CWH
  300. ;; It is buggy!  On the PDP10 it dies on (RATIONALIZE -1.16066076E-7) 
  301. ;; which calls FLOAT on some rather big numbers.  ($RATEPSILON is approx. 
  302. ;; 7.45E-9) - JPG
  303.  
  304. #-PDP10 (DEFUN FIXFLOAT (X)
  305.   (LET (($RATEPSILON #.(EXPT 2.0 (f- MACHINE-MANTISSA-PRECISION))))
  306.        (MAXIMA-RATIONALIZE X)))
  307.  
  308. ;; Takes a flonum arg and returns a rational number corresponding to the flonum
  309. ;; in the form of a dotted pair of two integers.  Since the denominator will
  310. ;; always be a positive power of 2, this number will not always be in lowest
  311. ;; terms.
  312.  
  313. ;; PDP-10 Floating Point number format:
  314. ;; 1 bit sign -- 0 = negative 1 = positive
  315. ;; 8 bit exponent -- If positive, excess 128 encoding used, i.e.
  316. ;;   -128 exponent = 0 and +127 exponent = 255.  If number is negative,
  317. ;;   ones complement of excess 128 is used.  This is done so that the
  318. ;;   representation of the negation of a floating point number is the twos
  319. ;;   complement of the number interpreted as an integer.  If x is the number in
  320. ;;   the 8 bit field, x-128 will yield the exponent if the sign bit is off, and
  321. ;;   127-x will yield the exponent if the sign bit is on.
  322. ;; 27 bit fraction -- If the number is normalized, this fraction will be
  323. ;;   between 1/2 and 1-2^-27 inclusive, i.e. the msb of the fraction will
  324. ;;   always be 1.  The fraction is stored in two's complement so the most
  325. ;;   negative flonum is (fsc (rot 3 -1) 0) and the most positive flonum is
  326. ;;   (fsc (lsh -1 -1) 0).
  327.  
  328. ;; Old definition which explicitly hacks floating point representations.
  329. #+PDP10 (PROGN 'COMPILE
  330.   (DECLARE (CLOSED T))
  331.   (DEFUN FIXFLOAT (X) 
  332.      (PROG (NEG NUM EXPONENT DENOM) 
  333.            (COND ((LESSP X 0.0) (SETQ NEG -1.) (SETQ X (MINUS X))))
  334.            (SETQ X (LSH X 0))
  335.            (SETQ EXPONENT (DIFFERENCE (LSH X -27.) 129.))
  336.            (SETQ NUM (LSH (LSH X 9.) -9.))
  337.            (SETQ DENOM #. (f* 1 (^ 2 26.)))        ;(^ 2 26)
  338.            (COND ((LESSP EXPONENT 0)
  339.               (SETQ DENOM (TIMES DENOM (EXPT 2 (MINUS EXPONENT)))))
  340.              (T (SETQ NUM (TIMES NUM (EXPT 2 EXPONENT)))))
  341.            (IF NEG (SETQ NUM (MINUS NUM)))
  342.            (RETURN (CONS NUM DENOM)))) 
  343.   (DECLARE (CLOSED NIL))
  344.  )
  345.  
  346. ;; Format of a floating point number on the Lisp Machine:
  347. ;; 
  348. ;; High 8 bits of mantissa (plus sign bit) ---------\
  349. ;; Exponent (excess 1024) --------------\           |
  350. ;; Type of extended number --\          |           |
  351. ;; DTP-HEADER (7) ---\       |          |           |
  352. ;; Not used --\      |       |          |           |
  353. ;;            |      |       |          |           |
  354. ;;         ------------------------------------------------
  355. ;;         |  3  |   5   |   5   |     11     |     8     |
  356. ;;         ------------------------------------------------
  357. ;;         ------------------------------------------------
  358. ;;         |  3  |   5   |             24                 |
  359. ;;         ------------------------------------------------
  360. ;;            |      |                  |
  361. ;; Not used --/      |                  |
  362. ;; DTP-FIX (5) ------/                  |
  363. ;; Low 24 bits of mantissa -------------/
  364.  
  365. ;; #+LISPM
  366. ;; (DEFUN FIXFLOAT (X)
  367. ;;   (LET ((EXPONENT (f- (%P-LDB-OFFSET #O 1013 X 0) #O 2000))
  368. ;;     (NUM (%P-LDB-OFFSET #O 0010 X 0))
  369. ;;     (DENOM 1_31.))
  370. ;;     ;;Extract the high portion of the mantissa and left justify it within
  371. ;;     ;;a fixnum.
  372. ;;     (SETQ NUM (LSH NUM 16.))
  373. ;;     ;;Then extract the high 16 bits of the low portion of the mantissa and
  374. ;;     ;;store into the fixnum.
  375. ;;     (SETQ NUM (LOGIOR NUM (%P-LDB-OFFSET #O 1020 X 1)))
  376. ;;     ;;Finally, convert what we've got into a bignum by shifting left by
  377. ;;     ;;8 bits and add in low 8 bits of the low portion of the mantissa.
  378. ;;     (SETQ NUM (LOGIOR (f* NUM 1_8) (%P-LDB-OFFSET #O 0010 X 1)))
  379. ;;     (COND ((< EXPONENT 0)
  380. ;;        (SETQ DENOM (f* DENOM (^ 2 (f- EXPONENT)))))
  381. ;;       (T (SETQ NUM (f* NUM (^ 2 EXPONENT)))))
  382. ;;     (CONS NUM DENOM)))
  383.  
  384. ;; The format of a floating point number on the H6180 is very similar
  385. ;; to that on the PDP-10.  There are 8 bits of exponent, 1 bit of sign,
  386. ;; and 27 bits of mantissa in that order.  The exponent is stored
  387. ;; in twos complement, and the low order 28 bits are the mantissa
  388. ;; in twos complement.
  389.  
  390. ;; #+H6180
  391. ;; (DEFUN FIXFLOAT (X) 
  392. ;;       (PROG (NEG NUM EXPONENT DENOM) 
  393. ;;          (WHEN (LESSP X 0.0) (SETQ NEG -1) (SETQ X (-$ X)))
  394. ;;          (SETQ X (LSH X 0))
  395. ;;          (SETQ EXPONENT (LSH X -28.))
  396. ;;          (AND (> EXPONENT 177) (SETQ EXPONENT (DIFFERENCE EXPONENT 400)))
  397. ;;          (SETQ NUM (BOOLE  BOOLE-AND X 1777777777))    ;2^29-1
  398. ;;          (SETQ DENOM
  399. ;;            (TIMES 1000000000    ;2^27
  400. ;;               (COND ((LESSP EXPONENT 0)
  401. ;;                  (EXPT 2 (MINUS EXPONENT)))
  402. ;;                 (T (SETQ NUM
  403. ;;                      (TIMES NUM (EXPT 2 EXPONENT)))
  404. ;;                    1))))
  405. ;;          (IF NEG (SETQ NUM (MINUS NUM)))
  406. ;;          (RETURN (CONS NUM DENOM)))) 
  407.  
  408. (DEFUN BCONS (S) `((BIGFLOAT SIMP ,FPPREC) . ,S)) 
  409.  
  410. (DEFMFUN $BFLOAT (X) 
  411.   (LET (Y)
  412.     (COND ((BIGFLOATP X))
  413.       ((OR (NUMBERP X) (MEMQ X '($%E $%PI)))
  414.        (BCONS (INTOFP X)))
  415.       ((OR (ATOM X) (MEMQ 'array (CDAR X)))
  416.        (IF (EQ X '$%PHI)
  417.            ($BFLOAT '((MTIMES SIMP)
  418.               ((RAT SIMP) 1 2)
  419.               ((MPLUS SIMP) 1 ((MEXPT SIMP) 5 ((RAT SIMP) 1 2)))))
  420.            X))
  421.       ((EQ (CAAR X) 'MEXPT)
  422.        (IF (EQUAL (CADR X) '$%E)
  423.            (*FPEXP (CADDR X))
  424.            (EXPTBIGFLOAT ($BFLOAT (CADR X)) (CADDR X))))
  425.       ((EQ (CAAR X) 'MNCEXPT)
  426.        (LIST '(MNCEXPT) ($BFLOAT (CADR X)) (CADDR X)))
  427.       ((SETQ Y (safe-GET (CAAR X) 'FLOATPROG))
  428.        (FUNCALL Y (MAPCAR #'$BFLOAT (CDR X))))
  429.       ((OR (TRIGP (CAAR X)) (ARCP (CAAR X)) (EQ (CAAR X) '$ENTIER))
  430.        (SETQ Y ($BFLOAT (CADR X)))
  431.        (IF ($BFLOATP Y)
  432.            (COND ((EQ (CAAR X) '$ENTIER) ($ENTIER Y))
  433.              ((ARCP (CAAR X))
  434.               (SETQ Y ($BFLOAT (LOGARC (CAAR X) Y)))
  435.               (IF (FREE Y '$%I)
  436.               Y (LET ($RATPRINT) (FPARCSIMP ($RECTFORM Y)))))
  437.              ((MEMQ (CAAR X) '(%COT %SEC %CSC))
  438.               (INVERTBIGFLOAT
  439.                ($BFLOAT (LIST (NCONS (safe-GET (CAAR X) 'RECIP)) Y))))
  440.              (T ($BFLOAT (EXPONENTIALIZE (CAAR X) Y))))
  441.            (SUBST0 (LIST (NCONS (CAAR X)) Y) X)))
  442.       (T (RECUR-APPLY #'$BFLOAT X))))) 
  443.  
  444. (DEFPROP MPLUS ADDBIGFLOAT FLOATPROG)
  445. (DEFPROP MTIMES TIMESBIGFLOAT FLOATPROG)
  446. (DEFPROP %SIN SINBIGFLOAT FLOATPROG)
  447. (DEFPROP %COS COSBIGFLOAT FLOATPROG)
  448. (DEFPROP RAT RATBIGFLOAT FLOATPROG)
  449. (DEFPROP %ATAN ATANBIGFLOAT FLOATPROG)
  450. (DEFPROP %TAN TANBIGFLOAT FLOATPROG)
  451. (DEFPROP %LOG LOGBIGFLOAT FLOATPROG)
  452. (DEFPROP MABS MABSBIGFLOAT FLOATPROG)
  453.  
  454. (DEFMFUN ADDBIGFLOAT (H)
  455.        (PROG (FANS TST R NFANS)
  456.          (SETQ FANS (SETQ TST BIGFLOATZERO) NFANS 0)
  457.          (DO ((L H (CDR L))) ((NULL L))
  458.          (COND ((SETQ R (BIGFLOATP (CAR L)))
  459.             (SETQ FANS (BCONS (FPPLUS (CDR R) (CDR FANS)))))
  460.                (T (SETQ NFANS (LIST '(MPLUS) (CAR L) NFANS)))))
  461.          (RETURN (COND ((EQUAL NFANS 0) FANS)
  462.                ((EQUAL FANS TST) NFANS)
  463.                (T (SIMPLIFY (LIST '(MPLUS) FANS NFANS))))))) 
  464.  
  465. (DEFMFUN RATBIGFLOAT (L) (BCONS (FPQUOTIENT (CDAR L) (CDADR L)))) 
  466.  
  467. (DEFUN DECIMALSIN (X) 
  468.  (DO ((I (QUOTIENT (TIMES 59. X) 196.)    ;log[10](2)=.301029
  469.      (f1+ I))) (NIL) (IF (> (HAULONG (EXPT 10. I)) X) (RETURN (f1- I))))) 
  470.  
  471. (DEFMFUN ATANBIGFLOAT (X) (*FPATAN (CAR X) (CDR X))) 
  472.  
  473. (DEFMFUN *FPATAN (A Y) 
  474.  (FPEND (LET ((FPPREC (PLUS 8. FPPREC)))
  475.          (IF (NULL Y)
  476.          (IF ($BFLOATP A) (FPATAN (CDR ($BFLOAT A)))
  477.                   (LIST '(%ATAN) A))
  478.          (FPATAN2 (CDR ($BFLOAT A))
  479.               (CDR ($BFLOAT (CAR Y))))))))
  480.  
  481. (DEFUN FPATAN (X)
  482.        (PROG (TERM X2 ANS OANS ONE TWO TMP)
  483.          (SETQ ONE (INTOFP 1) TWO (INTOFP 2))
  484.          (COND ((FPGREATERP (FPABS X) ONE)
  485.             (SETQ TMP (FPQUOTIENT (FPPI) TWO))
  486.             (SETQ ANS (FPDIFFERENCE TMP (FPATAN (FPQUOTIENT ONE X))))
  487.             (RETURN (COND ((FPGREATERP ANS TMP) (FPDIFFERENCE ANS (FPPI)))
  488.                   (T ANS))))
  489.            ((FPGREATERP (FPABS X) (FPQUOTIENT ONE TWO))
  490.             (SETQ TMP (FPQUOTIENT X (FPPLUS (FPTIMES* X X) ONE)))
  491.             (SETQ X2 (FPTIMES* X TMP) TERM (SETQ ANS ONE))
  492.             (DO ((N 0 (f1+ N))) ((EQUAL ANS OANS))
  493.             (SETQ TERM
  494.                   (FPTIMES* TERM (FPTIMES* X2 (FPQUOTIENT
  495.                                (INTOFP (f+ 2 (f* 2 N)))
  496.                                    (INTOFP (f+ (f* 2 N) 3))))))
  497.             (SETQ OANS ANS ANS (FPPLUS TERM ANS)))
  498.             (SETQ ANS (FPTIMES* TMP ANS)))
  499.            (T (SETQ ANS X X2 (FPMINUS (FPTIMES* X X)) TERM X)
  500.               (DO ((N 3 (f+ N 2))) ((EQUAL ANS OANS))
  501.               (SETQ TERM (FPTIMES* TERM X2))
  502.               (SETQ OANS ANS 
  503.                     ANS (FPPLUS ANS (FPQUOTIENT TERM (INTOFP N)))))))
  504.          (RETURN ANS)))
  505.  
  506. (DEFUN FPATAN2 (Y X) ; ATAN(Y/X) from -PI to PI
  507.        (COND ((EQUAL (CAR X) 0)       ; ATAN(INF), but what sign?
  508.           (COND ((EQUAL (CAR Y) 0) (MERROR "ATAN(0//0) has been generated."))
  509.             ((MINUSP (CAR Y))
  510.              (FPQUOTIENT (FPPI) (INTOFP -2)))
  511.             (T (FPQUOTIENT (FPPI) (INTOFP 2)))))
  512.          ((SIGNP G (CAR X))
  513.           (COND ((SIGNP G (CAR Y)) (FPATAN (FPQUOTIENT Y X)))
  514.             (T (FPMINUS (FPATAN (FPQUOTIENT Y X))))))
  515.          ((SIGNP G (CAR Y))
  516.           (FPPLUS (FPPI) (FPATAN (FPQUOTIENT Y  X))))
  517.          (T (FPDIFFERENCE (FPATAN (FPQUOTIENT Y X)) (FPPI))))) 
  518.  
  519. (DEFUN TANBIGFLOAT (A)
  520.  (SETQ A (CAR A)) 
  521.  (FPEND (LET ((FPPREC (PLUS 8. FPPREC)))
  522.          (COND (($BFLOATP A)
  523.             (SETQ A (CDR ($BFLOAT A)))
  524.             (FPQUOTIENT (FPSIN A T) (FPSIN A NIL)))
  525.            (T (LIST '(%TAN) A))))))     
  526.  
  527. ;; Returns a list of a mantissa and an exponent.
  528. (DEFUN INTOFP (L) 
  529.        (COND ((NOT (ATOM L)) ($BFLOAT L))
  530.          ((FLOATP L) (FLOATTOFP L))
  531.          ((EQUAL 0 L) '(0 0))
  532.          ((EQ L '$%PI) (FPPI))
  533.          ((EQ L '$%E) (FPE))
  534.          (T (LIST (FPROUND L) (PLUS *M FPPREC))))) 
  535.  
  536. ;; It seems to me that this function gets called on an integer
  537. ;; and returns the mantissa portion of the mantissa/exponent pair.
  538.  
  539. ;; "STICKY BIT" CALCULATION FIXED 10/14/75 --RJF
  540. ;; BASE must not get temporarily bound to NIL by being placed
  541. ;; in a PROG list as this will confuse stepping programs.
  542.  
  543. (DEFUN FPROUND (L &AUX #-cl (*print-base* 10.) #+cl (*print-base* 10.)
  544.                #-cl (*NOPOINT T)#+cl *print-radix*
  545.           )
  546.   (PROG () 
  547.     (COND
  548.      ((NULL *DECFP)
  549.       ;;*M will be positive if the precision of the argument is greater than
  550.       ;;the current precision being used.
  551.       (SETQ *M (f- (HAULONG L) FPPREC))
  552.       (COND ((= *M 0) (SETQ *CANCELLED 0) (RETURN L)))
  553.       ;;FPSHIFT is essentially LSH.
  554.       (SETQ ADJUST (FPSHIFT 1 (SUB1 *M)))
  555.       (COND ((MINUSP L) (SETQ ADJUST (MINUS ADJUST))))
  556.       (SETQ L (PLUS L ADJUST))
  557.       (SETQ *M (f- (HAULONG L) FPPREC))
  558.       (SETQ *CANCELLED (ABS *M))
  559.          
  560.       (COND (#+cl
  561.          (zerop (HIPART L (MINUS *M)))
  562.          #-cl
  563.          (SIGNP E (HIPART L (MINUS *M)))
  564.          ;ONLY ZEROES SHIFTED OFF
  565.          (RETURN (FPSHIFT (FPSHIFT L (DIFFERENCE -1 *M))
  566.                   1)))            ; ROUND TO MAKE EVEN
  567.         (T (RETURN (FPSHIFT L (MINUS *M))))))
  568.      (T
  569.       (SETQ *M (DIFFERENCE (FLATSIZE (ABS L)) FPPREC))
  570.       (SETQ ADJUST (FPSHIFT 1 (SUB1 *M)))
  571.       (COND ((MINUSP L) (SETQ ADJUST (MINUS ADJUST))))
  572.       (SETQ ADJUST (TIMES 5 ADJUST))
  573.       (SETQ *M
  574.         (DIFFERENCE (FLATSIZE (ABS (SETQ L (PLUS L ADJUST))))
  575.                 FPPREC))
  576.       (RETURN (FPSHIFT L (MINUS *M)))))))
  577.  
  578. ;; Compute (* L (expt d n)) where D is 2 or 10 depending on
  579. ;; *decfp. Throw away an fractional part by truncating to zero.
  580. (DEFUN FPSHIFT (L N) 
  581.   (COND ((NULL *DECFP)
  582.      (cond ((and (minusp n) (minusp l))
  583.         ;; Left shift of negative number requires some
  584.         ;; care. (That is, (truncate l (expt 2 n)), but use
  585.         ;; shifts instead.)
  586.         (- (ash (- l) n)))
  587.            (t
  588.         (ash l n))))
  589.     ((GREATERP N 0.)
  590.      (TIMES L (EXPT 10. N)))
  591.     ((LESSP N 0.)
  592.      (QUOTIENT L (EXPT 10. (MINUS N))))
  593.     (T L)))
  594.  
  595. ;; Bignum LSH -- N is assumed (and declared above) to be a fixnum.
  596. ;; This isn't really LSH, since the sign bit isn't propagated when
  597. ;; shifting to the right, i.e. (BIGLSH -100 -3) = -40, whereas
  598. ;; (LSH -100 -3) = 777777777770 (on a 36 bit machine).
  599. ;; This actually computes (TIMES X (EXPT 2 N)).  As of 12/21/80, this function
  600. ;; was only called by FPSHIFT.  I would like to hear an argument as why this
  601. ;; is more efficient than simply writing (TIMES X (EXPT 2 N)).  Is the
  602. ;; intermediate result created by (EXPT 2 N) the problem?  I assume that
  603. ;; EXPT tries to LSH when possible.
  604.  
  605. (DEFUN BIGLSH (X N)
  606.   (COND
  607.    ;; In MacLisp, the result is undefined if the magnitude of the
  608.    ;; second argument is greater than 36.
  609.    ((AND (NOT (BIGP X))
  610.      (< N #.(f- MACHINE-FIXNUM-PRECISION))) 0)
  611.    ;; Either we are shifting a fixnum to the right, or shifting
  612.    ;; a fixnum to the left, but not far enough left for it to become
  613.    ;; a bignum.
  614.    ((AND (NOT (BIGP X)) 
  615.      (OR (<= N 0)
  616.          (< (PLUS (HAULONG X) N) #.MACHINE-FIXNUM-PRECISION)))
  617.     ;; The form which follows is nearly identical to (ASH X N), however
  618.     ;; (ASH -100 -20) = -1, whereas (BIGLSH -100 -20) = 0.
  619.     (IF (>= X 0)
  620.     (LSH X N)
  621.     (f- (BIGLSH (MINUS X) N)))) ;(minus x) may be a bignum even is x is a fixnum.
  622.    ;; If we get here, then either X is a bignum or our answer is
  623.    ;; going to be a bignum.
  624.    ((< N 0)
  625.     (COND ((> (ABS N) (HAULONG X)) 0)
  626.       ((GREATERP X 0)
  627.        (HIPART X (PLUS (HAULONG X) N)))
  628.       (T (MINUS (HIPART X (PLUS (HAULONG X) N))))))
  629.    ((= N 0) X)
  630.    ;; Isn't this the kind of optimization that compilers are
  631.    ;; supposed to make?
  632.    ((< N #.(f1- MACHINE-FIXNUM-PRECISION)) (TIMES X (LSH 1 N)))
  633.    (T (TIMES X (EXPT 2 N)))))
  634.  
  635.  
  636. (DEFUN FPEXP (X)       
  637.   (PROG (R S)
  638.     (IF (NOT (SIGNP GE (CAR X)))
  639.         (RETURN (FPQUOTIENT (FPONE) (FPEXP (FPABS X)))))
  640.     (SETQ R (FPINTPART X))
  641.     (RETURN (COND ((LESSP R 2) (FPEXP1 X))
  642.               (T (SETQ S (FPEXP1 (FPDIFFERENCE X (INTOFP R))))
  643.              (FPTIMES* S
  644.                    (CDR (BIGFLOATP
  645.                      ((LAMBDA (FPPREC R) (BCONS (FPEXPT (FPE) R)))    ; patch for full precision %E
  646.                       (PLUS FPPREC (HAULONG R) -1)
  647.                       R)))))))))
  648.  
  649. (DEFUN FPEXP1 (X) 
  650.        (PROG (TERM ANS OANS) 
  651.          (SETQ ANS (SETQ TERM (FPONE)))
  652.          (DO ((N
  653.          1.
  654.          (ADD1 N)))
  655.          ((EQUAL ANS OANS))
  656.          (SETQ TERM (FPQUOTIENT (FPTIMES* X TERM) (INTOFP N)))
  657.          (SETQ OANS ANS)
  658.          (SETQ ANS (FPPLUS ANS TERM)))
  659.          (RETURN ANS))) 
  660.  
  661. ;; Does one higher precision to round correctly.
  662. ;; A and B are each a list of a mantissa and an exponent.
  663. (DEFUN FPQUOTIENT (A B) 
  664.        (COND ((EQUAL (CAR B) 0)
  665.           (MERROR "PQUOTIENT by zero"))
  666.          ((EQUAL (CAR A) 0) '(0 0))
  667.          (T (LIST (FPROUND (QUOTIENT (FPSHIFT (CAR A)
  668.                           (PLUS 3 FPPREC))
  669.                      (CAR B)))
  670.               (PLUS -3 (DIFFERENCE (CADR A) (CADR B)) *M))))) 
  671.  
  672. (DEFUN FPGREATERP (A B) (FPPOSP (FPDIFFERENCE A B))) 
  673.  
  674. (DEFUN FPLESSP (A B) (FPPOSP (FPDIFFERENCE B A))) 
  675.  
  676. (DEFUN FPPOSP (X) (GREATERP (CAR X) 0)) 
  677.  
  678. (DEFmfUN FPMIN NA
  679.  (PROG (MIN) 
  680.       (SETQ MIN (ARG 1))
  681.       (DO ((I 2 (f1+ I))) ((> I NA))
  682.           (IF (FPLESSP (ARG I) MIN) (SETQ MIN (ARG I))))
  683.       (RETURN MIN)))
  684.  
  685. ;; (FPE) RETURN BIG FLOATING POINT %E.  IT RETURNS (CDR BIGFLOAT%E) IF RIGHT
  686. ;; PRECISION.  IT RETURNS TRUNCATED BIGFLOAT%E IF POSSIBLE, ELSE RECOMPUTES.
  687. ;; IN ANY CASE, BIGFLOAT%E IS SET TO LAST USED VALUE. 
  688.  
  689. (DEFUN FPE NIL
  690.        (COND ((= FPPREC (CADDAR BIGFLOAT%E)) (CDR BIGFLOAT%E))
  691.          ((< FPPREC (CADDAR BIGFLOAT%E))
  692.           (CDR (SETQ BIGFLOAT%E (BIGFLOATP BIGFLOAT%E))))
  693.          ((< FPPREC (CADDAR MAX-BFLOAT-%E))
  694.           (CDR (SETQ BIGFLOAT%E (BIGFLOATP MAX-BFLOAT-%E))))
  695.          (T (CDR (SETQ MAX-BFLOAT-%E (SETQ BIGFLOAT%E (*FPEXP 1)))))))
  696.  
  697. (DEFUN FPPI NIL
  698.        (COND ((= FPPREC (CADDAR BIGFLOAT%PI)) (CDR BIGFLOAT%PI))
  699.          ((< FPPREC (CADDAR BIGFLOAT%PI))
  700.           (CDR (SETQ BIGFLOAT%PI (BIGFLOATP BIGFLOAT%PI))))
  701.          ((< FPPREC (CADDAR MAX-BFLOAT-%PI))
  702.           (CDR (SETQ BIGFLOAT%PI (BIGFLOATP MAX-BFLOAT-%PI))))
  703.          (T (CDR (SETQ MAX-BFLOAT-%PI (SETQ BIGFLOAT%PI (FPPI1)))))))
  704.  
  705. (DEFUN FPONE NIL 
  706.        (COND (*DECFP (INTOFP 1)) ((= FPPREC (CADDAR BIGFLOATONE)) (CDR BIGFLOATONE))
  707.          (T (INTOFP 1)))) 
  708.  
  709. ;; COMPPI computes PI to N bits.
  710. ;; That is, (COMPPI N)/(2.0^N) is an approximation to PI.
  711.  
  712. (DEFUN COMPPI (N) 
  713.        (PROG (A B C) 
  714.          (SETQ A (EXPT 2 N))
  715.          (SETQ C (PLUS (TIMES 3 A) (SETQ B (*QUO A 8.))))
  716.           (DO ((I 4 (f+ I 2)))
  717.          ((ZEROP B))
  718.          (SETQ B (*QUO (TIMES B (f1- I) (f1- I))
  719.                    (TIMES 4 I (f1+ I))))
  720.          (SETQ C (PLUS C B)))
  721.          (RETURN C))) 
  722.  
  723. (DEFUN FPPI1 NIL 
  724.        (BCONS (LIST (FPROUND (COMPPI (PLUS FPPREC 3))) (PLUS -3 *M)))) 
  725.  
  726. (DEFmfUN FPMAX NA
  727.  (PROG (MAX) 
  728.       (SETQ MAX (ARG 1))
  729.       (DO ((I 2 (f1+ I))) ((> I NA))
  730.           (IF (FPGREATERP (ARG I) MAX) (SETQ MAX (ARG I))))
  731.       (RETURN MAX)))
  732.  
  733. (DEFUN FPDIFFERENCE (A B) (FPPLUS A (FPMINUS B))) 
  734.  
  735. (DEFUN FPMINUS (X) (IF (EQUAL (CAR X) 0) X (LIST (MINUS (CAR X)) (CADR X)))) 
  736.  
  737. (DEFUN FPPLUS (A B) 
  738.        (PROG (*M EXP MAN STICKY) 
  739.          (SETQ *CANCELLED 0)
  740.     (COND ((EQUAL (CAR A) 0) (RETURN B))
  741.           ((EQUAL (CAR B) 0) (RETURN A)))
  742.     (SETQ EXP (DIFFERENCE (CADR A) (CADR B)))
  743.     (SETQ MAN (COND ((EQUAL EXP 0)
  744.               (SETQ STICKY 0)
  745.               (FPSHIFT (PLUS (CAR A) (CAR B)) 2))
  746.              ((GREATERP EXP 0)
  747.               (SETQ STICKY (HIPART (CAR B) (DIFFERENCE 1 EXP)))
  748.               (SETQ STICKY (COND ((SIGNP E STICKY) 0)
  749.                          ((SIGNP L (CAR B)) -1)
  750.                          (T 1)))
  751.                                        ; COMPUTE STICKY BIT
  752.               (PLUS (FPSHIFT (CAR A) 2)
  753.                                        ; MAKE ROOM FOR GUARD DIGIT & STICKY BIT
  754.                 (FPSHIFT (CAR B) (DIFFERENCE 2 EXP))))
  755.              (T (SETQ STICKY (HIPART (CAR A) (ADD1 EXP)))
  756.                 (SETQ STICKY (COND ((SIGNP E STICKY) 0)
  757.                            ((SIGNP L (CAR A)) -1)
  758.                            (T 1)))
  759.                 (PLUS (FPSHIFT (CAR B) 2)
  760.                   (FPSHIFT (CAR A) (PLUS 2 EXP))))))
  761.          (SETQ MAN (PLUS MAN STICKY))
  762.          (RETURN (COND ((EQUAL MAN 0) '(0 0))
  763.                (T (SETQ MAN (FPROUND MAN))
  764.                   (SETQ EXP
  765.                     (PLUS -2 *M (MAX (CADR A) (CADR B))))
  766.                   (LIST MAN EXP)))))) 
  767.  
  768. (DEFUN FPTIMES* (A B) 
  769.        (COND ((OR (EQUAL (CAR A) 0) (EQUAL (CAR B) 0)) '(0 0))
  770.          (T (LIST (FPROUND (TIMES (CAR A) (CAR B)))
  771.               (PLUS *M (CADR A) (CADR B) (MINUS FPPREC)))))) 
  772.  
  773. ;; Don't use the symbol BASE since it is SPECIAL.
  774.  
  775. (DEFUN FPINTEXPT (INT NN FIXPREC)            ;INT is integer
  776.        (SETQ FIXPREC (QUOTIENT FIXPREC (LOG2 INT)))    ;NN is pos
  777.        (LET ((BAS (INTOFP (EXPT INT (MIN NN FIXPREC))))) 
  778.         (COND ((GREATERP NN FIXPREC)
  779.            (FPTIMES* (INTOFP (EXPT INT (REMAINDER NN FIXPREC)))
  780.                  (FPEXPT BAS (QUOTIENT NN FIXPREC))))
  781.           (T BAS))))
  782.  
  783. ;; NN is positive or negative integer
  784.  
  785. (DEFUN FPEXPT (P NN) 
  786.        (COND ((EQUAL NN 0.) (FPONE))
  787.          ((EQUAL NN 1.) P)
  788.          ((LESSP NN 0.) (FPQUOTIENT (FPONE) (FPEXPT P (MINUS NN))))
  789.          (T (PROG (U) 
  790.               (COND ((ODDP NN) (SETQ U P))
  791.                 (T (SETQ U (FPONE))))
  792.               (DO ((II (QUOTIENT NN 2.) (QUOTIENT II 2.)))
  793.               ((ZEROP II))
  794.               (SETQ P (FPTIMES* P P))
  795.               (COND ((ODDP II) (SETQ U (FPTIMES* U P)))))
  796.               (RETURN U))))) 
  797.  
  798. (declare-top (NOTYPE N))
  799.  
  800. (DEFUN EXPTBIGFLOAT (P N) 
  801.   (COND ((EQUAL N 1) P)
  802.     ((EQUAL N 0) ($BFLOAT 1))
  803.     ((NOT ($BFLOATP P)) (LIST '(MEXPT) P N))
  804.     ((EQUAL (CADR P) 0) ($BFLOAT 0))
  805.     ((AND (LESSP (CADR P) 0) (RATNUMP N))
  806.      ($BFLOAT
  807.       ($EXPAND (LIST '(MTIMES)
  808.              ($BFLOAT ((LAMBDA ($DOMAIN $M1PBRANCH) (POWER -1 N))
  809.                    '$COMPLEX T))
  810.              (EXPTBIGFLOAT (BCONS (FPMINUS (CDR P))) N)))))
  811.     ((AND (LESSP (CADR P) 0) (NOT (INTEGERP N)))
  812.      (COND ((OR (EQUAL N 0.5) (EQUAL N BFHALF))
  813.         (EXPTBIGFLOAT P '((RAT SIMP) 1 2)))
  814.            ((OR (EQUAL N -0.5) (EQUAL N BFMHALF))
  815.         (EXPTBIGFLOAT P '((RAT SIMP) -1 2)))
  816.            (($BFLOATP (SETQ N ($BFLOAT N)))
  817.         (COND ((EQUAL N ($BFLOAT (FPENTIER N)))
  818.                (EXPTBIGFLOAT P (FPENTIER N)))
  819.               (T  ;; for P<0: P^N = (-P)^N*cos(pi*N) + i*(-P)^N*sin(pi*N)
  820.              (SETQ P (EXPTBIGFLOAT (BCONS (FPMINUS (CDR P))) N)
  821.                    N ($BFLOAT `((MTIMES) $%PI ,N)))
  822.              (ADD2 ($BFLOAT `((MTIMES) ,P ,(*FPSIN N NIL)))
  823.                    `((MTIMES SIMP) ,($BFLOAT `((MTIMES) ,P ,(*FPSIN N T)))
  824.                            $%I)))))
  825.            (T (LIST '(MEXPT) P N))))
  826.     ((AND (RATNUMP N) (LESSP (CADDR N) 10.))
  827.      (BCONS (FPEXPT (FPROOT P (CADDR N)) (CADR N))))
  828.     ((NOT (INTEGERP N))
  829.      (SETQ N ($BFLOAT N))
  830.      (COND
  831.       ((NOT ($BFLOATP N)) (LIST '(MEXPT) P N))
  832.       (T
  833.        ((LAMBDA (EXTRABITS) 
  834.          (SETQ 
  835.           P
  836.           ((LAMBDA (FPPREC) 
  837.                (FPEXP (FPTIMES* (CDR (BIGFLOATP N))
  838.                     (FPLOG (CDR (BIGFLOATP P))))))
  839.            (PLUS EXTRABITS FPPREC)))
  840.          (SETQ P (LIST (FPROUND (CAR P))
  841.                (PLUS (MINUS EXTRABITS) *M (CADR P))))
  842.          (BCONS P))
  843.         (MAX 1 (PLUS (CADDR N) (HAULONG (CADDR P))))))))
  844.                    ; The number of extra bits required 
  845.     ((LESSP N 0) (INVERTBIGFLOAT (EXPTBIGFLOAT P (MINUS N))))
  846.     (T (BCONS (FPEXPT (CDR P) N)))))
  847.  
  848. (defun fproot (a n) ; computes a^(1/n)  see Fitch, SIGSAM Bull Nov 74
  849.  (let* ((ofprec fpprec) (fpprec (f+ fpprec 2))        ;assumes a>0 n>=2
  850.         (bk (fpexpt (intofp 2)
  851.             (add1 (quotient (cadr (setq a (cdr (bigfloatp a)))) n)))))
  852.       (do ((x bk
  853.           (fpdifference
  854.            x (setq bk (fpquotient (fpdifference
  855.                        x (fpquotient a (fpexpt x n1))) n))))
  856.        (n1 (sub1 n))
  857.        (n (intofp n)))
  858.       ((or (equal bk '(0 0))
  859.            (greaterp (difference (cadr x) (cadr bk)) ofprec)) (setq a x))))
  860.  (list (fpround (car a)) (plus -2 *m (cadr a))))
  861.  
  862. (DEFUN TIMESBIGFLOAT (H) 
  863.   (PROG (FANS TST R NFANS) 
  864.     (SETQ FANS (SETQ TST (BCONS (FPONE))) NFANS 1)
  865.     (DO ((L H (CDR L))) ((NULL L))
  866.         (IF (SETQ R (BIGFLOATP (CAR L)))
  867.         (SETQ FANS (BCONS (FPTIMES* (CDR R) (CDR FANS))))
  868.         (SETQ NFANS (LIST '(MTIMES) (CAR L) NFANS))))
  869.     (RETURN (COND ((EQUAL NFANS 1) FANS)
  870.               ((EQUAL FANS TST) NFANS)
  871.               (T (SIMPLIFY (LIST '(MTIMES) FANS NFANS))))))) 
  872.  
  873. (DEFUN INVERTBIGFLOAT (A) 
  874.    (IF (BIGFLOATP A) (BCONS (FPQUOTIENT (FPONE) (CDR A)))
  875.              (SIMPLIFY (LIST '(MEXPT) A -1))))
  876.  
  877. (DEFUN *FPEXP (A) 
  878.   (FPEND (LET ((FPPREC (PLUS 8. FPPREC)))
  879.        (IF ($BFLOATP (SETQ A ($BFLOAT A)))
  880.            (FPEXP (CDR A))
  881.            (LIST '(MEXPT) '$%E A)))))
  882.  
  883. (DEFUN *FPSIN (A FL) 
  884.   (FPEND (LET ((FPPREC (PLUS 8. FPPREC)))
  885.        (COND (($BFLOATP A) (FPSIN (CDR ($BFLOAT A)) FL))
  886.          (FL (LIST '(%SIN) A))
  887.          (T (LIST '(%COS) A))))))
  888.  
  889. (DEFUN FPEND (A)
  890.   (COND ((EQUAL (CAR A) 0) (BCONS A))
  891.     ((NUMBERP (CAR A))
  892.      (SETQ A (LIST (FPROUND (CAR A)) (PLUS -8. *M (CADR A))))
  893.      (BCONS A))
  894.     (T A))) 
  895.  
  896. (DEFUN FPARCSIMP (E)  ; needed for e.g. ASIN(.123567812345678B0) with 
  897.               ; FPPREC 16, to get rid of the miniscule imaginary 
  898.               ; part of the a+bi answer.
  899.   (IF (AND (MPLUSP E) (NULL (CDDDR E))
  900.        (MTIMESP (CADDR E)) (NULL (CDDDR (CADDR E)))
  901.        ($BFLOATP (CADR (CADDR E)))
  902.        (EQ (CADDR (CADDR E)) '$%I)
  903.        (< (CADDR (CADR (CADDR E))) (f+ (f- FPPREC) 2)))
  904.       (CADR E)
  905.       E))
  906. (declare-top (FIXNUM N))
  907.  
  908. (DEFUN SINBIGFLOAT (X) (*FPSIN (CAR X) T)) 
  909.  
  910. (DEFUN COSBIGFLOAT (X) (*FPSIN (CAR X) NIL)) 
  911.  
  912. ;; THIS VERSION OF FPSIN COMPUTES SIN OR COS TO PRECISION FPPREC,
  913. ;; BUT CHECKS FOR THE POSSIBILITY OF CATASTROPHIC CANCELLATION DURING
  914. ;; ARGUMENT REDUCTION (E.G. SIN(N*%PI+EPSILON)) 
  915. ;; *FPSINCHECK WILL CAUSE PRINTOUT OF ADDITIONAL INFO WHEN
  916. ;; EXTRA PRECISION IS NEEDED FOR SIN/COS CALCULATION.  KNOWN
  917. ;; BAD FEATURES:  IT IS NOT NECESSARY TO USE EXTRA PRECISION FOR, E.G.
  918. ;; SIN(PI/2), WHICH IS NOT NEAR ZERO, BUT  EXTRA
  919. ;; PRECISION IS USED SIN IT IS NEEDED FOR COS(PI/2).
  920. ;; PRECISION SEEMS TO BE 100% SATSIFACTORY FOR LARGE ARGUMENTS, E.G.
  921. ;; SIN(31415926.0B0), BUT LESS SO FOR SIN(3.1415926B0).  EXPLANATION
  922. ;; NOT KNOWN.  (9/12/75  RJF)
  923. (declare-top (SPECIAL *FPSINCHECK))
  924. (SETQ *FPSINCHECK NIL)
  925.  
  926. (DEFUN FPSIN (X FL) 
  927.   (PROG (PIBY2 R SIGN RES K *CANCELLED) 
  928.     (SETQ SIGN (COND (FL (SIGNP G (CAR X))) (T)) X (FPABS X))
  929.     (COND ((EQUAL (CAR X) 0)
  930.            (RETURN (COND (FL (INTOFP 0)) (T (INTOFP 1))))))
  931.     (RETURN 
  932.      (CDR
  933.       (BIGFLOATP
  934.        ((LAMBDA (FPPREC XT *CANCELLED OLDPREC) 
  935.           (PROG (X) 
  936.                LOOP (SETQ X (CDR (BIGFLOATP XT)))
  937.                 (SETQ PIBY2 (FPQUOTIENT (FPPI) (INTOFP 2)))
  938.             (SETQ R (FPINTPART (FPQUOTIENT X PIBY2)))
  939.             (SETQ X
  940.               (FPPLUS X
  941.                   (FPTIMES* (INTOFP (MINUS R))
  942.                         PIBY2)))
  943.             (SETQ K *CANCELLED)
  944.             (FPPLUS X (FPMINUS PIBY2))
  945.             (SETQ *CANCELLED (MAX K *CANCELLED))
  946.             (COND (*FPSINCHECK
  947.                (PRINT `(*CANC= ,*CANCELLED FPPREC= ,FPPREC
  948.                        OLDPREC= ,OLDPREC))))
  949.  
  950.             (COND
  951.              ((NOT (GREATERP OLDPREC
  952.                      (DIFFERENCE FPPREC
  953.                          *CANCELLED)))
  954.               (SETQ R (REMAINDER R 4))
  955.               (SETQ RES
  956.                 (COND (FL (COND ((= R 0) (FPSIN1 X))
  957.                         ((= R 1) (FPCOS1 X))
  958.                         ((= R 2) (FPMINUS (FPSIN1 X)))
  959.                         ((= R 3) (FPMINUS (FPCOS1 X)))))
  960.                   (T (COND ((= R 0) (FPCOS1 X))
  961.                        ((= R 1) (FPMINUS (FPSIN1 X)))
  962.                        ((= R 2) (FPMINUS (FPCOS1 X)))
  963.                        ((= R 3) (FPSIN1 X))))))
  964.               (RETURN (BCONS (COND (SIGN RES) (T (FPMINUS RES))))))
  965.              (T (SETQ FPPREC (PLUS FPPREC *CANCELLED))
  966.             (GO LOOP)))))
  967.         (MAX FPPREC (PLUS FPPREC (CADR X)))
  968.         (BCONS X)
  969.         0
  970.         FPPREC)))))) 
  971.  
  972. (DEFUN FPCOS1 (X) (FPSINCOS1 X NIL))
  973.  
  974. ;; Compute SIN or COS in (0,PI/2).  FL is T for SIN, NIL for COS.
  975. (DEFUN FPSINCOS1 (X FL)
  976.        (PROG (ANS TERM OANS X2)
  977.          (SETQ ANS (IF FL X (INTOFP 1))
  978.            X2 (FPMINUS(FPTIMES* X X)))
  979.          (SETQ TERM ANS)
  980.          (DO ((N (COND (FL 3) (T 2)) (PLUS N 2))) ((EQUAL ANS OANS))
  981.          (SETQ TERM (FPTIMES* TERM (FPQUOTIENT X2 (INTOFP (f* N (SUB1 N))))))
  982.          (SETQ OANS ANS ANS (FPPLUS ANS TERM)))
  983.          (RETURN ANS)))
  984.  
  985. (DEFUN FPSIN1(X) (FPSINCOS1 X T)) 
  986.  
  987. (DEFUN FPABS (X) 
  988.        (COND ((SIGNP GE (CAR X)) X)
  989.          (T (CONS (MINUS (CAR X)) (CDR X))))) 
  990.  
  991. (DEFMFUN FPENTIER (F) (LET ((FPPREC (CADDAR F))) (FPINTPART (CDR F))))
  992.  
  993. (DEFUN FPINTPART (F) 
  994.        (PROG (M) 
  995.          (SETQ M (DIFFERENCE FPPREC (CADR F)))
  996.          (RETURN (COND ((GREATERP M 0)
  997.                 (QUOTIENT (CAR F) (EXPT 2 M)))
  998.                (T (TIMES (CAR F) (EXPT 2 (MINUS M)))))))) 
  999.  
  1000. (DEFUN LOGBIGFLOAT (A) 
  1001.  ((LAMBDA (MINUS)
  1002.    (SETQ A ((LAMBDA (FPPREC) 
  1003.          (COND (($BFLOATP (CAR A))
  1004.             (SETQ A ($BFLOAT (CAR A)))
  1005.             (COND ((ZEROP (CADR A)) (MERROR "LOG(0.0B0) has been generated"))
  1006.               ((MINUSP (CADR A))
  1007.                (SETQ MINUS T) (FPLOG (LIST (MINUS (CADR A)) (CADDR A))))
  1008.               (T (FPLOG (CDR A)))))
  1009.            (T (LIST '(%LOG) (CAR A)))))
  1010.         (PLUS 2 FPPREC)))
  1011.    (COND ((NUMBERP (CAR A))
  1012.       (SETQ A
  1013.         (if (zerop (car a))
  1014.             (list 0 0)
  1015.             (LIST (FPROUND (CAR A)) (PLUS -2 *M (CADR A)))))
  1016.       (SETQ A (BCONS A))))
  1017.    (COND (MINUS (ADD A (MUL '$%I ($BFLOAT '$%PI)))) (T A)))
  1018.   NIL)) 
  1019.  
  1020. ;;; Computes the log of a bigfloat number.
  1021. ;;;
  1022. ;;; Uses the series
  1023. ;;;
  1024. ;;; log(1+x) = sum((x/(x+2))^(2*n+1)/(2*n+1),n,0,inf);
  1025. ;;;
  1026. ;;;
  1027. ;;;                  INF      x   2 n + 1
  1028. ;;;                  ====  (-----)
  1029. ;;;                  \      x + 2
  1030. ;;;          =  2     >    --------------
  1031. ;;;                  /        2 n + 1
  1032. ;;;                  ====
  1033. ;;;                  n = 0
  1034. ;;;
  1035. ;;;
  1036. ;;; which converges for x > 0.
  1037. ;;;
  1038. ;;; Note that FPLOG is given 1+X, not X.
  1039. ;;;
  1040. ;;; However, to aid convergence of the series, we scale 1+x until 1/e
  1041. ;;; < 1+x <= e.
  1042. ;;;
  1043. (DEFUN FPLOG (X) 
  1044.        (PROG (OVER TWO ANS OLDANS TERM E sum) 
  1045.          (IF (NOT (GREATERP (CAR X) 0))
  1046.          (MERROR "Non-positive argument to FPLOG"))
  1047.          (SETQ E (FPE)
  1048.            OVER (FPQUOTIENT (FPONE) E)
  1049.            ANS 0)
  1050.          ;; Scale X until 1/e < X <= E.  ANS keeps track of how
  1051.          ;; many factors of E were used.  Set X to NIL if X is E.
  1052.          (DO () (NIL)
  1053.          (COND ((EQUAL X E) (SETQ X NIL) (RETURN NIL))
  1054.                ((AND (FPLESSP X E) (FPLESSP OVER X))
  1055.             (RETURN NIL))
  1056.                ((FPLESSP X OVER)
  1057.             (SETQ X (FPTIMES* X E))
  1058.             (SETQ ANS (SUB1 ANS)))
  1059.                (T (SETQ ANS (ADD1 ANS))
  1060.               (SETQ X (FPQUOTIENT X e)))))
  1061.          (COND ((NULL X) (RETURN (INTOFP (ADD1 ANS)))))
  1062.          ;; Prepare X for the series.  The series is for 1 + x, so
  1063.          ;; get x from our X.  TERM is (x/(x+2)).  X becomes
  1064.          ;; (x/(x+2))^2.
  1065.          (SETQ X (FPDIFFERENCE  X (FPONE))
  1066.            ANS (INTOFP ANS))
  1067.          (SETQ 
  1068.           X
  1069.           (FPEXPT (SETQ TERM
  1070.                 (FPQUOTIENT X (FPPLUS X (SETQ TWO (INTOFP 2)))))
  1071.               2))
  1072.          ;; Sum the series until the sum (in ANS) doesn't change
  1073.          ;; anymore.
  1074.          (setq sum (intofp 0))
  1075.          (DO ((N 1 (+ N 2)))
  1076.          ((EQUAL sum OLDANS))
  1077.           (SETQ OLDANS sum)
  1078.           (SETQ sum
  1079.             (FPPLUS
  1080.              sum
  1081.              (FPQUOTIENT TERM (INTOFP N))))
  1082.           (SETQ TERM (FPTIMES* TERM X)))
  1083.          
  1084.          (return (fpplus ANS (fptimes* two sum)))))
  1085.  
  1086. (DEFUN MABSBIGFLOAT (L) 
  1087.        (PROG (R) 
  1088.          (SETQ R (BIGFLOATP (CAR L)))
  1089.          (RETURN (COND ((NULL R) (LIST '(MABS) (CAR L)))
  1090.                (T (BCONS (FPABS (CDR R)))))))) 
  1091.  
  1092. #-NIL
  1093. (eval-when (load)
  1094. (FPPREC1 NIL $FPPREC)  ; Set up user's precision
  1095. )
  1096.  
  1097. ; Undeclarations for the file:
  1098. (declare-top (NOTYPE I N EXTRADIGS))
  1099.