2013-11-27

pari

pariの使い方すぐに忘れるのでメモ

/*
 コメント
*/
\\ コメント
\\ 特に指定しない場合は変数は全てグローバル扱いになる。
\\　関数定義内で局所変数を使う場合は localまたはmyを使う。
\\ 初期化されていない変数は、自動的に自由変数扱いになる。
\\ 自由変数を明示的に使う場合は 'x などとする。
\\ 例:
f(a)={
	local(x='x);
	print(x^2+x+a);
};
x=3;
f(5); \\ x^2+x+5と表示される。

\\ for文等で複数行に渡る文を書く場合には{}で囲む
{for(i1=0, 3,
	for(i2=0, 3, 
		print(i1, " ", i2, " ", i1*i2);
	);
);}

if(x==0, error("abc")); \\エラーを投げる

\\ モジュラー演算
lift(Mod(1/7,100)) = 43

\\ 有限体
T = ffinit(2,8);
y = ffgen(T, 'y);
h = y.pol; \\ hはyの持ち上げになる。h=lift(y)とかでは持ちあげられない。

\\ p進数
log(1+3+O(3^4))

\\ 多項式
Polrev([1,2,3]) => 1+2*x+3*x^2
Pol([1,2,3]) => x^2+2*x+3

\\ 二次体
w = quadgen(5);
w^10 => 34 + 55*w

\\ ベクトル
[1,2]*[3,4]~ => 11
[1,2]~*[3,4] => [3, 4] [6,8] 
v = [3,1,4,1,5,9,2,6,5]
v[1] => 1
v[3..5] => [4,1,5]
vector(5, i, 2^i) => [2,4,8,16,32]
Vecrev([(1+x+x^3)]) => [1,1,0,1]
primes(5) => [2, 3, 5, 7, 11]
[x+2 | x<-primes(x), x>2] => [5,7,9,13]

\\行列
Mat(2, 2, i, j, 2*i+j);
[3, 4; 5, 6]

\\1の冪根
sqrtn(Mod(1,13),4, &z) => Mod(1,13) \\ z=Mod(5,13)となる。

2011-12-08

ペル方程式の解の存在をCoqで証明してみた

初めに

この記事はTheorem Proving Advent Calendar 2011の８日目の記事です。
以前、ペル方程式の解の存在証明をCoqで書いたので、その紹介をしようと思います。
anarchy proofでPell's equationとして投稿されていた問題です。
ちなみにFormalizing 100 theorems in Coqというリストがあるのですが、
そこには、まだCoqで証明されてない命題として載っているようなので(39. Solutions to Pell’s Equation)、Coqで証明されたのは初めてかもしれないです。
この記事の内容は、証明方針と証明してみた感想等を述べるのが目的で、Coq自体についての解説等はないのでご了承下さい。

証明する主張

「nを0より大きい整数で平方数でないとすると、x*x-n*y*y=1となる整数x,yでy≠0となるものが存在する。」
具体例としては次のようなものとなります。
n=2のとき(x,y)=(3,2) など.
n=3のとき(x,y)=(2,1) など.
n=5のとき(x,y)=(9,4) など.

coqで主張を書くと

forall (n: Z), n > 0 -> { x:Z & { y:Z | x*x - n*y*y = 1 /\ y <> 0 }} + { m:Z | n = m * m }.

となります。

方針

ペル方程式の解の存在を証明する方法は色々あると思いますが、今回は比較的Coqで表現しやすいと思われる次の証明を使います。

証明

$f:Z^4 \rightarrow Z^4$
$M2T:Z^4 \rightarrow Z^3$
$h:Z^3 \rightarrow Z^3$
を以下で定義します。
$M2T(a,b,c,d)=(a^2-nb^2,ac-nbd,c^2-nd^2)$
$h(s,t,u)=\begin{cases}(s+2t+u,t+u,u)&(s+2t+u\gt 0)\\(s,s+t,s+2t+u)&(s+2t+u\lt 0)\end{cases}$
$f(a,b,c,d)=\begin{cases}(a+c,b+d,c,d)&(s+2t+u\gt 0)\\(a,b,a+c,b+d)&(s+2t+u\lt 0)\end{cases}$
ただし、最後の式のs,t,uはM2T(a,b,c,d)=(s,t,u)です。

このとき実はM2T(f(a,b,c,d))=h(M2T(a,b,c,d))となっていることが分かります。
ここで例として、n=10の場合にf^k(1,0,0,1)とM2T(f^k(1,0,0,1))=h^k(M2T(1,0,0,1))を計算してみます。
k,f^k(1,0,0,1),M2T(f^k(1,0,0,1))=h^k(1,0,-10)を各行ごとに書いています。
0, (1,0,0,1), (1,0,-10)
1, (1,0,1,1), (1,1,-9)
2, (1,0,2,1), (1,2,-6)
3, (1,0,3,1), (1,3,-1)
4, (4,1,3,1), (6,2,-1)
5, (7,2,3,1), (9,1,-1)
6, (10,3,3,1), (10,0,-1)
7, (13,4,3,1), (9,-1,-1)
8, (16,5,3,1), (6,-2,-1)
9, (19,6,3,1), (1,-3,-1)
10, (19,6,22,7), (1,-2,-6)
11, (19,6,41,13), (1,-1,-9)
12, (19,6,60,19), (1,0,-10)
13, (19,6,79,25), (1,1,-9)
...

一番右の数字の三つ組みに注目するとt*t-s*u=n, s>0, u<0 が成り立っていることが分かります。
このようなs,t,uの組みは有限通りしかありません。また、hはこのような(s,t,u)の組み全体の集合に対する自己全単射です。
よって、右側に現れる三つ組みは途中で最初の(1,0,-n)に戻ってきてループするということです。（今回の例ではk=12）。
(s,t,u)=(1,0,-n)のときの(a,b,c,d)の値はM2T(a,b,c,d)=(1,0,-n)を満たすので、特にa*a-n*b*b=1となります。
(今回の例では(a,b,c,d)=(19,6,60,19) )。
最後に示さなくてはいけないことはb!=0となることですが、左側の数字は大きくなり続けるので、k=0を除いてb=0とはなりえません（後で詳しく書きます）。

coqで証明を書く

実際に上の方針をためにCoqで証明する主な内容は必要があります。

Aが有限集合でxがAの元で、f:A→Aが単射ならあるk>0が存在してf^k(x)=x
t*t-s*u=n, s>0, u<0 を満たすs,t,uは有限個
k!=0ならb!=0

実は２と３が一番苦労した部分です。

まず、Coqで集合を現すのに、A:SetとP:A→Propの組みであらわすことにしました。Pが成り立つようなAの元全体を(A,P)で表すって事です。
１の主張を書くため、まず次の定義を書きます。

Definition injective_sub {A}{B} (P:A->Prop) (f:A->B) := forall (x y:A), P x -> P y -> f x = f y -> x=y. 
Definition finite_sub (X:Set) (P:X -> Prop) :={f:X -> nat &{n | injective_sub P f /\ forall x,  f x < n  } }.
Definition well_def{A}{B} (f:A ->B)(P:A -> Prop) (Q:B -> Prop) := forall x:A , P x -> Q (f x).

injecttive_sub P f はfが単射となることを表します。
finite_sub X P は(X,P)が有限集合となることを表します。
well_def f P Q はfがWell-definedな写像となることを表します。
このあたりの部分はStandard-LibraryのCoq.Setsあたりを使ったほうが賢明だったかもしれないです。

Hypothesis hy1: finite_sub X P .
Hypothesis hy2: well_def f P P.
Hypothesis hy3: injective_sub P f.
Hypothesis hy4: P x0.

↑の仮定の下で次を示します。

Fixpoint iterate{A} (n:nat) (f:A->A) (x:A) :=match n with
| 0 => x
| S m => f (iterate m f x)
end.

Lemma loop_sub: {n:nat | 0<n /\ iterate n f x0 = x0}.

↑を示すために使ったのが以下の定義と補題です。

Definition skip m n:=
if (le_lt_dec m n) then S n else n.

Definition skip_inv m n:=
if (le_lt_dec m n) then (match n with | 0 => 0 | S k => k end) else n.

Lemma skip_id:forall m n,skip_inv m (skip m n) = n.

Lemma skip_lem:forall m n,m <> skip m n.

Lemma php0:forall (n:nat) (f:nat -> nat) , {m | n <= f m} + {x:nat & {y| x<>y /\ f x=f y}}.

Lemma php:forall (n:nat) (f:nat -> nat) , {m | n <= f m} + {x:nat & {y| x<y /\ f x=f y}}.

Lemma loop_sub_iter: forall m, P (iterate m f x0).

Lemma loop_sub0:{m:nat &{ n | m<n /\ iterate m f x0 = iterate n f x0}}.

まず iterate m f x0 = iterate n f x0 となるm

Let reg_t (t:Tri):Prop:=match t with 
|tri s t u => s>0 /\ u<0 /\ t*t-s*u=n
end.

示したいことは次の補題です。

Lemma psub1:finite_sub Tri reg_t.

実はこのLemmaの証明だけに３００行近くかかっています。
このLemmaの証明は主に２パートに分かれていて、第一パートで有限集合の積の有限性を証明し(コードでいうとSection P2Nの内部)、第二パートでそれを使って、不等式を処理したりしながら有限性を示しています。証明自体はかなり冗長で気合で乗り切った部分です。
あまりうまい方法ではないですが、うまい方法を求めてぐだぐだ悩み続けるよりは、冗長な方針でもいいからとにかく書いてみるのがいいと思います。

最後に３の部分です。

右側の数字(s,t,u)が(1,0,-n)に戻ったとき、左側の数字(a,b,c,d)が(a,b)!=(1,0)となっていることを示す必要があります。
自明だと思えるような主張ですが、(a,b)!=(1,0)となる理由を書き出すと次のようなり、意外とステップ数が多いことが分かります。

(a,b,c,d)は常にa>=1,b>=0,c>=0,d>=1である。またk!=0を除いてc>=1である。
(a,b,c,d)->(a,b,a+c,b+d)の変換のみを繰り返していると、(s,t,u)のtが大きくなり続けて途中で矛盾する。
よって、途中で(a,b,c,d)->(a+c,b+d,c,d)の変換が必ず入る。
↑の変換後は常にb>=1を満たす。

というわけでb!=1を示すのは思ったより面倒です。一つ一つ丁寧に証明していけば証明が完了します。

まとめ

人に説明するときなら自明の一言で済ませたくなる部分でも、Coqで証明するには非常に面倒だったりする。
逆に、数学的に一番キーポイントだと思える部分がCoqで書くとほんの数行だったりする。

コード

最後に少し長いですが、Coqコードを載せておきます。anarchy proof に投稿したものに「若干の変更＋Show Scriptコマンドによる整形」を行ったものです。

Require Import ZArith Znumtheory Arith.
Open Scope Z_scope.

(* facts about finite set *)
Open Scope nat_scope.


Fixpoint iterate{A} (n:nat) (f:A->A) (x:A) :=match n with
| 0 => x
| S m => f (iterate m f x)
end.

Definition skip m n:=
if (le_lt_dec m n) then S n else n.
Definition skip_inv m n:=
if (le_lt_dec m n) then (match n with | 0 => 0 | S k => k end) else n.

Lemma skip_id:forall m n,skip_inv m (skip m n) = n.
intros.
unfold skip in |- *.
destruct le_lt_dec.
 unfold skip_inv in |- *.
 destruct le_lt_dec.
  auto.
  
  omega.
  
 unfold skip_inv in |- *.
 destruct le_lt_dec.
  omega.
  
  auto.
Qed.

Lemma skip_lem:forall m n,m <> skip m n.
intros.
unfold skip in |- *.
destruct le_lt_dec; omega.
Qed.

Lemma php0:forall (n:nat) (f:nat -> nat) , {m | n <= f m} + {x:nat & {y| x<>y /\ f x=f y}}.
induction n.
 left.
 exists 0.
 omega.
 
 intro.
 case (IHn f).
  intro.
  destruct s.
  destruct (eq_nat_dec (f x) n).
   destruct (IHn (fun k => f (skip x k))); destruct s.
    destruct (eq_nat_dec (f (skip x x0)) n).
     right.
     exists x.
     exists (skip x x0).
     split.
      apply skip_lem.
      
      omega.
      
     left.
     exists (skip x x0).
     omega.
     
    destruct s.
    destruct a.
    right.
    exists (skip x x0).
    exists (skip x x1).
    split.
     intro.
     assert (x0 = x1).
      rewrite <- (skip_id x x0) in |- *.
      rewrite <- (skip_id x x1) in |- *.
      rewrite H1 in |- *.
      auto.
      
      auto.
      
     auto.
     
   left.
   exists x.
   omega.
   
  auto.
Qed.

Lemma php:forall (n:nat) (f:nat -> nat) , {m | n <= f m} + {x:nat & {y| x<y /\ f x=f y}}.
intros.
destruct (php0 n f).
 auto.
 
 destruct s as (x, (y, (H1, H2))).
 right.
 destruct (le_lt_dec x y).
  exists x.
  exists y.
  omega.
  
  exists y.
  exists x.
  omega.
Qed.

Definition injective_sub {A}{B} (P:A->Prop) (f:A->B) := forall (x y:A), P x -> P y -> f x = f y -> x=y.
Definition finite_sub (X:Set) (P:X -> Prop) :={f:X -> nat &{n | injective_sub P f /\ forall x,  f x < n  } }.
Definition well_def{A}{B} (f:A ->B)(P:A -> Prop) (Q:B -> Prop) := forall x:A , P x -> Q (f x).

Section LOOP.

Variable (X:Set) (f:X->X) (x0:X) (P:X -> Prop).
Hypothesis hy1: finite_sub X P .
Hypothesis hy2: well_def f P P.
Hypothesis hy3: injective_sub P f.
Hypothesis hy4: P x0.

Lemma loop_sub_iter: forall m, P (iterate m f x0).
induction m; simpl in |- *; auto.
Qed.

Lemma loop_sub0:{m:nat &{ n | m<n /\ iterate m f x0 = iterate n f x0}}.
destruct hy1 as (inj, (upper, (H1, H2))).
clear hy1 hy2.
destruct (php upper (fun n => inj (iterate n f x0))).
 destruct s as (m, H3).
 assert (inj (iterate m f x0) < upper).
  apply H2.
  
  exists 0.
  exists 0.
  apply False_ind.
  omega.
  
 destruct s as (x, (y, (H3, H4))).
 exists x.
 exists y.
 split.
  auto.
  
  apply H1; try apply loop_sub_iter.
  auto.
Qed.

Lemma minus:forall m n,n<=m -> {k | m=n+k}.
intros.
exists (m - n).
omega.
Qed.

Lemma loop_sub: {n:nat | 0<n /\ iterate n f x0 = x0}.
destruct loop_sub0 as (m1, (m2, (H1, H2))).
assert (forall n, n <= m1 -> iterate (m1 - n) f x0 = iterate (m2 - n) f x0).
 induction n.
  intro.
  assert (forall mm, mm - 0 = mm).
   intro.
   omega.
   
   rewrite H0,  H0 in |- *.
   auto.
   
  intro.
  apply hy3; try apply loop_sub_iter.
  assert (forall mm, S n <= mm -> mm - n = S (mm - S n)).
   intro.
   omega.
   
   rewrite H0,  H0 in IHn; try omega.
   simpl in IHn.
   apply IHn.
   omega.
   
 exists (m2 - m1).
 split.
  omega.
  
  rewrite <- (H m1) in |- *.
   replace (m1 - m1) with 0  by omega.
   auto.
   
   omega.
Qed.

End LOOP.

Check loop_sub.

Lemma p_fin_inj:forall X Y P Q f ,finite_sub Y Q -> well_def f P Q -> injective_sub P f -> finite_sub X P.
unfold finite_sub in |- *.
unfold injective_sub in |- *.
unfold well_def in |- *.
intros.
destruct H.
exists (fun xx => x (f xx)).
destruct s.
exists x0.
destruct a.
split.
 intros.
 apply H1.
  auto.
  
  auto.
  
  apply H.
   apply H0.
   auto.
   
   apply H0.
   auto.
   
   auto.
   
 intro.
 apply H2.
Qed.

Definition prod_prop{A}{B}(P:A->Prop)(Q:B->Prop)(z:(A*B)%type):= let (x,y) := z in P x /\ Q y.

Section aP2N.
Variable xs:nat.
Definition P2N x y:=x+y*xs.

Lemma p_p2n1 x1 x2 y1 y2:x1<xs -> x2<xs -> y1<y2 -> P2N x1 y1 < P2N x2 y2.
intros.
unfold P2N in |- *.
assert {dy : _ | y2 = dy + 1 + y1}.
 exists (y2 - y1 - 1).
 omega.
 
 destruct H2.
 rewrite e in |- *.
 replace ((x + 1 + y1) * xs) with (x * xs + y1 * xs + xs)  by ring.
 assert (0 <= x * xs).
  apply le_O_n.
  
  omega.
Qed.

Lemma p_p2n3 x1 x2 y1 y2:x1<xs -> x2<xs -> P2N x1 y1 =P2N x2 y2 -> y1=y2.
intros.
destruct (lt_eq_lt_dec y1 y2).
 destruct s.
  assert (P2N x1 y1 < P2N x2 y2).
   apply p_p2n1.
    auto.
    
    omega.
    
    auto.
    
   omega.
   
  auto.
  
 assert (P2N x2 y2 < P2N x1 y1).
  apply p_p2n1.
   auto.
   
   omega.
   
   auto.
   
  omega.
Qed.

Lemma p_p2n4 x1 x2 y1 y2 :x1<xs -> x2 <xs ->
 P2N x1 y1 = P2N x2 y2 -> x1 =x2 /\ y1=y2.
intros.
assert (y1 = y2).
 apply (p_p2n3 x1 x2 y1 y2); omega.
 
 split.
  unfold P2N in H1.
  rewrite H2 in H1.
  omega.
  
  auto.
Qed.

Definition P2Nprod z:=match z with (x,y) => P2N x y end.
Definition regular_prod (z:nat*nat):=match z with (x,y) => x<xs end.

Lemma  p_p2n5 z1 z2:regular_prod z1 -> regular_prod z2 -> P2Nprod z1 = P2Nprod z2 -> z1=z2.
intros z1 z2 H1 H2 H.
destruct z1 as (x1, y1).
destruct z2 as (x2, y2).
unfold regular_prod in *.
assert (x1 = x2 /\ y1 = y2).
 apply p_p2n4; auto.
 
 destruct H0.
 auto.
Qed.

End aP2N.

Check P2Nprod.

Lemma p_fin_prod:forall X Y P Q , finite_sub X P -> finite_sub Y Q -> 
 finite_sub ( (X*Y)%type ) (prod_prop P Q).
intros.
destruct H as (jx, (xs, (HX1, HX2))).
destruct H0 as (jy, (ys, (HY1, HY2))).
exists (fun z => match z with
                 | pair x y => P2Nprod xs (jx x, jy y)
                 end).
exists (ys * xs).
split.
 unfold injective_sub in |- *.
 intros z1 z2 H1 H2 H.
 destruct z1 as (x1, y1).
 destruct z2 as (x2, y2).
 assert (x1 = x2 /\ y1 = y2).
  assert ((jx x1, jy y1) = (jx x2, jy y2)).
   apply (p_p2n5 xs); simpl in |- *; auto.
   
   inversion H0.
   split; [ apply HX1 | apply HY1 ]; destruct H1; destruct H2; auto.
   
  destruct H0.
  rewrite H0,  H3 in |- *.
  auto.
  
 intro.
 destruct x as (x, y).
 unfold P2Nprod in |- *.
 unfold P2N in |- *.
 assert (jx x < xs).
  apply HX2.
  
  assert (jy y < ys).
   apply HY2.
   
   assert (S (jy y) * xs <= ys * xs).
    apply mult_le_compat_r; omega.
    
    simpl in H1.
    omega.
Qed.

Open Scope Z_scope.

Inductive perfect_square : Z -> Prop :=
  ps_intro : forall n, perfect_square (n * n).

Lemma perfect_square_dec:
 forall n, { m | n=m*m} + {~perfect_square n}.
intros.
destruct (Z_eq_dec (Zsqrt_plain n * Zsqrt_plain n) n).
 left.
 rewrite <- e in |- *.
 exists (Zsqrt_plain n).
 auto.
 
 right.
 intro.
 destruct H.
 destruct (Z_lt_dec 0 n).
  rewrite Zsqrt_square_id in n0.
   omega.
   
   omega.
   
  replace (n * n) with (- n * - n) in n0  by ring.
  rewrite Zsqrt_square_id in n0.
   omega.
   
   omega.
Qed.

Lemma sqpos:forall n, 0<>n -> 0<n*n.
intros.
destruct (Z_lt_dec 0 n).
 apply Zmult_lt_0_compat; auto.
 
 replace (n * n) with (- n * - n)  by ring.
 apply Zmult_lt_0_compat; omega.
Qed.


Section A.
Variable n:Z.
Hypothesis pos:n>0.
Hypothesis psn:~(perfect_square n).

Inductive Tri: Set :=
| tri (s t u:Z) :Tri.

Inductive Mat: Set :=
|mat (a b c d:Z):Mat.

Let M2T (m:Mat) := match m with
|mat a b c d => tri (a*a-n*b*b) (a*c-n*b*d) (c*c-n*d*d)
end.

Let f (m:Mat) := match m with
|mat a b c d =>match M2T(m) with
|tri s t u =>
if (Z_lt_dec 0 (s+t+t+u)) then mat (a+c) (b+d) c d else mat a b (a+c) (b+d) 
end
end.

Let reg_t (t:Tri):Prop:=match t with 
|tri s t u => s>0 /\ u<0 /\ t*t-s*u=n
end.


Definition h t:=match t with
|tri s t u => if (Z_lt_dec 0 (s+t+t+u)) then tri (s+t+t+u) (t+u) u else tri s (s+t) (s+t+t+u)  
end.

Definition h_inv t:=match t with
|tri s t u => if (Z_lt_dec 0 (s-t-t+u)) then tri (s-t-t+u) (t-u) u else tri s (-s+t) (s-t-t+u)
end.

Definition Z2:=(Z*Z)%type.
Definition Z3:=(Z*Z*Z)%type.
Definition prZ z:=0<=z<=n+n.
Definition prZ2 z2:=match z2 with
| (a,b) => 0<=a<=n+n /\ 0<=b<=n+n
end.
Definition prZ3 z3:=match z3 with
| (a,b,c) => 0<a<=n /\ 0<=b<=n+n /\ 0<c<=n
end.
Definition prZ3em t :=match t with
| tri s t u =>(s,t+n,-u)
end.


Lemma p_fin_prod2:forall (X Y:Set) (P:X -> Prop) (Q:Y->Prop) (R:(X*Y)%type -> Prop),
  finite_sub X P -> finite_sub Y Q -> (forall x y , R (x,y) -> P x /\ Q y)  -> finite_sub (X*Y) R.
intros.
assert (finite_sub (X * Y) (prod_prop P Q)).
 apply p_fin_prod; auto.
 
 unfold finite_sub in |- *.
 unfold finite_sub in H2.
 destruct H2 as (f00).
 destruct s as (n00).
 exists f00.
 exists n00.
 split.
  destruct a.
  unfold injective_sub in |- *.
  unfold injective_sub in H2.
  unfold prod_prop in H2.
  intros r1 r2.
  intros.
  destruct r1 as (x1, y1).
  destruct r2 as (x2, y2).
  apply H2.
   apply H1; auto.
   
   apply H1; auto.
   
   auto.
   
  apply a.
Qed.

Lemma pzn:forall n,0<=n -> Z_of_nat (Zabs_nat n) = n.
intros.
destruct n0.
 simpl in |- *.
 auto.
 
 simpl in |- *.
 rewrite Zpos_eq_Z_of_nat_o_nat_of_P in |- *.
 auto.
 
 assert (Zneg p < 0).
  apply Zlt_neg_0.
  
  omega.
Qed.


Definition embzn zup z:=
if (Z_le_dec z zup) then (if Z_le_dec 0 z then Zabs_nat z else 0%nat )else 0%nat.

Definition ZNP:=Zpos_eq_Z_of_nat_o_nat_of_P.

Lemma psubz:finite_sub Z prZ.
unfold prZ in |- *.
unfold finite_sub in |- *.
assert (0 < n + n).
 omega.
 
 revert H.
 generalize (n + n).
 intros.
 exists (embzn z).
 exists (1 + Zabs_nat z)%nat.
 unfold injective_sub in |- *.
 split.
  intros.
  unfold embzn in H2.
  destruct Z_le_dec; destruct Z_le_dec; try omega.
  destruct Z_le_dec; try destruct Z_le_dec; try omega.
  rewrite <- pzn in |- *.
   rewrite <- (pzn x) in |- *.
    auto.
    
    omega.
    
   omega.
   
  unfold embzn in |- *.
  intro.
  destruct Z_le_dec; try destruct Z_le_dec; try omega.
  unfold Zabs_nat in |- *.
  destruct x.
   omega.
   
   destruct z.
    omega.
    
    apply inj_lt_rev.
    rewrite inj_plus in |- *.
    rewrite <- ZNP in |- *.
    rewrite <- ZNP in |- *.
    replace (Z_of_nat 1) with 1 .
     omega.
     
     auto.
     
    assert (Zneg p0 < 0).
     apply Zlt_neg_0.
     
     omega.
     
   assert (Zneg p < 0).
    apply Zlt_neg_0.
    
    omega.
Qed.


Lemma psub10:finite_sub Z3 prZ3.
apply (p_fin_prod2 Z2 Z prZ2 prZ).
 apply (p_fin_prod2 Z Z prZ prZ).
  exact psubz.
  
  exact psubz.
  
  unfold prZ2 in |- *.
  unfold prZ in |- *.
  intros.
  omega.
  
 exact psubz.
 
 destruct x.
 unfold prZ3 in |- *.
 unfold prZ2 in |- *.
 unfold prZ in |- *.
 intros.
 omega.
Qed.

Lemma sqpos2:forall x,0<=x*x.
intro.
destruct (Z_zerop x).
 rewrite e in |- *.
 simpl in |- *.
 omega.
 
 assert (0 < x * x).
  apply sqpos.
  omega.
  
  omega.
Qed.

Lemma sql0:forall x , 0<=x -> x*x<=n -> -n<=x<=n.
intros.
destruct (Z_zerop x).
 omega.
 
 destruct (Z_le_dec x n).
  omega.
  
  assert (x * n < x * x).
   apply Zmult_lt_compat_l; omega.
   
   assert (1 * n <= x * n).
    apply Zmult_le_compat_r; omega.
    
    omega.
Qed.

Lemma sql1:forall x , x*x<=n -> -n<=x<=n.
intros.
destruct (Z_lt_dec x 0).
 cut (- n <= - x <= n).
  intro.
  omega.
  
  apply sql0.
   omega.
   
   replace (- x * - x) with (x * x)  by ring.
   omega.
   
 apply sql0; omega.
Qed.

Lemma sql2:forall x y,0<x -> 0<y -> x*y<=n -> x<=n.
intros.
destruct (Z_le_dec x n).
 omega.
 
 assert (n * 1 <= n * y).
  apply Zmult_le_compat_l; omega.
  
  assert (n * y < x * y).
   apply Zmult_lt_compat_r; omega.
   
   omega.
Qed.


Lemma psub1:finite_sub Tri reg_t.
apply (p_fin_inj Tri Z3 reg_t prZ3 prZ3em).
 exact psub10.
 
 unfold well_def in |- *.
 intros.
 destruct x.
 unfold reg_t in H.
 unfold prZ3em in |- *.
 unfold prZ3 in |- *.
 destruct H.
 assert (0 <= t * t).
  apply sqpos2.
  
  assert (s * u < 0).
   assert (0 < s * - u).
    apply Zmult_lt_0_compat; omega.
    
    replace (s * u) with (- (s * - u))  by ring.
    omega.
    
   assert (- n <= t <= n).
    apply sql1.
    omega.
    
    assert (s * - u <= n).
     replace (s * - u) with (- (s * u))  by ring.
     omega.
     
     split.
      assert (s <= n).
       apply (sql2 s (- u)); omega.
       
       omega.
       
      split.
       omega.
       
       split.
        omega.
        
        apply (sql2 (- u) s).
         omega.
         
         omega.
         
         replace (- u * s) with (s * - u)  by ring.
         omega.
         
 unfold injective_sub in |- *.
 intros.
 unfold prZ3em in H1.
 destruct x.
 destruct y.
 injection H1.
 intros.
 assert (s = s0 /\ t = t0 /\ u = u0) by complete omega.
 destruct H5.
 destruct H6.
 rewrite H5,  H6,  H7 in |- *.
 auto.
Qed.


Lemma psub2:well_def h reg_t reg_t.
unfold well_def in |- *.
unfold reg_t in |- *.
intros.
destruct x.
unfold h in |- *.
destruct H.
destruct H0.
destruct Z_lt_dec; split.
 omega.
 
 split.
  omega.
  
  rewrite <- H1 in |- *.
  ring.
  
 omega.
 
 split.
  destruct (Z_zerop (s + t + t + u)).
   assert (s = - t - t - u).
    omega.
    
    rewrite H2 in H1.
    apply False_ind.
    apply psn.
    replace n with ((t + u) * (t + u)) .
     apply ps_intro.
     
     rewrite <- H1 in |- *.
     ring.
     
   omega.
   
  rewrite <- H1 in |- *.
  ring.
Qed.

Lemma pp_inj_tri :forall s1 s2 t1 t2 u1 u2 , s1=s2 -> t1=t2 -> u1=u2 -> tri s1 t1 u1 = tri s2 t2 u2.
intros.
rewrite H,  H0,  H1 in |- *.
auto.
Qed.

Lemma h_id:forall t, (reg_t t) -> h_inv (h t) = t.
intros.
destruct t.
unfold h in |- *.
unfold h_inv in |- *.
unfold reg_t in H.
destruct Z_lt_dec.
 destruct Z_lt_dec.
  apply pp_inj_tri; omega.
  
  apply False_ind; omega.
  
 destruct Z_lt_dec.
  apply pp_inj_tri; omega.
  
  apply pp_inj_tri; omega.
Qed.

Lemma psub3:injective_sub reg_t h.
unfold injective_sub in |- *.
intros.
rewrite <- h_id in |- *.
 rewrite <- (h_id x) in |- *.
  rewrite H1 in |- *.
  auto.
  
  auto.
  
 auto.
Qed.


Lemma loop_h:{ k:nat | (0<k)%nat /\ iterate k h (tri 1 0 (-n) ) = (tri 1 0 (-n))}.
apply loop_sub with reg_t.
 exact psub1.
 
 exact psub2.
 
 exact psub3.
 
 unfold reg_t in |- *.
 split.
  omega.
  
  split.
   omega.
   
   ring.
Qed.

Definition K:=(proj1_sig loop_h).
Definition Finalmat:=iterate K f (mat 1 0 0 1).

Definition X:=match Finalmat with
| mat a b c d => a
end.

Definition Y:=match Finalmat with
|mat  a b c d => b
end.

Definition X2:=match Finalmat with
|mat a b c d => c
end.

Definition Y2:=match Finalmat with
|mat a b c d => d
end.


Lemma loopK : (0<K)%nat /\ iterate K h (tri 1 0 (-n)) = tri 1 0 (-n).
exact (proj2_sig loop_h).
Qed.

Lemma loopK_1: (0<K)%nat.
apply loopK.
Qed.

Lemma loopK_2:iterate K h (tri 1 0 (-n)) = tri 1 0 (-n).
apply loopK.
Qed.

Lemma eq_Finalmat:mat X Y X2 Y2 = Finalmat.
unfold X, Y, X2, Y2 in |- *.
destruct Finalmat.
auto.
Qed.

Lemma pp1 m:M2T(f(m)) = h( M2T(m)).
destruct m.
unfold f in |- *.
unfold h in |- *.
unfold M2T in |- *.
destruct Z_lt_dec; apply pp_inj_tri; ring.
Qed.

Lemma pp2:M2T Finalmat = tri 1 0 (-n).
rewrite <- loopK_2 in |- *.
unfold Finalmat in |- *.
generalize K.
induction n0.
 unfold iterate in |- *.
 unfold M2T in |- *.
 replace (1 * 1 - n * 0 * 0) with 1  by ring.
 replace (1 * 0 - n * 0 * 1) with 0  by ring.
 replace (0 * 0 - n * 1 * 1) with (- n)  by ring.
 auto.
 
 simpl in |- *.
 rewrite <- IHn0 in |- *.
 apply pp1.
Qed.


Definition reg1_m m:=match m with
|mat a b c d =>a>0 /\ b>0 /\ c>=0 /\ d>0
end.

Definition reg2_m m:=match m with
|mat a b c d =>a>0 /\ b>=0 /\ c>0 /\ d>0
end.

Definition reg3_m m:={reg1_m m} + {reg2_m m}.

Lemma f_dec:forall a b c d, {f (mat a b c d) = mat (a+c) (b+d) c d} + {f (mat a b c d) = mat a b (a+c) (b+d)}.
intros.
unfold f in |- *.
destruct M2T.
destruct (Z_lt_dec 0 (s + t + t + u)).
 left.
 auto.
 
 right.
 auto.
Qed.
 
Lemma pp3_1: forall m,reg1_m m -> reg1_m (f m).
destruct m.
destruct (f_dec a b c d); rewrite e in |- *; unfold reg1_m in |- *; omega.
Qed.

Lemma pp3_2: forall m,reg2_m m -> reg2_m (f m).
destruct m.
destruct (f_dec a b c d); rewrite e in |- *; unfold reg2_m in |- *; omega.
Qed.

Lemma pp3_3:forall m,reg3_m m -> reg3_m (f m).
intros.
unfold reg3_m in |- *.
destruct H.
 left.
 apply pp3_1.
 auto.
 
 right.
 apply pp3_2.
 auto.
Qed.

Lemma pp3_4:reg3_m Finalmat.
unfold Finalmat in |- *.
assert (0 < K)%nat.
 exact loopK_1.
 
 induction K.
  left.
  apply False_ind.
  omega.
  
  destruct (zerop n0).
   rewrite e in |- *.
   unfold iterate in |- *.
   destruct (f_dec 1 0 0 1); rewrite e0 in |- *.
    left.
    unfold reg1_m in |- *.
    omega.
    
    right.
    unfold reg2_m in |- *.
    omega.
    
   simpl in |- *.
   apply pp3_3.
   auto.
Qed.


Lemma pp4:   1=(X*X - n*Y*Y) /\ 0 = (X*X2 -  n*Y*Y2) /\ (-n) =(X2*X2 - n*Y2*Y2) .
assert
 (tri 1 0 (- n) =
  tri (X * X - n * Y * Y) (X * X2 - n * Y * Y2) (X2 * X2 - n * Y2 * Y2)).
 rewrite <- pp2 in |- *.
 rewrite <- eq_Finalmat in |- *.
 unfold M2T in |- *.
 auto.
 
 injection H.
 auto.
Qed.

Lemma pp3_5:X<>1 /\ Y<>0.
destruct pp4.
destruct H0.
destruct pp3_4; rewrite <- eq_Finalmat in r.
 unfold reg1_m in r.
 split.
  destruct (Z_eq_dec X 1).
   rewrite e in H.
   replace (1 * 1) with 1 in H  by ring.
   assert (n * Y * Y = 0).
    omega.
    
    replace (n * Y * Y) with (n * (Y * Y)) in H2  by ring.
    assert (n = 0).
     apply (Zmult_integral_l (Y * Y)).
      assert (0 < Y * Y).
       apply sqpos.
       omega.
       
       omega.
       
      auto.
      
     omega.
     
   auto.
   
  omega.
  
 unfold reg2_m in r.
 destruct (Z_zerop Y).
  rewrite e in H0.
  replace (n * 0 * Y2) with 0 in H0  by ring.
  assert (X = 0).
   apply (Zmult_integral_l X2); omega.
   
   omega.
   
  split.
   destruct (Z_eq_dec X 1).
    rewrite e in H.
    replace (1 * 1) with 1 in H  by ring.
    assert (n * Y * Y = 0).
     omega.
     
     replace (n * Y * Y) with (n * (Y * Y)) in H2  by ring.
     assert (n = 0).
      apply (Zmult_integral_l (Y * Y)).
       assert (0 < Y * Y).
        apply sqpos.
        omega.
        
        omega.
        
       auto.
       
      omega.
      
    auto.
    
   auto.

Qed.


Lemma pell_eq:{ x:Z & { y:Z | x*x - n*y*y = 1 /\ x <> 1 /\ y <> 0 }}.
exists X.
exists Y.
split.
 symmetry  in |- *.
 apply pp4.
 
 apply pp3_5.
Qed.

End A.

Theorem Pell's_equation:
  forall (n: Z), n > 0 -> { x:Z & { y:Z | x*x - n*y*y = 1 /\ x <> 1 /\ y <> 0 }} + { m:Z | n = m * m }.
Proof.
intros.
destruct (perfect_square_dec n).
 auto.
 
 left.
 apply pell_eq.
  auto.
  
  auto.
Qed.

2011-11-22

数列の推測

Project Eulerなどの問題を解いていると、「mod による場合分けを含んだ多項式により定義された数列」なんだろうけど、一般項を求めるのは面倒だなって思うことが良くある気がします。そこで、そのようなタイプの数列の一般項をすばやく求めるためのプログラムをpythonで書いてみました。
使い方は簡単です。a_nを求めたい場合はf(a,n)がa_nを返してくれます。


class Est:
  def __init__(self):
    self.mz = None
    pass

  def build(self,x):
    n = len(x)
    for k in xrange(1,n):
      z=[x[i::k] for i in xrange(k)]
      mz =[ self.f(x2) for x2 in z ]
      if all( len(mz[i])+2 <= len(z[i])  for i in xrange(k)):
        break
    else:
      print "Failed!"
    self.mz=mz
    return

  def get(self,n):
    k= n % len(self.mz)
    n0 = n / len(self.mz)
    return self.g( self.mz[k] , n0 )

  def f(self,_x):
    x=list(_x)
    m=[]
    while True:
      if all(i==0 for i in x):break
      m.append(x[0])
      x=[ x[i+1]-x[i] for i in xrange(len(x)-1) ]
    return m

  def g(self, m, n):
    s=1
    ret=0
    for i,v in enumerate(m):
      ret += v*s
      s = s*(n-i)/(i+1)
    return ret

def f(x,n):
  est=Est()
  est.build(x)
  return est.get(n)

def test():
  z=[(i**3,i+1,i-1)[i%3] for i in xrange(20)]
  print f(z, 300) # 27000000 = 30**3
  print f(z,301) # 301+1
  print f(z,302) # 302-1

test()

2011-09-14

Codeforces86E(Sleeping)

http://www.codeforces.com/problemset/problem/113/E
参加しなかった回の問題ですが練習で解いてみました。
この手の問題ってどうやって解くのが一番いいのか悩みます。
とりあえず、再帰をフルに使って解いてみました。

案の定間違えまくって、４回目ぐらいでやっとAcceptになりました。
やっぱりtestを書かなきゃ通すのは不可能なのかもしれません。



def readints():
  return map(int, raw_input().split())


def dn(x):
  k=0
  t=x-1
  while t:
    t/=10
    k+=1
  return k



def _li(x, dn):
  ret=[]
  for _ in xrange(dn):
    ret.append(x%10)
    x/=10
  return ret

def li(h,m):
  return _li(h,dh)+_li(m,dm)

def dif(h1,m1,h2,m2):
  return sum(v1!=v2 for v1,v2 in zip(li(h1,m1), li(h2,m2)) )

pow10=[10**i for i in xrange(12,-1,-1)]

memo={}

def g(h,m):
  if (h,m) in memo:return memo[(h,m)]
  ret=f(0,1,h,m)+(1>=k)
  memo[(h,m)]=ret
  return ret


def f(h1,m1,h2,m2):
  if (h1,m1)==(h2, m2):ret = 0
  elif (h1,m1)>(h2,m2):
    ret = f(h1, m1, H-1, M-1) + f(0, 0, h2, m2) + ( dif(H-1,M-1,0,0) >= k)
  else:
    for hh in pow10:
      if (h1+hh, m1)<=(h2, m2) and h1/hh%10!=9:
        ret = g(hh,0)+f(h1+hh, m1, h2, m2)
        break
    else:
      for mm in pow10:
        if m1+mm<M and (h1,m1+mm)<=(h2,m2) and m1/mm%10!=9:
          ret = g(0,mm) + f(h1,m1+mm, h2, m2)
          break
      else:
        nt=h1*M+m1+1
        nh,nm=nt/M,nt%M
        ret = ( dif(h1,m1,nh,nm)>=k ) + f(nh,nm,h2,m2)
  #print h1,m1,h2,m2,ret
  return ret

def main():
  global H,M,k,dh,dm
  H,M,k=readints()
  h1,m1=readints()
  h2,m2=readints()
  """H,M,k=24,60,1
  h1,m1=0,0
  h2,m2=23,59
  H,M,k=24,60,3
  h1,m1=23,59
  h2,m2=23,59"""
  dh=dn(H)
  dm=dn(M)
  print f(h1,m1,h2,m2)

main()

2011-06-15

TopCoderOpen 2011 Marathon Roun2

今回は頑張りが足りなかったです。最終順位はおそらく５０位ぐらいですが、実は通過の１００位以内も結構危なかった。
以下、参加したときのメモ。特に纏めてはいないのでぐちゃぐちゃです。CodeChefのJune Long Contestにも参加してたので、それも混じってます。

２０１１・６・１

19:37
とりあえずCodeChefのCLONEを解いてる。これは楽に解けそう。
２１：０７
やっとCLONEを通した。最後に%modをするのを忘れて何回かWAを食らった。一番乗りでちょっと嬉しい。
１１：５２
CTEAMSを送信。pythonでかくよりCで書いたほうが楽だったかも。しかもRAをくらってる。
２２：０４
訳が分からないが、最初のN=input()でRuntimeErrorを食らってるらしい。謎だ。
２２：１６
N=int(raw_input())にしたら通った。いくらなんでも謎過ぎる。

疲れたのでちょっと休憩。

２０１１・６・２

１：０４
Maximal Score Pathをサブミット。無事一発で通過。これで三問。

２：１７
marathon matchの問題。第一印象。幾何か・・・。
sizeの定義が書いてない気がする・・・。面積のことでいいのだろうか。
面積だとしたら、あまり辺数の多い多角形を作っても仕方ない。とりあえずは正三角形を作るべき？
得点が相対性なのがちょっと不満。そこそこ奥が深そうなので一位続出はないと信じたい。テストケースの数を増やすかもしれないし。

２：２２
どうやって正三角形を探索しよう？頂点数５０００は結構厳しい。

２：２３
test Case Generationが結構ややこしそう。実際の図をみてもかなり特徴的な分布という感じがする。

２：２７
半径１とすると三角形より四角形のほうが有利っぽい。
三角形　ルート３/2=0.4380 正方形1/2=0.50...

２：５１
正方形の探索は、誤差を気にしなければOrder(P^2)で終わりそう。

３：００
眠くなってきたので、適当に何か実装、送信して、とりあえず寝よう。

４：０１
区間分割して、正方形を作るのが妥当な作戦？

１１：２０
起床。得点の予想つかないけど、とりあえず組んでみるか。今回の問題はむずい・・・。１００位以内も実は少し怖い。

２０１１・６・３

２：１８
SRMに参加してた。そろそろ、今日寝るまでに、何か送信しておかないとずるずる行ってまずい気がする。

３：１５
javaのsubprocessとして実行すると、デバッグしにくいので、ファイルに書き出すことに。

４：０９
choosePolygonsを呼び出しただけでsegmentaton faultを引き起こす謎のバグが発生。原因が分からなくてどうしようもない。

４：２８
原因が分かった気がする。どうやら全然関係なさそうな部分に原因があるらしい。

１０：４４
起きた。

１１：０９
バグをなんとか修正して、やっと正のスコアが取れた。かなり出遅れてる感があるが大丈夫かな・・・。

１１：１３
１６点ってどういうことだろう？いやな予感がする。

１１：２０
問題文よく見たら書いてた。・・・・。やばい。sizeは面積のことじゃなかった。頂点数の二乗らしい。これは非常にまずい。

１１：２１
かなり無駄に時間を過ごしてしまった。というか戦略を１から考え直すはめになった。

１１：２５
ふ、二日間ロスっただけだ。まだ、何とかなる。しかし、この得点の決め方だと、同点続出の可能性が確かにあるなあ。面積の方がよかった気がする。

１１：２６
やっぱり、すぐにプロトタイプを作らないと危険だということが分かった。恐ろしい・・・。

１１：３４
しかし、逆にちょっと考えやすい問題になったのは確かだ。基本的にでかい多角形から順に作っていけばいいだけの気がするし。

１１：４９
図形は小さくてもいいんだよなあ。それが一番重大な事実かもしれない。

１５：４２
Test Case Generationをちゃんと読んだ。螺旋状に点が配置されるのかと思ったらそんなことは全然無かった。
意図的に正多角形を作りやすくしてる感じ。

１５：４４
点が多い場合は中心位置を推理するのがまず重要になりそうだ。
中心の個数は最大５〜２０個。

１５：４６
点の個数が５０−４９９　と、５００−５０００で、色々パターンがあるので、場合わけをするべきかも。とりあえずは点の個数が多い場合の対策をするべき。

１５：５３
明確に円っぽい物が見えるときもあるけど、必ずしもそうではないなあ。

６/4

００：３４
ニコニコみて時間つぶしてた。何やってんだ・・・。

２２：０７
pygameで極座標表示をしてみたり。
相当だらだら過ごしてるから、このままのペースでいくとまずい。

２２：０８
辺数を増やしたとき、実はradiiDiffってあんまり関係ない気がする。凸性と、sidesDiffの条件からradiiDiffは自然に大丈夫になるんじゃないかと。

２２：０９
実験してみないと分からないけど。

６/6

００：４１
中心決める方法はあんまりうまく行かなさそうなきがするなあ。

１４：２０
正多角形の位置の頂点で、一番近い点を選択するという、普通の方法を実装することにする。何故かこれを思いつくのに時間がかかった。

６/８

4:30
あと一時間以内に方針を決めて、さらに１時間以内にsubmitしよう。

４：５０
やりかた候補
完全にランダムに中心と半径を決める。
二点決めて、その中心と距離で半径を決める。
３点選んで中心と半径を決める。

６/10

謎の雑務が終わったので、今から本気出して頑張る。とりあえず駄目コードで一回submitしている状態。
１８．０３点で１４６位。clock()関数があまりよく機能していないことに気付いて衝撃を受けている。これはまずいだろ。

6/11

1:21
起きたので開始。

１：５３
clock()関数で１秒で打ち切ったら、実行時間が３〜５秒ほど。なにこれ・・・。

２：２１
clock()関数の代わりにgettimesofdayを使ったら、うまく行った。何か変な気がするけど、これで戦える。

２：４８
timerをちゃんとセットして送信したら、28.03点。１２６位。

６：０５
exampleの結果をメモ
112.0
928.0
8352.0
13040.0
864.0
8576.0
736.0
512.0
3968.0
208.0

これは二点選んで、それを対点とする正方形を探索し続ける作戦。

６：２１
時間が<5のときは辺数６
＞５のときは辺数４で
テスト

116.0
940.0
10992.0
17872.0
916.0
9828.0
696.0
496.0
5292.0
200.0

逆にSCOREが落ちているケースもある。（３つのケース）。これは意外。

１５：２３
全く同じコードでテストしてみたら、結果が不安定。これはちょっとやる気をなくす。

6/12

０：２９
正三角形も真面目にサーチしたら、それだけで割りとましな点数になりそう。

163.0
1424.0
11754.0
18614.0
1225.0
11049.0
1118.0
762.0
5626.0
288.0

全て上昇した。

送信してみる。今は順位126->134らしい。

００：３４
三回目のsubmit。57.66点ゲット。134->70位。100位以内に入ったのがなんとなく安心材料。まだまだ油断できないけど。

００：３９
これ以上、先に進める場合は、実行時間によって最善の選択肢が変化する可能性があるので、最適化を先に行っておくほうがいいかもしれない。

1:57
正多角形を作ろうとする試みを何回行っているかメモするようにした。
グローバル変数をinfoを作った。あまりグローバル変数を使いたくないが仕方ない。

２：３０
何気なく結果を見てたけど、どう考えても処理に時間がかかりすぎている。例えば最初の５秒で、7480のバケットしか
調べられていない。多分bucketのサイズが大きくなりすぎている。（ひとつのバケットに点が入りすぎている）のが原因と思われる。
２：５０
変数がよほど大きくならない限りは、最も近い点を調べる作戦に大きな問題があるとは思われない。

３：１６
ボロノイ図的なものをつくるのがいいかな。

３：２９
時間を増やしても、SCOREの上昇は多少に留まる。もう点数は限界が近いかな。

１２：３１
起きた。なんかぼろぼろだ。

１９：００
一応profileをとったらsq()に時間がかかってるので、やっぱりbucket大きすぎ説は正しい。

１９：０１
二度寝してしまったので、

２１：２２
ボロノイ図作ってみたけど、まだ速度が甘いかな。

２１：２３
と思って、何も変更せずにもう一度送信してみたら、探索量が増えてる。謎。さっきたまたま遅かっただけか・・・。

２３：３４
探索量があまりに不安定だし、example testでも似たような感じだから、なにかおかしなことが起きてそうだ。調べたほうがいい。
177.0
1294.0
11077.0
17272.0
1117.0
10746.0
1064.0
675.0
5225.0
288.0

6/13

ボロノイ的なことをしたら、余計に得点が下がってしまったが、これは多分、近距離候補の１０点をを全て使い切ってしまっているため。対策はどうしよう・・・。

０：３８
きれいなやり方ではないけど、残りのpointの数が減ったら、全探索に切り替えるという手もある。

３：４４
次の手はどう打つべきか・・・。ボロノイを利用する作戦はあまりうまく行ってないようにみてる。色々なサイズのbucketを用意するのがいい気がする。

５：３８
一時間以内に次の作戦を決める。

６：２７
bucket法を拡張するのが妥当という結論に落ち着いた。

６：２８
点の数が多い場合は思った以上に煮詰まっているんじゃないかという気がしてきた。

７：３７
何か適当に改良してexample testしてみた。
177.0
1452.0
13136.0
21312.0
1272.0
12207.0
1108.0
821.0
6169.0
299.0

ローカルよりいい点数が出やすいようだ。これはちょっとやりにくい・・・。

７：３９
submitしてみるか・・・。現在５４．２９点で８３位。これは４回目のsubmit。

7:43
結果は65.64点で５４位だった。土日を乗り切ったから多分、もう１００以内からは多分落ちない。
前日ラッシュの可能性も無いわけではないけど。

８：１１
しまった。よく考えたら -O2オプションつけるの忘れてた。そりゃローカルの挙動がおかしくなって当然だ。あほだ・・・。

９：４６
頂点数が多い場合は遅かれ早かれ、最適化ゲーになる可能性が高い。とりあえず、まずは頂点数<500の対策を先にやっておこう。

１７：２９
まず、最後の三角形の選択の改良をどうにかしないと。

6/14

22:22
もう残り時間が少ない。結局だらするだけで時間が過ぎ去ってしまっているので、今からは本気でやる。
今は７０位なので非常に危ない。

まず、方針はひとつに絞ってやるべき。あんまりだらだらやってると、予選落ちもありうる。
現在は正多角形を、選んだら、それで固定してその点は使わないアルゴリズムになっている。
長所は、使われていない点による正多角形のみを、検索しているので効率がよい。
短所は、一般に解は最適とは限らない。
この短所を克服すれば多少の点数アップが望めると思われる。
というより、これをいかに行うかが高得点化の鍵だと思われる。
そのためには正多角形をリストアップしてから、どの多角形を使うか決める方法がいい。
ただし、この戦略のせいで、長所が死んでしまうと、逆に点数が落ちかねないからそれは気をつける必要がある。
もうひとつ、今までは、時間が来たら探索を終了していた。
正多角形をリストアップしてから、どれを使うかきめる方法をとったとき、
正多角形をきめるステップの実行時間を正しく評価できていないと、タイムオーバーの危険があったりで、
有効に計算時間を使うのが難しい。正多角形を見つけるたびに動的に使う多角形リストを変化させれば、単に時間いっぱい計算すればいいことになって、有効に時間を使える。どっちの作戦がいいかは分からない。
安全なのは後者なので、後者の作戦でいくべきだろうか？残りの時間は少ないので。
後者の作戦で行くとして、どのようにアルゴリズムを設計するべきだろう？
まず、辺数の大きいほうから決定することに大きな問題はない気がする。大きな辺数の多角形を犠牲にして逆に点数を稼げるパターンも無いわけではないけど、割合としては少ないだろう。より上位を目指すならそうするべきだが、とりあえず後回し。
で、動的に多角形リストを更新ってどうやってやるんだ？

２２：４９
ちょっと紙に書いて考える。いつもこのまま考えてる間に堕落するパターンなので、２３：１０までに切り上げて、状況を整理する。２３：００
前回のR２と同じようにラグランジュ緩和法をつかうのはひとつの手だ。多分前回よりうまく行きやすい状況だとは思う。
でも、正直なんでうまくいってるのか自分でもあまり理解していない方法なので、あまり使いたくはない。というのがある。

２３：０３
しかし、これはどうみても重みつきset packing problemだよなあ。納得はいかないけど、ラグランジュ緩和法をつかうのは仕方ない気がしてきた。

２３：０９
packing problemとして定式化する方法を詰めてみる。一つ目の問題点は、使われていない点による正多角形の検索を放棄していること。

２３：１１
現在の正多角形のリストアップ方法は、
１：ランダムに二点を選ぶ
２：その二点を対点とする正多角形を検索する。
３：そのさい頂点は、理想の位置に最も近い点で、まだ使われていない点を検索する。
という方法をとっている。

２３：１８
最悪、前のアルゴリズムを繰り返しまくる作戦でもそこまで大きな問題はない気がする。

２３：２３
とりあえず、ラグランジュ緩和的なことをすることで、方針は決定した。よく考えればこれは非常にいい方法だという気がしてくる。

２３：２５
まだ一時間ほどしか考えてないが疲れた。とりあえず後一時間ほどは双対的な方法のやり方を自分なりに整理する。

6/15

５：３９
単体法を今から実装する。

単体法は放棄しました。

２２：３８
メモをとる余裕もないぐらいあせりまくってます。順位は１１０位。やばい。
example testの結果。

179.0
1559.0
13339.0
21514.0
1333.0
12434.0
1158.0
830.0
6400.0
328.0

一応全て点数が上がったから何とかなるかな。何とかなってほしい。

２２：４１
最適化に取り組む。

２３：５３
送信した。71.23点。５０位。

example 15のresult
179.0
1575.0
14283.0
24640.0
1262.0
12240.0
1169.0
853.0
7120.0
328.0

2011-05-08

Google Code Jam 2011 Qualification Round

全てpythonでGolfして解きました。本当はbefungeとかで解きたかったけど、予想以上に問題がHARDだったのでやめました。
努力不足なのであんまり短くないです。

A(204 Byte)

for T in range(input()):
 x=raw_input().split();t=i=1;w=[1,1];d=[0,0]
 while i<len(x):
	c=x[i]<"F";v=int(x[i+1])
	if abs(w[c]-v)<=d[c]:i+=2;w[c]=v;d[c]=-1
	d[0]+=1;d[1]+=1;t+=1
 print"Case #%s:"%(T+1),t-1

B(309 Byte)

for T in range(input()):
 x=raw_input().split();i=int(x[0]);m={};d={};h=[]
 for s in x[1:i+1]:m[s[:2]]=m[s[1::-1]]=s[2]
 for s in x[i+2:-2]:d[s]=d[s[::-1]]=1
 for c in x[-1]:
  while h and h[-1]+c in m:c=m[h.pop()+c]
  if any(a+c in d for a in h):h=[]
  else:h+=[c]
 print"Case #%s:"%(T+1),`h`.replace("'","")

C(147 Byte)

for T in range(input()):n=input();x=sorted(map(int,raw_input().split()));z=i=0;exec"z^=x[i];i+=1;"*n; print"Case #%s:"%(T+1),[sum(x[1:]),"NO"][0<z]

D(112 Byte)

for T in range(input()):input();print"Case #%s:"%(T+1),sum(i+1!=int(v)for i,v in enumerate(raw_input().split()))

2011-03-29

■

D,Cを解きました。

D

def readints():
  return map(int,raw_input().split())

def main():
  global n,m
  n,m=readints()
  data=[raw_input().split() for j in xrange(n)]
  alldata=[d for d2 in reps(data,"J1") for d in reps(d2,"J2")]
  #print alldata
  for d in alldata:
    allplace=[(i,j) for i in xrange(n) for j in xrange(m) if okp(i,j,d)]
    for x1,y1 in allplace:
      for x2,y2 in allplace:
        if abs(x1-x2)<3 and abs(y1-y2)<3:continue
        print"Solution exists."
        joker_mes=[]
        for jk in ["J1","J2"]:
          for i in xrange(n):
            for j in xrange(m):
              if data[i][j]==jk:
                joker_mes.append((jk,d[i][j]))
        if len(joker_mes)==0:
          print"There are no jokers."
        if len(joker_mes)==1:
          jk,c=joker_mes[0]
          print "Replace %s with %s."%(jk,c)
        if len(joker_mes)==2:
          (j1,c1),(j2,c2)=joker_mes
          print "Replace J1 with %s and J2 with %s."%(c1,c2)
        mes((x1,y1),0)
        mes((x2,y2),1)
        return
  print "No solution."
  pass

def ins(x,y):
  return 0<x<n-1 and 0<y<m-1

def reps(data,s):
  rem = allcards-set(data[i][j]for i in xrange(n) for j in xrange(m))
  ret=[]
  for i in xrange(n):
    for j in xrange(m):
      if data[i][j]==s:
        for card in rem:
          ndata=deepcopy(data)
          ndata[i][j]=card
          ret.append(ndata)
        return ret
  return [data]

def mes(p,n):
  th=["first","second"][n]
  print"Put the %s square to (%s, %s)."%(th,p[0],p[1],)

def okp(x,y,dd):
  if not ins(x,y):return False
  adj=[(x+dx,y+dy) for dx in [-1,0,1] for dy in [-1,0,1]]
  return all(dd[x][y][1] == dd[nx][ny][1] for nx,ny in adj) or len(set(dd[nx][ny][0] for nx,ny in adj)) ==9


from copy import deepcopy
ranks="2 3 4 5 6 7 8 9 T J Q K A".split()
suits="C D H S".split()
allcards=set([ra+su for ra in ranks for su in suits])

if __name__=='__main__':
  main()

C

def readints():
  return map(int,raw_input().split())

def main():
  n=input()
  x=readints()
  for k in xrange(3,n+1):
    if n%k:continue
    m=n/k
    for i in xrange(m):
      if all(x[i+t*m] for t in xrange(k)):
        print "YES"
        return
  print "NO"
  pass