ヒッパルコスの弦の数表

このページは元々

プトレマイオスの弦の数表

の訳注に含めていましたが、長くなりすぎたので、現在の独立のページにしました。

ヒッパルコスの弦の表は現在は失われていますが、 7.5° ごとのものであったことが知られています。最初は、あまりに大ざっぱな表で大したことがないと思っていたのですが

大きさと距離について

を読んでいくうちに、ヒッパルコスは角度 α が小さい時に、α の弦 chord(α) を求めるのに次の近似式を利用していることに気がつきました。(これはトゥーマーによる推論だと思われます。)

chord(α)＝2 sin(α/2)≒ α×60/3438

ヒッパルコスの近似式は、円周率を 3.1416 とした時、角度 α が小さければ、中心角が α の弦の長さと円弧の長さがほぼ等しいとしたときに得られるものです：

sin(α)≒(α/180)×3.1416≒α×60/3438

最後の近似式が成立することは電卓で確認することができます。

従って、加法定理を併用することにより、ヒッパルコスは任意の角 θ の弦をかなり高い精度で求めていたことになります。従って、ヒッパルコスは加法定理を知っていたことも明白となります。

注意

角度 α の弦 chord(α) を

chord(α) = 2 sin(α/2)

と考える方が理解しやすいですが、ヒッパルコスもプトレマイオスも円の半径を R に選び

R chord(α)

を角度 α の弦と考えています。つまり半径 R の円の中心角 α の「弦の長さ」を「弦」と呼んでいるわけです。ヒッパルコスの場合は R = 3438/60 (60 進表示では 57; 18) に選び、プトレマイオスの場合には R = 60 に選んでいるのです。以上は「カッツ、数学の歴史」による。なお、「カッツ、数学の歴史」にはヒッパルコスが円周率を 60 進表示で 3;8,30 としており、この値はアルキメデスの 2 つの値 3¹⁰/₇₁ と 3¹/₇ の平均値に近いと書いています。

3;8,30 を 10 進表示すると 3.14166.. となり、3.1416 より少し大きくなります。円周率を 60 表示すると 3;8,29,44,0,.. となり小数点以下 3 桁を丸めると 3;8,30 となります。また 3.1416 を 60 進表示して丸めると 3;8,30 となるので、ヒッパルコスの円周率は 3.1416 であるとして良いと思います。

最初はヒッパルコスが誤差の評価までしているとは思わなかったのですが、プトレマイオスの弦の数表では至る所で厳密な誤差を評価しているように思うことになったので、ヒッパルコスも厳密に誤差を評価しているはずだと考えて、まとめたのがヒッパルコスの弦の表はどの程度正確か ?です。

ヒッパルコスの方法の一番問題となりそうなのが円周率 π が 3.1416 であることです。アルキメデスは円に内接する正角形の周囲の長さと、円に外接する正多角形の周囲の長さから、 n = 96 の場合に π = 3.14 程度の値を得ているそうです。ヒッパルコスは更に三角法を使用できるので、アルキメデスよりはずっと効率よく π = 3.1416 を決定することができたはずだと考えて正多角形で円周率 π = 3.1416 を求める方法をまとめました。

正多角形で円周率 π = 3.1416 を求める方法でしていることは、単にアルキメデスの方法に三角法を援用しただけのものです。しかしこれを実際に試してみると、この方法はプトレマイオスが弦の数表を計算するときに、実施した方法と類似のものとなります。即ち、ある角度 θ が与えられた時、cos の倍角公式を利用して、 cos(θ×2^-n) (n = 1, 2, .. ) を順に決定し、ある精度に達した時に、対応する角度の sin を求めることになります。

このような考察を進めていくうちに、プトレマイオスは弦の数表を作成する際に π=3;8,30 となることを再チェックし、これを弦の数表の検算に使用したことに気がつきました。

π=3;8,30 はヒッパルコスが発見したことではないかもしれませんが、円周率 π=3;8,30 はヒッパルコスの理論にとって重要な役割を果たしますから、恐らくヒッパルコスも π=3;8,30 を再チェックしたのに相違なく、正 768 角形を使用したのに相違ないと思われます。

以下では n が半整数でも負の数でも、t が微小量であれば

(1 + t)ⁿ ≒ 1 + nt

の近似が古代にも知られていたと仮定しています。

ヒッパルコスの弦の表はどの程度正確か ?
正多角形で円周率 π=3.1416 を求める方法

ヒッパルコスの弦の表はどの程度正確か ?

ヒッパルコスの近似式がどの程度正確であるかを、Java プログラミングで確認してみることにします。i = 1, 2,.., 10° に対して、sin(i) と i×60/3438 の計算をして、その差を出力してみます。ソースは次のようになります。プログラミングを知らなければ単に信じてください。最後にはプログラミングは不必要となります。

MyCalculation.java


public class MyCalculation {
    static double  pi = 3.1416;

    public static void main(String[] args) {
        int i;
        double x,y,z;
        String str;
        for(i = 1; i<=10; i++){
            x = Math.sin(toRadian((double)i ));
            y = ((double) i)*60.0/3438.0;
            z = x-y;
            str = String.format("%.5f, %.5f, %.5f", x, y, z);
            System.out.println(str);
        }

    }

    public static double toRadian(double a){
        return a*pi/180;
    }

}

注意

上のプログラムは古い Java では動きません。String.format を使用しています。

すると出力結果は次のようになります。小数点以下 5 桁未満を丸めて表示し、わかりやすいように行番号も付けてあります。

	sin(α)	α×60/3438	差
1	0.01745	0.01745	0.00000
2	0.03490	0.03490	-0.00000
3	0.05234	0.05236	-0.00002
4	0.06976	0.06981	-0.00005
5	0.08716	0.08726	-0.00010
6	0.10453	0.10471	-0.00018
7	0.12187	0.12216	-0.00029
8	0.13917	0.13962	-0.00044
9	0.15643	0.15707	-0.00063
10	0.17365	0.17452	-0.00087

出力結果から判断すれば

sin(α/2)≒(α/2)×(60/3438)

