プロジェクトオイラー Problem 237 「4 × n のゲーム盤上を進む順路」 †

Problem 237 「4 × n のゲーム盤上を進む順路」 †

T(n) を以下のルールに従い 4 × n のゲーム盤上を進む順路の数と定義する:

• 左上の角から始める
• 1マス分の上下左右の移動を繰り返す
• 各マスを全てちょうど1回ずつ通る
• 左下の角で終わる

下の図は 4 × 10 の盤上の順路の一例である:

T(10) は 2329 である. T(1012) を 108 で割った余りを求めよ.

解法
一列ずつの動的計画法でまず計算をしてみました。
次に1列を2つつなげて2列
2列を二つつなげて4列
4列を二つつなげて8列、、、
としてあとは１，２，4,8、、、列を結合して10＾12-２列を作ります。
最初の列とそれをつなげて最後の一列は実行結果を目視で確認して足し算。
最新バージョンのＳＷＩＰｒｏｌｏｇでは10＾12や10＾8でエラーが出るようです。

integer(10^12)-2などのようにして囲む必要があるようです。
おそらく高速化のためだと思いますが個人的には不便になっただけな気がします。
数千桁や数万桁をシームレスにあつかえてこそのＰｒｏｌｏｇだと思うのですが。



記述言語
Ｐｒｏｌｏｇ


datas(X):-
      X=[[[1,2,3,4],[1,2,3,4],1],
         [[1,4,3,2],[1,4,3,2],1],
         [[3,2,1,4],[3,2,1,4],1],
         [[1,0,0,2],[1,2,3,4],1],
         [[1,0,0,2],[1,2,0,0],1],
         [[1,4,3,2],[1,0,0,2],1],
         [[3,2,1,4],[1,0,0,2],1],
         [[1,2,0,0],[1,0,0,2],1],
         [[1,0,0,2],[0,1,2,0],1],
         [[1,0,0,2],[0,0,1,2],1],
         [[1,0,2,0],[0,1,0,2],1],
         [[0,1,2,0],[1,0,0,2],1],
         [[0,1,0,2],[1,0,2,0],1],
         [[0,0,1,2],[1,0,0,2],1],
         [[1,2,0,0],[1,4,3,2],1],
         [[0,0,1,2],[3,2,1,4],1],
         [[1,2,3,4],[0,0,1,2],1],
         [[1,2,3,4],[1,2,0,0],1],
         [[1,4,3,2],[1,0,0,2],1]].

union_sumB([],[E,E1,Count],[[E,E1,Count]]):-!.
union_sumB([[E,E1,Count]|Rest],[E,E1,Count1],Result):-
      !,
      Count2 is (Count+Count1) mod 10^8,
      union_sumB(Rest,[E,E1,Count2],Result).
union_sumB([E|Rest],E1,[E1|Result]):-
      !,
      union_sumB(Rest,E,Result).


union_sum([],[E,Count],[[E,Count]]):-!.
union_sum([[E,Count]|Rest],[E,Count1],Result):-
      !,
      Count2 is (Count+Count1) mod 10^8,
      union_sum(Rest,[E,Count2],Result).
union_sum([E|Rest],E1,[E1|Result]):-
      !,
      union_sum(Rest,E,Result).

next_calc(Base,Datas,[E1,Count3]):-
      member([E,Count1],Base),
      member([E,E1,Count2],Datas),
      Count3 is Count1*Count2.

next_calc_w(Base,Base3,Datas):-
      findall(E,next_calc(Base,Datas,E),Base1),
      msort(Base1,[Top|Base2]),
      union_sum(Base2,Top,Base3).



next_double(Datas,[E,E2,Count3]):-
      member([E1,E2,Count1],Datas),
      member([E,E1,Count2],Datas),
      Count3 is (Count1*Count2) mod 10^8.

next_double_w(Datas,Datas3):-
      findall(E,next_double(Datas,E),Datas1),
      msort(Datas1,[Top|Datas2]),
      union_sumB(Datas2,Top,Datas3).



dp(0,Base,_):-
      !,
      member([[1,0,0,2],C1],Base),
      member([[1,2,3,4],C2],Base),
      Ans is (C1+C2) mod 10^8,
      write([ans,Ans]).

dp(R,Base,Datas):-
      R mod 2=:=1,
      !,
      R1 is R//2,
      next_calc_w(Base,Base1,Datas),
      next_double_w(Datas,Datas1),
      dp(R1,Base1,Datas1).
dp(R,Base,Datas):-
      !,
      next_double_w(Datas,Datas1),
      R1 is R//2,
      dp(R1,Base,Datas1).

main237:-
      datas(X),
      sort(X,Datas),
      Seed=[[[1,2,0,0],1],[[1,2,3,4],1],
            [[0,1,2,0],1],[[0,0,1,2],1]],
      R is 10,
      R1 is R^12-2,
      dp(R1,Seed,Datas).

プロジェクトオイラー Problem 79 「パスコードの導出」 †

Problem 79 「パスコードの導出」 †

オンラインバンクで通常使われるsecurity methodは, パスコードからランダムに選んだ3文字をユーザーに要求するものである.

たとえば, パスコードが531278のとき, 2番目, 3番目, 5番目の文字を要求されるかもしれない. このとき, 期待される答えは: 317 である.

テキストファイルkeylog.txtには, ログインに成功した50回の試行が記録されている.

3つの文字が常に順番通りに要求されるとするとき, ファイルを分析して, 可能なパスコードのなかでもっとも短いものを見つけよ.

