;; mpfi.scm ; Copyright (c) 2010, Jeronimo C. Pellegrini ; All rights reserved. ; ; Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following ; conditions are met: ; ; Redistributions of source code must retain the above copyright notice, this list of conditions and the following ; disclaimer. ; Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following ; disclaimer in the documentation and/or other materials provided with the distribution. ; Neither the name of the author nor the names of its contributors may be used to endorse or promote ; products derived from this software without specific prior written permission. ; ; THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS ; OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY ; AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR ; CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR ; CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR ; SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY ; THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR ; OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ; POSSIBILITY OF SUCH DAMAGE. (module mpfi * (import scheme chicken foreign) (foreign-declare "#include ") ;; Makes an empty interval (mpfi_t in C) and ;; returns a reference to it. (define make-ia-interval-empty (foreign-lambda* (c-pointer "mpfi_t") () "mpfi_t *x; \n" "x = (mpfi_t *) malloc (sizeof (mpfi_t)); \n" "mpfi_init (*x); \n" "C_return (x);")) ;; Turns a mpfr (multi-precision real, used in the external lib) ;; into a double ;; FIXME: shouldn't this be named mpfr->flonum ? (define mpfr->double (lambda (mpfr) (let ((mpfi-tmp (make-ia-interval 0 0))) ((foreign-lambda int mpfi_init_set_fr c-pointer c-pointer) mpfi-tmp mpfr) ((foreign-lambda double mpfi_get_d c-pointer) mpfi-tmp)))) ;; Makes an empty interval [lo..hi], given f64 flonums lo and hi. (define make-ia-interval-f64 (foreign-lambda* (c-pointer "mpfi_t") ((double lo) (double hi)) "mpfi_t *x; \n" "x = (mpfi_t *) malloc (sizeof (mpfi_t));\n" "mpfi_init (*x);\n" "mpfi_interv_d (*x, lo, hi);\n" "C_return (x);")) ;; Makes an empty interval. (define make-ia-interval (foreign-lambda* (c-pointer "mpfi_t") (((c-pointer "mpfi_t") interval)) "mpfi_t *x; \n" "x = (mpfi_t *) malloc (sizeof (mpfi_t));\n" "mpfi_init_set (*x, *interval);\n" "C_return (x);")) ;; Makes an empty mpfr (multi-precision real, used in the external lib) (define make-mpfr (foreign-lambda* (c-pointer "mpfr_t") () "mpfr_t *x; \n" "x = (mpfr_t *) malloc (sizeof (mpfr_t));\n" "mpfr_init (*x);\n" "C_return (x);")) ;; ia-operation-text is used by ia-define-operation. ;; ;; Given the Scheme name and the C name of the operation, and ;; lists of argument names and types, will send back the text that should ;; be used in side the procedure body. ;; ;; (ia-operation-text '(scheme-op "op" ("t1" a1) ("t2" a2))) ;; => ;; (foreign-lambda* int (((c-pointer "t1") a1) ;; ((c-pointer "t2") a2)) ;; "op (*a1 , *a2 ); \n") (define-for-syntax ia-operation-text (lambda arguments (let ((c-op (cadar arguments)) (rest (cddar arguments))) (let ((arglist (map (lambda (arg) (list (list 'c-pointer (car arg)) (cadr arg))) rest))) `(foreign-lambda* int ,arglist ,(let loop ((str (format "~a (" c-op)) (args rest)) (let ((comma (if (> (length args) 1) ", " ""))) (if (null? args) (string-append str "); \n") (loop (string-append str (format "*~a ~a" (cadar args) comma)) (cdr args)))))))))) #| Defines a simple FFI call, given Scheme/C names and argument names and types. (define T "type") (expand `(ia-define-operation ia-ia "mpfi_sub" (,T c) (,T a) (,T b))) ==> (##core#set! ia-ia (foreign-lambda* int (((c-pointer "type") c) ((c-pointer "type") a) ((c-pointer "type") b)) "mpfi_sub (*c , *a , *b ); \n")) |# (define-syntax ia-define-operation (lambda (exp rename compare) (let ((scheme-name (cadr exp)) (args (cdr exp))) (display (format "Defining: ~a~%" args)) `(,(rename 'define) ,scheme-name ,(ia-operation-text args))))) ;; ia-define-nondestructive-operation: ;; ;; Uses their destructive counterparts... (define-syntax ia-define-nondestructive-operation (syntax-rules () ((_ op! op binary) (define op (lambda (arg1 arg2) (let ((res (make-ia-interval-empty))) (op! res arg1 arg2) res)))) ((_ op! op unary) (define op (lambda (arg1) (let ((res (make-ia-interval-empty))) (op! res arg1) res)))) ((_ op! op noargs) (define op (lambda () (let ((res (make-ia-interval-empty))) (op! res) res)))))) ;; Interval Operations ;; (ia-define-operation ia+! "mpfi_add" ("mpfi_t" ROP) ("mpfi_t" OP1) ("mpfi_t" OP2)) (ia-define-operation ia+f64! "mpfi_add_d" ("mpfi_t" ROP) ("mpfi_t" OP1) (double OP2)) (ia-define-operation ia-! "mpfi_sub" ("mpfi_t" ROP) ("mpfi_t" OP1) ("mpfi_t" OP2)) (ia-define-operation ia-f64! "mpfi_sub_d" ("mpfi_t" ROP) ("mpfi_t" OP1) (double OP2)) (ia-define-operation f64-ia! "mpfi_d_sub" ("mpfi_t" ROP) (double OP1) ("mpfi_t" OP2)) (ia-define-operation ia*! "mpfi_mul" ("mpfi_t" ROP) ("mpfi_t" OP1) ("mpfi_t" OP2)) (ia-define-operation ia*f64! "mpfi_mul_d" ("mpfi_t" ROP) ("mpfi_t" OP1) (double OP2)) (ia-define-operation ia/! "mpfi_div" ("mpfi_t" ROP) ("mpfi_t" OP1) ("mpfi_t" OP2)) (ia-define-operation ia/f64! "mpfi_div_d" ("mpfi_t" ROP) ("mpfi_t" OP1) (double OP2)) (ia-define-operation f64/ia! "mpfi_d_div" ("mpfi_t" ROP) (double OP1) ("mpfi_t" OP2)) (ia-define-nondestructive-operation ia+! ia+ binary) (ia-define-nondestructive-operation ia+f64! ia+f64 binary) (ia-define-nondestructive-operation ia-! ia- binary) (ia-define-nondestructive-operation ia-f64! ia-f64 binary) (ia-define-nondestructive-operation f64-ia! f64-ia binary) (ia-define-nondestructive-operation ia*! ia* binary) (ia-define-nondestructive-operation ia/! ia/ binary) (ia-define-nondestructive-operation ia/f64! ia/f64 binary) (ia-define-operation ia-neg! "mpfi_neg" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-sqr! "mpfi_sqr" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-inv! "mpfi_inv" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-sqrt! "mpfi_sqrt" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-abs! "mpfi_abs" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-log! "mpfi_log" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-exp! "mpfi_exp" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-exp2! "mpfi_exp2" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-nondestructive-operation ia-neg! ia-neg unary) (ia-define-nondestructive-operation ia-sqr! ia-sqr unary) (ia-define-nondestructive-operation ia-inv! ia-inv unary) (ia-define-nondestructive-operation ia-sqrt! ia-sqrt unary) (ia-define-nondestructive-operation ia-abs! ia-abs unary) (ia-define-nondestructive-operation ia-log! ia-log unary) (ia-define-nondestructive-operation ia-exp! ia-exp unary) (ia-define-nondestructive-operation ia-exp2! ia-exp2 unary) (ia-define-operation ia-cos! "mpfi_cos" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-sin! "mpfi_sin" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-tan! "mpfi_tan" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-nondestructive-operation ia-cos! ia-cos unary) (ia-define-nondestructive-operation ia-sin! ia-sin unary) (ia-define-nondestructive-operation ia-tan! ia-tan unary) (ia-define-operation ia-acos! "mpfi_acos" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-asin! "mpfi_asin" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-atan! "mpfi_atan" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-nondestructive-operation ia-acos! ia-acos unary) (ia-define-nondestructive-operation ia-asin! ia-asin unary) (ia-define-nondestructive-operation ia-atan! ia-atan unary) (ia-define-operation ia-cosh! "mpfi_cosh" ("mpfi_t" COP) ("mpfi_t" OP)) (ia-define-operation ia-sinh! "mpfi_sinh" ("mpfi_t" SOP) ("mpfi_t" OP)) (ia-define-operation ia-tanh! "mpfi_tanh" ("mpfi_t" TOP) ("mpfi_t" OP)) (ia-define-nondestructive-operation ia-cosh! ia-cosh unary) (ia-define-nondestructive-operation ia-sinh! ia-sinh unary) (ia-define-nondestructive-operation ia-tanh! ia-tanh unary) (ia-define-operation ia-acosh! "mpfi_acosh" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-asinh! "mpfi_asinh" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-atanh! "mpfi_atanh" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-nondestructive-operation ia-acosh! ia-acosh unary) (ia-define-nondestructive-operation ia-asinh! ia-asinh unary) (ia-define-nondestructive-operation ia-atanh! ia-atanh unary) (ia-define-operation ia-log1p! "mpfi_log1p" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-expm1! "mpfi_expm1" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-log2! "mpfi_log2" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-operation ia-log10! "mpfi_log10" ("mpfi_t" ROP) ("mpfi_t" OP)) (ia-define-nondestructive-operation ia-log1p! ia-log1p unary) (ia-define-nondestructive-operation ia-expm1! ia-expm1 unary) (ia-define-nondestructive-operation ia-log2! ia-log2 unary) (ia-define-nondestructive-operation ia-log10! ia-log10 unary) (ia-define-operation ia-const-log2! "mpfi_const_log2" ("mpfi_t" ROP)) (ia-define-operation ia-const-pi! "mpfi_const_pi" ("mpfi_t" ROP)) (ia-define-operation ia-const-e! "mpfi_const_euler" ("mpfi_t" ROP)) (ia-define-nondestructive-operation ia-const-log2! ia--log2 noargs) (ia-define-nondestructive-operation ia-const-pi! ia-pi noargs) (ia-define-nondestructive-operation ia-const-e! ia-e noargs) ;; Comparisons (define ia< (lambda (a b) (fx< ((foreign-lambda int mpfi_cmp c-pointer c-pointer) a b) 0))) (define ia> (lambda (a b) (fx> ((foreign-lambda int mpfi_cmp c-pointer c-pointer) a b) 0))) (define ia-overlap? (lambda (a b) (fx= ((foreign-lambda int mpfi_cmp c-pointer c-pointer) a b) 0))) (define ia-nan? (foreign-lambda bool mpfi_nan_p c-pointer)) (define ia-inf? (foreign-lambda bool mpfi_inf_p c-pointer)) (define ia-zero? (foreign-lambda bool mpfi_is_zero c-pointer)) (define ia-has-zero? (foreign-lambda bool mpfi_has_zero c-pointer)) (define ia-positive? (foreign-lambda bool mpfi_is_pos c-pointer)) (define ia-strictly-positive? (foreign-lambda bool mpfi_is_strictly_pos c-pointer)) (define ia-strictly-negative? (foreign-lambda bool mpfi_is_strictly_neg c-pointer)) (define ia-negative? (foreign-lambda bool mpfi_is_neg c-pointer)) (define ia-nonpositive? (foreign-lambda bool mpfi_is_nonpos c-pointer)) (define ia-nonnegative? (foreign-lambda bool mpfi_is_nonneg c-pointer)) (define ia-bounded? (foreign-lambda bool mpfi_bounded_p c-pointer)) ;; Set (define ia-empty? (foreign-lambda bool mpfi_is_empty c-pointer)) (define ia-inside? (foreign-lambda bool mpfi_is_inside c-pointer c-pointer)) (define ia-f64-inside? (foreign-lambda bool mpfi_is_inside_d double c-pointer)) (ia-define-operation ia-intersect "mpfi_intersect" ("mpfi_t" ROP) ("mpfi_t" OP1) ("mpfi_t" OP2)) (ia-define-operation ia-union "mpfi_union" ("mpfi_t" ROP) ("mpfi_t" OP1) ("mpfi_t" OP2)) ;; Other interval functions (define ia-lo (foreign-lambda* double (((c-pointer "mpfi_t") x)) "mpfr_t y; \n" "mpfi_t z; \n" "mpfr_init (y); \n" "mpfi_get_left (y, *x); \n" "mpfi_init_set_fr (z, y); \n" "C_return (mpfi_get_d (z) );")) (define ia-hi (foreign-lambda* double (((c-pointer "mpfi_t") x)) "mpfr_t y; \n" "mpfi_t z; \n" "mpfr_init (y); \n" "mpfi_get_right (y, *x); \n" "mpfi_init_set_fr (z, y); \n" "C_return (mpfi_get_d (z) );")) (define ia->string (lambda (x) (format "[~a__~a]" (ia-lo x) (ia-hi x)))) (define-syntax make-ia-proc-returning-double (syntax-rules () ((_ ret c-name scheme-name) (define scheme-name (lambda (interval) (let ((ROP (make-mpfr))) ((foreign-lambda ret c-name c-pointer c-pointer) ROP interval) (mpfr->double ROP))))))) ;;(define ia-diameter-absolute ;; (lambda (interval) ;; (let ((ROP (make-mpfr))) ;; ((foreign-lambda int mpfi_diam_abs c-pointer c-pointer) ROP interval) ;; (mpfr->double ROP)))) (make-ia-proc-returning-double int mpfi_diam ia-diameter) (make-ia-proc-returning-double int mpfi_diam_abs ia-diameter-absolute) (make-ia-proc-returning-double int mpfi_diam_rel ia-diameter-relative) (make-ia-proc-returning-double int mpfi_mag ia-magnitude) (make-ia-proc-returning-double int mpfi_mig ia-mignitude) (make-ia-proc-returning-double void mpfi_alea ia-random) (make-ia-proc-returning-double int mpfi_mid ia-middle) )