ccl-format.patch 4.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129
  1. diff -r -U1 ccl.orig/lib/format.lisp ccl/lib/format.lisp
  2. --- ccl.orig/lib/format.lisp 2015-11-07 02:10:10.000000000 +0600
  3. +++ ccl/lib/format.lisp 2015-11-20 22:51:51.736191995 +0600
  4. @@ -1296,5 +1296,2 @@
  5. -
  6. -
  7. -
  8. ;;; Given a non-negative floating point number, SCALE-EXPONENT returns a
  9. @@ -1305,41 +1302,74 @@
  10. -
  11. -(defconstant long-log10-of-2 0.30103d0)
  12. -
  13. -#|
  14. -(defun scale-exponent (x)
  15. - (if (floatp x )
  16. - (scale-expt-aux (abs x) 0.0d0 1.0d0 1.0d1 1.0d-1 long-log10-of-2)
  17. - (report-bad-arg x 'float)))
  18. -
  19. -#|this is the slisp code that was in the place of the error call above.
  20. - before floatp was put in place of shortfloatp.
  21. - ;(scale-expt-aux x (%sp-l-float 0) (%sp-l-float 1) %long-float-ten
  22. - ; %long-float-one-tenth long-log10-of-2)))
  23. -|#
  24. -
  25. -; this dies with floating point overflow (?) if fed least-positive-double-float
  26. -
  27. -(defun scale-expt-aux (x zero one ten one-tenth log10-of-2)
  28. - (let ((exponent (nth-value 1 (decode-float x))))
  29. - (if (= x zero)
  30. - (values zero 1)
  31. - (let* ((e (round (* exponent log10-of-2)))
  32. - (x (if (minusp e) ;For the end ranges.
  33. - (* x ten (expt ten (- -1 e)))
  34. - (/ x ten (expt ten (1- e))))))
  35. - (do ((d ten (* d ten))
  36. - (y x (/ x d))
  37. - (e e (1+ e)))
  38. - ((< y one)
  39. - (do ((m ten (* m ten))
  40. - (z y (* z m))
  41. - (e e (1- e)))
  42. - ((>= z one-tenth) (values x e)))))))))
  43. -|#
  44. -
  45. -(defun scale-exponent (n)
  46. - (let ((exp (nth-value 1 (decode-float n))))
  47. - (values (round (* exp long-log10-of-2)))))
  48. -
  49. +(defconstant single-float-min-e
  50. + (nth-value 1 (decode-float least-positive-single-float)))
  51. +(defconstant double-float-min-e
  52. + (nth-value 1 (decode-float least-positive-double-float)))
  53. +
  54. +;;; Adapted from CMUCL.
  55. +
  56. +;; This is a modified version of the scale computation from Burger and
  57. +;; Dybvig's paper "Printing floating-point quickly and accurately."
  58. +;; We only want the exponent, so most things not needed for the
  59. +;; computation of the exponent have been removed. We also implemented
  60. +;; the floating-point log approximation given in Burger and Dybvig.
  61. +;; This is very noticeably faster for large and small numbers. It is
  62. +;; slower for intermediate sized numbers.
  63. +(defun accurate-scale-exponent (v)
  64. + (declare (type float v))
  65. + (if (zerop v)
  66. + 1
  67. + (let ((float-radix 2) ; b
  68. + (float-digits (float-digits v)) ; p
  69. + (min-e
  70. + (etypecase v
  71. + (single-float single-float-min-e)
  72. + (double-float double-float-min-e))))
  73. + (multiple-value-bind (f e)
  74. + (integer-decode-float v)
  75. + (let ( ;; FIXME: these even tests assume normal IEEE rounding
  76. + ;; mode. I wonder if we should cater for non-normal?
  77. + (high-ok (evenp f)))
  78. + ;; We only want the exponent here.
  79. + (labels ((flog (x)
  80. + (declare (type (float (0.0)) x))
  81. + (let ((xd (etypecase x
  82. + (single-float
  83. + (float x 1d0))
  84. + (double-float
  85. + x))))
  86. + (ceiling (- (the (double-float -400d0 400d0)
  87. + (log xd 10d0))
  88. + 1d-10))))
  89. + (fixup (r s m+ k)
  90. + (if (if high-ok
  91. + (>= (+ r m+) s)
  92. + (> (+ r m+) s))
  93. + (+ k 1)
  94. + k))
  95. + (scale (r s m+)
  96. + (let* ((est (flog v))
  97. + (scale (the integer (10-to-e (abs est)))))
  98. + (if (>= est 0)
  99. + (fixup r (* s scale) m+ est)
  100. + (fixup (* r scale) s (* m+ scale) est)))))
  101. + (let (r s m+)
  102. + (if (>= e 0)
  103. + (let* ((be (expt float-radix e))
  104. + (be1 (* be float-radix)))
  105. + (if (/= f (expt float-radix (1- float-digits)))
  106. + (setf r (* f be 2)
  107. + s 2
  108. + m+ be)
  109. + (setf r (* f be1 2)
  110. + s (* float-radix 2)
  111. + m+ be1)))
  112. + (if (or (= e min-e)
  113. + (/= f (expt float-radix (1- float-digits))))
  114. + (setf r (* f 2)
  115. + s (* (expt float-radix (- e)) 2)
  116. + m+ 1)
  117. + (setf r (* f float-radix 2)
  118. + s (* (expt float-radix (- 1 e)) 2)
  119. + m+ float-radix)))
  120. + (scale r s m+))))))))
  121. @@ -1922,3 +1952,3 @@
  122. (format-error "incompatible values for k and d")))
  123. - (when (not exp) (setq exp (scale-exponent number)))
  124. + (when (not exp) (setq exp (accurate-scale-exponent (abs number))))
  125. AGAIN