This file is indexed.

/usr/share/slib/limit.scm is in slib 3b1-3.1.

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))))))))