|
23 | 23 |
|
24 | 24 | (define (taylor-coefficients batch brfs vars transforms-to-try) |
25 | 25 | (define expander (expand-taylor! batch)) |
26 | | - (define taylor-coeffs |
27 | | - (for*/list ([var (in-list vars)] |
28 | | - #:do [(define taylorer (taylor var batch))] |
29 | | - [transform-type transforms-to-try]) |
30 | | - (match-define (list name f finv) transform-type) |
31 | | - (define replacer (batch-replace-expression! batch var (f var))) |
32 | | - (for/list ([brf (in-list brfs)]) |
33 | | - (taylorer (expander (reducer (replacer brf))))))) |
34 | | - taylor-coeffs) |
| 26 | + (for*/list ([var (in-list vars)] |
| 27 | + #:do [(define taylorer (taylor var batch))] |
| 28 | + [transform-type transforms-to-try]) |
| 29 | + (match-define (list name f finv) transform-type) |
| 30 | + (define replacer (batch-replace-expression! batch var (f var))) |
| 31 | + (for/list ([brf (in-list brfs)]) |
| 32 | + (taylorer (expander (reducer (replacer brf))))))) |
35 | 33 |
|
36 | 34 | (define (approximate taylor-approxs |
37 | 35 | batch |
|
235 | 233 | (define (series-ref s n) |
236 | 234 | (define cache (series-cache s)) |
237 | 235 | (define builder (series-f s)) |
238 | | - (define fetch (λ (i) (dvector-ref cache i))) |
| 236 | + (define (fetch i) |
| 237 | + (dvector-ref cache i)) |
239 | 238 | (when (>= n (dvector-length cache)) |
240 | 239 | (for ([i (in-range (dvector-length cache) (add1 n))]) |
241 | 240 | (define value (reducer (adder (builder fetch i)))) |
242 | 241 | (dvector-set! cache i value))) |
243 | 242 | (dvector-ref cache n)) |
244 | 243 |
|
245 | | -(define (series-function s) |
246 | | - (λ (n) (series-ref s n))) |
| 244 | +(define ((series-function s) n) |
| 245 | + (series-ref s n)) |
247 | 246 |
|
248 | 247 | (define (taylor-add left right) |
249 | 248 | ;(-> term? term? term?) |
|
331 | 330 | (define offset (series-offset normalized)) |
332 | 331 | (define coeffs (series-function normalized)) |
333 | 332 | (define offset* (+ offset (modulo (- offset) n))) |
334 | | - (if (= offset offset*) |
335 | | - normalized |
336 | | - (let ([cache (make-dvector 2)]) ;; never called more than twice |
337 | | - (define (coeffs* i) |
338 | | - (unless (and (> (dvector-capacity cache) i) (dvector-ref cache i)) |
339 | | - (define res |
340 | | - (match i |
341 | | - [0 |
342 | | - (adder (make-sum (for/list ([j (in-range (modulo offset n))]) |
343 | | - `(* ,(coeffs j) (pow ,var ,(+ j (modulo (- offset) n)))))))] |
344 | | - [_ |
345 | | - #:when (< i n) |
346 | | - (adder 0)] |
347 | | - [_ (coeffs (+ (- i n) (modulo offset n)))])) |
348 | | - (dvector-set! cache i res)) |
349 | | - (dvector-ref cache i)) |
350 | | - (make-series offset* (λ (f i) (deref (coeffs* i))))))) |
| 333 | + (cond |
| 334 | + [(= offset offset*) normalized] |
| 335 | + [else |
| 336 | + (define cache (make-dvector 2)) ;; never called more than twice |
| 337 | + (define (coeffs* i) |
| 338 | + (unless (and (> (dvector-capacity cache) i) (dvector-ref cache i)) |
| 339 | + (define res |
| 340 | + (match i |
| 341 | + [0 |
| 342 | + (adder (make-sum (for/list ([j (in-range (modulo offset n))]) |
| 343 | + `(* ,(coeffs j) (pow ,var ,(+ j (modulo (- offset) n)))))))] |
| 344 | + [_ |
| 345 | + #:when (< i n) |
| 346 | + (adder 0)] |
| 347 | + [_ (coeffs (+ (- i n) (modulo offset n)))])) |
| 348 | + (dvector-set! cache i res)) |
| 349 | + (dvector-ref cache i)) |
| 350 | + (make-series offset* (λ (f i) (deref (coeffs* i))))])) |
351 | 351 |
|
352 | 352 | (define (taylor-sqrt var num) |
353 | 353 | ;(-> symbol? term? term?) |
|
423 | 423 | ;(-> (-> number? batchref?) term?) |
424 | 424 | (make-series 0 |
425 | 425 | (λ (f n) |
426 | | - (if (zero? n) |
427 | | - `(exp ,(coeffs 0)) |
428 | | - (let* ([coeffs* (list->vector (map coeffs (range 1 (+ n 1))))] |
429 | | - [nums (for/list ([i (in-range 1 (+ n 1))] |
430 | | - [coeff (in-vector coeffs*)] |
431 | | - #:unless (equal? (deref coeff) 0)) |
432 | | - i)]) |
433 | | - `(* (exp ,(coeffs 0)) |
434 | | - (+ ,@(for/list ([p (all-partitions n (sort nums >))]) |
435 | | - `(* ,@(for/list ([(count num) (in-dict p)]) |
436 | | - `(/ (pow ,(vector-ref coeffs* (- num 1)) ,count) |
437 | | - ,(factorial count)))))))))))) |
| 426 | + (cond |
| 427 | + [(zero? n) `(exp ,(coeffs 0))] |
| 428 | + [else |
| 429 | + (define coeffs* (list->vector (map coeffs (range 1 (+ n 1))))) |
| 430 | + (define nums |
| 431 | + (for/list ([i (in-range 1 (+ n 1))] |
| 432 | + [coeff (in-vector coeffs*)] |
| 433 | + #:unless (equal? (deref coeff) 0)) |
| 434 | + i)) |
| 435 | + `(* (exp ,(coeffs 0)) |
| 436 | + (+ ,@(for/list ([p (all-partitions n (sort nums >))]) |
| 437 | + `(* ,@(for/list ([(count num) (in-dict p)]) |
| 438 | + `(/ (pow ,(vector-ref coeffs* (- num 1)) ,count) |
| 439 | + ,(factorial count)))))))])))) |
438 | 440 |
|
439 | 441 | (define (taylor-sin coeffs) |
440 | 442 | ;(-> (-> number? batchref?) term?) |
441 | 443 | (make-series 0 |
442 | 444 | (λ (f n) |
443 | | - (if (zero? n) |
444 | | - 0 |
445 | | - (let* ([coeffs* (list->vector (map coeffs (range 1 (+ n 1))))] |
446 | | - [nums (for/list ([i (in-range 1 (+ n 1))] |
447 | | - [coeff (in-vector coeffs*)] |
448 | | - #:unless (equal? (deref coeff) 0)) |
449 | | - i)]) |
450 | | - `(+ ,@(for/list ([p (all-partitions n (sort nums >))]) |
451 | | - (if (= (modulo (apply + (map car p)) 2) 1) |
452 | | - `(* ,(if (= (modulo (apply + (map car p)) 4) 1) 1 -1) |
453 | | - ,@(for/list ([(count num) (in-dict p)]) |
454 | | - `(/ (pow ,(vector-ref coeffs* (- num 1)) ,count) |
455 | | - ,(factorial count)))) |
456 | | - 0)))))))) |
| 445 | + (cond |
| 446 | + [(zero? n) 0] |
| 447 | + [else |
| 448 | + (define coeffs* (list->vector (map coeffs (range 1 (+ n 1))))) |
| 449 | + (define nums |
| 450 | + (for/list ([i (in-range 1 (+ n 1))] |
| 451 | + [coeff (in-vector coeffs*)] |
| 452 | + #:unless (equal? (deref coeff) 0)) |
| 453 | + i)) |
| 454 | + `(+ ,@(for/list ([p (all-partitions n (sort nums >))]) |
| 455 | + (if (= (modulo (apply + (map car p)) 2) 1) |
| 456 | + `(* ,(if (= (modulo (apply + (map car p)) 4) 1) 1 -1) |
| 457 | + ,@(for/list ([(count num) (in-dict p)]) |
| 458 | + `(/ (pow ,(vector-ref coeffs* (- num 1)) ,count) |
| 459 | + ,(factorial count)))) |
| 460 | + 0)))])))) |
457 | 461 |
|
458 | 462 | (define (taylor-cos coeffs) |
459 | 463 | ;(-> (-> number? batchref?) term?) |
460 | 464 | (make-series 0 |
461 | 465 | (λ (f n) |
462 | | - (if (zero? n) |
463 | | - 1 |
464 | | - (let* ([coeffs* (list->vector (map coeffs (range 1 (+ n 1))))] |
465 | | - [nums (for/list ([i (in-range 1 (+ n 1))] |
466 | | - [coeff (in-vector coeffs*)] |
467 | | - #:unless (equal? (deref coeff) 0)) |
468 | | - i)]) |
469 | | - `(+ ,@(for/list ([p (all-partitions n (sort nums >))]) |
470 | | - (if (= (modulo (apply + (map car p)) 2) 0) |
471 | | - `(* ,(if (= (modulo (apply + (map car p)) 4) 0) 1 -1) |
472 | | - ,@(for/list ([(count num) (in-dict p)]) |
473 | | - `(/ (pow ,(vector-ref coeffs* (- num 1)) ,count) |
474 | | - ,(factorial count)))) |
475 | | - 0)))))))) |
| 466 | + (cond |
| 467 | + [(zero? n) 1] |
| 468 | + [else |
| 469 | + (define coeffs* (list->vector (map coeffs (range 1 (+ n 1))))) |
| 470 | + (define nums |
| 471 | + (for/list ([i (in-range 1 (+ n 1))] |
| 472 | + [coeff (in-vector coeffs*)] |
| 473 | + #:unless (equal? (deref coeff) 0)) |
| 474 | + i)) |
| 475 | + `(+ ,@(for/list ([p (all-partitions n (sort nums >))]) |
| 476 | + (if (= (modulo (apply + (map car p)) 2) 0) |
| 477 | + `(* ,(if (= (modulo (apply + (map car p)) 4) 0) 1 -1) |
| 478 | + ,@(for/list ([(count num) (in-dict p)]) |
| 479 | + `(/ (pow ,(vector-ref coeffs* (- num 1)) ,count) |
| 480 | + ,(factorial count)))) |
| 481 | + 0)))])))) |
476 | 482 |
|
477 | 483 | ;; This is a hyper-specialized symbolic differentiator for log(f(x)) |
478 | 484 |
|
|
522 | 528 | (define base |
523 | 529 | (make-series 0 |
524 | 530 | (λ (f n) |
525 | | - (if (zero? n) |
526 | | - `(log ,(maybe-negate (coeffs 0))) |
527 | | - (let ([tmpl (logcompute n)]) |
528 | | - `(/ (+ ,@(for/list ([term tmpl]) |
529 | | - (match-define `(,coeff ,k ,ps ...) term) |
530 | | - `(* ,coeff |
531 | | - (/ (* ,@(for/list ([i (in-naturals 1)] |
532 | | - [p ps]) |
533 | | - (if (= p 0) |
534 | | - 1 |
535 | | - `(pow (* ,(factorial i) ,(coeffs i)) ,p)))) |
536 | | - (exp (* ,(- k) ,(f 0))))))) |
537 | | - ,(factorial n))))))) |
| 531 | + (cond |
| 532 | + [(zero? n) `(log ,(maybe-negate (coeffs 0)))] |
| 533 | + [else |
| 534 | + (define tmpl (logcompute n)) |
| 535 | + `(/ (+ ,@(for/list ([term tmpl]) |
| 536 | + (match-define `(,coeff ,k ,ps ...) term) |
| 537 | + `(* ,coeff |
| 538 | + (/ (* ,@(for/list ([i (in-naturals 1)] |
| 539 | + [p ps]) |
| 540 | + (if (= p 0) |
| 541 | + 1 |
| 542 | + `(pow (* ,(factorial i) ,(coeffs i)) ,p)))) |
| 543 | + (exp (* ,(- k) ,(f 0))))))) |
| 544 | + ,(factorial n))])))) |
538 | 545 |
|
539 | 546 | (if (zero? shift) |
540 | 547 | base |
|
0 commit comments