■[俺用メモ]ニュートン力学

アルツってつらいわ。

力学

1. 位置、速度、加速度
   • 
   • 
   • 
   • 
2. 運動法則
   • （力の釣り合い）
   • （慣性）
   • （作用反作用の法則）
   • （天秤の釣り合い）
   • （重心位置）
   • （重心速度）
   • （重心の等速運動）
3. 重力
   • （万有引力の法則）
   • （重力場）
   • （重力質量と慣性質量の等価性）
   • （落下）
   • 
   • （位置エネルギー）
4. 等速円運動
   • 
   • 
   • 
   • （周期）
   • （バネの力）
   • （バネのエネルギー）
   • （単振動の角速度）
   • （単振り子の角速度）
5. 力学的保存量
   • （力積）
   • （運動量保存）
   • （運動エネルギー）
   • （力学エネルギー保存）
   • （仕事）
   • （仕事率、ワット）

波動

1. 波
   • 横波: 振動方向と進行方向が直交
   • 縦波: 振動方向と進行方向が同じ
   • （周波数）
   • （波長）
   • 
   • （干渉の振幅最大条件）
   • （干渉の振幅最小条件）
2. 波の屈折と反射
   • （媒体境界での周波数の連続）
   • （媒体境界での波長の変化）
   • （スネルの法則）
   • （相対的な屈折率）
   • 固定端反射: 端で常に振幅0、反射波は位相逆転
   • 開放端反射: 端で常に振幅最大、反射波は位相不変
3. 音
   • （波長のドップラー効果）
   • （周波数のドップラー効果）
   • （うなり振動数）
   • （弦、両閉、両開気柱の振動モード）
   • （片開気柱の振動モード）
   • （弦の周波数）
4. 光
   • 屈折率は波長によって変化する
   • 横波のため偏光が存在する
   • （絶対屈折率）
   • （全反射の条件）
   • （強めあう回折条件）
   • （弱めあう回折条件）
   • （真空中の光速）

Permalink | コメント(0) | トラックバック(0)

■[俺用メモ]相対性理論

天才アインシュタインが、当時、誰も思いつかなかった高度な理論を唐突に発表した。
当時は、誰もが彼の理論を疑ったが、その後、彼が正しいことが証明された。

Permalink | コメント(0) | トラックバック(0)

■[俺用メモ]C/C++ スタイルルール

とりあえずのスタイルルール

基本的にはその言語の開発者の流儀に従う。

C は Brian W. Kernighan and Dennis M. Ritchie, The C Programming Language, second edition (Prentice-Hall, 1988) のスタイル（K&R 流）に従う。

C++ も基本的には上と同じでいいが， Bjarne Stroustrup, The C++ Programming Language, third edition (Addison-Wesley, 1997) に従って char *p; ではなく char* p; と書く。 私の CGI で学ぶ C++ のコード例参照。

Java も基本的には上と同じでいいが， Java 独自の部分は Sun から出ている Java の本の流儀に従う。 私の Java 言語について のコード例参照。

インデントは K&R の本にならって4桁ごととする。 TAB コードは使わず，スペースで展開する^*。

^* 『C言語による最新アルゴリズム事典』 のソースコードでは VZ エディタで C を編集したときの既定値で TAB を4桁としたが， 環境によって見え方が違ってしまった。 ちなみに，この本は紙を節約するため詰め込みすぎたところがあるが， これは真似る必要はない。

右側のコメントは33バイト目から始める（これはあまり気にする必要はない）。 33という数字は下の indent や CC mode での既定値。

//345678901234567890123456789012345678901234567890
#include <iostream>             // 入出力の準備
using namespace std;            // std を使う
int main(int argc, char* argv[]) // メイン
{
    for (int i = 1; i < argc; i++) // 引数を出力
        cout << argv[i] << '\n';
}
</iostream>

indent を使う

GNU の indent を次のように使えば C ソースを（ほぼ）上のルールに則って書き直してくれます：

indent -kr -ts255 *.c

-kr は K&R 流，-ts255 は255文字のスペースをTABコードに置き換える （したがって実際にはTABの代わりにスペースを使う） という意味です。

たとえば次のようなコードを処理してみましょう。

int
foo( double x ){
    if( sin( x ) == 0 )
    {
        return 3;
    }
    else
    {
        return 5;
    }
}

indent -kr -ts255 の出力は次のようになります。

