# Ticket #869: fn-ticket.lisp

File fn-ticket.lisp, 9.7 KB (added by dfindlay, 5 years ago) |
---|

Line | |
---|---|

1 | (in-package :cl-user) |

2 | |

3 | ;;; This code is written in cl-user, not ccl |

4 | ;;; The code for new-sqrt, new-log and new-atan follows that for CCL's sqrt, log |

5 | ;;; and atan (in l0-float.lisp) in terms of the case structure; but where the branch |

6 | ;;; is not changed the CCL function is called (rather than copying the code in CCL's |

7 | ;;; branch, which probably wouldn't work here). I hope it shouldn't be too difficult |

8 | ;;; to cut and paste the appropriate code bits into CCL/s code. |

9 | |

10 | ;;; Use of new %hypot is assumed: |

11 | ;;; (but this version isn't as efficient as the *real* one) |

12 | (defun ccl::%hypot (x y) |

13 | (let ((a (abs x)) |

14 | (b (abs y))) |

15 | (when (> a b) |

16 | (psetq a b b a)) |

17 | (if (= b 0.0d0) |

18 | b |

19 | (let ((c (/ a b))) |

20 | (* b (sqrt (+ 1.0d0 (* c c)))))))) |

21 | |

22 | |

23 | (defmacro single-float-sqrt (x) |

24 | ;; usual expression for single float sqrt |

25 | ;; (referred to several times in new-sqrt) |

26 | `(sqrt ,x)) |

27 | |

28 | |

29 | ;; fix sqrt to return values for rationals whose square root can be represented |

30 | ;; as single float, even though they can't be themselves. |

31 | ;; also take care to avoid intermediate overflow for some complex arguments |

32 | (defun new-sqrt (x) |

33 | "Return the square root of NUMBER." |

34 | (cond ((zerop x) x) |

35 | ((complexp x) |

36 | ;; take care that intermediate results do not overflow/underflow: |

37 | ;; pre-scale argument by an appropriate power of two |

38 | ;; (there may be more efficient ways of doing this; and strictly |

39 | ;; if doesn't need to be done for 'reasonably sized' arguments) |

40 | (let ((m (max (abs (realpart x)) |

41 | (abs (imagpart x)))) |

42 | (s 1) |

43 | (s2 1)) |

44 | (typecase m |

45 | (integer (let ((expon (floor (integer-length m) 2))) |

46 | (setq s (ash 1 expon)) |

47 | (setq s2 (ash s expon)))) |

48 | (ratio (let ((expon (floor (- (integer-length (numerator m)) |

49 | (integer-length (denominator m))) |

50 | 2))) |

51 | (setq s (expt 2 expon)) |

52 | (setq s2 (* s s)))) |

53 | (float (let ((expon (truncate (nth-value 1 (decode-float m)) 2))) |

54 | (if (minusp expon) |

55 | (incf expon) |

56 | (decf expon)) |

57 | (setq s (scale-float (float 1 m) expon)) |

58 | (setq s2 (scale-float s expon))))) |

59 | (let ((z (/ x s2))) |

60 | (* s (* (sqrt (abs z)) (cis (/ (phase z) 2))))))) |

61 | ((minusp x) (complex 0 (new-sqrt (- x)))) |

62 | ((floatp x) |

63 | ;; usual floating point sqrt |

64 | (sqrt x)) |

65 | ((integerp x) |

66 | (let* ((a (isqrt x)) |

67 | (test (* a a))) |

68 | (cond ((eql x test) |

69 | a) |

70 | ((<= x most-positive-single-float) |

71 | (single-float-sqrt x)) |

72 | (t |

73 | ;; handle case where (sqrt x) can be represented as single float |

74 | ;; but x can't |

75 | (let ((error (/ x test))) |

76 | (* (float a 1.0) |

77 | (single-float-sqrt error))))))) |

78 | ((ratiop x) |

79 | ;; similar to integer case - be careful when (sqrt x) can be represented |

80 | ;; as single float but x can't |

81 | (let* ((n (numerator x)) |

82 | (a (isqrt n)) |

83 | (test1 (* a a))) |

84 | (unless (eql n test1) |

85 | (when (and (<= x most-positive-single-float) |

86 | (>= x least-positive-normalized-single-float)) |

87 | (return-from new-sqrt (single-float-sqrt x)))) |

88 | (let* ((d (denominator x)) |

89 | (b (isqrt d)) |

90 | (test2 (* b b))) |

91 | (cond ((and (eql d test2) |

92 | (eql n test1)) |

93 | (/ a b)) |

94 | ((and (<= x most-positive-single-float) |

95 | (>= x least-positive-normalized-single-float)) |

96 | (single-float-sqrt x)) |

97 | (t |

98 | (let ((error (/ x (/ test1 test2)))) |

99 | (* (float (/ a b) 1.0) |

100 | (single-float-sqrt error)))))))) |

101 | (t |

102 | (single-float-sqrt x)))) |

103 | |

104 | |

105 | ;; makes log work for small rationals (it already worked for big ones) |

106 | ;; correct wrong return type for negative bignums |

107 | ;; fix float/rational contagion in two argument case |

108 | ;; propagate rational behaviour to complex rationals as well |

109 | (defun new-log (x &optional (b nil b-p)) |

110 | "Return the logarithm of NUMBER in the base BASE, which defaults to e." |

111 | (if b-p |

112 | (if (zerop b) |

113 | (if (zerop x) |

114 | ;(report-bad-arg x '(not (satisfies zerop) )) |

115 | (error "NEW-LOG: zero arguments") |

116 | ;(if (floatp x) (float 0.0d0 x) 0)) ** CORRECT THE CONTAGION for complex args ** |

117 | ; (is this the best way to get a zero of the right type?) |

118 | (+ 0 (* x b))) |

119 | ;; ** Don't we have to do the float/rational contagion before the division? |

120 | ;; (is this the best way to do it?) |

121 | (/ (new-log-e (+ x (* b 0))) |

122 | (new-log-e (+ b (* x 0))))) |

123 | (new-log-e x))) |

124 | |

125 | |

126 | ;;; constant required: will need to be 'magic' number somewhere |

127 | (defconstant +log2+ #.(log 2.0d0)) ; = 0.6931471805599453D0 |

128 | |

129 | (defun new-log-e (x) |

130 | (cond |

131 | ((bignump x) |

132 | (if (minusp x) |

133 | (complex (new-log-e (- x)) #.(coerce pi 'single-float)) ; ** FIX RETURN TYPE ERROR ** |

134 | (let* ((base1 3) |

135 | (guess (floor (1- (integer-length x)) |

136 | (log base1 2))) |

137 | (guess1 (* guess (new-log-e base1)))) |

138 | (+ guess1 (new-log-e (/ x (expt base1 guess))))))) |

139 | ((and (ratiop x) ; Escape out small values as well as large |

140 | (if (minusp x) |

141 | (or (> x least-negative-normalized-short-float) |

142 | (< x most-negative-short-float)) |

143 | (or (> x most-positive-short-float) |

144 | (< x least-positive-normalized-short-float)))) |

145 | (- (new-log-e (numerator x)) (new-log-e (denominator x)))) |

146 | ((typep x 'complex) |

147 | ;; take care that intermediate results do not overflow/underflow: |

148 | ;; pre-scale argument by an appropriate power of two |

149 | ;; we only need to scale for very large and very small values - |

150 | ;; hence the various 'magic' numbers (values may not be too |

151 | ;; critical but do depend on the sqrt update to fix abs's operation) |

152 | (let ((m (max (abs (realpart x)) |

153 | (abs (imagpart x)))) |

154 | (log-s 0) |

155 | (s 1)) |

156 | (typecase m |

157 | (integer (let ((expon (integer-length m))) |

158 | (when (> expon 126) |

159 | (setq log-s (float (* expon +log2+) 1.0)) |

160 | (setq s (ash 1 expon))))) |

161 | (ratio (let ((expon (- (integer-length (numerator m)) |

162 | (integer-length (denominator m))))) |

163 | (when (or (> expon 125) |

164 | (< expon -122)) |

165 | (setq log-s (float (* expon +log2+) 1.0)) |

166 | (setq s (expt 2 expon))))) |

167 | (single-float (let ((expon (nth-value 1 (decode-float m)))) |

168 | (when (or (> expon 126) |

169 | (< expon -124)) |

170 | (if (minusp expon) |

171 | (incf expon) |

172 | (decf expon)) |

173 | (setq log-s (float (* expon +log2+) 1.0)) ; assumes (float-radix m) = 2 |

174 | (setq s (scale-float 1.0 expon))))) |

175 | (double-float (let ((expon (nth-value 1 (decode-float m)))) |

176 | (when (or (> expon 1022) |

177 | (< expon -1020)) |

178 | (if (minusp expon) |

179 | (incf expon) |

180 | (decf expon)) |

181 | (setq log-s (* expon +log2+)) ; assumes (float-radix m) = 2 |

182 | (setq s (scale-float 1.0d0 expon)))))) |

183 | (if (eql s 1) |

184 | (complex (new-log-e (abs x)) (phase x)) |

185 | (let ((z (/ x s))) |

186 | (complex (+ log-s (new-log-e (abs z))) (phase z)))))) |

187 | ((typep x 'double-float) |

188 | ;; usual double float log-e |

189 | (log x)) |

190 | (t |

191 | ;; usual single float log-e |

192 | (log x) |

193 | ))) |

194 | |

195 | |

196 | ;; atan needs to be more careful about small and large rational arguments |

197 | ;; (this will make phase as careful too for complex rationals) |

198 | (defun new-atan (y &optional (x nil x-p)) |

199 | "Return the arc tangent of Y if X is omitted or Y/X if X is supplied." |

200 | (cond (x-p |

201 | (when (or (rationalp x) (rationalp y)) |

202 | (unless (zerop y) ; don't tread on special zero handling |

203 | ;; rescale arguments so that the maximum absolute value is 1 |

204 | (cond ((> (abs y) (abs x)) |

205 | (setf x (/ x (abs y))) |

206 | (setf y (signum y))) |

207 | (t |

208 | (setf y (/ y (abs x))) |

209 | (setf x (signum x)))))) |

210 | ;; usual expression |

211 | (atan y x)) |

212 | (t |

213 | ;; usual expressions |

214 | ;; [but (sqrt -1) in the complex case could be just #C(0 1)] |

215 | (atan y)))) |

216 | |

217 | |

218 | |

219 | |

220 | ;;; FIXES |

221 | #| |

222 | |

223 | (sqrt (expt 10 47)) => 3.1622778E+23 |

224 | (sqrt (/ (expt 10 47) 3)) => 1.8257418E+23 |

225 | (sqrt (complex (expt 10 46) (expt 10 47))) => #C(2.3505187E+23 2.12719E+23) |

226 | (sqrt (complex most-positive-short-float most-positive-short-float)) => #C(2.0267142E+19 8.394926E+18) |

227 | |

228 | (log (expt 10 -66)) => -151.97063 |

229 | (log (- (expt 10 66))) => #C(151.97063 3.1415927) [corrects type] |

230 | (log (complex (expt 10 65) (expt 10 66))) => #C(151.9756 1.4711276) |

231 | (log (complex (expt 10 -65) (expt 10 -66))) => #C(-149.66307 0.09966865) |

232 | (log 8.0d0 2) => 3.0D0 [corrects accuracy] |

233 | (log #C(0.0 1.0) 0) => #C(0.0 0.0) [corrects type] |

234 | |

235 | (atan (expt 10 46) (expt 10 47)) => 0.09966865 |

236 | (atan (expt 10 -46) (expt 10 -47)) => 1.4711276 [corrects loss of accuracy] |

237 | |

238 | |# |

239 |