の近似式は 0＜α/2＜7.5/2 の範囲で、小数点以下 4 桁までは一致することがわかります。

従って、ヒッパルコスが近似式を 0＜α≦7.5 に適用し、7.5° ごとの弦の表と、加法公式を使用すれば、ほぼ小数点以下 4 桁までの数値が得られることになります。

大きさと距離のトゥーマーによる第 1 巻の復元の計算ではヒッパルコスは chord(10) を近似式を使用して求めていますから、必ずしも厳密にしていたわけではないようです。

最初に次を示します。

補題 1.1

C = π/180 とする (π は円周率)。0°＜β＜α＜90° であれば 0 ＜ Cβ - sin(β) ＜ Cα - sin(α)

注意

この補題は明らかです。中心角 2α に対応する円弧の長さが 2Cα で、弦の長さが 2 sin(α) となります。 2Cα の方が 2 sin(α) より常に大きく、その差は中心角が大きくなれば大きくなるからです。
角度だと、ややこしくなります。ラジアンで表示すれば
0 ＜ β - sin(β) ＜ α - sin(α)
この補題はグラフを使用しても明らかです。y = x - sin(x) のグラフは y = x のグラフから y = sin(x) のグラフを引くことによって描け、これは単調増加です。以下はグラフを使用しない証明です。(面倒なので詳しく説明していない個所があります。不足分は補ってください。)

証明： 三角法の歴史 / ギリシャの数学のアリスタルコスの不等式から

＜

そこで

= b,

= a

と置けば

b＜a, β = b sin β, α = a sin α

このとき

	(Cα - sin α) - (Cβ - sin β)
＝	(Ca sin α - sin α) - (Cb sin β - sin β)
＝	(sin α) (Ca - 1) - (sin β) (Cb - 1)
＞	(sin β) (Ca - 1) - (sin β) (Cb - 1)
＝	(sin β) C(a - b) ＞ 0

これが証明したいことでした。証明終

さて 0≦θ≦7.5°/2 とする。この時、円周率を π とし、π=3.1416 とすれば

0 ＜

θ - sin θ ＜

(7.5/2) - sin(7.5/2) = 0.000046...

従って

定理 1.2

0 ≦ θ ≦ 7.5°/2 であれば、ヒッパルコスの近似式

sin θ

≒

×θ ≒

θ×

は小数点以下 4 桁まで有効である。

注意

ここで与えた近似式は θ の大きさを十分小さくすれば、もう少し高い精度で成立することがわかります。これに関しては、定理 2.2.2 をご覧ください。

ヒッパルコスの 7.5° ごとの弦の表は 60 進数で表示され、半径が 3438/60 (= 57.3) の場合の弦の長さですから、恐らく各弦の長さは小数点以下 1 桁の 60 進数で表示されていたはずです。ヒッパルコスの半径の選び方がややこしいのですが、もしも半径を 1 に選べば、おおよそ 2 桁の 60 進数で表示したと考えることができます。

仮にヒッパルコスの考えを使い、10 進表示で弦の表を作れば、小数点以下が 4 桁まで有効となり、考え方としては非常に優れていたものと思われます。

但し、これ以上の精度は恐らくありえない。それは

7.5° ごとの表 + 近似式 + 加法公式

を組み合わせて使用するため、桁数が多いと現実的でなくなるためです。

正多角形で円周率 π=3.1416 を求める方法

アルキメデスは正 96 角形を利用して、円周率を決めたと言われています。ここでは最初にアルキメデスの方法に三角法を援用すればどうなるかを考えます。正 n 角形の各辺を 2 等分すれば、正 2n 角形ができます。正 6 角形を 4 回 2 等分すれば正 96 角形ができます。

ヒッパルコスは 7.5° の弦を決定するのに、cos の倍角公式 (半角公式というべきか)

(sin α)² = (1 - cos(2α))/2

を使用して、60° (正六角形) から開始し、順に 30°, 15°, 7.5° に至るまで計算したと思われています。 (以上は「カッツ、数学の歴史」)

7.5° は正 6 角形を 3 回 2 等分して得られる正 48 角形の一辺の中心角です。これを更に何度も 2 等分して得られる正多角形で円周率を決める方法を考えましょう。

π = 3.1416 を求めるには正 768 角形を考えればよい -- 計算機チェック

ヒッパルコスの近似式は円周率が 3.1416 であることに依存しています。アルキメデスは円周率を内接する正多角形と外接する正多角形から求めたそうですが、ヒッパルコスはその上に三角法が使用できます。但しヒッパルコスは弦しか使用できないので、記号がややこしくなるはずです。ここでは面倒なので tan も使用して、どのようにして π = 3.1416 を求めることができるか調べてみましょう。まず、先に補題が必要となります。

補題 2.1.1

0＜2α＜90° のとき t = sin(α), A = sin(2α) とすれば
0＜2α＜90° のとき t = tan(α), A = tan(2α) とすれば

注意

常識的には cos の倍角公式を使用して、半分の角度の cos を決定してから sin, tan を決定する方が単純ですが、このようにすると誤差が拡大する可能性があります。そこで、ここでは sin, tan の倍角公式に相当するものを導いています。上の式からわかるように sin の場合は 2 度平方根を取る必要があり、tan の場合には桁数の大きな数で割る必要があります。古代の割り算はとても大変な計算となるようなので、いずれの公式も実質的な意味を持ちません。しかし、ここでは以上の式を使用します。

証明： 1. に関して

α が 0 に近づく時、A も 0 に近づき、この時 t も 0 に近づかないといけないので

2. に関して

α が 0 に近づく時、A も 0 に近づき、この時 t も 0 に近づかないといけないので

これが証明したいことでした。証明終

n を自然数とし、 θ = 360°/2n とすれば、半径 1 の円に内接する正 n 角形の 1 辺の長さは 2 sin(θ) で与えられ、外接する正 n 角形の 1 辺の長さは 2 tan(θ) で与えられる。

そこで

n = 6, 　2・6, 　2²・6, 　2³・6,....

に対して、半径 1 の円に内接する正 n 角形、外接する正 n 角形を構成し、その各々の辺の長さの総和の列を出力するような Java プログラムを作成すると、例えば次のようになります。

MyCircumference.java

public class MyCircumference {

  public static void main(String[] args) {
    int i;
    int n = 6;
    double x = 0.5;     // 内接する正 n 角形の 1 辺の長さの半分
                        // x = sin(2*pi/(2*m))
    double L = n*2*x;   // 内接する正 n 角形の周囲の長さ
    double y = 1/Math.sqrt(3.0);  
                        // 外接する正 n 角形の 1 辺の長さの半分
                        // y = tan(2*pi/(2*n))
    double M = n*2*y;   // 外接する正 n 角形の周囲の長さ
    String str;
    
    for(i = 1; i<=10; i++){
      str = String.format("%d, %d, %f, %f", i, n, L, M);
      System.out.println(str);
      n = n*2;
      x = half_sin(x);
      y = half_tan(y);
      L = n*2*x;
      M = n*2*y;
    }
  }
  
  // sin(2α) が与えられたときに sin(α) を返す関数。
  static double half_sin(double x){
    double y;
    y = Math.sqrt((1.0- Math.sqrt(1.0-x*x))/2.0);
    return y;
  }
  // tan(2α) が与えられたときに tan(α) を返す関数
  static double half_tan(double x){
    double y;
    y = (-1.0 + Math.sqrt(1+ x*x))/x;
    return y;
  }
}

結果は次のようになります。

	n	内接正多角形	外接正多角形
1	6	6.000000	6.928203
2	12	6.211657	6.430781
3	24	6.265257	6.319320
4	48	6.278700	6.292172
5	96	6.282064	6.285429
6	192	6.282905	6.283746
7	384	6.283115	6.283325
8	768	6.283168	6.283220
9	1536	6.283181	6.283194
10	3072	6.283184	6.283187

従って、円周率 3.1416 が求められるのは 8 番目の正 768 角形である (2×3.1416 = 6.2832)。ヒッパルコスがこの方法を取ったかどうか不明ですが、ともかくも可能な方法です。

アルキメデスは円周率を計算するために正 96 角形を使用しているそうですが、ヒッパルコスが類似の方法を採用すれば、三角法を利用したはずですから、ここで述べたような計算に類似していた可能性があります。

正弦の値を 10 進数で小数点以下 5 桁まで正多角形の方法で求めようとすると、円周率を 3.14159 まで計算する必要がある。しかし、円周率 3.14159 を求めるには 10 番目の正 3072 角形まで求めないといけない (2×3.14159 = 6.28318)。

768 を手計算で求め、円周率を計算機で求める -- 必要な桁数の推測

前節では n_i = 3×2ⁱ によって n_i (i = 1, 2, 3,...) を定め、ここから正 n_i 角形から円周率を決めるには、n_i = 768 の場合まで、計算する必要があることを計算機で確かめました。この節では n_i = 768 を計算機を使用せずに確認することを目的とします。

n を自然数とし、 θ = 360°/2n とすれば、半径 1 の円に内接する正 n 角形の辺の長さの合計は 2n sin(θ) で与えられ、外接する正 n 角形の辺の長さの合計は 2n tan(θ) で与えられます。円周率 π を小数点以下 4 桁まで求めるためには、あらかじめどこまで n を計算するかわかっていないといけません。その理由は、辺の長さの合計は sin や tan の値を 2n 倍して求めるため、sin や tan の精度をどこまで上げて計算すればよいかがわからないからです。しかし、これは比較的簡単に示すことができます。先に次を示します。

補題 2.2.1

α が微小角の時に次が成立。 tan α - sin α≒(1/2)(πα/180)³

証明：

tan α - sin α	＝	sin α ((cos α)^-1 - 1)
	＝	sin α ((1 - (sin α)²)^-1/2 - 1)
	≒	sin α ((1 + (1/2) (sin α)²) - 1)　　(∵ sin α は微小量)
	≒	(1/2) (sin α)³
	≒	(1/2) (πα/180)³　　(∵ α は微小角)

証明終

定理 2.2.2

正 n 角形の一辺の中心角を 2θとする。すなわち θ = 180/n とする。このとき、おおよそ 0＜π - n sin θ＜(1/2) π³/n² が成立する。但し右端の π は大きさを評価するためのもので π≒3.14 程度の近似でよい。

注意

定理の式を書き直すと

0＜(π/n) - sin θ	＜(1/2)×(π/n)³
0＜(π/n) - sin(180/n)	＜(1/2)×(π/n)³

『(180/n)度 = (180/n)×(π/180) ラジアン』なので、角度をラジアン表示して π/n = x とおくと、上の式は次の式と同等になります。

0＜x - sin(x)＜(1/2)x³

従って、これはテーラー展開

sin(x) = x - (1/3!) x³ + (1/5!) x⁵ - ...

から得られる式 x - sin(x)≒(1/3!)x³ の結果です。

証明： 内接正 n 角形と外接正 n 角形の周囲の長さから、次が成立していた。

n sin θ＜π＜n tan θ

よって

0＜π - n sin θ＜n(tan θ - sin θ)

補題 2.2.1 を適用すれば、おおよそ次の不等式が成立する。

０＜π - n sin θ＜n (1/2) (πθ/180)³

ここで

n (1/2)(πθ/180)³ = n(1/2)((π/180)×(180/n))³ = (1/2)×π³/n²

以上により、定理を得る。証明終

この定理により、円周率の決定方法が基本的に変化してしまいます。確かに、正多角形から円周率を求めることにはなりますが、 sin の計算のみでよいことになりました。(内接正多角形のみでよい。)

この定理を使用せずに正多角形から円周率 π=3.1416 を決めようとすると、 tan の評価で精度が落ちてしまいます。この節では最終的に π=3.1416 を決めるのに 60 進数で小数点以下 5 桁の計算でよいことを見ますが、定理を使用しないと、60 進数で小数点以下 6 桁の計算が必要となり、計算がとてもわずらわしいことになる。このような点から、円周率 π=3.1416 を決める際に、この定理は知られていたであろうと思われます。