int foo(double x)
{
    if (sin(x) == 0) {
        return 3;
    } else {
        return 5;
    }
}

indent は C++ でも簡単な場合は使えますが，完全ではありません。 Java ではさらに悲惨なことになります。

CC mode を使う

XEmacs/Emacs/Mule ではファイルを開くと拡張子で自動的に言語を判断し， 適当な編集モードになります。 インデントのルールを指定しておけば， 自動的にそれに則ってインデントしてくれますので， スタイルをそろえるのに便利です。 最下行に C とか C++ とか Java とか出るはずですので，確認してください。 実際に働いているのは CC mode という emacs lisp プログラムです。

インデントのルールを K&R 流にするには， ホームディレクトリの .emacs に次のように書き加えます：

(require 'cc-mode)

;; Kernighan & Ritchie スタイルにする
(setq c-default-style "k&r")

;; BackSpace キーを「賢く」し，インデント幅は4桁，タブはスペースに展開
(add-hook 'c-mode-common-hook
     	  '(lambda ()
             (progn
               (c-toggle-hungry-state 1)
               (setq c-basic-offset 4 indent-tabs-mode nil))))

;; .hpp と .h を C++ の拡張子とする（.h の方は余計かも）
(setq auto-mode-alist
      (append
       '(("\\.hpp$" . c++-mode)
         ("\\.h$"   . c++-mode)
         ) auto-mode-alist))

これでも望みのモードにならない場合は，たとえば M-x c++-mode のように打ち込みます。 あるいは，ファイルの先頭に -*- C++ -*- を含む文字列を，たとえば

// This is a -*- C++ -*- header.

CC mode のよく使う機能をまとめておきます。

• TAB でインデントしてくれる
• C-j で改行してインデントしてくれる
• M-; でコメントが挿入できる

（付）CC mode のインストール

最新の CC mode のインストールは次のようにします。

まず CC Mode サイトから最新版をいただいてきて，Emacs 用のディレクトリ（下の例では /usr/local/share/emacs/site-lisp）で展開します。

cd /usr/local/share/emacs/site-lisp
tar xvzf .../cc-mode-5.30.8.tar.gz
ln -s cc-mode-5.30.8 cc-mode

Mule (Emacs 19.34) の場合は cc-make.el の

(defconst cc-path-to-the-custom-library
  ;; 中略
  nil
  )

(defconst cc-path-to-the-custom-library
  ;; 中略
  "/usr/local/share/emacs/site-lisp/cc-mode/"
  )

EMACS   = emacs

EMACS   = mule

そうしてから，

mule -batch -no-site-file -q -l cc-make.el cc-*.el

.emacs に上記のパスを設定しておきます。

(setq load-path (cons "/usr/local/share/emacs/site-lisp/cc-mode" load-path))

Permalink | コメント(0) | トラックバック(0)

■[俺用メモ]アルゴリズム

y² < x < (y + 1)² … (1)

x == (y + 1/2)²

x == y² + y

err(+1/4..+1/2 LSB) = 536872070
err(-1/4..+1/4 LSB) = 1073739508
err(-1/2..-1/4 LSB) = 536872070

g(x) = (22 * x + 11) / 32 … (2)

Fig.1 sqrt(x) の近似誤差 (初期値とニュートン法１回目)

/* NAME * sqrtl - square root function * SYNOPSYS * long * sqrtl(long x) * DISCRIPTIONS * The sqrtl() function compute the non-negative square root of x. * ERROR * Below 1/2 LSB. * SEE ALSO * sqrt(3), http://www.finetune.co.jp/~lyuka/fract/sqrt_hypot.html * COPYRIGHT * Copyright 2002, Takayuki HOSODA. All rights reserved. */ long _sqrtl(a) long a; { register long x; register long t; register long s; int scale; x = a; if (x > 0) { scale = 0; if (x < 0x8000) { x <<= 16; scale = 8; a = x; } x >>= 8; s = 8; for (t = 0x400000L; x < t; t >>= 2) s--; t = 88; t <<= s; x *= 22; s += 5; x >>= s; /* -3.1e-2 < err < +2.9e-2 */ s = a; t += x; x = s; s /= t; s += t; s >>= 1; /* -4.8e-4 < err <= 0 */ t = x; x /= s; x += s; x >>= 1; /* -1.2e-7 < err <= 0 */ s = x; s++; s *= x; if (t > s) /* adjust LSB */ x++; if (scale) { x += 127; x >>= 8; } } return x; }

