数値計算アルゴリズム

ユークリッドの互除法(Euclidean Algorithm)#

ユークリッドの互除法は，2つの自然数の最大公約数を求めるための古典的なアルゴリズムである． $a,b\in \Z,a>b$ とした場合， $\gcd(a,b)=\gcd(b, a \mod b)$ の関係を用いて下記コードのように求めることができる．

1
int _gcd(int a, int b)
2
{
3
    return b == 0 ? a : _gcd(b, a % b);
4
}

$\gcd(a, b)=\gcd(b,a)$ を考慮した実装( $a, b$ の大小比較)は[gcd.cpp]を参照されたい．

最小公倍数(Least common multiple)#

正整数 $a,b$ に対して，最大公約数 $\gcd(a,b)$ と最小公倍数 $\text{lcm}(a,b)$ の間には，

\gcd(a,b)\cdot\text{lcm}(a,b)=ab

の関係があるので，下記コードにて求めることができる．完全なコードは，[lcm.cpp]を参照されたい．

1
int lcm(int a, int b)
2
{
3
    return a * b / gcd(a, b);
4
}

ニュートン法(Newton’s method)#

ニュートン法は数値計算による方程式解法で，微分可能性と2次微分符号を条件に，反復法で求根を行う．問題を定式化すると，関数 $f:\R\to\R,x\in\R$ として， $f(x)=0$ となる $x$ を求めること限定する．このとき， $x$ の付近に適当な値 $x_0$ をとり，次の漸化式によって， $x$ に収束する数列を得ることができる場合が多い．

x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}

具体的には，以下の反復を行う．完全なコードは，[newton.cpp]を参照されたい．

1
double x0 = 1.0;
2
double tolerance = 1e-6;
3
int max_iter = 1000;
4
int iter = 0;
5

6
while (iter < max_iter)
7
{
8
    double x1 = x0 - f(x0) / df(x0);
9

10
    if (abs(x1 - x0) < tolerance)
11
        break;
12

13
    x0 = x1;
14
    iter++;
15
}

ガウスの消去法(Gaussian elimination)#

ガウスの消去法は，線形方程式系 $Ax = b$ を解くために係数行列 $A$ を前進消去（forward elimination）と後退代入（back substitution）によって上三角行列へ変換し，逐次的に解を求める直接法である．前進消去では，行基本変形により，対角成分を基準に下側の成分をゼロにする．数値的安定性を考慮する場合，部分ピボット選択（partial pivoting）を併用する．後退代入により，最下行から順に未知数を逐次決定し，最終的に解ベクトル $x$ を得る．完全なコードは，[gauss.cpp]を参照されたい．

1
vector<double> gauss(vector<vector<double>> A, vector<double> b)
2
{
3
    int n = A.size();
4
    for (int i = 0; i < n; i++)
5
    {
6
        int p = i;
7
        for (int j = i + 1; j < n; j++)
8
            if (abs(A[j][i]) > abs(A[p][i]))
9
                p = j;
10
        swap(A[i], A[p]), swap(b[i], b[p]);
11
        for (int j = i + 1; j < n; j++)
12
        {
13
            double f = A[j][i] / A[i][i];
14
            for (int k = i; k < n; k++)
15
                A[j][k] -= f * A[i][k];
16
            b[j] -= f * b[i];
17
        }
18
    }
19
    vector<double> x(n);
20
    for (int i = n - 1; i >= 0; i--)
21
    {
22
        x[i] = b[i];
23
        for (int j = i + 1; j < n; j++)
24
            x[i] -= A[i][j] * x[j];
25
        x[i] /= A[i][i];
26
    }
27
    return x;
28
}

素数判定(Primality test)#

素数判定とは，与えられた自然数が素数か合成数かを判定することである．素数判定を行うアルゴリズムを素数判定法という．完全なコードは，[check_prime.cpp]を参照されたい．

1
bool is_prime(unsigned int n)
2
{
3
    if (n < 2)
4
        return false;
5
    if (n == 2)
6
        return true;
7
    if (n % 2 == 0)
8
        return false;
9
    for (unsigned int i = 3; i * i <= n; i += 2)
10
    {
11
        if (n % i == 0)
12
            return false;
13
    }
14
    return true;
15
}

sql-hkr

algo

Waiting for api.github.com...

00K

Waiting...