[PYTHON] Projecet Euler 12 Find the number of divisors without division.

problem

The sequence of triangular numbers is represented by the sum of natural numbers, and the 7th triangular number is 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. The first 10 terms of the triangular number are:

1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ... Will be.

The divisors of the first 7 terms are listed below.

1: 1 3: 1,3 6: 1,2,3,6 10: 1,2,5,10 15: 1,3,5,15 21: 1,3,7,21 28: 1,2,4,7,14,28

From this, it can be seen that the seventh triangular number, 28, is the first triangular number with more than five divisors.

So, what are the first triangular numbers with divisors greater than 500? http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2012

Mathematical consideration

Triangular numbers are difficult, but the point is that triangular numbers are the sum of natural numbers from 1 to n. The sum of natural numbers from 1 to n is expressed by the following formula.

n * (n+1) /2

And since n and (n + 1) are relatively prime (have no common prime factors), n / 2 and n + 1 or n and (n + 1) / 2 are also relatively prime. Whether it should be n / 2 and n + 1 or n and (n + 1) / 2 depends on which is even. The number of divisors of the number obtained by multiplying the coprime natural numbers n1 and n2 is equal to the number obtained by multiplying the number of divisors of n1 and the number of divisors of n2.

Answer policy 1

Based on the above mathematical considerations, we will answer this question using the following algorithm.

Make a list of prime numbers using mymath. 　mymath.py 　http://qiita.com/cof/items/45d3823c3d71e7e22920
While the number of divisors is 500 or less, the variable n is incremented in order from 1 in the loop loop.
After determining whether n + 1 is even or not, perform n + 1 or (n + 1) / 2 prime factorization using the prime number list.
Calculate the number of divisors of n + 1 below the obtained prime factors, and multiply it by the number of divisors of n previously memorized.
Memorize the multiplied number as a divisor and return to 2.

code

import mymath

def cof():
  pri = mymath.get_primes(10**4)
  n = 1
  num_fac = 1
  lim = 500
  n1_fac =  mymath.get_number_of_factor(n,pri)
  while num_fac < lim:
    n += 1
    if (n+1)%2:
      n2_fac =  mymath.get_number_of_factor(n+1,pri)
    else:
      n2_fac =  mymath.get_number_of_factor((n+1)/2,pri)
    num_fac = n1_fac * n2_fac
    n1_fac = n2_fac
  print n * (n+1) / 2

Discovery of problems and new awareness

When I tried to execute it, it was 1.62 seconds per execution, which is slow at the poop level. When I wondered if I could speed it up somehow, I realized that factoring was the cause of the slowness in the first place. I noticed that the divisor of each number can be obtained without division by the same method as Eratosthenes.

Answer policy 2

For more details,

Make a list for each number. "" "," ", ..."
Loop from 1 to an appropriate number and add the number of targets as an element at a position that is a multiple of the number of targets. For example, if 2 is the target number [[]. [1], [1], [1], [1] ← Add 2 to the list corresponding to multiple 4, the same as 6,8,10 below.
Repeat 2 and, strangely, you will have a list of true divisors of the corresponding numbers. [[].[1],[1],[1],[1,2],[1],[1,2,3],[1],[1,2,4],[1,3],[1,2,5],...]
Based on this divisor list, find the number of divisors of n in response policy 1.

code

Added the following functions to mymath.

def factor_seq(max):
  ret = [[0]]
  for i in range(max):
    ret += [[1]]
  seq = range(max+1)
  for i in seq[2:max//2+1]:
    for j in seq[i*2::i]:
      ret[j] += [i]
  return ret

Coded using the above

def cof2():
  fq = mymath2.factor_seq(10**5)
  n = 1
  num_fac = 1
  lim = 500
  n1_fac = len(fq[n])
  while num_fac < lim:
    n += 1
    if (n+1)%2:
      n2_fac =  len(fq[n+1])
    else:
      n2_fac =  len(fq[(n+1)/2])
    num_fac = n1_fac * n2_fac
    n1_fac = n2_fac
  print n * (n+1) / 2

As a result, it was 640ms, which was a level that I could put up with. Probably, if you change ret [j] + = [i] to ret [j] + = 1 in factor_seq, you can find the number of divisors, so I think it will be faster if you use it.