# Ticket #1059: math-patches.lisp

File math-patches.lisp, 6.4 KB (added by dfindlay, 2 years ago) |
---|

Line | |
---|---|

1 | ;;; Math Patches |

2 | |

3 | ;;; phase (l0-numbers.lisp) |

4 | ;;; Returns wrong value for negative zero |

5 | ;;; - use %double-float-sign, %short-float-sign instead of minusp |

6 | |

7 | (defun phase (number) |

8 | "Return the angle part of the polar representation of a complex number. |

9 | For complex numbers, this is (atan (imagpart number) (realpart number)). |

10 | For non-complex positive numbers, this is 0. For non-complex negative |

11 | numbers this is PI." |

12 | (number-case number |

13 | (rational |

14 | (if (minusp number) |

15 | (%short-float pi) |

16 | 0.0f0)) |

17 | (double-float |

18 | |

19 | ;;; Changed code (1) |

20 | |

21 | (if (%double-float-sign number) |

22 | |

23 | ;;; End of change (1) |

24 | |

25 | (%double-float pi) |

26 | 0.0d0)) |

27 | (complex |

28 | (atan (%imagpart number) (%realpart number))) |

29 | (short-float |

30 | |

31 | ;;; Changed code (2) |

32 | |

33 | (if (%short-float-sign number) |

34 | |

35 | ;;; End of change (2) |

36 | |

37 | (%short-float pi) |

38 | 0.0s0)))) |

39 | |

40 | ;;; exp (l0-float.lisp) |

41 | ;;; Allow ridiculously large negative real rationals |

42 | ;;; - these just need to return 0.0 |

43 | |