g(x) = (53 * x + 129) / 128 … (3)

h(x) = 13 * x / 32 + 1 … (4)

Fig.2 sqrt(1+x) の近似誤差 (初期値とニュートン法１回目)

/* NAME * hypots - euclidean distance function * SYNOPSYS * unsigned short * hypots(short x, short y) * DESCRIPTION * The hypots() function compute the sqrt(x * x + y * y). * ERROR * Below 1/2 LSB * SEE ALSO * hypot(3), http://www.finetune.co.jp/~lyuka/fract/hypot.html * COPYRIGHT * Copyright 2002, Takayuki HOSODA. All rights reserved. */ unsigned short _hypots(x, y) short x; short y; { register unsigned long a; register unsigned long b; register unsigned long c; a = x; b = y; if(x < 0) a = -x; if(y < 0) b = -y; c = b; if(a != 0) { c = a; if (b != 0) { if(a < b) { a = b; b = c; } b *= b; c = b; b /= a; b *= 13; b >>= 5; b += a; /* 0 < error < +1.9e-2 */ a *= a; c += a; a = c; /* apply newton-raphson methode twice */ a /= b; a += b; a >>= 1; /* -1.7e-4 < error < 0 */ b = c; c /= a; c += a; c >>= 1; /* -1.3e-8 < error < 0 */ a = c; a++; a *= c; if (b > a) /* adjust LSB */ c++; } } return (unsigned short)c; }

/* NAME * hypot14 - euclidean distance function * SYNOPSYS * short * hypot14(short x, short y) * DESCRIPTION * The hypot14() functions compute the sqrt(x * x + y * y). * ERROR * Below 1/2 LSB * BUGS * Input x, y must be within range of -11,715 to 11,715 for 1/2 LSB error. * SEE ALSO * hypot(3), http://www.finetune.co.jp/~lyuka/fract/hypot.html * COPYRIGHT * Copyright 2002, Takayuki HOSODA. All rights reserved. */ #include <math.h> short hypot14(short x, short y) { short x; short y; { register long a; register long b; register long c; a = (long)abs(x); b = (long)abs(y); c = b; if (a != 0) { c = a; if (b != 0) { if(a < b) { a = b; b = c; } b *= b; c = b; b /= a; b *= 53; b += a; b >>= 7; b += a; /* |err| < 8.5e-3 */ a *= a; c += a; /* apply newton-raphson methode */ a = c; c /= b; c += b; c >>= 1; /* -7.2e-5 < error < 0, 1/sqrt(2^27) = 8.6e-5 */ b = c; b++; b *= c; if (a > b) /* adjust LSB */ c++; } } return c; }

g(x) = (39 * x + 29) / 64 … (5)

x_n+1 = x_n * (x_n³ + 2 a) / (2 x_n³ + a) … (6)

Fig.3 cbrt(x) の近似誤差

/* NAME * cbrtl - cube root function. * SYNOPSYS * long * cbrtl(long x) * DISCRIPTIONS * The cbrtl() function compute the cube root of x. * ERROR * Below 1/2 LSB. * SEE ALSO * cbrt(3), http://www.finetune.co.jp/~lyuka/fract/sqrt_hypot.html * COPYRIGHT * Copyright 2003, Takayuki HOSODA. All rights reserved. */ long cbrtl(a) long a; { long x, t, b; if (a == 0) return 0; if (a == 0x80000000L) return -1290; x = b = abs(a); for (t = 1; x > t; t <<= 1) x >>= 2; x = ((x * 39 + t * 29) >> 6) + 1; /* -0.06 < err < 0.06 */ t = b / (x * x); x = (x * (t + t + x)) / (x + x + t); /* |err| < 1.5e-4 */ t = x * x; if ((6 * t + 3 * x) < (b - x * t) * 4) x++; /* adjust LSB */ if (a < 0) return -x; return x; }

x' = x * ((n - 1) * xn + (n + 1) * a) / ((n + 1) * xn + (n - 1) * a)

Permalink | コメント(0) | トラックバック(0)

■[俺用メモ]ニュートン法

平方根

X_n+1 = (X_n + A / X_n) / 2

X_n+1 = X_n * (3 - A * X_n²) / 2

h_n = 1 - A * x_n²
x_n+1 = x_n * (1 + h_n * (4 + 3 * h_n) / 8)

