continued fractions for nth degree roots

Here's the code for continued fraction expansion of *arbitrary*
numerical-valued expressions, including transcendental functions. It is
pretty robust even for rational results, e.g. sin(2)^2+cos(2)^2 or
asin(1/sqrt(2))/%pi, and gives a warning when the result might look more
precise than it is.

Areas to improve:

   - It is *not* robust at complex singularities, so cfb(
   asin(sin(2)^2+cos(2)^2)/%pi , 5) gives an error.
   - It could do better at noticing how far it needs to expand (though
   the user can always increase cfb_maxcount).

Feedback is welcome.

                 -s

--------------------------------------------------------------

cfb_maxcount: 6$

cfb(ex,terms) :=
     block(
      [res,                     /* current CF */
       trunc,                   /* truncated CF */
       fpprec:16,               /* increasing precision */
       count:0],                /* iteration count */
      if terms<1 then error("cfb requires a positive number of terms"),
      if ratnump(ex) then res:cf(ex) else (
        if floatnump(ex) then ex: bfloat(ex),
           /* just convert floats to bfloat once */
        while ( length( res: cf(bfloat(ex)) ) )  <  terms + 3
                 and
              (count: count+1) < cfb_maxcount
          do fpprec:fpprec*2,
        if count = cfb_maxcount then
          (print("Warning: cfb may be off by 10^ -", fpprec),
           /* discard last term? */
           if ( bfloat(
               cfdisrep(res) -
               (if length(res)<=1 then (trunc:[0],0)
                else cfdisrep(trunc: reverse(rest(reverse(res))))))
              / 10^(1-fpprec) )
            < 1
           then res:trunc)),
      /* truncate to requested number of terms + 1
         then normalize final term */
      res: rest( reverse(res), max(0,length(res)-terms-1 )),
      reverse( if length(res) < 2 then res
                elseif res[1]=1 then cons(res[2]+1,rest(res,2))
                elseif count=cfb_maxcount then res
                else rest(res,1)
                ))$

Test cases

cfb(.25,5);
Warning:  Float to bigfloat conversion of 0.25
Warning: cfb may be off by 10^ - 512
  =>  [0, 4]

cfb(1/4,5);
  =>  [0, 4]

cfb(sqrt(7),3);
  =>  [2, 1, 2]

cfb(sqrt(7),100);
  =>  [2, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1,
1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1,
4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4,
1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1,
1, 1]

cfb(sin(%pi/10),10);
  =>  [0, 3, 4, 4, 4, 4, 4, 4, 4, 4]

cfb(sin(1)^2+cos(1)^2,20);
Warning: cfb may be off by 10^ - 512
  =>  [1]

cfb(sin(1)^2+cos(1)^2+10^-500,20);
Warning: cfb may be off by 10^ - 512
  =>  [1,
99999999999974775246987723107459958889016439203282138228671991035264810602432912576727037490136046301332696648835698988844793279690261282050403492720817487566483059823998582824631757396835957729325333565882924021109558696660456072021234688669172916757132142477460719221645589937039887851829233616435452465975791561135280819321112028808791096875355800150184829419265834651862286504605407797719239296022467521797487303629823365392127805180388589756624999779798816755036489606075204618469638428130129716]

cfb(sin(1)^2+cos(1)^2+10^-600,20);
Warning: cfb may be off by 10^ - 512
  =>  [1]

cfb(asin(sin(2)^2+cos(2)^2)/%pi,20);
%i - not a continued fraction
#0: cfb(ex=asin(sin(2)^2+cos(2)^2)/%pi,terms=20)
 -- an error.  To debug this try debugmode(true);

cfb((17/23)^(9/10),500);
  =>  [0, 1, 3, 5, 24, 1, 4, 5, 1, 35, 2, 1, 1, 3, 11, 2, 1, 1, 2, 1, 1, 2,
1, 2, 2, 1, 2, 2, 2, 1, 11, 1, 68, 1, 1, 2, 3, 2, 1, 2, 1, 1, 3, 5, 1, 1, 2,
59, 11, 1, 3, 1, 6, 3, 1, 23, 2, 17, 2, 1, 1, 1, 4, 4, 1, 4, 39, 3, 22, 1,
2, 3, 1, 2, 7, 1, 2, 2, 1, 2, 2, 2, 1, 4, 1, 3, 4, 3, 3, 1, 1, 2, 1, 8, 1,
4, 17, 1, 7, 1, 4, 1, 1, 1, 1, 12, 1, 11, 1, 1, 1, 7, 1, 1, 73, 1, 1, 1, 1,
3, 2, 1, 1, 1, 7, 1, 19, 4, 2, 1, 7, 1, 2, 3, 15, 1, 2, 1, 2, 1, 1, 1, 3, 2,
1, 8, 1, 4, 25, 3, 2, 1, 6, 1, 2, 1, 15, 2, 1, 2, 1, 130, 1, 22, 5, 7, 1, 4,
1, 1, 1, 7, 17, 1, 23, 1, 1, 1, 17, 11, 1, 3, 1, 1, 1, 20, 2, 3, 26, 4, 113,
2, 3, 1, 3, 2, 3, 111, 5, 4, 22, 5, 3, 2, 4, 7, 1, 1, 9, 1, 1, 67, 3, 2, 1,
1, 11, 1, 1, 2, 3, 21, 1, 1, 1, 1, 9, 3, 8, 2, 1, 1, 1, 1, 6, 6, 1, 1, 190,
4, 3, 1, 1, 1, 2, 46, 1, 16, 1, 3, 4, 1, 1, 1, 5, 8, 1, 1, 3, 14, 3, 1, 33,
7, 6, 5, 1, 1, 2, 1, 11, 47, 1, 4, 5, 2, 3, 1, 3, 1, 4, 1, 4, 1, 1, 2, 1, 1,
1, 4, 5, 1, 1, 4, 21, 1, 5, 2, 100, 1, 1, 1, 3, 1, 96, 2, 5, 1, 1, 1, 2, 3,
1, 4, 1, 1, 3, 1, 1, 16, 36, 3, 1, 3, 8, 1, 1, 2, 1, 6, 8, 1, 1, 3, 1, 2, 7,
1, 2, 2, 13, 2, 1, 19, 15, 1, 5, 3, 1, 9, 5, 1, 7, 1, 1, 1, 3, 37, 1, 2, 1,
1, 1, 2, 1, 1, 1, 3, 1, 2, 1, 1, 2, 6, 1, 7, 1, 2, 1, 1, 1, 4, 11, 168, 8,
7, 3, 16, 1, 1, 6, 1, 21, 4, 8, 1, 239, 3, 1, 2, 2, 3, 2, 3, 3, 1, 8, 2, 1,
1, 3, 1, 5, 3, 1, 2, 1, 52, 1, 5, 2, 2, 1, 3, 4, 1, 5, 2, 1, 13, 1, 23, 10,
1, 3, 343, 4, 1, 1, 17, 2, 17, 5, 5, 2, 8, 11, 1, 2, 2, 1, 4, 1, 35, 2, 1,
2, 3, 8, 1, 20, 1, 2, 1, 1, 1, 2, 2, 2, 4, 1, 1, 6, 1, 1, 1, 1, 9, 1, 28, 2,
7, 3, 2, 6, 1, 4, 2, 1, 9, 1, 28, 1, 1, 2, 1, 1, 5, 2, 8]

cfb(%e^(2/1001),50);
  =>  [1, 500, 6006, 2502, 1, 1, 3503, 18018, 5505, 1, 1, 6506, 30030, 8508,
1, 1, 9509, 42042, 11511, 1, 1, 12512, 54054, 14514, 1, 1, 15515, 66066,
17517, 1, 1, 18518, 78078, 20520, 1, 1, 21521, 90090, 23523, 1, 1, 24524,
102102, 26526, 1, 1, 27527, 114114, 29529, 2]

cfb((exp(2/5)-1)/(exp(2/5)+1),90)
  =>  [0, 5, 15, 25, 35, 45, 55, 65, 75, 85, 95, 105, 115, 125, 135, 145,
155, 165, 175, 185, 195, 205, 215, 225, 235, 245, 255, 265, 275, 285, 295,
305, 315, 325, 335, 345, 355, 365, 375, 385, 395, 405, 415, 425, 435, 445,
455, 465, 475, 485, 495, 505, 515, 525, 535, 545, 555, 565, 575, 585, 595,
605, 615, 625, 635, 645, 655, 665, 675, 685, 695, 705, 715, 725, 735, 745,
755, 765, 775, 785, 795, 805, 815, 825, 835, 845, 855, 865, 875, 885]

cfb(%e^(2/1001),100);
Warning: cfb may be off by 10^ - 512
  =>  [1, 500, 6006, 2502, 1, 1, 3503, 18018, 5505, 1, 1, 6506, 30030, 8508,
1, 1, 9509, 42042, 11511, 1, 1, 12512, 54054, 14514, 1, 1, 15515, 66066,
17517, 1, 1, 18518, 78078, 20520, 1, 1, 21521, 90090, 23523, 1, 1, 24524,
102102, 26526, 1, 1, 27527, 114114, 29529, 1, 1, 30530, 126126, 32532, 1, 1,
33533, 138138, 35535, 1, 1, 36536, 150150, 38538, 1, 1, 39539, 162162,
41541, 1, 1, 42542, 174174, 44544, 1, 1, 45545, 186186, 47547, 1, 1, 48548,
198198, 50550, 1, 1, 51551, 210210, 53553, 1, 1, 54554]

cfb((exp(2/100001)-1)/(exp(2/100001)+1),40);
Warning: cfb may be off by 10^ - 512
  =>  [0, 100001, 300003, 500005, 700007, 900009, 1100011, 1300013, 1500015,
1700017, 1900019, 2100021, 2300023, 2500025, 2700027, 2900029, 3100031,
3300033, 3500035, 3700037, 3900039, 4100041, 4300043, 4500045, 4700047,
4900049, 5100051, 5300053, 5500055, 5700057, 5900059, 6100061, 6300063,
6500065, 6700067, 6900069, 7100071, 7300073, 7500075, 7700077]