定理を適用すれば、π を小数点以下 4 桁まで求めるには、おおよそ

(1/2)×π³/n²	＜	(1/2)×10^-4
n²	＞	π³×10⁴

とすればよい。π = 3.14 の概算値を使用すれば、おおよそ

n＞556

とすればよい。

注意

以上の評価は「おおよそ」であって厳密ではありません。

0＜π - n sin θ＜(1/2)×π³/n²＜0.5×10^-4

であっても π と n sin θ の小数点以下 4 桁目が一致するとはいえないからです。実際は π = 3.14159 ですから、ほぼ上の条件で決まります。
ヒッパルコスもプトレマイオスも厳密には円周率を小数点以下 2 桁の 60 進数 3;8,30 としています。そこで小数点以下 2 桁の 60 進数として求めることを考えれば、条件はおおよそ

0＜π - n sin θ＜(1/2)×π³/n²＜0.5×60^-2

となり、n＞333 となります。プトレマイオスが弦の数表を構成した時に考察した正 480 角形 (360/480 = 3/4) はこの条件を満たしており、円周率が確定する候補となります。誤差を計算してみると

(1/2)×3.14³/480²≒6.71..×10^-5≒0.24..×60^-2≒14.5..×60^-3

となり円周率が 60 進数の小数点以下 2 桁が確定する可能性があります。実際の円周率が 3;8,29,44,0 であることを知っていればほぼ確定することがわかる。しかし運悪く『3;8,29 もしくは 3;8,30』の結論しかない可能性もある。

i = 1, 2, ... に対して

n_i = 3 × 2ⁱ, θ_i = 360 / (2n_i) (= 180/(3×2ⁱ))

と置くと、n₇ = 384, n₈ = 768 なので、定理 2.2.2 により、正 n₈ 角形で評価すればよい。このとき定理 2.2.2 で π=3.14 の概算値を適用すれば

0＜π - n sin θ＜(1/2) π³/n²≒(1/2) 3.14³/768²＝2.6..×10^-5

従って、n₈ sin(θ₈) の正確な計算のみで π の値がおおよそ小数点以下 4 桁の精度で確定するはずである。

中間の角度の sin を求める必要がないので、cos の倍角公式から θ_i (i = 2,3,.., 7) を決めた後で cos θ₇ から sin θ₈ を求めるのが最も簡単な方法となります。

どこまでの精度で決定すればよいか不明なので、仮に小数点以下 d 桁まで cos θ₇ を漸化式

cos θ_i = ((cos θ_i-1 + 1)/2)^1/2 (i≦7)

によって決定し, sin θ₈ を

sin θ₈ = ((1 - cos θ₇)/2)^1/2

によって決定しよう。

次の方針で Java プログラムを作成する。実際の計算は 10 進法を使用するが、数値の記録は 60 進法でする。記録する時の 60 進数の精度を小数点以下 d にし、これはプログラムで可変となるようにし、出力を見ながら調節する。

計算には切り上げ・切り捨ての「はさみうち」と「丸め」(四捨五入) の 2 種類の計算をする。また計算を開始するのは θ₄ = 7.5°/2 とし、cos(7.5°/2) が正確に決定できているとする。 (正 48 角形が正確に確定しているとする。) 計算をするのは次の値である。

cos θ₅, cos θ₆, cos θ₇, sin θ₈

ほとんど電卓の計算でできる計算であるが、桁を調節できる点がプログラムの利点である。

出力結果を整理したものとプログラムを掲載する。以下の表で、円周率の 60 進表示の 1 桁目の 188 が妙に見えるが、これは 3×60 + 8 のこと (つまり 3;8 のこと) で、プログラムで少し手を抜いたためこのような表示になった。精度が 4 桁では円周率を求めることは到底無理であることがわかる。精度が 5 桁では下限・上限の「はさみうち」方式では 3;8,30 も 3.1416 も決まっているとは言えないが、「丸め」では 3;8,30 も 3.1416 も決まっている。誤差の 2.6..×10^-5 (= 5.6..×60^-3) が加わってもちゃんと評価できている。

精度 4 桁
	下限	丸め	上限
cos(θ₅)	59 58 4 20	59 58 4 21	59 58 4 22
cos(θ₆)	59 59 31 4	59 59 31 5	59 59 31 6
cos(θ₇)	59 59 52 45	59 59 52 46	59 59 52 47
sin(θ₈)	0 14 42 50	0 14 43 51	0 14 44 53
円周率(60進)	188 20 15 59	188 33 16 48	188 46 30 24
円周率(10進)	3.13896288	3.14257777	3.14625185

精度 5 桁
	下限	丸め	上限
cos(θ₅)	59 58 4 21 3	59 58 4 21 3	59 58 4 21 4
cos(θ₆)	59 59 31 5 8	59 59 31 5 9	59 59 31 5 10
cos(θ₇)	59 59 52 46 16	59 59 52 46 17	59 59 52 46 18
sin(θ₈)	0 14 43 32 58	0 14 43 34 0	0 14 43 35 2
円周率(60進)	188 29 25 58 23	188 29 39 12 0	188 29 52 25 36
円周率(10進)	3.14150913	3.14157037	3.14163160

Java プログラム

