組み合わせの計算の高速化

Question

下記のリンクの問題を解いていたのですが，sampleケースが両方通って他のケースは一つがWA，残りはTLEとなってしまいました．while文も終了条件はきちんとできていますし，for文もそこまで計算量があるようには思えません．あるとしたら組み合わせを求めているところかなと思うのですが，ここはどのようにしたら高速化できますでしょうか．
問題

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

double nPm(int n,int m){
    if(m<=1)  return (double)n;
    return (n-m+1)*(nPm(n,m-1));
}

double nCm(int n,int m){
    return nPm(n,m)/nPm(m,m);
}


int main(void){
    int N,M;
    int *r;

    scanf("%d %d",&N,&M);

    r=(int *)malloc(sizeof(int)*M);

    for(int i=0;i<M;i++) scanf("%d",&r[i]);

    double p=1/(double)M;

    double P=1;

    double LOG10=N*log10(p);

    //printf("p=%f\n",p);

    int A=N;
    for(int i=0;i<M;i++){
        P=nCm(A,r[i]);
        LOG10+=log10(P);
        A-=r[i];
    }

    //printf("P=%f\n",P);

    int ans=0;
    double divide=1.0;
    while(1){
        if(LOG10>=log10(divide)) break;
        divide/=10;
        ans++;
    }

    printf("%d\n",ans);
    return 0;
}

Comments

[swordone]

間違った指摘でした。削除します。

Answer 1

swordoneさんの回答は、次のように修正すればいいということですよね。

 double nPm(int n,int m){
-    if(m<=1)  return (double)n;
-    return (n-m+1)*(nPm(n,m-1));
+    if(m<=1)  return log10(n);
+    return log10(n-m+1) + nPm(n,m-1);
 }
 
 double nCm(int n,int m){
-    return nPm(n,m)/nPm(m,m);
+    return nPm(n,m) - nPm(m,m);
 }
     for(int i=0;i<M;i++){
         P=nCm(A,r[i]);
-        LOG10+=log10(P);
+        LOG10 += P;
         A-=r[i];
     }

追記
選ばれなかった色がある場合、すなわち r[i] = 0 の場合、
double nPm(int n,int m) の m が 0 になりますよ。
次のような修正が必要なのではありませんか？

double nPm(int n,int m){
    if(m==0)  return 0;
    if(m<=1)  return log10(n);
    return log10(n-m+1) + nPm(n,m-1);
}

Comments

[kazuma-s]

お見事です。 次の 2個所を修正すればよいでしょう。

- logs=(double *)malloc(sizeof(double)*N);
+ logs = (double *)malloc(sizeof(double)*(N + 1));

- LOG10 -= -N*log10(M);
+ LOG10 -= N*log10(M);

[swordone]

これは私の回答へのコメントでしょうか？
ご指摘ありがとうございます。修正しました。

[kazuma-s]

すみません。コメントを書くところを間違えました。

Answer 2

常用対数で計算するなら、全部常用対数で計算すればいいのでは？
階乗が絡んでくると、掛け算ではちょっと数が大きくなればすぐオーバーフローしますよ。

条件を満たす投票パターンは「同じものを含む順列」の公式から
N!/(r1!r2!...rM!)通りとなるので、
これを全投票パターンのM^N通りで割れば確率pになります。
p≧10^(-x) ⇔ log10(p)≧-x ⇔ -log10(p)≦xなので、
pの常用対数を取って、-1倍して天井関数を取れば、条件を満たす最小の整数xが求められます。
pの常用対数の計算において、対数の性質から
log10(p)=log10(N!)-{log10(r1!)+log10(r2!)+...+log10(rM!)}-Nlog10(M)
であり、階乗の定義と対数の性質から自然数nに対し、
log10(n!)=log10(1)+log10(2)+...+log10(n)
なので、掛け算を使うことなく計算が可能です。
log10(N!)はループで計算することになりますが、その過程でN以下のすべての自然数nに対するlog10(n!)の値が出てくることになるので、それを配列に保存しておけば、log10(ri!)に再利用できます。
なお、0!=1に注意してください。

Cのコードは書いたことがないのですが、たぶんこんな感じ？
ポインタなども理解が十分でないため、見様見真似で書いてます。

int main(void){
    int N,M;
    double *logs;

    scanf("%d %d",&N,&M);

    logs=(double *)malloc(sizeof(double)*(N+1));
    logs[0]=logs[1]=0.0; //log10(0!)=log10(1!)=0

    double LOG10 = 0.0;

    for(int i=2;i<=N;i++){ //log10(1)=0なのでスキップして2から
        LOG10+=log10(i);
        logs[i] = LOG10;
    }

    int r;
    for(int i=0;i<M;i++) {
        scanf("%d",&r);
        LOG10 -= logs[r];
    }

    LOG10 -= N*log10(M);
    int ans = (int)ceil(-LOG10);

    printf("%d\n",ans);
    return 0;
}

Comments

[grape_ll]

確かにこの考え方なら常用対数の加法減法で表すことができるうえ，計算量も削減できるように思えます．
また，慣れない言語であるのに参考コードを書いてくださりありがとうございます．きれいに問題が解決しました．

Answer 3

nCmの値を愚直に計算すると倍精度でも表せません
C言語でこういう値がどう扱われるのかわかりませんが、おそらくそのせいでLOG10がNaNか何かになってbreakされず無限ループになってるんでしょう。

Comments

[grape_ll]

組み合わせを使用すると精度が足りないということは知りませんでした．
ですが，私の考えでは計算の過程にどうしても組み合わせを用いる必要が出てきてしまいます．なにか他のうまい方法はありますでしょうか．

[swordone]

たぶんn,mが大きくなるとnPm(n,m)もnPm(m,m)もINFINITYになると思います。
INFINITY同士の割り算だから、NaNになるのでしょう。