44 | (defun exp (x) |

45 | "Return e raised to the power NUMBER." |

46 | (typecase x |

47 | (complex (* (exp (realpart x)) (cis (imagpart x)))) |

48 | (double-float (%double-float-exp! x (%make-dfloat))) |

49 | (t |

50 | |

51 | ;;; Changed code |

52 | |

53 | (if (and (typep x 'rational) |

54 | (< x -104)) |

55 | 0.0s0 |

56 | #+32-bit-target |

57 | (target::with-stack-short-floats ((sx x)) |

58 | (%single-float-exp! sx (%make-sfloat))) |

59 | #+64-bit-target |

60 | (%single-float-exp (%short-float x)))))) |

61 | ;;; End of change |

62 | |

63 | ;;; tan (l0-float.lisp) |

64 | ;;; Complex numbers with very large imaginary parts were overflowing in (* 2 (abs i)) |

65 | ;;; - recode to avoid this |

66 | |

67 | (defun tan (x) |

68 | "Return the tangent of NUMBER." |

69 | (cond ((complexp x) |

70 | (let ((r (realpart x)) |

71 | (i (imagpart x))) |

72 | (if (zerop i) |

73 | (complex (tan r) i) |

74 | (let* ((tx (tan r)) |

75 | (ty (tanh i)) |

76 | (tx2 (* tx tx)) |

77 | (d (1+ (* tx2 (* ty ty)))) |

78 | |

79 | ;;; Changed code |

80 | |

81 | (c (if (> (abs i) 20) |

82 | (* 2 (exp (- (abs i)))) |

83 | (/ (cosh i))))) |

84 | (complex (/ (* (* c c) tx) d) |

85 | (/ (* ty (1+ tx2)) d)))))) |

86 | |

87 | ;;; End of change |

88 | |

89 | ((or (typep x 'ratio) |

90 | (> (abs x) two^23)) |

91 | (let ((c (%extended-cis x))) |

92 | (if (typep x 'double-float) |

93 | (/ (imagpart c) (realpart c)) |

94 | (%short-float (/ (imagpart c) (realpart c)))))) |

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

96 | (%double-float-tan! x (%make-dfloat))) |

97 | (t |

98 | #+32-bit-target |

99 | (target::with-stack-short-floats ((sx x)) |

100 | (%single-float-tan! sx (%make-sfloat))) |

101 | #+64-bit-target |

102 | (%single-float-tan (%short-float x)) |

103 | ))) |

104 | |

105 | ;;; %complex-atanh (l0-float.lisp) |

106 | ;;; The expression (/ 2 y) was overflowing for very small y |

107 | ;;; - recode log(2/y) as log(2) - log(y) |

108 | ;;; replace atan(2/y, -1) with atan(2, -y) of atan(-2, y) depending on sign of y |

109 | ;;; Argument to log1+ becoming -1 when x is near 1 and y is small |

110 | ;;; - add additional cond clause to handle this |

111 | |

112 | (defun %complex-atanh (z) |

113 | (let* ((rx (realpart z)) |

114 | (ix (imagpart z)) |

115 | (sign (typecase rx |

116 | (double-float (%double-float-sign rx)) |

117 | (short-float (%short-float-sign rx)) |

118 | (t (minusp rx)))) |

119 | (x rx) |

120 | (y ix) |

121 | (y1 (abs y)) |

122 | ra ia) |

123 | ;; following code requires non-negative x |

124 | (when sign |

125 | (setf x (- x)) |

126 | (setf y (- y))) |

127 | (cond ((> (max x y1) 1.8014399e+16) |

128 | ;; large value escape |

129 | (setq ra (if (> x y1) |

130 | (let ((r (/ y x))) |

131 | (/ (/ x) (1+ (* r r)))) |

132 | (let ((r (/ x y))) |

133 | (/ (/ r y) (1+ (* r r)))))) |

134 | (setq ia (typecase y |

135 | (double-float (float-sign y double-float-half-pi)) |

136 | (single-float (float-sign y single-float-half-pi)) |

137 | (t (if (minusp y) #.(- single-float-half-pi) single-float-half-pi))))) |

138 | ((= x 1) |

139 | |

140 | ;;; Changed code |

141 | |

142 | (cond ((< y1 1e-9) |

143 | (setq ra (/ (- (if (typep y 'double-float) double-float-log2 single-float-log2) |

144 | (log-e y1)) |

145 | 2)) |

146 | (setq ia (/ (if (minusp y) (atan -2 y) (atan 2 (- y))) 2))) |

147 | (t |

148 | (setq ra (/ (log1+ (/ 4 (* y y))) 4)) |

149 | (setq ia (/ (atan (/ 2 y) -1) 2))))) |

150 | ((and (< y1 1) |

151 | (< 0.5 x 2)) |

152 | (let ((x-1 (- x 1)) |

153 | (x+1 (+ x 1)) |

154 | (y2 (* y y))) |

155 | (setq ra (/ (log-e (/ (+ (* x-1 x-1) y2) (+ (* x+1 x+1) y2))) -4)) |

156 | (setq ia (/ (atan (* 2 y) (- 1 (+ (* x x) y2))) 2)))) |

157 | |

158 | ;;; End of change |

159 | |

160 | (t |

161 | (let ((r2 (+ (* x x) (* y y)))) |

162 | (setq ra (/ (log1+ (/ (* -4 x) (1+ (+ (* 2 x) r2)))) -4)) |

163 | (setq ia (/ (atan (* 2 y) (- 1 r2)) 2))))) |

164 | ;; fixup signs, with special case for real arguments |

165 | (cond (sign |

166 | (setq ra (- ra)) |

167 | (when (typep z 'complex) |

168 | (setq ia (- ia)))) |

169 | (t |

170 | (unless (typep z 'complex) |

171 | (setq ia (- ia))))) |

172 | (complex ra ia))) |

173 | |

174 | ;;; tanh (numbers.lisp) |

175 | ;;; Complex numbers with very large real parts were overflowing in (* 2 (abs r)) |

176 | ;;; - recode to avoid this |

177 | ;;; (analogous fix to that in tan) |

178 | |

179 | (defun tanh (x) |

180 | "Return the hyperbolic tangent of NUMBER." |

181 | (cond ((complexp x) |

182 | (let ((r (realpart x)) |

183 | (i (imagpart x))) |

184 | (if (zerop r) |

185 | (complex r (tan i)) |

186 | (let* ((tx (tanh r)) |

187 | (ty (tan i)) |

188 | (ty2 (* ty ty)) |

189 | (d (1+ (* (* tx tx) ty2))) |

190 | |

191 | ;;; Changed code |

192 | |

193 | (c (if (> (abs r) 20) |

194 | (* 2 (exp (- (abs r)))) |

195 | (/ (cosh r))))) |

196 | (complex (/ (* tx (1+ ty2)) d) |

197 | (/ (* (* c c) ty) d)))))) |

198 | |

199 | ;;; End of change |

200 | |

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

202 | (%double-float-tanh! x (%make-dfloat))) |

203 | ((and (typep x 'rational) |

204 | (> (abs x) 12)) |

205 | (if (plusp x) 1.0s0 -1.0s0)) |

206 | (t |

207 | #+32-bit-target |

208 | (target::with-stack-short-floats ((sx x)) |

209 | (%single-float-tanh! sx (%make-sfloat))) |

210 | #+64-bit-target |

211 | (%single-float-tanh (%short-float x))))) |

212 |