#!/usr/bin/env python3
"""Number of n < N such that n and n+1 have the same number of divisors."""
import sys
from itertools import islice, starmap, tee
def ndiv(n, prime_factors):
"""Return number of divisors of `n`.
prime_factors: prime factors of `n`.
>>> from itertools import starmap
>>> list(starmap(ndiv, ((1, []), (12, [2, 3]))))
[1, 6]
"""
assert n > 0
phi = 1
for prime in prime_factors:
alpha = 0 # multiplicity of `prime` in `n`
q, r = divmod(n, prime)
while r == 0: # `prime` is a factor of `n`
n = q
alpha += 1
q, r = divmod(n, prime)
phi *= alpha + 1 # see http://e...content-available-to-author-only...a.org/wiki/Divisor_function
return phi
def prime_factors_gen():
"""Yield prime factors for each natural number.
Based on
http://stackoverflow.com/questions/567222/simple-prime-generator-in-python/568618#568618
>>> from itertools import islice
>>> list(islice(prime_factors_gen(), 20)) #doctest:+NORMALIZE_WHITESPACE
[(1, []), (2, [2]), (3, [3]), (4, [2]), (5, [5]), (6, [3, 2]),
(7, [7]), (8, [2]), (9, [3]), (10, [5, 2]), (11, [11]), (12, [3, 2]),
(13, [13]), (14, [7, 2]), (15, [5, 3]), (16, [2]), (17, [17]),
(18, [3, 2]), (19, [19]), (20, [5, 2])]
"""
D = {1:[]} # nonprime -> prime factors of `nonprime`
#NOTE: dictionary could be replaced by priority queue
q = 1
while True: # Sieve of Eratosthenes algorithm
if q not in D: # `q` is a prime number
D[q + q] = [q]
yield q, [q]
else: # q is a composite
for p in D[q]: # `p` is a factor of `q`
# therefore `p` is a factor of `p + q` too
D.setdefault(p + q, []).append(p)
yield q, D[q]
del D[q]
q += 1
def pairwise(it):
a, b = tee(it)
next(b, None)
return zip(a, b)
N = int(input('Enter max n: '))
taus = islice(starmap(ndiv, prime_factors_gen()), N)
print(sum(x == y for x, y in pairwise(taus)))
IyEvdXNyL2Jpbi9lbnYgcHl0aG9uMwoiIiJOdW1iZXIgb2YgbiA8IE4gc3VjaCB0aGF0IG4gYW5kIG4rMSBoYXZlIHRoZSBzYW1lIG51bWJlciBvZiBkaXZpc29ycy4iIiIKaW1wb3J0IHN5cwpmcm9tIGl0ZXJ0b29scyBpbXBvcnQgaXNsaWNlLCBzdGFybWFwLCB0ZWUKCmRlZiBuZGl2KG4sIHByaW1lX2ZhY3RvcnMpOgogICAgIiIiUmV0dXJuIG51bWJlciBvZiBkaXZpc29ycyBvZiBgbmAuCgogICAgcHJpbWVfZmFjdG9yczogcHJpbWUgZmFjdG9ycyBvZiBgbmAuCgogICAgPj4+IGZyb20gaXRlcnRvb2xzIGltcG9ydCBzdGFybWFwCiAgICA+Pj4gbGlzdChzdGFybWFwKG5kaXYsICgoMSwgW10pLCAoMTIsIFsyLCAzXSkpKSkKICAgIFsxLCA2XQogICAgIiIiCiAgICBhc3NlcnQgbiA+IDAKICAgIHBoaSA9IDEKICAgIGZvciBwcmltZSBpbiBwcmltZV9mYWN0b3JzOgogICAgICAgIGFscGhhID0gMCAjIG11bHRpcGxpY2l0eSBvZiBgcHJpbWVgIGluIGBuYAogICAgICAgIHEsIHIgPSBkaXZtb2QobiwgcHJpbWUpCiAgICAgICAgd2hpbGUgciA9PSAwOiAjIGBwcmltZWAgaXMgYSBmYWN0b3Igb2YgYG5gCiAgICAgICAgICAgIG4gPSBxCiAgICAgICAgICAgIGFscGhhICs9IDEKICAgICAgICAgICAgcSwgciA9IGRpdm1vZChuLCBwcmltZSkKICAgICAgICBwaGkgKj0gYWxwaGEgKyAxICMgc2VlIGh0dHA6Ly9lLi4uY29udGVudC1hdmFpbGFibGUtdG8tYXV0aG9yLW9ubHkuLi5hLm9yZy93aWtpL0Rpdmlzb3JfZnVuY3Rpb24KICAgIHJldHVybiBwaGkKCmRlZiBwcmltZV9mYWN0b3JzX2dlbigpOgogICAgIiIiWWllbGQgcHJpbWUgZmFjdG9ycyBmb3IgZWFjaCBuYXR1cmFsIG51bWJlci4KCiAgICBCYXNlZCBvbgogICAgaHR0cDovL3N0YWNrb3ZlcmZsb3cuY29tL3F1ZXN0aW9ucy81NjcyMjIvc2ltcGxlLXByaW1lLWdlbmVyYXRvci1pbi1weXRob24vNTY4NjE4IzU2ODYxOAoKICAgID4+PiBmcm9tIGl0ZXJ0b29scyBpbXBvcnQgaXNsaWNlCiAgICA+Pj4gbGlzdChpc2xpY2UocHJpbWVfZmFjdG9yc19nZW4oKSwgMjApKSAjZG9jdGVzdDorTk9STUFMSVpFX1dISVRFU1BBQ0UKICAgIFsoMSwgW10pLCAoMiwgWzJdKSwgKDMsIFszXSksICg0LCBbMl0pLCAoNSwgWzVdKSwgKDYsIFszLCAyXSksCiAgICAoNywgWzddKSwgKDgsIFsyXSksICg5LCBbM10pLCAoMTAsIFs1LCAyXSksICgxMSwgWzExXSksICgxMiwgWzMsIDJdKSwKICAgICgxMywgWzEzXSksICgxNCwgWzcsIDJdKSwgKDE1LCBbNSwgM10pLCAoMTYsIFsyXSksICgxNywgWzE3XSksCiAgICAoMTgsIFszLCAyXSksICgxOSwgWzE5XSksICgyMCwgWzUsIDJdKV0KICAgICIiIgogICAgRCA9IHsxOltdfSAjIG5vbnByaW1lIC0+IHByaW1lIGZhY3RvcnMgb2YgYG5vbnByaW1lYAogICAgI05PVEU6IGRpY3Rpb25hcnkgY291bGQgYmUgcmVwbGFjZWQgYnkgcHJpb3JpdHkgcXVldWUKICAgIHEgPSAxCiAgICB3aGlsZSBUcnVlOiAjIFNpZXZlIG9mIEVyYXRvc3RoZW5lcyBhbGdvcml0aG0KICAgICAgICBpZiBxIG5vdCBpbiBEOiAjIGBxYCBpcyBhIHByaW1lIG51bWJlcgogICAgICAgICAgICBEW3EgKyBxXSA9IFtxXQogICAgICAgICAgICB5aWVsZCBxLCBbcV0KICAgICAgICBlbHNlOiAjIHEgaXMgYSBjb21wb3NpdGUKICAgICAgICAgICAgZm9yIHAgaW4gRFtxXTogIyBgcGAgaXMgYSBmYWN0b3Igb2YgYHFgCiAgICAgICAgICAgICAgICAjIHRoZXJlZm9yZSBgcGAgaXMgYSBmYWN0b3Igb2YgYHAgKyBxYCB0b28KICAgICAgICAgICAgICAgIEQuc2V0ZGVmYXVsdChwICsgcSwgW10pLmFwcGVuZChwKQogICAgICAgICAgICB5aWVsZCBxLCBEW3FdCiAgICAgICAgICAgIGRlbCBEW3FdCiAgICAgICAgcSArPSAxCgpkZWYgcGFpcndpc2UoaXQpOgogICAgYSwgYiA9IHRlZShpdCkKICAgIG5leHQoYiwgTm9uZSkKICAgIHJldHVybiB6aXAoYSwgYikKCgpOID0gaW50KGlucHV0KCdFbnRlciBtYXggbjogJykpCnRhdXMgPSBpc2xpY2Uoc3Rhcm1hcChuZGl2LCBwcmltZV9mYWN0b3JzX2dlbigpKSwgTikKcHJpbnQoc3VtKHggPT0geSBmb3IgeCwgeSBpbiBwYWlyd2lzZSh0YXVzKSkpCg==