弓一夫一付派人共四千打的急但是大刀小斧不合又多欠工才怯下了兵威。
各位同學:上面這一段半通不通的文字, 是我在十幾年前寫的。雖然文字不十分通順, 但寫的時候卻花費了一番心血, 因為其中暗含著一些玄機:你如果數一數文字中每一個字的筆畫, 你便會發現弓字三筆, 一字一筆, 夫字四筆, 等等。連貫起來, 便得到$$3.141592653589793238462643383279,$$ 這包括著$\pi$的整數部分和其小數部分的前三十位。
現在的數學家已經知道, $\pi$是一個無理數, 不能用有盡或循環小數表示。$\pi$的計算, 自古便是數學家極感興趣的問題。 例如公元前三世紀, 阿基米德曾算出 $$\frac{{221}}{{71}} \lt \pi \lt \frac{{22}}{7}$$ 即$3.1126760563\lt \pi \lt 3.1428571429$,
公元後三世紀三國的魏人劉徽計算出$\pi$的近似值為3.14159。1706年英人馬秦(John Machin)將$\pi$算到 100 位小數。 1974年利用電子計算機更將$\pi$算到小數一百萬位。時至今日每當新的電腦或新的編譯軟體(compiler)設計完成後, 常用它計算$\pi$來試測其速度和精確度。
今天的講演是想介紹給大家一個算$\pi$的方法和程式, 並利用這程式將$\pi$算到小數 1,000 位。 我們的方法基於下面的兩則公式和一條定理。
反正切的加法公式： $$\arctan x + \arctan x = \arctan \frac{{x + y}}{{1 - xy}}$$ 反正切級數展開公式 $$\arctan x = \sum\limits_{n = 1}^\infty {{{( - 1)}^{n + 1}}} \frac{{{x^{2n - 1}}}}{{2n - 1}}$$ 交錯級數定理:設{$a_n$}為遞減正項敘列,即$a_1\ge a_2\ge\dots\ge a_n\ge\dots\gt 0$
且$\lim a_n=0$, 則級數 $$\sum\limits_{n = 1}^\infty {{{( - 1)}^{n + 1}}{a_n}} $$ 必收斂。設其和為s , 則 $$|s - \sum\limits_{n = 1}^m {{{( - 1)}^{n + 1}}{a_n}} | \le {a_{m + 1}}$$ 反正切的加法公式見於高中數學課本。反正切級數展開公式和交錯級數定理見於徵積分課本。 如果我們能把$\pi$利用反正切函數表示出來, 那麼至少在理論點我們可以利用反正切級數展式的部分和將$\pi$算出, 而交錯級數定理可以用來估計計算的誤差。
利用反正切計算$\pi$值, 最直接的公式似乎是 $$\pi = 4\arctan 1 = 4\sum\limits_{n = 1}^\infty {\frac{{{{( - 1)}^{n + 1}}}}{{2n - 1}}} $$ 但如果我們利用這公式計算部分和至一萬項, 則得到的近似值為3.141693, 其精確度僅至小數點後的第三位。事實上, 若用交錯級數定理估計誤差, 取$n$項的部分和, 希望正確到小數點後的1,000位, 則必須 $$\frac{4}{{2n + 1}} \lt {10^{ - 1000}}$$ 即$n \gt 2 \times {10^{1000}}$, - 2 後加上一千個0。這數字比天文數字還大了許多, 絕非以有涯的生命可以對付得了的。
雖然這方法失敗了, 但從這方法我們可以得到一些啟示:若$\theta$相當小, 則$\arctan \theta$的計算並不那麼困難。 因此我們要尋覓的是用較小的有理數的反正切表示$\pi$的公式。
英國天文學家馬秦便找到了一個這樣的公式。令 $$x = \arctan \frac{1}{5}$$ 則由反正切的加法公式得 $$2x = \arctan \frac{{\frac{2}{5}}}{{1 - \frac{1}{{25}}}} = \arctan \frac{5}{{12}}$$ $$4x = \arctan \frac{{\frac{5}{6}}}{{1 - \frac{{25}}{{144}}}} = \arctan \frac{{120}}{{119}}$$ 因$\frac{120}{119}$和1很接近, 所以$4x$和$\arctan 1=\frac{\pi}4$也很接近。它們的差是 $$4x - \frac{\pi }{4} = \frac{{\frac{{120}}{{119}} - 1}}{{1 + \frac{{120}}{{119}}}} = \arctan \frac{1}{{239}}$$ 換言之, $$\frac{\pi }{4} = 4\arctan \frac{1}{5} - \arctan \frac{1}{{239}}$$ 這公式很好用, 因為第二項的級數收斂很快, 而第一項的級數各項的計算又很容易。馬秦用這公式把$\pi$算到了一百位, 他死後後人把他算出的$\pi$值刻到他的墓碑上。這墓碑也是人類進步的一個里程碑。
第一次用電腦算出$\pi$, 是用ENIAC機器作的。這機器花了七十小時的時間, 將$\pi$計算到2037位。計算的方法, 便是利用馬秦公式。
但馬秦公式並不是最好的工具, 後來人用電腦算$\pi$, 慣用的是1896年斯妥摩(Stormer)所發現的公式。 $$\pi = 24\arctan \frac{1}{8} + 8\arctan \frac{1}{{57}} + 4\arctan \frac{1}{{239}}$$ 但馬秦公式和斯妥摩公式算$\pi$的速度都不如下述的公式好。這公式是 $$\pi = 32\arctan \frac{1}{{10}} - 4\arctan \frac{1}{{239}} - 16\arctan \frac{1}{{515}}$$
這公式我最初在一本1922年出版的德文小冊子O.Th.Bürklen:Mathematische Formelsammlung, W.de Gruyter上看到。 小冊子稱這公式為邁瑟爾(Meisel)公式, 但對邁瑟爾是何許人並沒有介紹。後來又在W. L. Schaaf 著的The History, Nature and Computation of $\pi$, Stanford, 1967(中央書局出版有賴建業的中譯本:$\pi$的歷史、性質和計算) 看到這公式, 說它是由S. Klingenstierna 在 1730 年得到的。但我已經把這式叫慣了邁瑟爾公式, 以下仍稱它為邁瑟爾公式吧。
邁瑟爾公式的導法如下:由反正切的加法公式得 $$2\arctan \frac{1}{{10}} = \arctan \frac{{\frac{1}{5}}}{{1 - \frac{1}{{100}}}} = \arctan \frac{{20}}{{99}}$$ 因為$ \frac{20}{99}$和$\frac{1}{5}$很接近, 我們願意算它們的反正切的差 ： $$2\arctan \frac{1}{{10}} - \arctan \frac{1}{5} = \arctan \frac{{\frac{{20}}{{99}} - \frac{1}{5}}}{{1 + \frac{1}{5}.\frac{{20}}{{99}}}} = \arctan \frac{1}{{515}}$$ 故得 $$\arctan \frac{1}{5} = 2\arctan \frac{1}{{10}} - \arctan \frac{1}{{515}}$$ 將這結果代入馬泰公式, 便得邁瑟爾公式。
我們現在的目標是將$\pi$計算到1000位小數, 即包含整數部分的1001位有效數字。在馬秦、斯妥摩和邁瑟爾的三個公式中, 若僅計算到1000位小數, 恐怕還是以馬秦的公式為最快。但若計算更多位小數, 邁瑟爾的公式便會大大地超前了。 因此今天向諸位介紹的便是這方法：——先用反正切的級數展開計算 $\arctan \frac{1}{{10}},\arctan \frac{1}{{239}}\arctan \frac{1}{{515}}$, 再用邁瑟爾公式求得$\pi$的值。
這三個反正切級數該計算到多少項才能保證$\pi$的值。值正確到1,000位小數呢？注意三個級數都呈$\arctan \frac{1}{n}$的形式, 其中$n$是正整數。因 $$\arctan \frac{1}{n} = \sum\limits_{k = 0}^\infty {\frac{{{{( - 1)}^k}}}{{2k + 1}}{{\left( {\frac{1}{n}} \right)}^k}} $$ 滿足交錯級數定理的條件, 故若該級數的前$m$項的和作為$\arctan \frac{1}{n}$的值, 其誤差的絕對值一定不超過$\frac{1}{{2m + 1}}\frac{1}{{{n^m}}}$。所以要想這個和正確到小數點後$a$位, 則必須 $$\frac{1}{{2m + 1}}\frac{1}{{{n^m}}} \lt \frac{5}{{{{10}^{a + 1}}}}$$ 即 $$ - \log (2m + 1) - m\log n \lt \log 5 - (a + 1)$$ 或 $$\log (2m + 1) + m\log n \gt a + 1 - \log 5$$ 式中 log 表常用對數。找尋這個$m$的方法可略述如下:我們先找最小的$m$, 使 $$m\log n \gt a + 1 - \log 5$$ 只要用除法便可以找到這個$m$了, 但找到的$m$可能比需要的大。因此我們將$m$的值遞減, 試測不等式是否仍然正確, 直到不能再減為止。
其次我們要問, 到底我們要算到多少項呢？假定每個級數近似反正切時的誤差都不超過$\varepsilon$, 那麼邁瑟爾公式近似$\pi$的誤差便不會超過$(32 + 4 + 16)\varepsilon = 52\varepsilon $。 若計算每個反正切級數至$a$位小數能使$\pi$的值正確至1,000 位小數, 則須 $$52 \times 5 \times {10^{ - a - 1}} \lt 5 \times {10^{ - 1001}}$$
可以看出, 若取$a=1003$ 便可以達成這要求。但為了防止下文所述由四捨五入引進的另一種誤差, 我們取$a=1008$。 如是 1000 位小數的結果可以完全保證其正確。
我的計算程式是使用Pascal 語言撰寫, 在詮腦電腦上編譯執行的。電腦中裝有8087 晶片, 可以迅速且精確地處理實數。 但這個 Pascal 的編譯軟體所認識的實數, 僅有十五位有效數字, 而整數也必須在-32768和32767之間。 我們如何能在這種情況下計算$\pi$至一千位小數呢?答案是我們利用多重精密度的算術。
諸位從小學起所學的算術都是十進位的, 在中學對其他進位的算術也略有涉獵。等到學習計算機科學, 則著重於2進位, 8進位和16進位的算術。在我們算$\pi$的程式中, 我們採用萬進位的四則運算, 每種運算都必須寫入程式之中。
每位萬進位的數字, 都可以用四位十進位的數字表示, 所以用一個 Pascal 整數表示它便足夠了, 我們稱一個具有1008位小數的實數為長實數。長實數都可以用引數(index)在0和252間的整數敍列(array)來表示, 第0項表示其整數部分, 以後每項依遞降次序表示其各位的萬進位數字。計算長實數的四則運算如下:因為兩位萬進位數字的和與差, 都在-9999 和 19998 之間, 我們仍可以用整數來表示它, 然後再用進位和借位的方法得到結果。至於乘法就沒有這樣容易了。 兩位萬進位數字的乘積很容易超出 Pascal 語言中整數的範圈。因此我們可以藉助於實數的運算, 再在進位後化回整數來求得乘積。 在我們的$\pi$的計算中, 只用到一個一位的萬進位數字和一個長實數的乘法, 和一個長實數被一個一位的萬進位數字除的除法, 所以乘式並不難寫, 當然我們也可以用組合語言直接寫多重精確度的四則運算。但組合語言依附於機種, 不適合今天我們公開講演, 所以不予採用。
在我們的程式中, 長實數有 252 位萬進位小數, 也就是 1008 位十進位小數。 對第 1009 位小數, 我們不採四捨五入的辦法, 而是一概捨棄。因此計算級數時每加一項都可能在最後一位上造成1的誤差。若反覆計算一千項, 那麼最後的三位數字便不甚可靠。 好在我們保留了八位數字作為謹慎的預防, 這才有1000位小數正確的把握。實際上八位十進位的數字只不過是兩位萬進位的數字而已。
我們把所寫的 Pascal 程式及其計算結果附在後面。這程式全部執行的時間不到十五分鐘。以下讓我們略述程式的內容。 先看宣告部分: 我們用 Ten Thousand 表示一萬, 是由於要用萬進位的緣故。 Qd Places = 252, 我們要計算到 252 位萬進位數字。 長實數的型(type) 為 long, 這表示以 0 至 252 為指標的整數敍列, 第 0 元為整數部分, 第 $k$ 元為第 $k$ 位萬進位小數 。
其次我們依次看程式中的子程式, 諸位可參閱所附程式的報表。Write Where 模組決定在螢幕上, 列表機上或檔案中輸出。 Outp 將一個萬進位的長實數寫成通常小數的形式。注意我們每行 15 個萬進位數字, 即 60 個十進位數字。
Add, Sub, Mul 和 Divl 諸模組作萬進位小數的四則計算, 我們只作下文所需的運算, 因此沒有計算兩個長實數的積和商。 這種計算可以留給諸位作為程式設計的線習。
函數 log 求正整數$n$的常用對數, 這是為下一個函數 TermNo作準備的, 而TermNo的目的是求$\arctan\frac1n$正確至$a$(=1008)位小數需要在其級數展式中計算的項數。 下一個子程式便是$\arctan\frac1n$的計算了。最後的主程式將以上的資料連慣起來, 利用邁瑟爾公式, 求出$\pi$的值至 1000 位。作了這麼多的計算, 若只印出$\pi$實在有點不甘心。所以我們也把 $\arctan \frac{1}{10} \arctan \frac{1}{239}和\arctan \frac{1}{515}$的值也印了出來, 全都印到 1000 位小數。
THE PROGRAM
program Meise1 ( input,tx)；
{ Meise1's algorithm to compute pi:
pi = 8 arctan 1/10 - arctan 1/239 -4 arctan 1/515
Programmed by Ju-kwei Wang. }
const TenThousand=10000；
QdPlaces-252;
type long=array[0..QdPlaces] of integer；
var tx:text;
one, refsult,summand, pi: long;
i:iriteger;
procedure WriteWhere(var t:text);
{This procedure decides where to place the output}
var ans: char;
Filespec: string[14];
begin
write('P(rinter, V(ideo or F(i]e output?');
repeat read(kbd, ans) unti1 ans in ['p','v','f','P','V','F'];
case ans of
'P','P': assign(t, 'LST:');
'v','V': begin
assign(t,'CON:');
clrscr;
writeln
end;
'f','F'： begin
write('Filespec:');
readln(Fi1espec)；
assign(t,Filespec)
end
end; {case}
rewrite(t)
end;
procedure outp(a:long);
var i:integer;
begin
write(tx,a[0]:4,'.');
for i:=l to QdPlaces-2 do begin
if a[i]>999 then write(tx,a[i]:4,' ')
else if a[i]>99 then write(tx,'0',a[i]:3,' ')
else if a[i]>9 then write(tx,'00',a[i]:2,' ')
else write(tx,'000',a[i]:1,' ');
if (i+1) mod 15=0 then writeln(tx)
end;
writeln(tx)；
writeln(tx)
end;
procedure add(a,b:long; var sum:long)；
var i:integer；
begin
for i:=0 to QdPlaces do sum[i]:=a[i]+b[i]；
for i:=QdPlaces downto 1 do if sum[i]>=TenThousand then begin
sum[i]:=sum[i]-TenThousand;
sum[i-1]：=sum[i-1]+1
end
end;
procedure sub(a,b:long；var diff:long);
var i :integer;
begin
for i:=0 to QdPlaces do diff[i]:=a[i]-b[i];
for i:=QdPiaces downto 1 do if diff[ij <0 then begin
diff[i]:=diff[i]+TenThousand;
diff[i-1]:=diff[i-1]-1
end
end;
procedure mul(a:long; n:integer; var prod:long);
var i,advance:integer;
x:real;
begin
advance:=0;
for i:=QdPlaces downto 1 do begin
x:=1.0*a[i]*n+advance;
advance:=trune(x/TenThousand);
prod[i]:=trunc(x-1.0*advance*TenThousand)
end;
prod[0]:=a[0]+advance
end;
procedure divl(a:long; n:integer; var qo:long)；
var i:integer;
x:real;
begin
x:=a[0];
for i:=0 to QdPlaces-1 do begin
qot[i]:=trunc(x/n);
x:=(x-1*0*qot[i]*n)*TenThousand+a[i+1]
end;
qot[QdPlaces]:=trunc(x/n)
end;
function log(n:integer):regl;
begin
log:=ln(1.0*n)/1n(10.0)
end;
function TermNo(n,a：integer)：integer;
{ How many terms does one need to compute arctan 1/n so that the result
is correct to a decimal places? }
var m:integer;
x:real;
begin
m:=trunc((a+l)/log(n))+1；
while log(n)*(2*m+1)+log(2*m+1) > a+1-log(5) do m: =m-1;
TermNo:=m+1
end;
procedure invtan (n:integer； var atn: long);
{atn = arctan 1/n by MacLaurin series to m terms}
var m,nSq,i:integer;
pwr,trm:long;
begin
divl(one,n,pwr);
atn:=pwr;
m:=TermNo(n,4*QdPlaces);
if n>108 then for i:=1 to m do begin
divl(pwr,n,trm);
divl(trm,n,pwr);
divl(pwr, 2*i+1,trm);
if odd(i) then sub (atn, trm, atn) else add(atn ,trm, atn)
end else begin
nSq:=n*n;
for i:=1 to m do begin
divl(pwr,nSq,trm);
pwr:=trm;
divl(pwr,2*i+l,trm);
if odd(i) then sub (atn, trm atn) else add(atn, trm, atn)
end
end
end;
begin
one[0]:=1;
for i:=1 to QdPlaces do one[i]:=0;
TextBackground(blue);
TextColor(LightMagenta);
clrscr;
WriteWhere(tx);
invtan(10,result);
writeln (tx,'arctan 1/10 =')；outp(result);
mul(result,32,pi);
invtan(239,result);
writeln(tx,'arctan 1/239 =');outp(result);
mul(result,4,summand);
sub(pi,summand,pi);
invtan(515,result);
writeln(tx,'arctan 1/515 =');outp(result);
mul(result,16,summand);
sub(pi,summand,pi);
writeln(tx,'pi'='); outp(pi);
close(tx)
end・
THE OUTPUT
arctan 1/10 =
0.0396 6865 2491 1620 2737 8446 1198 7802 0590 2432 7832 2504 3146 4801
5508 7768 1002 7747 4475 5065 4420 6126 2443 4286 3715 7955 8386 4088 2739
8969 5679 2706 6563 1569 1279 0302 0720 8528 3902 8312 4322 3413 1452 0899
9120 3117 7645 1880 7542 7113 8281 1288 2359 3605 1731 1328 0153 5327 4276
6156 3882 5823 5491 4853 2458 7272 1476 0662 8140 5298 9463 7174 9699 7574
7101 4359 6227 3736 4575 2379 9778 0795 2648 9347 2492 9628 1854 9098 8357
0529 6856 4858 4916 2407 1550 0922 7014 3997 1745 5311 3784 8904 7538 9119
6646 5570 9523 9507 2692 0497 5355 0270 0220 9240 7194 9393 9323 5771 9184
0717 8695 9221 5984 8828 3363 3310 2491 1444 5127 5118 7074 0166 4303 8293
1083 0089 1901 7716 5321 2275 3913 1671 4214 0417 0455 9491 7296 5694 8052
2593 2178 2056 1144 7695 7767 9157 2075 2408 4060 9141 3778 0643 1803 1790
7809 3812 3279 1593 8266 4054 1140 3984 6766 8319 8180 3959 2521 6005 4106
5270 9824 3563 7314 3136 7097 6752 4114 4531 5933 0006 7415 6728 7565 8425
9304 6214 6327 3018 3990 0018 6663 9306 2097 5055 5237 6351 2328 3681 8403
8223 7057 7903 5281 4738 1287 5664 3645 8264 5611 9169 5171 1686 9735 3318
8378 0600 0792 6636 7362 3526 8794 4206 0269 9244 3332 0755 7332 0047 0228
2987 2557 4646 0128 5061 4702 6067 0701 4040 4795 2889
arctan 1/239 =
0.0041 8407 6002 0747 2386 4538 2149 5928 5452 7410 4806 5307 6319 5082
7019 6128 8718 1778 3414 2289 3273 7826 0581 3622 9094 5497 5450 6664 4486
3756 0524 5839 4789 3118 6505 8922 1288 3309 2800 8462 7196 2330 7733 7594
7634 6033 1847 3414 5703 3198 6015 4548 1480 5992 4498 3021 1460 3912 5394
9527 6077 9688 1558 8812 7339 7953 3465 1804 5742 5481 3586 7464 4751 9791
0232 8309 7700 2064 6528 2763 4653 2969 1048 1838 6543 5607 8919 5914 5123
2220 9446 3686 2766 1552 0831 6796 4264 6574 6551 1032 5103 4352 6282 4451
2693 5567 0499 6844 4524 7904 3317 7283 9307 0863 1401 9386 9519 5037 0586
4107 7085 5855 4045 2235 5388 1423 7677 0836 5156 9182 5270 2002 2930 8954
4950 0435 8544 0934 4964 4014 2418 7249 5092 2838 6239 5455 3335 6511 7197
3747 0202 3494 7597 7909 7469 5011 1888 5476 6739 7957 3153 7093 0327 8211
3089 8425 8308 3677 1909 1008 3909 8516 5510 4192 2416 7809 2053 2686 4916
2667 4027 1684 4424 4773 1579 6452 0275 4957 4158 8258 2909 4058 5090 3820
7331 7590 8431 9977 8432 7604 2858 6383 7347 4688 2493 4807 6992 8733 1405
2238 9423 0869 0507 3039 5440 6912 5801 0248 6061 6072 3699 5039 10S4 7488
1561 7625 8645 5122 5901 5542 9569 7941 6157 1069 4111 3503 6698 7619 4315
5931 8681 7207 4076 9131 0279 7191 2912 2049 3392 1211
arctan 1/515 =
0.0019 4174 5132 4432 9638 6842 4745 6125 0887 0389 7154 1220 7771 9451
3328 6133 7902 1524 9168 5752 3108 5274 2083 1284 6990 4648 8654 7439 6343
7834 6793 7424 6332 0744 6983 3038 6489 5426 4535 6102 7569 9824 7259 1644
2234 4949 3435 1098 1828 1040 7281 5932 5750 6591 3933 8597 7713 9530 6037
0984 2451 1723 1271 6472 9539 2364 6867 9559 1450 0072 4197 0383 3673 8644
5315 0563 9360 9043 8571 1642 9962 7397 6779 8058 0066 2281 1769 6489 1104
6106 6344 2979 5981 6413 4976 5059 7870 7901 0508 8108 7337 3141 9528 6129
0622 6567 5259 7467 0545 0195 2811 5519 9689 6111 7551 9173 9419 3189 6435
5863 5063 4301 7322 0215 3197 5649 9517 1664 6424 7539 8542 8495 1033 4377
4382 4325 2278 8958 1311 5404 9062 9519 7303 1800 8881 7715 4083 0719 3871
392b 9305 4087 7443 5124 3913 2825 0071 2555 4075 0287 6932 4689 4323 4850
9789 9752 6022 2884 1962 0447 5612 6595 0408 0012 8929 3828 5782 6269 1177
6002 5601 8991 2387 5196 6114 3468 6091 3070 5907 0941 8166 4942 9875 9783
8654 36Z 2223 2741 3311 6145 0424 1729 3262 2257 2162 7484 7764 2522 0552
0572 9721 1019 8675 4335 1195 3485 6620 8425 7987 6065 7033 5825 5963 4333
4659 2C44 7839 6341 0277 1488 4498 9199 8601 3074 2019 0674 5991 9575 5266
5708 7837 4822 6324 4713 9786 0548 1177 3136 0216 2851
pi =
3.1415 9265 3589 7932 3846 2643 3832 7950 2884 1971 6939 9375 1058 2097
4944 5923 0781 6406 2862 0899 8628 0348 2534 2117 0679 8214 8086 5132 8230
6647 0938 4460 9550 5822 3172 5359 4081 2848 1117 4502 8410 2701 9385 2110
5559 6446 2294 8954 9303 8196 4428 8109 7566 5933 4461 2847 5648 2337 8678
3165 2712 0190 9145 6485 6692 3460 3486 1045 4326 6482 1339 3607 2602 4914
1273 7245 8700 6606 3155 8817 4881 5209 2096 2829 2540 9171 5364 3678 9259
0360 0113 3053 0548 8204 6652 1384 1469 5194 1511 6094 3305 7270 3657 5959
1953 0921 8611 7381 9326 1179 3105 1185 4807 4462 3799 6274 9567 3518 8575
2724 8912 2793 8183 0119 4912 9833 6733 6244 0656 6430 8602 1394 9463 9522
4737 1907 0217 9860 9437 0277 0539 2171 7629 3176 7523 8467 4818 4676 6940
5132 0005 6812 7145 2635 6082 7785 7713 4275 7789 6091 7363 7178 7214 6844
0901 2249 5343 0146 5495 8537 1050 7922 7968 9258 9235 4201 9956 1121 2902
1960 8640 3441 8159 8136 2977 4771 3099 6051 8707 2113 4999 9998 3729 7804
9951 0597 3173 2816 0963 1859 5024 4594 5534 6908 3026 4252 2308 2533 4468
5035 2619 3118 8171 0100 0313 7838 7528 8658 7533 2083 8142 0617 1776 6914
7303 5982 5349 0428 7554 6873 1159 5628 6388 2353 7875 9375 1957 7818 5778
0532 1712 2680 6613 0019 2787 6611 1959 0921 6420 1989
---本文作者現任教於國立中央大學數學系---