I am implementing [A064604 - OEIS][1] for positive integers up to 10 billion.

I am finding the divisors in \$O(\sqrt N)\$. So, the overall time complexity of running the formula right now is \$O(N\sqrt N)\$. How do I improve on this?

[Pastebin][2]

    import math

    def factors(n):    
        return set(reduce(list.__add__, 
                    ([i, n//i] for i in range(1, int(n**0.5) + 1) if n % i == 0)))

    def sigma_4(n):
        l = factors(n)
        ans=0
        for factor in l:
            ans += (pow(factor,4))
        return ans

    n=int(raw_input())
    ans=0
    for i in xrange(1,n+1):
        ans+=sigma_4(i)
    print ans


  [1]: https://oeis.org/A064604
  [2]: http://pastebin.com/tYxRkqp2