diff options
Diffstat (limited to 'Transcen.scm')
| -rw-r--r-- | Transcen.scm | 136 |
1 files changed, 124 insertions, 12 deletions
diff --git a/Transcen.scm b/Transcen.scm index b0d1a2b..dd869a7 100644 --- a/Transcen.scm +++ b/Transcen.scm @@ -1,4 +1,4 @@ -;; Copyright (C) 1992, 1993, 1995, 1997 Free Software Foundation, Inc. +;; Copyright (C) 1992, 1993, 1995, 1997, 2005 Free Software Foundation, Inc. ;; ;; This program is free software; you can redistribute it and/or modify ;; it under the terms of the GNU General Public License as published by @@ -40,6 +40,8 @@ ;;;; "Transcen.scm", Complex trancendental functions for SCM. ;;; Author: Jerry D. Hedden. +;;;; 2005-05 SRFI-70 extensions. +;;; Author: Aubrey Jaffer (define compile-allnumbers #t) ;for HOBBIT compiler @@ -63,17 +65,6 @@ ($sqrt z)) (make-polar ($sqrt (magnitude z)) (/ (angle z) 2)))) -(define expt - (let ((integer-expt integer-expt)) - (lambda (z1 z2) - (cond ((zero? z1) (if (zero? z2) 1 0)) - ((exact? z2) - (integer-expt z1 z2)) - ((and (real? z2) (real? z1) (>= z1 0)) - ($expt z1 z2)) - (else - (exp (* z2 (log z1)))))))) - (define (sinh z) (if (real? z) ($sinh z) (let ((x (real-part z)) (y (imag-part z))) @@ -137,3 +128,124 @@ (if (real? z) ($atan z) (/ (log (/ (- +i z) (+ +i z))) +2i)) ($atan2 z (car y)))) + +;;;; SRFI-70 +(define expt + (let ((integer-expt integer-expt)) + (lambda (z1 z2) + (cond ((and (exact? z2) (not (and (zero? z1) (not (positive? z2))))) + (integer-expt z1 z2)) + ((and (real? z2) (real? z1) (positive? z1)) + ($expt z1 z2)) + (else + (exp (* (if (zero? z1) (real-part z2) z2) (log z1)))))))) + +(set! quotient + (let ((integer-quotient quotient)) + (lambda (x1 x2) + (if (and (exact? x1) (exact? x2)) + (integer-quotient x1 x2) + (truncate (/ x1 x2)))))) + +(set! remainder + (let ((integer-remainder remainder)) + (lambda (x1 x2) + (if (and (exact? x1) (exact? x2)) + (integer-remainder x1 x2) + (- x1 (* x2 (quotient x1 x2))))))) + +(set! modulo + (let ((integer-modulo modulo)) + (lambda (x1 x2) + (if (and (exact? x1) (exact? x2)) + (integer-modulo x1 x2) + (- x1 (* x2 (floor (/ x1 x2)))))))) + +(define (infinite? z) (and (= z (* 2 z)) (not (zero? z)))) +(define (finite? z) (not (infinite? z))) + +(define (invintp f1 f2 f3) + (define f1^2 (* f1 f1)) + (define f2^2 (* f2 f2)) + (define f3^2 (expt f3 2)) + (let ((c (+ (* -3 f1^2 f2) + (* 3 f1 f2^2) + (* (- (* 2 f1^2) f2^2) f3) + (* (- f2 (* 2 f1)) f3^2))) + (b (+ (- f1^2 (* 2 f2^2)) f3^2)) + (a (- (* 2 f2) f1 f3))) + (define disc (- (* b b) (* 4 a c))) + (if (negative? (real-part disc)) + (/ b -2 a) + (let ((sqrt-disc (sqrt disc))) + (define root+ (/ (- sqrt-disc b) 2 a)) + (define root- (/ (+ sqrt-disc b) -2 a)) + (if (< (magnitude (- root+ f1)) (magnitude (- root- f1))) + root+ + root-))))) + +(define (extrapolate-0 fs) + (define n (length fs)) + (define (choose n k) + (do ((kdx 1 (+ 1 kdx)) + (prd 1 (/ (* (- n kdx -1) prd) kdx))) + ((> kdx k) prd))) + (do ((k 1 (+ 1 k)) + (lst fs (cdr lst)) + (L 0 (+ (* -1 (expt -1 k) (choose n k) (car lst)) L))) + ((null? lst) L))) + +(define (sequence->limit proc sequence) + (define lval (proc (car sequence))) + (if (finite? lval) + (let ((val (proc (cadr sequence)))) + (define h_n*nsamps (* (length sequence) (magnitude (- val lval)))) + (if (finite? val) + (let loop ((sequence (cddr sequence)) + (fxs (list val lval)) + (trend #f) + (ldelta (- val lval)) + (jdx (+ -1 (length sequence)))) + (cond ((null? sequence) + (case trend + ((diverging) (and (real? val) (* ldelta 1/0))) + ((bounded) (invintp val lval (caddr fxs))) + (else (cond ((zero? ldelta) val) + ((not (real? val)) #f) + (else (extrapolate-0 fxs)))))) + (else + (set! lval val) + (set! val (proc (car sequence))) + (if (finite? val) + (let ((delta (- val lval))) + (define h_j (/ h_n*nsamps jdx)) + (cond ((case trend + ((converging) (<= (magnitude delta) h_j)) + ((bounded) (<= (magnitude ldelta) (magnitude delta))) + ((diverging) (>= (magnitude delta) h_j)) + (else #f)) + (loop (cdr sequence) (cons val fxs) trend delta (+ -1 jdx))) + (trend #f) + (else + (loop (cdr sequence) (cons val fxs) + (cond ((> (magnitude delta) h_j) 'diverging) + ((< (magnitude ldelta) (magnitude delta)) 'bounded) + (else 'converging)) + delta (+ -1 jdx))))) + (and (eq? trend 'diverging) val))))) + (and (real? val) val))) + (and (real? lval) lval))) + +(define (limit proc x1 x2 . k) + (set! k (if (null? k) 8 (car k))) + (cond ((not (finite? x2)) (slib:error 'limit 'infinite 'x2 x2)) + ((not (finite? x1)) + (or (positive? (* x1 x2)) (slib:error 'limit 'start 'mismatch x1 x2)) + (limit (lambda (x) (proc (/ x))) 0.0 (/ x2) k)) + ((= x1 (+ x1 x2)) (slib:error 'limit 'null 'range x1 (+ x1 x2))) + (else (let ((dec (/ x2 k))) + (do ((x (+ x1 x2 0.0) (- x dec)) + (cnt (+ -1 k) (+ -1 cnt)) + (lst '() (cons x lst))) + ((negative? cnt) + (sequence->limit proc (reverse lst)))))))) |