problem 7
今日はもうやらないつもりだったのだけど，ちらっと見たらまた素数だったので problem 3 のコードを流用した．
#include <stdio.h> #include <stdlib.h> #define SIZE 10001 int main() { unsigned long long int *primes = (long long int *)malloc(sizeof(long long int)*SIZE); int n=5,i=2,j; primes[0]=2; primes[1]=3; for(;i<=10001;n+=2) { long long int f=1; for(j=1;j<i&&f;j++) { f = 1; f*=n%primes[j]; } if(f) primes[i++]=n; } printf("%lld\n",primes[10000]); free(primes); return 0; }
追記: problem 6 かと思ったら 7 だったので訂正．
problem 5
Problem 3 のコードの使い回しを行った．
#include <stdio.h> #define SIZE 20 int main() { unsigned int multiples[SIZE][2]; unsigned int primes[SIZE]; int result=1,max=20,n=4,i=2,j=0,k; primes[0]=2; primes[1]=3; for(;n <= max;n+=1) { int f=1; for(k=0;k<i&&f;k++) { f=1; f*=n%primes[k]; } if(f) primes[i++]=n; else { f=1; for(k=i-1;k>=0&&f;k--) { f=1; f*=n%primes[k]; } multiples[j][0]=n; multiples[j++][1]=n/primes[k+1]; } } for(k=0;k<i;k++) { result *= primes[k]; } for(k=0;k<j;k++) { result *= (result%multiples[k][0])?multiples[k][1]:1; } printf("%d\n", result); return 0; }
problem 1-4
Problem 1:
#include <stdio.h> int main() { int n=0,i=0,j; while(i<997)n+=(i+=3); for(i=5;i<1000;i+=15)n+=((j=i+5)<1000)?i+j:i; printf("%d\n", n); }
Problem2:
#include <stdio.h> int main() { int n=0,t1=1,t2=2; while(t2 < 4000000) { n+=(t2%2==0)?t2:0; t2+=t1; t1=t2-t1; } printf("%d\n", n); }
Problem3:
#include <stdio.h> #include <stdlib.h> #define SIZE 10000 int main() { unsigned long long int *result=(long long int *)malloc(sizeof(long long int)*SIZE); unsigned long long int *primes=(long long int *)malloc(sizeof(long long int)*SIZE); long long int max=600851475143LL,n=5,i=2,j,k; primes[0]=2; primes[1]=3; for(;n*n < max;n+=2) { long long int f=1; for(k=1;k<i&&f;k++,f=1) { f*=n%primes[k]; } if(f) primes[i++]=n; } for(j=0,k=0;j<i;j++) { if((max%primes[j])==0) result[k++]=primes[j]; } for(j=0;j<k;j++) { printf("%lld ", result[j]); } printf("\n"); free(result); free(primes); return 0; }
Problem4:
#include <stdio.h> #include <math.h> int digits(int n) { int i=1,j=0; if(n==0) return 1; while(n>=i) i*=10,j+=1; return j; } int nth(int n, int i) { if(i==1) return n-(n/10)*10; return nth(n/10,i-1); } int isPalindrome(int n, int d_prev) { int i,j,d; if((d=digits(n))<d_prev-2) { if(n==0) return 1; for(i=0;i<d_prev-2-d;i++) { if(nth(n,1)!=0) return 0; n /= 10; } return isPalindrome(n,digits(n)+2); } if(d==1) return 1; if(d==2) { for(i=1;i<10;i++) if(i*11==n) return 1; return 0; } return ((i=nth(n,d))==(j=nth(n,1)))?isPalindrome((n-i*(int)pow(10,d-1)-j)/10,d):0; } int main() { int i,j,max=0; for(i=999;i>100;i--) { for(j=999;j>100;j--){ int n = i*j; if(isPalindrome(n,digits(n)+2)) max=(n>max)?n:max; } } printf("%d\n",max); }
優先順位つきキューを用いたエラトステネスの篩の変形
Haskellでエラトステネスの篩 - 簡潔なQ
上の記事に出てくる「優先順位つきキューを用いたエラトステネスの篩の変形」が C++ で書かれたコードを読んでもよく分からなかったので本人に尋ねてみた所，そもそも優先順位つきキュー自体がよく分かっていないことが判明し，それも含めて解説してもらった．
Sieve of Eratosthenes using a priority queue · GitHub
それでとりあえず分かった気になっていたのだけど，なんかひっかかるなあという感じになったので少し考えてみたら割と単純な話だった:
- i の初期値は素数 2 である．
- キューには i を越える最小の偶数 (素数 2 の倍数) が必ず含まれている．
- キューから pop した値を x としたとき 1 ≦ x - i ≦ 2 である．つまり i と x は隣り同士，または奇数 q を挟む形になっている．
- i と x が隣り同士で且つ i が偶数の場合，x は奇数である．このとき x は 2 以外で i 未満の素数 p の倍数である．従って i も x も素数ではないので，i に x + 1 を代入し，キューに x + 2p (x + p は偶数であるため既にキュー内に存在する) を追加する．
- i と x が奇数 q を挟む場合，q は素数．なぜなら，i 未満の素数において i を越える最小の倍数を y とすると x ≦ y が成り立ち q ≠ y となるからである．i に x を代入，キューに q^2 を追加して処理を続ける．
- (4.), (5.) より i は奇数になり得ないことに注意．
ちなみに 4. の
i 未満の素数において i を越える最小の倍数を y とすると x ≦ y が成り立ち
は
// 該当部分を抜き出したコード int p = i+1; pq.push(make_pair(-p*p, p)); pq.push(make_pair(-x-q, q)); i = x;
から分かる．
SKI combinator calculus
*Main> termToString $ ski (E2 (E3 S (E2 K (Var 'x')) (E3 S I I)) (Var 'y')) "xyy"
という感じで動作するプログラム．某 Skype チャットに SKI という単語が出てきてふと思い出したので書いてみた．例によって間違ってる可能性ありありなので適当に指摘してもらえるとありがたいのです．
data Term = Var Char | I | K | S | E2 Term Term | E3 Term Term Term | E4 Term Term Term Term ski :: Term -> Term ski (E2 t1 t2) = eval (E2 (ski t1) (ski t2)) ski (E3 t1 t2 t3) = eval (E3 (ski t1) (ski t2) (ski t3)) ski (E4 t1 t2 t3 t4) = eval (E4 (ski t1) (ski t2) (ski t3) (ski t4)) ski t = t eval :: Term -> Term eval (E2 I t) = eval t eval (E3 K t1 t2) = eval t1 eval (E4 S t1 t2 t3) = eval (E2 (eval (E2 (eval t1) (eval t3))) (eval (E2 (eval t2) (eval t3)))) eval (E2 (E2 t1 t2) t3) = eval (E3 t1 t2 t3) eval (E2 (E3 t1 t2 t3) t4) = eval (E4 t1 t2 t3 t4) eval (E2 (E2 t1 t2) (E2 t3 t4)) = eval (E3 t1 t2 (E2 t3 t4)) eval (E3 (E2 t1 t2) t3 t4) = eval (E4 t1 t2 t3 t4) eval t = t termToString :: Term -> String termToString S = "S" termToString K = "K" termToString I = "I" termToString (Var c) = [c] termToString (E2 t1 t2) = termToString t1 ++ termToString t2 termToString (E3 t1 t2 t3) = termToString t1 ++ termToString t2 ++ termToString t3 termToString (E4 t1 t2 t3 t4) = termToString t1 ++ termToString t2 ++ termToString t3 ++ termToString t4
エラトステネスの篩
よく haskell コードの例として出される
primes = sieve [2..] sieve (p:xs) = p : sieve [x | x <- xs, x `mod` p /= 0]
は遅いよ，という割と有名な話を昨日読んでみたんだけど，なんだか腑に落ちない部分があったのでここに書いてみる．参考にした記事は以下の 3 つ:
- http://lowlife.jp/mft/weblog/2006/04/21.html
- http://haskell.g.hatena.ne.jp/nobsun/20060618
- http://d.hatena.ne.jp/kazu-yamamoto/20100624/1277348961
まず上の 2 つを読んで「篩に掛けられた後のリストの先頭を n とすると n^2 未満の (それまでの篩に引っかからなかった) 数は素数だからそれを利用して計算量を減らそう」ていうことだろうと適当に理解したつもりだったんですが，3 番目を読むと
エラトステネスの方法で、たとえば 5 の倍数を消すとき、5 の次は 10 を対象とする。7 は相手にしない。(6、8、9 はすでに消えてる。)
しかし、このコードでは 7 を 5 で割ってしまう。これでは誇大広告と言われてもしかたたないだろう。
という風に書いてあって，なんか自分の理解と違うなあという感じで頭が混乱してるなうというわけです．それで，もやもやな疑問を放置するのもあれなのでとりあえずこの記事の意図するのはどういうことだろうと考えてみたのが，
- 等差数列ということで (この例だと) 5 の位置から 4 個飛ばしで値を除いていけば速いんじゃね? という話
でもこれも間違いで，一見すると正しいようにみえるけどその前の篩でいくつか数が抜け落ちているため 4 個飛ばしで値を除くという処理自体が不可能．
そういうわけで Haskell さんむずかしくて意味が分からない感じなので誰か助言ください…．
isTaut 2
抜けている条件があったので補ったら先程の論理式も正しく判定するようになった．やりましたね!
*Main> isTaut (Not (And (And (Or (Var 'x') (Var 'y')) (Or (Not (Var 'x')) (Var 'y'))) (And (Or (Var 'x') (Not (Var 'y'))) (Or (Not (Var 'x')) (Not (Var 'y')))))) Const True
data Prop = Const Bool | Var Char | Not Prop | Or Prop Prop | And Prop Prop | Imply Prop Prop deriving (Eq, Ord, Show) isTaut :: Prop -> Prop isTaut (Const b) = (Const b) isTaut (Var x) = (Var x) isTaut (Not p) = rule (Not (isTaut p)) isTaut (Or p1 p2) = rule (Or (isTaut p1) (isTaut p2)) isTaut (And p1 p2) = rule (And (isTaut p1) (isTaut p2)) isTaut (Imply p1 p2) = rule (Or (isTaut (Not p1)) (isTaut p2)) rule :: Prop -> Prop rule (Const b) = (Const b) rule (Not (Const b)) = (Const (not b)) rule (Or (Const b1) (Const b2)) = (Const (b1 || b2)) rule (And (Const b1) (Const b2)) = (Const (b1 && b2)) rule (Imply (Const b1) (Const b2)) = (Const (not b1 && b2)) rule (Var x) = (Var x) rule (Not (Not p)) = rule p rule (Not (Or p1 p2)) = rule (And (rule (Not (rule p1))) (rule (Not (rule p2)))) rule (Not (And p1 p2)) = rule (Or (rule (Not (rule p1))) (rule (Not (rule p2)))) rule (Or (Const True) p) = (Const True) rule (Or p (Const True)) = (Const True) rule (Or (Const False) p) = p rule (Or p (Const False)) = p rule (Or (Not p1) p2) | (rule p1) == (rule p2) = (Const True) rule (Or p1 (Not p2)) | p1 == p2 = (Const True) rule (Or p1 p2) | p1 == p2 = p1 rule (Or (Or p1 p2) p3) | p1 == p3 = rule (Or p2 p3) | p2 == p3 = rule (Or p1 p3) rule (Or (Or p1 p2) p3) | (p1 == (Not p3) || p2 == (Not p3)) = (Const True) rule (Or (Or p1 p2) p3) | ((Not p1) == p3 || (Not p2) == p3) = (Const True) rule (Or p1 (Or p2 p3)) | p1 == p2 = rule (Or p1 p3) | p1 == p3 = rule (Or p1 p2) rule (Or p1 (Or p2 p3)) | ((Not p1) == p2 || (Not p1) == p3) = (Const True) rule (Or p1 (Or p2 p3)) | (p1 == (Not p2) || p1 == (Not p3)) = (Const True) rule (Or (And p1 p2) p3) | (p1 == p3 || p2 == p3) = p3 rule (Or (And p1 p2) p3) | p1 == (Not p3) = rule (And p2 p3) | p2 == (Not p3) = rule (And p1 p3) rule (Or p1 (And p2 p3)) | (p1 == p2 || p1 == p3) = p1 rule (Or p1 (And p2 p3)) | (Not p1) == p3 = rule (And p1 p2) | (Not p2) == p3 = rule (And p1 p3) rule (Or (Or p1 p2) (Or p3 p4)) | (p1 == p3 || p2 == p3) = rule (Or (Or p1 p2) p4) | (p1 == p4 || p2 == p4) = rule (Or (Or p1 p2) p3) | ((Not p1) == p3 || (Not p2) == p3 || p1 == (Not p3) || p2 == (Not p3)) = (Const True) | ((Not p1) == p4 || (Not p2) == p4 || p1 == (Not p4) || p2 == (Not p4)) = (Const True) rule (Or (And p1 p2) (And p3 p4)) | p1 == p3 = rule (And p1 (Or p2 p4)) | p1 == p4 = rule (And p1 (Or p2 p3)) | p2 == p3 = rule (And p2 (Or p1 p4)) | p2 == p4 = rule (And p2 (Or p1 p3)) rule (Or (And p1 p2) (Or p3 p4)) | p1 == p3 = rule (Or p1 p4) | p1 == p4 = rule (Or p1 p3) | p2 == p3 = rule (Or p2 p4) | p2 == p4 = rule (Or p2 p3) | ((Not p1) == p3 || (Not p1) == p4) = rule (Or (Or p3 p2) p4) | ((Not p2) == p3 || (Not p2) == p4) = rule (Or (Or p3 p1) p4) | (p1 == (Not p3) || p1 == (Not p4)) = rule (Or (Or p3 p2) p4) | (p2 == (Not p3) || p2 == (Not p4)) = rule (Or (Or p3 p1) p4) rule (Or (Or p3 p4) (And p1 p2)) | p1 == p3 = rule (Or p1 p4) | p1 == p4 = rule (Or p1 p3) | p2 == p3 = rule (Or p2 p4) | p2 == p4 = rule (Or p2 p3) | ((Not p1) == p3 || (Not p1) == p4) = rule (Or (Or p3 p2) p4) | ((Not p2) == p3 || (Not p2) == p4) = rule (Or (Or p3 p1) p4) | (p1 == (Not p3) || p1 == (Not p4)) = rule (Or (Or p3 p2) p4) | (p2 == (Not p3) || p2 == (Not p4)) = rule (Or (Or p3 p1) p4) rule (And (Const False) p) = (Const False) rule (And p (Const False)) = (Const False) rule (And (Const True) p) = p rule (And p (Const True)) = p rule (And (Not p1) p2) | p1 == p2 = (Const False) rule (And p1 (Not p2)) | p1 == p2 = (Const False) rule (And p1 p2) | p1 == p2 = p1 rule (And (And p1 p2) p3) | p1 == p3 = rule (And p2 p3) | p2 == p3 = rule (And p1 p3) rule (And (And p1 p2) p3) | (p1 == (Not p3) || p2 == (Not p3)) = (Const False) rule (And (And p1 p2) p3) | ((Not p1) == p3 || (Not p2) == p3) = (Const False) rule (And p1 (And p2 p3)) | p1 == p2 = rule (And p1 p3) | p1 == p3 = rule (And p1 p2) rule (And p1 (And p2 p3)) | ((Not p1) == p2 || (Not p1) == p3) = (Const False) rule (And p1 (And p2 p3)) | (p1 == (Not p2) || p1 == (Not p3)) = (Const False) rule (And (Or p1 p2) p3) | (p1 == p3 || p2 == p3) = p3 rule (And (Or p1 p2) p3) | p1 == (Not p3) = rule (Or p2 p3) | p2 == (Not p3) = rule (Or p1 p3) rule (And p1 (Or p2 p3)) | (p1 == p3 || p1 == p2) = p1 rule (And p1 (Or p2 p3)) | (Not p1) == p3 = rule (Or p1 p2) | (Not p2) == p3 = rule (Or p1 p3) rule (And (And p1 p2) (And p3 p4)) | (p1 == p3 || p2 == p3) = rule (And (And p1 p2) p4) | (p1 == p4 || p2 == p4) = rule (And (And p1 p2) p3) | ((Not p1) == p3 || (Not p2) == p3 || p1 == (Not p3) || p2 == (Not p3)) = (Const False) | ((Not p1) == p4 || (Not p2) == p4 || p1 == (Not p4) || p2 == (Not p4)) = (Const False) rule (And (Or p1 p2) (Or p3 p4)) | p1 == p3 = rule (Or p1 (And p2 p4)) | p1 == p4 = rule (Or p1 (And p2 p4)) | p2 == p3 = rule (Or p2 (And p1 p4)) | p2 == p4 = rule (Or p2 (And p1 p3)) rule (And (Or p1 p2) (And p3 p4)) | p1 == p3 = rule (And p1 p4) | p1 == p4 = rule (And p1 p3) | p2 == p3 = rule (And p2 p4) | p2 == p4 = rule (And p2 p3) | ((Not p1) == p3 || (Not p1) == p4) = rule (And (And p2 p4) p3) | ((Not p2) == p3 || (Not p2) == p4) = rule (And (And p1 p4) p3) | (p1 == (Not p3) || p1 == (Not p4)) = rule (And (And p3 p2) p4) | (p2 == (Not p3) || p2 == (Not p4)) = rule (And (And p3 p1) p4) rule (And (And p3 p4) (Or p1 p2)) | p1 == p3 = rule (And p1 p4) | p1 == p4 = rule (And p1 p3) | p2 == p3 = rule (And p2 p4) | p2 == p4 = rule (And p2 p3) | ((Not p1) == p4 || (Not p1) == p3) = rule (And (And p2 p4) p3) | ((Not p2) == p4 || (Not p2) == p3) = rule (And (And p1 p4) p3) | (p1 == (Not p3) || p1 == (Not p4)) = rule (And (And p3 p2) p4) | (p2 == (Not p3) || p2 == (Not p4)) = rule (And (And p3 p1) p4) rule p = p
