12 Haziran 2018

Euler Projesi 241. Soru

Mükemmellik Bölümleri

Pozitif bir n tam sayı için σ(n), n'nin tüm bölenlerinin toplamı olsun. Örneğin σ(6) = 1 + 2 + 3 + 6 = 12.

Muhtemelen bildiğiniz üzere mükemmel bir sayı için σ(n) = 2n dir.

Pozitif bir tam sayının mükemmellik bölümünü p(n) = σ(n)/n olarak tanımlayalım.

k bir tam sayı olmak üzere p(n)'nin k + 1⁄2 formuna sahip olduğu tüm pozitif n ≤ 1018 tam sayılarının toplamını bulunuz.
Cevap: 482316491800641154

Python:
import math
                               
MAX_PRIME = 10 ** 3

def isprime(x) :  
  if x <= 1 :
    return False
  if x < 4 :
    return True
  if x % 2 == 0 :
    return False
  if x < 9 :
    return True
  if x % 3 == 0 :
    return False
  
  r = int(math.sqrt(x))
  f = 5
  while f <= r :
    if x % f == 0 :
      return False
    if x % (f+2) == 0 :
      return False
    f += 6
  return True 


p = [x for x in range(MAX_PRIME + 1) if isprime(x)]
 
def factor(x):
    factors = {}
    for q in p:
        if q * q > x:
            break
        i = 0
        while x % q == 0:
            x //= q
            i += 1
        if i:
            factors[q] = i
    if x > 1:
        factors[x] = 1
    return factors


def sigma(prime, power) :
  res = 0
  m = 1
  for x in xrange(0, power + 1) :
    res += m
    m *= prime
  return factor(res)




MAX =  10 ** 18 

total = 0

def walk(current, factors, balance, number) :

  # check end  
  if balance.count(0) == len(balance) :
    print "found: ", number, factors
    global total
    total += number
    return
  
  if len(factors) > 30 :
    return

  if current == len(balance) :
    return

  # try powers
  start = balance[current]
  if start < 1 :
    walk(current + 1, factors, balance, number) # for zero
    start = 1

  p = factors[current]

  for power in xrange(start, 30) :

    f = sigma(p, power)   

    number_new = number * p ** power

    #check if applicable
    if number_new > MAX :
      break

    apply = True
    for v in f :           
      if v in factors :
        i = factors.index(v)
        if i < current and balance[i] + f[v] > 0 :
          apply = False
        break    

    # apply
    if apply :
      balance_new = balance[:]
      factors_new = factors[:]

      balance_new[current] -= power

      for v in f :           
        if not v in factors_new :
          factors_new.append(v)
          balance_new.append(f[v])
        else :
          i = factors.index(v)        
          balance_new[i] += f[v]

      walk(current + 1, factors_new, balance_new, number_new)
        
   
for x in xrange(3, 15, 2) :
  print "Looking: ", x / 2.0
  if x == 9 :
    balance = [1, -2]
    factors = [2, 3]
  else :
    balance = [1, -1]
    factors = [2, x]
  walk(0, factors, balance, 1)


print "Total: ", total

Maple:
pb241:=proc(a,b,N,p,primes)
local res,p1,aa,bb,nn,alpha,g,sp,p2,f,pp2;

if a=b and a=1 then return([1]) fi:

if p>min(N,200) then return([]) fi:

if a=1
then
  p1:=p:alpha:=1:
  p2:=nextprime(p):
  while member(p2,primes) do p2:=nextprime(p2) od:
  res:=op(pb241(a,b,N,p2,{op(primes),p1}))
else
  f:=ifactors(a)[2]:
  if f[1][1]=2 then p1:=f[1][1]:alpha:=f[1][2] 
  else p1:=f[nops(f)][1]:alpha:=f[nops(f)][2] fi: # choose the biggest
  if f[1][1]
p then p2:=p else p2:=nextprime(p):while member(p2,primes) do p2:=nextprime(p2) od: fi: res:=NULL fi: while p1**alpha<=N do sp:=(p1**(alpha+1)-1)/(p1-1):g:=igcd(a*sp,b*p1**alpha): res:=res,op(p1**alpha*pb241(a*sp/g,b*p1**alpha/g,N/p1**alpha,p2,{op(primes),p1})): alpha:=alpha+1 od: return([res]): end:

C++:
#include 

/*
ans cand: 2
ans cand: 24
ans cand: 4680
ans cand: 26208
ans cand: 17428320
ans cand: 4320
ans cand: 8910720
ans cand: 91963648
ans cand: 197064960
ans cand: 20427264
ans cand: 88898072401645056
ans cand: 8583644160
ans cand: 57629644800
ans cand: 206166804480
ans cand: 17116004505600
ans cand: 75462255348480000
ans cand: 57575890944
ans cand: 10200236032
ans cand: 21857648640
ans cand: 1416963251404800
ans cand: 15338300494970880
ans cand: 301183421949935616
result: 482316491800641154

*/

#define N 1000000000000000000LL
#define SQRN 1000000000LL

int prime_t[100000];
int prime_f[100000];
int prime_search[SQRN/10];
int primen;

int myprime_t[200];
int myprime_f[200];
int myprime_i[200];
int myprimen;

#define ANSWERN 100000
long long int answer[ANSWERN];
int answern = 0;
long long int sum = 0;

int target2;
int sum2;
long long int deno;
long long int pow2_deno = 1;

void make_primetable();
long long int mypow_sum(int b, int n);
void prime_deno();
int answer_check(long long int ans);

int main(int argc, char **argv) {
  int i;
  int j, k;
  int s;
  int max;
  
  long long int pow2_add = 1;
  int primemax;

  long long int target;
  make_primetable();

  int s2count;


  // main routine: check 2^n
  for (i=1;; i++) {
    pow2_deno *= 2;
    pow2_add = pow2_add * 2+1;
    if (pow2_deno > N /5)
      break;
    //printf ("%d %lld %lld %lld\n", i, pow2_deno, pow2_add, pow2_add*pow2_deno);

    for (j=0; j target2)
    return;
  if (sum2 == target2)
    answer_check(deno);


  myprimen_l = myprimen;
  deno_l = deno;

  for (i=0; i N)
 break;

      // numerator
      target = mypow_sum(myprime_t[i], myprime_i[i]);


      // prime decompostion
      mymyprimen = 0;
      sqrt_t = sqrt(target);
      for (j=0; j< primen; j++) {
 
 if (sqrt_t < prime_t[j])
   break;
 if (!(target % prime_t[j])) {
   target /= prime_t[j];
   mymyprime_t[mymyprimen] = prime_t[j];
   mymyprime_f[mymyprimen++] = 1;
   while (!(target % prime_t[j])) {
     mymyprime_f[mymyprimen-1]++;
     target /= prime_t[j];
   }
   sqrt_t = sqrt(target);
 }
      }
      if (target != 1) {
 mymyprime_t[mymyprimen] = target;
 mymyprime_f[mymyprimen++] = 1;
      }
 

      for (j=0; j myprime_f[i])
      return -1;


  answer[answern++] = ans;
  sum += ans;

  printf ("ans cand: %lld\n", ans);

  if (answern >= ANSWERN)
    printf ("overflow\n");

  return 1;
}