Euler Projesi 234. Soru

Yarıbölünebilir Sayılar

Bir n ≥ 4 tam sayısı için n'nin lps(n) ile gösterilen alt asal karekökünü en büyük asal ≤ √n ve n'nin ups(n) ile gösterilen üst asal karekökünü en küçük asal sayı ≥ √n olarak tanımlarız.

Yani, örneğin, lps (4) = 2 = up (4), lps (1000) = 31, ups (1000) = 37.
lps(n) ve ups(n)'den sadece biri n'yi bölüyorsa n ≥ 4 tam sayısına yarıbölünebilir sayı diyelim.

15'i aşmayan yarıbölünebilir sayıların toplamı 30'dur; bunlar 8, 10 ve 12'dir.
15, hem lps(15) = 3 hem de ups(15) = 5'in bir katı olduğu için yarıbölünebilir değildir.
Başka bir örnek olarak, 1000'e kadar olan 92 yarıbölünebilir sayının toplamı 34825'tir.

999966663333'ü aşmayan tüm yarıbölünebilir sayıların toplamı nedir?

Cevap: 1259187438574927161

C++:

lint Sum(lint M1, lint M2, lint d) {
    M1 = (M1 - 1) / d + 1;
    M2 = M2 / d;
    lint ans = (M2-M1+1)*(M2+M1);
    ans /= 2;
    ans *= d;
    return ans;
}

lint Problem() {
 lint ans = 0;

    lint N = 999966663333;
    for (int i = 0; ; ++i) {
        lint p1 = viPrimes[i];
        lint p2 = viPrimes[i+1];
        if (p1 * p1 > N)
            break;
        lint M1 = p1*p1+1;
        lint M2 = p2*p2-1;
        if (M2 > N)
            M2 = N;
        ans += Sum(M1, M2, p1);
        ans += Sum(M1, M2, p2);
        ans -= 2*Sum(M1, M2, p1*p2);
    }

    return ans;
}

Python:

import math

nmax = 1000000

primes = nmax*[True]

primes[0:2] = [False,False]

for x in range(2, int(math.sqrt(len(primes)) + 1)):
    if primes[x]:
        for y in range(2 * x, len(primes), x):
            primes[y] = False

p = [x for x in range(len(primes)) if primes[x]]

s = 0

nmax = 999966663333

for i in range(len(p) - 1):
    n = p[i] ** 2
    if n > nmax:
        break
    n += p[i]
    while n < p[i+1] ** 2:
        if n <= nmax and n % p[i+1]:
            s += n
        n += p[i]
    n = p[i+1] ** 2 - p[i+1]
    while n > p[i] ** 2:
        if n <= nmax and n % p[i]:
            s += n
        n -= p[i+1]

print(s)

Java:

public class p234 {
 public static int CalcPrimes( int [] primes, int max ) {
  int  cnt;
  int  j, k;
  int          maxK;
  int  maxKPrimeSq;
  boolean         isPrime;
  
  primes[0] = 2;
  primes[1] = 3;
  cnt = 2;
  maxK = 0;
  maxKPrimeSq = primes[maxK] * primes[maxK];
  for( j=5; ; j+=2 ){
   if( j > max ) break;
   while( maxKPrimeSq < j ) {
    maxK ++;
    maxKPrimeSq = primes[maxK] * primes[maxK];
   }
   isPrime = true;
   for( k=0; k <= maxK; k++ ) {
    if( j % primes[k] == 0 ) {
     isPrime = false;
     break;
    }
   }
   if( isPrime ) {
    primes[ cnt ] = j;
    cnt ++;
    if( cnt >= primes.length ) {
     System.out.print( "prime array too small\n" );
     break;
    }
   }
  }
  return( cnt );
 }
 public static void main(String[] args) {
  long   start, stop;
  int          primes[] = new int[250000]; 
  int   primeCnt;
  int   j;
  long   min, max;
  long   test;
  long   sum   = 0;
  long   MAX_VAL  = 999966663333l;
  
  start = System.nanoTime();

  primeCnt = CalcPrimes( primes, 1000100 );

  
  for( j=0; j < primeCnt - 1; j++ ) {
   min = (long)primes[j] * (long)primes[j];
   max = (long)primes[ j+1 ] * (long)primes[ j+1 ];
   if( max > MAX_VAL + 1 ) max = MAX_VAL + 1;
   for( test = min + primes[j]; test < max; test += primes[j] ) {
    if( test % primes[j+1] != 0 ) {
     sum += test;
    }
   }
   for( test = min - ( min % primes[j+1] ) + primes[j+1]; 
                             test < max; test += primes[j+1] ) {
    if( test % primes[j] != 0 ) {
     sum += test;
    }
   }

  }
  stop = System.nanoTime();
  System.out.print( "Run time: " 
                          + ( String.valueOf( ((double)(stop-start))/ 1000000000)) 
                          + "\n");
  System.out.print( "Sum: " + String.valueOf(sum) + "\n");
 }
}

Matlab:

m = 999966663333;
BreakFlag = 0;

pp = primes(1000003);

num_semidivisible = 0;
sum_semidivisible = 0;

% the LPS and UPS are changing only around squares of primes. iterating on
% consecutive primes:

for n = 1:numel(pp)-1
    first = (pp(n))^2 + 1;
    last = (pp(n+1))^2 - 1;
    
    if last >= m
        last = m;
        BreakFlag = 1;
    end
    
    lps = pp(n);
    ups = pp(n+1);
    
    % the first and last are common with adj. sections, and both can't be
    % semidevisible since they are squares of their lps or ups
    
    a = ceil(first/lps);
    b = floor(last/lps);
    
    num = (b-a + 1);
    s = lps*(a+b)*num/2;
    num_semidivisible = num_semidivisible + num;
    
    a = ceil(first/ups);
    b = floor(last/ups);
 
    num = (b-a + 1);
    s = s + ups*(a+b)*num/2;
    num_semidivisible = num_semidivisible + num;

    % since LPS and UPS are primes, gcd(LPS,UPS)=1 so a number which is
    % divisble by both has to be divisible by LPS*UPS.
    % moreover, since first-1 = lps^2 and last-1 = ups^2, there's at least
    % one number which equals lps*ups, and therefore is divisible by it.
    
    a = ceil(first/ups/lps);
    b = floor(last/ups/lps);

    num = (b-a + 1);
    s = s - 2*lps*ups*(a+b)*num/2;
    num_semidivisible = num_semidivisible - 2*num;
    
    sum_semidivisible = BigNumAdd(sum_semidivisible,s);
    
    if BreakFlag
        break
    end
end