#!/usr/bin/env python
"""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(taus)
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)))
IyEvdXNyL2Jpbi9lbnYgcHl0aG9uCiIiIk51bWJlciBvZiBuIDwgTiBzdWNoIHRoYXQgbiBhbmQgbisxIGhhdmUgdGhlIHNhbWUgbnVtYmVyIG9mIGRpdmlzb3JzLiIiIgppbXBvcnQgc3lzCmZyb20gaXRlcnRvb2xzIGltcG9ydCBpc2xpY2UsIHN0YXJtYXAsIHRlZQoKZGVmIG5kaXYobiwgcHJpbWVfZmFjdG9ycyk6CiAgICAiIiJSZXR1cm4gbnVtYmVyIG9mIGRpdmlzb3JzIG9mIGBuYC4KCiAgICBwcmltZV9mYWN0b3JzOiBwcmltZSBmYWN0b3JzIG9mIGBuYC4KCiAgICA+Pj4gZnJvbSBpdGVydG9vbHMgaW1wb3J0IHN0YXJtYXAKICAgID4+PiBsaXN0KHN0YXJtYXAobmRpdiwgKCgxLCBbXSksICgxMiwgWzIsIDNdKSkpKQogICAgWzEsIDZdCiAgICAiIiIKICAgIGFzc2VydCBuID4gMAogICAgcGhpID0gMQogICAgZm9yIHByaW1lIGluIHByaW1lX2ZhY3RvcnM6CiAgICAgICAgYWxwaGEgPSAwICMgbXVsdGlwbGljaXR5IG9mIGBwcmltZWAgaW4gYG5gCiAgICAgICAgcSwgciA9IGRpdm1vZChuLCBwcmltZSkKICAgICAgICB3aGlsZSByID09IDA6ICMgYHByaW1lYCBpcyBhIGZhY3RvciBvZiBgbmAKICAgICAgICAgICAgbiA9IHEKICAgICAgICAgICAgYWxwaGEgKz0gMQogICAgICAgICAgICBxLCByID0gZGl2bW9kKG4sIHByaW1lKQogICAgICAgIHBoaSAqPSBhbHBoYSArIDEgIyBzZWUgaHR0cDovL2UuLi5jb250ZW50LWF2YWlsYWJsZS10by1hdXRob3Itb25seS4uLmEub3JnL3dpa2kvRGl2aXNvcl9mdW5jdGlvbgogICAgcmV0dXJuIHBoaQoKZGVmIHByaW1lX2ZhY3RvcnNfZ2VuKCk6CiAgICAiIiJZaWVsZCBwcmltZSBmYWN0b3JzIGZvciBlYWNoIG5hdHVyYWwgbnVtYmVyLgoKICAgIEJhc2VkIG9uCiAgICBodHRwOi8vc3RhY2tvdmVyZmxvdy5jb20vcXVlc3Rpb25zLzU2NzIyMi9zaW1wbGUtcHJpbWUtZ2VuZXJhdG9yLWluLXB5dGhvbi81Njg2MTgjNTY4NjE4CgogICAgPj4+IGZyb20gaXRlcnRvb2xzIGltcG9ydCBpc2xpY2UKICAgID4+PiBsaXN0KGlzbGljZShwcmltZV9mYWN0b3JzX2dlbigpLCAyMCkpICNkb2N0ZXN0OitOT1JNQUxJWkVfV0hJVEVTUEFDRQogICAgWygxLCBbXSksICgyLCBbMl0pLCAoMywgWzNdKSwgKDQsIFsyXSksICg1LCBbNV0pLCAoNiwgWzMsIDJdKSwKICAgICg3LCBbN10pLCAoOCwgWzJdKSwgKDksIFszXSksICgxMCwgWzUsIDJdKSwgKDExLCBbMTFdKSwgKDEyLCBbMywgMl0pLAogICAgKDEzLCBbMTNdKSwgKDE0LCBbNywgMl0pLCAoMTUsIFs1LCAzXSksICgxNiwgWzJdKSwgKDE3LCBbMTddKSwKICAgICgxOCwgWzMsIDJdKSwgKDE5LCBbMTldKSwgKDIwLCBbNSwgMl0pXQogICAgIiIiCiAgICBEID0gezE6W119ICMgbm9ucHJpbWUgLT4gcHJpbWUgZmFjdG9ycyBvZiBgbm9ucHJpbWVgCiAgICAjTk9URTogZGljdGlvbmFyeSBjb3VsZCBiZSByZXBsYWNlZCBieSBwcmlvcml0eSBxdWV1ZQogICAgcSA9IDEKICAgIHdoaWxlIFRydWU6ICMgU2lldmUgb2YgRXJhdG9zdGhlbmVzIGFsZ29yaXRobQogICAgICAgIGlmIHEgbm90IGluIEQ6ICMgYHFgIGlzIGEgcHJpbWUgbnVtYmVyCiAgICAgICAgICAgIERbcSArIHFdID0gW3FdCiAgICAgICAgICAgIHlpZWxkIHEsIFtxXQogICAgICAgIGVsc2U6ICMgcSBpcyBhIGNvbXBvc2l0ZQogICAgICAgICAgICBmb3IgcCBpbiBEW3FdOiAjIGBwYCBpcyBhIGZhY3RvciBvZiBgcWAKICAgICAgICAgICAgICAgICMgdGhlcmVmb3JlIGBwYCBpcyBhIGZhY3RvciBvZiBgcCArIHFgIHRvbwogICAgICAgICAgICAgICAgRC5zZXRkZWZhdWx0KHAgKyBxLCBbXSkuYXBwZW5kKHApCiAgICAgICAgICAgIHlpZWxkIHEsIERbcV0KICAgICAgICAgICAgZGVsIERbcV0KICAgICAgICBxICs9IDEKCmRlZiBwYWlyd2lzZShpdCk6CiAgICBhLCBiID0gdGVlKHRhdXMpCiAgICBuZXh0KGIsIE5vbmUpCiAgICByZXR1cm4gemlwKGEsIGIpCgoKTiA9IGludChpbnB1dCgnRW50ZXIgbWF4IG46ICcpKQp0YXVzID0gaXNsaWNlKHN0YXJtYXAobmRpdiwgcHJpbWVfZmFjdG9yc19nZW4oKSksIE4pCnByaW50KHN1bSh4ID09IHkgZm9yIHgsIHkgaW4gcGFpcndpc2UodGF1cykpKQo=