Ask Your Question
0

how accelerate this simple program

asked 2017-07-27 17:37:44 -0500

wisher gravatar image

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 flag offensive close merge delete

2 answers

Sort by ยป oldest newest most voted
1

answered 2017-07-27 21:36:30 -0500

dan_fulea gravatar image

updated 2017-07-27 21:43:38 -0500

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 addtocounter:
            counter += addtocounter
            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
edit flag offensive delete link more

Comments

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

wisher gravatar imagewisher ( 2017-07-28 17:14:16 -0500 )edit
1

answered 2017-07-28 04:08:50 -0500

ndomes gravatar image

updated 2017-07-28 11:49:32 -0500

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)
edit flag offensive delete link more

Comments

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?

wisher gravatar imagewisher ( 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.

ndomes gravatar imagendomes ( 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 ...?

wisher gravatar imagewisher ( 2017-07-30 17:15:17 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2017-07-27 17:37:44 -0500

Seen: 57 times

Last updated: Jul 28