# Ticket #869: ccl-math-reference.lisp

File ccl-math-reference.lisp, 53.2 KB (added by dfindlay, 5 years ago) |
---|

Line | |
---|---|

1 | ;;; Reference math function implementations |

2 | ;;; ======================================= |

3 | ;;; Hopefully, all cases where the ANSI function definitions are "terrible methods |

4 | ;;; for floating point computation" have been addressed :-) |

5 | ;;; Care has been taken for rational (and complex rational) arguments to maintain |

6 | ;;; precision by not converting to (short) float too early. |

7 | ;;; |

8 | ;;; In most cases (exp and exp-like functions excepted) the common lisp functions |

9 | ;;; are only called for known floating point arguments. It is assumed that they |

10 | ;;; return accurate values in these cases! |

11 | ;;; |

12 | ;;; This lot should be applicable in any Lisp that (just) uses IEEE single and |

13 | ;;; double floating point formats. |

14 | ;;; (This code does assume that short-float = single-float) |

15 | ;;; |

16 | ;;; It should be "obvious" how this lot can be (back) ported to CCL... |

17 | ;;; (Replace calls to common lisp functions by #| ... |# commented-out blocks.) |

18 | |

19 | (cl:defpackage "MATH-REFERENCE" |

20 | (:use "COMMON-LISP") |

21 | (:nicknames "MATH") |

22 | (:shadow "SIGNUM" "ABS" "PHASE" "CIS" "SIN" "COS" "TAN" "ASIN" "ACOS" "ATAN" |

23 | "LOG" "EXP" "EXPT" "SQRT" "SINH" "COSH" "TANH" "ASINH" "ACOSH" "ATANH") |

24 | (:export "SIGNUM" "ABS" "PHASE" "CIS" "SIN" "COS" "TAN" "ASIN" "ACOS" "ATAN" |

25 | "LOG" "EXP" "EXPT" "SQRT" "SINH" "COSH" "TANH" "ASINH" "ACOSH" "ATANH")) |

26 | |

27 | (in-package :math) |

28 | |

29 | ;;; === Stubs for low level CCL functions |

30 | |

31 | (defmacro number-case (keyform &body clauses) |

32 | (let ((value (gensym))) |

33 | `(let ((,value ,keyform)) |

34 | (typecase ,value |

35 | ,@clauses |

36 | (otherwise |

37 | (error "value ~S is not of the expected type NUMBER." ,value)))))) |

38 | |

39 | (defmacro %realpart (x) |

40 | `(realpart ,x)) |

41 | |

42 | (defmacro %imagpart (x) |

43 | `(imagpart ,x)) |

44 | |

45 | (defmacro %numerator (x) |

46 | `(numerator ,x)) |

47 | |

48 | (defmacro %denominator (x) |

49 | `(denominator ,x)) |

50 | |

51 | (defmacro %short-float (x) |

52 | `(float ,x 1s0)) |

53 | |

54 | (defmacro %double-float (x) |

55 | `(float ,x 1d0)) |

56 | |

57 | (defmacro with-stack-double-floats (vars &body body) |

58 | `(let ,(mapcar (lambda (item) |

59 | `(,(first item) (%double-float ,(second item)))) |

60 | vars) |

61 | ,@body)) |

62 | |

63 | (defmacro report-bad-arg (arg type) |

64 | `(error "~S: not of type ~S" ,arg ,type)) |

65 | |

66 | (defmacro ratiop (x) |

67 | `(typep ,x 'ratio)) |

68 | |

69 | (defmacro bignump (x) |

70 | `(typep ,x 'bignum)) |

71 | |

72 | (defmacro fsqrt (x) |

73 | `(cl:sqrt ,x)) |

74 | |

75 | |

76 | (defconstant IEEE-single-float-bias 126) |

77 | (defconstant IEEE-double-float-bias 1022) |

78 | |

79 | |

80 | (defun %double-float-sign (x) |

81 | (or (minusp x) (eql x -0.0d0))) |

82 | |

83 | (defun %short-float-sign (x) |

84 | (or (minusp x) (eql x -0.0s0))) |

85 | |

86 | ;;; These next two aren't quite accurate, but I don't think that matters :-) |

87 | (defun %short-float-exp (x) |

88 | (+ (nth-value 1 (decode-float x)) IEEE-single-float-bias)) |

89 | |

90 | (defun %double-float-exp (x) |

91 | (+ (nth-value 1 (decode-float x)) IEEE-double-float-bias)) |

92 | |

93 | |

94 | (defun %integer-power (b e) |

95 | (declare (type unsigned-byte e)) |

96 | (if (zerop e) |

97 | (+ 1 (* b 0)) |

98 | (if (eql b 2) |

99 | (ash 1 e) |

100 | (do* ((next (ash e -1) (ash e -1)) |

101 | (oddexp (oddp e) (oddp e)) |

102 | (total (if oddexp b 1) (if oddexp (* b total) total))) |

103 | ((zerop next) total) |

104 | (declare (type unsigned-byte next)) |

105 | (setq b (* b b) e next))))) |

106 | |

107 | (defun %hypot (x y) |

108 | (let ((a (abs x)) |

109 | (b (abs y))) |

110 | (when (> a b) |

111 | (psetq a b b a)) |

112 | (if (= b 0.0d0) |

113 | 0.0d0 |

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

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

116 | |

117 | |

118 | ;;; === Unmodified CCL function definitions |

119 | ;;; (Copies needed to pick up shadowed versions of called functions) |

120 | |

121 | (defun abs (number) |

122 | "Return the absolute value of the number." |

123 | (number-case number |

124 | (real (cl:abs number)) |

125 | #| |

126 | (fixnum |

127 | (locally (declare (fixnum number)) |

128 | (if (minusp number) (- number) number))) |

129 | (double-float |

130 | (%double-float-abs number)) |

131 | (short-float |

132 | (%short-float-abs number)) |

133 | (bignum |

134 | (if (bignum-minusp number)(negate-bignum number) number)) |

135 | (ratio |

136 | (if (minusp number) (- number) number)) |

137 | |# |

138 | (complex |

139 | (let ((rx (%realpart number)) |

140 | (ix (%imagpart number))) |

141 | (number-case rx |

142 | (rational |

143 | (sqrt (+ (* rx rx) (* ix ix)))) |

144 | (short-float |

145 | (%short-float (%hypot (%double-float rx) |

146 | (%double-float ix)))) |

147 | (double-float |

148 | (%hypot rx ix))))))) |

149 | |

150 | (defun phase (number) |

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

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

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

154 | numbers this is PI." |

155 | (number-case number |

156 | (rational |

157 | (if (minusp number) |

158 | (%short-float pi) |

159 | 0.0f0)) |

160 | (double-float |

161 | (if (minusp number) |

162 | (%double-float pi) |

163 | 0.0d0)) |

164 | (complex |

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

166 | (short-float |

167 | (if (minusp number) |

168 | (%short-float pi) |

169 | 0.0s0)))) |

170 | |

171 | (defun exp (x) |

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

173 | (typecase x |

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

175 | #| |

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

177 | (t |

178 | #+32-bit-target |

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

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

181 | #+64-bit-target |

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

183 | |# |

184 | (t (cl:exp x)))) |

185 | |

186 | ;;; === New function definitions === |

187 | ;;; 1. in l0-float.lisp |

188 | |

189 | (eval-when (:execute :compile-toplevel) |

190 | (defconstant two^23 (ash 1 23)) |

191 | ) |

192 | |

193 | (defun sin (x) |

194 | "Return the sine of NUMBER." |

195 | (cond ((complexp x) |

196 | (let* ((r (realpart x)) |

197 | (i (imagpart x))) |

198 | (complex (* (sin r) (cosh i)) |

199 | (* (cos r) (sinh i))))) |

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

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

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

203 | (imagpart (%extended-cis x)) |

204 | (%short-float (imagpart (%extended-cis x))))) |

205 | #| |

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

207 | (%double-float-sin! x (%make-dfloat))) |

208 | |# |

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

210 | (cl:sin x)) |

211 | #| |

212 | (t |

213 | #+32-bit-target |

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

215 | (%single-float-sin! sx (%make-sfloat))) |

216 | #+64-bit-target |

217 | (%single-float-sin (%short-float x))) |

218 | |# |

219 | (t (cl:sin (%short-float x))))) |

220 | |

221 | (defun cos (x) |

222 | "Return the cosine of NUMBER." |

223 | (cond ((complexp x) |

224 | (let* ((r (realpart x)) |

225 | (i (imagpart x))) |

226 | (complex (* (cos r) (cosh i)) |

227 | (- (* (sin r) (sinh i)))))) |

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

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

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

231 | (realpart (%extended-cis x)) |

232 | (%short-float (realpart (%extended-cis x))))) |

233 | #| |

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

235 | (%double-float-cos! x (%make-dfloat))) |

236 | |# |

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

238 | (cl:cos x)) |

239 | #| |

240 | (t |

241 | #+32-bit-target |

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

243 | (%single-float-cos! sx (%make-sfloat))) |

244 | #+64-bit-target |

245 | (%single-float-cos (%short-float x))) |

246 | |# |

247 | (t (cl:cos (%short-float x))))) |

248 | |

249 | (defun tan (x) |

250 | "Return the tangent of NUMBER." |

251 | (cond ((complexp x) |

252 | (let ((r (realpart x)) |

253 | (i (imagpart x))) |

254 | (if (zerop i) |

255 | (complex (tan r) i) |

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

257 | (ty (tanh i)) |

258 | (tx2 (* tx tx)) |

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

260 | (n (if (> (abs i) 20) |

261 | (* 4 (exp (* -2 (abs i)))) |

262 | (let ((c (cosh i))) |

263 | (/ (* c c)))))) |

264 | (complex (/ (* n tx) d) |

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

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

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

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

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

270 | (/ (imagpart c) (realpart c)) |

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

272 | #| |

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

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

275 | |# |

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

277 | (cl:tan x)) |

278 | #| |

279 | (t |

280 | #+32-bit-target |

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

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

283 | #+64-bit-target |

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

285 | ) |

286 | |# |

287 | (t (cl:tan (%short-float x))))) |

288 | |

289 | ;;; Helper function for sin/cos/tan for ratio or large arguments |

290 | ;;; (Will become inaccurate for ridiculously large arguments though) |

291 | ;;; Does not assume that float library returns accurate values for large arguments |

292 | ;;; (many don't) |

293 | |

294 | ;;; hexdecimal representations of pi at various precisions |

295 | (defconstant pi-vector |

296 | #(#x3243F6A8885A308D313198A2E0 |

297 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D008 |

298 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D |

299 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC |

300 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B5470 |

301 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310B |

302 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB |

303 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045 |

304 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70 |

305 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871 |

306 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D9 |

307 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD588215 |

308 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D |

309 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085 |

310 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085F0CA417918B8DB38EF8E79DCB |

311 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085F0CA417918B8DB38EF8E79DCB0603A180E6C9E0E8BB01E8A3E |

312 | #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085F0CA417918B8DB38EF8E79DCB0603A180E6C9E0E8BB01E8A3ED71577C1BD314B2778AF2FDA5 |

313 | )) |

314 | |

315 | (defun %extended-cis (x) |

316 | (let (size pi-size) |

317 | (typecase x |

318 | (integer (setq size (1- (integer-length (abs x))))) |

319 | (ratio (setq size (- (integer-length (abs (numerator x))) |

320 | (integer-length (denominator x))))) |

321 | (float (multiple-value-bind (mantissa exponent sign) |

322 | (integer-decode-float x) |

323 | (setq size (+ (1- (integer-length mantissa)) exponent)) |

324 | (setq x (* sign mantissa (expt 2 exponent)))))) |

325 | (setq pi-size (ceiling (+ size 64) 100)) |

326 | (cond ((< pi-size 1) (setq pi-size 1)) |

327 | ((> pi-size 17) (setq pi-size 17))) |

328 | (let* ((2pi-approx (/ (aref pi-vector (1- pi-size)) |

329 | (ash 1 (1- (* 100 pi-size))))) |

330 | (reduced-x (rem x 2pi-approx)) |

331 | (x0 (float reduced-x 1.0d0)) |

332 | (x1 (- reduced-x (rational x0)))) |

333 | (* (cis x0) (cis (float x1 1.0d0)))))) |

334 | |

335 | ;;; Multiply by i (with additional optional scale factor) |

336 | ;;; Does the "right thing" with minus zeroes (see CLTL2) |

337 | (defun i* (number &optional (scale 1)) |

338 | (complex (* (- scale) (imagpart number)) |

339 | (* scale (realpart number)))) |

340 | |

341 | |

342 | |

343 | ;;; >>> move these values of constants up from below <<< |

344 | (eval-when (:execute :compile-toplevel) |

345 | (defconstant single-float-pi (coerce pi 'single-float)) |

346 | (defconstant double-float-half-pi (/ pi 2)) |

347 | (defconstant single-float-half-pi (coerce (/ pi 2) 'single-float)) |

348 | ) |

349 | |

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

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

352 | (cond (x-p |

353 | (cond ((or (typep x 'double-float) |

354 | (typep y 'double-float)) |

355 | (with-stack-double-floats ((dy y) |

356 | (dx x)) |

357 | #| |

358 | (%df-atan2 dy dx) |

359 | |# |

360 | (cl:atan (%double-float dy) (%double-float dx)))) |

361 | (t |

362 | (when (and (rationalp x) (rationalp y)) |

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

364 | (let ((x1 (abs x)) (y1 (abs y))) |

365 | (cond ((> y1 x1) |

366 | (setf x (/ x y1)) |

367 | (setf y (signum y))) |

368 | ((not (zerop x)) |

369 | (setf y (/ y x1)) |

370 | (setf x (signum x)))))) |

371 | #| |

372 | #+32-bit-target |

373 | (target::with-stack-short-floats ((sy y) |

374 | (sx x)) |

375 | (%sf-atan2! sy sx)) |

376 | #+64-bit-target |

377 | (%sf-atan2 (%short-float y) (%short-float x)) |

378 | |# |

379 | (cl:atan (%short-float y) (%short-float x))))) |

380 | ((typep y 'float) |

381 | (cl:atan y)) |

382 | #| |

383 | ((typep y 'double-float) |

384 | (%double-float-atan! y (%make-dfloat))) |

385 | ((typep y 'single-float) |

386 | #+32-bit-target |

387 | (%single-float-atan! y (%make-sfloat)) |

388 | #+64-bit-target |

389 | (%single-float-atan y)) |

390 | |# |

391 | ((typep y 'rational) |

392 | (cond ((<= (abs y) most-positive-short-float) |

393 | #| |

394 | #+32-bit-target |

395 | (target::with-stack-short-floats ((sy y)) |

396 | (%single-float-atan! sy (%make-sfloat))) |

397 | #+64-bit-target |

398 | (%single-float-atan (%short-float y)) |

399 | |# |

400 | (cl:atan (%short-float y))) |

401 | ((minusp y) |

402 | #.(- single-float-half-pi)) |

403 | (t |

404 | single-float-half-pi))) |

405 | (t |

406 | (let ((r (realpart y)) |

407 | (i (imagpart y))) |

408 | (if (zerop i) |

409 | (complex (atan r) i) |

410 | (i* (%complex-atanh (complex (- i) r)) -1)))) |

411 | )) |

412 | |

413 | |

414 | ;;; useful constants |

415 | (eval-when (:execute :compile-toplevel) |

416 | (defconstant single-float-log2 0.6931472) ; (log 2) |

417 | (defconstant double-float-log2 0.6931471805599453d0) ; (log 2.0d0) |

418 | (defconstant double-float-log2^23 15.942385152878742d0) ; (log (expt 2 23)) |

419 | ) |

420 | |

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

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

423 | (if b-p |

424 | (if (zerop b) |

425 | (if (zerop x) |

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

427 | ;(if (floatp x) (float 0.0d0 x) 0) |

428 | ; ** CORRECT THE CONTAGION for complex args ** |

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

430 | ; do the float/rational contagion before the division |

431 | ; but take care with negative zeroes |

432 | (let ((x1 (realpart x)) |

433 | (b1 (realpart b))) |

434 | (if (and (typep x1 'float) |

435 | (typep b1 'float)) |

436 | (/ (log-e (* x (float 1.0 b1))) |

437 | (log-e (* b (float 1.0 x1)))) |

438 | (let ((r (/ (cond ((typep x 'rational) |

439 | (%rational-log x 1.0d0)) |

440 | ((typep x1 'rational) |

441 | (%rational-complex-log x 1.0d0)) |

442 | (t |

443 | (log-e (* x 1.0d0)))) |

444 | (cond ((typep b 'rational) |

445 | (%rational-log b 1.0d0)) |

446 | ((typep b1 'rational) |

447 | (%rational-complex-log b 1.0d0)) |

448 | (t |

449 | (log-e (* b 1.0d0))))))) |

450 | (cond ((or (typep x1 'double-float) |

451 | (typep b1 'double-float)) |

452 | r) |

453 | ((complexp r) |

454 | (complex (%short-float (realpart r)) |

455 | (%short-float (imagpart r)))) |

456 | (t |

457 | (%short-float r))))))) |

458 | (log-e x))) |

459 | |

460 | ;;; Cases rearranged to put more esoteric cases later |

461 | (defun log-e (x) |

462 | (cond |

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

464 | (if (%double-float-sign x) |

465 | (complex (cl:log (- x)) pi) |

466 | (cl:log x))) |

467 | ((typep x 'short-float) |

468 | (if (%short-float-sign x) |

469 | (complex (cl:log (- x)) single-float-pi) |

470 | (cl:log x))) |

471 | #| |

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

473 | (if (%double-float-sign x) |

474 | (with-stack-double-floats ((dx x)) |

475 | (complex (%double-float-log! (%%double-float-abs! dx dx) (%make-dfloat)) pi)) |

476 | (%double-float-log! x (%make-dfloat)))) |

477 | ((typep x 'short-float) |

478 | #+32-bit-target |

479 | (if (%short-float-sign x) |

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

481 | (complex (%single-float-log! (%%short-float-abs! sx sx) (%make-sfloat)) |

482 | single-float-pi)) |

483 | (%single-float-log! x (%make-sfloat))) |

484 | #+64-bit-target |

485 | (if (%short-float-sign x) |

486 | (complex (%single-float-log (%short-float-abs (%short-float x))) |

487 | single-float-pi) |

488 | (%single-float-log (%short-float x)))) |

489 | |# |

490 | ((typep x 'complex) |

491 | (if (typep (realpart x) 'rational) |

492 | (%rational-complex-log x 1.0s0) |

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

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

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

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

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

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

499 | (abs (imagpart x)))) |

500 | (log-s 0) |

501 | (s 1)) |

502 | (if (typep m 'short-float) |

503 | (let ((expon (- (%short-float-exp m) IEEE-single-float-bias))) |

504 | (cond ((> expon 126) |

505 | (setq log-s double-float-log2^23) |

506 | (setq s #.(ash 1 23))) |

507 | ((< expon -124) |

508 | (setq log-s #.(- double-float-log2^23)) |

509 | (setq s #.(/ 1.0s0 (ash 1 23)))))) |

510 | (let ((expon (- (%double-float-exp m) IEEE-double-float-bias))) |

511 | (cond ((> expon 1022) |

512 | (setq log-s double-float-log2^23) |

513 | (setq s #.(ash 1 23))) |

514 | ((< expon -1020) |

515 | (setq log-s #.(- double-float-log2^23)) |

516 | (setq s #.(/ 1.0d0 (ash 1 23))))))) |

517 | (if (eql s 1) |

518 | (complex (log-abs x) (phase x)) |

519 | (let ((temp (log-abs (/ x s)))) |

520 | (complex (float (+ log-s temp) temp) (phase x))))))) |

521 | (t |

522 | (%rational-log x 1.0s0)))) |

523 | |

524 | ;;; (log1+ x) = (log (1+ x)) |

525 | ;;; but is much more precise near x = 0 |

526 | (defun log1+ (x) |

527 | ;(cond ((typep x 'complex) |

528 | ; (let ((r (realpart x)) |

529 | ; (i (imagpart x))) |

530 | ; (if (and (< (abs r) 0.5) |

531 | ; (< (abs i) 3)) |

532 | ; (let* ((n (+ (* r (+ 2 r)) (* i i))) |

533 | ; (d (1+ (sqrt (1+ n))))) |

534 | ; (complex (log1+ (/ n d)) (atan i (1+ r)))) |

535 | ; (log (1+ x))))) |

536 | ; (t |

537 | (if (and (typep x 'ratio) |

538 | (< -0.5 x 2)) |

539 | (setq x (%short-float x))) |

540 | (let ((y (1+ x))) |

541 | (if (eql y x) |

542 | (log-e y) |

543 | (let ((e (1- y))) |

544 | (if (zerop e) |

545 | (* x 1.0) |

546 | (- (log-e y) (/ (- e x) y))))))) |

547 | ;)) |

548 | |

549 | ;;; helper for complex log |

550 | ;;; uses more careful approach when (abs x) is near 1 |

551 | (defun log-abs (x) |

552 | (let ((a (abs x))) |

553 | (if (< 0.5 a 3) |

554 | (let* ((r (realpart x)) |

555 | (i (imagpart x)) |

556 | (n (if (> (abs r) (abs i)) |

557 | (+ (* (1+ r) (1- r)) (* i i)) |

558 | (+ (* r r) (* (1+ i) (1- i)))))) |

559 | (log1+ (/ n (1+ a)))) |

560 | (log-e a)))) |

561 | |

562 | ;;; helpers for rational log |

563 | (defun %rational-log (x prototype) |

564 | (cond ((minusp x) |

565 | (complex (%rational-log (- x) prototype) |

566 | (if (typep prototype 'short-float) |

567 | single-float-pi |

568 | pi))) |

569 | ((bignump x) |

570 | ;(let* ((base1 3) |

571 | ; (guess (floor (1- (integer-length x)) |

572 | ; (log base1 2))) |

573 | ; (guess1 (* guess (log-e base1)))) |

574 | ; (+ guess1 (log-e (/ x (expt base1 guess))))) |

575 | ; Using base1 = 2 is *much* faster. Was there a reason for 3? |

576 | (let* ((guess (1- (integer-length x))) |

577 | (guess1 (* guess double-float-log2))) |

578 | (float (+ guess1 (log-e (float (/ x (ash 1 guess)) 1.0d0))) prototype))) |

579 | ((and (ratiop x) |

580 | ;; Rational arguments near +1 can be specified with great precision: don't lose it |

581 | (cond ((< 0.5 x 3) |

582 | (log1+ (float (- x 1) prototype))) |

583 | ( |

584 | ;; Escape out small values as well as large |

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

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

587 | ;; avoid loss of precision due to subtracting logs of numerator and denominator |

588 | (let* ((n (%numerator x)) |

589 | (d (%denominator x)) |

590 | (sn (1- (integer-length n))) |

591 | (sd (1- (integer-length d)))) |

592 | (float (+ (* (- sn sd) double-float-log2) |

593 | (- (log1+ (float (1- (/ n (ash 1 sn))) 1.0d0)) |

594 | (log1+ (float (1- (/ d (ash 1 sd))) 1.0d0)))) |

595 | prototype)))))) |

596 | ((typep prototype 'short-float) |

597 | #| |

598 | #+32-bit-target |

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

600 | (%single-float-log! sx (%make-sfloat))) |

601 | #+64-bit-target |

602 | (%single-float-log (%short-float x)) |

603 | |# |

604 | (cl:log (%short-float x))) |

605 | (t |

606 | #| |

607 | (with-stack-double-floats ((dx x)) |

608 | (%double-float-log! dx (%make-dfloat))) |

609 | |# |

610 | (cl:log (%double-float x))))) |

611 | |

612 | (defun %rational-complex-log (x prototype &aux ra ia) |

613 | (let* ((rx (realpart x)) |

614 | (ix (imagpart x)) |

615 | (x (abs rx)) |

616 | (y (abs ix))) |

617 | (if (> y x) |

618 | (let ((r (float (/ rx y) 1.0d0))) |

619 | (setq ra (+ (%rational-log y 1.0d0) |

620 | (/ (log1+ (* r r)) 2))) |

621 | (setq ia (atan (if (minusp ix) -1.0d0 1.0d0) r))) |

622 | (let ((r (float (/ ix x) 1.0d0))) |

623 | (setq ra (+ (%rational-log x 1.0d0) |

624 | (/ (log1+ (* r r)) 2))) |

625 | (setq ia (atan r (if (minusp rx) -1.0d0 1.0d0))))) |

626 | (if (typep prototype 'short-float) |

627 | (complex (%short-float ra) (%short-float ia)) |

628 | (complex ra ia)))) |

629 | |

630 | |

631 | |

632 | ;;; (It may be possible to do something with rational exponents, e.g. so that |

633 | ;;; (expt x 1/2) == (sqrt x) |

634 | ;;; (expt x 3/2) == (expt (sqrt x) 3) ** NOT (sqrt (expt x 3)) !! ** |

635 | ;;; or (* x (sqrt x)) |

636 | ;;; (expt x 1/8) == (sqrt (sqrt (sqrt x))) |

637 | ;;; even, in the case of rational x, returning a rational result if possible.) |

638 | ;;; |

639 | (defun expt (b e) |

640 | "Return BASE raised to the POWER." |

641 | (cond ((zerop e) (1+ (* b e))) |

642 | ((integerp e) |

643 | (if (minusp e) (/ 1 (%integer-power b (- e))) (%integer-power b e))) |

644 | ((zerop b) |

645 | (if (plusp (realpart e)) (* b e) (report-bad-arg e '(number (0) *)))) |

646 | ((and (realp b) (plusp b) (realp e) |

647 | ; escape out very small or very large rationals |

648 | ; - avoid problems converting to float |

649 | (typecase b |

650 | (bignum (<= b most-positive-short-float)) |

651 | (ratio (cond ((< b 0.5) |

652 | (>= b least-positive-normalized-short-float)) |

653 | ((> b 3) |

654 | (<= b most-positive-short-float)))) |

655 | (t t))) |

656 | ;; assumes that the library routines are accurate |

657 | ;; (if not, just excise this whole clause) |

658 | #| |

659 | (if (or (typep b 'double-float) |

660 | (typep e 'double-float)) |

661 | (with-stack-double-floats ((b1 b) |

662 | (e1 e)) |

663 | (%double-float-expt! b1 e1 (%make-dfloat))) |

664 | #+32-bit-target |

665 | (target::with-stack-short-floats ((b1 b) |

666 | (e1 e)) |

667 | (%single-float-expt! b1 e1 (%make-sfloat))) |

668 | #+64-bit-target |

669 | (%single-float-expt (%short-float b) (%short-float e)) |

670 | ) |

671 | |# |

672 | (if (or (typep b 'double-float) |

673 | (typep e 'double-float)) |

674 | (cl:expt (%double-float b) (%double-float e)) |

675 | (cl:expt (%short-float b) (%short-float e)))) |

676 | ((typep b 'rational) |

677 | (let ((r (exp (* e (%rational-log b 1.0d0))))) |

678 | (cond ((typep (realpart e) 'double-float) |

679 | r) |

680 | ((typep r 'complex) |

681 | (complex (%short-float (realpart r)) (%short-float (imagpart r)))) |

682 | (t |

683 | (%short-float r))))) |

684 | ((typep (realpart b) 'rational) |

685 | (let ((r (exp (* e (%rational-complex-log b 1.0d0))))) |

686 | (if (typep (realpart e) 'double-float) |

687 | r |

688 | (complex (%short-float (realpart r)) (%short-float (imagpart r)))))) |

689 | (t |

690 | ;; type upgrade b without losing -0.0 ... |

691 | (let ((r (exp (* e (log-e (* b 1.0d0)))))) |

692 | (cond ((or (typep (realpart b) 'double-float) |

693 | (typep (realpart e) 'double-float)) |

694 | r) |

695 | ((typep r 'complex) |

696 | (complex (%short-float (realpart r)) (%short-float (imagpart r)))) |

697 | (t |

698 | (%short-float r))))))) |

699 | |

700 | |

701 | |

702 | ;;; For complex square root, |

703 | ;;; if (a + i*b)^2 = rx + i * ix |

704 | ;;; then, the basic solution is |

705 | ;;; a = sqrt((rx + sqrt(rx^2 + ix^2))/2) |

706 | ;;; b = ix / (2 * a) |

707 | ;;; but, in practice, there are various cases that must be handled specially |

708 | ;;; Note: for rational complex arguments, we can search for (and return) exact results |

709 | (defun sqrt (x &aux a b) |

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

711 | (cond ((zerop x) x) |

712 | ((complexp x) |

713 | (let ((rx (realpart x)) |

714 | (ix (imagpart x))) |

715 | (cond ((rationalp rx) |

716 | (if (zerop rx) |

717 | (let ((s (sqrt (/ (abs ix) 2)))) |

718 | (complex s (if (minusp ix) (- s) s))) |

719 | (let* ((s (+ (* rx rx) (* ix ix))) |

720 | (d (if (ratiop s) |

721 | (/ (isqrt (%numerator s)) |

722 | (isqrt (%denominator s))) |

723 | (isqrt s)))) |

724 | (unless (eql s (* d d)) |

725 | (setf d (%hypot (%double-float rx) |

726 | (%double-float ix)))) |

727 | (cond ((minusp rx) |

728 | (setq b (sqrt (/ (- d rx) 2))) |

729 | (when (minusp ix) |

730 | (setq b (- b))) |

731 | (setq a (/ ix (+ b b)))) |

732 | (t |

733 | (setq a (sqrt (/ (+ rx d) 2))) |

734 | (setq b (/ ix (+ a a))))) |

735 | (if (rationalp a) |

736 | (complex a b) |

737 | (complex (%short-float a) (%short-float b)))))) |

738 | ((minusp rx) |

739 | (if (zerop ix) |

740 | (complex 0 (float-sign ix (sqrt (- rx)))) |

741 | (let ((shift (cond ((< rx -1) -3) |

742 | ((and (> rx -5.9604645E-8) (< (abs ix) 5.9604645E-8)) 25) |

743 | (t -1)))) |

744 | (setq rx (scale-float rx shift)) |

745 | (let ((s (fsqrt (- (abs (complex rx (scale-float ix shift))) rx)))) |

746 | (setq b (scale-float s (ash (- -1 shift) -1))) |

747 | (when (minusp ix) |

748 | (setq b (- b))) |

749 | (setq a (/ ix (scale-float b 1))) |

750 | (complex a b))))) |

751 | (t |

752 | (if (zerop ix) |

753 | (complex (sqrt rx) ix) |

754 | (let ((shift (cond ((> rx 1) -3) |

755 | ((and (< rx 5.9604645E-8) (< (abs ix) 5.9604645E-8)) 25) |

756 | (t -1)))) |

757 | (setq rx (scale-float rx shift)) |

758 | (let ((s (fsqrt (+ rx (abs (complex rx (scale-float ix shift))))))) |

759 | (setq a (scale-float s (ash (- -1 shift) -1))) |

760 | (setq b (/ ix (scale-float a 1))) |

761 | (complex a b)))))))) |

762 | ((minusp x) (complex 0 (sqrt (- x)))) |

763 | ((floatp x) |

764 | (fsqrt x)) |

765 | ((and (integerp x) (eql x (* (setq a (isqrt x)) a))) a) |

766 | ((and (ratiop x) |

767 | (let ((n (numerator x)) |

768 | d) |

769 | (and (eql n (* (setq a (isqrt n)) a)) |

770 | (eql (setq d (denominator x)) |

771 | (* (setq b (isqrt d)) b))))) |

772 | (/ a b)) |

773 | (t |

774 | ;; use double float (sometimes x can't be represented in single) |

775 | (with-stack-double-floats ((f1 x)) |

776 | (%short-float (fsqrt f1)))))) |

777 | |

778 | |

779 | |

780 | |

781 | (defun asin (x) |

782 | "Return the arc sine of NUMBER." |

783 | (cond ((and (typep x 'double-float) |

784 | (locally (declare (type double-float x)) |

785 | (and (<= -1.0d0 x) |

786 | (<= x 1.0d0)))) |

787 | #| |

788 | (%double-float-asin! x (%make-dfloat)) |

789 | |# |

790 | (cl:asin x)) |

791 | ((and (typep x 'single-float) |

792 | (locally (declare (type single-float x)) |

793 | (and (<= -1.0s0 x) |

794 | (<= x 1.0s0)))) |

795 | #| |

796 | #+32-bit-target |

797 | (%single-float-asin! x (%make-sfloat)) |

798 | #+64-bit-target |

799 | (%single-float-asin x) |

800 | |# |

801 | (cl:asin x)) |

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

803 | (<= (abs x) 1)) |

804 | (if (integerp x) |

805 | (case x |

806 | (0 0.0s0) ; or simply 0 ?? |

807 | (1 single-float-half-pi) |

808 | (-1 #.(- single-float-half-pi))) |

809 | (atan x (sqrt (- 1 (* x x)))))) |

810 | (t |

811 | (%complex-asin/acos x nil)) |

812 | )) |

813 | |

814 | |

815 | (defun acos (x) |

816 | "Return the arc cosine of NUMBER." |

817 | (cond ((and (typep x 'double-float) |

818 | (locally (declare (type double-float x)) |

819 | (and (<= -1.0d0 x) |

820 | (<= x 1.0d0)))) |

821 | #| |

822 | (%double-float-acos! x (%make-dfloat)) |

823 | |# |

824 | (cl:acos x)) |

825 | ((and (typep x 'short-float) |

826 | (locally (declare (type short-float x)) |

827 | (and (<= -1.0s0 x) |

828 | (<= x 1.0s0)))) |

829 | #| |

830 | #+32-bit-target |

831 | (%single-float-acos! x (%make-sfloat)) |

832 | #+64-bit-target |

833 | (%single-float-acos x) |

834 | |# |

835 | (cl:acos x)) |

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

837 | (<= (abs x) 1)) |

838 | (if (integerp x) |

839 | (case x |

840 | (0 single-float-half-pi) |

841 | (1 0.0s0) ; or simply 0 ?? |

842 | (-1 single-float-pi)) |

843 | (atan (sqrt (- 1 (* x x))) x))) |

844 | (t |

845 | (%complex-asin/acos x t)) |

846 | )) |

847 | |

848 | ;;; combined complex asin/acos routine |

849 | ;;; argument acos is true for acos(z); false for asin(z) |

850 | ;;; |

851 | ;;; based on Hull, Fairgrieve & Tang, ACM TMS 23, 3, 299-335 (Sept. 1997) |

852 | (defun %complex-asin/acos (z acos) |

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

854 | (ix (imagpart z)) |

855 | (x (abs rx)) |

856 | (y (abs ix)) |

857 | (m (max x y))) |

858 | (if (> m 1.8014399E+16) |

859 | ;; Large argument escape |

860 | (let ((log-s 0)) |

861 | (if (typep m 'double-float) |

862 | (if (> m #.(/ most-positive-double-float 2)) |

863 | (setq log-s double-float-log2) |

864 | (setq z (* 2 z))) |

865 | (if (> m #.(/ most-positive-short-float 2)) |

866 | (setq log-s single-float-log2) |

867 | (setq z (* 2 z)))) |

868 | (if acos |

869 | (i* (+ log-s (log-e z)) |

870 | (if (minusp ix) +1 -1)) |

871 | (if (minusp ix) |

872 | (i* (+ log-s (log-e (i* z 1))) -1) |

873 | (i* (+ log-s (log-e (i* z -1))) 1)))) |

874 | (let ((qrx rx) |

875 | (qx x) |

876 | x-1 y2 s) |

877 | (cond ((rationalp rx) |

878 | (setq x-1 (float (abs (- x 1)))) |

879 | (setq rx (float rx)) |

880 | (setq x (abs rx)) |

881 | (setq y (float y)) |

882 | (setq y2 (* y y)) |

883 | (setq s (cond ((zerop x-1) |

884 | y) |

885 | ((> y x-1) |

886 | (let ((c (/ x-1 y))) |

887 | (* y (sqrt (1+ (* c c)))))) |

888 | (t |

889 | (let ((c (/ y x-1))) |

890 | (* x-1 (sqrt (1+ (* c c))))))))) |

891 | (t |

892 | (setq x-1 (abs (- x 1))) |

893 | (setq y2 (* y y)) |

894 | (setq s (if (zerop x-1) |

895 | y |

896 | (sqrt (+ (* x-1 x-1) y2)))))) |

897 | (let* ((x+1 (+ x 1)) |

898 | (r (sqrt (+ (* x+1 x+1) y2))) |

899 | (a (/ (+ r s) 2)) |

900 | (b (/ rx a)) |

901 | (ra (if (<= (abs b) 0.6417) |

902 | (if acos (acos b) (asin b)) |

903 | (let* ((r+x+1 (+ r x+1)) |

904 | (s+x-1 (+ s x-1)) |

905 | (a+x (+ a x)) |

906 | (ry (if (<= qx 1) |

907 | (let ((aa (+ (/ y2 r+x+1) s+x-1))) |

908 | (sqrt (/ (* a+x aa) 2))) |

909 | (let ((aa (+ (/ a+x r+x+1) (/ a+x s+x-1)))) |

910 | (* y (sqrt (/ aa 2))))))) |

911 | (if acos (atan ry rx) (atan rx ry))))) |

912 | (ia (if (<= a 1.5) |

913 | (let* ((r+x+1 (+ r x+1)) |

914 | (s+x-1 (+ s x-1)) |

915 | (ll (if (< qx 1) |

916 | (let* ((aa (/ (+ (/ 1 r+x+1) (/ 1 s+x-1)) 2))) |

917 | (+ (* aa y2) (* y (sqrt (* aa (1+ a)))))) |

918 | (let* ((a-1 (/ (+ (/ y2 r+x+1) s+x-1) 2))) |

919 | (+ a-1 (sqrt (* a-1 (1+ a)))))))) |

920 | (log1+ ll)) |

921 | (log (+ a (sqrt (1- (* a a)))))))) |

922 | ;; final fixup of signs |

923 | (if acos |

924 | (if (complexp z) |

925 | (if (typep ix 'float) |

926 | (setq ia (float-sign (- ix) ia)) |

927 | (if (plusp ix) |

928 | (setq ia (- ia)))) |

929 | (if (< qrx -1) |

930 | (setq ia (- ia)))) |

931 | (if (complexp z) |

932 | (if (typep ix 'float) |

933 | (setq ia (float-sign ix ia)) |

934 | (if (minusp ix) |

935 | (setq ia (- ia)))) |

936 | (if (> qrx 1) |

937 | (setq ia (- ia))))) |

938 | (complex ra ia)))))) |

939 | |

940 | ;;; complex atanh |

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

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

943 | (ix (imagpart z)) |

944 | (sign (typecase rx |

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

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

947 | (t (minusp rx)))) |

948 | (x rx) |

949 | (y ix) |

950 | (y1 (abs y)) |

951 | ra ia) |

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

953 | (when sign |

954 | (setf x (- x)) |

955 | (setf y (- y))) |

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

957 | ;; large value escape |

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

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

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

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

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

963 | (setq ia (typecase y |

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

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

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

967 | ((= x 1) |

968 | (setq ra (if (< y1 1e-9) |

969 | (/ (log-e (/ 2 y1)) 2) |

970 | (/ (log1+ (/ 4 (* y y))) 4))) |

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

972 | (t |

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

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

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

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

977 | (cond (sign |

978 | (setq ra (- ra)) |

979 | (when (typep z 'complex) |

980 | (setq ia (- ia)))) |

981 | (t |

982 | (unless (typep z 'complex) |

983 | (setq ia (- ia))))) |

984 | (complex ra ia))) |

985 | |

986 | |

987 | ;;; |

988 | ;;; 2. In l0-numbers.lisp |

989 | |

990 | (defun cis (theta) |

991 | "Return cos(Theta) + i sin(Theta), i.e. exp(i Theta)." |

992 | (cond ((complexp theta) |

993 | (error "Argument to CIS is complex: ~S" theta)) |

994 | ((or (typep theta 'ratio) |

995 | (> (abs theta) #.(ash 1 23))) |

996 | (if (typep theta 'double-float) |

997 | (%extended-cis theta) |

998 | (coerce (%extended-cis theta) '(complex single-float)))) |

999 | (t |

1000 | (complex (cos theta) (sin theta))))) |

1001 | |

1002 | ;;; |

1003 | ;;; 3. In numbers.lisp |

1004 | |

1005 | (defun signum (x) |

1006 | "If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))." |

1007 | (cond ((complexp x) (if (zerop x) |

1008 | x |

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

1010 | (abs (imagpart x))))) |

1011 | (cond ((rationalp m) |

1012 | ;; rescale to avoid intermediate under/overflow |

1013 | (setq x (/ x m))) |

1014 | ((> m #.(ash 1 23)) |

1015 | ;; ensure no overflow for abs |

1016 | (setq x (/ x 2)))) |

1017 | (/ x (abs x))))) |

1018 | ((rationalp x) (if (plusp x) 1 (if (zerop x) 0 -1))) |

1019 | ((zerop x) (float 0.0 x)) |

1020 | (t (float-sign x)))) |

1021 | |

1022 | (defun sinh (x) |

1023 | "Return the hyperbolic sine of NUMBER." |

1024 | (if (complexp x) |

1025 | (let ((r (realpart x)) |

1026 | (i (imagpart x))) |

1027 | (complex (* (sinh r) (cos i)) |

1028 | (* (cosh r) (sin i)))) |

1029 | #| |

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

1031 | (%double-float-sinh! x (%make-dfloat)) |

1032 | #+32-bit-target |

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

1034 | (%single-float-sinh! sx (%make-sfloat))) |

1035 | #+64-bit-target |

1036 | (%single-float-sinh (%short-float x))) |

1037 | |# |

1038 | (cl:sinh x))) |

1039 | |

1040 | (defun cosh (x) |

1041 | "Return the hyperbolic cosine of NUMBER." |

1042 | (if (complexp x) |

1043 | (let ((r (realpart x)) |

1044 | (i (imagpart x))) |

1045 | (complex (* (cosh r) (cos i)) |

1046 | (* (sinh r) (sin i)))) |

1047 | #| |

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

1049 | (%double-float-cosh! x (%make-dfloat)) |

1050 | #+32-bit-target |

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

1052 | (%single-float-cosh! sx (%make-sfloat))) |

1053 | #+64-bit-target |

1054 | (%single-float-cosh (%short-float x))) |

1055 | |# |

1056 | (cl:cosh x))) |

1057 | |

1058 | (defun tanh (x) |

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

1060 | (cond ((complexp x) |

1061 | (let ((r (realpart x)) |

1062 | (i (imagpart x))) |

1063 | (if (zerop r) |

1064 | (complex r (tan i)) |

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

1066 | (ty (tan i)) |

1067 | (ty2 (* ty ty)) |

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

1069 | (n (if (> (abs r) 20) |

1070 | (* 4 (exp (- (* 2 (abs r))))) |

1071 | (let ((c (cosh r))) |

1072 | (/ (* c c)))))) |

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

1074 | (/ (* n ty) d)))))) |

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

1076 | (cl:tanh x)) |

1077 | #| |

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

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

1080 | |# |

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

1082 | (> (abs x) 12)) |

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

1084 | #| |

1085 | (t |

1086 | #+32-bit-target |

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

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

1089 | #+64-bit-target |

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

1091 | |# |

1092 | (t (cl:tanh (%short-float x))))) |

1093 | |

1094 | (defun asinh (x) |

1095 | "Return the hyperbolic arc sine of NUMBER." |

1096 | (cond ((typep x 'double-float) |

1097 | #| |

1098 | (%double-float-asinh! x (%make-dfloat)) |

1099 | |# |

1100 | (cl:asinh x)) |

1101 | ((typep x 'short-float) |

1102 | #| |

1103 | #+32-bit-target |

1104 | (%single-float-asinh! x (%make-sfloat)) |

1105 | #+64-bit-target |

1106 | (%single-float-asinh (%short-float x)) |

1107 | |# |

1108 | (cl:asinh x)) |

1109 | ((typep x 'rational) |

1110 | (if (< (abs x) most-positive-short-float) |

1111 | #| |

1112 | #+32-bit-target |

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

1114 | (%single-float-asinh! sx (%make-sfloat))) |

1115 | #+64-bit-target |

1116 | (%single-float-asinh (%short-float x)) |

1117 | |# |

1118 | (cl:asinh (%short-float x)) |

1119 | (* (signum x) (log-e (* 2 (abs x)))))) |

1120 | (t |

1121 | (i* (%complex-asin/acos (i* x) nil) -1)))) |

1122 | |

1123 | ;;; for complex case, use acos and post-fix the branch cut |

1124 | (defun acosh (x) |

1125 | "Return the hyperbolic arc cosine of NUMBER." |

1126 | (cond ((and (typep x 'double-float) |

1127 | (locally (declare (type double-float x)) |

1128 | (<= 1.0d0 x))) |

1129 | #| |

1130 | (%double-float-acosh! x (%make-dfloat)) |

1131 | |# |

1132 | (cl:acosh x)) |

1133 | ((and (typep x 'short-float) |

1134 | (locally (declare (type short-float x)) |

1135 | (<= 1.0s0 x))) |

1136 | #| |

1137 | #+32-bit-target |

1138 | (%single-float-acosh! x (%make-sfloat)) |

1139 | #+64-bit-target |

1140 | (%single-float-acosh (%short-float x)) |

1141 | |# |

1142 | (cl:acosh x)) |

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

1144 | (<= 1 x)) |

1145 | (cond ((< x 2) |

1146 | (log1+ (+ (- x 1) (sqrt (- (* x x) 1))))) |

1147 | ((<= x most-positive-short-float) |

1148 | #| |

1149 | #+32-bit-target |

1150 | (target::with-stack-short-floats ((x1 x)) |

1151 | (%single-float-acosh! x1 (%make-sfloat))) |

1152 | #+64-bit-target |

1153 | (%single-float-acosh (%short-float x)) |

1154 | |# |

1155 | (cl:acosh (%short-float x))) |

1156 | (t |

1157 | (log-e (* 2 x))))) |

1158 | (t |

1159 | (let ((sign (and (typep x 'complex) |

1160 | (let ((ix (imagpart x))) |

1161 | (typecase ix |

1162 | (double-float (%double-float-sign ix)) |

1163 | (single-float (%short-float-sign ix)) |

1164 | (t (minusp ix))))))) |

1165 | (i* (%complex-asin/acos x t) (if sign -1 1)))))) |

1166 | |

1167 | (defun atanh (x) |

1168 | "Return the hyperbolic arc tangent of NUMBER." |

1169 | (cond ((and (typep x 'double-float) |

1170 | (locally (declare (type double-float x)) |

1171 | (and (<= -1.0d0 x) |

1172 | (<= x 1.0d0)))) |

1173 | #| |

1174 | (%double-float-atanh! x (%make-dfloat)) |

1175 | |# |

1176 | (cl:atanh x)) |

1177 | ((and (typep x 'short-float) |

1178 | (locally (declare (type short-float x)) |

1179 | (and (<= -1.0s0 x) |

1180 | (<= x 1.0s0)))) |

1181 | #| |

1182 | #+32-bit-target |

1183 | (%single-float-atanh! x (%make-sfloat)) |

1184 | #+64-bit-target |

1185 | (%single-float-atanh x) |

1186 | |# |

1187 | (cl:atanh x)) |

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

1189 | (<= (abs x) 1)) |

1190 | (let ((n (numerator x)) |

1191 | (d (denominator x))) |

1192 | (/ (log-e (/ (+ d n) (- d n))) 2))) |

1193 | (t |

1194 | (let ((r (realpart x))) |

1195 | (if (zerop r) |

1196 | (complex r (atan (imagpart x))) |

1197 | (%complex-atanh x)))))) |

1198 | |

1199 | |

1200 | (defun approx= (x y) |

1201 | (labels ((ulp (x) |

1202 | (multiple-value-bind (m e s) (integer-decode-float x) |

1203 | (if (eql m 0) |

1204 | (values x x) |

1205 | (let ((x2 (scale-float (float (1+ m) x) e)) |

1206 | (x1 (scale-float (float (1- m) x) e))) |

1207 | (if (eql s 1) |

1208 | (values x1 x2) |

1209 | (values (- x2) (- x1))))))) |

1210 | (match (x y) |

1211 | (handler-case |

1212 | (or (and (zerop x) (zerop y)) |

1213 | (multiple-value-bind (x1 x2) (ulp x) |

1214 | (multiple-value-bind (y1 y2) (ulp y) |

1215 | (and (> y2 x1) (> x2 y1))))) |

1216 | (error () nil)))) |

1217 | (and (numberp x) |

1218 | (numberp y) |

1219 | (cond ((eql x y) t) |

1220 | ((or (and (typep x 'single-float) (typep y 'single-float)) |

1221 | (and (typep x 'double-float) (typep y 'double-float))) |

1222 | (and (match x y) :approx)) |

1223 | ((or (and (typep x '(complex single-float)) (typep y '(complex single-float))) |

1224 | (and (typep x '(complex double-float)) (typep y '(complex double-float)))) |

1225 | (and (match (realpart x) (realpart y)) |

1226 | (match (imagpart x) (imagpart y)) |

1227 | :approx)))))) |

1228 | |

1229 | (defun test (list-of-forms) |

1230 | (flet ((eval-form (fn form) |

1231 | (handler-case |

1232 | (eval (cons fn (rest form))) |

1233 | (error () :Error)))) |

1234 | (format t "~&; Expression~54TReference Value~108TCL Value~%") |

1235 | (format t "; ----------~54T---------------~108T--------~%") |

1236 | (dolist (form list-of-forms) |

1237 | (let* ((label (princ-to-string form)) |

1238 | (fn-name (symbol-name (car form))) |

1239 | (ref-fn (find-symbol fn-name :math)) |

1240 | (cl-fn (find-symbol fn-name :cl)) |

1241 | (ref-r (eval-form ref-fn form)) |

1242 | (cl-r (eval-form cl-fn form)) |

1243 | (mark (case (approx= ref-r cl-r) |

1244 | ((nil) #\#) |

1245 | (:approx #\~) |

1246 | (t #\Space)))) |

1247 | (format t |

1248 | "~&;~A~A~A~54T~A~108T~A~159T~A~%" |

1249 | mark |

1250 | (if (> (length label) 50) |

1251 | (subseq label 0 50) |

1252 | label) |

1253 | (if (> (length label) 50) #\\ #\Space) |

1254 | ref-r |

1255 | cl-r |

1256 | mark)))) |

1257 | (values)) |

1258 | |

1259 | (test '( |

1260 | (math:expt 1.0 0) |

1261 | (math:expt (complex 2.0 3.0) 0) |

1262 | (math:expt 0 1.2) |

1263 | (math:expt 0 (complex 2.0d0 3.0d0)) |

1264 | |

1265 | (math:log (expt 10 -66)) |

1266 | (math:log (expt 10 66)) |

1267 | (math:log (- (expt 10 -66))) |

1268 | (math:log (- (expt 10 66))) |

1269 | (math:log (complex (expt 10 66) (expt 10 68))) |

1270 | (math:log (complex (expt 10 -66) (expt 10 -68))) |

1271 | (math:log 8.0d0 2) |

1272 | (math:log #c(2 3) 0) |

1273 | (math:log #c(2.5 3.6) 0) |

1274 | (math:log #c(2.5d0 3.6d0) 0) |

1275 | (math:log (complex most-positive-single-float most-positive-single-float)) |

1276 | (math:log (complex least-positive-single-float most-positive-single-float)) |

1277 | (math:log (complex least-positive-single-float least-positive-single-float)) |

1278 | (math:log (complex most-positive-double-float most-positive-double-float)) |

1279 | (math:log (complex least-positive-double-float least-positive-double-float)) |

1280 | (math:log (rational pi)) |

1281 | (math:log (/ (expt 400 66) (expt 27 66))) |

1282 | (math:log (/ (expt 27 66) (expt 400 66))) |

1283 | |

1284 | (math:sqrt (expt 10 47)) |

1285 | (math:sqrt (/ (expt 10 47) 3)) |

1286 | (math:sqrt (complex (expt 10 46) (expt 10 47))) |

1287 | (math:sqrt (complex (expt 10 -46) (expt 10 -47))) |

1288 | (math:sqrt (complex (expt 10 65) (expt 10 67))) |

1289 | (math:sqrt (complex (expt 10 -65) (expt 10 -67))) |

1290 | (math:sqrt (complex (/ (expt 10 65) 3) (expt 10 67))) |

1291 | (math:sqrt (complex most-positive-single-float most-positive-single-float)) |

1292 | (math:sqrt (complex least-positive-single-float most-positive-single-float)) |

1293 | (math:sqrt (complex least-positive-single-float least-positive-single-float)) |

1294 | (math:sqrt (complex most-positive-double-float most-positive-double-float)) |

1295 | (math:sqrt (complex least-positive-double-float least-positive-double-float)) |

1296 | |

1297 | (math:expt (expt 10 47) 1/2) |

1298 | (math:expt (/ (expt 10 47) 3) 1/2) |

1299 | (math:expt (complex (expt 10 46) (expt 10 47)) 1/2) |

1300 | (math:expt (complex (expt 10 -46) (expt 10 -47)) 1/2) |

1301 | (math:expt (complex (expt 10 65) (expt 10 67)) 1/2) |

1302 | (math:expt (complex (expt 10 -65) (expt 10 -67)) 1/2) |

1303 | (math:expt (complex (/ (expt 10 65) 3) (expt 10 67)) 1/2) |

1304 | (math:expt (expt 10 600) 0.5d0) |

1305 | |

1306 | (math:expt (expt 10 47) 0.5) |

1307 | (math:expt (/ (expt 10 47) 3) 0.5) |

1308 | (math:expt (complex (expt 10 46) (expt 10 47)) 0.5) |

1309 | (math:expt (complex (expt 10 -46) (expt 10 -47)) 0.5) |

1310 | (math:expt (complex (expt 10 65) (expt 10 67)) 0.5) |

1311 | (math:expt (complex (expt 10 -65) (expt 10 -67)) 0.5) |

1312 | (math:expt (complex (/ (expt 10 65) 3) (expt 10 67)) 0.5) |

1313 | |

1314 | (math:expt (complex -2.0 -0.0) (complex 0.5d0 0.0d0)) |

1315 | |

1316 | (math:atan (expt 10 46) (expt 10 47)) |

1317 | (math:atan (expt 10 -46) (expt 10 -47)) |

1318 | (math:atan most-positive-single-float most-positive-single-float) |

1319 | (math:atan least-positive-single-float least-positive-single-float) |

1320 | |

1321 | (math:sin (expt 3 20)) |

1322 | (math:cos (expt 3 20)) |

1323 | (math:tan (expt 3 20)) |

1324 | (math:sin (expt 3 40)) |

1325 | (math:cos (expt 3 40)) |

1326 | (math:tan (expt 3 40)) |

1327 | (math:sin (expt 10 47)) |

1328 | (math:cos (expt 10 47)) |

1329 | (math:tan (expt 10 47)) |

1330 | (math:tan (complex 3 4)) |

1331 | (math:tan (complex 1.2 40)) |

1332 | |

1333 | (math:sinh (complex 1e-5 1e-6)) |

1334 | (math:cosh (complex 1e-5 1e-6)) |

1335 | (math:tanh (complex 1e6 10.0)) |

1336 | |

1337 | (math:tanh (complex 3 4)) |

1338 | (math:tanh (complex 40 1.2)) |

1339 | |

1340 | (math:asin (complex 1.0 0.0)) |

1341 | (math:asin (complex 0.1 1e-23)) |

1342 | (math:asin (complex 0.99 1e-23)) |

1343 | (math:asin (complex 1.01 1e-23)) |

1344 | (math:asin (complex 10.0 1e-23)) |

1345 | (math:asin (complex 1e-23 1.0)) |

1346 | (math:asin (complex 1e-24 1e16)) |

1347 | (math:asin (+ 1 (expt 10 -10))) |

1348 | (math:asin (- 1 (expt 10 -10))) |

1349 | (math:asin (complex 1 (expt 10 -10))) |

1350 | (math:asin (complex 1 (- (expt 10 -10)))) |

1351 | (math:asin -10000) |

1352 | (math:asin (complex most-positive-single-float most-positive-single-float)) |

1353 | |

1354 | (math:acos (complex -1.0 0.0)) |

1355 | (math:acos (complex 0.1 1e-23)) |

1356 | (math:acos (complex 0.99 1e-23)) |

1357 | (math:acos (complex 1.01 1e-23)) |

1358 | (math:acos (complex 10.0 1e-23)) |

1359 | (math:acos (complex 1e-23 1.0)) |

1360 | (math:acos (complex 1e-24 1e16)) |

1361 | (math:acos (+ 1 (expt 10 -10))) |

1362 | (math:acos (- 1 (expt 10 -10))) |

1363 | (math:acos (complex 1 (expt 10 -10))) |

1364 | (math:acos (complex 1 (- (expt 10 -10)))) |

1365 | (math:acos -10000) |

1366 | (math:acos (complex most-positive-single-float most-positive-single-float)) |

1367 | |

1368 | (math:acosh (+ 1 (expt 10 -10))) |

1369 | (math:acosh (- 1 (expt 10 -10))) |

1370 | (math:log (+ 1 (expt 10 -10))) |

1371 | (math:log (- 1 (expt 10 -10))) |

1372 | (math:log 2 (+ 1 (expt 10 -10))) |

1373 | (math:log 2.0 (+ 1 (expt 10 -10))) |

1374 | (math:log 2.0d0 (+ 1 (expt 10 -10))) |

1375 | |

1376 | (math:atan (complex 0 (+ 1 (expt 10 -15)))) |

1377 | (math:atan (complex 0 (- 1 (expt 10 -15)))) |

1378 | (math:atan (complex (expt 10 -10) 1)) |

1379 | (math:atan (complex (- (expt 10 -10)) 1)) |

1380 | (math:atanh (+ 1 (expt 10 -10))) |

1381 | (math:atanh (- 1 (expt 10 -10))) |

1382 | |

1383 | (math:asinh (expt 10 46)) |

1384 | (math:acosh (expt 10 46)) |

1385 | (math:atanh (expt 10 46)) |

1386 | (math:asinh (complex 1e-5 1e-6)) |

1387 | (math:acosh (complex 1e-5 1e-6)) |

1388 | (math:atanh (complex 1e-5 1e-6)) |

1389 | (math:asinh (complex 1e-15 1e-6)) |

1390 | (math:acosh (complex 1e-15 1e-6)) |

1391 | (math:atanh (complex 1e-15 1e-6)) |

1392 | (math:asinh (complex 1e-3 -1e6)) |

1393 | (math:acosh (complex 1e-3 -1e6)) |

1394 | (math:atanh (complex 1e-3 -1e6)) |

1395 | (math:asinh (complex 1d-30 1d-12)) |

1396 | (math:acosh (complex 1d-30 1d-12)) |

1397 | (math:atanh (complex 1d-30 1d-12)) |

1398 | (math:asinh (complex 1d-6 -1d12)) |

1399 | (math:acosh (complex 1d-6 -1d12)) |

1400 | (math:atanh (complex 1d-6 -1d12)) |

1401 | (math:asinh (complex 1e-5 1e-16)) |

1402 | (math:acosh (complex 1e-5 1e-16)) |

1403 | (math:atanh (complex 1e-5 1e-16)) |

1404 | (math:atanh (complex 1.0 1e-32)) |

1405 | )) |

1406 |