# how accelerate this simple program

Hello everybody. Here is my question We seek integers of which the decimal part of their square root begin by 2017. Sample: sqrt(10858) = 104.201727.... Here is my program:

for n in range (1, 10^5):
if int(  N(sqrt(n)).frac() * 10000 )  == 2017:
cpt += 1


the problem is it takes too much time. For the bound 10^5 it takes about 55 sec on my iMac, for 10^6 it takes 620 sec, for 10^7 it takes almost 6000 sec ... The question is to obtain the result for the bound 10^10 !

edit retag close merge delete

Sort by » oldest newest most voted

Let us put some more notation in there, that will support the code. So we seek integers $M\le 10^{10}$ of which the decimal part $K$ of their square root begin by $2017$. Mathematically written: $$K+\frac{2017}{10\ 000} \le \sqrt M < K+\frac{2018}{10\ 000}$$

Sample: sqrt(10858) = 104.201727... . So $M=10\, 858$ is a "good $M$-value".

Proposal: Let us equivalently look for the $K$-values instead of the $M$ values. So we rewrite the statement:

we seek integers $K< 10^{10/2}$ so that between the two numbers $$a(K)=\left(K+\frac{2017}{10\ 000}\right)^2 \ ,\ b(K)=\left(K+\frac{2018}{10\ 000}\right)^2$$ there exist at least an integer $M$.

Sample: 104.2017^2 = 10857.99428289 and 104.2018^2 = 10858.01512324 . So $M=10858$ is a good $M$-value.

An other sample.

sage: "%.4f" % ( 8888.2017^2 )
'79000129.4599'
sage: "%.4f" % ( 8888.2018^2 )
'79000131.2375'


So the two numbers $M=79000129$ and $M=79000130$ are both good values. Explicitly:

sage: for K in range( 79000124, 79000135 ):    print "sqrt(%s) ~ %.4f" % ( K, sqrt(RR(K)) )
sqrt(79000124) ~ 8888.2014
sqrt(79000125) ~ 8888.2014
sqrt(79000126) ~ 8888.2015
sqrt(79000127) ~ 8888.2016
sqrt(79000128) ~ 8888.2016
sqrt(79000129) ~ 8888.2017
sqrt(79000130) ~ 8888.2017
sqrt(79000131) ~ 8888.2018
sqrt(79000132) ~ 8888.2018
sqrt(79000133) ~ 8888.2019
sqrt(79000134) ~ 8888.2020


In such a case, we have to count the integer $K$ twice.

Since $K^2\in\mathbb Z$, we can simplify the above condition. We seek equivalently integers $K< 10^{5}$ so that between the two numbers $$A(K)=2K\cdot \frac{2017}{10\ 000} +\left(\frac{2017}{10\ 000}\right)^2$$ $$B(K)=2K\cdot \frac{2018}{10\ 000} +\left(\frac{2018}{10\ 000}\right)^2$$ there exist at least an integer $M$.

The posted code seems to only count the "good" integers $M$. Then the equivalent counting program would be:

def count_M_values( bound, list_M_vales=False ):
s = 2017 / 10000
t = 2018 / 10000
ss = s^2
tt = t^2

counter = 0
mList   = []
for K in xrange( bound ):
addtocounter = floor( 2*K*t + tt ) - floor( 2*K*s + ss )
if list_M_vales:
mList . extend( [ floor( (K+s)^2 ) + 1  .. floor( (K+t)^2 ) ] )
if list_M_vales:
return counter, mList
return counter


This gives for instance:

sage: count_M_values( 10**1 )
0
sage: count_M_values( 10**2 )
0
sage: count_M_values( 200 )
3
sage: count_M_values( 200, True )
(3, [10858, 24399, 25986])
sage: count_M_values( 300, True )
(9, [10858, 24399, 25986, 45455, 47612, 49819, 73009, 75736, 78513])
sage: count_M_values( 10**3 )
99
sage: count_M_values( 10**4 )
9998
sage: count_M_values( 10**5 )
999980
sage: count_M_values( 10**6 )
99999800
sage: count_M_values( 10**7 )
9999998000


(We have a nice count pattern that we can now investigate mathematically.)

Note that the bound 10**7 above applies for $K$, we are thus going up to an $M$--bound equal to $10^{14}$. The post wanted the answer for the $K$-bound 10**5.

One more check for the $K$-bound $300$:

count, mList = count_M_values( 300, True )
for M in mList:
print "sqrt( %s ) ~ %.7f" % ( M, RR( sqrt(M) ) )


We get:

sqrt( 10858 ) ~ 104.2017274
sqrt( 24399 ) ~ 156.2017926
sqrt( 25986 ) ~ 161.2017370
sqrt( 45455 ) ~ 213.2017824
sqrt( 47612 ) ~ 218.2017415
sqrt( 49819 ) ~ 223.2017025
sqrt( 73009 ) ~ 270.2017765
sqrt( 75736 ) ~ 275.2017442
sqrt( 78513 ) ~ 280.2017131

more

Good evening. A very big thank you for your long explanation that I find extraordinary. I would never have thought of doing that. I read and reread your text several times and I think I understood the principle. The result is found in just over a second! It's awesome

( 2017-07-28 17:14:16 -0500 )edit

One bottleneck comes from working with the symbolic ring. The following code is some hundred times faster (may be still not fast enough ;-) )

import math
def count_2017(bound=100000):
cpt = 0
for n in range (1,bound):
if int(RR(math.sqrt(n)).frac()*10000) == 2017:
cpt += 1
return cpt


Update: Using cython makes the code another ten times faster, for example:

fits_2017 = cython_lambda("double x", "bool(int(10000*(x-int(x))) == 2017)")
f = cython_lambda("int n", "sqrt(n)")

def count_2017(bound=100000):
cpt = 0
for n in range(1,bound):
if fits_2017(f(n)):
cpt += 1
return cpt

count_2017(10^7)

more

Indeed your code works more than 100 times faster. I do not understand why RR (math.sqrt (n)) is 100 times faster than N (sqrt (n)). Where can I find an explanation (as simple as possible)?

I tried your cython code and I get the following error.

compilation terminated. error: command 'gcc' failed with exit status 1

Is not Sagemat already written in cython?

( 2017-07-28 16:41:36 -0500 )edit

math.sqrt returns a numerical value, sage sqrt returns a symbolic expression.

See an example: http://sagecell.sagemath.org/?q=bgdeqj

Exact symbolic operations are much slower than numerical algorithms. Converting a symbolic expression to a numerical value takes additional time.

( 2017-07-29 07:43:35 -0500 )edit

thank you for your explanations. I see that I still have a lot to learn. It's frustrating to be stuck in solving a problem because I do not know all the tricks that make a program much more effective. Where can I find answers? Which book or site or ...?

( 2017-07-30 17:15:17 -0500 )edit