としおの読書生活

田舎に住む社会人の読書記録を綴ります。 主に小説や新書の内容紹介と感想を書きます。 読書の他にもワイン、紅茶、パソコン関係などの趣味を詰め込んだブログにしたいです。

カテゴリ: C++

C++で微分を実装してみたので忘備録として記録しておきます。


数値微分を実装する


f(x)=f(x+Δx)f(x)Δx

まずは典型的な上記の式の前進差分近似を使った方法で微分をしてみます。

#include <iostream>

double func(double x){
    return x*x;
}

double numericalDiff(double (*f)(double), double xdouble eps=0.001){
    return (f(x + eps) - f(x)) / eps;
}

int main(void){
    double (*pfunc)(double);
    pfunc = func;
    printf("%lf"numericalDiff(pfunc2));
    return 0;
}

4.001000

多少は誤差があるものの微分ができました。

関数ポインタを使ってnumericalDiffに微分したい関数を渡しているため、他の式に対しても微分を行うことが可能です。

また、引数epsの値を小さくすることでより精度の高い微分を多なうことができます。



誤差の少ない数値微分


前進差分近似を使うよりも中心差分近似を使った微分の方が誤差が少なくなるみたいですのでそちらも実装していきます。

先ほどのコードからnumericalDiff関数だけ以下の式に合わせて変更していきます。

f(x)=f(x+Δx)f(x - Δx)x


double numericalDiff(double (*f)(double), double xdouble eps=0.001){
    return (f(x + eps) - f(x - eps)) / (2*eps);
}

4.000000

中心差分近似を使うことで上記のサンプルに対しては誤差なく微分を行うことができました。



3753817_s

C++でモンテカルロ法を使って円周率を求めるプログラムを作ったので概要と実装を忘備録として残します。


モンテカルロ法とは


モンテカルロ法とは乱数を用いてシミュレーションや数値計算を行う方法の一つです。

少し前には将棋や囲碁のAIでモンテカルロ法を使うのがブームになりました。

今回はモンテカルロ法を使った例題としてポピュラーな円周率の近似値を求めるという問題を解いていきます。



円周率を求める手順


①正方形内にランダムに点を打つ。今回は正方形の大きさを1×1として(x, y)座標のx, yを0~1の範囲で乱数を使って点を打っていきます。

②生成した点が円の内側にあるか、外側にあるか判定する。内側にある場合はポイントを加算する。

キャプチャ

点が上記の図の青色の範囲内ならポイントを加算するということです。

③上記の①と②の行為をn回繰り返す。試行回数が多いほど円周率が正しい値へと近づいていきます。

④円周率の近似値をもとめる。円周率の近似値は以下の式で求めることができます。

π = 4 × ポイント






実装


#include <iostream>
#include <ctime>
#include <cstdlib>

int main(void){
    int n = 0;
    std::cout << "試行回数を入力してください" << std::endl;
    std::cin >> n;

    // 試行回数が0以下ならエラー
    if (n <= 0){
        return -1;
    }

    // 試行回数分ランダムに点を描画
    std::srand(time(NULL));
    int point = 0;
    double x = 0;
    double y = 0;
    for(int i=0i<n; ++i){
        // ①x, yは0~1の範囲
        x = (double)rand() / 32767;
        y = (double)rand() / 32767;
        // ②
        if (x*x+y*y < 1.0){
            // 円内に点が打たれた場合
            point++;
        }
    }

    // ④円周率を求める
    float pi = 4.0 * point / n;
    std::cout << "pi = " << pi << std::endl;
    
    return 0;
}




実行結果


1000
pi = 3.064

1000000
pi = 3.14275

試行回数が多いほど近似した値になっていることが分かりますね。





C++でStackを使用して逆ポーランド記法で入力された式の結果を求めるというAIZU ONLINE JUDGEのALDS_1_3_Aを解きました。



逆ポーランド記法とは


逆ポーランド記法とは、演算子をオペランドの後に記述するプログラムを記述する記法です。

例えば、(1 + 2) * (5 + 4)という中間記法で書かれた数式があったとする。

これを逆ポーランド記法では次のように書くことができます。

1 2 + 5 4 + *



問題


逆ポーランド記法で与えられた数式の計算結果を出力してください。

入力例


1つの数式が1行に与えられます。連続するシンボル(オペランドあるいは演算子)は1つの空白で区切られて与えられます。
1 2 +


出力例


計算結果を1行に出力してください。

3

https://onlinejudge.u-aizu.ac.jp/courses/lesson/1/ALDS1/3/ALDS1_3_Aより



Stackの実装


#include <iostream>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>

using namespace std;

struct Stack{
int S[1000];
int Top;
};

void push(Stack* stack, int num){
stack->S[stack->Top] = num;
stack->Top++;
}

