How can I find the entries of the continued fraction of $2^{1/3}$ more efficiently?

Question

I wanted to detect large entries in the continued fraction of $$2^{1/3}$$

I found two ways to find the continued fraction of an irrational number like $$2^{1/3}$$

in PARI/GP. The most obvious is the contfrac-command which requires a very large precision, if I want to calculate, lets say , $10^7$ entries.

If the continued fraction expansion of $2^{1/3}$ is $[a_1,a_2,a_3,\cdots]$, then define $x_k:=[a_k,a_{k+1},a_{k+2},\cdots]$

I had the idea that I can calculate the minimal polynomials of the numbers $x_k$ (See the above definition). I wrote the program

? maxi=0;j=0;f=x^3-2;while(j<10^7,j=j+1;s=0;while(subst(f,x,s)<0,s=s+1);s=s-1;if
(s>maxi,maxi=s;print(j,"  ",s));f=subst(f,x,x+s);f=subst(f,x,1/x)*x^3;if(pollead
(f)<0,f=-f))
1  1
2  3
4  5
10  8
12  14
32  15
36  534
572  7451
1991  12737 

which works very well. The output shows the first few successive record entries in the continued fraction expansion. The program does not need floating-point calculations. It only calculates the new minimal polynomial and determines the truncated part of the (unique) positive root.

But the coefficients of the minimal polynomial get soon extremely large, hence the program needs already much longer for the first $2\cdot 10^4$ entries than the cont-frac-command.

Does anyone know a faster method with which I can calculate $10^7$ or $10^8$ entries in a reasonable time ? The contfrac-command requires much memory , so a way to avoid large arrays would be very appreciated.

Answer 1

I've computed the first 9.6 million coefficients and get the following maximum values. There is only one additional large value added to your list.

       n      a_n
       1        1
       2        3
       4        5
      10        8
      12       14
      32       15
      36      534
     572     7451
    1991    12737
   20857    22466
   27432    68346
   28763   148017
  155122   217441
  190271   320408
  288108   533679
  484709  4156269
 1395499  4886972
 9370521 10253793

Here is my Ruby code.

The big improvement over your code is that I do a binary search, first finding the smallest $n$ such that $p(2^n)<0<p(2^{n+1})$. This should reduce the number of times you have to evaluate $p$, and $p(2^n)$ can be computed faster with bit-shifts rather than multiplications. We are dealing with numbers with millions of digits. so the cost of multiplication can be high.

Once we have the values $2^n$ and $2^{n+1}$, this algorithm performs a binary search between those values to find $m$ so that $p(m)<0<p(m+1)$.

Even this could be optimized, since the binary search is pure, so we are always doing $p(N+2^k)$ for some $N$ where $p(N)$ is already known. We get $$p(N+2^k)=p(N)+p'(N)2^k + p''(N)2^{2k}/2 + p'''(N)2^{3k},$$ so we might be able to get away with speeding up those calculations. The second code example implements that. It is a bit more opaque.

# Evaluate the polynomial at a value
def evaluate(v,poly)
    total = 0
    poly.each do |c|
       total = total*v + c
    end
    total
end

# Performance hack - evaluate the polynomial at 2**n
# Uses bit-shifting to improve performance in this cases
def evaluate_bits(n,coefficients)
    total = 0
    coefficients.each do |c|
       total = (total << n) + c
    end
    total
end

# Given the current polynomial, and a value, compute
# y^3p(v+1/y)
#
# p(x)=ax^3+bx^2+cx+d 
# y^3p(v+1/y)=a(yv+1)^3 + by(yv+1)^2+cy^2(yv+1)+dy^3
#            = p(v)y^3+(3av^2+2bv+c)y^2+(3av+b)y +a
def nextPoly(v,poly)
   a,b,c,d = poly
   v2= v*v
   v3 = v2*v
   a1 = a*v3+b*v2+c*v+d
   b1 = 3*a*v2 + 2*b*v + c
   c1 = 3*a*v + b
   d1 = a
   # Take the negative, so the new polynomial is negative at 1
   [-a1,-b1,-c1,-d1]
