π DouKaku?

円周率を計算してください。

積分を計算するも、素朴な方法も、速さを目指すも、LLで計算する意味を問うもあるでしょう。

http://ja.wikipedia.org/wiki/%E5%86%86%E5%91%A8%E7%8E%87

とりあえず、参考ページにあった数式で素直に計算。

use strict;
use warnings;

use bignum;

use List::Util qw/reduce/;
use Memoize;

sub fract
{
  $_[0] == 0 || $_[0] == 1 ? 1 :
  $_[0] * fract($_[0]-1);
}
memoize(q/fract/);

my $N = 32;

print 6 * reduce {
  no warnings qw/once/;
  $a + $b
} map {
  fract(2 * $_) /
  ((2 ** (4 * $_ + 1)) * (fract($_) ** 2) * (2 * $_ + 1))
} 0 .. $N;

ビジュアル的にはこっちのほうがおもしろい

use strict;
use warnings;

use bignum;

use Memoize;

sub fract
{
  my $n = shift;
  $n == 0 || $n == 1 ? 1 :
  $n * fract($n-1);
}
sub f
{
  my $n = shift;
  $n < 0 ? 0 :
  fract(2 * $n) /
  ((2 ** (4 * $n + 1)) * (fract($n) ** 2) * (2 * $n + 1)) +
  f($n-1);
}
memoize(q/fract/);
memoize(q/f/);

my $N = 1024;

print "@{[6 * f($_)]}\n" for 1 .. $N;

ここでいう「計算」とはどういう意味なのでしょうか。

遅い、精度が低い

import java.util.Random;
/**
 *　教科書的なモンテカルロ法
 */
public class Pi {
    public static void main(String args[]) {
    Random r= new Random();
    int m=1000000;
    double j=0;
    for(int i=0;i<m;i++)
        if(Math.hypot(r.nextDouble(),r.nextDouble())<1.0)
          j++;
    System.out.println(4*j/m);
    }
}

問題文の趣旨をまるで無視。

1	N[Pi,10000]

Rにはpiという組み込み変数がありますが、計算であればRらしくモンテカルロで。

1
2
3

pi.montecarlo <- function(n){
  mean(replicate(n, sum(runif(2)^2) < 1)) * 4
}

Squeak Smalltalk で。

| pi n next |
pi := 0.0.
n := 0.
[   next := 6 * (2 * n) factorial 
        / ((2 raisedTo: 4 * n + 1) * n factorial squared * (2 * n + 1)).
    next < 1e-16] whileFalse: [pi := pi + next. n := n + 1].
^pi   "=> 3.141592653589793 "

昔のメモにのこっていたコード。出典がわからないのだがなにかのpaperにあったコードだと思う。
桁数を指定する。各桁の数のリストで出力

*Main> :main 30
[3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3,3,8,3,2,7]

import System.Environment
main = print . flip take  (g (1,180,60,2)) . read . head =<< getArgs
g (q,r,t,i) = y : g (10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
  where (u,y) = (3*(3*i+1)*(3*i+2),div (q*(27*i-12)+5*r) (5*t))

言語指定しわすれた Haskell です

すみません再投稿します

import System.Environment
main = print . flip take  (g (1,180,60,2)) . read . head =<< getArgs
g (q,r,t,i) = y : g (10*q*i*(2*i-1),10*u*(q*(5*i-2)+r-y*t),t*u,i+1)
  where (u,y) = (3*(3*i+1)*(3*i+2),div (q*(27*i-12)+5*r) (5*t))

せっかくmathematicaなんだから積分計算してくださいよ。

それから何桁まで埋め込みでもっているのでしょうか？

Mathematica の Pi は組み込みの無理数オブジェクトで、N は引数を任意の精度の実数に変換する関数なのだそうです。

see: Mathematica への誘い - ２．数値計算

前半は半径rの1/4円の面積aを求めています。
後半では面積aを単純に4倍し、r*rで割ってpiを求めたかったのですが、
bashでは小数を持てないので、小数部分を10倍しながら表示しています。

$ ./pi.sh 100
3.1416
$ ./pi.sh 241
3.1400974
$ ./pi.sh 1000
3.141548

#/bin/bash

r="$1"
a=0

for ((x=1;x<=r;x++)) {
    for ((y=0;y<r;y++)) {
        if ((x*x+y*y<=r*r));then
            a=$((a+1))
        else break
        fi
    }
}

z=$((a*4));m=$((z/(r*r)))
echo -n $m
z=$((z-m*r*r))
if ((z));then
    echo -n .
fi
for ((i=1;i<8&&z>0;i++)) {
    n=$((z*10/(r*r)))
    z=$((z*10-n*r*r))
    echo -n $n
}
echo

１万回ループの無限級数(1-1/3+1/5-1/7…)で三桁の精度です。

# 組み込みPI
PIを表示

# 無限級数
a=1
nで1から10000まで繰り返す
    a=a+(-1^n)*(1/(2*n+1))
a*4を表示

# モンテカルロ法
t=10000
m=10000
a=0
nで0からtまで繰り返す
    もし(HYPOT(mの乱数/m,mの乱数/m))<1ならば
        a=a+1
a*4/tを表示

単純に円周率の数字がほしいだけならHaskellには最初から用意されています。

1	main = print pi

Gauche で (use math.const) (print pi) というのも風情がないので三角関数で。

1 2	(display (* (asin 1) 2)) (newline)

なるほど。どういう実装なんでしょうね。

  Mathematicaにはいくつかの数学定数が組み込まれています。これらのうち興味深いものは
  無理数で、それらはMathematicaではいかなる精度でも計算できるように定義されています。
  RAMの容量と計算を求めるのにかかる時間だけが、唯一計算可能な精度を制限しています。
  数学定数の中でも最も良く知られているのはPiとEでしょう。

Wikipediaにあったarcsinのテーラー展開を利用してます。

整数だけで計算してみました。
またよくわからない計算になってしまいました。

#include <stdio.h>

/*
       ∞
    pi=Σ(3*(2n)!/(2^4n*(n!)^2*(2n+1)))
       n=0
*/
void put_pi(){
    unsigned long long pi=0;
    int n;

    unsigned long long factor=1;
    int f2_m=2;
    int f2_m_a=10;
    int pow2_4n=0;
    int p2f_a=1;
    int add2n_1=1;
    
    for(n=0;n<16;n++){
        pi+=(3*(1000000000000000000>>pow2_4n)*(factor)/add2n_1);
        factor=factor*f2_m/p2f_a;
        f2_m+=f2_m_a;
        f2_m_a+=8;
        pow2_4n+=4;
        add2n_1+=2;
        p2f_a+=add2n_1;
        
        printf("PI=%Ld\n",pi);
    }
}
int main(){
    put_pi();
}

arcsin(1/2)のテイラー展開により円周率を一万桁まで求めるプログラムです。

最後の20桁は05600101655256367493です（9995桁まで正確に求められているようです）。

import java.math.*;

public class Pi10000 {
    static final BigInteger TWO = BigInteger.valueOf(2);
    static final BigInteger SIXTEEN = BigInteger.valueOf(16);
    static final int PRECISION = 10000;
    
    public static void main(String[] args) {
        BigInteger pi = BigInteger.ZERO;
        BigInteger a = BigInteger.valueOf(3).multiply(BigInteger.TEN.pow(PRECISION));
        BigInteger d = BigInteger.ONE;
        BigInteger e;
        long i = 0, j = 0;
        do {
            e = a.divide(d);
            pi = pi.add(e);
            a = a.multiply(BigInteger.valueOf(++i)).multiply(BigInteger.valueOf(++i));
            a = a.divide(SIXTEEN).divide(BigInteger.valueOf(++j).pow(2));
            d = d.add(TWO);
        } while(!e.equals(BigInteger.ZERO));
        System.out.println(new BigDecimal(pi, PRECISION));
    }

}

せっかくなので数値積分させてみました．

Abs[Pi - npi]を評価すると2x10^(-100)でした．

この値は積分範囲の逆数程度になっている様子．

1	npi = NIntegrate[1/(1 + x^2), {x, -10^100, 10^100}, WorkingPrecision -> 200, MaxRecursion -> 50]

VBA for Excel (2003)

Sub PI()
    Range("A2").FormulaR1C1 = "=ROW(RC)-2"
    Range("B2").FormulaR1C1 =  _
        "=FACT(RC[-1]*2)/FACT(RC[-1])/FACT(RC[-1])" & _
        "/POWER(2,4*RC[-1]+1)/(2*RC[-1]+1)"
    Range("A2:B22").Select
    Selection.FillDown
    Range("A1").Select
    ActiveCell.FormulaR1C1 = "=SUM(R[1]C[1]:R[21]C[1])*6"
End Sub

横に長いと醜いので改行とインデントを入れている。

see: Google ヘルプ:電卓機能

string s = new System.Net.WebClient().DownloadString
  ("http://www.google.co.jp/search?hl=en&q=pi");
System.Console.WriteLine
  (s.Substring(s.IndexOf("pi =") + 5, 10));

今更ながらfactorialのつもりのサブルーチン名がfractになってることに気が付いた。
恥さらしもいいとこだ orz...

arcsinのテイラー展開を写経する的な

see:

import (rnrs hashtables (6)))

(define (memoize f)
  (let ((cache (make-eqv-hashtable)))
    (lambda args
      (let ((pvalue (hashtable-ref cache args #f)))
    (or    pvalue
        (let ((result (apply f args)))
          (hashtable-set! cache args result)
          result))))))

(define fact
   (lambda (n)
     (if (<= n 1)
     1
     (* n (fact (- n 1))))))

(define (pi freq)
  (define (iter n value)
    (let ((value2 (+ value (taylor-step n))))
      (if (= n freq)
      value2
      (iter (+ n 1) value2))))
  (define (taylor-step n)
    (/ (fact (* 2 n))
       (*
    (expt 2 (+ (* 4 n) 1))
    (expt (fact n) 2)
    (+ (* 2 n) 1))))
  (* (iter 0 0) 6))

(define result (pi 30))
(display result)
(newline)
(display (inexact result))

1 2	import math print math.asin(1)*2

高野喜久雄の公式を使って高速化してみました。

π/4 = 12 arctan (1/49) + 32 arctan (1/57) - 5 arctan(1/239) + 12 arctan(1/110443)

arcsin(1/2)に比べて3.4倍高速になっています（小数点以下10万桁求めています。99,998桁まで正確に求められているようです）。

import java.math.*;

public class Pi100000 {
    static final BigInteger TWO = BigInteger.valueOf(2);
    static final int PRECISION = 100000;
    
    public static void main(String[] args) {
        BigInteger pi = BigInteger.ZERO;
        pi = arctan(pi, 4 * 12, 49);
        pi = arctan(pi, 4 * 32, 57);
        pi = arctan(pi, -4 * 5, 239);
        pi = arctan(pi, 4 * 12, 110443);
        System.out.println(new BigDecimal(pi, PRECISION));
    }
    
    public static BigInteger arctan(BigInteger pi, int n, int x) {
        BigInteger a = BigInteger.valueOf(n).multiply(BigInteger.TEN.pow(PRECISION));
        a = a.divide(BigInteger.valueOf(x));
        BigInteger x2 = BigInteger.valueOf(x).pow(2);
        BigInteger d = BigInteger.ONE;
        for (;;) {
            BigInteger e = a.divide(d);
            if (e.signum() == 0)
                break;
            pi = pi.add(e);
            a = a.divide(x2);
            d = d.add(TWO);
            e = a.divide(d);
            pi = pi.subtract(e);
            a = a.divide(x2);
            d = d.add(TWO);
        }
        return pi;
    }

}

面積から攻めるのはよく知られていますね。

「ビュフォンの針」を実装するのもありでしょう。

手でがんばった人

http://portal.nifty.com/2008/10/01/c/

applet

http://www.f.waseda.jp/takezawa/math/number/buffon/buffon.htm

なるほど。googleに計算させるのですね。

http://3.141592653589793238462643383279502884197169399375105820974944592.jp/ というのもあるらしいですが、サーバが死んでいるみたいです。

print math.atan(1)*4 でもいいですな。

IOCCC に円周率を計算するプログラムが投稿されたことがありますね http://www0.us.ioccc.org/1988/westley.c

see:

IOCCCって、AA好きですよね。＃フライトシミュレーターでソースが飛行機のAAだったのがあったな。

common-lisp:pi と比べると意外と精度が低い？

1	(princ (* (acos 0) 2))

「ビュフォンの針」は初めて知りました。
というわけで実装。

class Needle
  LENGTH = 5
  RAND_MAX = 10000
  attr_reader :x1, :x2

  def initialize()
    @x1 = rand(RAND_MAX) + rand
    @x2 = LENGTH * Math.cos(rand * 2 * Math::PI) + @x1
  end

end

n = 10000
crossed = 0
n.times {
  needle = Needle.new()
  k1 = (needle.x1 / (2 * Needle::LENGTH)).to_i
  k2 = (needle.x2 / (2 * Needle::LENGTH)).to_i
  crossed += 1 if k1 != k2
}

pi = n.to_f / crossed.to_f
puts "pi = #{pi}"

> 昔のメモにのこっていたコード。出典がわからないのだがなにかのpaperにあったコードだと思う。

コードの出典はわかりませんが計算方法は多分
「円周率の公式集（そのた・不明）」のページ
http://www.pluto.ai.kyutech.ac.jp/plt/matumoto/pi_small/node14.html
の下の方にある

t[n] = n(2n-1)((5n+3)+t[n+1])/((3n+1)(3n+2)3)
π[n] = 3 + t[1]
π = lim π[n]
というやつですね。

> 各桁の数のリストで出力

収束の関係上いつかは誤差で間違った数がでてくるんでしょうね
（後の桁で１０以上の数がでてきて全体としては正しい値に調整されると思いますが）
その桁にたどりつく前に計算機資源が足りなくなるか。

DouKaku?

π

Posted feedbacks - Flatten

Index

Feed

Other

Link

Pathtraq