プロジェクトオイラー Problem 240 「上位のサイコロ」 †

Problem 240 「上位のサイコロ」 †

6面のサイコロ(各面は 1 から 6)を 5 個振って, 上位 3 個の合計が 15 となる場合は 1111 通りある. いくつか例を挙げる:

D1,D2,D3,D4,D5 = 4,3,6,3,5
D1,D2,D3,D4,D5 = 4,3,3,5,6
D1,D2,D3,D4,D5 = 3,3,3,6,6
D1,D2,D3,D4,D5 = 6,6,3,3,3

12面のサイコロ(各面は 1 から 12)を 20 個振って, 上位 10 個の合計が 70 となる場合は何通りあるか.

http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20240


解法
数字の大きなサイコロから探索で10個試します。
その10個の中の一番小さなサイと同じものが何個あるかとそれより小さなサイコロが何個あるかで組み合わせを計算するだけです。
20！を分子として同じ目のサイコロが何個あるかで分母の割り算が決定されるだけです。
書くのは疲れましたが整理して考えれば簡単な問題。

もっとさいころの数が増えて高速化することを考えた場合、数千桁を普通に扱える言語で動的計画法が必要になるかと思います。

コンパイラー　ＢＣＣ5.5
言語Ｃ++
64ビット整数型でもギリギリの答えですが計算は一秒もかかりません。
コード製作者　堀江伸一


#include<stdio.h>
#include<queue>
#include<iostream>
#include<map>
#include<math.h>

const int XAI_LIMIT=12;//n面ダイス
const int XAI_TOP_10 =10;//上位10個のサイコロの数
const int XAI_COUNT_LIMIT=20;
const int ANS_SUM=70;//答えとなる数

struct S{
      int sum;
      unsigned __int64 allXaiCount,nowXaiCount,nowDown,div;
      bool operator<(const S& s)const{
            if(sum!=s.sum)return sum<s.sum;
            if(allXaiCount!=s.allXaiCount)return allXaiCount<s.allXaiCount;
            if(nowXaiCount!=s.nowXaiCount)return nowXaiCount<s.nowXaiCount;
            if(nowDown!=s.nowDown)return nowDown<s.nowDown;
            return div<s.div;
      }
};

unsigned __int64 fact(unsigned __int64 n){
      unsigned __int64 result=1;
      while(n>0){
            result*=n;
            n--;
      }
      return result;
}
unsigned __int64 nCr(unsigned __int64 n,unsigned __int64 r){
      unsigned __int64 n1=fact(n),div1=fact(r),div2=fact(n-r);
      return (n1/div1)/div2;
}
unsigned __int64 calc_perm(S s1){
      unsigned __int64 all=fact(XAI_COUNT_LIMIT);
      unsigned __int64 div,perm;
      unsigned __int64 result=0;
      unsigned __int64 down=s1.nowDown-1;
      int start=0;
      if(down==0)start=XAI_COUNT_LIMIT-s1.allXaiCount;
      for(int i=start;i+s1.allXaiCount<=XAI_COUNT_LIMIT;i++){
            perm=all/(fact(i+s1.nowXaiCount)*fact(XAI_COUNT_LIMIT-XAI_TOP_10 -i));
            result+=perm/s1.div*(down==0?1:pow(down,XAI_COUNT_LIMIT-i-s1.allXaiCount));          
      }
      return result;
}

int main(){
      unsigned __int64 all;
      all=fact(XAI_TOP_10 );
      std::map<S,unsigned __int64> dp,nextDP;
      std::map<S,unsigned __int64>::iterator it;
      S s,s2;
      s.sum=0;
      s.allXaiCount=s.nowXaiCount=s.nowDown=0;
      s.div=1;
      dp[s]=1;
      unsigned __int64 ans=0;
      for(unsigned __int64 i=XAI_LIMIT;i>=1;i--){
            nextDP.clear();
            for(it=dp.begin();it!=dp.end();it++){
                  s=(*it).first;
                  for(int j=0;j<=XAI_TOP_10 ;j++){
                        s2=s;
                        s2.allXaiCount+=j;
                        if(s2.allXaiCount>XAI_TOP_10 )break;
                        s2.nowXaiCount=j;
                        s2.sum+=j*i;
                        if(s2.sum>ANS_SUM)break;
                     
                        if(j>0)s2.nowDown=i;
                     
                        if(s2.sum==ANS_SUM && s2.allXaiCount==XAI_TOP_10 &&j>0){
                              ans+=calc_perm(s2)*(*it).second;
                        }else{
                              s2.div*=fact(j);
                              if(nextDP.find(s2)==nextDP.end())nextDP[s2]=0;
                              nextDP[s2]+=(*it).second;
                        }
                  }
            }
            dp.clear();
            dp.insert(nextDP.begin(),nextDP.end());
      }
      std::cout<<ans;
}