1;;; a first stab at a numerics timing test 2(set! (*s7* 'heap-size) (* 6 1024000)) 3 4(define dolph-1 5 (let ((+documentation+ "(dolph-1 n gamma) produces a Dolph-Chebyshev FFT data window of 'n' points using 'gamma' as the window parameter.")) 6 (lambda (N gamma) 7 (let ((vals (make-vector N))) 8 (let ((alpha (cosh (/ (acosh (expt 10.0 gamma)) N)))) 9 (do ((den (/ 1.0 (cosh (* N (acosh alpha))))) 10 (freq (/ pi N)) 11 (mult -1 (- mult)) 12 (i 0 (+ i 1)) 13 (phase (* -0.5 pi))) 14 ((= i N)) 15 (set! (vals i) (* mult den (cos (* N (acos (* alpha (cos phase))))))) 16 (set! phase (+ phase freq)))) 17 ;; now take the DFT 18 (let ((pk 0.0) 19 (w (make-vector N)) 20 (c (* 2.0 0+1.0i pi))) 21 (do ((i 0 (+ i 1)) 22 (sum 0.0 0.0)) 23 ((= i N)) 24 (do ((k 0 (+ k 1))) 25 ((= k N)) 26 (set! sum (+ sum (* (vals k) (exp (/ (* c k i) N)))))) 27 (set! pk (max pk (vector-set! w i (magnitude sum))))) 28 ;; scale to 1.0 (it's usually pretty close already, that is pk is close to 1.0) 29 (do ((i 0 (+ i 1))) 30 ((= i N)) 31 (vector-set! w i (/ (vector-ref w i) pk))) 32 w))))) 33 34(dolph-1 (expt 2 10) 0.5) 35 36 37(define src-duration 38 (let ((+documentation+ "(src-duration envelope) returns the new duration of a sound after using 'envelope' for time-varying sampling-rate conversion")) 39 (lambda (e) 40 (let ((len (- (length e) 2))) 41 (do ((all-x (- (e len) (e 0))) ; last x - first x 42 (dur 0.0) 43 (i 0 (+ i 2))) 44 ((>= i len) dur) 45 (let ((area (let ((x0 (e i)) 46 (x1 (e (+ i 2))) 47 (y0 (e (+ i 1))) ; 1/x x points 48 (y1 (e (+ i 3)))) 49 (if (< (abs (real-part (- y0 y1))) .0001) 50 (/ (- x1 x0) (* y0 all-x)) 51 (/ (* (log (/ y1 y0)) 52 (- x1 x0)) 53 (* (- y1 y0) all-x)))))) 54 (set! dur (+ dur (abs area))))))))) 55 56 57(define (src-test) 58 (src-duration (do ((env (make-float-vector 7000)) 59 (i 0 (+ i 2))) 60 ((= i 7000) env) 61 (set! (env i) i) 62 (set! (env (+ i 1)) (+ .1 (random 1.0)))))) 63(src-test) 64 65 66(define invert-matrix 67 (let ((+documentation+ "(invert-matrix matrix b (zero 1.0e-7)) inverts 'matrix'")) 68 (lambda* (matrix b (zero 1.0e-7)) 69 ;; translated from Numerical Recipes (gaussj) 70 (call-with-exit 71 (lambda (return) 72 (let ((n (car (vector-dimensions matrix)))) 73 (let ((cols (make-int-vector n 0)) 74 (rows (make-int-vector n 0)) 75 (pivots (make-int-vector n 0))) 76 (do ((i 0 (+ i 1)) 77 (col 0 0) 78 (row 0 0)) 79 ((= i n)) 80 (do ((biggest 0.0) 81 (j 0 (+ j 1))) 82 ((= j n) 83 (if (< biggest zero) 84 (return #f))) ; this can be fooled (floats...): (invert-matrix (subvector (float-vector 1 2 3 3 2 1 4 5 6) 0 9 (list 3 3))) 85 (if (not (= (pivots j) 1)) 86 (do ((k 0 (+ k 1))) 87 ((= k n)) 88 (if (= (pivots k) 0) 89 (let ((val (abs (matrix j k)))) 90 (when (> val biggest) 91 (set! col k) 92 (set! row j) 93 (set! biggest val))) 94 (if (> (pivots k) 1) 95 (return #f)))))) 96 (set! (pivots col) (+ (pivots col) 1)) 97 (if (not (= row col)) 98 (let ((temp (if (sequence? b) (b row) 0.0))) 99 (when (sequence? b) 100 (set! (b row) (b col)) 101 (set! (b col) temp)) 102 (do ((k 0 (+ k 1))) 103 ((= k n)) 104 (set! temp (matrix row k)) 105 (set! (matrix row k) (matrix col k)) 106 (set! (matrix col k) temp)))) 107 (set! (cols i) col) 108 (set! (rows i) row) 109 ;; round-off troubles here 110 (if (< (abs (matrix col col)) zero) 111 (return #f)) 112 (let ((inverse-pivot (/ 1.0 (matrix col col)))) 113 (set! (matrix col col) 1.0) 114 (do ((k 0 (+ k 1))) 115 ((= k n)) 116 (set! (matrix col k) (* inverse-pivot (matrix col k)))) 117 (if b (set! (b col) (* inverse-pivot (b col))))) 118 (do ((k 0 (+ k 1))) 119 ((= k n)) 120 (if (not (= k col)) 121 (let ((scl (matrix k col))) 122 (set! (matrix k col) 0.0) 123 (do ((m 0 (+ 1 m))) 124 ((= m n)) 125 (set! (matrix k m) (- (matrix k m) (* scl (matrix col m))))) 126 (if b (set! (b k) (- (b k) (* scl (b col))))))))) 127 (do ((i (- n 1) (- i 1))) 128 ((< i 0)) 129 (if (not (= (rows i) (cols i))) 130 (do ((k 0 (+ k 1))) 131 ((= k n)) 132 (let ((temp (matrix k (rows i)))) 133 (set! (matrix k (rows i)) (matrix k (cols i))) 134 (set! (matrix k (cols i)) temp))))) 135 (list matrix b)))))))) 136 137(define (matrix-solve A b) 138 (cond ((invert-matrix A b) => cadr) (else #f))) 139 140(define (invert-test) 141 (matrix-solve (do ((A (make-float-vector '(100 100))) 142 (i 0 (+ i 1))) 143 ((= i 100) A) 144 (do ((j 0 (+ j 1))) 145 ((= j 100)) 146 (set! (A i j) (random 1.0)))) 147 (do ((b (make-float-vector 100)) 148 (i 0 (+ i 1))) 149 ((= i 100) b) 150 (set! (b i) (random 1.0))))) 151(invert-test) 152 153 154(define* (cfft data n (dir 1)) ; complex data 155 (unless n (set! n (length data))) 156 (do ((i 0 (+ i 1)) 157 (j 0)) 158 ((= i n)) 159 (if (> j i) 160 (let ((temp (data j))) 161 (set! (data j) (data i)) 162 (set! (data i) temp))) 163 (do ((m (/ n 2) (/ m 2))) 164 ((not (<= 2 m j)) 165 (set! j (+ j m))) 166 (set! j (- j m)))) 167 (do ((ipow (floor (log n 2))) 168 (prev 1) 169 (lg 0 (+ lg 1)) 170 (mmax 2 (* mmax 2)) 171 (pow (/ n 2) (/ pow 2)) 172 (theta (complex 0.0 (* pi dir)) (* theta 0.5))) 173 ((= lg ipow)) 174 (do ((wpc (exp theta)) 175 (wc 1.0) 176 (ii 0 (+ ii 1))) 177 ((= ii prev) 178 (set! prev mmax)) 179 (do ((jj 0 (+ jj 1)) 180 (i ii (+ i mmax)) 181 (j (+ ii prev) (+ j mmax)) 182 (tc 0.0)) 183 ((>= jj pow) 184 (set! wc (* wc wpc))) 185 (set! tc (* wc (data j))) 186 (set! (data j) (- (data i) tc)) 187 (set! (data i) (+ (data i) tc))))) 188 data) 189 190(define (cfft-test) 191 (let* ((size (expt 2 16)) 192 (data (make-vector size))) 193 (do ((i 0 (+ i 1))) 194 ((= i size)) 195 (set! (data i) (complex (- 1.0 (random 2.0)) (- 1.0 (random 2.0))))) 196 (cfft data size))) 197 198(cfft-test) 199 200 201;; these Bessel functions are from Numerical Recipes 202(define (bes-j0-1 x) ;returns J0(x) for any real x 203 (if (< (abs x) 8.0) 204 (let ((y (* x x))) 205 (let ((ans1 (+ 57568490574.0000 (* y (- (* y (+ 651619640.7 (* y (- (* y (+ 77392.33017 (* y -184.9052456))) 11214424.18)))) 13362590354.0)))) 206 (ans2 (+ 57568490411.0 207 (* y (+ 1029532985.0 208 (* y (+ 9494680.718 209 (* y (+ 59272.64853 210 (* y (+ 267.8532712 y))))))))))) 211 (/ ans1 ans2))) 212 (let* ((ax (abs x)) 213 (z (/ 8.0 ax)) 214 (y (* z z))) 215 (let ((xx (- ax 0.785398164)) 216 (ans1 (+ 1.0 (* y (- (* y (+ 2.734510407e-05 (* y (- (* y 2.093887211e-07) 2.073370639e-06)))) 0.001098628627)))) 217 (ans2 (- (* y (+ 0.0001 (* y (- (* y (+ 7.621095160999999e-07 (* y -9.34945152e-08))) 6.911147651000001e-06)))) 0.0156))) 218 (* (sqrt (/ 0.636619772 ax)) 219 (- (* ans1 (cos xx)) 220 (* z (sin xx) ans2))))))) 221 222(define bes-j1-1 223 (let ((signum (lambda (x) (if (= x 0.0) 0 (if (< x 0.0) -1 1))))) 224 (lambda (x) ;returns J1(x) for any real x 225 (if (< (abs x) 8.0) 226 (let ((y (* x x))) 227 (let ((ans1 (* x (+ 72362614232.0000 (* y (- (* y (+ 242396853.1 (* y (- (* y (+ 15704.4826 (* y -30.16036606))) 2972611.439)))) 7895059235.0))))) 228 (ans2 (+ 144725228442.0 229 (* y (+ 2300535178.0 230 (* y (+ 18583304.74 231 (* y (+ 99447.43394 232 (* y (+ 376.9991397 y))))))))))) 233 (/ ans1 ans2))) 234 (let* ((ax (abs x)) 235 (z (/ 8.0 ax)) 236 (y (* z z))) 237 (let ((xx (- ax 2.356194491)) 238 (ans1 (+ 1.0 (* y (+ 0.00183105 (* y (- (* y (+ 2.457520174e-06 (* y -2.40337019e-07))) 3.516396496e-05)))))) 239 (ans2 (+ 0.0469 (* y (- (* y (+ 8.449199096000001e-06 (* y (- (* y 1.05787412e-07) 8.8228987e-07)))) 0.0002002690873))))) 240 (* (signum x) 241 (sqrt (/ 0.636619772 ax)) 242 (- (* ans1 (cos xx)) 243 (* z (sin xx) ans2))))))))) 244 245(define (bes-jn nn x) 246 (let ((besn (let ((n (abs nn))) 247 (cond ((= n 0) (bes-j0-1 x)) 248 ((= n 1) (bes-j1-1 x)) 249 ((= x 0.0) 0.0) 250 (else 251 (let ((iacc 40) 252 (ans 0.0000) 253 (bigno 1.0e10) 254 (bigni 1.0e-10)) 255 (if (> (abs x) n) 256 (do ((tox (/ 2.0 (abs x))) 257 (bjm (bes-j0-1 (abs x))) 258 (bj (bes-j1-1 (abs x))) 259 (j 1 (+ j 1)) 260 (bjp 0.0)) 261 ((= j n) (set! ans bj)) 262 (set! bjp (- (* j tox bj) bjm)) 263 (set! bjm bj) 264 (set! bj bjp)) 265 (let ((tox (/ 2.0 (abs x))) 266 (m (* 2 (floor (/ (+ n (sqrt (* iacc n))) 2)))) 267 (jsum 0) 268 (bjm 0.0000) 269 (sum 0.0000) 270 (bjp 0.0000) 271 (bj 1.0000)) 272 (do ((j m (- j 1))) 273 ((= j 0)) 274 (set! bjm (- (* j tox bj) bjp)) 275 (set! bjp bj) 276 (set! bj bjm) 277 (when (> (abs bj) bigno) 278 (set! bj (* bj bigni)) 279 (set! bjp (* bjp bigni)) 280 (set! ans (* ans bigni)) 281 (set! sum (* sum bigni))) 282 (if (not (= 0 jsum)) 283 (set! sum (+ sum bj))) 284 (set! jsum (- 1 jsum)) 285 (if (= j n) (set! ans bjp))) 286 (set! ans (/ ans (- (* 2.0 sum) bj))))) 287 (if (and (< x 0.0) (odd? n)) 288 (- ans) 289 ans))))))) 290 (if (and (< nn 0) 291 (odd? nn)) 292 (- besn) 293 besn))) 294 295(define (bes-i0 x) 296 (if (< (abs x) 3.75) 297 (let ((y (expt (/ x 3.75) 2))) 298 (+ 1.0 299 (* y (+ 3.5156229 300 (* y (+ 3.0899424 301 (* y (+ 1.2067492 302 (* y (+ 0.2659732 303 (* y (+ 0.360768e-1 304 (* y 0.45813e-2))))))))))))) 305 (let* ((ax (abs x)) 306 (y (/ 3.75 ax))) 307 (* (/ (exp ax) (sqrt ax)) 308 (+ 0.3989 (* y (+ 0.0133 309 (* y (+ 0.0023 310 (* y (- (* y (+ 0.0092 311 (* y (- (* y (+ 0.02635537 312 (* y (- (* y 0.00392377) 0.01647633)))) 313 0.02057706)))) 314 0.0016))))))))))) 315 316(define (bes-i1 x) ;I1(x) 317 (if (< (abs x) 3.75) 318 (let ((y (expt (/ x 3.75) 2))) 319 (* x (+ 0.5 320 (* y (+ 0.87890594 321 (* y (+ 0.51498869 322 (* y (+ 0.15084934 323 (* y (+ 0.2658733e-1 324 (* y (+ 0.301532e-2 325 (* y 0.32411e-3)))))))))))))) 326 (let ((ax (abs x))) 327 (let ((ans2 (let* ((y (/ 3.75 ax)) 328 (ans1 (+ 0.02282967 (* y (- (* y (+ 0.01787654 (* y -0.00420059))) 0.02895312))))) 329 (+ 0.39894228 (* y (- (* y (- (* y (+ 0.00163801 (* y (- (* y ans1) 0.01031555)))) 0.00362018)) 0.03988024))))) 330 (sign (if (< x 0.0) -1.0 1.0))) 331 (/ (* (exp ax) ans2 sign) (sqrt ax)))))) 332 333(define (bes-in n x) 334 (cond ((= n 0) (bes-i0 x)) 335 ((= n 1) (bes-i1 x)) 336 ((= x 0.0) 0.0) 337 (else 338 (let ((bigno 1.0e10) 339 (bigni 1.0e-10) 340 (ans 0.0000) 341 (tox (/ 2.0 (abs x))) 342 (bip 0.0000) 343 (bi 1.0000) 344 (m (* 2 (+ n (truncate (sqrt (* 40 n)))))) ; iacc=40 345 (bim 0.0000)) 346 (do ((j m (- j 1))) 347 ((= j 0)) 348 (set! bim (+ bip (* j tox bi))) 349 (set! bip bi) 350 (set! bi bim) 351 (when (> (abs bi) bigno) 352 (set! ans (* ans bigni)) 353 (set! bi (* bi bigni)) 354 (set! bip (* bip bigni))) 355 (if (= j n) (set! ans bip))) 356 (if (and (< x 0.0) (odd? n)) 357 (set! ans (- ans))) 358 (* ans (/ (bes-i0 x) bi)))))) 359 360 361(define (fm-complex-component freq-we-want wc wm a b interp using-sine) 362 (let ((sum 0.0) 363 (mxa (ceiling (* 7 a))) 364 (mxb (ceiling (* 7 b)))) 365 (do ((k (- mxa) (+ k 1))) 366 ((>= k mxa)) 367 (do ((j (- mxb) (+ j 1))) 368 ((>= j mxb)) 369 (if (< (abs (- freq-we-want wc (* k wm) (* j wm))) 0.1) 370 (let ((curJI (* (bes-jn k a) 371 (bes-in (abs j) b) 372 (expt 0.0+1.0i j)))) 373 (set! sum (+ sum curJI)))))) 374 (list sum 375 (+ (* (- 1.0 interp) (real-part sum)) 376 (* interp (imag-part sum)))))) 377 378(define (fm-cascade-component freq-we-want wc wm1 a wm2 b) 379 (let ((sum 0.0) 380 (mxa (ceiling (* 7 a))) 381 (mxb (ceiling (* 7 b)))) 382 (do ((k (- mxa) (+ k 1))) 383 ((>= k mxa)) 384 (do ((j (- mxb) (+ j 1))) 385 ((>= j mxb)) 386 (if (< (abs (- freq-we-want wc (* k wm1) (* j wm2))) 0.1) 387 (let ((curJJ (* (bes-jn k a) 388 (bes-jn j (* k b))))) 389 (set! sum (+ sum curJJ)))))) 390 sum)) 391 392(define (test-fm-components) 393 (do ((i 0.0 (+ i .1))) 394 ((>= i 3.0)) 395 (do ((k 1.0 (+ k 1.0))) 396 ((>= k 10.0)) 397 (fm-complex-component 1200 1000 100 i k 0.0 #f) 398 (fm-cascade-component 1200 1000 100 i 50 k)))) 399 400(test-fm-components) 401 402 403(define (pisum) ; from Julia microbenchmarks 404 (let ((sum 0.0)) 405 (do ((j 0 (+ j 1))) 406 ((= j 500) 407 sum) 408 (set! sum 0.0) 409 (do ((k 1 (+ k 1))) 410 ((> k 10000)) 411 (set! sum (+ sum (/ (* k k)))))))) 412 413(define (quicksort a lo hi) ; from Julia microbenchmarks (slightly optimized for s7) 414 (do ((i lo) 415 (j hi hi)) 416 ((>= i hi)) 417 (do ((t (a 0)) 418 (pivot (a (floor (/ (+ lo hi) 2))))) 419 ((> i j)) 420 (set! i (do ((k i (+ k 1))) 421 ((>= (a k) pivot) k))) 422 (set! j (do ((k j (- k 1))) 423 ((<= (a k) pivot) k))) 424 (when (<= i j) 425 (set! t (a i)) 426 (set! (a i) (a j)) 427 (set! (a j) t) 428 (set! i (+ i 1)) 429 (set! j (- j 1)))) 430 (if (< lo j) 431 (quicksort a lo j)) 432 (set! lo i))) 433 434(define (jtests) 435 (let-temporarily (((*s7* 'equivalent-float-epsilon) 1e-12)) 436 (unless (equivalent? (pisum) 1.6448340718480652) 437 (format *stderr* "pisum: ~S~%" (pisum)))) 438 439 (do ((n 0 (+ n 1)) ; Julia original only runs it once! 440 (a (make-float-vector 5000) (make-float-vector 5000))) 441 ((= n 10)) 442 (do ((i 0 (+ i 1))) 443 ((= i 5000)) 444 (set! (a i) (random 5000.0))) 445 (quicksort a 0 4999) 446 (when (zero? n) 447 ;; make sure it worked... 448 (do ((i 1 (+ i 1))) 449 ((= i 5000)) 450 (if (< (a i) (a (- i 1))) 451 (display ".")))))) 452 453(jtests) 454 455 456(define (mean v) 457 (let ((len (length v)) 458 (sum 0)) 459 (do ((i 0 (+ i 1))) 460 ((= i len) 461 (/ sum len)) 462 (set! sum (+ sum (v i)))))) 463 464(define (medium v) 465 (let ((len (length v))) 466 (quicksort v 0 (- len 1)) 467 (if (odd? len) 468 (v (ceiling (/ len 2))) 469 (/ (+ (v (/ len 2)) (v (+ (/ len 2) 1))) 2)))) 470 471(define (mtests) 472 (let ((v (make-float-vector 5000))) 473 (do ((n 0 (+ n 1))) 474 ((= n 10)) 475 (do ((i 0 (+ i 1))) 476 ((= i 5000)) 477 (set! (v i) (- (random 2.0) 1.0))) 478 (mean v) 479 (medium v)))) 480 481(mtests) 482 483 484(define (gammln xx) ;Ln(gamma(xx)), xx>0 485 (let* ((stp 2.5066282746310005) 486 (x xx) 487 (tmp (+ x 5.5)) 488 (tmp1 (- tmp (* (+ x 0.5) (log tmp)))) 489 (ser (+ 1.000000000190015 490 (/ 76.18009172947146 (+ x 1.0)) 491 (/ -86.50532032941677 (+ x 2.0)) 492 (/ 24.01409824083091 (+ x 3.0)) 493 (/ -1.231739572450155 (+ x 4)) 494 (/ 0.1208650973866179e-2 (+ x 5.0)) 495 (/ -0.5395239384953e-5 (+ x 6.0))))) 496 (- (log (/ (* stp ser) x)) tmp1))) 497 498(define (gser a x) ;P(a,x) evaluated as series, also Ln(gamma) 499 (if (< x 0.0) 500 (format #f "~F is less than 0" x) 501 (if (= x 0.0) 502 0.0 503 (let* ((gln (gammln a)) 504 (itmax 100) 505 (eps 3.0e-7) 506 (ap a) 507 (sum (/ 1.0 a)) 508 (del sum)) 509 (do ((n 1 (+ n 1))) 510 ((or (> n itmax) 511 (< (abs del) (* (abs sum) eps))) 512 (* sum (exp (- (* a (log x)) x gln)))) 513 (set! ap (+ ap 1)) 514 (set! del (* del (/ x ap))) 515 (set! sum (+ sum del))))))) 516 517(define (gcf a x) ;Q(a,x) evaluated as continued fraction 518 (let ((itmax 100) 519 (eps 3.0e-7) 520 (gln (gammln a)) 521 (gold 0.0) ;previous value of g, tested for convergence 522 (a0 1.0) 523 (a1 x) 524 (b0 0.0) 525 (b1 1.0) ;setting up continued fraction 526 (fac 1.0) 527 (ana 0.0) (g 0.0) (anf 0.0)) 528 (call-with-exit 529 (lambda (return) 530 (do ((n 1 (+ n 1))) 531 ((> n itmax) 532 (* g (exp (- (* a (log x)) x gln)))) 533 (set! ana (- n a)) 534 (set! a0 (* fac (+ a1 (* a0 ana)))) 535 (set! b0 (* fac (+ b1 (* b0 ana)))) 536 (set! anf (* fac n)) 537 (set! a1 (+ (* x a0) (* anf a1))) 538 (set! b1 (+ (* x b0) (* anf b1))) 539 (unless (= 0.0 a1) ;renormalize? 540 (set! fac (/ 1.0 a1)) 541 (set! g (* b1 fac)) 542 (if (< (abs (/ (- g gold) g)) eps) 543 (return (* g (exp (- (* a (log x)) x gln))))))))))) 544 545(define (gammp a x) ; incomplete gamma function P(a,x) 546 (if (or (<= a 0.0) 547 (< x 0.0)) 548 (format #f "Invalid argument to gammp: ~F" (if (<= 0.0 a) a x)) 549 (if (< x (+ a 1.0)) 550 (gser a x) ;use series 551 (- 1.0 (gcf a x))))) ;use continued fraction 552 553(define (erf x) ; error function erf(x) 554 (if (< x 0.0) 555 (- (gammp 0.5 (* x x))) 556 (gammp 0.5 (* x x)))) 557 558(define (gammq a x) ;incomplete gamma function Q(a,x) = 1 - P(a,x) 559 (- 1.0 (gammp a x))) 560 561(define (erfc x) ;complementary error function erfc(x) 562 (if (< x 0.0) 563 (+ 1.0 (gammp 0.5 (* x x))) 564 (gammq 0.5 (* x x)))) 565 566(define (test-erf) 567 (let-temporarily (((*s7* 'equivalent-float-epsilon) 1e-12)) 568 (unless (equivalent? (erf 0) 0.0) (format *stderr* "erf 0: ~S~%" (erf 0))) 569 (unless (equivalent? (erf 1) 0.8427007900291826) (format *stderr* "erf 1: ~S~%" (erf 1))) 570 (unless (equivalent? (erfc 1) 0.15729920997081737) (format *stderr* "erfc 1: ~S~%" (erfc 1))) 571 (unless (equivalent? (erf 0.5) 0.5204998760683841) (format *stderr* "erf 0.5: ~S~%" (erf 0.5))) 572 (unless (equivalent? (erf 2) 0.9953222650189529) (format *stderr* "erf 2: ~S~%" (erf 2))) 573 (unless (equivalent? (erf 0.35) 0.3793820529938486) (format *stderr* "erf 0.35: ~S~%" (erf 0.35))) 574 575 (do ((i 0 (+ i 1)) 576 (x 0.0 (+ x .001))) 577 ((= i 3000)) 578 (unless (equivalent? (+ (erf x) (erfc x)) 1.0) 579 (format *stderr* "erf: ~S trouble\n" x))))) 580 581(test-erf) 582 583 584(define show-digits-of-pi-starting-at-digit 585 ;; piqpr8.c 586 ;; This program implements the BBP algorithm to generate a few hexadecimal 587 ;; digits beginning immediately after a given position id, or in other words 588 ;; beginning at position id + 1. On most systems using IEEE 64-bit floating- 589 ;; point arithmetic, this code works correctly so long as d is less than 590 ;; approximately 1.18 x 10^7. If 80-bit arithmetic can be employed, this limit 591 ;; is significantly higher. Whatever arithmetic is used, results for a given 592 ;; position id can be checked by repeating with id-1 or id+1, and verifying 593 ;; that the hex digits perfectly overlap with an offset of one, except possibly 594 ;; for a few trailing digits. The resulting fractions are typically accurate 595 ;; to at least 11 decimal digits, and to at least 9 hex digits. 596 ;; David H. Bailey 2006-09-08 597 598 (let ((ihex (lambda (x nhx chx) 599 ;; This returns, in chx, the first nhx hex digits of the fraction of x. 600 (do ((y (abs x)) 601 (hx "0123456789ABCDEF") 602 (i 0 (+ i 1))) 603 ((= i nhx) chx) 604 (set! y (* 16.0 (- y (floor y)))) 605 (set! (chx i) (hx (floor y)))))) 606 (series (lambda (m id) 607 ;; This routine evaluates the series sum_k 16^(id-k)/(8*k+m) using the modular exponentiation technique. 608 (let ((expm (let ((ntp 25)) 609 (let ((tp1 0) 610 (tp (make-vector ntp))) 611 (lambda (p ak) 612 ;; expm = 16^p mod ak. This routine uses the left-to-right binary exponentiation scheme. 613 ;; If this is the first call to expm, fill the power of two table tp. 614 (when (= tp1 0) 615 (set! tp1 1) 616 (set! (tp 0) 1.0) 617 (do ((i 1 (+ i 1))) 618 ((= i ntp)) 619 (set! (tp i) (* 2.0 (tp (- i 1)))))) 620 621 (if (= ak 1.0) 622 0.0 623 (let ((pl -1)) 624 ;; Find the greatest power of two less than or equal to p. 625 (do ((i 0 (+ i 1))) 626 ((or (not (= pl -1)) 627 (= i ntp))) 628 (if (> (tp i) p) 629 (set! pl i))) 630 631 (if (= pl -1) (set! pl ntp)) 632 (let ((pt (tp (- pl 1))) 633 (p1 p) 634 (r 1.0)) 635 ;; Perform binary exponentiation algorithm modulo ak. 636 637 (do ((j 1 (+ j 1))) 638 ((> j pl) r) 639 (when (>= p1 pt) 640 (set! r (* 16.0 r)) 641 (set! r (- r (* ak (floor (/ r ak))))) 642 (set! p1 (- p1 pt))) 643 (set! pt (* 0.5 pt)) 644 (when (>= pt 1.0) 645 (set! r (* r r)) 646 (set! r (- r (* ak (floor (/ r ak)))))))))))))) 647 (eps 1e-17) 648 (s 0.0)) 649 (do ((k 0 (+ k 1))) 650 ((= k id)) 651 (let* ((ak (+ (* 8 k) m)) 652 (t (expm (- id k) ak))) 653 (set! s (+ s (/ t ak))) 654 (set! s (- s (floor s))))) 655 656 ;; Compute a few terms where k >= id. 657 (do ((happy #f) 658 (k id (+ k 1))) 659 ((or (> k (+ id 100)) happy) s) 660 (let ((t (/ (expt 16.0 (- id k)) (+ (* 8 k) m)))) 661 (set! happy (< t eps)) 662 (set! s (+ s t)) 663 (set! s (- s (floor s))))))))) 664 (lambda (id) 665 ;; id is the digit position. Digits generated follow immediately after id. 666 (let ((chx (make-string 10 #\space)) 667 (pid (let ((s1 (series 1 id)) 668 (s2 (series 4 id)) 669 (s3 (series 5 id)) 670 (s4 (series 6 id))) 671 (- (+ (* 4.0 s1) (* -2.0 s2)) s3 s4)))) 672 (ihex (- (+ 1.0 pid) (floor pid)) 10 chx) 673 chx)))) 674 675(define (test-digits) 676 (do ((i 0 (+ i 1))) 677 ((= i 5)) 678 (show-digits-of-pi-starting-at-digit (* i 1000)))) 679 680(test-digits) 681 682 683(exit) 684 685