尤拉計劃699:三腳數

Triffle Numbers

Problem 699 Let σ(n) be the sum of all the divisors of the positive integer n, for example: σ(10)=1+2+5+10=18.

Define T(N) to be the sum of all numbers n≤N such that when the fraction σ(n)/n is written in its lowest form a/b, the denominator is a power of 3 i.e. b=3^k,k>0.

You are given T(100)=270 and T(10^6)=26089287.

Find T(10^14)

測試程式

def sieve(n):
 primes=[]
 h  = [True] * int(n**0.5+1)
 h[:2] = [False, False]
 for i in range(2, int(n ** 0.25) + 1):
    if h[i]:
        h[i*i::i] = [False] * len(h[i*i::i])
 for i, e in enumerate(h):
    if e:
        primes.append(i)
 return primes

def gcd(a,b):
 if a<b: a,b=b,a
 while(b>0):a,b=b,a%b
 return a

import math

def divs(n):
 d=[]
 c=0
 st=math.sqrt(n)
 for i in p:
  if i>st:break
  while n % i ==0: n//=i;c+=1  #;print(i,c);
  else:
   if(c>0):d.append([i,c]);c=0
 if n>st: d.append([n,1])
 return d

import itertools

def tdiv(n):
 t=divs(n)
 s=0
 z=[[i[0]**j for j in range(i[1]+1)] for i in t]
 #[[for i in range(len(item)]for item in itertools.product(*z)]
 s=0
 for item in itertools.product(*z):
  s1=1
  for i in item: s1*=i
  s+=s1
 return s

def three(n):
 if n==1:return 0
 while n%3==0:n//=3
 return n

n=10**6
p=sieve(n)
def t(n):
 s=0
 for i in range(3,n+1,3):
  deno=i/gcd(tdiv(i),i)
  #if 3**int(math.log(deno,3))==deno:s+=i
  if three(deno)==1:s+=i #;print(i,s)
 return s

>>> t(10**6)
26089287