목표는 간단하다 : 계산 totient 기능을 에서 당신이 할 수 많은 번호 등을 위해 십초 와 숫자의 합계.
마지막에 결과를 인쇄해야하며 실제로 계산해야합니다. 자동 강의 기능은 허용되지 않지만 bignum 라이브러리는 허용됩니다. 1부터 시작하여 모든 정수를 연속적으로 계산해야합니다. 당신은되어 있지 숫자를 건너 뛸 수있었습니다.
귀하의 점수는 귀하의 프로그램이 귀하의 컴퓨터에서 계산할 수있는 숫자 / 내 프로그램이 귀하의 컴퓨터에서 계산할 수있는 숫자입니다 . 내 코드는 C ++ (최적화 해제)의 간단한 프로그램이므로 실행할 수 있기를 바랍니다.
사용할 수있는 중요한 속성!
- 만약
gcd(m,n) = 1, phi(mn) = phi(m) * phi(n) - 경우
p소수,phi(p) = p - 1(에 대한p < 10^20) - 경우
n도있다phi(2n) = 2 phi(n) - 첫 번째 링크에 나열된 기타
내 코드
#include <iostream>
using namespace std;
int gcd(int a, int b)
{
while (b != 0)
{
int c = a % b;
a = b;
b = c;
}
return a;
}
int phi(int n)
{
int x = 0;
for (int i=1; i<=n; i++)
{
if (gcd(n, i) == 1)
x++;
}
return x;
}
int main()
{
unsigned int sum = 0;
for (int i=1; i<19000; i++) // Change this so it runs in 10 seconds
{
sum += phi(i);
}
cout << sum << endl;
return 0;
}답변
님로드 : ~ 38,667 (580,000,000 / 15,000)
이 답변은 매우 간단한 접근법을 사용합니다. 이 코드는 복합 숫자 (프라임의 경우 0)에 대해 각 슬롯에서 가장 작은 소수의 소수를 저장 한 다음 소수 범위의 동적 함수를 사용하여 동일한 범위에 대한 배신자 기능을 구성한 다음 결과를 합산하는 간단한 소수 체를 사용합니다. 이 프로그램은 체를 구성하는 데 거의 모든 시간을 소비 한 다음, 짧은 시간 안에 계류 기능을 계산합니다. 효율적인 체를 구성하는 것으로 보입니다 (결과에서 복합 수의 주요 요소를 추출 할 수 있어야하고 메모리 사용을 합리적인 수준으로 유지 해야하는 약간의 왜곡이 있음).
업데이트 : 메모리 공간을 줄이고 캐시 동작을 개선하여 성능을 향상시킵니다. 성능을 5 % -10 % 향상시킬 수는 있지만 코드 복잡성의 증가는 그만한 가치가 없습니다. 궁극적으로이 알고리즘은 주로 CPU의 von Neumann 병목 현상을 발생 시키며,이를 해결할 수있는 알고리즘 조정은 거의 없습니다.
또한 C ++ 코드가 모든 최적화와 함께 컴파일되도록 설계되지 않았으므로 다른 누구도 수행하지 않았기 때문에 성능 분할기가 업데이트되었습니다. 🙂
업데이트 2 : 메모리 액세스 개선을위한 최적화 된 체 작업. memcpy () (~ 5 % speedup)를 통해 작은 소수를 일괄 처리하고 더 큰 소수를 체질 할 때 2, 3, 5의 배수를 건너 뜁니다 (~ 10 % 속도 향상).
C ++ 코드 : 9.9 초 (g ++ 4.9 포함)
님로드 코드 : 9.9 초 (-d : release, gcc 4.9 백엔드 사용)
proc handleSmallPrimes(sieve: var openarray[int32], m: int) =
# Small primes are handled as a special case through what is ideally
# the system's highly optimized memcpy() routine.
let k = 2*3*5*7*11*13*17
var sp = newSeq[int32](k div 2)
for i in [3,5,7,11,13,17]:
for j in countup(i, k, 2*i):
sp[j div 2] = int32(i)
for i in countup(0, sieve.high, len(sp)):
if i + len(sp) <= len(sieve):
copyMem(addr(sieve[i]), addr(sp[0]), sizeof(int32)*len(sp))
else:
copyMem(addr(sieve[i]), addr(sp[0]), sizeof(int32)*(len(sieve)-i))
# Fixing up the numbers for values that are actually prime.
for i in [3,5,7,11,13,17]:
sieve[i div 2] = 0
proc constructSieve(m: int): seq[int32] =
result = newSeq[int32](m div 2 + 1)
handleSmallPrimes(result, m)
var i = 19
# Having handled small primes, we only consider candidates for
# composite numbers that are relatively prime with 31. This cuts
# their number almost in half.
let steps = [ 1, 7, 11, 13, 17, 19, 23, 29, 31 ]
var isteps: array[8, int]
while i * i <= m:
if result[i div 2] == 0:
for j in 0..7: isteps[j] = i*(steps[j+1]-steps[j])
var k = 1 # second entry in "steps mod 30" list.
var j = 7*i
while j <= m:
result[j div 2] = int32(i)
j += isteps[k]
k = (k + 1) and 7 # "mod 30" list has eight elements.
i += 2
proc calculateAndSumTotients(sieve: var openarray[int32], n: int): int =
result = 1
for i in 2'i32..int32(n):
var tot: int32
if (i and 1) == 0:
var m = i div 2
var pp: int32 = 2
while (m and 1) == 0:
pp *= 2
m = m div 2
if m == 1:
tot = pp div 2
else:
tot = (pp div 2) * sieve[m div 2]
elif sieve[i div 2] == 0: # prime?
tot = i - 1
sieve[i div 2] = tot
else:
# find and extract the first prime power pp.
# It's relatively prime with i/pp.
var p = sieve[i div 2]
var m = i div p
var pp = p
while m mod p == 0 and m != p:
pp *= p
m = m div p
if m == p: # is i a prime power?
tot = pp*(p-1)
else:
tot = sieve[pp div 2] * sieve[m div 2]
sieve[i div 2] = tot
result += tot
proc main(n: int) =
var sieve = constructSieve(n)
let totSum = calculateAndSumTotients(sieve, n)
echo totSum
main(580_000_000)
답변
Java, ~ 24,000 점 (360,000,000 / 15,000)
아래의 자바 코드는 totient 함수와 주요 체를 함께 계산합니다. 컴퓨터에 따라 초기 / 최대 힙 크기를 늘려야합니다 (느린 노트북의 경우 최대 올라갔습니다 -Xmx3g -Xms3g).
public class Totient {
final static int size = 360000000;
final static int[] phi = new int[size];
public static void main(String[] args) {
long time = System.currentTimeMillis();
long sum = 0;
phi[1] = 1;
for (int i = 2; i < size; i++) {
if (phi[i] == 0) {
phi[i] = i - 1;
for (int j = 2; i * j < size; j++) {
if (phi[j] == 0)
continue;
int q = j;
int f = i - 1;
while (q % i == 0) {
f *= i;
q /= i;
}
phi[i * j] = f * phi[q];
}
}
sum += phi[i];
}
System.out.println(System.currentTimeMillis() - time);
System.out.println(sum);
}
}답변
님로드 : ~ 2,333,333 (42,000,000,000 / 18,000)
이것은 이전 답변과 완전히 다른 접근법을 사용합니다. 자세한 내용은 의견을 참조하십시오. longint모듈을 찾을 수 있습니다 여기에 .
import longint
const max = 500_000_000
var ts_mem: array[1..max, int]
# ts(n, d) is defined as the number of pairs (a,b)
# such that 1 <= a <= b <= n and gcd(a,b) = d.
#
# The following equations hold:
#
# ts(n, d) = ts(n div d, 1)
# sum for i in 1..n of ts(n, i) = n*(n+1)/2
#
# This leads to the recurrence:
# ts(n, 1) = n*(n+1)/2 - sum for i in 2..n of ts(n, i)
#
# or, where ts(n) = ts(n, 1):
# ts(n) = n*(n+1)/2 - sum for i in 2..n of ts(n div i)
#
# Note that the large numbers that we deal with can
# overflow 64-bit integers.
proc ts(n, gcd: int): int =
if n == 0:
result = 0
elif n == 1 and gcd == 1:
result = 1
elif gcd == 1:
result = n*(n+1) div 2
for i in 2..n:
result -= ts(n, i)
else:
result = ts(n div gcd, 1)
# Below is the optimized version of the same algorithm.
proc ts(n: int): int =
if n == 0:
result = 0
elif n == 1:
result = 1
else:
if n <= max and ts_mem[n] > 0:
return ts_mem[n]
result = n*(n+1) div 2
var p = n
var k = 2
while k < n div k:
let pold = p
p = n div k
k += 1
let t = ts(n div pold)
result -= t * (pold-p)
while p >= 2:
result -= ts(n div p)
p -= 1
if n <= max:
ts_mem[n] = result
proc ts(n: int128): int128 =
if n <= 2_000_000_000:
result = ts(n.toInt)
else:
result = n*(n+1) div 2
var p = n
var k = 2
while k < n div k:
let pold = p
p = n div k
k += 1
let t = ts(n div pold)
result = result - t * (pold-p)
while p >= 2:
result = result - ts(n div p)
p = p - 1
echo ts(42_000_000_000.toInt128)
답변
C # : 49,000 (980,000,000 / 20,000)
/codegolf//a/26800 “하워드 코드”.
그러나 수정 된 phi 값은 홀수 정수에 대해 계산됩니다.
using System;
using sw = System.Diagnostics.Stopwatch;
class Program
{
static void Main()
{
sw sw = sw.StartNew();
Console.Write(sumPhi(980000000) + " " + sw.Elapsed);
sw.Stop(); Console.Read();
}
static long sumPhi(int n) // sum phi[i] , 1 <= i <= n
{
long s = 0; int[] phi;
if (n < 1) return 0; phi = buildPhi(n + 1);
for (int i = 1; i <= n; i++) s += getPhi(i, phi);
return s;
}
static int getPhi(int i, int[] phi)
{
if ((i & 1) > 0) return phi[i >> 1];
if ((i & 3) > 0) return phi[i >> 2];
int z = ntz(i); return phi[i >> z >> 1] << z - 1;
}
static int[] buildPhi(int n) // phi[i >> 1] , i odd , i < n
{
int i, j, y, x, q, r, f; int[] phi;
if (n < 2) return new int[] { 0 };
phi = new int[n / 2]; phi[0] = 1;
for (j = 2, i = 3; i < n; i *= 3, j *= 3) phi[i >> 1] = j;
for (x = 4, i = 5; i <= n >> 1; i += x ^= 6)
{
if (phi[i >> 1] > 0) continue; phi[i >> 1] = i ^ 1;
for (j = 3, y = 3 * i; y < n; y += i << 1, j += 2)
{
if (phi[j >> 1] == 0) continue; q = j; f = i ^ 1;
while ((r = q) == i * (q /= i)) f *= i;
phi[y >> 1] = f * phi[r >> 1];
}
}
for (; i < n; i += x ^= 6) // primes > n / 2
if (phi[i >> 1] == 0)
phi[i >> 1] = i ^ 1;
return phi;
}
static int ntz(int i) // number of trailing zeros
{
int z = 1;
if ((i & 0xffff) == 0) { z += 16; i >>= 16; }
if ((i & 0x00ff) == 0) { z += 08; i >>= 08; }
if ((i & 0x000f) == 0) { z += 04; i >>= 04; }
if ((i & 0x0003) == 0) { z += 02; i >>= 02; }
return z - (i & 1);
}
}새로운 점수 : 61,000 (1,220,000,000 / 20,000)
“App.config”에서 “gcAllowVeryLargeObjects enabled = true”를 추가해야했습니다.
static long sumPhi(int n)
{
int i1, i2, i3, i4, z; long s1, s2, s3, s4; int[] phi;
if (n < 1) return 0; phi = buildPhi(n + 1); n -= 4; z = 2;
i1 = 1; i2 = 2; i3 = 3; i4 = 4; s1 = s2 = s3 = s4 = 0;
if (n > 0)
for (; ; )
{
s1 += phi[i1 >> 1];
s2 += phi[i2 >> 2];
s3 += phi[i3 >> 1];
s4 += phi[i4 >> z >> 1] << z - 1;
i1 += 4; i2 += 4; i3 += 4; i4 += 4;
n -= 4; if (n < 0) break;
if (z == 2)
{
z = 3; i4 >>= 3;
while ((i4 & 3) == 0) { i4 >>= 2; z += 2; }
z += i4 & 1 ^ 1;
i4 = i3 + 1;
}
else z = 2;
}
if (n > -4) s1 += phi[i1 >> 1];
if (n > -3) s2 += phi[i2 >> 2];
if (n > -2) s3 += phi[i3 >> 1];
if (n > -1) s4 += phi[i4 >> z >> 1] << z - 1;
return s1 + s2 + s3 + s4;
}
static int[] buildPhi(int n)
{
int i, j, y, x, q0, q1, f; int[] phi;
if (n < 2) return new int[] { 0 };
phi = new int[n / 2]; phi[0] = 1;
for (uint u = 2, v = 3; v < n; v *= 3, u *= 3) phi[v >> 1] = (int)u;
for (x = 4, i = 5; i <= n >> 1; i += x ^= 6)
{
if (phi[i >> 1] > 0) continue; phi[i >> 1] = i ^ 1;
for (j = 3, y = 3 * i; y < n; y += i << 1, j += 2)
{
if (phi[j >> 1] == 0) continue; q0 = j; f = i ^ 1;
while ((q1 = q0) == i * (q0 /= i)) f *= i;
phi[y >> 1] = f * phi[q1 >> 1];
}
}
for (; i < n; i += x ^= 6)
if (phi[i >> 1] == 0)
phi[i >> 1] = i ^ 1;
return phi;
}답변
파이썬 3 : ~ 24000 (335,000,000 / 14,000)
내 버전은 Howard 알고리즘의 Python 포트입니다 . 내 원래 기능은이 블로그 게시물에 소개 된 알고리즘을 수정 한 것입니다 .
Numpy 및 Numba 모듈을 사용하여 실행 시간을 단축하고 있습니다. 일반적으로 Numba를 사용할 때 로컬 변수의 유형을 선언 할 필요는 없지만이 경우 가능한 한 실행 시간을 짜고 싶었습니다.
편집 : constructsieve 및 summarum 함수를 단일 함수로 결합했습니다.
C ++ : 9.99 초 (n = 14,000); 파이썬 3 : 9.94 (n = 335,000,000)
import numba as nb
import numpy as np
import time
n = 335000000
@nb.njit("i8(i4[:])", locals=dict(
n=nb.int32, s=nb.int64, i=nb.int32,
j=nb.int32, q=nb.int32, f=nb.int32))
def summarum(phi):
s = 0
phi[1] = 1
i = 2
while i < n:
if phi[i] == 0:
phi[i] = i - 1
j = 2
while j * i < n:
if phi[j] != 0:
q = j
f = i - 1
while q % i == 0:
f *= i
q //= i
phi[i * j] = f * phi[q]
j += 1
s += phi[i]
i += 1
return s
if __name__ == "__main__":
s1 = time.time()
a = summarum(np.zeros(n, np.int32))
s2 = time.time()
print(a)
print("{}s".format(s2 - s1))답변
다음은 10 초 동안 ~ 60000 개의 숫자를 크랭크 할 수있는 파이썬 구현입니다. pollard rho 알고리즘을 사용하여 숫자를 분해하고 Rabin miller primality 테스트를 사용하고 있습니다.
from Queue import Queue
import random
def gcd ( a , b ):
while b != 0: a, b = b, a % b
return a
def rabin_miller(p):
if(p<2): return False
if(p!=2 and p%2==0): return False
s=p-1
while(s%2==0): s>>=1
for _ in xrange(10):
a=random.randrange(p-1)+1
temp=s
mod=pow(a,temp,p)
while(temp!=p-1 and mod!=1 and mod!=p-1):
mod=(mod*mod)%p
temp=temp*2
if(mod!=p-1 and temp%2==0): return False
return True
def pollard_rho(n):
if(n%2==0): return 2;
x=random.randrange(2,1000000)
c=random.randrange(2,1000000)
y=x
d=1
while(d==1):
x=(x*x+c)%n
y=(y*y+c)%n
y=(y*y+c)%n
d=gcd(x-y,n)
if(d==n): break;
return d;
def primeFactorization(n):
if n <= 0: raise ValueError("Fucked up input, n <= 0")
elif n == 1: return []
queue = Queue()
factors=[]
queue.put(n)
while(not queue.empty()):
l=queue.get()
if(rabin_miller(l)):
factors.append(l)
continue
d=pollard_rho(l)
if(d==l):queue.put(l)
else:
queue.put(d)
queue.put(l/d)
return factors
def phi(n):
if rabin_miller(n): return n-1
phi = n
for p in set(primeFactorization(n)):
phi -= (phi/p)
return phi
if __name__ == '__main__':
n = 1
s = 0
while n < 60000:
n += 1
s += phi(n)
print(s)
답변
φ (2 n ) = 2 n − 1
Σ φ (2 i ) = 2 i − 1 ~ 1에서 n까지
먼저, 시간을 찾을 무언가 :
import os
from time import perf_counter
SEARCH_LOWER = -1
SEARCH_HIGHER = 1
def integer_binary_search(start, lower=None, upper=None, big_jump=1):
if lower is not None and lower == upper:
raise StopIteration # ?
result = yield start
if result == SEARCH_LOWER:
if lower is None:
yield from integer_binary_search(
start=start - big_jump,
lower=None,
upper=start - 1,
big_jump=big_jump * 2)
else:
yield from integer_binary_search(
start=(lower + start) // 2,
lower=lower,
upper=start - 1)
elif result == SEARCH_HIGHER:
if upper is None:
yield from integer_binary_search(
start=start + big_jump,
lower=start + 1,
upper=None,
big_jump=big_jump * 2)
else:
yield from integer_binary_search(
start=(start + upper) // 2,
lower=start + 1,
upper=upper)
else:
raise ValueError('Expected SEARCH_LOWER or SEARCH_HIGHER.')
search = integer_binary_search(start=1000, lower=1, upper=None, big_jump=2500)
n = search.send(None)
while True:
print('Trying with %d iterations.' % (n,))
os.spawnlp(
os.P_WAIT,
'g++', 'g++', '-Wall', '-Wextra', '-pedantic', '-O0', '-o', 'reference',
'-DITERATIONS=%d' % (n,),
'reference.cpp')
start = perf_counter()
os.spawnl(os.P_WAIT, './reference', './reference')
end = perf_counter()
t = end - start
if t >= 10.1:
n = search.send(SEARCH_LOWER)
elif t <= 9.9:
n = search.send(SEARCH_HIGHER)
else:
print('%d iterations in %f seconds!' % (n, t))
break
참조 코드의 경우 다음과 같습니다.
…
14593 반복 시도.
9.987747
초에 64724364 14593 반복!
이제 Haskell :
import System.Environment (getArgs)
phiSum :: Integer -> Integer
phiSum n = 2 ^ n - 1
main :: IO ()
main = getArgs >>= print . phiSum . (2^) . read . head
0.718 초 안에 2525224 자리의 무언가를 만듭니다. 그리고 지금은 @Howard의 의견을 주목하고 있습니다.