NextFloor

すきな物やモノに囲まれて生きていく!

「ディビジョン・サム」

実行結果
863367
863367
717091
619112705397
638209
395123300598713973
566966
224021477247250466422784
74047
16588118325727155287303913472
263052
4363537701819179652567390088069120
449898
1963146884973045287293214903948919439360
511579
1004304720267625535004699227821413354746937344
820681
824213802133955191719266839633940485471141808308224
610913
503522926503060968017985550085287381179521128334729674752
445459
224298819317127035751542749517470882407276058817326651226980352
388807


参考

  1. codeiq.jp
  2. 約数 - Wikipedia
  3. 冪剰余 - Wikipedia
  4. 等比数列 - Wikipedia
  5. 素因数分解 - Wikipedia
  6. エラトステネスの篩 - Wikipedia


感想
問題解説を読んで、理解できた範囲で組んでみたけど、n=32までしか計算できない。
しかも計算結果が違う。n=4まで確認したけど、n=12ですでに計算が狂っている。
考えられるのは、とても大きな値を小さい値で割っているので、そこに誤差があり、
誤差があるまま乗算を繰り返しているので、計算が狂っているのかな。


ソース

/*
 入力 n (1 <= n <= 1000)
 (n!)^n の約数の和から1000003で割った余を求める
 手順
 
 ・約数の和
 ある自然数Nの素因数がpのみなら、pを0からk個選べばいいので、
 N=p^k
 (p^0 + p^1 + p^2 + ... p^k)
 
 ある自然数Nの素因数がpとqのみなら、pを0からk個、qを0からm個選べばいいので、
 N=p^k*q^m
 (p^0 + p^1 + p^2 + ... p^k)(q^0 + q^1 + q^2 + ... q^m)
 
 ある自然数Nの素因数がp0,p1,p2...pnの場合、p0を0からk0個、p1を0からk1個...pnを0からkn個選べばいいので、
 (p0^0 + ... + p0^k0)(p1^0 + ... + p1^k1)(pn^0 + ... + pn^kn)
 
 
 ・(n!)^nの素因数分解
 n=6なら、6!= 1*2*3*4*5*6
 
 1から6までの各自然数を素因数分解すると、
 1=1
 2=2^1
 3=3^1
 4=2^2
 5=5^1
 6=2^1*3^1
 
 各素因数を集めると、
  6!=1*2^4*3^2*5^1
 
 (6!)^6=(1*2^4*3^2*5^1)^6=1^6*2^24*3^12*5^6
 
 (n!)^nは、0からnまでの各自然数を素因数分解して、各素因数の出現個数を数えて、
 その各素因数の出現個数をn倍した値が、(n!)^nの各素因数の指数になる。
 
 (n!)^n=(p0^i0*p1^i1...*pk^ik)^n=p0^i0n*p1^i1n...*pk^ikn
 
 ・(n!)^nの各素因数の指数を元にして約数の和の式に当てはめて、
 (n!)^nの約数の和を求める
 ・求めた約数の和を1000003で割った余を求める

 */

#include<stdio.h>
#include<math.h>


#define N 1000
int numlist[N+1], plist = 2;

/* 探索リストから先頭の素数を取得する */
int gethedplist(void)
{
    while(numlist[plist++])
        ;
    return(plist - 1);
}

/* 探索リストから、ある素数の倍数をふるい落とす */
void eratosthenes(void)
{
    int i, j;
    int root = sqrt(N);
    while((i = gethedplist()) <= root){
        for(j = i + i; j <= N; j += i){
            numlist[j] = 1;
        }
    }
    return;
}

int isPrime(int n){ return(numlist[n]); };

/* 探索リストから素数リストを作成する エラトステネス前提 */
struct PLIST{
    int p;   //素数
    int pow; //素数の指数
}structPlist[N/2];
int plistSize;
int newPlist(void)
{
    int i;
    for(i = 2; i <= N; i++){
        if(isPrime(i) == 0){
            structPlist[plistSize].p = i;
            structPlist[plistSize].pow = 0;
            plistSize++;
        }
    }
    return(plistSize);
}
int getPlistSize(void){ return(plistSize); };

/* (n!)^nに出現する各値を素因数分解して出現する素数の個数を数える */
void factorization(int n)
{
    // n! 2,3,4,...n を素因数分解する
    int i, j, k;
    for(i = 2; i <= n; i++){
        j = i;
        k = 0;
        while(1 < j && k < getPlistSize()){
            if((j % structPlist[k].p) == 0){
                structPlist[k].pow++;
                j /= structPlist[k].p;
            }else{
                k++;
            }
        }
    }
    
    /* 
       n!で求めた各値の素因数分解の素数の指数にnを掛けることによって
       (n!)^nを素因数分解したときの素数の指数を求めることができる。
     */
    for(i = 0; i < getPlistSize(); i++){
        structPlist[i].pow *= n;
    }
    return;
}

/* 約数関数から約数の和を(等比数列)10^6+3で割った余を求める */
/*
    S1S2S3...Sn mod M
    ((S1 mod M)(S2 mod M)(S3 mod M)...(Sn mod M)) mod M
 
    S1 = (P1^(n+1) - 1)/(P1 - 1)
    P1 = 2
 */
long double ans2(void)
{
    long double ans2 = 1, mod = 1000003, e=1, p=1, t;
    int i;
    for(i = 0; i < getPlistSize(); i++){
        if(0 < structPlist[i].pow){
            p = pow(structPlist[i].p, structPlist[i].pow + 1) - 1;
            e = structPlist[i].p - 1;
            t = fmodl(p/e, mod);
            printf("%.0Lf\n", t);
            ans2 *= t;
            printf("%.0Lf\n", ans2);
        }
    }
    return(fmodl(ans2, mod));
}

void printtest()
{
    int i;
    for(i = 2; i <= 1000; i++){
        if(numlist[i] == 0){
            printf("%d,", i);
        }
    }
}

void printn(void)
{
    int i;
    printf("\n");
    for(i = 0; i < getPlistSize(); i++){
        printf("%d::%d\n",structPlist[i].p, structPlist[i].pow);
    }
}


int main(int argc, const char * argv[])
{
    eratosthenes();
    newPlist();
    factorization(32);
    printf("%.0Lf\n", ans2());
    
    return 0;
}