[C++] 프로젝트 오일러 #58 : 나선형 소수들

#58 문제도 난이도는 5% 문제입니다.

 

문제 자체는 복잡해 보여도 사실 푸는데 전혀 어려움이 없는 문제예요.

 

 

문제는 1부터 차례대로 나선형으로 수들을 배치하게 되면, 대각선쪽에는 홀수가 항상 배치되게 됩니다.  그럴 경우, 대각선에 있는 홀수 중 소수의 비율이 얼마일것인가가 이 문제를 푸는데 중요한 부분입니다.

 

이 문제는 소수의 비율이 10% 미만으로 떨어질 때, 한 변에 위치하는 수들의 갯수를 답하라는 것입니다.

 

배치하는 방법이 일정하기 때문에, 대각선에 위치한 숫자들의 순열만 파악하면 됩니다.

 

일단 첫번째 대각선(가장 작은 수)의 수열은 다음과 같이 표현할 수 있습니다.  s를 한변 길이의 절반이라고 할 때(즉, 1이 위치한 s는 0, 그 다음은 1, 그 다음은 2.. ),

 

$$f_s=4s^2-2s+1$$

 

이 것을 이용하여, 다음 수는 2s, 4s, 6s 만큼 더해주면 되겠죠.

 

이것을 단순하게 프로그램하면 아래와 같습니다.

//------------------------------------------------
//    Project Euler #58 - Spiral Primes
//        - by Aubrey Choi
//        - created at 2015-09-28
//------------------------------------------------
#include <stdio.h>
#include <time.h>
#include "isprime.h"

void solve1()
{
    long long s, pcount = 0;
    for( s = 1 ;  ; s++ )
    {
        int64_t f = 4*s*s - 2*s + 1;
        if( isPrime(f) ) pcount++;
        if( isPrime(f+2*s) ) pcount++;
        if( isPrime(f+4*s) ) pcount++;
        if( isPrime(f+6*s) ) pcount++;
        if( pcount*10 < s*4+1 ) break;
    }
    printf("Ans = %lld\n", s*2+1);
}

int main()
{
    clock_t t;

    t = clock();

    solve1();

    printf("Elapsed time is %.3f seconds \n", (float)(clock()-t)/CLOCKS_PER_SEC);
}

finding spiral primes

 

isPrime(...) 함수는 밀러-라빈 소수 판별법을 이용했습니다.  어느정도 범위의 소수에 대해서는 판별하는 숫자를 늘려놓으면 소수 검사 확률이 100%에 가깝습니다.  하지만 숫자의 범위가 커지면 소수가 아님에도 불구하고 소수로 판별할 수 있습니다.  소수를 소수가 아니라고 판정하는 경우는 없습니다.

 

#include <stdio.h>
#include <stdint.h>

bool testPrime(uint64_t n, int s, uint64_t d, int v, int a);
uint64_t MultiplyAndMod(uint64_t a, uint64_t b, uint64_t n);

//    Miller-Rabin primality test
bool isPrime(uint64_t n)
{
    int a[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 0 };
    int s = 0;
    uint64_t n1 = n-1;
    int v = 0;

    for( int i = 0 ; a[i] ; i++ )
    {
        if( a[i] == n ) return true;
        if( n%a[i] == 0 ) return false;
    }
    if( n <= 50 ) return false;

    for( s = 0 ; !(n1&1) ; s++, n1>>=1 ) ;

    for( v = 62 ; ((n1>>v)&1) == 0 ; v-- ) ;

    if( n < 2047 )
    {
        if( testPrime(n, s, n1, v, 2) == false ) return false;
        return true;
    }
    if( n < 1373653 )
    {
        if( testPrime(n, s, n1, v, 2) == false ) return false;
        if( testPrime(n, s, n1, v, 3) == false ) return false;
        return true;
    }
    if( n < 9080191 )
    {
        if( testPrime(n, s, n1, v, 31) == false ) return false;
        if( testPrime(n, s, n1, v, 73) == false ) return false;
        return true;
    }
    if( n < 25326001 )
    {
        if( testPrime(n, s, n1, v, 2) == false ) return false;
        if( testPrime(n, s, n1, v, 3) == false ) return false;
        if( testPrime(n, s, n1, v, 5) == false ) return false;
        return true;
    }
    if( n < 4759123141LL )
    {
        if( testPrime(n, s, n1, v, 2) == false ) return false;
        if( testPrime(n, s, n1, v, 7) == false ) return false;
        if( testPrime(n, s, n1, v, 61) == false ) return false;
        return true;
    }
    if( n < 1122004669633LL )
    {
        if( testPrime(n, s, n1, v, 2) == false ) return false;
        if( testPrime(n, s, n1, v, 13) == false ) return false;
        if( testPrime(n, s, n1, v, 23) == false ) return false;
        if( testPrime(n, s, n1, v, 1662803) == false ) return false;
        return true;
    }
    for( int i = 0 ; a[i] ; i++ )
        if( testPrime(n, s, n1, v, a[i]) == false ) return false;

    return true;
}

bool testPrime(uint64_t n, int s, uint64_t d, int v, int a)
{
    uint64_t ad = 1;
    uint64_t an[64];
    int i;

    for( an[0] = a, i = 1 ; i <= v ; i++ )
        an[i] = MultiplyAndMod(an[i-1], an[i-1], n);

    for( ad = 1, i = 0 ; i <= v ; i++ )
        if( (d>>i)&1 ) ad = MultiplyAndMod(ad, an[i], n);

    if( ad == 1 || ad == n-1 ) return true;

    for( i = 1 ; i < s ; i++ )
    {
        ad = MultiplyAndMod(ad, ad, n);
        if( ad == n-1 ) return true;
    }

    return false;
}

uint64_t MultiplyAndMod(uint64_t a, uint64_t b, uint64_t n)
{
    if( !(a>>32) && !(b>>32) ) return (a*b)%n;

    uint64_t x = 0, y = a;

    if( a < b ) { int64_t t = a; a = b; b = t; y = a; }

    int k = 1;
    for( int i = 2 ; i < 64 ; i++ )
        if( (b>>i)&1 ) k = i;

    for( int i = 0 ; i <= k ; i++ )
    {
        if( b&1 ) x = (x+y)%n;
        y = (y<<1)%n;
        b >>= 1;
    }
    return x;
}