http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2079
詳細はリンク先で。





解法
トポロジカルソートで片が付きます。
トポロジカルでない場合答えが複数になるのでトポロジカルになるということですかね。

e(E1,E2,_,E1,E2).
e(_,E2,E3,E2,E3):-!.

e1(Top,_,Top).
e1(_,Tail,Tail).

change_G(Datas,[Top,Tail]):-
      member([E1,E2,E3],Datas),
      e(E1,E2,E3,Top,Tail).


dells(G,Top,[Top1,Tail1]):-
      member([Top1,Tail1],G),
      not(Top=:=Top1).

search_next([],[Num]):-
      !,
      Num1 is Num-48,
      write(Num1).

search_next(G,Nums):-
      select(Top,Nums,Nums1),
      not(member([_,Top],G)),
      !,
      findall(E,dells(G,Top,E),G1),
      Top1 is Top-48,
      write(Top1),
      search_next(G1,Nums1).

numbers1(G,E):-
      member([Top,Tail],G),
      e1(Top,Tail,E).
numbers(G,Nums2):-
      findall(E,numbers1(G,E),Nums1),
      sort(Nums1,Nums2).

main79:-
      open('pe79.txt',read,IS),
      read_term(IS,Datas,[]),
      close(IS),
      findall(E,change_G(Datas,E),G),
      numbers(G,Nums),
      search_next(G,Nums).

プロジェクトオイラー 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;
}

プロジェクトオイラー Problem 244 「スライダー」 †

Problem 244 「スライダー」 †

おそらく15パズルはご存知だろう. ここでは, 数字の書かれたタイルの代わりに, 7枚の赤いタイルと8枚の青いタイルを使用する.

移動は, タイルの動いた方向(Left, Right, Up, Down)の大文字のイニシャルで表す. すなわち, 最初の状態 (S) から文字列 LULUR を経て状態 (E) にたどり着く.

(S)  , (E) 

各経路は, 以下に示す擬似コードによってチェックサムが計算される:

checksum = 0
checksum = (checksum × 243 + m1) mod 100 000 007
checksum = (checksum × 243 + m2) mod 100 000 007
…
checksum = (checksum × 243 + mn) mod 100 000 007

mk は移動の文字列の k 番目の文字のアスキーコードの値である. アスキーコードは以下の通り:

L	76
R	82
U	85
D	68

上で例に挙げた文字列 LULUR の場合, チェックサムは 19761398 になるだろう.

では, 状態 (S) から始めて, 状態 (T) に到達する最短の経路を全て求めよ.

(S)  , (T) 

最小の長さとなる経路のチェックサム全ての合計を求めよ.

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

解法
盤面を一マス2ビットで表現します。
隙間　11
赤　10
青　01
右上から左氏へビットであらわしてあとは同じ手数で同じ盤面のものを統合しながら動的計画法でアセプトです。
最短手数を残すたびに一手先を計算するたびに一手前の盤面と同じになったらそれを削除します。
そして解の手が見つかったらその時点で計算を終了します。
コードのほとんどの処理がＤＰが苦手なＰｒｏｌｏｇ言語にＤＰを行わせるための処理になっており計算本体はかなりコンパクトにまとめたつもりです。


コード実行結果
% 6,090,943 inferences, 1.469 CPU in 1.488 seconds (99% CPU, 4147025 Lips)
出てきた答えのＬｉｓｔ2つ目が答え。


modN(100000007):-!.
seed([[2779096487,0,1]]):-!.
goal([1721329307,_,_]):-!.


move_next_pos(P1,[P2,Mk]):-
      member([D,Mk],[[-1,82],[1,76]]),
      P2 is P1+D,
      0=<P2,
      P2=<15,
      P1//4=:=P2//4.
move_next_pos(P1,[P2,Mk]):-
      member([D,Mk],[[4,85],[-4,68]]),
      P2 is P1+D,
      0=<P2,
      P2=<15.

dells([],List,List):-!.
dells(_,[],[]):-!.
dells([[E,_,_]|Rest],[[E,_,_]|Rest1],Result):-
      !,
      dells(Rest,Rest1,Result).
dells([[E,S1,C1]|Rest],[[E1,CheckSum,Count]|Rest1],[[E1,CheckSum,Count]|Result]):-
      E>E1,
      !,
      dells([[E,S1,C1]|Rest],Rest1,Result).
dells([_|Rest],Es,Result):-
      !,
      dells(Rest,Es,Result).

calc_next(Boards,[E1,CheckSum1,Count]):-
      member([E,CheckSum,Count],Boards),
      between(0,15,P1),
      Bit1 is 3<<(P1*2),
      (E/\Bit1)=:=Bit1,
      move_next_pos(P1,[P2,Mk]),
      Bit2 is 3<<(P2*2),
      Bit3 is (E /\ (3<<(P2*2))),
      Bit4 is (Bit3>>(P2*2))<<(P1*2),
      E1 is E-Bit1-Bit3+Bit2+Bit4,
      modN(Mod),
      CheckSum1 is (CheckSum*243+Mk*Count) mod Mod.

union_sum([],E,[E]):-!.
union_sum([[E,CheckSum,Count]|Rest],[E,CheckSum1,Count1],Result):-
      !,
        modN(Mod),
      CheckSum2 is (CheckSum+CheckSum1 ) mod Mod,
      Count2 is Count+Count1,
      union_sum(Rest,[E,CheckSum2,Count2],Result).
union_sum([[E,CheckSum,Count]|Rest],E1,[E1|Result]):-
      !,
      union_sum(Rest,[E,CheckSum,Count],Result).

dp(_,NowList):-
      goal(X),
      member(X,NowList),
      !,
      write(X).
dp(OldList,NowList):-
      write(s),nl,
      findall(E,calc_next(NowList,E),NextList1),
      msort(NextList1,[Top|NextList2]),
      union_sum(NextList2,Top,NextList3),
      dells(OldList,NextList3,NextList4),
      dp(NowList,NextList4).


main244:-
      seed(Seed),
      dp([],Seed).

Problem 78　「コインの分割」 †

n 枚のコインを異なった方法で山に分ける場合の数を p(n) と表わす. 例えば, 5枚のコインを山に分ける異なったやり方は７通りなので p(5)=7 となる.

OOOOO

OOOO   O

OOO   OO

OOO   O   O

OO   OO   O

OO   O   O   O

O   O   O   O   O

p(n) が100万で割り切れる場合に最小となる n を求めよ.

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

解法
問76のコードを少しだけ改良して流用します。
計算はmod100万と取った時0になればよいので合同式の中で組み合わせを計算します。

#include<stdio.h>
#include<iostream>
#include<vector>

const int MOD=1000*1000;

int pos(int n){
      n++;
      int n1=n/2*(n%2==0?1:-1);
      return (n1*(3*n1-1))/2;
}


int main(){
      std::vector<int> dp;
      dp.push_back(1);
      int num=1;
      while(1){
            dp.push_back(0);
            for(int j=1;j<=num;j++){
                  int p=pos(j);
                  if(num-p<0)break;
                  __int64 mult=((j-1)%4)/2==0?1:-1;
                  dp[num]=(dp[num]+dp[num-p]*mult)%MOD;
            }
            if(num>3&&dp[num]==0)break;
            num++;
      }
      std::cout<<num;
}

プロジェクトオイラー Problem 77 「素数の和」 †

Problem 77 「素数の和」 †

10は素数の和として5通りに表すことができる:

7 + 3
5 + 5
5 + 3 + 2
3 + 3 + 2 + 2
2 + 2 + 2 + 2 + 2

素数の和としての表し方が5000通り以上になる最初の数を求めよ.

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

解法
素数だけを対象にした動的計画法で計算がすみます。
5000通りは少ないので
50000000000通りを超える場合を計算させていますが動的計画法なのでもちろん一瞬で答えが出ます。
コード中の
 if(sum>=50000000000)break;
を
 if(sum>=5000)break;
とすれば5000通りを超える場合が計算できます。




#include<stdio.h>
#include<vector>
#include<map>
#include<iostream>
bool is_prime(int num){
      if(num<2)return false;
      for(int i=2;i*i<=num;i++){
            if(num%i==0)return false;
      }
      return true;
}

int main(){
      std::map<int, std::map<int,__int64> > dp;
      std::map<int,__int64> map;

      std::map<int, __int64>::iterator it;
      std::vector<int> primes;
      dp[0][0]=0;
      int num=2;
      while(1){
            __int64 sum=0;
            if(is_prime(num)){
                  primes.push_back(num);
            }
            for(int i=0;i<primes.size();i++){
                  int p=primes[i];
                  if(dp.find(num-p)!=dp.end()){
                        it=dp[num-p].upper_bound(p);
                        it--;
                        sum+=(*it).second;
                  }
                  dp[num][p]=sum;
            }
            if(sum>=50000000000)break;
            if(is_prime(num)){
                  if(dp.find(num)==dp.end())dp[num]=map;
                  dp[num][num]=sum+1;
                  dp[num][0]=0;
            }
           
            num++;
      }
      printf("%d\n",num);
}

プロジェクトオイラー Problem 76 「和の数え上げ」 †

Problem 76 「和の数え上げ」 †

5は数の和として6通りに書くことができる:

4 + 1
3 + 2
3 + 1 + 1
2 + 2 + 1
2 + 1 + 1 + 1
1 + 1 + 1 + 1 + 1

2つ以上の正整数の和としての100の表し方は何通りか.

解法
素朴な動的計画法ではmemo[101][101]の配列が必要ですが。
Ｗｉｋｉに掲載されていた分割数の漸化式をそのまま実装します。
100は
Ｃ++


#include<stdio.h>
#include<iostream>

int pos(int n){
      n++;
      int n1=n/2*(n%2==0?1:-1);
      return (n1*(3*n1-1))/2;
}


int main(){
      __int64 memo[101]={0};
      memo[0]=1;
      for(int i=1;i<=100;i++){
            for(int j=1;j<=i;j++){
                  int p=pos(j);
                  if(i-p<0)break;
                  __int64 mult=((j-1)%4)/2==0?1:-1;
                  memo[i]+=memo[i-p]*mult;
            }
      }
      std::cout<<memo[100];
}

プロジェクトオイラー Problem 75 「1通りの整数直角三角形」 †

Problem 75 「1通りの整数直角三角形」 †

ある長さの鉄線を折り曲げたときに1通りの直角三角形を作る最短の長さは12 cmである. 他にも沢山の例が挙げられる.

12 cm: (3,4,5)
24 cm: (6,8,10)
30 cm: (5,12,13)
36 cm: (9,12,15)
40 cm: (8,15,17)
48 cm: (12,16,20)

それとは対照的に, ある長さの鉄線 (例えば20cm) は整数長さで折ることができない. また2つ以上の折り曲げ方があるものもある. 2つ以上ある例としては, 120cmの長さの鉄線を用いた場合で, 3通りの折り曲げ方がある.

120 cm: (30,40,50), (20,48,52), (24,45,51)

Lを鉄線の長さとする. 直角三角形を作るときに1通りの折り曲げ方しか存在しないような L ≤ 1,500,000 の総数を答えよ.

注: この問題は最近変更されました. あなたが正しいパラメータを使っているか確認してください.

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

解法
原始ピタゴラス数とその定数倍が150万以下のものをカウントするだけです。


#include<stdio.h>
const int LIMIT=1500*1000;
int counts[LIMIT+1]={0};
int gcd ( int a, int b ){
      int c;
      while ( a != 0 ) {
            c = a; a = b%a;  b = c;
      }
      return b;
}



int main(){
      int m=2;
      while(2*m*(m+1)<=LIMIT){
            int n=1;
            while(2*m*(m+n)<=LIMIT&&n<m){
                  if(((m-n)%2==0)||(gcd(m,n)!=1)){
                        //何もしない
                  }else{
                        int d=1;
                        while(2*m*(m+n)*d<=LIMIT){
                              counts[2*m*(m+n)*d]++;
                              d++;
                        }
                  }
                  n++;
            }
            m++;
      }
      int ans=0;
      for(int i=1;i<=LIMIT;i++){
            ans+=(counts[i]==1);
      }
      printf("%d\n",ans);
}

プロジェクトオイラー Problem 247 「双曲線下の正方形」 †

Problem 247 「双曲線下の正方形」 †

1 ≤ x, 0 ≤ y ≤ 1/x の領域について考える.

S1 をこの曲線の下に入る最大の正方形とする.
S2 を残りの空間に入る最大の正方形とし, これを繰り返す.
Sn のインデックスを (left, below) とする. left は Sn の左にある正方形の数を, below は Sn の下にある正方形の数を表す.

これらの正方形に番号を記したものを上の図に示す.
S2 は左に 1 個, 下に 0 個の正方形があるので, S2 のインデックスは (1,0) である.
S32 のインデックスは (1,1) であることがわかる. S50 のインデックスも同じである.
50 は (1,1) をインデックスに持つ Sn の中で, 最大の n である.

(3,3) をインデックスに持つ Sn の中で, 最大の n を求めよ.

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




解法
正方形の左下をx1,y1とすると
y=x-x1+y1とy=1/x交点のＸ座標とｙ１、交点のｙ座標とx1が次の正方形の位置となります。
交点のｘ座標と右下のｘ1の差が正方形のサイズです。
これを優先順位付きキューに入れて３、３がでてきた回数を調べます。
後は３，３のインデックスを持つ箱の積み上げかたは再帰関数で考えて
1  1  1  1
1  2  3.  4
1  3  6 10
1 4 10  20
で20通りしかありえないので20個めの3,3インデックス正方形が答えです。



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

struct S{
      double x1,y1,size;
      int down,left;
      bool operator<(const S& s)const{
            return size<s.size;
      }
};

double calcSize(S& s1){
      double x1=s1.x1;
      double y1=s1.y1;
      return (-y1+x1+sqrt((y1-x1)*(y1-x1)+4))/2.0-s1.x1;
}

int main(){
      S s1,s2;
      s1.x1=1;
      s1.y1=0;
      s1.size=calcSize(s1);
      s1.down=0;
      s1.left=0;
      int count=0,count33=0;
      std::priority_queue<S> pq;
      pq.push(s1);
      while(1){
            s2=s1=pq.top();
            pq.pop();
            count++;
            if(s1.down==3&&s1.left==3)count33++;
            if(count33==20)break;
            s2.x1+=s1.size;
            s2.size=calcSize(s2);
            s2.left++;
            pq.push(s2);

            s2=s1;
            s2.y1+=calcSize(s2);
            s2.down++;
            s2.size=calcSize(s2);
            pq.push(s2);
      }
      printf("%d\n",count);
}

プロジェクトオイラー 248 Numbers for which Euler’s totient function equals 13!

Numbers for which Euler’s totient function equals 13!

Problem 248

The first number n for which φ(n)=13! is 6227180929.
Find the 150,000^th such number.

http://projecteuler.net/problem=248


オイラーのファイ関数を通した結果が13！になるｎを小さいほうから見たとき15万番目となるｎを答えよ。



解法
ちょっとダメな感じです。
一分ルールを守れてるコードは書けたのですが。
無駄組み合わせもいくらか計算していてそれをファイ関数で力づくではじいています。
無駄組み合わせが一つも出ないのが理想なのですが？


コードの発想
オイラーのphi関数を通して13！と同じ素因数を持つものを生成します。
13！＝13*11*7*5＾２*3＾5*2＾10であることに注目します。

Ａ
ｎが13より大きな素因数ｐを持つ場合、その素因数ｐがｐ＾２などの二つ以上ｎに含まれる場合
phi関数を通しても13より大きな素因数を持つ数が生成されて不適。
よって13より大きな素因数ｐは一つでなくてはいけない。
またｐ＞13の素数の時ｐ-1が13より大きな素因数をもつならそれは13！の素因数でないものがphi関数を通して残ることになるのでｐ-1の素因数は13！の素因数、２，３、5,7,11,13でなくてはいけない。

Ｂ
ｎの13以下の素因数ｐはそのままｐが残っても問題ありません。
ｐ＾１でなくｐ＾２以上ならｐがphi関数を通してもその素因数が生き残ります。
ｎの素数ｐが13以下の素数なら、ｐ-1の２，３、5,7,11の素因数を考慮しておきます。

ＡとＢを合成して作ったphi関数の結果が
13！＝13*11*7*5＾２*3＾5*2＾10
となればそれは条件を満たします。


後はＡとＢの組み合わせを無駄なく全探索するだけです。

この発想で大丈夫なはずなのですが私のコードは何かミスがあるらしく。
無駄なパターンをいくつか余分に検討しています。
無駄なｎを生成した場合はph(n)iの結果が13!にならないならはじいています。
一応正答はしてるのですが何か微妙に間違ってるようです。


#include<stdio.h>
#include<set>
#include<iostream>

int counts[]={1,1,1,2,5,10};
int dells[6][6]={      {0,0,0,0,1,2},
                  {0,0,0,1,0,1},
                  {0,0,0,0,1,1},
                  {0,0,0,0,0,2},
                  {0,0,0,0,0,1},
                  {0,0,0,0,0,0}};
__int64 base[]={13,11,7,5,3,2};

std::set<__int64> all;

void saiki1(int,__int64);
void saiki2(int,int,__int64,__int64);


__int64 phi( __int64 n )
{
  __int64 res, p;

  res = n;
  if ( n % 2 == 0 ) {
    res /= 2;
    do { n /= 2; } while( n % 2 == 0 );
  }
  p = 3;
  while ( p <= n / p ) {
    if ( n % p == 0 ) {
      res = res / p * ( p - 1 );
      do { n /= p; } while( n % p == 0 );
    }
    p += 2;
  }
  if ( n > 1 ) { res = res / n * (n - 1); }

  return res;
}

bool isPrime(__int64 num){
      if(num<2)return false;
      for(__int64 i=2;i*i<=num;i+=1+(i%2)){
            if(num%i==0)return false;
      }
      return true;
}

void saiki2(int p1,int p2,__int64 m1,__int64 m2,__int64 loopCount){
      if(p2==6){
            if(m2>1&&isPrime(m2+1)==false)return ;
            if(loopCount!=0&&m2>13){
                  m2=m2;
            }else if(loopCount==0){
                  m2=0;
            }else{
                  return ;
            }
            bool ok=true;
            for(int i=0;i<6;i++){
                  if(counts[i]<dells[p1][i])ok=false;
            }
            if(ok==true){
                  for(int i=0;i<6;i++){
                        counts[i]-=dells[p1][i];
                  }
                  __int64 m3=1;
                  for(int i=0;i<=counts[p1];i++)m3*=base[p1];
                  saiki1(p1+1,m1*(m2+1)*m3);
                  for(int i=0;i<6;i++){
                        counts[i]+=dells[p1][i];
                  }
            }
            if(counts[p1]==0){
                  saiki1(p1+1,m1*(m2+1));
            }
            if(m2>13){
                  saiki2(p1,p1,m1*(m2+1),1,loopCount);
            }
      }else{
   
            __int64 m3=1;
            int temp=counts[p2];
            int start=0;
            if(p1==p2){
                  if(counts[p1]<loopCount)return ;
                  start=1;
                  if(start<loopCount)start=loopCount;
                  for(int i=0;i<start;i++){
                        m3*=base[p1];
                        counts[p1]--;
                  }
            }
         
            while(counts[p2]>=0){
                  if(p1==p2&&loopCount<start)loopCount=start;
                  int temp2=counts[p2];
                  saiki2(p1,p2+1,m1,m2*m3,loopCount);
                  counts[p2]=temp2-1;
                  m3*=base[p2];
                  start++;
            }
            counts[p2]=temp;
      }
}
void saiki1(int p1,__int64 m1){
      if(p1==6){
            if(phi(m1)==6227020800){
                  all.insert(m1);
            }else{
                  //std::cout<<"bad"<<m1<<" "<<phi(m1)<<"\n";
            }
      }else{
            saiki2(p1,p1,m1,1,1);
            saiki2(p1,6,m1,1,0);
      }

}

int main(){
      saiki1(0,1);
   
      std::set<__int64>::iterator it;
      std::cout<<all.size()<<"\n";
      int c=1;
      for(it=all.begin();it!=all.end();it++){
            if(c==150000){
                  std::cout<<"ans="<<(*it);
                  break;
            }
            c++;
      }
}

プロジェクトオイラー Problem 249 「素数の部分集合和」 †

Problem 249 「素数の部分集合和」 †

S = {2, 3, 5, ..., 4999} を 5000 より小さい素数の集合とする.

要素の合計が素数となるような, S の部分集合の数を求めよ.
最下位の 16 桁を答えとして入力せよ.

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





解法
これも問250と同じサービス問題のようで難易度は低いです。
ただ私の発想が悪いのか動的計画法の本体で速度がでてません。
オーソドックスな動的計画法を書いたつもりなのですが？
たぶんmod演算が遅いのでmod演算を実行回数へらしたらいいのかもしれません。
そう考えてmod演算を引き算に置き換えてみたところ驚くほど高速化しました。
Ｍｏｄ演算が遅いのはＣＰＵが悪いのかもしれません。



#include<stdio.h>
#include<string.h>
#include<iostream>
#include<vector>
__int64 dp[1600000];

bool is_prime(int n){
      if(n<2)return false;
      for(int i=2;i*i<=n;i+=1+(i&1)){
            if(n%i==0)return false;
      }
      return true;
}

int main(){
   
      __int64 MOD_LIMIT=100000000;
      MOD_LIMIT*=100000000;
      memset(dp,0,sizeof(dp));
      int now,next,sum=0;
      std::vector<int> primes;
      for(int i=2;i<5000;i++){
            if(is_prime(i)){
                  primes.push_back(i);
            }
      }
      for(int i=0;i<primes.size();i++){
            int p=primes[i];
            sum+=p;
            for(int j=sum;j>=p;j--){
                  dp[j]=(dp[j]+dp[j-p]);
                  if(dp[j]>=MOD_LIMIT)dp[j]-=MOD_LIMIT;
            }
            dp[p]+=1;
      }
      __int64 ans=0;
      for(int j=2;j<=sum;j++){
            if(is_prime(j)==true){
                  ans=(ans+dp[j])%MOD_LIMIT;
            }
      }
      std::cout<<ans<<"\n";
}

プロジェクトオイラー Problem 250 「250250」 †

Problem 250 「250250」 †

要素の合計が 250 で割り切れるような, {11, 22, 33,..., 250250250250} の空でない部分集合の数を求めよ. 最下位の 16 桁を答えとして入力せよ.

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






解法
難易度の低いサービス問題のようです。
何も考えない動的計画法で片が付きます。
難易度低すぎで他の難問の中で浮いてますねこの問題。
プログラム　Ｃ++　ＢＣＣ5.5


#include<stdio.h>
#include<string.h>
#include<iostream>

const int MOD=250;

int mod_pow(int n){
      int r=n,baseM=n % MOD,result=1;
      while(r>0){
            if(r%2==1){
                  result=(result*baseM) % MOD;
            }
            r/=2;
            baseM=(baseM*baseM) % MOD;
      }
      return result;
}
int main(){
      __int64 dp[2][MOD];
      __int64 MOD_LIMIT=100000000;
      MOD_LIMIT*=100000000;
      memset(dp,0,sizeof(dp));
      int now,next;
      for(int i=1;i<=250250;i++){
            now=i%2;
            next=(i+1)%2;
            memcpy(dp[next],dp[now],sizeof(dp[next]));
            int add=mod_pow(i);
            dp[next][add]+=1;
            for(int i=0;i<MOD;i++){
                  dp[next][(i+add) % MOD]=(dp[next][(i+add)%MOD]+dp[now][i])%MOD_LIMIT;
            }
      }
      std::cout<<dp[next][0]<<"\n";
}

プロジェクトオイラー Problem 74 「桁の階乗による連鎖」 †

Problem 74 「桁の階乗による連鎖」 †

145は各桁の階乗の和が145と自分自身に一致することで有名である.

1! + 4! + 5! = 1 + 24 + 120 = 145

169の性質はあまり知られていない. これは169に戻る数の中で最長の列を成す. このように他の数を経て自分自身に戻るループは3つしか存在しない.

169 → 363601 → 1454 → 169
871 → 45361 → 871
872 → 45362 → 872

どのような数からスタートしてもループに入ることが示せる.

例を見てみよう.

69 → 363600 → 1454 → 169 → 363601 (→ 1454)
78 → 45360 → 871 → 45361 (→ 871)
540 → 145 (→ 145)

69から始めた場合, 列は5つの循環しない項を持つ. また100万未満の数から始めた場合最長の循環しない項は60個であることが知られている.

100万未満の数から開始する列の中で, 60個の循環しない項を持つものはいくつあるか?

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

解法
数字の中身を辞書順に並べ変えたものを一纏めにして検討します。
3345,3354,3435,3453,3534,3543,4335,4353,4533,5334,5343,5433,
などは全部次の数が同じになります。

ということでこれを辞書順の3345で代表させて長さを調べます。
辞書順で昇順になってるものだけ生成してから検討します。

最後に長さ60であることが判明したら、それの並べ替え数を計算して集計してアセプトです。
昇順なら20000個程度の数字を検討するだけで答えが出ます。





コード実行結果
100万までの場合
10 ?- time(main74).
402
% 17,362,257 inferences, 2.500 CPU in 2.494 seconds (100% CPU, 6944903 Lips)


1000万までの場合
11 ?- time(main74).
11322
% 43,154,394 inferences, 6.281 CPU in 6.281 seconds (100% CPU, 6870351 Lips)

メモ化のための配列や木を使ってないことを考えたらまあまあの速度だと思います。









以下この問いのＰｒｏｌｏｇ言語による記述


sum([],0):-!.
sum([X|Xs],Result):-
      !,
      sum(Xs,Re),
      Result is Re+X.

to_num([],[]):-!.
to_num([X|Xs],[Y|Result]):-
      !,
      Y is X-48,
      to_num(Xs,Result).

fact(0,1):-!.
fact(N,Result):-
      !,
      N1 is N-1,
      fact(N1,Re),
      Result is Re*N.

factC(48,1):-!.
factC(N,Result):-
      !,
      N1 is N-1,
      factC(N1,Re),
      Result is Re*(N-48).

calc_next(List,Result):-
      findall(F,(member(E,List),factC(E,F)),List1),
      sum(List1,Num),
      number_codes(Num,Result).



count([],[N,Count],[[N,Count]]):-!.
count([N|Rest],[N,Count],Result):-
      !,
      Count1 is Count+1,
      count(Rest,[N,Count1],Result).
count([N|Rest],[N1,Count],[[N1,Count]|Result]):-
      !,
      count(Rest,[N,1],Result).

count_w([Top|List],Result):-
      count(List,[Top,1],Result).

calc_perm3([],Perm,Perm):-!.
calc_perm3([[_,X]|Xs],Perm,Result):-
      !,
      fact(X,F),
      Perm1 is Perm//F,
      calc_perm3(Xs,Perm1,Result).

calc_perm2(Counts,AllPerm,Result):-
      between(1,9,N),
      select([N,Count],Counts,Counts1),
      Count1 is Count-1,
      calc_perm3([[N,Count1]|Counts1],AllPerm,Result)
      .


calc_perm(List,Result):-
      !,
      length(List,Len),
      Len1 is Len-1,
      fact(Len1,AllPerm),
      to_num(List,List1),
      count_w(List1,Counts),
      findall(P,calc_perm2(Counts,AllPerm,P),Perms),
      sum(Perms,Result).

seed(_,48,6):-!,fail.
seed([],_,6):-!.
seed([],Down,Keta):-Keta>0,Down>48.
seed([Down1|Result],Down,Keta):-
      !,
      between(Down,57,Down1),
      Keta1 is Keta+1,
      seed(Result,Down1,Keta1).

calc(List,Len,Old):-
      number_codes(Num,List),
      member(Num,Old),
      !,
      Len=:=60.
calc(List,Len,[_,Old2,Old3]):-
      !,
      number_codes(Num,List),
      Len1 is Len+1,
      calc_next(List,List1),
      calc(List1,Len1,[Old2,Old3,Num]).


search(Perm):-
      seed(List,48,0),
      calc(List,0,[-1,-1,-1]),
      calc_perm(List,Perm).

main74:-
      findall(E,search(E),List),
      sum(List,Sum),
      write(Sum).

プロジェクトオイラー Problem 73 「ある範囲内の分数の数え上げ」 †

nとdを正の整数として, 分数 n/d を考えよう. n<d かつ HCF(n,d)=1 のとき, 真既約分数と呼ぶ.

d ≤ 8 について既約分数を大きさ順に並べると, 以下を得る:

1/8, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3, 5/7, 3/4, 4/5, 5/6, 6/7, 7/8

1/3と1/2の間には3つの分数が存在することが分かる.

では, d ≤ 12,000 について真既約分数をソートした集合では, 1/3 と 1/2 の間に何個の分数があるか?

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

解法
ファレイ数列を元に高速化できないかいろいろ数式変形を考えてみましたが。
纏めて計算しようとしてもΣ計算が複雑につみあがるばかりでどう考えたらいいのわからず。
結局全探索。


search(L,R):-
      L+R>12000,
      !,
      fail.
search(_,_).

search(L,R):-
      L1 is L+R,
      search(L1,R).
search(L,R):-
      !,
      R1 is L+R,
      search(L,R1).
start(1):-
      search(3,2).
main73:-
      findall(E,start(E),List),
      length(List,Ans),
      write(Ans).

プロジェクトオイラー Problem 72 「分数の数え上げ」 †

Problem 72 「分数の数え上げ」 †

nとdを正の整数として, 分数 n/d を考えよう. n<d かつ HCF(n,d)=1 のとき, 真既約分数と呼ぶ.

d ≤ 8について真既約分数を大きさ順に並べると, 以下を得る:

1/8, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3, 5/7, 3/4, 4/5, 5/6, 6/7, 7/8

この集合は21個の要素をもつことが分かる.

d ≤ 1,000,000について, 真既約分数の集合は何個の要素を持つか?

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

解法
２～100万=n Σオイラーのファイ関数（ｎ）が答えです。
エラトステネスの篩もどきで計算が終わります。


#include<stdio.h>
#include<iostream>
const int LIMIT=1000*1000;
int memo[LIMIT+1]={0};

int main(){
      __int64 ans=0;
     
      for(int i=2;i<=LIMIT;i++){
            if(memo[i]!=0){
                  ans+=memo[i];
                  continue;
            }else{
                  ans+=i-1;
            }
            for(int j=i+i;j<=LIMIT;j+=i){
                  if(memo[j]==0)memo[j]=j;
                  memo[j]=(memo[j]/i)*(i-1);
            }
      }
      std::cout<<ans<<"\n";
}

プロジェクトオイラー Problem 71 「順序分数」 †

Problem 71 「順序分数」 †

nとdを正の整数として, 分数 n/d を考えよう. n<d かつ HCF(n,d)=1 のとき, 真既約分数と呼ぶ.

d ≤ 8について既約分数を大きさ順に並べると, 以下を得る:

1/8, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3, 5/7, 3/4, 4/5, 5/6, 6/7, 7/8

3/7のすぐ左の分数は2/5である.

d ≤ 1,000,000について真既約分数を大きさ順に並べたとき, 3/7のすぐ左の分数の分子を求めよ.

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

解法

Ｐｒｏｌｏｇ言語

ファレイ数列より

X is (1000*1000-5)//7,Ans is X*3+2.

が答え。

プロジェクトオイラー Problem 70 「トーティエント関数の置換」 †

Problem 70 「トーティエント関数の置換」 †

オイラーのトーティエント関数 φ(n) (ファイ関数とも呼ばれる) とは, n 未満の正の整数で n と互いに素なものの個数を表す. 例えば, 1, 2, 4, 5, 7, 8 は9未満で9と互いに素であるので, φ(9) = 6 となる.
1 は全ての正の整数と互いに素であるとみなされる. よって φ(1) = 1 である.

面白いことに, φ(87109)=79180 であり, 87109は79180を置換したものとなっている.

1 < n < 107 で φ(n) が n を置換したものになっているもののうち, n/φ(n) が最小となる n を求めよ.

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



解法
素数は答えになりません。
nとｎ-1各桁の数字の個数が絶対合わないからです。

次の候補として
2素因数で1000万に最も近い数が答えの候補となります。
3素因数の場合もあり得ますが確率は大幅に下がります。
そして2素因数で正解が見つかります。


not_prime(N):-N<2,!.
not_prime(N):-
      between(2,N,D),
      (N<D*D->!,fail;true),
      N mod D=:=0,
      !.
is_prime(N):-not(not_prime(N)).

prime_list(N):-
      between(2000,5000,N),
      is_prime(N).

calc(Primes,[D,P3]):-
      member(P1,Primes),
      member(P2,Primes),
      P1<P2,
      P3 is P1*P2,
      P3<1000*10000,
      P4 is (P1-1)*(P2-1),
      number_codes(P3,List3),
      number_codes(P4,List4),
      msort(List3,List33),
      msort(List4,List33),
      D is P3/P4.

main70:-
        findall(P,prime_list(P),Primes),
      findall(E,calc(Primes,E),AnsList),
      msort(AnsList,[Ans|_]),
      write(Ans)
      .

プロジェクトオイラー Problem 69 「トーティエント関数の最大値」 †

Problem 69 「トーティエント関数の最大値」 †

オイラーのトーティエント関数, φ(n) [時々ファイ関数とも呼ばれる]は, n と互いに素な n 未満の数の数を定める. たとえば, 1, 2, 4, 5, 7, そして8はみな9未満で9と互いに素であり, φ(9)=6.

n	互いに素な数	φ(n)	n/φ(n)
2	1	1	2
3	1,2	2	1.5
4	1,3	2	2
5	1,2,3,4	4	1.25
6	1,5	2	3
7	1,2,3,4,5,6	6	1.1666...
8	1,3,5,7	4	2
9	1,2,4,5,7,8	6	1.5
10	1,3,7,9	4	2.5

n ≤ 10 では n/φ(n) の最大値は n=6 であることがわかる.

n ≤ 1,000,000で n/φ(n) が最大となる値を見つけよ.

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

解法

答え

ファイ関数の定義より答えは

2*3*5*7*11*13*17

プロジェクトオイラー Problem 68 「Magic 5-gon ring」 †

Problem 68 「Magic 5-gon ring」 †

下に示す図のようなものを"magic" 3-gon ringという. これは1～6の数字を当てはめて, 各列の数字の和が9となっている. これを例として説明する.

外側のノードのうち一番小さいものの付いた列(例では4,3,2)から時計回りに回ってそれぞれ列の数字を3つ連ねて説明する. 例えば例のものは4,3,2; 6,2,1; 5,1,3という組で説明することができる.

1～6の数字を当てはめて, 各列の数字の和が等しくなるものは次の8通りある.

合計	組
9	4,2,3; 5,3,1; 6,1,2
9	4,3,2; 6,2,1; 5,1,3
10	2,3,5; 4,5,1; 6,1,3
10	2,5,3; 6,3,1; 4,1,5
11	1,4,6; 3,6,2; 5,2,4
11	1,6,4; 5,4,2; 3,2,6
12	1,5,6; 2,6,4; 3,4,5
12	1,6,5; 3,5,4; 2,4,6

この組の各数字を連結して, 9桁の数字で表すことができる. 例えば, 上の図のものは4,3,2; 6,2,1; 5,1,3であるので432621513である.

さて, 下の図に1～10の数字を当てはめ, 各列の数字の和が等しくなる"magic" 5-gon ringを作って, それを表す16桁または17桁の数字のうち, 16桁のものの最大の数字を答えよ.

(注, 3つの場合の例を見ても分かる通り, 列の始まりの数字を比べた時一番小さい数字で始まる列から時計回りに繋げるという条件のもとで文字列を生成する必要があります. この条件下で最大となる数字を答えてください. )

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


解法
条件を満たすか調べながら全探索するだけでも解けます。
Ｐｒｏｌｏｇ言語

select_num(N,Nums,Nums):-integer(N),!.
select_num(N,Nums,Nums1):-!,select(N,Nums,Nums1).

ok(-1,Xa,Xa):-!.
ok(X,Xa,X):-X<Xa,!.

sum_ok(-1,Xa,Xb,Xc,Sum):-!,Sum is Xa+Xb+Xc.
sum_ok(Sum,Xa,Xb,Xc,Sum):-!,Sum=:=Xa+Xb+Xc.

calc([],_,_,[]):-
      !.
calc([[Xa,Xb,Xc]|Board],Sum,FirstNum,Nums):-
      select_num(Xa,Nums,Nums1),
      ok(FirstNum,Xa,FirstNum1),
      select_num(Xb,Nums1,Nums2),
      select_num(Xc,Nums2,Nums3),
      Xb<10,
      Xc<10,
      sum_ok(Sum,Xa,Xb,Xc,Sum1),
      calc(Board,Sum1,FirstNum1,Nums3).

main68:-
      Board=[[_,X2,X3],[_,X3,X5],[_,X5,X7],[_,X7,X9],[_,X9,X2]],
      calc(Board,-1,-1,[1,2,3,4,5,6,7,8,9,10]),
      write(Board).