/[cmucl]/src/code/irrat-dd.lisp
ViewVC logotype

Diff of /src/code/irrat-dd.lisp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.18 by rtoy, Wed Mar 18 01:24:52 2009 UTC revision 1.18.14.1 by rtoy, Thu Feb 25 20:34:50 2010 UTC
# Line 19  Line 19 
19    
20  (in-package "KERNEL")  (in-package "KERNEL")
21    
22    (intl:textdomain "cmucl")
23    
24  ;;;; Random constants, utility functions.  ;;;; Random constants, utility functions.
25    
26  (defconstant max-log  (defconstant max-log
27    7.0978271289338399678773454114191w2    7.0978271289338399678773454114191w2
28    "log(most-positive-double-double-float)")    _N"log(most-positive-double-double-float)")
29    
30  (defconstant min-log  (defconstant min-log
31    -7.4444007192138126231410729844608w2    -7.4444007192138126231410729844608w2
32    "log(least-positive-double-double-float")    _N"log(least-positive-double-double-float")
33    
34    
35  (defconstant loge2  (defconstant loge2
36    0.6931471805599453094172321214581765680755001w0    0.6931471805599453094172321214581765680755001w0
37    "log(2)")    _N"log(2)")
38    
39  (defconstant log2e  (defconstant log2e
40    1.442695040888963407359924681001892137426646w0    1.442695040888963407359924681001892137426646w0
41    "Log base 2 of e")    _N"Log base 2 of e")
42    
43  (defconstant log2ea  (defconstant log2ea
44    4.4269504088896340735992468100189213742664595w-1    4.4269504088896340735992468100189213742664595w-1
45    "log2(e)-1")    _N"log2(e)-1")
46    
47  (defconstant dd-pi  (defconstant dd-pi
48    3.141592653589793238462643383279502884197169w0    3.141592653589793238462643383279502884197169w0
49    "Pi")    _N"Pi")
50    
51  (defconstant dd-pi/2  (defconstant dd-pi/2
52    1.570796326794896619231321691639751442098585w0    1.570796326794896619231321691639751442098585w0
53    "Pi/2")    _N"Pi/2")
54    
55  (defconstant dd-pi/4  (defconstant dd-pi/4
56    0.7853981633974483096156608458198757210492923w0    0.7853981633974483096156608458198757210492923w0
57    "Pi/4")    _N"Pi/4")
58    
59  ;; log2-c1 and log-c2 are log(2) arranged in such a way that log2-c1 +  ;; log2-c1 and log-c2 are log(2) arranged in such a way that log2-c1 +
60  ;; log2-c2 is log(2) to an accuracy greater than double-double-float.  ;; log2-c2 is log(2) to an accuracy greater than double-double-float.
# Line 64  Line 66 
66    
67  (defconstant sqrt-1/2  (defconstant sqrt-1/2
68    0.7071067811865475244008443621048490392848w0    0.7071067811865475244008443621048490392848w0
69    "Sqrt(1/2)")    _N"Sqrt(1/2)")
70    
71  ;; Evaluate polynomial  ;; Evaluate polynomial
72  (declaim (maybe-inline poly-eval poly-eval-1))  (declaim (maybe-inline poly-eval poly-eval-1))
# Line 144  Line 146 
146    ;; log(2)/2, where the coefficients of P and Q are given Pn and Qn    ;; log(2)/2, where the coefficients of P and Q are given Pn and Qn
147    ;; above.  Theoretical peak relative error = 8.1e-36.    ;; above.  Theoretical peak relative error = 8.1e-36.
148    (defun dd-%expm1 (x)    (defun dd-%expm1 (x)
149      "exp(x) - 1"      _N"exp(x) - 1"
150      (declare (type double-double-float x)      (declare (type double-double-float x)
151               (optimize (speed 3) (space 0)               (optimize (speed 3) (space 0)
152                         (inhibit-warnings 3)))                         (inhibit-warnings 3)))
# Line 1206  pi/4    11001001000011111101101010100010 Line 1208  pi/4    11001001000011111101101010100010
1208                  #x91615E #xE61B08 #x659985 #x5F14A0 #x68408D #xFFD880                  #x91615E #xE61B08 #x659985 #x5F14A0 #x68408D #xFFD880
1209                  #x4D7327 #x310606 #x1556CA #x73A8C9 #x60E27B #xC08C6B                  #x4D7327 #x310606 #x1556CA #x73A8C9 #x60E27B #xC08C6B
1210                  ))                  ))
1211    "396 (hex) digits of 2/pi")    _N"396 (hex) digits of 2/pi")
1212    
1213    
1214  (let ((y (make-array 3 :element-type 'double-float))  (let ((y (make-array 3 :element-type 'double-float))
# Line 1506  pi/4    11001001000011111101101010100010 Line 1508  pi/4    11001001000011111101101010100010
1508                 (setf s (* loge2 e))))                 (setf s (* loge2 e))))
1509          (when (> s max-log)          (when (> s max-log)
1510            ;; Overflow.  What to do?            ;; Overflow.  What to do?
1511            (error "Overflow"))            (error _"Overflow"))
1512          (when (< s min-log)          (when (< s min-log)
1513            (return-from dd-%powil 0w0))            (return-from dd-%powil 0w0))
1514    
# Line 1638  pi/4    11001001000011111101101010100010 Line 1640  pi/4    11001001000011111101101010100010
1640                 (values rho 0)))))))                 (values rho 0)))))))
1641    
1642  (defun dd-complex-sqrt (z)  (defun dd-complex-sqrt (z)
1643    "Principle square root of Z    _N"Principle square root of Z
1644    
1645  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1646    (declare (number z))    (declare (number z))
# Line 1680  Z may be any number, but the result is a Line 1682  Z may be any number, but the result is a
1682          (complex eta nu)))))          (complex eta nu)))))
1683    
1684  (defun dd-complex-log-scaled (z j)  (defun dd-complex-log-scaled (z j)
1685    "Compute log(2^j*z).    _N"Compute log(2^j*z).
1686    
1687  This is for use with J /= 0 only when |z| is huge."  This is for use with J /= 0 only when |z| is huge."
1688    (declare (number z)    (declare (number z)
# Line 1715  This is for use with J /= 0 only when |z Line 1717  This is for use with J /= 0 only when |z
1717                   (atan y x))))))                   (atan y x))))))
1718    
1719  (defun dd-complex-log (z)  (defun dd-complex-log (z)
1720    "Log of Z = log |Z| + i * arg Z    _N"Log of Z = log |Z| + i * arg Z
1721    
1722  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1723    (declare (number z))    (declare (number z))
# Line 1727  Z may be any number, but the result is a Line 1729  Z may be any number, but the result is a
1729  ;; never 0 since we have positive and negative zeroes.  ;; never 0 since we have positive and negative zeroes.
1730    
1731  (defun dd-complex-atanh (z)  (defun dd-complex-atanh (z)
1732    "Compute atanh z = (log(1+z) - log(1-z))/2"    _N"Compute atanh z = (log(1+z) - log(1-z))/2"
1733    (declare (number z))    (declare (number z))
1734    (cond ((and (realp z) (< z -1))    (cond ((and (realp z) (< z -1))
1735           ;; ATANH is continuous with quadrant III in this case.           ;; ATANH is continuous with quadrant III in this case.
# Line 1788  Z may be any number, but the result is a Line 1790  Z may be any number, but the result is a
1790                        (- (* beta nu))))))))                        (- (* beta nu))))))))
1791    
1792  (defun dd-complex-tanh (z)  (defun dd-complex-tanh (z)
1793    "Compute tanh z = sinh z / cosh z"    _N"Compute tanh z = sinh z / cosh z"
1794    (declare (number z))    (declare (number z))
1795    (let ((x (float (realpart z) 1.0w0))    (let ((x (float (realpart z) 1.0w0))
1796          (y (float (imagpart z) 1.0w0)))          (y (float (imagpart z) 1.0w0)))
# Line 1868  Z may be any number, but the result is a Line 1870  Z may be any number, but the result is a
1870    (complex (+ (realpart z) 1) (imagpart z)))    (complex (+ (realpart z) 1) (imagpart z)))
1871    
1872  (defun dd-complex-acos (z)  (defun dd-complex-acos (z)
1873    "Compute acos z = pi/2 - asin z    _N"Compute acos z = pi/2 - asin z
1874    
1875  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1876    (declare (number z))    (declare (number z))
# Line 1894  Z may be any number, but the result is a Line 1896  Z may be any number, but the result is a
1896                                              sqrt-1-z)))))))))                                              sqrt-1-z)))))))))
1897    
1898  (defun dd-complex-acosh (z)  (defun dd-complex-acosh (z)
1899    "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))    _N"Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1900    
1901  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1902    (declare (number z))    (declare (number z))
# Line 1920  Z may be any number, but the result is a Line 1922  Z may be any number, but the result is a
1922    
1923    
1924  (defun dd-complex-asin (z)  (defun dd-complex-asin (z)
1925    "Compute asin z = asinh(i*z)/i    _N"Compute asin z = asinh(i*z)/i
1926    
1927  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1928    (declare (number z))    (declare (number z))
# Line 1949  Z may be any number, but the result is a Line 1951  Z may be any number, but the result is a
1951                                                   sqrt-1+z)))))))))))                                                   sqrt-1+z)))))))))))
1952    
1953  (defun dd-complex-asinh (z)  (defun dd-complex-asinh (z)
1954    "Compute asinh z = log(z + sqrt(1 + z*z))    _N"Compute asinh z = log(z + sqrt(1 + z*z))
1955    
1956  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1957    (declare (number z))    (declare (number z))
# Line 1960  Z may be any number, but the result is a Line 1962  Z may be any number, but the result is a
1962               (- (realpart result)))))               (- (realpart result)))))
1963    
1964  (defun dd-complex-atan (z)  (defun dd-complex-atan (z)
1965    "Compute atan z = atanh (i*z) / i    _N"Compute atan z = atanh (i*z) / i
1966    
1967  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1968    (declare (number z))    (declare (number z))
# Line 1971  Z may be any number, but the result is a Line 1973  Z may be any number, but the result is a
1973               (- (realpart result)))))               (- (realpart result)))))
1974    
1975  (defun dd-complex-tan (z)  (defun dd-complex-tan (z)
1976    "Compute tan z = -i * tanh(i * z)    _N"Compute tan z = -i * tanh(i * z)
1977    
1978  Z may be any number, but the result is always a complex."  Z may be any number, but the result is always a complex."
1979    (declare (number z))    (declare (number z))

Legend:
Removed from v.1.18  
changed lines
  Added in v.1.18.14.1

  ViewVC Help
Powered by ViewVC 1.1.5