(module math.flonum.fpvector (inline-build-fpvector build-fpvector (inline-fpvector-map check-fpvector-lengths!) fpvector-map ;; Construction unsafe-fpvector-copy! fpvector-copy! ;; Conversion list->fpvector fpvector->list vector->fpvector fpvector->vector ;; Pointwise operations fpvector-scale fpvector-sqr fpvector-sqrt fpvector-abs fpvector+ fpvector* fpvector- fpvector/ fpvector-min fpvector-max ;; Sum fpvector-sum fpvector-sums fpsum) (import scheme chicken.flonum chicken.fixnum srfi-4 chicken.type (only chicken.base include-relative assert case-lambda error exact-integer? unless sub1 void when let*-values) (only chicken.format format) (only miscmacros define-syntax-rule) math.flonum.functions math.racket-shim math.flonum.utils) (define-syntax inline-build-fpvector (syntax-rules () [(_ n-expr f-expr) (let* ([xs (make-f64vector n-expr)] [n (f64vector-length xs)]) ;; (define-syntax new-f ;; (syntax-rules () ;; [(_ i) ;; (let ([x (f-expr i)]) ;; (if (flonum? x) x (error 'build-flvector "arg is not a flonum" x)))])) (let loop ((i 0)) (if (fx< i n) (begin (f64vector-set! xs i (f-expr i)) (loop (fx+ i 1)))) xs))])) ;; (define-syntax fpvector-mapfill! ;; (syntax-rules () ;; [(_ xs n f-expr) ;; (let loop ((i 0)) ;; (if (fx< i n) ;; (begin (f64vector-set! xs i (f-expr (f64vector-ref xs i))) ;; (loop (fx+ i 1)))) ;; xs)])) (define-syntax inline-fpvector-map (syntax-rules () [(_ fx-expr xs-expr) (let ([n (f64vector-length xs-expr)]) (assert (f64vector? xs-expr)) (define ys (make-f64vector n)) (let loop ((i 0)) (if (fx< i n) (begin (f64vector-set! ys i (fx-expr (f64vector-ref xs-expr i))) (loop (fx+ i 1)))) ys))] [(_ f-expr xs-expr xss-expr ...) (let ([n (f64vector-length xs-expr)]) (let-syntax ((assert-vector (syntax-rules () ((_ v) (assert (f64vector? v) 'fpvector-map "not an f64vector" v))))) (assert-vector xs-expr) (assert-vector xss-expr) ...) (check-fpvector-lengths! 'fpvector-map n xss-expr ...) (define ys (make-f64vector n)) (let loop ((i 0)) (if (fx< i n) (begin (f64vector-set! ys i (f-expr (f64vector-ref xs-expr i) (f64vector-ref xss-expr i) ...)) (loop (fx+ i 1)))) ys) )])) (: build-flvector (fixnum (float -> float) -> f64vector)) (define (build-fpvector n f) (inline-build-fpvector n f)) (: fpvector-map (procedure f64vector #!optional f64vector #!rest (list-of f64vector) -> f64vector)) (define fpvector-map (case-lambda [(f xs) (inline-fpvector-map f xs)] [(f xs ys) (inline-fpvector-map f xs ys)] [(f xs ys . yss) (define n (f64vector-length xs)) (apply check-fpvector-lengths! 'fpvector-map n ys yss) (inline-build-fpvector n (lambda (i) (apply f (f64vector-ref xs i) (f64vector-ref ys i) (map (lambda (ys) (f64vector-ref ys i)) yss))))])) ;; =================================================================================================== ;; flvector-copy! (: unsafe-fpvector-copy! (f64vector integer f64vector integer integer -> void)) (define (unsafe-fpvector-copy! dest dest-start src src-start src-end) (let loop ([i dest-start] [j src-start]) (when (fx< j src-end) (f64vector-set! dest i (f64vector-ref src j)) (loop (fx+ i 1) (fx+ j 1))))) (: fpvector-copy! (f64vector integer f64vector #!optional integer integer -> void)) (define fpvector-copy! (case-lambda [(dest dest-start src) (fpvector-copy! dest dest-start src 0 (f64vector-length src))] [(dest dest-start src src-start) (fpvector-copy! dest dest-start src src-start (f64vector-length src))] [(dest dest-start src src-start src-end) (define dest-len (f64vector-length dest)) (define src-len (f64vector-length src)) (cond [(or (< dest-start 0) (> dest-start dest-len)) (error 'fpvector-copy! (format "expected index <= ~a" dest-len) dest-start)] [(or (< src-start 0) (> src-start src-len)) (error 'fpvector-copy! (format "expected index <= ~a" src-len) src-start)] [(or (< src-end 0) (> src-end src-len)) (error 'fpvector-copy! (format "expected index <= ~a" src-len) src-end)] [(< src-end src-start) (error 'fpvector-copy! "ending index is smaller than starting index")] [(< (- dest-len dest-start) (- src-end src-start)) (error 'fpvector-copy! "not enough room in target vector")] [else (unsafe-fpvector-copy! dest dest-start src src-start src-end)])])) ;; =================================================================================================== ;; Conversion (define list->fpvector list->f64vector) (define fpvector->list f64vector->list) (: vector->fpvector (vector -> f64vector)) (define (vector->fpvector vs) (let* ((len (vector-length vs)) (xs (make-f64vector len))) (do ((i 0 (fx+ i 1))) ((= i len)) (f64vector-set! xs i (vector-ref vs i))))) (: fpvector->vector (f64vector -> vector)) (define (fpvector->vector xs) (let* ((len (f64vector-length xs)) (vs (make-vector len))) (do ((i 0 (fx+ i 1))) ((= i len)) (vector-set! vs i (f64vector-ref xs i))))) ;; =================================================================================================== ;; Pointwise operations (define-syntax-rule (lift1 f) (lambda (arr) (inline-fpvector-map f arr))) (define-syntax-rule (lift2 f) (lambda (arr0 arr1) (inline-fpvector-map f arr0 arr1))) (: fpvector-scale (f64vector float -> f64vector)) (define (fpvector-scale arr y) (inline-fpvector-map (lambda (x) (fp* x y)) arr)) (: fpvector-sqr (f64vector -> f64vector)) (define fpvector-sqr (lift1 (lambda (x) (fp* x x)))) (: fpvector-sqrt (f64vector -> f64vector)) (define fpvector-sqrt (lift1 fpsqrt)) (: fpvector-abs (f64vector -> f64vector)) (define fpvector-abs (lift1 fpabs)) (: fpvector+ (f64vector f64vector -> f64vector)) (define fpvector+ (lift2 fp+)) (: fpvector* (f64vector f64vector -> f64vector)) (define fpvector* (lift2 fp*)) (: fpvector- (f64vector #!optional f64vector -> f64vector)) (define fpvector- (case-lambda [(arr0) (inline-fpvector-map (lambda (x) (fp- 0.0 x)) arr0)] [(arr0 arr1) (inline-fpvector-map fp- arr0 arr1)])) (: fpvector/ (f64vector #!optional f64vector -> f64vector)) (define fpvector/ (case-lambda [(arr0) (inline-fpvector-map (lambda (x) (fp/ 1.0 x)) arr0)] [(arr0 arr1) (inline-fpvector-map fp/ arr0 arr1)])) (: fpvector-min (f64vector f64vector -> f64vector)) (define fpvector-min (lift2 fpmin)) (: fpvector-max (f64vector f64vector -> f64vector)) (define fpvector-max (lift2 fpmax)) ;; =================================================================================================== ;; Summation #| Algorithm adapted from: J R Shewchuk. Adaptive Precision Floating-Point Arithmetic and Fast Geometric Predicates. Discrete & Computational Geometry, 1996, vol 18, pp 305--363. |# (: flvector-sum (f64vector -> float)) ;; Returns the sum of the elements in xs in a way that incurs rounding error only once (define (fpvector-sum xs) (define n (f64vector-length xs)) ;; Vector of remainders (define rs (make-f64vector n)) ;; Loop over `xs' (let i-loop ([i 0] ;; p = Number of valid remainders in `rs' [p 0]) (cond [(fx< i n) ;; Add x consecutively to each remainder, storing the remainder of *those* additions in `rs' (let j-loop ([j 0] ;; q = Number of remainders generated by this j-loop: [q 0] [x (f64vector-ref xs i)]) (cond [(fx< j p) ;; Get the largest of x and r, or x if it's not rational (let ((r (f64vector-ref rs j))) (let*-values ([(x r) (if (fp< (fpabs x) (fpabs r)) (values r x) (values x r))] ;; Add with remainder [(z) (fp+ x r)] [(hi lo) (cond [(fprational? z) (values z (fp- r (fp- z x)))] [else (values x r)])]) (cond [(fp= lo 0.0) ;; No remainder: don't store (makes average case O(n*log(n))) (j-loop (fx+ j 1) q hi)] [else ;; Store the remainder, increment the counter (f64vector-set! rs q lo) (j-loop (fx+ j 1) (fx+ q 1) hi)])))] [else ;; Store the sum so far as the last remainder (f64vector-set! rs q x) (i-loop (fx+ i 1) (fx+ q 1))]))] [else ;; Add all the remainders (let j-loop ([j 0] [acc 0.0]) (cond [(fx< j p) (j-loop (fx+ j 1) (fp+ acc (f64vector-ref rs j)))] [else acc]))]))) (: flvector-sums (f64vector -> f64vector)) ;; Returns the partial sums of the elements in xs in a way that incurs rounding error only once ;; for each ;; This function works just like `flvector-sum', but keeps track of partial sums instead of ;; summing all the remainders at the end (define (fpvector-sums xs) (define n (f64vector-length xs)) (define rs (make-f64vector n)) (define ss (make-f64vector n)) (let i-loop ([i 0] [p 0]) (cond [(fx< i n) (let j-loop ([j 0] [q 0] [x (f64vector-ref xs i)] [s 0.0]) (cond [(fx< j p) (let ((r (f64vector-ref rs j))) (let*-values ([(x r) (if (fp< (fpabs x) (fpabs r)) (values r x) (values x r))] [(z) (fp+ x r)] [(hi lo) (cond [(fprational? z) (values z (fp- r (fp- z x)))] [else (values x r)])]) (cond [(fp= lo 0.0) (j-loop (fx+ j 1) q hi s)] [else (f64vector-set! rs q lo) (j-loop (fx+ j 1) (fx+ q 1) hi (fp+ s lo))])))] [else (f64vector-set! rs q x) (f64vector-set! ss i (fp+ s x)) (i-loop (fx+ i 1) (fx+ q 1))]))] [else ss]))) (: flsum ((list-of float) -> float)) (define (fpsum xs) (fpvector-sum (list->fpvector xs))))