Revision history [back]

The best performance I was able to archive is by the following code, which avoids recursion and uses a few tricks to save memory:

def divbelow_10(n,B):
    if B*B >= n:
        return number_of_divisors(n) - divbelow_10(n,n//(B+1))

    res = 0

    bs = {B:1}   # dict mapping bounds to their multiplicities

    for p,e in reversed(factor(n)):

        n //= p^e
        nd = number_of_divisors(n)
        n2 = 0                       # number of bounds >= n

        ext = bs
        for i in range(e):
            ext2 = dict()
            for b,c in ext.items():
                r = b//p
                if r >= 1:
                    ext2[r] = ext2.get(r,0) + c
            ext = ext2

            # extending bs with ext
            for b,c in ext.items():
                if b>=n:
                    n2 += c
                else:
                    bs.setdefault(b,0)
                    bs[b] += c

        #print('Extended:', len(bs))


        for b in [b for b in bs if b*b>=n]:
            c = bs.pop(b)
            n2 += c
            b_ = n//(b+1)
            if b_ >= 1:
                bs.setdefault(b_,0)
                bs[b_] -= c

        res += bs.pop(1,0)
        res += n2*nd

        #print('Trimmed:', len(bs))

    return res

Its performance is very good:

sage: %time divbelow_10(47987366728745697656468743117440000, 61448226992991993)
CPU times: user 67.5 ms, sys: 75 µs, total: 67.6 ms
Wall time: 66.9 ms
10074213

The best performance I was able to archive is by the following code, which avoids recursion and uses a few tricks to save memory:

def divbelow_10(n,B):
    if B*B >= n:
        return number_of_divisors(n) - divbelow_10(n,n//(B+1))

    res = 0

    bs = {B:1}   # dict mapping bounds to their multiplicities

    for p,e in reversed(factor(n)):

        n //= p^e
        nd = number_of_divisors(n)
        n2 = 0                       # number of bounds >= n

        ext = bs
        for i in range(e):
            ext2 = dict()
            for b,c in ext.items():
                r = b//p
                if r >= 1:
                    ext2[r] = ext2.get(r,0) + c
            ext = ext2

            # extending bs with ext
            for b,c in ext.items():
                if b>=n:
                    n2 += c
                else:
                    bs.setdefault(b,0)
                    bs[b] += c

        #print('Extended:', len(bs))


        for b in [b for b in bs if b*b>=n]:
            c = bs.pop(b)
            n2 += c
            b_ = n//(b+1)
            if b_ >= 1:
                bs.setdefault(b_,0)
                bs[b_] -= c

        res += bs.pop(1,0)
        res += n2*nd
n2 * number_of_divisors(n)

        #print('Trimmed:', len(bs))

    return res

Its performance is very good:

sage: %time divbelow_10(47987366728745697656468743117440000, 61448226992991993)
CPU times: user 67.5 ms, sys: 75 µs, total: 67.6 ms
Wall time: 66.9 ms
10074213