Skip to content


Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
tree: 1314d478fa
Fetching contributors…

Cannot retrieve contributors at this time

138 lines (131 sloc) 5.827 kb
(defun make-poly (coefs &optional round)
(declare (type simple-vector coefs))
(when (zerop (length coefs))
(return-from make-poly (lambda (x) x
(if round
(ecase *float-mode*
(let ((coefs (map '(simple-array double-float 1)
(lambda (x)
(float x 1d0))
(lambda (x)
(let ((x (float x 1d0))
(acc 0d0)
(x^i 1d0))
(declare (type double-float x acc x^i))
(map nil (lambda (coef)
(incf acc (* x^i coef))
(setf x^i (* x x^i)))
(rational acc)))))
(let ((coefs (map '(simple-array single-float 1)
(lambda (x)
(float x 1s0))
(lambda (x)
(let ((x (float x 1s0))
(acc 0s0)
(x^i 1s0))
(declare (type single-float x acc x^i))
(map nil (lambda (coef)
(incf acc (* x^i coef))
(setf x^i (* x x^i)))
(rational acc))))))
(let* ((lcm (reduce #'lcm coefs :key #'denominator))
(coefs (map 'simple-vector (lambda (x)
(* x lcm))
(lambda (x)
(let* ((x (rational x))
(num (numerator x))
(denom (denominator x))
(lb-denom (integer-length (1- denom)))
(acc 0)
(num^i 1))
(assert (= denom (ash 1 lb-denom)))
(map nil (lambda (coef)
(setf acc (ash acc lb-denom))
(incf acc (* num^i coef))
(setf num^i (* num num^i)))
(/ acc (ash lcm (* lb-denom (1- (length coefs))))))))))
(defun make-dpoly (coefs &optional (diff 1) round)
(when (> diff (length coefs))
(return-from make-dpoly (make-poly #() round)))
(flet ((diff (coefs)
(let ((coefs (subseq coefs 1)))
(loop for i upfrom 0
for c across coefs
do (setf (aref coefs i) (* c (1+ i)))
finally (return coefs)))))
(loop repeat diff do (setf coefs (diff coefs)))
(make-poly coefs round)))
(defun poly-approx-error (f df d2f coefs &optional round)
(values (let ((f~ (make-poly coefs round)))
(lambda (x)
(computable-reals:-r (funcall f x)
(funcall f~ x))))
(let ((df~ (make-dpoly coefs 1 nil)))
(lambda (x)
(computable-reals:-r (funcall df x)
(funcall df~ x))))
(let ((d2f~ (make-dpoly coefs 2 nil)))
(lambda (x)
(computable-reals:-r (funcall d2f x)
(funcall d2f~ x))))))
(defun %approximation-error-extrema (f~ ddelta-f d2delta-f points)
(assert (plusp (length points)))
(setf points (sort (copy-seq points) #'< :key #'point-loc))
(let ((*precision* 64)
(worst-diff 0)
(new-extrema '())
(roots 0)
(maximin-distance 0))
(labels ((delta (x)
(abs (- (funcall f~ (rational x))
(pull-bits (funcall *loc-value* x)))))
(new-root (root lb ub)
(incf roots)
(let ((root (round-to-float root)))
(setf maximin-distance (max maximin-distance
(min (abs (- (float-bits lb)
(float-bits root)))
(abs (- (float-bits ub)
(float-bits root))))))
(flet ((maybe-push (x &aux (x (round-to-float x)))
(unless (or (<= x lb)
(>= x ub))
(setf worst-diff (max worst-diff (delta x)))
(pushnew x new-extrema))))
(maybe-push root)
(maybe-push (float-from-bits
(1+ (float-bits root))))
(maybe-push (float-from-bits
(1- (float-bits root))))))))
(setf worst-diff (delta (point-loc (aref points 0))))
(loop for i from 0
for j from 1 below (length points)
for lb = (point-loc (aref points i))
for ub = (point-loc (aref points j))
do (setf worst-diff (max worst-diff (delta ub)))
when (> (- (float-bits ub) (float-bits lb)) 1)
do (let ((root (and (> (- (float-bits ub)
(float-bits lb))
(newton lb ub ddelta-f d2delta-f))))
(when root
(new-root (round-to-float root)
lb ub)))))
(values worst-diff new-extrema maximin-distance roots)))
;; first derivative
(defvar *loc-dvalue* #'computable-reals:exp-r)
(defvar *loc-d2value* #'computable-reals:exp-r)
(defun find-error-extrema (coefs points &optional round)
(multiple-value-bind (delta-f ddelta-f d2delta-f)
(poly-approx-error *loc-value* *loc-dvalue* *loc-d2value* coefs round)
(declare (ignore delta-f))
(%approximation-error-extrema (make-poly coefs) ddelta-f d2delta-f points)))
Jump to Line
Something went wrong with that request. Please try again.