continued fractions for nth degree roots



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