;;; -*- Mode:Common-Lisp; Package:ari; Base:10 -*-
;; Affine Real Interval (ARI arithmetic)
;; Based on ninterval.lisp in this directory,
;; See ../../papers/interval.tex for info on AAI.
(defpackage :ari ;uses generic arithmetic for endpoints
(:use :ga :cl)
(:shadowing-import-from
:ga
+ - / *
= /= > < <= >=
sin cos tan
atan asin acos atan2
sinh cosh tanh
asinh acosh atanh
expt log exp sqrt
min max
1- 1+ abs incf decf
numerator denominator
tocl re-intern
;; true false ;used by simpsimp
)
(:shadowing-import-from
:cl
union intersection)
(:export ari ri union intersection))
(eval-when (compile eval load)(require "ga") (provide "ari") (in-package :ari)
(in-package :ari)
(eval-when (compile eval load)
(defstruct ri mid exp)) ;structure for affine real interval
(defun ri(lo hi &key exp)
;; (format t "~% lo=~s, hi=~s, index=~s exp=~s" lo hi index exp)
(cond ((< hi lo)(error "attempt to create invalid affine real interval [~s, ~s]" lo hi))
(t (let* ((mid (/ (+ lo hi) 2))
(del (- hi mid)))
(setf exp (cons (cons (gensym) del) nil))
(make-ri :mid (/ (+ lo hi) 2):exp exp)))))
;;or, widen interval.
(defun wri(lo hi &key exp)
(cond ((< hi lo)(error "attempt to create invalid affine real interval [~s, ~s]" lo hi))
(t (let* ((mid (/ (+ lo hi) 2))
(del (ri::bu (- hi mid)))) ; bump up.
(setf exp (cons (cons (gensym) del) nil))
(make-ri :mid (/ (+ lo hi) 2):exp exp)))))
(defun into(lo hi)(ri lo hi)) ;; standard name would be ari::into
(defmethod print-object ((a ri) stream)
(format stream "[~a~a~a]" (ri-mid a) #.(code-char #x00b1) ;; unicode +-
(reduce #'cl::+ (mapcar #'(lambda(q)(abs(cdr q))) (ri-exp a))))
)
#+ignore ;; alternative
(defmethod print-object ((a ri) stream)
(let ((delta (reduce #'cl::+ (mapcar #'(lambda(q)(abs(cdr q))) (ri-exp a)))))
(if (> delta 0)
(format stream "[~a,~a]" (- (ri-mid a) delta)(+ (ri-mid a) delta ))
(format stream "[~a,~a]" (+ (ri-mid a) delta ) (- (ri-mid a) delta)))))
(defun mkstr (&rest args)
(with-output-to-string (s)(dolist (a args) (princ a s))))
(defun symb (&rest args) (values (intern (apply #'mkstr args))))
(defmacro with-ri2 (struct1 struct2 names1 names2 &body body)
(let ((gs1 (gensym))
(gs2 (gensym)))
`(let ((,gs1 ,struct1)
(,gs2 ,struct2))
(let ,(append
(mapcar #'(lambda (f field)
`(,f (,(symb "ri-" field) ,gs1)))
names1
'(mid exp))
(mapcar #'(lambda (f field)
`(,f (,(symb "ri-" field) ,gs2)))
names2
'(mid exp)))
,@body))))
(defmacro with-ri (struct1 names1 &body body)
(let ((gs1 (gensym)))
`(let ((,gs1 ,struct1))
(let ,(mapcar #'(lambda (f field)
`(,f (,(symb "ri-" field) ,gs1)))
names1
'(mid exp))
,@body))))
) ; end eval-when
(defmethod ga::two-arg-+ (r (s ari::ri)) ;adding num+interval
(with-ri s (mid exp)
(make-ri :mid (+ r mid) :exp exp)))
(defmethod ga::two-arg-+ ((s ari::ri) r) ;adding interval + num
(+ r s))
(defmethod ga::two-arg-- (r (s ari::ri))
(with-ri s (mid exp)
(make-ri :mid (- r mid) :exp exp)))
(defmethod ga::two-arg-- ( (s ari::ri) r)
(with-ri s (mid exp)
(make-ri :mid (- mid r) :exp exp)))
(defmethod ga::two-arg-- ((ri1 ri)(ri2 ri))
(+ ri1 (negate ri2)))
(defmethod negate ((r ri))
(make-ri :mid (- (ri-mid r))
:exp (mapcar #'(lambda(h)(cons (car h)(- (cdr h)))) (ri-exp r))))
(defmethod ga::two-arg-* ((s ari::ri) r) ;multiply interval * num
(with-ri s (mid exp)
(make-ri :mid (* r mid)
:exp (mapcar #'(lambda(q)(cons (car q)(* r (cdr q)))) exp))))
(defmethod ga::two-arg-* ( r(s ari::ri)) (ga::two-arg-* s r)) ;reverse the args
(defvar intindex 0)
;;new
(defmethod ga::two-arg-+ ((ri1 ri)(ri2 ri))
(with-ri2 ri1 ri2 (mid1 exp1)(mid2 exp2)
(let ((m (+ mid1 mid2))
;; implement a kind of sort-merge
(e (mapcar #'copy-list exp1)))
(loop for q in exp2 do
(let ((h (assoc (car q) e)))
; (format t "~%q,h=~s,~s" q h)
(if (null h)(setf e (cons q e))
(let ((sum (+ (cdr h)(cdr q))))
(if (= sum 0) (setf e (delete h e))
(setf (cdr h) sum))))))
;; (describe e)
(make-ri :mid m :exp e))))
;;#+ignore
(defmethod ga::two-arg-* ((r ri)(s ri))
;; a crude version.
;; (x0 + a1*e1 + ... +an*en) *( y0 + b1*e1 + ... +bm*em)
;; let A=|a1|+...+|an|
;; let B=|b1|+...+|bm|
;; then a bound could be (x0+ A*e_xx)*(y0+B*e_yy)
;; x0*y0+ (x0*B*e_yy)+(y0*A*e_xx)+ AB*e_zz
;; so we bound the result by
;; --> (x0*y0 + (A*y0+B*x0+A*B)e_new.
(with-ri2 r s (mid1 exp1)(mid2 exp2)
(let* ((m (* mid1 mid2))
(A (reduce #'cl::+ (mapcar #'(lambda(q)(abs(cdr q)))
exp1)))
(B (reduce #'cl::+ (mapcar #'(lambda(q)(abs(cdr q)))
exp2))) )
(make-ri :mid m
:exp `( (,(gensym) . ,(+ (abs(* mid1 B))(abs(* mid2 A))(* A B))))))))
;; this works, but by translating to conventional interval arithmetic.
;;not doing anything clever with dependencies
#+ignore
(defmethod ga::two-arg-* ((r ri)(s ri))
(with-ri2 r s (mid1 exp1)(mid2 exp2)
(let* ((A (reduce #'cl::+ (mapcar #'(lambda(q)(abs(cdr q)))
exp1)))
(B (reduce #'cl::+ (mapcar #'(lambda(q)(abs(cdr q)))
exp2)))
(lo1 (- mid1 A))
(hi1 (+ mid1 A))
(lo2 (- mid2 B))
(hi2 (+ mid2 B)) )
(if (and (> lo1 0)(> lo2 0))
(ri (* lo1 lo2)(* hi1 hi2))
(let ((all (list (* lo1 lo2)(* lo1 hi2)
(* hi1 lo2)(* hi1 hi2))))
;; (format t "~%A=~s, B=~s, all=" A B all)
(ri (apply #'cl::min all)(apply #'cl::max all)))))))
;; I think this is right, but not necessarily as tight as regular intervals sometimes.
;; Sometimes it is tighter, as shown by Stolfi.
(defmethod ga::two-arg-* ((r ri)(s ri))
;; maybe we can be tighter here...
;; since ei*ej is also in [-1,1], and ei*ei is in [0,1], we can do
;; this.. partial terms like a1*e1 * b1*e1 ==> (a1*b1)[0,1] ==> add a new
;; item with mid=(a1*b1)/2, del= (a1*b1)/2. We add to the middle
;; term, a1*b1/2, and introduce a variable (a1*b1)/2*e_new.
;; otherwise, accumulate A and B as above, if terms don't match.
;;
;; a1*b2 --> (a1*b1)*e_new
(with-ri2 r s (mid1 exp1)(mid2 exp2)
(let ((m (* (- mid1) mid2))
(e nil)
(e1 0)
(v1 0)
(leftover 0))
(loop for p in exp1 do ; compute quadratic terms
(setf e1 (car p) v1 (cdr p))
(loop for q in exp2 do
(let ((e2 (car q))
(v2 (cdr q)))
; (format t "~%v1=~s, v2=~s" v1 v2)
(cond ((eq e1 e2)
(let ((vv (* 1/2 (* v1 v2))))
; (format t "~%vv=~s, m=~s" vv m)
(setf e `( (,(gensym). ,vv),@ e))
(setf m (+ m vv))))
(t
(incf leftover (abs (* v1 v2)))))
)))
(unless (= leftover 0) (setf e `( (,(gensym). ,leftover),@ e)))
(+ ;; finally, add these intervals
(make-ri :mid m :exp e)
(* mid1 s) ;linear terms
(* mid2 r) ;linear terms
))));; works for squaring (ri -1/2 1/2). not (ri 1 3)
#|
[1,3] = 2 +- 1. to square it, compute
m= -4
e1=xx
v1=1
e2=xx
v2=1
vv=1/2
e = (( yy . 1/2))
m = -4+1/2 = -7/2
leftover=0
(* mid1 s) = (* mid2 r) = 2* [1,3] = [2,6] or 4+-1
add these two together to make [4,12] or 8 +-4
final result is
mid = -7/2 + 8 = -7/2+16/2 = 9/2
exp = ((yy . 1/2) (zz . 4))
By this computation,
low end should be 9/2 -1/2 -4 =0.
hi end should be 9/2+1/2+4 = 9.
A tighter bound would be [1,9].
Why can't we make it better??
We know that when zz =-1, yy=1, so low end would be 9/2+1/2-4 = 1.
But the representation has lost that information.
|#
;;hm.. from text on stolfi's site, he suggests
;; x_0*y_0 + sum( x_0*y_i+y_0*x_i)*e_i +rad(X)*rad(Y)e_new.
;; rad is the half-width. let's try it.
;; program below has bug? or maybe almost as much work and not as tight as above.
#+ignore
(defmethod ga::two-arg-* ((r ri)(s ri))
(with-ri2 r s (mid1 exp1)(mid2 exp2)
(let ((m (* mid1 mid2))
(e nil)
(e1 0)
(v1 0)
(rad1 (reduce #'cl::+ (mapcar #'(lambda(q)(abs(cdr q)))
exp1)))
(rad2 (reduce #'cl::+ (mapcar #'(lambda(q)(abs(cdr q)))
exp2)))
(mv 0))
(loop for p in exp1 do
(setf e1 (car p) v1 (cdr p) mv (* v1 mid2))
;; (setf e `((,e1 . ,(* v1 mid2)),@ e))
(loop for q in exp2 do
(let ((e2 (car q))
(v2 (cdr q)))
; (format t "~%v1=~s, v2=~s" v1 v2)
(cond ((eq e1 e2)
(let ((vv (+ (* mid1 v2) mv)))
(format t "~%vv=~s, m=~s" vv m)
(setf e `( (,e1 . ,vv),@ e)) ))))))
(format t "~%rad1=~s, rad2=~s" rad1 rad2)
(setf e `( (,(gensym) . ,(* rad1 rad2
)),@ e))
(make-ri :mid m :exp e)
)))
(defmethod ga::two-arg-expt ((s ri) (r fixnum))
;; can we do better than this?
(with-ri s (mid exp)
(let* ((del (reduce #'cl::+ (mapcar #'cdr exp)))
(lo (- mid del))
(hi (+ mid del)))
(if (and (evenp r)(<= lo 0 hi))
(ri 0 (expt (ga::two-arg-max (abs lo) (abs hi)) r))
(let ((e1 (expt lo r))
(e2 (expt hi r)))
(if (< e1 e2)(ri e1 e2 )
(ri e2 e1)))))))
;; make an affine into ordinary interval.
(defmethod aa2ia((a ri))
(let ((delta (reduce #'cl::+ (mapcar #'(lambda(q)(abs(cdr q))) (ri-exp a)))))
(if (> delta 0)
(ri::into (- (ri-mid a) delta)(+ (ri-mid a) delta )) (ri-mid a))))
(defmethod aa2ia(b) b)
;; make an ordinary interval into an affine interval
(defmethod ia2aa((a ri::ri))
(ari::ri (ri::ri-lo a)(ri::ri-hi a)))
(defmethod ia2aa(b) b)
(defmethod ga::invert((k ri))
(ia2aa(ri::invert (aa2ia k))))
(defmethod ga::two-arg-/((n ri) d) (* (/ 1 d) n))
(defmethod ga::two-arg-/(n (d ri)) (* (ga::invert d) n))
;;(defmethod ga::two-arg-/((n ri) (d ri)) (ia2aa (/ (aa2ia n)(aa2ia d))))
;; Every other method for IA may be transferred to AA by this technique,
;; though there may be other ways that involve less calculation or shorter program
;; as the version below.
(defmethod ga::two-arg-/((n ri) (d ri)) (* (ga::invert d) n))