最終更新日
最終更新日 |
|
|
円周率πの計算には、マチンの方法 16*ARCTAN (1/5) - 4*ARCTAN (1/239) というのが知られている。まず、double型で計算してみる。doubleは符号つき64ビットだから、その値をまず調べてみた。1-3で述べるが、BigIntegerクラスを使ったら、63bit= 9223372036854775808=9.2*10^18 とでた。
まずは、Mathクラスのatan2(アークタンジェント)を使って、マチンの公式そのもので計算し、続いて、テイラー展開したもので、やや複雑な繰り返し処理で計算。さらに、1-4で、それをBigIntegerでやってみる。
class Pai{
public static void main(String args[]) {
double pai=16.0*Math.atan2(1.0,5.0)-4.0*Math.atan2(1.0,239.0);
System.out.println("pai = " + Math.PIpai);
System.out.println("Math.PI = " + Math.PI);
}
}
pai = 3.1415926535897936
Math.PI = 3.141592653589793 (Mathクラスの静的変数)
アークタンジェントの展開形は、
ARCTAN (X) = X (1 - 1 / 3X^2 + 1 / 5 X^4- 1 / 7 X^6 +........)
だから、マチンの公式は、次の2つになる。
16 * ARCTAN (1/5) = 16 / 5 (1 - 1 / (3 * 25) + 1 / (5 * 25^2) -.....)
4 * ARCTAN (1/239)= 4/239 (1 - 1 /(3*57121) + 1 /(5*57121^2) - ........)
(注 239*239=57121)
2項目を double n1=-1.0/(3.0*25.0); とすると、
3項目はi=2として、n1=n1*(-(2.0*i-1)/((2.0*i+1)*25.0));
という関係があり、以下もiを1ずつ増やしていくことでこの関係は保たれるので、whileによる繰り
返しにこの関係を使う。18桁がいいところなので、1/10^18をsとして、そこまでくりかえす。
class Pai{
public static void main(String args[]) {
double n1=-1.0/(3.0*25.0);
double nn1=1.0+n1;
double n2=-1.0/(3.0*57121.0);
double nn2=1.0+n2;
System.out.println("n1= " + n1 + " n2= " + n2);
double s=0.000000000000000001;
double i=2.0;
while (Math.abs(n1) >s) {
n1=n1*(-(2.0*i-1)/((2.0*i+1)*25.0));
nn1 += n1;
i +=1.0;
}
System.out.println("i = " + i);
i=2.0;
while (Math.abs(n2) > s) {
n2=n2*(-(2.0*i-1)/((2.0*i+1)*57121.0));
nn2 += n2;
i +=1.0;
}
System.out.println("i = " + i);
double pai=(16.0/5.0*nn1-4.0/239.0*nn2);
System.out.println("pai = " + pai);
System.out.println("Math.PI = " + Math.PI);
}
}
i = 13.0 (マチンの式前半の繰り返し回数は2から12まで)
i = 5.0 (マチンの式後半の繰り返し回数は2から4まで)
pai = 3.1415926535897936 (1-1 atan2を使ったものと同じ結果)
Math.PI = 3.141592653589793
import java.math.*;
class Pai{
public static void main(String args[]) {
BigInteger a10 = new BigInteger("10");
BigInteger b2 = new BigInteger("2");
BigInteger n = b2.pow(63);
BigInteger m1 = a10.pow(20);
BigInteger n1 = m1.divide(n);
System.out.println("n = " + n + " n1 = " + n1);
}
}
結果は、 63bit= 9223372036854775808 で約9.2*10^18
10^20をこれで割ると10となる。
import java.math.*;
class Pai{
public static void main(String args[]) {
BigInteger a1 = new BigInteger("1");
BigInteger a2 = new BigInteger("2");
BigInteger a3 = new BigInteger("3");
BigInteger a4 = new BigInteger("4");
BigInteger a5 = new BigInteger("5");
BigInteger a16 = new BigInteger("16");
BigInteger a25 = new BigInteger("25");
BigInteger a239 = new BigInteger("239");
BigInteger a57121 = new BigInteger("57121");
BigInteger m1 = new BigInteger("-1");
BigInteger a10 = new BigInteger("10");
BigInteger n=a10.pow(25);
//式前半 16 * ARCTAN (1/5) = 16 / 5 (1 - 1 / (3 * 25) + 1 / (5 * 25^2) -.....)
//1項と2項の部分 (1 - 1 / (3 * 25))のみ
BigInteger n1=m1.multiply(n).divide(a3.multiply(a25));
BigInteger nn1=n.add(n1);
//式前半の3項以下を前項を使って求め求め、足し込む。
BigInteger i=new BigInteger("2");
while ((n1.abs()).doubleValue() > 1.0) {
BigInteger nc=a2.multiply(i).subtract(a1);
BigInteger nm=a2.multiply(i).add(a1);
n1 = m1.multiply(n1).multiply(nc).divide(nm.multiply(a25));
nn1=nn1.add(n1);
i=i.add(a1);
}
System.out.println("i = " + i); //繰り返し回数
//式後半 4 * ARCTAN (1/239)= 4/239 (1 - 1 /(3*57121) + 1 /(5*57121^2) - ....)
//1項と2項の部分 (1 - 1 / (3 * 57121))のみ
BigInteger n2=m1.multiply(n).divide(a3.multiply(a57121));
BigInteger nn2=n.add(n2);
//式後半の3項以下を前項を使って求め、足し込む。
i=new BigInteger("2");
while ((n2.abs()).doubleValue() > 1.0) {
BigInteger nc=a2.multiply(i).subtract(a1);
BigInteger nm=a2.multiply(i).add(a1);
n2=m1.multiply(n2).multiply(nc).divide(nm.multiply(a57121));
nn2 = nn2.add(n2);
i=i.add(a1);
}
System.out.println("i = " + i); //繰り返し回数
//式前半
BigInteger pai1=a16.multiply(n).multiply(nn1).divide(a5);
//式後半
BigInteger pai2=a4.multiply(n).multiply(nn2).divide(a239);
//全体
BigInteger pai=pai1.subtract(pai2);
System.out.println("pai = " + pai);
}
}
25桁 i = 18 i = 6
pai = 314159265358979323846264361974895397489539748953975
[参考] 有効桁数 25桁(Logo2writerでやったもの)
314159265358979323846264338327950288419716939937510582097494459237816406
2862898862834825332116679
100桁 i = 71 i = 22
pai = 3141592653589793238462643383279502884197169399375105820974944592307816406
2862089986280348253421170677774058577405857740585774058577405857740585774
0585774058577405857740585774058577405857740585774058578
10000桁
pai = 3141592653589793238462643383279502884197169399375105820974944592307816406
2862089986280348253421170679821480865132823066470938446095505822317253594
0812848111745028410270193852110555964462294895493038196442881097566593344
6128475648......以下略
1-2のdouble型をBigIntegerで多桁の整数計算するよう、BigIntegerクラスの変数をインスタンス化し、メソッドを使って置き換えたにすぎない。(クラスライブラリ概要(ネット上のJava2 SE5.0 API仕様 日本語版 参照)。
100桁なら、n=a10.pow(100) として、n倍しておいて計算結果の整数のみをとるものである。25桁、100桁、10000桁とやってみて、指定桁数まで一致しているといえない。どこかおかしいかもしれないが、繰り返し計算の中での誤差の累積であろう。[参考]としてだした有効桁数 25桁とした、Logo2writerでやったものは、余りを次々と計算して加えていくやり方をとったが、そこまでやる必要があるのだろう。
繰り返しの打ち切り判定の while ((n1.abs()).doubleValue() > 1.0) は、異様だが、BigIntegerクラスに、比較のメソッドがないため、doubleValue()でdoubleに変換してやっている。
二項定理、二項式のN 乗(a + b)^ N の展開式の係数を求めるもの。1行前の2つの係数を加えたものが
次の行の各項目になるという、不思議な面白さがある。
さらに、Nケの中からKケ選び出す場合の数にも相当し、「コアJava2 Vol1.1基礎編」改定版(ASCII 2001/12)
の P108では逆数が宝くじのあたる確率として紹介され、さらにP128では二次元配列の事例として紹介され
ている。このプログラムの一部は、2-3で紹介し、関連についても考察する。
[プログラム]
import java.awt.*;
import java.util.*;
import java.applet.Applet;
/*
<applet code="Pascal.class" width=800 height=500>
</applet>
*/
public class Pascal extends Applet {
public void paint(Graphics g) {
final int NMAX = 15;
int[][] odds = new int[NMAX + 1][];
for (int n = 0; n <= NMAX; n++)
odds[n] = new int[n + 1];
for (int n = 0; n < odds.length; n++){
odds[n][0] = 1;
if (n gt;= 1)
odds[n][odds[n].length-1] = 1;
}
for (int n = 2; n < odds.length; n++)
for (int k = 1; k < odds[n].length-1; k++) {
odds[n][k]=odds[n-1][k-1] + odds[n-1][k];
}
setBackground(Color.yellow);
g.setFont(new Font("Serif", Font.BOLD, 20));
g.setColor(Color.blue);
for (int n = 0; n < odds.length; n++) {
for (int k = 0; k < odds[n].length; k++) {
String output =(new Integer(odds[n][k])).toString();
g.drawString(output,400 - 25 * n + 45 * k, 30 + 30 * n);
}
}
}
}
[ソース]
Pascal.java
[実行結果]
[コメント]
1 前文にも書いたように、二項式のN 乗(a + b)^N の展開式の係数を求めている。
2 1 に始まり、(a + b) 、(a + b)^2------(a + b)^Nの係数だけを次々に図示し、数字だけで作った
三角形のようになり、パスカルの三角形といわれる。
3 各列は(a + b)^N の展開形で、a^N とb^N の係数はいずれも１、1 で始まり、1 で終わる。
その中は一つ上の列でその数の左右にある二つの数の和となる。
(a + b)^2 の場合を考えると、一つ前の(a + b) に(a + b) を掛けたもの、(a + b) の各項にa を掛け、
次にb を掛けて加え合わせたものになる。abの係数は、ab+baで2になるというわけ。
N 乗の場合も(a + b) ^ N-1 の各項に対して同様に考えればよい。
4 「コアJava2 Vol1.1基礎編」P126を参考にして、2次元配列を使った。
行数(配列の最初の[])を、int[][] odds = new int[NMAX + 1][];で決めた後、
行ごとの大きさ(列方向)は、各行のnの値を越えないとの制約で、
for (int n = 0; n <= NMAX; n++)
odds[n] = new int[n + 1];
としているのは、「コアJava2」P128による。0の選択も含めるので、n + 1(2-3のコメント1参照)。
5 両端の1は、前列の2項の和といかないため、それだけのためのforブロックを作り、強制的に1としている。
後述する、Logoのプログラム例は、列の両端に0を加えておいて、前列の2つの和を例外なくスムーズに
計算できるようにしている。その0は処理の便宜上、追加しただけのものなので、印字はしない。
[プログラム 一部省略]
//import文,appletタグ、2-1に同じ、省略
public class Pascal5 extends Applet {
public void paint(Graphics g) {
final int NMAX = 15;
int[][] odds = new int[NMAX +1][];
for (int n = 0; n <= NMAX; n++)
odds[n] = new int[n + 3]; // [n + 1]を変更
odds[0][1] = 1; //最初の項だけ1を入れておく。
for (int n = 0; n < odds.length; n++){
odds[n][0] = 0; //列の左端にはすべて0を設定
odds[n][odds[n].length-1] = 0; //列の右端もすべて0を設定
}
for (int n = 1; n < odds.length; n++) //n=1の場合もn=0をもとに計算できる
for (int k = 1; k < odds[n].length-1; k++) {
// System.out.println("n=" + n + " k=" + k);
odds[n][k]=odds[n-1][k-1] + odds[n-1][k];
}
// 文字、色設定、2-1に同じ(省略)
for (int n = 0; n < odds.length; n++) {
for (int k = 1; k < odds[n].length-1; k++) { //列の両端の0は印字しない
String output =(new Integer(odds[n][k])).toString();
g.drawString(output,400 - 25 * n + 45 * k, 30 + 30 * n);
}
}
}
}
[ソース]
Pascal5.java
2-1の前文にもふれた、Nケの中からKケ選び出す場合の数、「コアJava2 Vol1.1基礎編」P128にある
Application版リストの一部を紹介する。
[プログラム一部]
...........
for (int n = 0; n < odds.length; n++)
for (int k = 0; k < odds[n].length; k++) {
/*
compute binomial coefficient
n * (n - 1) * (n - 2) * . . . * (n - k + 1)
-------------------------------------------
1 * 2 * 3 * . . . * k
*/
int lotteryOdds = 1;
for (int i = 1; i <= k; i++)
lotteryOdds = lotteryOdds * (n - i + 1) / i;
odds[n][k] = lotteryOdds;
}
...........
[ソース]
Pascal2.java (「コアJava2」P128をApplet版に書き換えたもの)
[コメント]
1 (a + b)^n の展開式の係数と、上記のコメントの部分との関係を考えてみた。
a,b,cの3ケからなる、n=3 をみる。0ケ選ぶというのも1としているから、1 3 3 1。
k=1が3(a または bまたはc の3通り)
k=2も3(ab または bc またはca の3通り)
a,b,c,dの4ケからなる、n=4 をみる。1 4 6 4 1。
2 n=4,k=3の場合が、n=3,k=1の場合とn=3,k=2の場合の和になるというのが、パスカルの三角形。
その関係を、蛇足ながら、上記コメント内の式で記述してみると、
第1項、n-1,k-1の場合、
(n - 1) * (n - 2) * . . * (n - k + 1)
-------------------------------------------
1 * 2 * 3 * . . . * (k-1)
第2項、n-1,kの場合、
(n - 1) * (n - 2) * .. * (n - k + １) * (n - k)
------------------------------------------------
1 * 2 * 3 * . . . * k
この2つの和、第1項の分母、分子に、kを乗じて整理すると、
{(n - 1) * (n - 2) * . . . * (n - k + 1)} * (k + n - k)
--------------------------------------------------------
1 * 2 * 3 * . . . * k
結局、n,kの場合と同じになる。
2-2でふれたLogoのプログラム例は、2-2に書いた方法をもっと巧妙に使っている。列の先頭に0を
加えたリストNL1と、列の末尾に0を加えたリストNL2を作る。
それぞれの先頭の和をとって、先頭に加えた上で、リストの先頭を除いて、再帰で繰り返していく
というものである。
なお、Logoについては、Javaプログラミングメモ3(クィックソート) の2章参照。
[参考ーLogoのプログラム]
TO START
CG PU SETPOS [-50 130] SETC 3 LABEL [パスカルの三角形]
SETC 4 MAKE "X 0 MAKE "Y 110
PASCAL [1]
END
TO PASCAL :NL
IF (COUNT :NL) = 14 [STOP]
PR :NL
MAKE "NL1 FPUT 0 :NL
MAKE "NL2 LPUT 0 :NL
MAKE "NL PASCAL1 :NL1 :NL2
PASCAL :NL
END
TO PASCAL1 :NL1 :NL2
IF :NL1 = [ ] [OP [ ]]
OP FPUT (FIRST :NL1) + (FIRST :NL2)
PASCAL1 BF :NL1 BF :NL2
END
TO 三角形を描く :NL
PU SETPOS LIST :X :Y
REPEAT COUNT :NL
[LABEL FIRST :NL RT 90 FD 35 LT 90 MAKE "NL BF :NL]
WAIT 50
MAKE "X :X - 20 MAKE "Y :Y - 20
END
できたPascalの三角形のそれぞれの数字を、たとえば7で割った余りの数に各種の色を割り当て、小さな
三角形を塗りつぶすなどすると、結構、きれいな図形ができあがる(2-6参照)。
さて、Javaのリスト処理、再帰用法は、Logoのようにスマートにはいかないが、トライしてみた。
[LogoのプログラムーApplication版]
import java.util.*;
public class Pascal10 {
public static void main(String[] args) {
final int NMAX = 15;
ArrayList alt = new ArrayList();
alt.add(0,new Integer(1)); //Integerをいれる
ArrayList alt1 = new ArrayList();
ArrayList alt2 = new ArrayList();
for (int i = 0; i < NMAX; i++) {
System.out.println("alt=" + alt); //結果の打ち出し
alt1=(ArrayList)alt.clone();
alt2=(ArrayList)alt.clone();
alt1 = fput(new Integer(0),alt1); //リスト処理を巧みにする工夫
alt2 = lput(new Integer(0),alt2); //リスト処理を巧みにする工夫
alt=pascal20(alt1,alt2);
}
}
public static ArrayList pascal20(ArrayList nl1,ArrayList nl2) {
if (nl1.isEmpty()) return nl1;
Integer n1=(Integer)nl1.get(0);
Integer n2=(Integer)nl2.get(0);
int f = n1.intValue()+ n2.intValue();
return fput (new Integer(f),pascal20(bf(nl1),bf(nl2)));
}
public static ArrayList bf (Object a1) {
//リストa1の最初の要素を取り去る。
ArrayList aa = (ArrayList)a1;
aa.remove(0);
return aa;
}
public static ArrayList fput (Object a1,Object a2) {
//リストa1の最初にa1を追加する
ArrayList aa = (ArrayList)a2;
aa.add(0,a1);
return aa;
}
public static ArrayList lput (Object a1,Object a2) {
//リストa1の最後にa1を追加する
ArrayList aa = (ArrayList)a2;
aa.add(a1);
return aa;
}
}
[ソース]
Pascal10.java
[コンバイル 警告]
javac Pascal10.java -Xlint:unchecked としないと、エラーになる。そして、こうすると、下記の警告。
Pascal10.java:26: 警告: [unchecked] raw 型 java.util.ArrayList のメンバとしてのadd(int,E)
への無検査呼び出しです。 alt.add(0,new Integer(1));
[実行結果]
alt=[1]
alt=[1, 1]
alt=[1, 2, 1]
alt=[1, 3, 3, 1]
alt=[1, 4, 6, 4, 1]
alt=[1, 5, 10, 10, 5, 1]
alt=[1, 6, 15, 20, 15, 6, 1]
alt=[1, 7, 21, 35, 35, 21, 7, 1]
alt=[1, 8, 28, 56, 70, 56, 28, 8, 1]
alt=[1, 9, 36, 84, 126, 126, 84, 36, 9, 1]
alt=[1, 10, 45, 120, 210, 252, 210, 120, 45, 10, 1]
alt=[1, 11, 55, 165, 330, 462, 462, 330, 165, 55, 11, 1]
alt=[1, 12, 66, 220, 495, 792, 924, 792, 495, 220, 66, 12, 1]
alt=[1, 13, 78, 286, 715, 1287, 1716, 1716, 1287, 715, 286, 78, 13, 1]
alt=[1, 14, 91, 364, 1001, 2002, 3003, 3432, 3003, 2002, 1001, 364, 91, 14, 1]
2-4の中のLogoライクなfput,bfメソッドをリスト処理のパッケージとしてまとめ、それを呼び出し
て使ってみた。
同様のことを、Javaプログラミングメモ3(クィックソート) の6章でも、やっているが、一歩前進させた。
以前は、fput,bfなど、メソッドですむところをわざわざ、独立のクラスにしていたが、この場合はその必要はない。
バッケージ名を慣例に従い、小文字にした。バッケージ化したクラスを納めるディレクトリ(フォルダ)名を
バッケージ名と同じ にするだけのこと(「独習Java」P214)。
バッケージを使う側でも、package宣言が必要というのは納得できないとの以前の疑問はset CLASSPATH指定を綿密に
することで解決できた。また、バッケージと、それを使うプログラムを別のディレクトリにおくこともできた。
(コメント 3,4参照)。
[プログラムの一部]
import java.util.*;
import listp.*;
public class Pascal100 {
static ListP lt = new ListP(); //List処理への参照変数、static要
public static void main(String[] args) {
final int NMAX = 15;
// ListP lt = new ListP(in,alt2);
ArrayList alt = new ArrayList();
alt.add(0,new Integer(1)); //第1項 [1]
for (int i = 0; i < NMAX; i++) {
System.out.println("alt=" + alt);
ArrayList alt1=(ArrayList)alt.clone(); //clone()使用
ArrayList alt2=(ArrayList)alt.clone();
alt1 = lt.fput(new Integer(0),alt1);
alt2 = lt.lput(new Integer(0),alt2);
alt=pascal20(alt1,alt2);
}
}
public static ArrayList pascal20(ArrayList nl1,ArrayList nl2) {
if (nl1.isEmpty()) return nl1;
Integer n1=(Integer)nl1.get(0);
Integer n2=(Integer)nl2.get(0);
int f=n1.intValue()+ n2.intValue();
/* Integer ff=new Integer(f);
ListP lf = new ListP(ff,nl1);
ListP lt1 = new ListP(nl1);
return lf.fput (new Integer(f),pascal20(lt1.bf(nl1),lt1.bf(nl2)));
*/
return lt.fput (new Integer(f),pascal20(lt.bf(nl1),lt.bf(nl2)));
}
}
[ソース]
Pascal100.java
package listp(listpディレクトリにある)
ListP.java
[コメント]
1 実行結果は、当然、2-4 と同じである。
2 当初は、ListPクラスのさまざまな引数に対応したコンストラクタを使ったオブジェクト参照変数を用意
した。 上記のコメントの部分である。
それは煩雑なので、引数なしのコンストラクタで一本化をトライしたら、あっさりできてしまった。
冒頭の static ListP lt = new ListP(); である。staticが必要になるのは、嬉しくないが。
3 以上は、パッケージlistpもそれを使うPascal100.javaも同じディレクトリPascalにおいた場合だが、
listpはPascalに、Pascal200.javaはAppletにおいた場合、set CLASSPATH=Pascal として、コンパイル
はできるが、実行すると、下記のエラーメッセージが出る。
Exception in thread "main" java.lang.NoClassDefFoundError: Applet/Pascal200
4 しかし、set CLASSPATH=Pascal;Applet; とすることで、この問題も解決した(「独習Java」P217参照)。
2-4にちらっと書いた図形化を試みた。2で割った余りが0か1で、色を変えて菱形を塗りつぶすだけ
の簡単なもの。円もトライはしてみた。
Lineart.java
[プログラム]
//import文,appletタグ、2-1に同じ、省略
public class Lineart extends Applet {
public void paint(Graphics g) {
//バスカルの三角形 背景色、すべて2-1に同じ、省略
//ただし、文字を描くことはないので、文字の大きさ、色指定は削除
int x;
int y;
for (int n = 0; n < odds.length; n++) {
for (int k = 0; k < odds[n].length; k++) {
x= 400 - 15 * n + 30 * k;
y=40 + 30 * n;
if (odds[n][k] % 2 == 0)
g.setColor(Color.red);
else
g.setColor(Color.blue);
int[] xdata={x,x-15,x,x+15};
int[] ydata={y,y+10,y +20,y+10};
g.fillPolygon(xdata,ydata,4);
}
}
}
}
[ソース]
Lineart.java
未完
 |
|
|