int pop(Stack *stack){
stack->Top--;
return stack->S[stack->Top];
}

int main(void){
Stack stStack = {0};
char s[100];
int num1, num2;
while(scanf("%s", s) != EOF){
if (s[0] == '+'){
num1 = pop(&stStack);
num2 = pop(&stStack);
push(&stStack, num2+num1);
}
else if(s[0] == '-'){
num1 = pop(&stStack);
num2 = pop(&stStack);
push(&stStack, num2-num1);
}
else if(s[0] == '*'){
num1 = pop(&stStack);
num2 = pop(&stStack);
push(&stStack, num2*num1);
}
else{
push(&stStack, atoi(s));
}
}

num1 = pop(&stStack);
cout << num1 << endl;

return 0;
}




C++でバブルソートの問題であるAIZU ONLINE JUDGEのALDS_1_2_Dを解きました。



問題


1  insertionSort(A, n, g)
2      for i = g to n-1
3          v = A[i]
4          j = i - g
5          while j >= 0 && A[j] > v
6              A[j+g] = A[j]
7              j = j - g
8              cnt++
9          A[j+g] = v
10
11 shellSort(A, n)
12     cnt = 0
13     m = ?
14     G[] = {?, ?,..., ?}
15     for i = 0 to m-1
16         insertionSort(A, n, G[i])

shellSort(A, n) は、一定の間隔 gg だけ離れた要素のみを対象とした挿入ソートである insertionSort(A, n, g) を、最初は大きい値から gg を狭めながら繰り返します。これをシェルソートと言います。

上の疑似コードの ? を埋めてこのプログラムを完成させてください。nn と数列 AA が与えられるので、疑似コード中の mmmm 個の整数 Gi(i=0,1,...,m1)Gi(i=0,1,...,m1)、入力 AAを昇順にした列を出力するプログラムを作成してください。




シェルソートを実装する


ソートを行う間隔gは最後が1で終わればどんな数列でも良さそうだったので今回は、gn+1 = 3gn + 1という数列を使用します。

この数列に2のべき乗のような間隔が短すぎるものを入れると無駄な計算が発生しやすくなり、
シェルソートで実装する意味がなくなるみたいです。

#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

int cnt = 0;

void insertionSort(int *A, int n, int g){
for(int i=g; i<n; ++i){
int v = A[i];
int j = i-g;
while ((j >= 0) && (A[j] > v)){
A[j+g] = A[j];
j = j - g;
cnt++;
}
A[j+g] = v;
}
}

void shellSort(int* A, int n){
cnt = 0;
int m = 1;
vector<int> G;
G.push_back(1);
while(3*G[m-1]+1 <= n){
G.push_back(3*G[m-1]+1);
m++;
}

cout << m << endl;
for(int i=G.size()-1;i>=0; --i){
cout << G[i];
if (i != 0){
cout << " ";
}
}
cout << endl;

for(int i = G.size()-1; i >= 0; --i){
insertionSort(A, n, G[i]);
}
}

int main(void){
int n;
cin >> n;
int *A = new int[n];
for(int i=0; i<n; ++i){
cin >> A[i];
}
shellSort(A, n);
cout << cnt << endl;
for(int i=0; i<n; ++i){
cout << A[i] << endl;
}

return 0;
}



C++でバブルソートの問題であるAIZU ONLINE JUDGEのALDS_1_2_Aを解きました。



問題


数列 A を読み込み、バブルソートで昇順に並び変え出力するプログラムを作成してください。また、バブルソートで行われた要素の交換回数も報告してください。


入力例


入力の最初の行に、数列の長さを表す整数 N が与えられます。2行目に、N 個の整数が空白区切りで与えられます。

5
5 3 2 4 1


出力例


出力は 2 行からなります。1行目に整列された数列を 1 行に出力してください。数列の連続する要素は1つの空白で区切って出力してください。2 行目に交換回数を出力してください。
1 2 3 4 5
8

https://onlinejudge.u-aizu.ac.jp/courses/lesson/1/ALDS1/2/ALDS1_2_Aより



バブルソートを実装する


#include <iostream>
#include <cmath>

using namespace std;

void bubbleSort(int *A, int n){
int count = 0;
int flag = 1;
while(flag){
flag = 0;
for(int i=n-1; i>=0; --i){
if (A[i] < A[i-1]){
swap(A[i], A[i-1]);
count++;
flag = 1;
}
}
}
for(int i=0; i<n-1; ++i){
cout << A[i] << " ";
}
cout << A[n-1] << endl;
cout << count << endl;
}

int main(void){
int n;
cin >> n;
int* A = new int[n];
for (int i=0; i<n; ++i){
cin >> A[i];
}
bubbleSort(A, n);
delete[] A;
return 0;
}


↑このページのトップヘ