ratepsilon
- Subject: ratepsilon
- From: Henry Baker
- Date: Sun, 01 Jul 2012 12:51:20 -0700
Here's some "toy" code to play with.
(complex-rationalize (complex 1 1e-10) 1e-19) returns 1.
(complex-rationalize (complex 1 1e-10) 1e-20) returns #c(1 1e-10).
So setting ratepsilon to 1e-20 with this sort of rationalize function will do what I want.
;;; Implement Common Lisp/Maxima "rationalize" function for complex numbers.
;;; Many functions taken from "Complex Gaussian Integers for 'Gaussian Graphics' "
;;; http://home.pipeline.com/~hbaker1/Gaussian.html
(defun gaussian-integerp (z)
;;; Inexplicably left out of Common Lisp.
(assert (numberp z))
(and (integerp (realpart z)) (integerp (imagpart z))))
(defun norm (z)
;;; Inexplicably left out of Common Lisp.
;;; norm is multiplicative, i.e., norm(a*b)=norm(a)*norm(b).
;;; Almost as useful as abs, but less expensive to compute.
(realpart (* z (conjugate z))))
(defun complex-floor (z &optional (d 1))
;;; Inexplicably left out of Common Lisp.
;;; This can be optimized by utilizing the "floor" remainders directly.
(let* ((dc (conjugate d)) (dn (* d dc)) (zdc (* z dc))
(q (complex (floor (realpart zdc) dn) (floor (imagpart zdc) dn)))
(r (- z (* d q))))
(values q r)))
(defun complex-round (z &optional (d 1))
;;; Inexplicably left out of Common Lisp.
;;; This can be optimized by utilizing the "round" remainders directly.
(let* ((dc (conjugate d)) (dn (* d dc)) (zdc (* z dc))
(q (complex (round (realpart zdc) dn) (round (imagpart zdc) dn)))
(r (- z (* d q))))
(values q r)))
(defun gaussian-gcd1 (z1 z2)
;;; Compute gcd using Euclidean algorithm and least abs remainders.
(let* ((nz1 (norm z1)) (nz2 (norm z2)))
(if (< nz1 nz2) (gaussian-gcd1 z2 z1)
(if (= nz2 0) z1
(gaussian-gcd1 z2 (- z1 (* z2 (complex-round z1 z2))))))))
(defun gaussian-gcd (z1 z2)
;;; Inexplicably left out of Common Lisp.
;;; Fix up result of gaussian-gcd1 using units.
(let* ((g (gaussian-gcd1 z1 z2))
(gr (realpart g)) (gi (imagpart g)))
(cond ((plusp gr) (if (minusp gi) (* (complex 0 1) g) g))
((minusp gr) (if (plusp gi) (* (complex 0 -1) g) (- g)))
(t (abs gi)))))
(defun gaussian-evenp (z)
;;; Inexplicably left out of Common Lisp.
(evenp (+ (realpart z) (imagpart z))))
(defun gaussian-oddp (z)
;;; Inexplicably left out of Common Lisp.
(not (gaussian-evenp z)))
(defun gaussian-numerator (z)
;;; Inexplicably left out of Common Lisp.
(let* ((x (realpart z)) (y (imagpart z))
(xd (denominator x)) (yd (denominator y))
(zn (complex (* (numerator x) yd) (* (numerator y) xd)))
(zd (* xd yd)) (g (gaussian-gcd zn zd)) (r (/ zn g)))
(if (zerop (imagpart g)) r (* (complex 0 1) r))))
(defun gaussian-denominator (z)
;;; Inexplicably left out of Common Lisp.
(let* ((x (realpart z)) (y (imagpart z))
(xd (denominator x)) (yd (denominator y))
(zn (complex (* (numerator x) yd) (* (numerator y) xd)))
(zd (* xd yd)) (g (gaussian-gcd zn zd)) (r (/ zd g)))
(if (zerop (imagpart g)) r (* (complex 0 1) r))))
(defun newt1 (n ap)
(/ (+ ap (/ n ap)) 2))
(defun rat-to-cf (r &optional (floatr (float r)))
(if (integerp r) (list r)
(multiple-value-bind
(i f)
(floor r)
(print `(rat-to-cf r ,r i ,i f ,f))
(setq foo (read))
(cons i (rat-to-cf (/ 1 f) (/ 1 (- floatr i)))))))
(defun mrat-to-cf (r &optional (eps 0.5))
;;; Returns 3 values: the cf list and the rational approximation,
;;; and whether abs(r-ratio)<eps.
(if (gaussian-integerp r) (values (list r) r t)
(multiple-value-bind
(i f)
(complex-round r)
(if (>= eps 1)
(values (list i) i t)
(let* ((if (/ f))
(epsif (* eps if))
(neweps (abs (/ (* if epsif) (- 1 epsif)))))
(multiple-value-bind
(lis ratio)
(mrat-to-cf if neweps)
(let* ((rapprox (+ i (/ ratio))))
(values
(cons i lis)
rapprox
(< (abs (- rapprox r)) eps)))))))))
(defun complex-rationalize (z &optional (eps 1.0e-10))
;;; Inexplicably left out of Common Lisp.
;;; Returns 2 values: the rational approximation and the cf list.
;;; During the recursion, eps grows until approximation is good enough.
(if (>= eps 2) (values 'infinity nil)
(if (gaussian-integerp z) (values z (list z))
(multiple-value-bind
(i f)
(complex-round z)
; (print `(crationalize z ,z eps ,eps i ,i f ,f))
(if (zerop f) (values i (list i))
(let* ((if (/ f))
(epsif (* eps if))
(neweps (abs (/ (* if epsif) (- 1 epsif)))))
(multiple-value-bind
(ratio lis)
(complex-rationalize if neweps)
(let* ((iratio (if (eq ratio 'infinity) 0 (/ ratio)))
(zapprox (+ i iratio)))
(values
zapprox
(cons i lis))))))))))
(defun cf-to-rat (cf)
(if (cdr cf) (+ (car cf) (/ 1 (cf-to-rat (cdr cf))))
(car cf)))
(defparameter cfsqrt2 '(1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2))
(defparameter sqrt2 (cf-to-rat cfsqrt2))