continued fractions for nth degree roots
- Subject: continued fractions for nth degree roots
- From: Robert Dodier
- Date: Tue, 6 Feb 2007 08:44:25 -0700
Hello Andrei,
I've committed your code for continued fractions of nth degree roots,
from your message to the mailing list dated 2005/05/07.
E.g.
cflength : 8;
cf (2^(1/3));
=> [1, 3, 1, 5, 1, 1, 5]
float (cfdisrep (cf (2^(1/3))));
=> 1.259927797833935
float (2^(1/3));
=> 1.259921049894873
Thanks a lot for your help! The wheels of Maxima grind slowly,
but exceedingly fine.
best,
Robert Dodier
PS.
--- src/combin.lisp 3 Jan 2007 09:03:08 -0000 1.14
+++ src/combin.lisp 3 Feb 2007 04:17:10 -0000
@@ -591,6 +591,17 @@
(t (cons '(mlist cf simp)
(apply 'find-cf (cf-back-recurrence (cdr a)))))))
+; Code to expand nth degree roots of integers into continued fraction
+; approximations. E.g. cf(2^(1/3))
+; Courtesy of Andrei Zorine (feniy at mail.nnov.ru) 2005/05/07
+
+(defun cfnroot(b)
+ (let ((ans (list '(mlist xf))) ent ($algebraic $true))
+ (dotimes (i $cflength (nreverse ans))
+ (setq ent (meval `(($floor) ,b))
+ ans (cons ent ans)
+ b ($ratsimp (m// (m- b ent)))))))
+
(defun cfeval (a)
(let (temp $ratprint)
(cond ((integerp a) (list '(mlist cf) a))
@@ -622,6 +633,7 @@
(cftimes (cfsqrt (cfeval (cadr a)))
(cfexpt (cfeval (cadr a))
(m- (caddr a) '((rat) 1 2)))))
+ ((integerp (cadr a)) (cfnroot a)) ; <=== new case x
((cfexpt (cfeval (cadr a)) (caddr a)))))
((setq temp (assq (caar a) '((mplus . cfplus) (mtimes .
cftimes) (mquotient . cfquot)
(mdifference . cfdiff) (mminus
. cfminus))))