from math import isqrt
PHI = [0, 1, 2, 4, 6, 10, 12, 18, 22, 28, 32, 42, 46, 58, 64, 72, 80, 96, 102, 120, 128, 140, 150, 172, 180, 200, 212, 230, 242, 270, 278, 308, 324, 344, 360, 384, 396, 432, 450, 474, 490, 530, 542, 584, 604, 628, 650, 696, 712, 754, 774, 806, 830, 882, 900, 940, 964, 1000, 1028, 1086, 1102, 1162, 1192, 1228, 1260, 1308, 1328, 1394, 1426, 1470, 1494, 1564, 1588, 1660, 1696, 1736, 1772, 1832, 1856, 1934, 1966, 2020, 2060, 2142, 2166, 2230, 2272, 2328, 2368, 2456, 2480, 2552, 2596, 2656, 2702, 2774, 2806, 2902, 2944, 3004, 3044, 3144, 3176, 3278, 3326, 3374, 3426, 3532, 3568, 3676, 3716, 3788, 3836, 3948, 3984, 4072, 4128, 4200, 4258, 4354, 4386, 4496, 4556, 4636, 4696, 4796, 4832, 4958, 5022, 5106, 5154, 5284, 5324, 5432, 5498, 5570, 5634, 5770, 5814, 5952, 6000, 6092, 6162, 6282, 6330, 6442, 6514, 6598, 6670, 6818, 6858, 7008, 7080, 7176, 7236, 7356, 7404, 7560, 7638, 7742, 7806, 7938, 7992, 8154, 8234, 8314, 8396, 8562, 8610, 8766, 8830, 8938, 9022, 9194, 9250, 9370, 9450, 9566, 9654, 9832, 9880, 10060, 10132, 10252, 10340, 10484, 10544, 10704, 10796, 10904, 10976, 11166, 11230, 11422, 11518, 11614, 11698, 11894, 11954, 12152, 12232, 12364, 12464, 12632, 12696, 12856, 12958, 13090, 13186, 13366, 13414, 13624, 13728, 13868, 13974, 14142, 14214, 14394, 14502, 14646, 14726, 14918, 14990, 15212, 15308, 15428, 15540, 15766, 15838, 16066, 16154, 16274, 16386, 16618, 16690, 16874, 16990, 17146, 17242, 17480, 17544, 17784, 17894, 18056, 18176, 18344, 18424, 18640, 18760, 18924, 19024, 19274, 19346, 19566, 19692, 19820, 19948, 20204, 20288, 20504, 20600, 20768, 20898, 21160, 21240, 21448, 21556, 21732, 21864, 22132, 22204, 22474, 22602, 22746, 22882, 23082, 23170, 23446, 23584, 23764, 23860, 24140, 24232, 24514, 24654, 24798, 24918, 25158, 25254, 25526, 25638, 25830, 25974, 26266, 26350, 26582, 26726, 26906, 27054, 27318, 27398, 27650, 27800, 28000, 28144, 28384, 28480, 28786, 28906, 29110, 29230, 29540, 29636, 29948, 30104, 30248, 30404, 30720, 30824, 31104, 31232, 31444, 31576, 31864, 31972, 32212, 32374, 32590, 32750, 33026, 33106, 33436, 33600, 33816, 33982, 34246, 34342, 34678, 34834, 35058, 35186, 35486, 35594, 35888, 36056, 36232, 36404, 36750, 36862, 37210, 37330, 37546, 37706, 38058, 38174, 38454, 38630, 38822, 39000, 39358, 39454, 39796, 39976, 40196, 40340, 40628, 40748, 41114, 41290, 41530, 41674, 41986, 42106, 42478, 42638, 42838, 43022, 43358, 43466, 43844, 43988, 44240, 44430, 44812, 44940, 45180, 45372, 45624, 45816, 46204, 46300, 46652, 46820, 47080, 47276, 47588, 47708, 48104, 48302, 48518, 48678, 49078, 49210, 49570, 49770, 49986, 50154, 50514, 50642, 51050, 51210, 51482, 51686, 52034, 52166, 52494, 52686, 52962, 53142, 53560, 53656, 54076, 54286, 54562, 54770, 55090, 55230, 55590, 55802, 56042, 56210, 56640, 56784, 57216, 57396, 57620, 57836, 58232, 58376, 58814, 58974, 59226, 59418, 59860, 60004, 60356, 60578, 60874, 61066, 61514, 61634, 62034, 62258, 62558, 62784, 63072, 63216, 63672, 63900, 64188, 64364, 64824, 64944, 65406, 65630, 65870, 66102, 66568, 66712, 67108, 67292, 67604, 67836, 68256, 68412, 68772, 68964, 69276, 69514, 69992, 70120, 70552, 70792, 71056, 71276, 71660, 71822, 72308, 72548, 72872, 73040, 73530, 73690, 74138, 74354, 74594, 74834, 75254, 75418, 75916, 76116, 76448, 76698, 77200, 77344, 77744, 77964, 78276, 78528, 79036, 79164, 79596, 79852, 80176, 80432, 80840, 81008, 81468, 81684, 82028, 82220, 82740, 82908, 83430, 83690, 83930, 84192, 84672, 84832, 85338, 85546, 85894, 86110, 86590, 86766, 87190, 87454, 87810, 88078, 88498, 88642, 89182, 89452, 89812, 90068, 90500, 90644, 91190, 91462, 91822, 92022, 92526, 92702, 93170, 93446, 93734, 94010, 94566, 94746, 95250, 95442, 95762, 96042, 96604, 96788, 97236, 97518, 97842, 98122, 98690, 98834, 99404, 99644, 100024, 100264, 100704, 100896, 101472, 101744, 102128, 102352, 102844, 103036, 103556, 103844, 104132, 104424, 105010, 105178, 105718, 105950, 106342, 106630, 107222, 107402, 107786, 108082, 108478, 108742, 109340, 109500, 110100, 110352, 110748, 111048, 111488, 111688, 112294, 112582, 112918, 113158, 113710, 113902, 114514, 114820, 115140, 115380, 115996, 116200, 116818, 117058, 117454, 117764, 118292, 118484, 118984, 119296, 119656, 119968, 120544, 120688, 121318, 121630, 122050, 122366, 122870, 123078, 123582, 123862, 124282, 124538, 125178, 125390, 126032, 126296, 126632, 126920, 127566, 127782, 128362, 128602, 128962, 129286, 129938, 130154, 130674, 130994, 131426, 131702, 132360, 132520, 133180, 133510, 133894, 134222, 134654, 134870, 135486, 135818, 136262, 136526, 137126, 137318, 137990, 138326, 138686, 138998, 139674, 139898, 140474, 140730, 141182, 141482, 142164, 142380, 142924, 143218, 143674, 144010, 144634, 144810, 145500, 145844, 146204, 146550, 147102, 147326, 147966, 148314, 148778, 149018, 149718, 149934, 150582, 150902, 151270, 151622, 152222, 152454, 153162, 153442, 153910, 154262, 154922, 155114, 155594, 155950, 156426, 156784, 157502, 157694, 158306, 158648, 159128, 159488, 160048, 160268, 160994, 161282, 161768, 162056, 162728, 162968, 163700, 164066, 164402, 164754, 165414, 165654, 166392, 166680, 167112, 167424, 168166, 168406, 168998, 169370, 169862, 170182, 170818, 171018, 171768, 172136, 172636, 172972, 173572, 173788, 174544, 174922, 175362, 175650, 176410, 176662, 177310, 177690, 178074, 178456, 179152, 179408, 180176, 180416, 180928, 181312, 182084, 182336, 182936, 183320, 183752, 184140, 184860, 185052, 185752, 186104, 186608, 186944, 187568, 187828, 188614, 189006, 189530, 189842, 190514, 190754, 191474, 191870, 192286, 192682, 193478, 193694, 194430, 194750, 195278, 195678, 196398, 196662, 197190, 197550, 198086, 198486, 199294, 199510, 200320, 200656, 201196, 201556, 202204, 202460, 203216, 203624, 204056, 204376, 205196, 205468, 206290, 206698, 207098, 207446, 208272, 208536, 209364, 209692, 210244, 210628, 211300, 211576, 212240, 212600, 213140, 213558, 214396, 214588, 215400, 215820, 216380, 216800, 217424, 217700, 218360, 218776, 219340, 219660, 220452, 220732, 221584, 221944, 222376, 222800, 223656, 223896, 224754, 225090, 225570, 226000, 226862, 227150, 227838, 228270, 228814, 229174, 229954, 230178, 230970, 231402, 231978, 232374, 232974, 233262, 234138, 234576, 235160, 235480, 236360, 236612, 237494, 237878, 238342, 238784, 239670, 239958, 240714, 241066, 241606, 242050, 242878, 243174, 243886, 244270, 244798, 245246, 246086, 246326, 247158, 247558, 248062, 248510, 249230, 249530, 250436, 250888, 251488, 251776, 252686, 252974, 253794, 254250, 254730, 255186, 255966, 256254, 257172, 257524, 258136, 258596, 259436, 259676, 260396, 260858, 261470, 261918, 262846, 263086, 263842, 264306, 264926, 265392, 266032, 266320, 267256, 267652, 268276, 268644, 269584, 269896, 270776, 271240, 271672, 272092, 273038, 273350, 274214, 274574, 275206, 275590, 276542, 276854, 277614, 278090, 278650, 279128, 279944, 280200, 281130, 281562, 282198, 282678, 283446, 283710, 284676, 285116, 285692, 286076, 287046, 287370, 288198, 288684, 289164, 289644, 290620, 290944, 291824, 292160, 292808, 293298, 294280, 294600, 295384, 295832, 296384, 296816, 297740, 297980, 298970, 299450, 300110, 300530, 301322, 301650, 302646, 303144, 303792, 304192]
n = 65536
from functools import cache
print(isqrt(n))
@cache
def phi(n):
    c = isqrt(n)    
    ans = n * (n + 1) // 2
    for i in range(2, n // c + 1):
        ans -= PHI[n // i] if n // i < 1000 else phi(n // i)
    for j in range(1, c):
        ans -= (n // j - n // (j + 1)) * (PHI[j] if j < 1000 else phi(j))
    return ans

print(phi(10**6))
from time import time
from math import log, isqrt, inf
from itertools import takewhile, compress, count




def primegen(limit=inf):
    """
    Generates primes strictly less than limit almost-lazily by a segmented
    sieve of Eratosthenes.  Memory usage depends on the sequence of prime
    gaps; on Cramer's conjecture, it is O(sqrt(p) * log(p)^2), where p is
    the most-recently-yielded prime.
    
    Input: limit -- a number (default = inf)
    
    Output: sequence of integers
    
    Examples:
    
    >>> list(islice(primegen(), 19))
    [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67]
    
    >>> list(primegen(71))
    [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67]
    """
    # We do not sieve 2, so we ought to be able to get sigificant savings by halving the length of the sieve.
    # But the tiny extra computation involved in that seems to exceed the savings.
    yield from takewhile(lambda x: x < limit, (2,3,5,7,11,13,17,19,23,29,31,37,41,43,47))
    pl, pg = [3,5,7], primegen()
    for p in pl: next(pg)
    while True:
        lo = pl[-1]**2
        if lo >= limit: break
        pl.append(next(pg))
        hi = min(pl[-1]**2, limit)
        sieve = bytearray([True]) * (hi - lo)
        for p in pl: sieve[(-lo)%p::p] = bytearray([False]) * ((hi-1)//p - (lo-1)//p)
        yield from compress(range(lo,hi,2), sieve[::2])




def introot(n, r=2):    # TODO Newton iteration?
    """
    Returns the rth root of n, rounded to the nearest integer in the
    direction of zero.  Returns None if r is even and n is negative.
    
    Input:
        n -- an integer
        r -- a natural number or None
    
    Output: An integer
    
    Examples:
    
    >>> [introot(-729, 3), introot(-728, 3)]
    [-9, -8]
    
    >>> [introot(1023, 2), introot(1024, 2)]
    [31, 32]
    """
    if n < 0: return None if r%2 == 0 else -introot(-n, r)
    if n < 2: return n
    if r == 1: return n
    if r == 2: return isqrt(n)
    if r % 2 == 0: return introot(isqrt(n), r//2)
    lower = upper = 1 << (n.bit_length() // r)
    while lower ** r >  n: lower >>= 2
    while upper ** r <= n: upper <<= 2
    while lower != upper - 1:
        mid = (lower + upper) // 2
        m = mid**r
        if   m == n: return  mid
        elif m <  n: lower = mid
        elif m >  n: upper = mid
    return lower




def mobiussieve(limit=inf):
    """
    Uses a segmented sieve to compute the Mobius function for all positive
    integers strictly less than the input.
    
    The time- and space-complexities to iterate over the first n terms
    are within logarithmic factors of O(n) and O(sqrt(n)), respectively.
    
    Input: limit -- an integer.  Default == inf.
    
    Output: Sequence of integers
    
    Example:
    
    >>> list(mobiussieve(21))
    [1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0, -1, 1, 1, 0, -1, 0, -1, 0]
    """
    if limit <= 1: return
    yield 1
    pg = primegen()
    primes = [next(pg)]
    nextp = next(pg)
    lo, hi = 2, min(nextp**2, limit)
    # We can sieve up to hi - 1.
    while lo < limit:
        mobs = [1] * (hi - lo)
        for p in primes:
            for n in range((-lo) %   p  , hi - lo,  p ): mobs[n] *= -p
            for n in range((-lo) % (p*p), hi - lo, p*p): mobs[n]  =  0
        for n in range(hi - lo):
            m = mobs[n]
            if m == 0: continue
            if -lo-n < m < lo+n:
                if m > 0: mobs[n] = -1
                if m < 0: mobs[n] =  1
            else:
                if m > 0: mobs[n] =  1
                if m < 0: mobs[n] = -1
        
        yield from mobs
        
        primes.append(nextp)
        nextp = next(pg)
        lo, hi = hi, min(nextp**2, limit)




def totientsum(n):
    if n <= 10: return 0 if n < 0 else (0,1,2,4,6,10,12,18,22,28,32)[n]
    
    a = introot(int((n / log(log(n)))**2), 3)
    b = n // a
    nr = isqrt(n)
    if verbose:
        print("ln(ln(n)):", log(log(n)))
        print("sqrt(n)/b:", nr//b)
        print("  sqrt(a):", isqrt(a))
        print("        b:", b)
        print("  sqrt(n):", nr)
        print("        a:", a)
        starttime = time()
        p2batchnum = 0
    Mover = [0] * (b + 1)  # Mover[n//x] will store Mertens(x) for large x.
    Mblock = []
    
    mert = 0
    X, Y, Z = 0, 0, 0
    
    s = nr - (nr == n//nr)
    chi = n // s
    
    d = b
    xp = isqrt(n//d)
    
    for (x, mu) in enumerate(mobiussieve(a+1), start=1):
        #if x == 100: exit()
        v = int(str(n // x))    # The int(str( ... )) pushes us back down into the 64-bit data types, when applicable.
        mert += mu
        X += mu * (v * (v+1) // 2)
        
        if x <= nr:
            Mblock.append(mert)
            
            if x > 1 and mu != 0:
                if verbose and (x<10*b or x%b<2): print("\b"*80, " Phase 1:", x, "%f%%" % (100*x/nr), end='    ', flush=True)
                if mu > 0:
                    for y in range(1, min(b, v//x) + 1):
                        Mover[y] -= v // y
                else:
                    for y in range(1, min(b, v//x) + 1):
                        Mover[y] += v // y
            
            while x == xp:
                Mover[d] += 1 - (n//d) + x * mert
                d = (d - 1) if d > 1 else n
                xp = isqrt(n//d)
            
            if x % b == 0 or x == nr:
                if verbose: print("\b"*80, " P1 Batch", x//b, "%f%%" % (100*x/nr), end='    ', flush=True)
                A = 1 + (b * (x//b) if x % b != 0 else (x - b))
                for t in range(1, b+1):
                    #if verbose and x < 20*b: print("\b"*80, "         ", x//b, t, "%f%%" % (100*x/a), end='    ', flush=True)
                    nt = n // t
                    lmin = 1 + n // (t * (x+1))
                    lmax = min(isqrt(n//t), nt // A)
                    for l in range(lmin, lmax + 1):
                        k = nt // l
                        assert A <= k <= x
                        Mover[t] -= Mblock[k - A]
                Mblock.clear()
                if verbose and x == nr:
                    phase1time = time()
                    print("\b"*80 + ("Phase 1 took %f seconds.    " % (phase1time - starttime)))
        
        elif x == chi:
            if verbose and len(Mblock) % 100 == 0: print("\b"*80, " Phase 2:", x, "%f%%" % (100*x/a), end='    ', flush=True)
            if v != b:
                if len(Mblock) == 0: A = v
                Mblock.append(mert)
                B = v
            s -= 1
            chi = n // s
        
        if (x == a and len(Mblock) > 0) or (x > nr and len(Mblock) == b):
            if verbose:
                p2batchnum += 1
                print("\b"*80, " P2 Batch", p2batchnum, "%f%%" % (100*x/a), end='    ', flush=True)
            BnBn = B * n // (B + n)
            for y in range(1, b+1):
                ny = n // y
                tmin = max(2, BnBn // y)
                tmax = min(isqrt(ny), (A + 1) // y)
                for t in range(tmin, tmax+1):
                    nty = ny // t
                    if nr < nty:
                        if B <= n // nty <= A:
                            Mover[y] -= Mblock[A - t*y]
                    else: break
            Mblock.clear()
        
    Z = mert * (b * (b+1) // 2)
    
    if verbose:
        phase2time = time()
        print("\b"*80 + ("Phase 2 took %f seconds.    " % (phase2time - phase1time)))
    
    for y in range(b, 0, -1):
        v = n // y
        vr = isqrt(v)
        Mv = 0
        for t in count(2):
            if n >= (b+1) * (v//t): break
            Mv -= Mover[n//(v//t)]
        # Mv is now Mertens(v).
        Mover[y] += Mv
        Y += y * Mover[y]
    
    if verbose:
        phase3time = time()
        print("\b"*80 + ("Phase 3 took %f seconds." % (phase3time - phase2time)))
    
    return X + Y - Z




# verbose = False
# for n in (10**5,10**6,10**7):
#     print(n, end=' ', flush=True)
#     print(totientsum(n))

Embed on website

To embed this program on your website, copy the following code and paste it into your website's HTML: