π

1. プログラムの概要
- 1.1. Machinの公式
- 1.2. 相加相乗平均
2. C
- 2.1. Machinの公式
- 2.2. 相加相乗平均
3. Emacs Lisp
- 3.1. Machinの公式
- 3.2. 相加相乗平均
4. Bash
- 4.1. Machinの公式
- 4.2. 相加相乗平均
5. Fish
- 5.1. Machinの公式
- 5.2. 相加相乗平均

Code of Poem

Unknown programmer's programming note.

<2022-09-26 月>

1. プログラムの概要

1.1. Machinの公式

pi1 - Machinの公式を使って円周率を求める

Machinの公式：

\(\frac{\pi}{4} = 4 \arctan{\frac{1}{5}} - \arctan{\frac{1}{239}}\)

arctanは、次の級数展開を使えばよい：

\(\arctan{x} = x - \frac{x^3}{3} + \frac{x^5}{5} - \frac{x^7}{7} + ...\)

1.2. 相加相乗平均

pi2 - 相加相乗平均を使って円周率を求める

\(\pi_n = \frac{2 a_{n+1}^2}{1 - \sum_{k=0}^{n} {2^k (a_k^2 - b_k^2)}} = \frac{2 a_{n+1}^2}{1 - \sum_{k=0}^{n}{2^k (a_k - a_{k-1})^2}}\)

2. C

2.1. Machinの公式

/*
 * 出典元「C言語による最新アルゴリズム事典」奥村晴彦著 - P16
 * インデントは変更してあります。
 * 他にも、メッセージや変数名や処理を変更している場合があります。
 */

#include <stdio.h>

int main()
{
    int k;
    long double p, t, last;

    p = 0;
    k = 1;
    t = 16.0L / 5.0L;

    do
    {
        last = p;
        p += t / k;
        t /= -5.0L * 5.0L;
        k += 2;
    } while (p != last);

    k = 1;
    t = 4.0L / 239.0L;

    do
    {
        last = p;
        p -= t / k;
        t /= -239.0L * 239.0L;
        k += 2;
    } while (p != last);

    printf("%.30Lf\n", p);
}

2.2. 相加相乗平均

/*
 * 出典元「C言語による最新アルゴリズム事典」奥村晴彦著 - P16
 * インデントは変更してあります。
 * 他にも、メッセージや変数名や処理を変更している場合があります。
 */

#include <math.h>
#include <stdio.h>

int main()
{
    int i;
    double a, b, s, t, last;

    a = 1;
    b = 1 / sqrt(2.0);
    s = 1;
    t = 4;

    for (i = 0; i < 3; i++)
    {
        last = a;
        a = (a + b) / 2;
        b = sqrt(last * b);
        s -= t * (a - last) * (a - last);
        t *= 2;
        printf("%16.14f\n", (a + b) * (a + b) / s);
    }
}

3. Emacs Lisp

3.1. Machinの公式

(let ((k 1)
      (p 0)
      (x (/ 16.0 5.0))
      (last))

  (while (progn
           (setq last p
                 p (+ p (/ x k))
                 x (/ x (* -5.0 5.0))
                 k (+ k 2))
           (not (= p last))
           ))

  (setq k 1
        x (/ 4.0 239.0))

  (while (progn
           (setq last p
                 p (- p (/ x k))
                 x (/ x (* -239.0 239.0))
                 k (+ k 2))
           (not (= p last))
           ))

  (message "%.20f" p))

whileが奇妙な形に見えるかもしれません。これは do-while あるいは repeat-until を真似るためのものです。

do-while emulation | Iteration (GNU Emacs Lisp Reference Manual)

3.2. 相加相乗平均

(let ((a 1)
      (b (/ 1.0 (sqrt 2.0)))
      (s 1)
      (x 4))

  (message (with-output-to-string
             (dotimes (i 3)
               (setq last a
                     a (/ (+ a b) 2)
                     b (sqrt (* last b))
                     s (- s (* x (* (- a last) (- a last))))
                     x (* x 2))
               (princ (format "%16.14f\n" (/ (* (+ a b) (+ a b)) s)))))))

4. Bash

bashは浮動小数点数の計算ができません。 bcにお願いするのが代表的な解決策です。

4.1. Machinの公式

#!/bin/bash

function fmath {
    echo "scale=20; $1" | bc
}

p=0
k=1
t=$(fmath "16.0 / 5.0") 

while
    last=$p
    p=$(fmath "$p + $t / $k")
    t=$(fmath "$t / (-5.0 * 5.0)")
    k=$(( $k + 2 ))
    test $(fmath "$p != $last") -eq 1
do true; done

k=1
t=$(fmath "4.0 / 239.0")

while
    last=$p
    p=$(fmath "$p - $t / $k")
    t=$(fmath "$t / (-239.0 * 239.0)")
    k=$(( $k + 2 ))
    test $(fmath "$p != $last") -eq 1
do true; done

printf "%.30Lf\n" $p

whileループが異様な形をしています。これは、必ず一回は本体を実行してから、最後に条件を判定する do-while あるいは rpeat-until を真似るためにこのような形になっています。

Is there a do-while loop in bash? - Stack Overflow

4.2. 相加相乗平均

#!/bin/bash

function fmath {
    echo "scale=20; $1" | bc
}

a=1
b=$(fmath "1.0 / sqrt(2.0)")
s=1
t=4

for i in {1..3}
do
    last=$a
    a=$(fmath "($a + $b) / 2.0")
    b=$(fmath "sqrt($last * $b)")
    s=$(fmath "$s - $t * ($a - $last) * ($a - $last)")
    t=$(fmath "$t * 2")
    printf "%16.14f\n" $(fmath "($a + $b) * ($a + $b) / $s")
done

5. Fish

fishのmathは、デフォルトで浮動小数点数の計算を行います。 bcに頼らなくても大丈夫でした。

5.1. Machinの公式

#!/usr/bin/fish

function fmath
    math -s15 $argv
end

set p 0
set k 1
set t (fmath 16.0 / 5.0)

while true
    set last $p
    set p (fmath $p + $t / $k)
    set t (fmath "$t / (-5.0 * 5.0)")
    set k (math $k + 2)
    if test $p -eq $last
        break
    end
end

set k 1
set t (fmath 4.0 / 239.0)

while true
    set last $p
    set p (fmath $p - $t / $k)
    set t (fmath "$t / (-239.0 * 239.0)")
    set k (math $k + 2)
    if test $p -eq $last
        break
    end
end

printf "%.30Lf\n" $p

fishでは do-while あるいは repeat-until をエミュレートする手段がないようです。代替として、ループの条件をwhileで判定するのではなく、本体の最後でifで判定して、breakするようにしました。

5.2. 相加相乗平均

#!/usr/bin/fish

function fmath
    math -s15 $argv
end

set a 1
set b (fmath "1.0 / (sqrt 2.0)")
set s 1
set t 4

for i in (seq 3)
    set last $a
    set a (fmath "($a + $b) / 2.0")
    set b (fmath sqrt $last x $b)
    set s (fmath "$s - $t * ($a - $last) * ($a - $last)")
    set t (fmath $t x 2)
    printf "%16.14f\n" (fmath "($a + $b) * ($a + $b) / $s")
end

π

目次