suzanne.soy/scribble-math
1
0
Fork 0
Code Issues 8 Pull Requests 3 Actions Packages Projects Releases Wiki Activity

Adding old changes

Browse Source
This commit is contained in:
Jens Axel Søgaard 2013-12-06 18:51:16 +01:00
parent b063d54f6e
commit 6e41c86f54
6 changed files with 984 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
bracket/info.rkt Normal file
View File

@ -0,0 +1,5 @@
#lang setup/infotab
; (define scribblings '(["scribblings/datalog.scrbl" (multi-page) (language)]))
(define compile-omit-paths '("tests"))

4
bracket/undefined.rkt Normal file
View File

@ -0,0 +1,4 @@
#lang racket/base
(provide undefined undefined?)
(define undefined 'undefined)
(define (undefined? e) (eq? e 'undefined))

66
numeric/adaptive-integration-simpson.rkt Normal file
View File

@ -0,0 +1,66 @@
#lang racket
(require rackunit)
; Adaptive Simpson's Method
; See Felleisen's post:
; https://groups.google.com/forum/?fromgroups#!topic/plt-scheme/JXwe9nZmCnQ
(provide/contract
[adaptive-simpson
(->i ((f (-> real? real?)) (L real?) (R (L) (and/c real? (>/c L))))
(#:epsilon (ε real?))
(r real?))])
(define (adaptive-simpson f L R #:epsilon [ε .000000001])
(define f@L (f L))
(define f@R (f R))
(define-values (M f@M whole) (simpson-1call-to-f f L f@L R f@R))
(asr f L f@L R f@R ε whole M f@M))
;; computationally efficient: 2 function calls per step
(define (asr f L f@L R f@R ε whole M f@M)
(define-values (leftM f@leftM left*) (simpson-1call-to-f f L f@L M f@M))
(define-values (rightM f@rightM right*) (simpson-1call-to-f f M f@M R f@R))
(define delta* (- (+ left* right*) whole))
(cond
[(<= (abs delta*) (* 10 ε)) (+ left* right* (/ delta* 15.0))]
[else (define epsilon1 (/ ε 2))
(+ (asr f L f@L M f@M epsilon1 left* leftM f@leftM)
(asr f M f@M R f@R epsilon1 right* rightM f@rightM))]))
(define (simpson-1call-to-f f L f@L R f@R)
(define M (mid L R))
(define f@M (f M))
(values M f@M (* (/ (abs (- R L)) 6.0) (+ f@L (* 4.0 f@M) f@R))))
(define (mid L R) (/ (+ L R) 2.))
;; simplistic prototype: many calls to f per step
;; (asr f L R ε whole)
;; Simpson's rule for approximating an integral
#;(define (simpson f L R)
(* (/ (- R L) 6.0) (+ (f L) (* 4.0 (f (mid L R))) (f R))))
#;(define (asr.v0 f L R ε whole)
(define M (mid L R))
(define left* (simpson f L M))
(define right* (simpson f M R))
(define delta* (- (+ left* right*) whole))
(cond
[(<= (abs delta*) (* 15 ε)) (+ left* right* (/ delta* 15))]
[else (define epsilon1 (/ ε 2))
(+ (asr f L M epsilon1 left*) (asr f M R epsilon1 right*))]))
(module* test #f
(define const=5 (lambda (x) 5))
(check-equal? (adaptive-simpson const=5 0 1) 5.0)
(check-equal? (adaptive-simpson const=5 0 10) 50.0)
(check-equal? (adaptive-simpson values 0 1) .5)
(define step (lambda (x) (if (< x 1) (sin x) 1)))
(adaptive-simpson step 0 2)
(adaptive-simpson (λ (x) (* (/ 100. (* x x)) (sin (/ 10. x)))) 1.0 3.0)
)

107
numeric/integration.rkt Normal file
View File

@ -0,0 +1,107 @@
#lang racket
;;;
;;; DEFINITE INTEGRAL
;;;
;;; References
; [NA] Numerical Analysis 7th ed. Burden and Faires
; Composite Simpson's Rule
; Purpose: Approximate the integral I of f(x) from a to b.
; Input: a,b endpoints
; n even positive integer
; Output: Approximation to I
; Notes: If f is relatively smooth over [a,b] (e.g. when [a,b] small)
; then a good approximation is found.
; URL: http://en.wikipedia.org/wiki/Simpson's_rule
;
; THEOREM
; Let f is C^4 on [a,b], n even, h=(b-a)/n and x_j=a+jh for i=0,1,...n.
; Then there exists μ in ]a,b[ for which the Composite Simpson's rule
; for n subintervals can be written with its error term as:
;
; b n/2-1 n/2
; int f(x) dx = h/3 [ f(a) + 2 sum f(x_2j) + 4 sum f(x_{2j-1}) + f(b) ] - error
; a j=1 j=1
;
; b-a (4)
; where error = ----- h^4 f (μ).
; 180
(define (simpson f a b [n 100])
(let* ([a (exact->inexact a)]
[b (exact->inexact b)]
[h (exact->inexact (/ (- b a) n))]
[sum-even 0.0]
[sum-odd 0.0])
(for ([i (in-range 1 n)])
(let* ([x (+ a (* i h))]
[fx (f x)])
(if (even? i)
(+= sum-even fx)
(+= sum-odd fx))))
(/ (* h (+ (f a) (f b) (* 2 sum-even) (* 4 sum-odd))) 3.0)))
(define (adaptive f a b tol n)
(let* ([a (* 1.0 a)]
[b (* 1.0 b)]
[tol (* 1.0 tol)]
[app 0.0])
(let* ([i 1]
[toli (* 10.0 tol)]
[ai a]
[hi (/ (- b a) 2.0)]
[fai (f a)]
[fci (f (+ a hi))]
[fbi (f b)]
[si (/ (* hi (+ fai (* 4.0 fci fbi))) 3.0)]
[li 1.0])
(let loop ([i i])
(when (> i 0)
(:= fd (f (+ ai (/ hi 2.0))))
(:= fe (f (+ ai (/ (* 3.0 hi) 2.0))))
(:= s1 (/ (* hi (+ fai (* 4.0 fd) fci)) 6.0))
(:= s2 (/ (* hi (+ fci (* 4.0 fe) fbi)) 6.0))
(define-values (v1 v2 v3 v4 v5 v6 v7 v8)
(values ai fai fci fbi hi toli si li))
(+= i -1)
(if (< (abs (- (+ s1 s2) v7)) v6)
(+= app (+ s1 s2))
(cond
[(>= v8 n)
'error-level-exceeded]
[else
(+= i 1)
(:= ai (+ v1 v5))
(:= fai v3)
(:= fci fe)
(:= fbi v4)
(:= hi (/ v5 2.0))
(:= toli (/ v6 2.0))
(:= si s2)
(:= li (+ v8 1.0))
;;; Convenient summation
(define-syntax (+= stx)
(syntax-case stx ()
[(_ id expr)
(identifier? #'id)
(syntax/loc stx (set! id (+ id expr)))]
[_ (raise-syntax-error #f "Expected (+= id expr)" stx)]))
(module* test #f
(require rackunit)
(define (error f F a b [n 100])
(- (simpson f a b n) (- (F b) (F a))))
(define-values (a b) (values 0 1))
(check-equal? (error (λ (x) 1) (λ (x) x) a b) 0.0)
(check-pred (λ (err) (<= err 0.00002)) (error (λ (x) (sin x)) (λ (x) (- (cos x))) 0 pi 20)))

545
numeric/qags.rkt Normal file
View File

@ -0,0 +1,545 @@
#lang typed/racket
(require "tagged-begin.rkt"
racket/flonum)
; Search for TODO
;;;
;;; FORTRAN RUNTIME
;;;
(define epsilon
(let ()
(define (compute-machine-epsilon)
(let loop ([n 1.0])
(define new-n (fl/ n 2.0))
(cond
[(fl= 1.0 (fl+ 1.0 (fl/ new-n 2.0)))
new-n]
[else
(loop new-n)])))
(define machine-epsilon (compute-machine-epsilon))
(λ (x) machine-epsilon)))
;;;
;;; QUADPACK ROUTINES
;;;
; subroutine qags ( f, a, b, epsabs, epsrel, result, abserr, neval, ier )
; http://orion.math.iastate.edu/burkardt/f_src/quadpack/quadpack.f90
; QAGS estimates the integral of a function.
; Discussion:
; The routine calculates an approximation RESULT to a definite integral
; I = integral of F over (A,B),
; hopefully satisfying
; || I - RESULT || <= max ( EPSABS, EPSREL * ||I|| ).
;
; Reference:
;
; R Piessens, E de Doncker-Kapenger, C W Ueberhuber, D K Kahaner,
; QUADPACK, a Subroutine Package for Automatic Integration,
; Springer Verlag, 1983
;
; Parameters:
;
; Input, external real F, the name of the function routine, of the form
; function f ( x )
; real f
; real x
; which evaluates the integrand function.
;
; Input, real A, B, the limits of integration.
;
; Input, real EPSABS, EPSREL, the absolute and relative accuracy requested.
;
; Output, real RESULT, the estimated value of the integral.
;
; Output, real ABSERR, an estimate of || I - RESULT ||.
;
; Output, integer NEVAL, the number of times the integral was evaluated.
;
; Output, integer IER, error flag.
; ier = 0 normal and reliable termination of the
; routine. it is assumed that the requested
; accuracy has been achieved.
; ier > 0 abnormal termination of the routine
; the estimates for integral and error are
; less reliable. it is assumed that the
; requested accuracy has not been achieved.
; = 1 maximum number of subdivisions allowed
; has been achieved. one can allow more sub-
; divisions by increasing the data value of
; limit in qags (and taking the according
; dimension adjustments into account).
; however, if this yields no improvement
; it is advised to analyze the integrand
; in order to determine the integration
; difficulties. if the position of a
; local difficulty can be determined (e.g.
; singularity, discontinuity within the
; interval) one will probably gain from
; splitting up the interval at this point
; and calling the integrator on the sub-
; ranges. if possible, an appropriate
; special-purpose integrator should be used,
; which is designed for handling the type
; of difficulty involved.
; = 2 the occurrence of roundoff error is detec-
; ted, which prevents the requested
; tolerance from being achieved.
; the error may be under-estimated.
; = 3 extremely bad integrand behavior occurs
; at some points of the integration
; interval.
; = 4 the algorithm does not converge. roundoff
; error is detected in the extrapolation
; table. it is presumed that the requested
; tolerance cannot be achieved, and that the
; returned result is the best which can be
; obtained.
; = 5 the integral is probably divergent, or
; slowly convergent. it must be noted that
; divergence can occur with any other value
; of ier.
; = 6 the input is invalid, because
; epsabs < 0 and epsrel < 0,
; result, abserr and neval are set to zero.
;
; Local Parameters:
;
; alist - list of left end points of all subintervals
; considered up to now
; blist - list of right end points of all subintervals
; considered up to now
; rlist(i) - approximation to the integral over
; (alist(i),blist(i))
; rlist2 - array of dimension at least limexp+2 containing
; the part of the epsilon table which is still
; needed for further computations
; elist(i) - error estimate applying to rlist(i)
; maxerr - pointer to the interval with largest error
; estimate
; errmax - elist(maxerr)
; erlast - error on the interval currently subdivided
; (before that subdivision has taken place)
; area - sum of the integrals over the subintervals
; errsum - sum of the errors over the subintervals
; errbnd - requested accuracy max(epsabs,epsrel*
; abs(result))
; *****1 - variable for the left interval
; *****2 - variable for the right interval
; last - index for subdivision
; nres - number of calls to the extrapolation routine
; numrl2 - number of elements currently in rlist2. if an
; appropriate approximation to the compounded
; integral has been obtained it is put in
; rlist2(numrl2) after numrl2 has been increased
; by one.
; small - length of the smallest interval considered
; up to now, multiplied by 1.5
; erlarg - sum of the errors over the intervals larger
; than the smallest interval considered up to now
; extrap - logical variable denoting that the routine is
; attempting to perform extrapolation i.e. before
; subdividing the smallest interval we try to
; decrease the value of erlarg.
; noext - logical variable denoting that extrapolation
; is no longer allowed (true value)
;
;
; integer, parameter :: limit = 500
; subroutine qags ( f, a, b, epsabs, epsrel, result, abserr, neval, ier )
(define: (qags [f : (Float -> Float)] [a : Float] [b : Float]
[epsabs : Float] [epsrel : Float]
[result : Float] [abserr : Float]
[neval : Integer] [ier : Integer]) : (U Float Integer)
(let/ec: return : (U Float Integer)
(define: limit : Integer 500)
(: abseps Float)
(: alist (Vectorof Float)) ; real alist(limit)
(: area Float)
(: area1 Float)
(: area12 Float)
(: area2 Float)
(: a1 Float)
(: a2 Float)
(: blist (Vectorof Float)) ; real blist(limit)
(: b1 Float)
(: b2 Float)
(: correc Float)
(: defabs Float)
(: defabs1 Float)
(: defabs2 Float)
(: dres Float)
(: elist (Vectorof Float)) ; real elist(limit)
(: erlarg Float)
(: erlast Float)
(: errbnd Float)
(: errmax Float)
(: error1 Float)
(: error2 Float)
(: erro12 Float)
(: errsum Float)
(: errtest Float)
(: extrap Float)
; TODO: real, external :: f
(: id Integer)
(: irro Integer)
(: iord (Vectorof Integer))
(define iord (make-vector limit 0))
(: iroff1 Integer)
(: iroff2 Integer)
(: iroff3 Integer)
(: jupbnd Integer)
(: k Integer)
(: ksgn Integer)
(: ktmin Integer)
(: last Integer)
(: noext Boolean) ; logical noext
(: maxerr Integer)
(: nres Integer)
(: nrmax Integer)
(: numr12 Integer)
(: resabs Float)
(: reseps Float)
(: res31a (Vectorof Float)) ; real res3la(3)
(: rlist (Vectorof Float)) ; real rlist(limit)
(: rlist2 (Vectorof Float)) ; real rlist2(52)
(: small Real)
; The dimension of rlist2 is determined by the value of
; limexp in QEXTR (rlist2 should be of dimension
; (limexp+2) at least).
;
; Test on validity of parameters.
;
(define ier 0)
(define neval 0)
(define last 0)
(define result 0.0e+00)
(define abserr 0.0e+00)
(define alist (make-vector limit 0.0))
(vector-set! alist 1 a) ; alist(1) = a
(define blist (make-vector limit 0.0))
(vector-set! blist 1 b) ; blist(1) = b
(define rlist (make-vector limit 0.0))
(vector-set! rlist 1 0.0e+00) ; rlist(1) = 0.0e+00
(define elist (make-vector limit 0.0))
(vector-set! elist 1 0.0e+00) ; elist(1) = 0.0e+00
(when (and (< epsabs 0.0) (< epsrel 0.0))
(set! ier 6)
(return ier))
;;
;; First approximation to the integral.
;;
(define ierro 0)
; TODO TODO
; call qk21 ( f, a, b, result, abserr, defabs, resabs )
;;
;; Test on accuracy.
;;
(define dres (abs result))
(define errbnd (max epsabs (* epsrel dres)))
(set! last 1)
(vector-set! rlist 1 result) ; rlist(1) = result
(vector-set! elist 1 abserr) ; elist(1) = abserr
(vector-set! iord 1 1)
; if ( abserr <= 1.0e+02 * epsilon ( defabs ) * defabs .and. &
; abserr > errbnd ) then
; ier = 2
; end if
(when (and (<= abserr (* 1.0e+02 (epsilon defabs)))
(> abserr errbnd))
(set! ier 2))
(when (= limit 1)
(set! ier 1))
(when (or (not (= ier 0))
(and (<= abserr errbnd)
(not (= abserr resabs)))
(= abserr 0.0e+00))
(goto-140))
;;
;; Initialization.
;;
; rlist2(1) = result
; errmax = abserr
; maxerr = 1
; area = result
; errsum = abserr
; abserr = huge ( abserr )
; nrmax = 1
; nres = 0
; numrl2 = 2
; ktmin = 0
; extrap = .false.
; noext = .false.
; iroff1 = 0
; iroff2 = 0
; iroff3 = 0
;
; if ( dres >= (1.0e+00-5.0e+01* epsilon ( defabs ) )*defabs ) then
; ksgn = 1
; else
; ksgn = -1
; end if
;
; do last = 2, limit
;;
;; Bisect the subinterval with the nrmax-th largest error estimate.
;;
; a1 = alist(maxerr)
; b1 = 5.0e-01*(alist(maxerr)+blist(maxerr))
; a2 = b1
; b2 = blist(maxerr)
; erlast = errmax
; call qk21 ( f, a1, b1, area1, error1, resabs, defab1 )
; call qk21 ( f, a2, b2, area2, error2, resabs, defab2 )
;;
;; Improve previous approximations to integral and error
;; and test for accuracy.
;;
; area12 = area1+area2
; erro12 = error1+error2
; errsum = errsum+erro12-errmax
; area = area+area12-rlist(maxerr)
;
; if ( defab1 == error1 .or. defab2 == error2 ) go to 15
;
; if ( abs ( rlist(maxerr) - area12) > 1.0e-05 * abs(area12) &
; .or. erro12 < 9.9e-01 * errmax ) go to 10
;
; if ( extrap ) then
; iroff2 = iroff2+1
; else
; iroff1 = iroff1+1
; end if
;
;10 continue
;
; if ( last > 10.and.erro12 > errmax ) iroff3 = iroff3+1
;
;15 continue
;
; rlist(maxerr) = area1
; rlist(last) = area2
; errbnd = max ( epsabs,epsrel*abs(area))
;;
;; Test for roundoff error and eventually set error flag.
;;
; if ( iroff1+iroff2 >= 10 .or. iroff3 >= 20 ) then
; ier = 2
; end if
;
; if ( iroff2 >= 5 ) ierro = 3
;;
;; Set error flag in the case that the number of subintervals
;; equals limit.
;;
; if ( last == limit ) then
; ier = 1
; end if
;;
;; Set error flag in the case of bad integrand behavior
;; at a point of the integration range.
;;
; if ( max ( abs(a1),abs(b2)) <= (1.0e+00+1.0e+03* epsilon ( a1 ) )* &
; (abs(a2)+1.0e+03* tiny ( a2 ) ) ) then
; ier = 4
; end if
;;
;; Append the newly-created intervals to the list.
;;
; if ( error2 <= error1 ) then
; alist(last) = a2
; blist(maxerr) = b1
; blist(last) = b2
; elist(maxerr) = error1
; elist(last) = error2
; else
; alist(maxerr) = a2
; alist(last) = a1
; blist(last) = b1
; rlist(maxerr) = area2
; rlist(last) = area1
; elist(maxerr) = error2
; elist(last) = error1
; end if
;;
;; Call QSORT to maintain the descending ordering
;; in the list of error estimates and select the subinterval
;; with nrmax-th largest error estimate (to be bisected next).
;;
; call qsort ( limit, last, maxerr, errmax, elist, iord, nrmax )
;
; if ( errsum <= errbnd ) go to 115
;
; if ( ier /= 0 ) then
; exit
; end if
;
; if ( last == 2 ) go to 80
; if ( noext ) go to 90
;
; erlarg = erlarg-erlast
;
; if ( abs(b1-a1) > small ) then
; erlarg = erlarg+erro12
; end if
;
; if ( extrap ) go to 40
;;
;; Test whether the interval to be bisected next is the
;; smallest interval.
;;
; if ( abs(blist(maxerr)-alist(maxerr)) > small ) go to 90
; extrap = .true.
; nrmax = 2
;
;40 continue
;;
;; The smallest interval has the largest error.
;; Before bisecting decrease the sum of the errors over the
;; larger intervals (erlarg) and perform extrapolation.
;;
; if ( ierro /= 3 .and. erlarg > ertest ) then
;
; id = nrmax
; jupbnd = last
;
; if ( last > (2+limit/2) ) then
; jupbnd = limit+3-last
; end if
;
; do k = id, jupbnd
; maxerr = iord(nrmax)
; errmax = elist(maxerr)
; if ( abs(blist(maxerr)-alist(maxerr)) > small ) go to 90
; nrmax = nrmax+1
; end do
;
; end if
;;
;; Perform extrapolation.
;;
;60 continue
;
; numrl2 = numrl2+1
; rlist2(numrl2) = area
; call qextr ( numrl2, rlist2, reseps, abseps, res3la, nres )
; ktmin = ktmin+1
;
; if ( ktmin > 5 .and. abserr < 1.0e-03 * errsum ) then
; ier = 5
; end if
;
; if ( abseps < abserr ) then
;
; ktmin = 0
; abserr = abseps
; result = reseps
; correc = erlarg
; ertest = max ( epsabs,epsrel*abs(reseps))
;
; if ( abserr <= ertest ) then
; exit
; end if
;
; end if
;;
;; Prepare bisection of the smallest interval.
;;
; if ( numrl2 == 1 ) then
; noext = .true.
; end if
;
; if ( ier == 5 ) then
; exit
; end if
;
; maxerr = iord(1)
; errmax = elist(maxerr)
; nrmax = 1
; extrap = .false.
; small = small*5.0e-01
; erlarg = errsum
; go to 90
;
;80 continue
;
; small = abs(b-a)*3.75e-01
; erlarg = errsum
; ertest = errbnd
; rlist2(2) = area
;
;90 continue
;
; end do
;;
;; Set final result and error estimate.
;;
; if ( abserr == huge ( abserr ) ) go to 115
; if ( ier+ierro == 0 ) go to 110
;
; if ( ierro == 3 ) then
; abserr = abserr+correc
; end if
;
; if ( ier == 0 ) ier = 3
; if ( result /= 0.0e+00.and.area /= 0.0e+00 ) go to 105
; if ( abserr > errsum ) go to 115
; if ( area == 0.0e+00 ) go to 130
; go to 110
;
;105 continue
;
; if ( abserr/abs(result) > errsum/abs(area) ) go to 115
;;
;; Test on divergence.
;;
;110 continue
;
; if ( ksgn == (-1).and.max ( abs(result),abs(area)) <= &
; defabs*1.0e-02 ) go to 130
;
; if ( 1.0e-02 > (result/area) .or. (result/area) > 1.0e+02 &
; .or. errsum > abs(area) ) then
; ier = 6
; end if
;
; go to 130
;;
;; Compute global integral sum.
;;
;115 continue
;
; result = sum ( rlist(1:last) )
;
; abserr = errsum
;
;130 continue
;
; if ( ier > 2 ) ier = ier-1
;
;140 continue
; neval = 42*last-21
(define (goto-140)
(set! neval (- (* 42 last) 21)))
;
;
; return
;end
42.0
))

257
numeric/tagged-begin.rkt Normal file
View File

@ -0,0 +1,257 @@
#lang racket
(provide tagged-begin)
;;; INTRODUCTION
; This is a little macro that resembles the Common Lisp tagbody construct
; <http://www-2.cs.cmu.edu/Groups/AI/html/hyperspec/HyperSpec/Body/speope_tagbody.html#tagbody>
; See also "Applications of Continuations" of Daniel P. Friedman.
;;; MOTIVATION
; Many algorithms is specified in an imperative manner
; in the literature [See Example 5 from Knuth]. For a no-brain-
; conversion to Scheme tagged-begin is convenient.
;;; SYNTAX
; (tagged-begin
; (<tag> | <expression>)* )
; where <tag> is a symbol and duplicate tags are not allowed.
;;; SEMANTICS
; The form evaluates the expressions in a lexical environment
; that provides functions go and return both of one argument to
; transfer control.
; The expressions in tagged-begin are evaluated sequentially.
; If no expressions are left (void) is returned.
; If an expression evaluates (go tag) then control is transfered
; to the expression following tag. The tags have lexical scope.
; The dynamic extent of tag is indefinite. An (go tag) is allowed
; to tranfer control to an outer tagbody. The call (go tag) has the
; proper tail recursive property, even in situation where the call
; syntactically is not in tail position.
; If (return <expression>) is evaluted, the value of <expression> is
; the value of the entire tagged-begin form.
;;; EXAMPLES
; See below the implementation.
;;; IMPLEMENTATION
; Tagged begin is here implemented as a syntax-case macro.
; The rewrite rule is taken from Daniel P. Friedmans
; "Applications of Continuations".
; (tagged-begin
; tag_1 e1 ... ; If the body doesn't begin with a tag
; ... ; the macro inserts a fresh one
; tag_n-1 e_n-1 ...
; tag_n en ...)
; expands to
; ((let/cc go
; (let ([return (lambda (v) (go (lambda () v)))])
; (letrec
; ([tag_1 (lambda () e1 ... (tag2))]
; ...
; [tag_n-1 (lambda () e_n-1 ... (tag_n))]
; [tag_n (lambda () e_n ... (return (void)))]
; (tag_1))))))
; where (let/cc k expr ...) is short for (call/cc (lambda (k) expr ...)))]))
(require
(for-syntax
racket
(only-in (lib "list.ss" "srfi" "1")
drop-right take-while)))
(define-syntax (tagged-begin stx)
(define tag? identifier?)
(define (non-tag? o) (not (tag? o)))
(define (generate-binding tag-exprs next-tag)
(match tag-exprs
[(list tag exprs)
(quasisyntax/loc stx [#,tag (lambda () #,@exprs (#,next-tag))])]))
(define (generate-last-binding tag-exprs return)
(match tag-exprs
[(list tag exprs)
(quasisyntax/loc stx [#,tag (lambda () #,@exprs (#,return (void)))])]))
(syntax-case stx ()
[(tagged-begin . tag/exprs-stx)
(let ([tes (syntax->list #'tag/exprs-stx)])
; introduce a dummy start-tag, if the tagged-begin starts with an expression
(when (not (tag? (car tes)))
(set! tes (cons #'start tes)))
(let* ([first-tag (car tes)]
[tag-exprs-list (let loop ([tes tes]
[rev-result '()])
(if (null? tes)
(reverse rev-result)
(let ([p tes])
(if (tag? (car p))
(loop (cdr tes)
(cons (list (car p) (take-while non-tag? (cdr p)))
rev-result))
(loop (cdr tes)
rev-result)
))))
#;(list-ec (:pairs p tes)
(if (tag? (car p)))
(list (car p) (take-while non-tag? (cdr p))))
]
[tags (map car tag-exprs-list)])
; tag-exprs-list = ( (tag_1 (e1 ...)) (tag_2 (e2 ...)) ... )
(with-syntax ([go (syntax-local-introduce (syntax/loc stx go))]
[return (syntax-local-introduce (syntax/loc stx return))])
#`((let/cc go
(let ([return (lambda (v) (go (lambda () v)))])
(letrec
(#,@(map generate-binding
(drop-right tag-exprs-list 1)
(cdr tags))
#,(generate-last-binding (last tag-exprs-list) #'return))
(#,first-tag))))))))]))
(module* test #f
; Example 1 (tagged-begin returns (void))
(displayln
(let ([i 0])
(tagged-begin
loop (set! i (+ i 1))
(when (< i 41) (go loop)))
i))
; Example 2 (tagged-begin returns 42)
(displayln
(let ([i 0])
(tagged-begin
loop (set! i (+ i 1))
(when (< i 42) (go loop))
(return i))))
; Example 3 (tagged-begin returns 43)
(displayln
(let ([i 0])
(tagged-begin
loop (set! i (+ i 1))
(go b)
a (when (< i 43) (go loop))
(return i)
b (go a))))
; Example 4 ( <http://www.emacswiki.org/cgi-bin/wiki.pl?StateMachine> )
(let ((a 0))
(tagged-begin
start
(set! a 0)
part-1
(set! a (+ a 1))
(displayln a)
(cond
((>= a 9) (go end))
((even? a) (go part-1))
(else (go part-2)))
part-2
(set! a (+ a 1))
(go part-1)
end
(displayln "We're done printing the odd numbers between 0 and 10")))
; Example 5 ( Knuth: "The Art of Computer Programming", vol1, p.176)
; Inplace inversion of a permutation represented as a vector.
(define permutation (vector 'dummy 6 2 1 5 4 3)) ; (Knuth counts from 1 not 0 :-) )
(define n (- (vector-length permutation) 1))
(define (X i) (vector-ref permutation i))
(define (X! i j) (vector-set! permutation i j))
(let ([m 0] [i 0] [j 0])
(tagged-begin
I1 ; Initialize
(set! m n)
(set! j -1)
I2 ; Next element
(set! i (X m))
(when (< i 0) (go I5))
I3 ; Invert one
(X! m j)
(set! j (- m))
(set! m i)
(set! i (X m))
I4 ; End of cycle?
(when (> i 0) (go I3))
(set! i j)
I5 ; Store final value
(X! m (- i))
I6 ; Loop on m
(set! m (- m 1))
(when (> m 0) (go I2))))
(displayln permutation)
; Example 6 (The CommonLisp Hyper Spec examples of tagbody)
(define val 'foo)
(tagged-begin
(set! val 1)
(go a)
c (set! val (+ val 4))
(go b)
(set! val (+ val 32))
a (set! val (+ val 2))
(go c)
(set! val (+ val 64))
b (set! val (+ val 8)))
(displayln val)
(define (f1 flag)
(let ((n 1))
(tagged-begin
(set! n (f2 flag (lambda () (go out))))
out
(display n))))
(define (f2 flag escape)
(if flag (escape) 2))
(displayln (f1 #f))
(displayln (f1 #t))
; Example 7
; Demonstrates lexical scoping of tagged-begins,
; and that an inner tagged-begin can use an outer tag.
(tagged-begin
a (tagged-begin
(go b))
b (return 'hello-world))
; Demonstrates that tags are lexically shadowed.
(tagged-begin
a (tagged-begin
(go b)
(return 'wrong)
b (go c))
b (return 'wrong)
c (return 'correct))
)