|
298 | 298 | so we extract that case out." |
299 | 299 | (match-define (cons offset b) (normalize-series term)) |
300 | 300 | (define cache (make-dvector 10)) |
301 | | - (dvector-set! cache 0 (reducer (adder `(/ 1 ,(b 0))))) |
302 | 301 |
|
303 | 302 | (make-series n |
304 | 303 | (- offset) |
305 | 304 | cache |
306 | 305 | n* |
307 | | - `(neg (+ ,@(for/list ([i (range n*)]) |
308 | | - `(* ,(dvector-ref cache i) (/ ,(b (- n* i)) ,(b 0)))))))) |
| 306 | + (if (zero? n*) |
| 307 | + `(/ 1 ,(b 0)) |
| 308 | + `(neg (+ ,@(for/list ([i (range n*)]) |
| 309 | + `(* ,(dvector-ref cache i) (/ ,(b (- n* i)) ,(b 0))))))))) |
309 | 310 |
|
310 | 311 | (define (taylor-quotient num denom) |
311 | 312 | ;(-> term? term? term?) |
|
315 | 316 | (match-define (cons noff a) (normalize-series num)) |
316 | 317 | (match-define (cons doff b) (normalize-series denom)) |
317 | 318 | (define cache (make-dvector 10)) |
318 | | - (dvector-set! cache 0 (reducer (adder `(/ ,(a 0) ,(b 0))))) |
319 | | - |
320 | 319 | (make-series n |
321 | 320 | (- noff doff) |
322 | 321 | cache |
323 | 322 | n* |
324 | | - `(- (/ ,(a n*) ,(b 0)) |
325 | | - (+ ,@(for/list ([i (range n*)]) |
326 | | - `(* ,(dvector-ref cache i) (/ ,(b (- n* i)) ,(b 0)))))))) |
| 323 | + (if (zero? n*) |
| 324 | + `(/ ,(a 0) ,(b 0)) |
| 325 | + `(- (/ ,(a n*) ,(b 0)) |
| 326 | + (+ ,@(for/list ([i (range n*)]) |
| 327 | + `(* ,(dvector-ref cache i) (/ ,(b (- n* i)) ,(b 0))))))))) |
327 | 328 |
|
328 | 329 | (define (modulo-series var n series) |
329 | 330 | ;(-> symbol? number? term? term?) |
|
349 | 350 | ;(-> symbol? term? term?) |
350 | 351 | (match-define (cons offset* coeffs*) (modulo-series var 2 num)) |
351 | 352 | (define cache (make-dvector 10)) |
352 | | - (dvector-set! cache 0 (reducer (adder `(sqrt ,(coeffs* 0))))) |
353 | | - (dvector-set! cache 1 (reducer (adder `(/ ,(coeffs* 1) (* 2 (sqrt ,(coeffs* 0))))))) |
354 | | - |
355 | 353 | (make-series n |
356 | 354 | (/ offset* 2) |
357 | 355 | cache |
358 | 356 | n* |
359 | 357 | (cond |
| 358 | + [(zero? n*) `(sqrt ,(coeffs* 0))] |
| 359 | + [(= n* 1) `(/ ,(coeffs* 1) (* 2 (sqrt ,(coeffs* 0))))] |
360 | 360 | [(even? n*) |
361 | 361 | `(/ (- ,(coeffs* n*) |
362 | 362 | (pow ,(dvector-ref cache (/ n* 2)) 2) |
|
375 | 375 | ;(-> symbol? term? term?) |
376 | 376 | (match-define (cons offset* coeffs*) (modulo-series var 3 num)) |
377 | 377 | (define cache (make-dvector 10)) |
378 | | - (dvector-set! cache 0 (reducer (adder `(cbrt ,(coeffs* 0))))) |
379 | | - (dvector-set! cache |
380 | | - 1 |
381 | | - (reducer (adder `(/ ,(coeffs* 1) |
382 | | - (* 3 (cbrt (* ,(dvector-ref cache 0) ,(dvector-ref cache 0)))))))) |
383 | | - |
384 | | - (make-series n |
385 | | - (/ offset* 3) |
386 | | - cache |
387 | | - n* |
388 | | - `(/ (- ,(coeffs* n*) |
389 | | - ,@(for*/list ([terms (n-sum-to 3 n*)] |
390 | | - #:unless (set-member? terms n*)) |
391 | | - (match-define (list a b c) terms) |
392 | | - `(* ,(dvector-ref cache a) ,(dvector-ref cache b) ,(dvector-ref cache c)))) |
393 | | - (* 3 ,(dvector-ref cache 0) ,(dvector-ref cache 0))))) |
| 378 | + (make-series |
| 379 | + n |
| 380 | + (/ offset* 3) |
| 381 | + cache |
| 382 | + n* |
| 383 | + (cond |
| 384 | + [(zero? n*) `(cbrt ,(coeffs* 0))] |
| 385 | + [(= n* 1) `(/ ,(coeffs* 1) (* 3 (cbrt (* ,(dvector-ref cache 0) ,(dvector-ref cache 0)))))] |
| 386 | + [else |
| 387 | + `(/ (- ,(coeffs* n*) |
| 388 | + ,@(for*/list ([terms (n-sum-to 3 n*)] |
| 389 | + #:unless (set-member? terms n*)) |
| 390 | + (match-define (list a b c) terms) |
| 391 | + `(* ,(dvector-ref cache a) ,(dvector-ref cache b) ,(dvector-ref cache c)))) |
| 392 | + (* 3 ,(dvector-ref cache 0) ,(dvector-ref cache 0)))]))) |
394 | 393 |
|
395 | 394 | (define (taylor-pow coeffs n) |
396 | 395 | ;(-> term? number? term?) |
|
423 | 422 | (define (taylor-exp coeffs) |
424 | 423 | ;(-> (-> number? batchref?) term?) |
425 | 424 | (define cache (make-dvector 10)) |
426 | | - (dvector-set! cache 0 (reducer (adder `(exp ,(coeffs 0))))) |
427 | | - |
428 | 425 | (make-series n |
429 | 426 | 0 |
430 | 427 | cache |
431 | 428 | n* |
432 | | - (let* ([coeffs* (list->vector (map coeffs (range 1 (+ n* 1))))] |
433 | | - [nums (for/list ([i (in-range 1 (+ n* 1))] |
434 | | - [coeff (in-vector coeffs*)] |
435 | | - #:unless (equal? (deref coeff) 0)) |
436 | | - i)]) |
437 | | - `(* (exp ,(coeffs 0)) |
438 | | - (+ ,@(for/list ([p (all-partitions n* (sort nums >))]) |
439 | | - `(* ,@(for/list ([(count num) (in-dict p)]) |
440 | | - `(/ (pow ,(vector-ref coeffs* (- num 1)) ,count) |
441 | | - ,(factorial count)))))))))) |
| 429 | + (if (zero? n) |
| 430 | + `(exp ,(coeffs 0)) |
| 431 | + (let* ([coeffs* (list->vector (map coeffs (range 1 (+ n 1))))] |
| 432 | + [nums (for/list ([i (in-range 1 (+ n 1))] |
| 433 | + [coeff (in-vector coeffs*)] |
| 434 | + #:unless (equal? (deref coeff) 0)) |
| 435 | + i)]) |
| 436 | + `(* (exp ,(coeffs 0)) |
| 437 | + (+ ,@(for/list ([p (all-partitions n (sort nums >))]) |
| 438 | + `(* ,@(for/list ([(count num) (in-dict p)]) |
| 439 | + `(/ (pow ,(vector-ref coeffs* (- num 1)) ,count) |
| 440 | + ,(factorial count))))))))))) |
442 | 441 |
|
443 | 442 | (define (taylor-sin coeffs) |
444 | 443 | ;(-> (-> number? batchref?) term?) |
445 | 444 | (define cache (make-dvector 10)) |
446 | | - (dvector-set! cache 0 (adder 0)) |
447 | | - |
448 | 445 | (make-series n |
449 | 446 | 0 |
450 | 447 | cache |
451 | 448 | n* |
452 | | - (let* ([coeffs* (list->vector (map coeffs (range 1 (+ n* 1))))] |
453 | | - [nums (for/list ([i (in-range 1 (+ n* 1))] |
454 | | - [coeff (in-vector coeffs*)] |
455 | | - #:unless (equal? (deref coeff) 0)) |
456 | | - i)]) |
457 | | - `(+ ,@(for/list ([p (all-partitions n* (sort nums >))]) |
458 | | - (if (= (modulo (apply + (map car p)) 2) 1) |
459 | | - `(* ,(if (= (modulo (apply + (map car p)) 4) 1) 1 -1) |
460 | | - ,@(for/list ([(count num) (in-dict p)]) |
461 | | - `(/ (pow ,(vector-ref coeffs* (- num 1)) ,count) |
462 | | - ,(factorial count)))) |
463 | | - 0)))))) |
| 449 | + (if (zero? n*) |
| 450 | + 0 |
| 451 | + (let* ([coeffs* (list->vector (map coeffs (range 1 (+ n* 1))))] |
| 452 | + [nums (for/list ([i (in-range 1 (+ n* 1))] |
| 453 | + [coeff (in-vector coeffs*)] |
| 454 | + #:unless (equal? (deref coeff) 0)) |
| 455 | + i)]) |
| 456 | + `(+ ,@(for/list ([p (all-partitions n* (sort nums >))]) |
| 457 | + (if (= (modulo (apply + (map car p)) 2) 1) |
| 458 | + `(* ,(if (= (modulo (apply + (map car p)) 4) 1) 1 -1) |
| 459 | + ,@(for/list ([(count num) (in-dict p)]) |
| 460 | + `(/ (pow ,(vector-ref coeffs* (- num 1)) ,count) |
| 461 | + ,(factorial count)))) |
| 462 | + 0))))))) |
464 | 463 |
|
465 | 464 | (define (taylor-cos coeffs) |
466 | 465 | ;(-> (-> number? batchref?) term?) |
467 | 466 | (define cache (make-dvector 10)) |
468 | | - (dvector-set! cache 0 (adder 1)) |
469 | | - |
470 | 467 | (make-series n |
471 | 468 | 0 |
472 | 469 | cache |
473 | 470 | n* |
474 | | - (let* ([coeffs* (list->vector (map coeffs (range 1 (+ n* 1))))] |
475 | | - [nums (for/list ([i (in-range 1 (+ n* 1))] |
476 | | - [coeff (in-vector coeffs*)] |
477 | | - #:unless (equal? (deref coeff) 0)) |
478 | | - i)]) |
479 | | - `(+ ,@(for/list ([p (all-partitions n* (sort nums >))]) |
480 | | - (if (= (modulo (apply + (map car p)) 2) 0) |
481 | | - `(* ,(if (= (modulo (apply + (map car p)) 4) 0) 1 -1) |
482 | | - ,@(for/list ([(count num) (in-dict p)]) |
483 | | - `(/ (pow ,(vector-ref coeffs* (- num 1)) ,count) |
484 | | - ,(factorial count)))) |
485 | | - 0)))))) |
| 471 | + (if (zero? n*) |
| 472 | + 1 |
| 473 | + (let* ([coeffs* (list->vector (map coeffs (range 1 (+ n* 1))))] |
| 474 | + [nums (for/list ([i (in-range 1 (+ n* 1))] |
| 475 | + [coeff (in-vector coeffs*)] |
| 476 | + #:unless (equal? (deref coeff) 0)) |
| 477 | + i)]) |
| 478 | + `(+ ,@(for/list ([p (all-partitions n* (sort nums >))]) |
| 479 | + (if (= (modulo (apply + (map car p)) 2) 0) |
| 480 | + `(* ,(if (= (modulo (apply + (map car p)) 4) 0) 1 -1) |
| 481 | + ,@(for/list ([(count num) (in-dict p)]) |
| 482 | + `(/ (pow ,(vector-ref coeffs* (- num 1)) ,count) |
| 483 | + ,(factorial count)))) |
| 484 | + 0))))))) |
486 | 485 |
|
487 | 486 | ;; This is a hyper-specialized symbolic differentiator for log(f(x)) |
488 | 487 |
|
|
0 commit comments