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?
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
原文：http://blog.csdn.net/zhangzhengyi03539/article/details/43371483