/usr/share/slib/limit.scm is in slib 3b1-5.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 | ;;; "limit.scm" Scheme Implementation of one-side limit algorithm.
;Copyright 2005 Aubrey Jaffer
;
;Permission to copy this software, to modify it, to redistribute it,
;to distribute modified versions, and to use it for any purpose is
;granted, subject to the following restrictions and understandings.
;
;1. Any copy made of this software must include this copyright notice
;in full.
;
;2. I have made no warranty or representation that the operation of
;this software will be error-free, and I am under no obligation to
;provide any services, by way of maintenance, update, or otherwise.
;
;3. In conjunction with products arising from the use of this
;material, there shall be no use of my name in any advertising,
;promotional, or sales literature without prior written consent in
;each case.
;;@code{(require 'limit)}
(define (inv-root f1 f2 f3 prec)
(define f1^2 (* f1 f1))
(define f2^2 (* f2 f2))
(define f3^2 (expt f3 2))
(require 'root) ; SLIB
(newton:find-root (lambda (f0)
(+ (- (* (expt f0 2) f1))
(* f0 f1^2)
(* (- (* 2 (expt f0 2)) (* 3 f1^2)) f2)
(* (+ (- (* 2 f0)) (* 3 f1)) f2^2)
(* (- (+ (- (expt f0 2)) (* 2 f1^2)) f2^2)
f3)
(* (+ (- f0 (* 2 f1)) f2) f3^2)))
(lambda (f0)
(+ (- (+ (* -2 f0 f1) f1^2 (* 4 f0 f2))
(* 2 f2^2)
(* 2 f0 f3))
f3^2))
f1
prec))
(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)))
;;(printf "discriminant: %g\n" disc)
(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 0.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)))
;;(printf "f(%12g)=%12g; delta=%12g hyp=%12g j=%3d %s\n" (car sequence) val (- val lval) (/ h_n*nsamps jdx) jdx (or trend ""))
(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))))))))
|