end

# Find 2^n so that p(2^n)<0<p(2^{n+1})
def next_value_logarithmic(poly)
   bits = 0
   while evaluate_bits(bits+1,poly)<0 do
      bits = bits+1
   end
   1 << bits
end

def next_value(poly)
   #puts "#{poly}"
   # Find the smallest 2^n so that p(2^n)<0<p(2^(n+1))
   range1 = next_value_logarithmic(poly)
   range2 = range1*2
   # Do a binary search on the range 2^n.. 2^{n+1} to find value
   # v where p(v)<0<p(v+1)
   while range2-range1>1 do
       #STDERR.puts "#{range1} #{range2}"
       middle = (range1+range2)/2
       midval = evaluate(middle,poly)
       if midval < 0 then
           range1 = middle
       else
           range2 = middle
       end
   end
   range1
end

p = [1,0,0,-2]
denom1 = 1
denom2 = 0
numer1 = 0
numer2 = 1
(1..500000).each do |i|
   c = next_value(p)
   p = nextPoly(c,p)
   puts "#{i} #{c}"
end
puts '#{p}'

The more opaque script using polynomial shifting to try to avoid most multiplications.

require 'csv'

# Performance hack - evaluate a polynomial at 2**n
# Uses bit-shifting to improve performance
def evaluate_bits(n,coefficients)
    total = 0
    coefficients.each do |c|
       total = (total << n) + c
    end
    total
end

# p(x)=x^3-2
p = [1,0,0,-2]
max_index = 0
pn_1=0
qn_1=1
pn=1
qn=0
if ARGV.size>0
   # Read the CSV file for the first coefficients
   CSV.foreach(ARGV[0],col_sep: ' ') do |row|
     index = row[0].to_i
     coefficient = row[1].to_i
     max_index = index
     pn_1,pn = pn, pn*coefficient +pn_1
     qn_1,qn = qn, qn*coefficient +qn_1
     if index % 10000 == 0 
         STDERR.puts "Finished coefficient #{index}"
     end
     puts "#{index} #{coefficient}"
   end
   # p(x)=(x*p_n + p_(n-1))^3 - 2*(x*q_n+q_(n-1))^3
   a = pn*pn*pn - 2*qn*qn*qn
   b = 3*(pn*pn*pn_1 - 2*qn*qn*qn_1)
   c = 3*(pn*pn_1*pn_1 - 2*qn*qn_1*qn_1)
   d = pn_1*pn_1*pn_1 - 2*qn_1*qn_1*qn_1
   if (a>0) 
       p = [a,b,c,d]
   else
       p = [-a,-b,-c,-d]
   end
end

# Find n so that 2^n so that p(2^n)<0<p(2^{n+1})
# Returns pair n,p(2^n)
def next_value_logarithmic(poly)
   bits = 0
   while (value = evaluate_bits(bits+1,poly))<0 do
      bits = bits+1
      last_value = value
   end
   if bits == 0
     last_value = poly.inject(0){|sum,x| sum + x }    
   end
   [bits,last_value]
end

# Given cubic p(x) and n and p(2^n), compute polynomial q(x)=p(2^n+x)
# p(2^n+x)=p(2^n)+xp'(2^n)+x^2p''(2^{n})/2+x^3p'''(2^{n})/6
# No multiplication, only bit shifts and addition; 
# Total of eight operations
def shift_poly_bits(n,p_2ton,poly)
   a,b,c,d= poly
   d1 = p_2ton
   # Precompute 2^n*3*a
   a3 = (a+(a<<1))<<n
   # c1=3a*2^{2n}+2b*2^n + c
   c1 = (a3<<n) + (b<<(n+1)) + c 
   # b1 = 3a2^n + b
   b1 = a3 + b
   [a,b1,c1,d1]
end  

# Find integer m\geq 1 so that p(m)<0<p(m+1)    
def next_value(poly)
   # Find the smallest 2^n so that p(2^n)<0<p(2^(n+1))
   bits,value = next_value_logarithmic(poly)
   # At all points, shifted(x)=poly(x+minimum)
   shifted_poly = shift_poly_bits(bits,value,poly)
   minimum = 1<<bits
   # Do a binary search 
   while bits>0 do
       bits = bits-1
       midval = evaluate_bits(bits,shifted_poly)
       if midval < 0 then
           shifted_poly = shift_poly_bits(bits,midval,shifted_poly)
           minimum = minimum + (1<<bits)
       end
   end
   # returns value c wher poly(c)<0<poly(c+1)
   # and the polynomial poly(c+x)
   [minimum,shifted_poly]
end

((max_index+1)..10000000).each do |i|
   c,p2 = next_value(p)
   # p2(x) is p(c+x)
   p = p2.reverse.map {|v| -v }
   puts "#{i} #{c}"
   if i%10000== 0 then
       max = p.map { |c| c.abs() }.max
       STDERR.puts "Coefficient #{i} max bits=#{max.bit_length}"
   end
end

Closing in on n=10,000,000, the program is dealing with cubic polynomials with coefficients of about 16,500,000 bits each. In base 10, these numbers would have more than 5,000,000 digits. My laptop takes about 6.7 minutes to compute each additional 10,000 coefficients, or about 25 additional coefficients per second.

Answer 2

I'm not sure if you are still interested in this problem. Continued fractions of roots of integers have a very regular format, but unfortunately only as a general continued fraction and anyway the terms grow without bound. Below, $\alpha$ is the nth integer root of $x$ and $\beta$ is the difference between x and $\alpha^n$. This continued fraction I found in "General Methods for Extracting Roots" by Manny Sardina.

$$\sqrt[n]{x^m} = (\alpha^n + \beta)^{\frac{m}{n}} = \alpha^m + \cfrac{\beta m}{n\alpha^{n-m} + \cfrac{\beta (n-m)}{2\alpha^m + \cfrac{\beta (n+m)}{3n\alpha^{n-m} + \cfrac{\beta (2n-m)}{2\alpha^m + \cfrac{\beta (2n+m)}{5n\alpha^{n-m} + \cdots}}}}}$$

Converting this to matrix multiplication is straightforward.

$$\begin{pmatrix} \alpha^m & \beta \cdot m \\ 1 & 0 \end{pmatrix} \times \begin{pmatrix} n\cdot \alpha^{n-m} & \beta \cdot (n - m) \\ 1 & 0 \end{pmatrix} \times \cdots$$

Since matrix multiplication is associative and the terms of this continued fraction have such a regular form you can program a slightly more efficient representation by doing pairs of multiplications up front.

At any point from left to right as you gobble up terms you will have some state $\begin{pmatrix} a & b \\ c & d \end{pmatrix}$ and you can emit a term of a simple continued fraction $t$ when $\lfloor \frac{a}{c} \rfloor = \lfloor \frac{b}{d} \rfloor \rightarrow t$ whereupon the state matrix becomes $\begin{pmatrix} c & d \\ a-tc & b-td \end{pmatrix}$.

You will observe that in general the state matrix terms will grow without bound for all irrationals, but some peace may be had since these terms will definitely not grow as fast as polynomial coefficients in Vincent's Theorem. I don't know what PARI/GP uses to calculate these continued fractions but as a guess they would experience some similar kind of overhead in any case.

Answer 3

The continued fraction of $\sqrt[3]{2}$ are conjectured to never repeat themeselves (to not be periodic). That's quite a statement. This is discussed in the paper by Lang and Trotter Continued Fractions of Some Algebraic Numbers. This is also sequence A002945 in the OEIS.

Here's the first 100 continued fraction digits: