6
$\begingroup$

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.

$\endgroup$
18
  • 1
    $\begingroup$ So far, the largest coefficient I've found is 4,156,269, at position 484709 (counting from 1.) $\endgroup$ Feb 6, 2017 at 17:12
  • 2
    $\begingroup$ I managed to compute 5769963 coefficients without finding a larger one than 4,886,972. At this point I was doing computations with numbers that had close to 10,000,000 bits. $\endgroup$ Feb 8, 2017 at 14:29
  • 1
    $\begingroup$ Yes, it is free and open source. It comes with most Macs, and is available for Linux and Windows, as well. I've never used it on Windows, personally. It isn't really known for its numeric packages - it is not great at crunching big data sets - but it does have arbitrary-precision integers that perform reasonably well. $\endgroup$ Feb 8, 2017 at 14:58
  • 1
    $\begingroup$ I've gotten 8.5 million coefficients, still with none bigger than 4,886,972. At this point, there were coefficients of the polynomial with 13 million bits, and it was taking about 6 minutes to compute each 10,000 coefficients additionally. I have a file with all the coefficients to that point. $\endgroup$ Feb 12, 2017 at 15:56
  • 2
    $\begingroup$ I've found a new record: the 9,370,521th coefficient is 10,253,793 $\endgroup$ Feb 13, 2017 at 14:27

3 Answers 3

3
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ This is essentially what I have already done. $\endgroup$
    – Peter
    Feb 6, 2017 at 14:31
  • $\begingroup$ All right, real Ruby code present. $\endgroup$ Feb 6, 2017 at 18:41
2
$\begingroup$

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.

$\endgroup$
2
  • 2
    $\begingroup$ This is 'a' continued fraction, but not 'the' continued fraction. I believe OP's interest is very specifically in the continued fraction itself (and specifically in properties of its coefficients, such as how well they fit the Khinchin distribution), and not how to evaluate $\sqrt[3]{2}$ using it. $\endgroup$ Apr 14, 2017 at 0:53
  • $\begingroup$ Hi @stevenstadnicki, I understood the question to be how to get at the terms of the simple continued fraction in some more efficient manner than he was doing. I know of no direct method, only by consuming terms of the general continued fraction into a state matrix and then extracting the terms of the simple continued fraction from this state matrix (the $t$ terms as above). $\endgroup$ Apr 14, 2017 at 1:10
-1
$\begingroup$

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:

enter image description here

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .