2011-08-24

Project Euler 87

python

問題は、Problem 87 - Project Euler。
50 * 10^6以下の数の中で、素数の二乗+三乗+四乗でかけるものの数を求めよというものです。

解答している人が前後の問題と比べて二倍程多いだけあって単純な総当たりでも数秒程で答えがでてしまいます。

from number_theory import *

p1 = primes(int(pow(50 * 10**6, 0.5)))
p2 = primes(int(pow(50 * 10**6, 1.0/3.0)))
p3 = primes(int(pow(50 * 10**6, 1.0/4.0)))

candidates = []

for x in p1:
    for y in p2:
        for z in p3:
            candidates += [x**2 + y**3 + z**4]

print len([x for x in set(candidates) if x < 50 * 1000000])

2011-08-15

レビ・チビタの記号についての2,3の公式

数学 memo

3次元のレビ・チビタの記号についての公式をまとめておきます。

以下縮約ルールを用いる。(同じ添字が現れたら1〜３までの和を取る.)

同じ添字が1つ現れる場合

$\epsilon_{iab}\epsilon_{inm}=\delta_{an}\delta_{bm}-\delta_{am}\delta_{nb}$

同じ添字が2つ現れる場合

$\epsilon_{ijb}\epsilon_{ijm}=\delta_{jj}\delta_{bm}-\delta_{jm}\delta_{jb}$
$= 3\delta_{bm}-\delta_{mb}$
$= 2\delta_{bm}$

すべて同じ添字の場合

$\epsilon_{ijk}\epsilon_{ijk}=2\delta_{kk}=2\times 3 = 3!$

追記

d次元の場合にも類似した公式を導くことが出来る。
例えば、

すべて同じ添字の場合

$\epsilon_{i_1i_2\ldots i_d}\epsilon_{i_1i_2\ldots i_d}=d!$

2011-08-13

ピタゴラス三体問題

物理

ピタゴラス3体問題 (Burrau’s problem of three bodies) | YABUKI Taro’s Home Page
のを見て興味を持ったのでグラフを書いて見ました。

重力多体系の数値計算を行なう場合、力が距離の逆二乗に比例する為粒子間の距離が近いときには誤差がとても大きくなってしまいます。
これを防ぐには、ポテンシャルの距離の部分に適当な定数を加えてポテンシャルを丸めたり、距離が近いところでは時間刻みを細かくする等の工夫が必要になります。

今回は距離が近いところでの時間刻みを細かくすることによって誤差の増大を抑えることにしました。

粒子の軌跡

粒子2つが組になり残り1つは弾き飛ばされてしまいます。

エネルギーと運動量は次の様になりました。

多体系の全エネルギー

1%程度の誤差が出ています。

多体系のx方向の全運動量

エネルギーに比べて運動量は比較的良く保存しています。なんでだろ?

2011-08-12

理想気体の断熱関係式を統計力学から導く

物理統計力学

熱力学の教科書を見ると大抵導出が乗っていますがここでは、統計力学を使ってそれを求めて見ます。

理想気体のエントロピー

カノニカル分布関数は、
$Z_N = \frac{1}{h^3}\int d^3pd^3q \exp \left[ -\beta \frac{p^2}{2m} \right$ ]
$= \frac{V^N}{h^{3N}}\left(\frac{2\pi m}{\beta} \right)^{3N/2}$
これを用いてエントロピーは次の様に求めることができる。

$S = -\frac{\partial}{\partial T}\left( -kT \log \left( \frac{V}{h^3} \left( 2\pi m kT \right)^{3/2} \right) \right)$
$= N\left( k\log \left( V \left( \frac{2\pi m k T}{h^2} \right)^{3/2} \right)-\frac{3}{2}k\right)$

ここで N は定数とする。

エントロピーの変化が0になることを用いて断熱関係式を導出する。

$dS \propto \frac{1}{V} dV + \frac{3}{2}\frac{1}{T}dT$

であるが特に断熱変化ではこれが $0$ にならなければならならないので、
$V T^{3/2} = const$
が成立する。

はてなのtex記法ちょっと使い難いです。

2011-08-11

cで木構造を書く2

memo c++ tree

http://serennz.sakura.ne.jp/sb/log/eid120.html
を読んだまとめ２。

前回作った木構造に作用する関数を書いてみます。
具体的には、
1. 自分の階層とその下の階層そのさらに下の階層と再帰的に作用する関数。
2. 自分の階層とその上の階層そのさらに上の階層と再帰的に作用する関数。
を作ります。

ここでは各nodeに作用する関数を、hogeとしましょう。
実装には再帰を用います。

まず1.の場合は

void re_hoge1 (node *n) {
  if(n == NULL) return;
  hoge(&n);
  re_hoge1(n->next); 
  re_hoge1(n->child);
  return;
}

この場合の終了条件は二行目

if(n==NULL) return;

で、引数がNULLになることです。
３行目４行目では自分自身を呼び出しています。

 re_hoge1(n->next); 
 re_hoge1(n->child);

前回例で出した構造の場合にどのように動くのかを考えてみます。
re_hoge1(&root);
とすると、これはNULLではないので、一行目はスルー。
二行目でrootにhogeが作用します。
３行目のre_hoge1には,root.next=NULLなのでNULLが入ります。
しかしここで呼び出された関数は、一行目の条件を満たすので、何もしないで終了します。
4行目のre_hoge1には,root.child = n1が入ります。
以下n1.nextからn2,n1.childからn3に作用してと再帰的にhogeが作用して行きます。

2.の場合は

void re_hoge2 (node *n) {
  if(n == NULL) return;
  hoge(%n);
  re_hoge2(n->next);
  re_hoge2(n->parent);
  return;
}

とします。

2011-08-08

cで木構造を書く1

c++ memo tree

http://serennz.sakura.ne.jp/sb/log/eid120.html
このページを参考に。
自分用に調べたことをまとめておきます。
院試が終わって時間が空けば自分で実装してみる予定。

実装にはポインターを用います。

struct node {
  struct node * parent;
  struct node * child;
  struct node * next;
  char *name;
  int depth;
};

nodeは自身と同じ型への３つのポインターparent,child,nextを持つ構造体として定義されます。
ここで,parentは今のnodeの一つ上のnodeを,nextはひとつ隣の,childはひとつ下のnodeを表します。
少し例を見てみましょう。

例えばroot以下、n1,n2,n1以下n3,n4が続くような構造を表したいのであれば、

- root
  - n1
    - n3
  - n2

まずrootを次の用に定義します。

node root, n1, n2, n3, n4;
root.parent = NULL;
root.child = &n1;
root.next = NULL;

rootは木構造の頂点に位置するのですからこれより上(parent)はありません(NULL)。
また今の例では頂点は一つですから同じ階層にこれに隣接するnodeもありません。(root.next = NULL)
rootの下にはn1が位置します(root.child = &n1）

次に,n1を定義します。

n1.parent = &root;
n1.child = &n3;
n1.next = &n2;

n1の上にはrootがあります。(n1.parent = &root)
n1の下にはn3が、同じ階層の隣にはn2があります。(n1.child = &n3, n1.next = &n2)

最後にn3,n2を定義しましょう。

n2.parent = &root;
n2.child = NULL;
n2.next = NULL;

n3.parent = &n1;
n3.child = NULL;
n3.next = NULL;

これで図の木構造を構成できたことになります。

2011-08-04

ランダムに配置された粒子を格子で区切る

c++ tree

粒子をランダムに配置した後、空間を粒子が一つか0こしか含まないような格子に区切るプログラムを書きました。
クイックソートの要領で再帰を使って実装しました。

実はこのプログラムは別の数値計算に用いる予定で書いたのですが予想外に時間が掛かりそうなので単体で記録として残しておくことにしました。

#include <iostream>
#include <vector>
#include <cmath>
using namespace std;
const int N = 1000;
const int x = 0, y = 1, z = 2;
typedef vector<double> Vec;
typedef vector<Vec> Mat;
typedef vector<Mat> Ten;
Mat X(N, Vec(3, 0));
Mat Dv;
void init() {
  for(int i = 0; i < N; i++) {
    for(int d = 0; d < 3; d++) {
      X[i][d] = drand48() * 20 - 10;
    }
  }
  return;
}
int how_many(double xmin, double xmax,
	     double ymin, double ymax,
	     double zmin, double zmax) {
  int count = 0;
  for(int i = 0; i < N; i++) {
    if((xmin < X[i][x] and X[i][x] < xmax) and
       (ymin < X[i][y] and X[i][y] < ymax) and
       (zmin < X[i][z] and X[i][z] < zmax)) {
      count ++;
    }
  }
  return count;
}
      
void dvide(double xmin, double xmax,
	   double ymin, double ymax,
	   double zmin, double zmax) {
  Vec tmp;
  tmp.push_back(xmin);
  tmp.push_back(xmax);
  tmp.push_back(ymin);
  tmp.push_back(ymax);
  tmp.push_back(zmin);
  tmp.push_back(zmax);
  Dv.push_back(tmp);
  if( 2 > how_many(xmin, xmax,
		   ymin, ymax,
		   zmin, zmax)) {
    return;
  } else {
    double dx = (xmax - xmin) / 2.0;
    double dy = (ymax - ymin) / 2.0;
    double dz = (zmax - zmin) / 2.0;
    dvide(xmin, xmin + dx,
	  ymin, ymin + dy,
	  zmin, zmin + dz);
    dvide(xmin, xmin + dx,
	  ymin, ymin + dy,
	  zmin + dx, xmax);
    dvide(xmin, xmin + dx,
	  ymin + dy, ymax,
	  zmin, zmin + dz);
    dvide(xmin, xmin + dx,
	  ymin + dy, ymax,
	  zmin + dz, zmax);
    dvide(xmin + dx, xmax,
	  ymin, ymin + dy,
	  zmin, zmin + dz);
    dvide(xmin + dx, xmax,
	  ymin, ymin + dy,
	  zmin + dz, zmax);
    dvide(xmin + dx, xmax,
	  ymin + dy, ymax,
	  zmin, zmin + dz);
    dvide(xmin + dx, xmax,
	  ymin + dy, ymax,
	  zmin + dz, zmax);
  }
}
void split() {
  Dv.clear();
  dvide(-10, 10,
	-10, 10,
	-10, 10);
}

int main(int argc, char **argv) {
  init();
  split();
  cout << Dv.size() << endl;
  return 0;
}