首页 > 代码库 > 欧拉计划(python) problem 64

欧拉计划(python) problem 64

2024-11-24 22:06:01   203人阅读

Odd period square roots

Problem 64

All square roots are periodic when written as continued fractions and can be written in the form:

√N = a0 +
1
 a1 +
1
  a2 +
1
   a3 + ...

For example, let us consider √23:

√23 = 4 + √23 — 4 = 4 + 
1
 = 4 + 
1
 
1
√23—4
 1 + 
√23 – 3
7

If we continue we would get the following expansion:

√23 = 4 +
1
 1 +
1
  3 +
1
   1 +
1
    8 + ...

The process can be summarised as follows:

a0 = 4, 
1
√23—4
 = 
√23+4
7
 = 1 + 
√23—3
7
a1 = 1, 
7
√23—3
 = 
7(√23+3)
14
 = 3 + 
√23—3
2
a2 = 3, 
2
√23—3
 = 
2(√23+3)
14
 = 1 + 
√23—4
7
a3 = 1, 
7
√23—4
 = 
7(√23+4)
7
 = 8 + √23—4
a4 = 8, 
1
√23—4
 = 
√23+4
7
 = 1 + 
√23—3
7
a5 = 1, 
7
√23—3
 = 
7(√23+3)
14
 = 3 + 
√23—3
2
a6 = 3, 
2
√23—3
 = 
2(√23+3)
14
 = 1 + 
√23—4
7
a7 = 1, 
7
√23—4
 = 
7(√23+4)
7
 = 8 + √23—4

It can be seen that the sequence is repeating. For conciseness, we use the notation √23 = [4;(1,3,1,8)], to indicate that the block (1,3,1,8) repeats indefinitely.

The first ten continued fraction representations of (irrational) square roots are:

√2=[1;(2)], period=1
√3=[1;(1,2)], period=2
√5=[2;(4)], period=1
√6=[2;(2,4)], period=2
√7=[2;(1,1,1,4)], period=4
√8=[2;(1,4)], period=2
√10=[3;(6)], period=1
√11=[3;(3,6)], period=2
√12= [3;(2,6)], period=2
√13=[3;(1,1,1,1,6)], period=5

Exactly four continued fractions, for N ≤ 13, have an odd period.

How many continued fractions for N ≤ 10000 have an odd period?


Answer:
1322
Completed on Sun, 1 Feb 2015, 06:17

Go to the thread for problem 64 in the forum.


Python code :

import math
sqrt=math.sqrt


def func(x):
    index=0
    b=[0.0,0.0]
    c=[0.0]
    intcoeff={}
    fractioncoeff={}
    a=int(sqrt(x))
    if sqrt(x)-a==0:
        return 0
    b[0]=1.0
    b[1]=-a*1.0
    c=1.0
    index+=1
    intcoeff[index]=a
    tempstr=str(b[0])+‘_‘+str(b[1])+‘_‘+str(c)
    fractioncoeff[tempstr]=index
    while 1:
        result=rationalization(b,c,x)
        a=result[‘a‘]
        b=result[‘b‘]
        c=result[‘c‘]
        tempstr=str(b[0])+‘_‘+str(b[1])+‘_‘+str(c)
        k=fractioncoeff.get(tempstr)
        if k==None:
            index+=1
            intcoeff[index]=a
            fractioncoeff[tempstr]=index
        else:
            return index-k+1
# 分母包含无理数，分子是有理数
def rationalization(b,c,x):
    result={}
    tc=b[0]*b[0]*x-b[1]*b[1]
    b[0]=int(b[0]*c)
    b[1]=int(-b[1]*c)
    a=int((b[0]*sqrt(x)+b[1])/tc)
    b[1]-=int(a*tc)
    k=gcd(gcd(b[0],abs(b[1])),tc)
    b[0]/=k
    b[1]/=k
    tc/=k
    result[‘a‘]=a
    result[‘b‘]=b
    result[‘c‘]=tc
    return result
#求最大公约数
def gcd(m,n):
    k=min(m,n)
    for i in range(k,0,-1):
        if m%i==0 and n%i==0:
            return i
N=10000
count=0
for i in range(2,N+1):
    k=int(func(i))
    if k%2==1:
        count+=1
print(count)


time: 3s

欧拉计划(python) problem 64


联系
我们