多倍長Ｃ言語インタープリタ

ＴＣ

Ｔｉｎｙ Ｃ

福田隆 (日本大学生産工学部)

nekoma@math.cit.nihon-u.ac.jp

第 3 回 JANT シンポジウム２００１年１月２７日

速度機能プログラミング
Mathematica 遅い高機能容易ではない

KASH 遅い高機能 FORTRAN or BASIC 好

きに向いている

PARI  速い高機能難しい
UBASIC 速い低水準 BASIC 好きに向いている
TC  速い低水準 C 好きに向いている

TC の特徴

• C で書かれており、 UNIX マシン上で動作

• 多倍長演算部はアセンブラ (Alpha, Sparc, Pentium)

• C 準拠なのでプログラミングが容易

• 3000 桁程度までの高精度数値計算を得意とする

• int, real, rational, complex が変数の主な型

どうすれば入手できるか

• TC はフリーソフトです。バイナリパッケージはありませんがコンパイルは容易です。

• 東京都立大学の ftp サーバー

ftp://tnt.math.metro-u.ac.jp/pub/math-packs/tc/

に置かせて頂いています。

• ねこまのホームページ

http://www.math.cit.nihon-u.ac.jp/~nekoma/

にも置いてあります。

何故コンパイラでなくインタープリタなのか

• コンパイラを作るのは大変

• インタープリタの方が TC 上でのプログラム作成とバグとりが容易

• 電卓代わりに使える

• 計算の精度が上がればインタープリトにかかる時間は多倍長計算に吸収され無視できる。

Ｃと似ている点

• 大域変数と局所変数を区別する • 配列、ポインタが使える

• 関数の再帰呼出しができる • 配列の動的確保ができる

• if, while, for, goto が使える

Ｃと異なる点

• 構造体がない

• 関数へのポインタが使えない • char がない

• do while がない

• switch 文がない

ＴＣの演算子

演 算 子結合規則

() [] {}  左から右
! ++ -- * & (type)  右から左

^右から左
+ - (符号)  右から左

* / // %  左から右

+ - (加減)  左から右

<< >>  左から右
< <= > >=  左から右
== !=  左から右

& |  左から右

&& ||  左から右

?:  左から右

= ^= *= /= //= \%= += -= <<= >>= &= |= 右から左

,  左から右

コンパイラとの比較

a(n) = ∑_{k=1}^{n} k

C コンパイラは 64 bit long で計算

n  TC C コンパイラ TC/C
10^{8} 126 秒  0.2 秒  630

b(n) = ∏_{k=1}^{n} 2 sin(\sqrt{2}k)

1000 桁で計算

n  TC C コンパイラ TC/C
10^{4} 28.7 秒  28.6 秒  1.00

tc-3.10, pari-2.1, kash-2.1 on Alpha 21264 667MHz (ミリ秒)

 桁数 exp(\sqrt{2}) sinh(\sqrt{2}) sin(\sqrt{2}) log(\sqrt{2}) tan^{−1}(\sqrt{2})
200 0.11  0.13  0.16  0.20  0.37
TC  500 0.47  0.54  0.70  1.3  1.9
1000  1.9  2.0  2.7  5.4  7.0
2000  8.5  8.8  11.7  15  35
200 0.42  0.42  0.37  0.83  1.5
PARI 500  1.6  1.7  1.6  3.0  5.7
1000  6.1  6.2  6.2  10  20
2000  26  26  28  40  77
200  3.2  10  9.2
KASH 500  19  63  48
1000  92  330  200
2000  460  1700  920

tc-3.10, pari-2.1, kash-2.1 on Pentium II 333MHz (ミリ秒)

 桁数 exp(\sqrt{2}) sinh(\sqrt{2}) sin(\sqrt{2}) log(\sqrt{2}) tan^{−1}(\sqrt{2})
200  0.7  0.8  1.0  1.2  2.2
TC  500  3.2  3.5  4.7  5.6  8.8
1000  12  13  17  24  36
2000  42  44  60  110  147
200 0.85  0.88  0.85  1.5  2.9
PARI 500  4.8  4.9  5.4  7.6  16
1000  23  24  26  33  73
2000  120  120  140  160  360
200  12  34  22
KASH 500  75  230  120
1000  350  1200  530
2000 1700  6300  2500