public class Pi {
  static int d = 5; // 桁数、これを調節する
  public static void main(String[] args) {
    int i,j;
    int[] u = new int[d];
    int[][] lcos = new int[3][d]; // cos の下限
    int[][] rcos = new int[3][d]; // cos の丸め
    int[][] ucos = new int[3][d]; // cos の上限
    int[] lsin = new int[d];      // sin の下限
    int[] rsin = new int[d];      // sin の丸め
    int[] usin = new int[d];      // sin の上限
    int[] lpi = new int[d];       // pi の下限
    int[] rpi = new int[d];       // pi の丸め
    int[] upi = new int[d];       // pi の上限
    
    //---------cos 下限-----------
    myprint("cos 下限\n");
    u = lower_cos(7.5/2.0);
    for(i=0; i<3; i++){
      u = lower_half_cos(u);
      for(j=0; j<d; j++){
        lcos[i][j] = u[j];
      }
      myprint(u, "\n");
    }
    //----------cos 丸め-----------
    myprint("cos 丸め\n");
    u = round_cos(7.5/2.0);
    for(i=0; i<3; i++){
        u = round_half_cos(u);
        for(j=0; j<d; j++){
            rcos[i][j] = u[j];
        }
        myprint(u, "\n");
    }
    //----------cos 上限-----------
    myprint("cos 上限\n");
    u = upper_cos(7.5/2.0);
    for(i=0; i<3; i++){
      u = upper_half_cos(u);
      for(j=0; j<d; j++){
        ucos[i][j] = u[j];
      }
      myprint(u, "\n");
    }

// ------ π 下限 ------------
    // cos は上限を使用
    for(i=0; i<d; i++){
        u[i] = ucos[2][i];
    }
    lsin = lower_sin(u);
    lpi = toSexa(768.0*toDecimal(lsin));
    
    // ----- π 丸め-------------
    // cos は丸めを使用
    for(i=0; i<d; i++){
        u[i] = rcos[2][i];
    }
    rsin = round_sin(u);
    rpi = toRoundSexa(768.0*toDecimal(rsin));

//------ π 上限-------------
    // cos は下限を使用
    for(i=0; i<d; i++){
        u[i] = lcos[2][i];
    }
    
    usin = upper_sin(u);
    upi = myceil(toSexa(768.0*toDecimal(usin)));
    
    // -----------------------
    myprint("\t\t下限  \t\t丸め\t\t 上限\n");
    myprint("sin 　　");
    myprint(lsin, " ");
    myprint(rsin, " ");
    myprint(usin, "\n");
    myprint("円周率　");
    myprint(lpi, " ");
    myprint(rpi, " ");
    myprint(upi, "\n");
    System.out.print("　　　　");
    System.out.print(toDecimal(lpi) + " ");
    System.out.print(toDecimal(rpi) + " ");
    System.out.println(toDecimal(upi));
  }

// sin θ = sqrt((1 - cos(2θ))/2)
  public static int[] lower_sin(int[] u){
      double w, z;
      w = toDecimal(u);
      z = Math.sqrt((1.0 - w)/2.0);
      return toSexa(z); 
  }
  public static int[] upper_sin(int[] u){
      double w, z;
      w = toDecimal(u);
      z = Math.sqrt((1.0 - w)/2.0);
      return myceil(toSexa(z)); 
  }
  public static int[] round_sin(int[] u){
      double w, z;
      w = toDecimal(u);
      z = Math.sqrt((1.0 - w)/2.0);
      return toRoundSexa(z); 
  }
  
  // 角度をラジアンに変換
  public static double toRad(double u){
    return ((Math.PI)*u/180.0);
  }
  // 10 進数を 60 進数 (配列) に変換。切り捨て
  public static int[] toSexa(double u){
    int[] x = new int[d];
    double v = u, w;
    for(int i=0; i<d; i++){
      v = 60.0*v;
      w = Math.floor(v);
      x[i]=(int)w;
      v = v - w;
    }
    return x;
  }
  //  10 進数を 60 進数 (配列) に変換, 丸め
  public static int[] toRoundSexa(double u){
        int[] x = new int[d];
        double v = u, w;
        for(int i=0; i<d; i++){
          v = 60.0*v;
          w = Math.floor(v);
          x[i]=(int)w;
          v = v - w;
        }
        v = 60.0*v;
        w = Math.floor(v);
        if (w < 30.0 ){
            return x;
        }
        return myceil(x);     
  }
  // cos(u°) を 60 進数 (配列) に変換、切り捨て
  public static int[] lower_cos(double u){
    int[] x = new int[d];
    x = toSexa(Math.cos(toRad(u)));
    return x;
  }
  // cos(u°) を 60 進数 (配列) に変換、切り上げ
  public static int[] upper_cos(double u){
    int[] x;
    x = toSexa(Math.cos(toRad(u)));
    x = myceil(x);
    return x;
  }
  // cos(u°) を 60 進数 (配列) に変換、丸め
  public static int[] round_cos(double u){
      int[] x;
      x = toRoundSexa(Math.cos(toRad(u)));
      return x;
  }
  //60 進数 (配列) を切り上げ
  public static int[] myceil(int[] u){
    int i;
    int[] x = new int[d];
    for(i = 0; i<d; i++){
        x[i] = u[i];
    }
    x[d-1]++;
    i = d-1;
    while(i>0){
      if(x[i]<60) return x;
      x[i]=0;
      x[i-1]++;
      i--;
    }
    return x;   
  }
  // 60 進数 (配列) を 10 進数へ
  public static double toDecimal(int[] u){
    double v = 0.0, t = 60.0;
    // 60 進数 (配列) を 10 進数に変換
    for(int i = 0; i<d; i++){
      v = v + ((double)u[i])/t;
      t = t*60.0;
    }
    return v;
  }
  // cos θ = sqrt((cos 2θ + 1)/2)
  public static int[] lower_half_cos(int[] u){
    double v;
    v = toDecimal(u);
    return toSexa(Math.sqrt((v + 1.0)/2.0));
  }
  public static int[] upper_half_cos(int[] u){
    return myceil(lower_half_cos(u)); 
  }
  public static int[] round_half_cos(int[] u){
      double v;
      v = toDecimal(u);
      return toRoundSexa(Math.sqrt((v + 1.0)/2.0));
  }
  
  
  public static void myprint(int[] u, String str){
    for(int i=0; i<d; i++){
      System.out.print(u[i] + " ");
    }
    System.out.print(str);
  }
  public static void myprint(String str){
      System.out.print(str);
  }
}

計算機の出力結果から cos(θ_i) の計算にはあまり誤差が出ていないことが分かるが、 sin(θ) の計算にはかなり誤差があることが分かる。次節でこれを評価してみよう。

誤差の評価

