<p>Offhand I don't see a slick way to make this truly efficient, but it's easy to make it a whole lot faster.  If you view your examples as matrices, you're summing them a row at a time.  This requires, for each <code>i</code>, finding all the divisors of <code>i</code> and summing their cubes.  In all, this requires a number of operations proportional to <code>x**2</code>.</p>

<p>You can easily cut that to a number of operations proportional to <code>x</code>, by summing the matrix by <em>columns</em> instead.  Given an integer <code>j</code>, how many integers in <code>1..x</code> are divisible <em>by</em> <code>j</code>?  That's easy:  there are <code>x//j</code> multiples of <code>j</code> in the range, so divisor <code>j</code> contributes <code>j**3 * (x // j)</code> to the grand total.</p>

<pre><code>def better(x):
    return sum(j**3 * (x // j) for j in range(1, x+1))
</code></pre>

<p>That runs much faster, but still takes time proportional to <code>x</code>.</p>

<p>There are lower-level tricks you can play to speed that in turn by constant factors, but they still take <code>O(x)</code> time overall.  For example, note that <code>x // j == 1</code> for all <code>j</code> such that <code>x // 2 &lt; j &lt;= x</code>.  So about half the terms in the sum can be skipped, replaced by closed-form expressions for a sum of consecutive cubes:</p>

<pre><code>def sum3(x):
    """Return sum(i**3 for i in range(1, x+1))"""
    return (x * (x+1) // 2)**2

def better2(x):
    result = sum(j**3 * (x // j) for j in range(1, x//2 + 1))
    result += sum3(x) - sum3(x//2)
    return result
</code></pre>

<p><code>better2()</code> is about twice as fast as <code>better()</code>, but to get faster than <code>O(x)</code> would require deeper insight.</p>

<h2>Quicker</h2>

<p>Thinking about this in spare moments, I still don't have a truly clever idea.  But the last idea I gave can be carried to a logical conclusion:  don't just group together divisors with only one multiple in range, but also those with two multiples in range, and three, and four, and ... That leads to <code>better3()</code> below, which does a number of operations roughly proportional to the square root of <code>x</code>:</p>

<pre><code>def better3(x):
    result = 0
    for i in range(1, x+1):
        q1 = x // i
        # value i has q1 multiples in range
        result += i**3 * q1
        # which values have i multiples?
        q2 = x // (i+1) + 1
        assert x // q1 == i == x // q2
        if i &lt; q2:
            result += i * (sum3(q1) - sum3(q2 - 1))
        if i+1 &gt;= q2: # this becomes true when i reaches roughly sqrt(x)
            break
    return result
</code></pre>

<p>Of course <code>O(sqrt(x))</code> is an enormous improvement over the original <code>O(x**2)</code>, but for very large arguments it's still impractical.  For example <code>better3(10**6)</code> appears to complete instantly, but <code>better3(10**12)</code> takes a few seconds, and <code>better3(10**16)</code> is time for a coffee break ;-)</p>

<p>Note:  I'm using Python 3.  If you're using Python 2, use <code>xrange()</code> instead of <code>range()</code>.</p>

<h2>One more</h2>

<p><code>better4()</code> has the same <code>O(sqrt(x))</code> time behavior as <code>better3()</code>, but does the summations in a different order that allows for simpler code and fewer calls to <code>sum3()</code>.  For "large" arguments, it's about 50% faster than <code>better3()</code> on my box.</p>

<pre><code>def better4(x):
    result = 0
    for i in range(1, x+1):
        d = x // i
        if d &gt;= i:
            # d is the largest divisor that appears `i` times, and
            # all divisors less than `d` also appear at least that
            # often.  Account for one occurence of each.
            result += sum3(d)
        else:
            i -= 1
            lastd = x // i
            # We already accounted for i occurrences of all divisors
            # &lt; lastd, and all occurrences of divisors &gt;= lastd.
            # Account for the rest.
            result += sum(j**3 * (x // j - i)
                          for j in range(1, lastd))
            break
    return result
</code></pre>

<p>It may be possible to do better by extending the algorithm in <a href="https://arxiv.org/abs/1206.3369" rel="nofollow noreferrer">"A Successive Approximation Algorithm for Computing the Divisor Summatory Function"</a>.  That takes <code>O(cube_root(x))</code> time for the possibly simpler problem of summing the number of divisors.  But it's much more involved, and I don't care enough about this problem to pursue it myself ;-)</p>

<h2>Subtlety</h2>

<p>There's a subtlety in the math that's easy to miss, so I'll spell it out, but only as it pertains to <code>better4()</code>.</p>

<p>After <code>d = x // i</code>, the comment claims that <code>d</code> is the largest divisor that appears <code>i</code> times.  But is that true?  The actual number of times <code>d</code> appears is <code>x // d</code>, which we did <em>not</em> compute.  How do we know that <code>x // d</code> in fact equals <code>i</code>?</p>

<p>That's the purpose of the <code>if d &gt;= i:</code> guarding that comment.  After <code>d = x // i</code> we know that</p>

<pre><code>x == d*i + r
</code></pre>

<p>for some integer <code>r</code> satisfying <code>0 &lt;= r &lt; i</code>.  That's essentially what floor division <em>means</em>.  But since <code>d &gt;= i</code> is also known (that's what the <code>if</code> test ensures), it must also be the case that <code>0 &lt;= r &lt; d</code>.  And that's how we know <code>x // d</code> is <code>i</code>.</p>

<p>This can break down when <code>d &gt;= i</code> is <em>not</em> true, which is why a different method needs to be used then.  For example, if <code>x == 500</code> and <code>i == 51</code>, <code>d</code> (<code>x // i</code>) is 9, but it's certainly not the case that 9 is the largest divisor that appears 51 times.  In fact, 9 appears <code>500 // 9 == 55</code> times.  While for positive real numbers</p>

<pre><code>d == x/i
</code></pre>

<p>if and only if</p>

<pre><code>i == x/d
</code></pre>

<p>that's <em>not</em> always so for floor division.  But, as above, the first does imply the second if we also know that <code>d &gt;= i</code>.</p>

<h2>Just for Fun</h2>

<p><code>better5()</code> rewrites <code>better4()</code> for about another 10% speed gain.  The real pedagogical point is to show that it's easy to compute all the loop limits in advance.  Part of the point of the odd code structure above is that it magically returns 0 for a 0 input without needing to test for that.  <code>better5()</code> gives up on that:</p>

<pre><code>def isqrt(n):
    "Return floor(sqrt(n)) for int n &gt; 0."
    g = 1 &lt;&lt; ((n.bit_length() + 1) &gt;&gt; 1)
    d = n // g
    while d &lt; g:
        g = (d + g) &gt;&gt; 1
        d = n // g
    return g

def better5(x):
    assert x &gt; 0
    u = isqrt(x)
    v = x // u
    return (sum(map(sum3, (x // d for d in range(1, u+1)))) +
            sum(x // i * i**3 for i in range(1, v)) -
            u * sum3(v-1))
</code></pre>
How to find sum of cubes of the divisors for every number from 1 to input number x in python where x can be very large
<pre><code>Examples, 
1.Input=4
Output=111

Explanation,
1 = 1³(divisors of 1)
2 = 1³ + 2³(divisors of 2)
3 = 1³ + 3³(divisors of 3)
4 = 1³ + 2³ + 4³(divisors of 4)
------------------------
sum = 111(output)

1.Input=5
Output=237

Explanation,
1 = 1³(divisors of 1)
2 = 1³ + 2³(divisors of 2)
3 = 1³ + 3³(divisors of 3)
4 = 1³ + 2³ + 4³(divisors of 4)
5 = 1³ + 5³(divisors of 5)
-----------------------------
sum = 237 (output)

x=int(raw_input().strip())
tot=0
for i in range(1,x+1):
    for j in range(1,i+1):
        if(i%j==0):
            tot+=j**3
print tot
</code></pre>

<p>Using this code I can find the answer for small number less than one million. 
   But I want to find the answer for very large numbers. Is there any algorithm 
   for how to solve it easily for large numbers?</p>