ratepsilon



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