この節の内容は前節の誤差の理由を考えるもので、本当は計算の下準備で必要があったことです。一応他の節とは独立なので、読み飛ばしても前後関係が不明になることはない。

前節で計算機を使用して見いだした誤差の変動を証明してみることにしましょう。

以下の議論は微分の概念があればただちにわかることです。簡単にその理由を説明しておきます。関数 y = f(x) が与えられた時、x = a における y の値を決定することを考えます。今 x の値に誤差 ε がある時、y の値に誤差 ε' が生じるとします。この時

ε'/ε≒(dy/dx)_{x = a}

の関係が成り立つ。これは微分係数の意味から明らかなことです。つまり ε'/ε は y = f(x) の点 a における接線の傾きです。

古代ギリシャには微分の概念はありませんが、恐らく次は知られていた：

n が半整数 (= 整数の 1/2) でも負の数でも、 t が微小量であれば (1 + t)ⁿ≒1 + nt

以下の議論では以上を前提にします。

cos の倍角公式を利用して cos θ から cos(θ/2) を求める時に、 θ が余り小さくなければ誤差に余り変動がないことが想定される。そこで、θ が微小角であれば cos θ から cos(θ/2) を求める際に誤差が拡大しないことを示そう。

補題 2.3.1

角 θ に対して cos(θ/2) は cos(θ/2) = ((cos θ + 1)/2)^1/2 によって求められる。さて θ が微小角の時に cos θ の近似値 X が与えられ、誤差が ε とする。そのとき cos(θ/2) = ((X + 1)/2)^1/2 によって cos(θ/2) を定めた時の誤差 ε' は ε/4 程度である。

注意

少しあいまいな表現であるが、ここで述べている誤差とは |ε|＜10^-4 のような条件のことで、この場合には大体 |ε'|＜10^-4/4 となることを述べている。ここでは誤差の大きさ (オーダー) を問題としている。

証明： θ は微小角であるから

cos θ = 1 - A

と書ける。但し A は微小量である。 cos θ の近似値 X を

X = 1 - A + ε

とする。その時 cos(θ/2) の近似値は次で与えられる。

cos(θ/2)	≒ ((X+1)/2)^1/2
	＝ ((1 - A + ε + 1)/2)^1/2
	＝ ((1 - A/2) + ε/2)^1/2
	＝ (1 - A/2)^1/2 (1 + (ε/2)/(1 - A/2))^1/2
	≒ (1 - A/2)^1/2 (1 + (ε/4)/(1 - A/2))
	＝ (1 - A/2)^1/2 + (ε/4)/(1 - A/2)^1/2

cos(θ/2) の正確な値は (1 - A/2)^1/2 であるから、誤差 ε' はおおよそ

ε' = （ε/4)/(1 - A/2)^1/2

である。(1 - A/2)^1/2 は 1 に近いので ε' はおおよそ ε/4 に等しい。証明終

次に θ が微小角の時に cos θ から sin θ を求める時に誤差が拡大することを示したい。そのためには次が必要となるが、これは明らかである。

補題 2.3.2

θ が微小角の時に

sin θ	≒ πθ/180
cos θ	＝ (1 - (sin θ)²)^1/2
	≒ 1 - (sin θ)²/2
	≒ 1 - (πθ/180)² / 2

示したいことは次である。sin をこの方法では決定しなかったが、この方法で sin を決定しようとすると誤差が大きくなりすぎることが判明したからである。

補題 2.3.3

sin θ は sin θ = (1 - cos² θ)^1/2 によって求められる。θ が微小角の時 cos θ の近似値 X が与えられ、誤差が εであるとする。その時 sin θ = (1 - X²)^1/2 によって sin θ を定めた時の誤差 ε' はおおよそ ε' ≒ 180ε/(πθ) である。つまり cos から sin を求める時に誤差は 180/πθ 倍される。

証明： θ は微小角であるから

cos θ = 1 - A

と書ける。但し A は微小量である。cos θ の近似値 X を

X = 1 - A + ε

とする。このとき sin θ の近似値は次で与えられる。

sin θ	≒ (1 - X²)^1/2
	＝ (1 - (1 - A + ε)²)^1/2
	＝ (1 - (1 - A)² - 2ε(1 - A) - ε²)^1/2
	≒ ((2A - A²) - 2ε(1 - A))^1/2
	＝ (2A - A²)^1/2 (1 - 2ε(1 - A)/(2A - A²))^1/2
	≒ (2A - A²)^1/2 (1 - ε(1 - A)/(2A - A²))
	＝ (2A - A²)^1/2 - ε(1 - A)/(2A - A²)^1/2
	≒ (2A - A²)^1/2 - ε/(2A)^1/2

となる。A は 1 に比べて微小量であることに注意する。従って ε'≒ε/(2A)^1/2 である。ここで補題 2.3.2 により、

A ≒ (sin θ)²/2 ≒ (πθ/180)²/2

であるから

ε'≒ε/(2A)^1/2≒ε/(πθ/180)＝180ε/(πθ)

を得る。証明終

sin を決める時に cos の倍角公式を利用した。こちらの方が精度で優れていることを示そう。補題 2.3.3 と比べると、誤差がわずか 1/4 しか低減しないが、これによりかなり良い評価になる。

補題 2.3.4

sin θ は sin θ = ((1 - cos(2θ))/2)^1/2 によって求められる。θ が微小角の時 cos(2θ) の近似値 X が与えられ、誤差が ε であるとする。その時 sin θ = ((1 - X)/2)^1/2 によって sin θ を定めた時の誤差 ε' はおよそ ε' ≒ (ε/4) × 180/(πθ) である。

証明： cos(2θ) を

cos(2θ) = 1 - A

と書く。但し A = 2 (sin θ)² である。cos(2θ) の近似値 X を

X = 1 - A + ε

とする。このとき sin θ の近似値は次で与えられる。

sin θ	＝((1 - X)/2)^1/2
	＝((A - ε)/2)^1/2
	＝(A/2)^1/2(1 - ε/A)^1/2
	≒(A/2)^1/2(1 - ε/(2A))
	＝(A/2)^1/2 - ε/(2^3/2×A^1/2)