h_n は、x_n の誤差と同じである。繰り返しごとに有効桁は３倍になる。

h_n = 1 - A * x_n²
x_n+1 = x_n * (1 + h_n * (8 + h_n * (6 + 5 * h_n)) / 16)

h_n = 1 - A * x_n²
x_n+1 = x_n * (1 + (64 * h_n + h_n² * (48 + 40 * h_n + 35 * h_n²)) / 128)

h_n = 1 - A * x_n²
x_n+1 = x_n * (1 + h_n * (1/2 + h_n * (3/8 + h_n * (5/16 + h_n * (35/128 + h_n * 63/256)))))

h_n = 1 - A / x_n²
x_n+1 = x_n * (1 - h_n * (1/2 + h_n * (1/8 + h_n * (1/16 + h_n * (5/128 + h_n * 7/256)))))

h_n = 1 - A * x_n³
x_n+1 = x_n + x_n * h_n * (1/3 + h_n * (2/9 + h_n * (14/81 + h_n * (35/243 + h_n * 91/729))))

h_n = 1 - A * x_n⁴
x_n+1 = x_n + x_n * h_n * (1/4 + h_n * (5/32 + h_n * (15/128 + h_n * (195/2048 + h_n * 663/8192))))

#include <stdio.h>
#include <math.h>
#ifndef DBL_EPSILON
#define DBL_EPSILON     2.2204460492503131E-16
#endif

double
_sqrtinv(a)
    double a;
{
    double x, h, g;
    int    e;

    /* enhalf the exponent for a half digit initial accuracy */
    frexp(a, &e);
    x = ldexp(1.0, -e >> 1);

    /* 1/sqrt(a) 4th order convergence */
    g = 1.0;
    while(fabs( h = 1.0 - a * x * x) < fabs(g)) {
        x += x * (h * (8.0 + h * (6.0 + 5.0 * h)) / 16.0);
        g = h;
    }
    return(x);
}

double
sqrt(a)
    double(a);
{
    if(a < 0){
        fprintf(stderr, "sqrt: domain error\n");
        return(0);
    }
    return(a * _sqrtinv(a));
}

Permalink | コメント(0) | トラックバック(0) | 18:55

■[俺用メモ]収束公式-GaussLegendre-

very_long_float
pi2(iter)
    int             iter;
{
    int             i;
    very_long_float a, b, t, x, y;

    a = 1;
    b = 1 / sqrt(2);
    t = 1;
    x = 4;
    for (i = 0; i < iter; i++) {
        y = a;
        a = (a + b) / 2;
        b = sqrt(y * b);
        t -= x * sqr(y - a);
        x *= 2;
    }
    return (sqr(a + b) / t);
}

(但し、sqr(x) は x², pow4(x) は x⁴, sqrt(x) は x^1/2, qroot(x) は x^1/4)

err( 1) = -3.2257622e-4
err( 2) = -2.3479336e-9
err( 3) = -5.8292283e-20
err( 4) = -1.7418264e-41
err( 5) = -7.6589246e-85
err( 6) = -7.3484406e-172
err( 7) = -3.3697129e-346
err( 8) = -3.5362794e-695
err( 9) = -1.9454428e-1393
err(10) = -2.9425892e-2790

very_long_float
pi(iter)
    int             iter;
{
    int             i;
    very_long_float a, b, y, t;

    y = sqrt(2) - 1;
    a = 6 - sqrt(32);
    b = 8;
    for (i = 0; i < iter; i++) {
        t = qroot(1 - pow4(y));
        y = (1 - t) / (1 + t);
        t = sqr(1 + y);
        a = sqr(t) * a - y * (t - y) * b;
        b *= 4;
    }
    return (1 / a);
}

err(1) = -2.3479336e-9
err(2) = -1.7418264e-41
err(3) = -7.3484406e-172
err(4) = -3.5362794e-695
err(5) = -2.9425892e-2790

very_long_float h;
very_long_float a0 = 355.0/113.0;


/* ai+1 = ai + sin(ai) は π に収束する。
 * ai+1 = ai - tan(ai) でも同じ。
 */ 

while(abs(h = sin(a)) > epsilon) {
    a += h;
}

err(1) = 1.007e-21
err(2) = 1.680e-63
err(3) = 7.804e-189

Permalink | コメント(0) | トラックバック(0)