summaryrefslogtreecommitdiffstats
path: root/Transcen.scm
diff options
context:
space:
mode:
Diffstat (limited to 'Transcen.scm')
-rw-r--r--Transcen.scm136
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))))))))