一方

A	＝ 2(sin θ)² ≒ 2 (πθ/180)²
A^1/2	≒ 2^1/2 (πθ/180)

よって

ε/(2^3/2×A^1/2)≒(ε/2^3/2)×(1/2^1/2)×(180/(πθ)) = (ε/4)×(180/(πθ))

を得る。以上により補題が示された。証明終

この節における結果は、実際に円周率を計算する前に必ずしたはずである。円周率の精度を上げようとして何度も失敗をしたはずで、そのあげく、この節におけるような評価をする必要があったはずである。これは前節の計算機による円周率 π = 3.1416 の決定の際にも必要であったことであり、誤差の評価なしに計算を遂行することは恐らく不可能である。 (定理 2.2.2 は誤差の評価をしている最中に気が付いた。)

どの正多角形から開始すればよいか

定理 2.2.2 を前提にすれば、正多角形から円周率を決める際にアルキメデスがしたように正六角形から開始する必要はない。定理 2.2.2 の直後に記述したことから、正 n 角形から円周率 π = 3.1416 を決めるには n＞556 であればよいだけである。

それではどの正 n 角形を考えればもっとも効率がよいであろうか ? この正 n 角形は定規とコンパスで作図可能でないといけない。そうでないと、一辺の長さを平方根で記述できないからである。そうすれば

n = 3, 4, 5, 15

から始め、これを 2 等分してできる正多角形しかない。

注意

正 n 角形が定規とコンパスで作図可能であるためには n が

n = 2^r p₁... p_s

の形でないといけない。但し p_i は互いに相異なる奇素数で 2^{2^t} + 1 の形をしている必要がある (フェルマー素数)。17, 257, 65537 もこのような形をした素数であるが、古代ギリシャでは知られていなかったと思う。

つまり n は次の形をしている。

n = 3×2ⁱ, 4×2ⁱ, 5×2ⁱ, 15×2ⁱ (i = 0, 1, 2, 3,... )

列挙してみると次のようになる。

3, 6, 12, 24, 48, 96, 192, 384, 768
4, 8, 16, 32, 64, 128, 256, 512, 1024
5, 10, 20, 40, 80, 160, 320, 640
15, 30, 60, 120, 240, 480, 960

従って、正 768 角形以外に、π = 3.1416 を求めることが容易である可能性のあるのが正 640 角形ですが、これは恐らく似たような計算になります。

60 進数で円周率が小数点以下 2 桁で確定するには n＞333 が候補になり、正 480 角形もその 1 つです。実際に計算してみることにしましょう。次を実行します。(この計算は基本的にはアルマゲストに記述されている chord(3/4) を計算する方法です。)

n_i = 60×2ⁱ, θ_i = 180/n_i (i = 0, 1, 2, 3) と置き

cos(θ₀) を決める
cos(θ_i+1) = ((cos(θ_i) + 1)/2)^1/2 (i = 1, 2) により、θ₁, θ₂ を決める
sin(θ₃) = ((1 - cos(θ₂))/2)^1/2 により sin(θ₃) を決める。
480×sin(θ₃) を出力する。 (注意：sin(θ₃) = sin(3/8) = (1/2)×chord(3/4))

プログラムは正 768 角形から円周率を決めたプログラムの定数を変更するだけの、微修正で十分です。結果を整理すると次のようになります。円周率(60進) と書いてある行は sin(θ₃)×480 の値で、これに誤差の最大である 0;0,0,30 を加えても 60 進数の小数点以下 2 桁に変化が見られないので、これで円周率は 60 進数の小数点以下 2 桁である 3;8,30 が確定したことになります。 10 進表示の場合も 3.1416 が決定されているかのように見えますが、誤差が大きいのでこれは確定していません。定理 2.2.2 はもう少し高い精度で成立する (定理 2.2.2 の注意を参照) ので、実際は 3.1416 が決定されると思いますが、古代では精度を上げることが少し無理だと思う。

精度 5 桁
	下限	丸め	上限
cos(θ₁)	59 58 45 58 55	59 58 45 58 56	59 58 45 58 57
cos(θ₂)	59 59 41 29 40	59 59 41 29 41	59 59 41 29 42
sin(θ₃)	0 23 33 41 50	0 23 33 42 28	0 23 33 43 7
円周率(60進)	188 29 34 40 0	188 29 39 44 0	188 29 44 56 1
円周率(10進)	3.14154938	3.14157283	3.14159691

微修正したプログラムも掲載しておきます。