\underline{自然対数の底 e の近似計算}

f(n) = ∑_{k=0}^{n} \frac{1}{k!}

TC
f(n)
int n;
{    real x,y;int k;
     for(y=x=k=1;k<=n;k++) y+=x/=k;
     return(y);}
PARI
f(n)=x=1.0;y=1.0;for(k=1,n,x=x/k;y=y+x);y

KASH
F:=function(n)
local x,y,k;
x:=1.0;y:=1.0;
for k in [1..n] do x:=x/k;y:=y+x; od;
return(y);
end;

1000 桁の精度での計算時間

Alpha 21264 667MHz

tc-3.10 pari-2.1 kash-2.1 KASH/TC
f(1000) 0.01 秒 0.03 秒 41 秒  4100
f(10000) 0.34 秒 0.35 秒 388 秒  1141

\underline{Theta 関数の近似値の計算}

Θ(a, b, c, d; n) = ∑_{x^{2}+y^{2}≤n} exp(2πiA(a, b, c, d, x, y))
A(a, b, c, d, x, y) = \frac{1}{10q^{2}} ((a + qx)^{2}(2 + ζ − ζ^{3} − 2ζ^{4})

+2(2 − ζ + ζ^{2} − 2ζ^{3})(a + qx)(b + qy) +(ζ + 2ζ^{2} − 2ζ^{3} − ζ^{4})(b + qy)^{2})

+\frac{1}{q^{2}} (c(a + qx) + d(b + qy))

但し q = 18, ζ = exp(2πi/5) で x, y は整数を動く。

Theta(n)         /* TC */
int n;
{    int x,y,m;complex z;
     z=0;m=isqrt(n);
     for(x=-m;x<=m;x++) for(y=-m;y<=m;y++)
        if(x*x+y*y<=n) z+=exp(2*_pi*_i*A(1,2,3,4,x,y));
     return(z);}
A(a,b,c,d,x,y)
int a,b,c,d,x,y;
{    complex z;int q;
     q=18;z=exp(2*_pi*_i/5);
     return(((a+q*x)^2*(2+z-z^3-2*z^4)
             +2*(2-z+z^2-2*z^3)*(a+q*x)*(b+q*y)+
             (z+2*z^2-2*z^3-z^4)*(b+q*y)^2)/(10*q^2)+
             (c*(a+q*x)+d*(b+q*y))/(q^2));}

Theta:=function(n)     /* KASH */
local x,y,z,m,I;
I:=Comp(0,1);z:=0;m:=Floor(Sqrt(n));
for x in [-m..m] do for y in [-m..m] do
  if x*x+y*y <= n then z:=z+Exp(2*pi*I*A(1,2,3,4,x,y)); fi;
od; od;
return(z);
end;
A:=function(a,b,c,d,x,y)
local q,z,I;
I:=Comp(0,1);q:=18;z:=Exp(2*pi*I/5);
return(((a+q*x)^2*(2+z-z^3-2*z^4)
       +2*(2-z+z^2-2*z^3)*(a+q*x)*(b+q*y)
       +(z+2*z^2-2*z^3-z^4)*(b+q*y)^2)/(10*q^2)
       +(c*(a+q*x)+d*(b+q*y))/(q^2));
end;

PARI
Theta(n,x,y,m,z)=z=0;m=sqrtint(n);\
for(x=-m,m,for(y=-m,m,\
if(x*x+y*y<=n,z=z+exp(2*Pi*I*A(1,2,3,4,x,y)),)));z
A(a,b,c,d,x,y,z,q)=q=18;z=exp(2*Pi*I/5);
((a+q*x)^2*(2+z-z^3-2*z^4)+2*(2-z+z^2-2*z^3)*(a+q*x)*(b+q*y)\
+(z+2*z^2-2*z^3-z^4)*(b+q*y)^2)/(10*q^2)\
+(c*(a+q*x)+d*(b+q*y))/(q^2)

500 桁での計算時間 Alpha 21264 667MHz

tc-3.10 pari-2.1 kash-2.1 KASH/TC
Theta(100) 0.81 秒 1.9 秒 296 秒  365
Theta(1000) 8.1 秒 19 秒 2959 秒  365
sin(\sqrt{2})  0.70 ms 1.6 ms 63 ms  90