;;; rational.el --- Rational number in factorial-base format. ;; Copyright (C) 2000-2019 Vinicius Jose Latorre ;; Author: Vinicius Jose Latorre ;; Maintainer: Vinicius Jose Latorre ;; Keywords: data, curiosity, mathematics ;; Time-stamp: <2019/10/10 01:18:56 vinicius> ;; Version: 1.1.1 ;; X-URL: https://www.emacswiki.org/cgi-bin/wiki/ViniciusJoseLatorre ;; This file is *NOT* (yet?) part of GNU Emacs. ;; This program is free software; you can redistribute it and/or modify ;; it under the terms of the GNU General Public License as published by ;; the Free Software Foundation; either version 2, or (at your option) ;; any later version. ;; This program is distributed in the hope that it will be useful, ;; but WITHOUT ANY WARRANTY; without even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ;; GNU General Public License for more details. ;; You should have received a copy of the GNU General Public License along with ;; GNU Emacs; see the file COPYING. If not, write to the Free Software ;; Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. ;;; Commentary: ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; Introduction ;; ------------ ;; ;; This package was adapted from article: ;; ;; Error-Free Fractions ;; Peter Wayner ;; BYTE, june 1988, pages 289-298 ;; ;; This package provides a way to represent a rational number exactly with the ;; use of factorial-base format. ;; ;; To use rational, insert in your ~/.emacs: ;; ;; (require 'rational) ;; ;; For good performance, be sure to byte-compile rational.el, e.g. ;; ;; M-x byte-compile-file ;; ;; This will generate rational.elc, which will be loaded instead of ;; rational.el. ;; ;; rational was tested with GNU Emacs 20.4.1. ;; ;; Please, read the article for mathematical proof and references. ;; ;; ;; Some Experimental Results ;; ------------------------- ;; ;; As the article cited above explains (page 292): ;; ;; ...... ;; [Consider the Pascal-like code:] ;; ;; x := 1 / n; ;; for i := 1 to 30 do ;; x := (n + 1) * x - 1; ;; ;; Mathematically, the function f(x) = (n + 1)x - 1 is invariant at the point ;; x = 1/n; that is, f(x) = x for x = 1/n. On paper, then, you would expect ;; the variable x to remain unchanged after 30 interations of the loop. This ;; is the case when I used factorial-base numbers. But the standard ;; floating-point system failed badly and returned 286,331,161.6 instead of ;; 0.33333 when n was set to 3. ;; From previous discussions, you might expect the floating-point software to ;; find the correct answer at least for n = 10 because 1/10 = 0.1 exactly in ;; base 10. This is quite far from the truth: x should be equaled 0.1 but ;; turned into 2.36378547759e21 after 30 loops. All the calculations are, of ;; course, done in binary. The floating-point software finds the correct ;; answer only when n is 2. ;; The only negative aspect of the factorial-base system is the slowness of ;; the calculations. ;; ...... ;; ;; Try to execute the test: ;; ;; M-x rational-test RET ;; ;; It's displayed: ;; ;; *** Rational Test *** ;; ;; function: f(x) = (n + 1)x - 1 with x = 1/n and n = 3 ;; (after 30 interactions of the loop) ;; ;; Using rational package (factorial-base): 0 . 0 2 0 0 0 0 0 0 0 0 0 ;; Using rational package (decimal-base) : 0.3333333333333333 ;; Using floating-point : -21.0 ;; ;; ********************* ;; ;; The floating-point result may differ depending on which machine/environment ;; you run, but the result will not be 0.333333. ;; ;; ;; Using rational ;; -------------- ;; ;; Use `rational-customize' to customize rational options by typing: ;; ;; M-x rational-customize RET ;; ;; You can also bind `rational-customize' to some key, like: ;; ;; (global-set-key "\C-c\C-c" 'rational-customize) ;; ;; There are the following predefined rational constants: `rational-zero' and ;; `rational-one'. ;; ;; To create a rational number you can use `make-rational', ;; `integer-to-rational' or `rational-convert'. ;; ;; To check if an object is a rational factorial-base number, use `rationalp'. ;; ;; To set a rational number to zero or one, use `rational-zero' or ;; `rational-one' functions, respectively. ;; ;; To set a rational number to any integer, use `rational-set'. ;; ;; To copy a rational, use `rational-copy'. ;; ;; To compare two rationals, use `rational-lessequal', `rational-less' or ;; `rational-equal'. ;; ;; There are the following rational operations: `rational-add', ;; `rational-subtract', `rational-absolute', `rational-negative', ;; `rational-multbyint', `rational-divbyint', `rational-multiply' and ;; `rational-divide'. ;; ;; To translate a rational to another representation, use `rational-to-string' ;; or `rational-to-float'. ;; ;; ;; Options ;; ------- ;; ;; The variable `rational-max-size' specifies the maximum array storage for ;; factorial-base number. ;; ;; ;; Things To Change ;; ---------------- ;; ;; . At moment, nothing. Any idea? Send it! ;; ;; ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; code: ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; User Variables: ;;; Interface to the command system (defgroup rational nil "Rational numbers in factorial-base format" :link '(emacs-library-link :tag "Source Lisp File" "rational.el") :prefix "rational-" :group 'data) (defcustom rational-max-size 12 "*Maximum array storage for factorial-base number." :type 'integer :group 'rational) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Customization ;;;###autoload (defun rational-customize () "Customize rational options." (interactive) (customize-group 'rational)) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; User functions ;;;###autoload (defun make-rational (&optional int) "Create a rational factorial-base vector. The coefficients of a factorial-base number. Slot i contains a(2 - i). Slot 1 is used to carry the integer part of the number." (let ((rational (make-vector (1+ rational-max-size) 0))) (aset rational 0 'rational) (and int (aset rational 1 int)) rational)) ;;;###autoload (defsubst rationalp (x) "Return t if X is a rational factorial-base vector." (and (vectorp x) (= (length x) (1+ rational-max-size)) (eq (aref x 0) 'rational))) (defconst rational-zero (make-rational) "Constant zero.") (defconst rational-one (make-rational 1) "Constant one.") ;;;###autoload (defun rational-set (x int) "Set rational number X to integer INT." (let ((i rational-max-size)) (while (> i 1) (aset x i 0) (setq i (1- i))) (aset x 1 int) x)) ;;;###autoload (defun rational-zero (x) "Set rational number X to 0." (rational-set x 0)) ;;;###autoload (defun rational-one (x) "Set rational number X to 1." (rational-set x 1)) ;;;###autoload (defun rational-copy (x &optional result) "Create a copy from rational number X." (let ((copy (if (and result (not (eq x result))) result (make-rational))) (i rational-max-size)) (while (> i 0) (aset copy i (aref x i)) (setq i (1- i))) copy)) ;;;###autoload (defun integer-to-rational (p &optional result) "Convert integer P into a factorial-base number." (if result (rational-set result p) (make-rational p))) ;;;###autoload (defun rational-convert (p q &optional result) "Convert (integers) P and Q into a factorial-base number." (let ((rational (or result (make-rational))) (i 2)) (aset rational 1 (/ p q)) (and (< p 0) (/= (% p q) 0) (aset rational 1 (1- (aref rational 1)))) (setq p (% p q)) (while (<= i rational-max-size) (setq p (* p i)) (aset rational i (/ p q)) (setq p (% p q) i (1+ i))) rational)) ;;;###autoload (defun rational-lessequal (x y) "Return t if X and Y are rational and X <= Y." (let ((i (rational-differ x y))) (<= (aref x i) (aref y i)))) ;;;###autoload (defun rational-less (x y) "Return t if X and Y are rational and X < Y." (let ((i (rational-differ x y))) (< (aref x i) (aref y i)))) ;;;###autoload (defun rational-equal (x y) "Return t if X and Y are rational and X = Y." (let ((i (rational-differ x y))) (= (aref x i) (aref y i)))) ;;;###autoload (defun rational-add (x y &optional result) "Add rational numbers X and Y." (let ((temp (or result (make-rational))) (i 1)) (while (<= i rational-max-size) (aset temp i (+ (aref x i) (aref y i))) (setq i (1+ i))) (rational-smooth temp))) ;;;###autoload (defun rational-subtract (x y &optional result) "Subtract rational number X from Y." (let ((temp (or result (make-rational))) (i 1)) (while (<= i rational-max-size) (aset temp i (- (aref x i) (aref y i))) (setq i (1+ i))) (rational-smooth temp))) ;; used by `rational-absolute' (defvar rational-abs (make-rational)) ;;;###autoload (defun rational-absolute (x &optional result) "Return the absolute value of rational number X. This is trickier than flipping a bit because the sign bit is attached to x[1]. So if x[1]<0 then compute -x[1] subtract the rest of the terms x[2..n]." (cond ((< (aref x 1) 0) (rational-zero rational-abs) (aset rational-abs 1 (- (aref x 1))) (aset x 1 0) (rational-subtract rational-abs x result)) ((eq x result) x) (t (rational-copy x result)) )) ;;;###autoload (defun rational-negative (x &optional result) "Convert a positive rational number X to negative form." (let ((temp (- -1 (aref x 1)))) (aset x 1 0) (setq result (rational-subtract rational-one x result)) (aset result 1 temp) result)) ;;;###autoload (defun rational-multbyint (x int &optional result) "Multiply rational number X by an integer INT." (let ((temp (or result (make-rational))) (i 1)) (while (<= i rational-max-size) (aset temp i (* (aref x i) int)) (setq i (1+ i))) (rational-smooth temp))) ;;;###autoload (defun rational-divbyint (x int &optional result) "Divide rational number X by an integer INT." (let ((temp (or result (make-rational))) (negative (< (aref x 1) 0)) (i 1) (carry 0)) (and negative (setq x (rational-absolute x))) (while (<= i rational-max-size) (let ((part (+ (aref x i) (* carry i)))) (setq carry (% part int)) (aset temp i (/ part int))) (setq i (1+ i))) (rational-smooth temp) (if negative (rational-negative temp temp) temp))) ;; used by `rational-multiply' and `rational-divide' (defvar rational-partial (make-rational)) (defvar rational-temp (make-rational)) (defvar rational-x (make-rational)) (defvar rational-y (make-rational)) ;;;###autoload (defun rational-multiply (x y &optional result) "Multiply ratinal numbers X and Y." (let ((i 1)) (rational-copy y rational-y) (rational-zero rational-partial) (while (<= i rational-max-size) ;; shift y over one decimal place (rational-divbyint rational-y i rational-y) ;; now temp contains y*(x[i]/i!) (rational-multbyint rational-y (aref x i) rational-temp) ;; add it and continue (rational-add rational-partial rational-temp rational-partial) (setq i (1+ i))) (rational-copy rational-partial result))) ;;;###autoload (defun rational-divide (x y &optional result) "Divide ratinal numbers X by Y. Begins by scaling the numbers to find an easy, accurate way of computing the first value. After that it proceeds to use long division." (let ((negative (or (and (< (aref x 1) 0) (> (aref y 1) 0)) (and (> (aref x 1) 0) (< (aref y 1) 0)))) (i 1) (posit 1) (denomfact 1) (denom (aref y 1))) (rational-absolute x rational-x) (rational-absolute y rational-y) (while (and (< denom 100) (< i 7)) ;; get approximately 3 significant figures (setq i (1+ i) denom (+ (* denom i) (aref rational-y i)) denomfact (* denomfact i))) ;; do loop until rational-x = 0 for best accuracy (rational-zero rational-partial) (while (and (not (rational-equal rational-x rational-zero)) (<= posit rational-max-size)) (let ((estimate (/ (* (aref rational-x 1) denomfact) denom))) (while (progn (rational-multbyint rational-y estimate rational-temp) (setq estimate (1- estimate)) (rational-lessequal rational-temp rational-x))) ;; record the result (aset rational-partial posit (1+ estimate))) ;; calculate the remainder (rational-subtract rational-x rational-temp rational-x) ;; move over one notch (setq posit (1+ posit)) ;; shift the numerator over one notch (rational-multbyint rational-x posit rational-x)) (rational-smooth rational-partial) (if negative (rational-negative rational-partial result) (rational-copy rational-partial result)))) ;;;###autoload (defun rational-to-string (x) "Convert rational number X to string." (let ((i 2) str) (rational-copy x rational-x) (when (< (aref rational-x 1) 0) ;; the numbers are stored in the form sum x[i]/i! so convert to an ;; equivalent positive number by subtracting the x[2..max] terms from ;; -x[1], the integer part of the number. (rational-zero rational-y) (aset rational-y 1 (- (aref rational-x 1))) (aset rational-x 1 0) (rational-subtract rational-y rational-x rational-x) (setq str "-")) (setq str (concat str (number-to-string (aref rational-x 1)) " .")) (while (<= i rational-max-size) (setq str (concat str " " (number-to-string (aref rational-x i))) i (1+ i))) str)) ;;;###autoload (defun rational-to-float (x) "Convert rational number X to float." (let ((float (aref x 1)) (fact 1.0) (i 2) (negative (< (aref x 1) 0))) (rational-copy x rational-x) (when negative ;; the numbers are stored in the form sum x[i]/i! so convert to an ;; equivalent positive number by subtracting the x[2..max] terms from ;; -x[1], the integer part of the number. (rational-zero rational-y) (aset rational-y 1 (- (aref rational-x 1))) (aset rational-x 1 0) (rational-subtract rational-y rational-x rational-x)) (while (<= i rational-max-size) (setq fact (* fact i) float (+ float (/ (aref rational-x i) fact)) i (1+ i))) (if negative (- float) float))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Internal functions (defun rational-smooth (x) "Take rational number X and does all the carrying and the borrowing." (let ((i rational-max-size)) (while (>= i 2) (aset x (1- i) (+ (aref x (1- i)) (/ (aref x i) i))) (and (< (aref x i) 0) (/= (% (aref x i) i) 0) (aset x (1- i) (1- (aref x (1- i))))) (aset x i (% (aref x i) i)) (setq i (1- i))) x)) (defun rational-differ (x y) "Return index where X and Y differ." (let ((i 1)) (while (and (< i rational-max-size) (= (aref x i) (aref y i))) (setq i (1+ i))) i)) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Test (defun rational-test (&optional no-display) "Rational test. Consider the Pascal-like code: x := 1 / n; for i := 1 to 30 do x := (n + 1) * x - 1; Mathematically, the function f(x) = (n + 1)x - 1 is invariant at the point x = 1/n; that is, f(x) = x for x = 1/n. On paper, then, you would expect the variable x to remain unchanged after 30 interations of the loop. This routine does the above test using rational package and floating-point. Returns a string with the results." (interactive) (let ((float (/ 1.0 3.0)) (a (rational-divbyint rational-one 3)) (i 30)) (while (> i 0) (rational-multbyint a 4 a) (rational-subtract a rational-one a) (setq float (- (* 4.0 float) 1.0) i (1- i))) (let ((str (concat "\n*** Rational Test ***\n" "\nfunction: f(x) = (n + 1)x - 1 with x = 1/n and n = 3" "\n (after 30 interactions of the loop)\n" "\nUsing rational package (factorial-base): " (rational-to-string a) "\nUsing rational package (decimal-base) : " (number-to-string (rational-to-float a)) "\nUsing floating-point : " (number-to-string float) "\n\n*********************\n"))) (or no-display ;; Display the data (with-output-to-temp-buffer " *Rational Test*" (princ str) (print-help-return-message))) str))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (provide 'rational) ;;; rational.el ends here