public class Pi2 {
  static int d = 5; // 桁数、これを調節する
  public static void main(String[] args) {
    int i,j;
    int[] u = new int[d];
    int[][] lcos = new int[2][d]; // cos の下限
    int[][] rcos = new int[2][d]; // cos の丸め
    int[][] ucos = new int[2][d]; // cos の上限
    int[] lsin = new int[d];      // sin の下限
    int[] rsin = new int[d];      // sin の丸め
    int[] usin = new int[d];      // sin の上限
    int[] lpi = new int[d];       // pi の下限
    int[] rpi = new int[d];       // pi の丸め
    int[] upi = new int[d];       // pi の上限
    
    //---------cos 下限-----------
    myprint("cos 下限\n");
    u = lower_cos(3.0);
    for(i=0; i<2; i++){
      u = lower_half_cos(u);
      for(j=0; j<d; j++){
        lcos[i][j] = u[j];
      }
      myprint(u, "\n");
    }
    //----------cos 丸め-----------
    myprint("cos 丸め\n");
    u = round_cos(3.0);
    for(i=0; i<2; i++){
        u = round_half_cos(u);
        for(j=0; j<d; j++){
            rcos[i][j] = u[j];
        }
        myprint(u, "\n");
    }
    //----------cos 上限-----------
    myprint("cos 上限\n");
    u = upper_cos(3.0);
    for(i=0; i<2; i++){
      u = upper_half_cos(u);
      for(j=0; j<d; j++){
        ucos[i][j] = u[j];
      }
      myprint(u, "\n");
    }

// ------ π 下限 ------------
    // cos は上限を使用
    for(i=0; i<d; i++){
        u[i] = ucos[1][i];
    }
    lsin = lower_sin(u);
    lpi = toSexa(480.0*toDecimal(lsin));
    
    // ----- π 丸め-------------
    // cos は丸めを使用
    for(i=0; i<d; i++){
        u[i] = rcos[1][i];
    }
    rsin = round_sin(u);
    rpi = toRoundSexa(480.0*toDecimal(rsin));

//------ π 上限-------------
    // cos は下限を使用
    for(i=0; i<d; i++){
        u[i] = lcos[1][i];
    }
    
    usin = upper_sin(u);
    upi = myceil(toSexa(480.0*toDecimal(usin)));
    
    // -----------------------
    myprint("\t\t下限  \t\t丸め\t\t 上限\n");
    myprint("sin 　　");
    myprint(lsin, " ");
    myprint(rsin, " ");
    myprint(usin, "\n");
    myprint("円周率　");
    myprint(lpi, " ");
    myprint(rpi, " ");
    myprint(upi, "\n");
    System.out.print("　　　　");
    System.out.print(toDecimal(lpi) + " ");
    System.out.print(toDecimal(rpi) + " ");
    System.out.println(toDecimal(upi));
  }

以上の計算から、プトレマイオスは弦の数表の計算をする際に、円周率 π が 60 進表示で 3;8,30 であることを再確認できたことがわかります。

プトレマイオスの弦の数表との関連

前節でみたように、プトレマイオスが弦の数表を作成した時に円周率 π が 60 進表示で 3;8,30 であることを再確認できたことがわかりました。これは単に、ついでにしたことでしょうか ? 恐らくこれは必要性に駆られたからだと思います。定理 2.2.2 の式を書き直します。

正 n 角形の一辺の中心角を 2 θとする。すなわり θ = 180/n とする。このとき、おおよそ 0＜π/n - sin θ＜(1/2) (π/n)³ が成立する。但し右端の π は大きさを評価するためのもので π≒3.14 程度の近似でよい。

この式を使用すれば、プトレマイオスが論証した chord(1°) = 0;1,2,50 をチェックできます。単に上の式で n = 360 で適用すればよいだけです。そうすると

0＜π/360 - sin(1/2)	＜(1/2) (π/360)³
0＜π/180 - chord(1)	＜ (π/360)³≒0;0,0,0,8,.. (60進)

従って chord(1) は π/180 にほぼ等しくなります。 π/60 は π = 3;8,30 の小数点を左にずらして得られる数値 0;3,8,30 に一致し、従って π/180 はそれを 3 で割った数 0;1,2,50 に一致します。しかも 3 で割っていますから、誤差が最大 0;0,0,0,30 であったものが 0;0,0,10 に減ってしまいました。つまり

π/180 = 0;1,2,50,0 ± 0;0,0,0,10

だからプトレマイオスは完全に chord(1) = 0;1,2,50 をチェックできたのです。これは極めて重大です。つまりプトレマイオスは π=3;8,30 の数値をうのみにして、このチェックをしたのではなく、 π=3;8,30 であることを再確認した上で、このチェックをしているのです。

なお、これと同じ方法で chord(3/2) = 0;1,34,15 及び chord(3/4) = 0;0,47,7 をチェックできます。プトレマイオスの弦の数表の訳注で、アルマゲストに chord(3/4) = 0;0,47,8 と書いてあるのは、誰かが写しまちがえたのに相違ないと書きましたが、これが裏付けられることになりました。プトレマイオスには検算方法があるので、間違えるはずがないのです。ここまでをまとめると

結論

プトレマイオスは弦の数表を作成するときに、円周率 π = 3;8,30 (60進) の再確認をして
これを使用して、弦の数表の検算に利用した。

このように結論すると色々と憶測できることがあります。

ひょっとすると π=3.1416 (ないしは 60 進数で π=3;8,30) を最初に決定した人はヒッパルコスでないかもしれません。しかし、この値はヒッパルコスにとって致命的な値で、これが正確でないと彼の理論はすべて無益となります。だから、これを最初に決定したのはヒッパルコスでないかもしれませんが、その場合でも必ず再チェックをしたはずです。

プトレマイオスにとっても π=3;8,30 は弦の数表の検算をするのに極めて重要な役目をし、そのためプトレマイオスは π=3;8,30 の再チェックをしているのです。だからヒッパルコスも同様なことをしたのに相違ないのです。

ヒッパルコスが π=3;8,30 を決定ないしは再チェックするのに使用したのは正 768 角形かもしれません。

プトレマイオスも最初は π=3;8,30 を再チェックしようとして、ヒッパルコスよりも易しくする方法を考え、正 480 角形に到達したのかもしれません。その過程で、chord(1°) を決定できることに気がつき、弦の数表の作成に着手したのかもしれません。

もう一つ付け加えておくべきでしょう。それはプトレマイオスが定理 2.2.2 を知っていたのであれば、何故それを残さなかったのかという点です。定理 2.2.2 の式にせよ、誤差の評価式にせよ、これはおよその評価です。厳密な評価もありますが、大体は大きさ (オーダー) を見ているだけです。従って、これは厳密なギリシャの数学の言葉にはなりません。 (現代の解析学の入門書には登場する言葉です。) だから定理 2.2.2 からわかる事実であれ、あえて定理 2.2.2 を登場させずに、代替の論法を取ったのであろうと思われます。プトレマイオスは chord(1) を決定するときには、ギリシャ数学における厳密な論法を取っているのです。

ヒッパルコスは α が微小角の時、現代の言葉では sin(α) = πα/180 となる近似式を使用しているのです。これは高い精度で一致する式で、よしんば定理 2.2.2 に気がつかなくても、chord(1) を求めた時点で近似式を適用してみれば chord(1) の値が正確にわかることに気が付いたはずで、必ず定理 2.2.2 のような事実に気がつき、またこの理由を示すことができたはずです。表面にこれを登場させなかったのは、これが単に厳密なギリシャ数学でなかったからです。