fork(1) download
  1. #!/usr/bin/env python3
  2. """Number of n < N such that n and n+1 have the same number of divisors."""
  3. import sys
  4. from itertools import islice, starmap, tee
  5.  
  6. def ndiv(n, prime_factors):
  7. """Return number of divisors of `n`.
  8.  
  9. prime_factors: prime factors of `n`.
  10.  
  11. >>> from itertools import starmap
  12. >>> list(starmap(ndiv, ((1, []), (12, [2, 3]))))
  13. [1, 6]
  14. """
  15. assert n > 0
  16. phi = 1
  17. for prime in prime_factors:
  18. alpha = 0 # multiplicity of `prime` in `n`
  19. q, r = divmod(n, prime)
  20. while r == 0: # `prime` is a factor of `n`
  21. n = q
  22. alpha += 1
  23. q, r = divmod(n, prime)
  24. phi *= alpha + 1 # see http://e...content-available-to-author-only...a.org/wiki/Divisor_function
  25. return phi
  26.  
  27. def prime_factors_gen():
  28. """Yield prime factors for each natural number.
  29.  
  30. Based on
  31. http://stackoverflow.com/questions/567222/simple-prime-generator-in-python/568618#568618
  32.  
  33. >>> from itertools import islice
  34. >>> list(islice(prime_factors_gen(), 20)) #doctest:+NORMALIZE_WHITESPACE
  35. [(1, []), (2, [2]), (3, [3]), (4, [2]), (5, [5]), (6, [3, 2]),
  36. (7, [7]), (8, [2]), (9, [3]), (10, [5, 2]), (11, [11]), (12, [3, 2]),
  37. (13, [13]), (14, [7, 2]), (15, [5, 3]), (16, [2]), (17, [17]),
  38. (18, [3, 2]), (19, [19]), (20, [5, 2])]
  39. """
  40. D = {1:[]} # nonprime -> prime factors of `nonprime`
  41. #NOTE: dictionary could be replaced by priority queue
  42. q = 1
  43. while True: # Sieve of Eratosthenes algorithm
  44. if q not in D: # `q` is a prime number
  45. D[q + q] = [q]
  46. yield q, [q]
  47. else: # q is a composite
  48. for p in D[q]: # `p` is a factor of `q`
  49. # therefore `p` is a factor of `p + q` too
  50. D.setdefault(p + q, []).append(p)
  51. yield q, D[q]
  52. del D[q]
  53. q += 1
  54.  
  55. def pairwise(it):
  56. a, b = tee(it)
  57. next(b, None)
  58. return zip(a, b)
  59.  
  60.  
  61. N = int(input('Enter max n: '))
  62. taus = islice(starmap(ndiv, prime_factors_gen()), N)
  63. print(sum(x == y for x, y in pairwise(taus)))
  64.  
Success #stdin #stdout 0.23s 28776KB
stdin
100000
stdout
Enter max n: 10585