mpmath 速度比べ

import mpmath
import numpy as np
import time

sample=500
iter=500

x = np.array([1.0/3.0]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + 1.0/3.0
    else:
        x = x - 1.0/3.0
goal=time.time()
print("---- numpy array -----")
print("%.20f" % init)
print("%.20f" % x[0])
print("time:", goal-start)

mpmath.mp.dps = 18
x = np.array([mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + mpmath.mpf("1.0")/mpmath.mpf("3.0")
    else:
        x = x - mpmath.mpf("1.0")/mpmath.mpf("3.0")
goal=time.time()
print("---- mpmath with numpy array (dps=%d) -----" % mpmath.mp.dps)
print("%.20s" % init)
print("%.20s" % x[0])
print("time:", goal-start)

mpmath.mp.dps = 100
x = np.array([mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + mpmath.mpf("1.0")/mpmath.mpf("3.0")
    else:
        x = x - mpmath.mpf("1.0")/mpmath.mpf("3.0")
goal=time.time()
print("---- mpmath with numpy array (dps=%d) -----" % mpmath.mp.dps)
print("%.20s" % init)
print("%.20s" % x[0])
print("time:", goal-start)

mpmath.mp.dps = 18
x = mpmath.matrix([mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + mpmath.mpf("1.0")/mpmath.mpf("3.0")
    else:
        x = x - mpmath.mpf("1.0")/mpmath.mpf("3.0")
goal=time.time()
print("---- mpmath matrix (dps=%d)-----" % mpmath.mp.dps)
print("%.20s" % init)
print("%.20s" % x[0])
print("time:", goal-start)

mpmath.mp.dps = 100
x = mpmath.matrix([mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample)
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = x + mpmath.mpf("1.0")/mpmath.mpf("3.0")
    else:
        x = x - mpmath.mpf("1.0")/mpmath.mpf("3.0")
goal=time.time()
print("---- mpmath matrix (dps=%d)-----" % mpmath.mp.dps)
print("%.20s" % init)
print("%.20s" % x[0])
print("time:", goal-start)


x = [mpmath.mpf("1.0")/mpmath.mpf("3.0")]*sample
init=x[0]
start = time.time()
for i in range(iter):
    if i<iter/2:
        x = [xx + mpmath.mpf("1.0")/mpmath.mpf("3.0") for xx in x]
    else:
        x = [xx - mpmath.mpf("1.0")/mpmath.mpf("3.0") for xx in x]
goal=time.time()
print("---- list -----")
print("%.20f" % init)
print("%.20f" % x[0])
print("time:", goal-start)

実行結果

$ python speed_test_array.py 
---- numpy array -----
0.33333333333333331483
0.33333333333332632042
('time:', 0.0022330284118652344)
---- mpmath with numpy array (dps=18) -----
0.333333333333333333
0.333333333333333326
('time:', 1.7795300483703613)
---- mpmath with numpy array (dps=100) -----
0.333333333333333333
0.333333333333333333
('time:', 2.18923282623291)
---- mpmath matrix (dps=18)-----
0.333333333333333333
0.333333333333333326
('time:', 5.815176010131836)
---- mpmath matrix (dps=100)-----
0.333333333333333333
0.333333333333333333
('time:', 6.364114046096802)
---- list -----
0.33333333333333331483
0.33333333333333331483
('time:', 12.696530103683472)

gccの4倍精度ライブラリー

gcc 4.6 から4倍精度計算ができる様になりました．昔は4倍精度で計算できるマシンがあったみたいだけどね．

今回は、c, (c++はやってない), fortranの言語で実装するとしてどの言語が最も実装しやすいか比較する．
gcc のversionは4.8 を用いる．結論としてはfortranが一番楽だと思う．
cおよびc++は
http://homepage3.nifty.com/maikaze/float128.html
で解説されているようにc(, c++)では幾分実装が不自然な所もあるので，一応比較してみる．
計算する内容は1/3, 4*arctan(1)=pi,sqrt(2)である．

まずfortranから

program fort_quadmath
 implicit none
 real (kind(0D0)) dx, dy, dpi,dsqrt2
 real (kind(0Q0)) qx, qy, qpi,qsqrt2
 dx = 1.0D0;
 dy = 3.0D0;
 dsqrt2 = sqrt(2.0D0)
 dpi = 4.0D0*atan(1.0D0)
 qx = 1.0Q0;
 qy = 3.0Q0;
 qsqrt2 = sqrt(2.0Q0)
 qpi = 4.0Q0*atan(1.0Q0)
 print *, "---- double precision ----"
 print *, "1/3     =", dx/dy
 print *, "pi      =", dpi
 print *, "sqrt(2) =", dsqrt2
 print *, "----  quad  precision ----" 
 print *, "1/3     =", qx/qy
 print *, "pi      =", qpi
 print *, "sqrt(2) =", qsqrt2
end program fort_quadmath 

コンパイルは

 gfortran -o fort_quadmath fort_quadmath.f03 

であり，実行結果は

$ ./fort_quadmath 
 ---- double precision ----
 1/3     =  0.33333333333333331     
 pi      =   3.1415926535897931     
 sqrt(2) =   1.4142135623730951     
 ----  quad  precision ----
 1/3     =  0.333333333333333333333333333333333317      
 pi      =   3.14159265358979323846264338327950280      
 sqrt(2) =   1.41421356237309504880168872420969798   

 pi           3.1415926535897932384626433832795028841971693993751
 sqrt(2)      1.4142135623730950488016887242096980785696718753769

厳密解と比較して32桁くらいまでは遜色無く一致している．
実装に関して変数の型の宣言と、リテラル(0Q0)さえ気をつけていれば，関数の型を意識する必要がない(というか関数が代表名で与えられているので，変数の型によって関数が勝手に陰的に変わるんだと理解している．)

一方c言語に関して

#include <stdio.h>
#include <quadmath.h>
#include <math.h>

int main(void)
{
 char c1[80],c2[80],c3[80],c4[80],c5[80];
 double dx,dy,dpi;
 long double lx,ly,lpi;
 __float128 x,y, pi1,pi2,sqrt1,sqrt2;
 dx = 1.0;
 dy = 3.0;
 dpi = 4.0*atan(1.0);
 lx = 1.0L;
 ly = 3.0L;
 lpi = 4.0L*atan(1.0L);
 x = 1.0Q;
 y = 3.0Q;
 pi1 = 4.0Q*atan(1.0Q);
 pi2 = 4.0Q*atanq(1.0Q); 
 sqrt1 = sqrt(2.0Q);
 sqrt2 = sqrtq(2.0Q);  
 
 printf("---- double precision ----\n");
 printf("1/3 = %.30lf\n", dx/dy);
 printf("pi  = %lf\n", dpi); 
 printf("---- long double prescision ----\n");
 printf("1/3 = %.30Lf\n", lx/ly);
 printf("pi  = %.30Lf\n", lpi);
 printf("---- quad precision ----\n");
 quadmath_snprintf(c1,80,"%.40Qf", x/y);
 quadmath_snprintf(c2,80,"%.40Qf", pi1); 
 quadmath_snprintf(c3,80,"%.40Qf", pi2);  
 quadmath_snprintf(c4,80,"%.40Qf", sqrt1);  
 quadmath_snprintf(c5,80,"%.40Qf", sqrt2);    
 printf("1/3 = %s\n", c1);
 printf("pi1  = %s\n", c2);
 printf("pi2  = %s\n",c3);
 printf("sqrt1  = %s\n",c4);
 printf("sqrt2  = %s\n",c5);  
 return 0;
}

のコンパイルは

gcc -o c_quadmath c_quadmath.c -lquadmath

で行う．実行結果は

$ ./c_quadmath 
---- double precision ----
1/3 = 0.333333333333333314829616256247
pi  = 3.141593
---- long double prescision ----
1/3 = 0.333333333333333333342368351437
pi  = 3.141592653589793115997963468544
---- quad precision ----
1/3 = 0.3333333333333333333333333333333333172839
pi1  = 3.1415926535897931159979634685441851615906
pi2  = 3.1415926535897932384626433832795027974791
sqrt1  = 1.4142135623730951454746218587388284504414
sqrt2  = 1.4142135623730950488016887242096981769402

 
pi     3.1415926535897932384626433832795028841971693993751
sqrt(2)  1.4142135623730950488016887242096980785696718753769

の厳密解と比較して pi1とsqrt1がdouble程度の精度しかない． これは代入する変数の型と関数の返り値の型が一致していないため，関数の方の型で丸められてしまったせいだと考えられる．(c++でも？)
既存のプログラムの精度を倍精度から4倍精度に書き直すにはc言語では関数も書き換えなければならない． (gccがversion upしたらなんとかなるんか？) 精度の仕様変更に関しては圧倒的にfortranの方がいじりやすいと個人的に思う．
 (fortran はあんまり使ってないけどね・・・)

gnuplot animation in C

　とある需要で，c言語でreal time animationを作る事になった．
gnuplotで描画させるのが一番簡単だろうなーと思ったので，基本的な操作方法をまとめる．
gnuplotを立ち上げて

gnuplot> plot "-"

と打つとデータの入力待ちになるので，

gnuplot> plot[-1:1][-1:1] "-" w p ps 7 pt 7
input data ('e' ends) > 0.0 0.0 
input data ('e' ends) > e

とかしてやると，データを手動で入力することができる．gnuplotを日常的に使っている人ならば直感的に，

gnuplot> plot[-1:1][-1:1] "-" w p ps 7 pt 7,"-" w p ps 7 pt 7
input data ('e' ends) > 0.5 0.0
input data ('e' ends) > e
input data ('e' ends) > -0.5 0.0
input data ('e' ends) > e

と二つ以上のデータ点を入力することができることに気がつくだろう．

結果は上図の様になる．つまり，cからgnuplotに上記のような命令でデータを入力してやれば，real time にデータを更新しながら描画することが可能になる．
(note:multiplotを利用するとデータ更新毎に再描画してしまうので，ちらつきの原因となる)
さて，ここからが本題．C言語からgnuplotに命令するため，ここではpopen関数をつかう．(多分やり方は色々有るはず)

FILE *gp;
gp = popen("gnuplot", "w");

でgnuplotにパイプを通すことができる．もちろんgnuplotのPATHが通ってなければならない．絶対PATHでも構わない． もし描画終了後にgnplotのウィンドウ画面が閉じるのが嫌であればgnuplotの起動オプション

gp = popen("gnuplot -persist", "w"); 

を付けてやれば良い．

以下簡単なサンプル

#include <stdio.h>
#include <stdlib.h>
#include <tmath.h>
#include <unistd.h>

void gnuplot_config(FILE *gp);
void plot2d(FILE *gp,const int num,const double *x,const double *y);
void move(const int num, double const t, double *x, double *y);
int main(void){
    FILE *gp;

    //gp = popen("gnuplot -persist", "w");    /*プログラムが終了後GnuplotのWindowを閉じない*/
    gp = popen("gnuplot", "w");               /*プログラムが終了後GnuplotのWindowを閉じる*/
    gnuplot_config(gp);  /* gnuplot configuration */


    int i;          // loop variable
    int step=100;
    int sample_num=10;
    double t;
    double dt=0.1;
    double x[sample_num], y[sample_num];
    
    printf("途中で終了する場合は，Ctrl-C をターミナルに入力\n");
    
    for (i=0;i<step;i++){
        fprintf(gp,"set title 'animated-%d/%d' font 'Times, 20'\n",i,step);    
        t = i*dt;
        move(sample_num, t, x, y);
        plot2d(gp,sample_num, x, y); 
        usleep(30000); //指定された(マイクロ?)秒プログラムを停止させます (Unix system only)
//        while(getchar() != '\n');       /*Enter keyが入力されるまで，計算停止*/        
    }
    printf("計算終了\n");


    return 0;
}

void gnuplot_config(FILE *gp)
{
/*
    fprintf(gp,"set terminal gif animate optimize size \n");
    fprintf(gp,"set output 'anime.gif'\n");
*/
    /*アニメーションをgif動画として保存する場合は，上のコメントを解除*/    
    fprintf(gp,"set xrange [-1.5:1.5]\n");
    fprintf(gp,"set yrange [-1.5:1.5]\n");
    fprintf(gp,"unset tics\n");
//    fprintf(gp,"set xtics 0.5\n");
//    fprintf(gp,"set ytics 0.5\n");
//    fprintf(gp,"set xlabel 'x' font 'Times, 20'\n");
//    fprintf(gp,"set ylabel 'y' font 'Times, 20'\n");
//    fprintf(gp,"set grid\n");
    fprintf(gp,"set border 0\n");
    fprintf(gp,"set size square\n");
}

void plot2d(FILE *gp, const int num, const double *x, const double *y)
{
    int i;
    fprintf(gp,"plot");
    for(i=0;i<num;i++){
        if (i!=num-1)
            fprintf(gp,"'-' with p pt 7 ps 10 lc %d tit '',", i);//plotデータの読み込み開始
        else
            fprintf(gp,"'-' with p pt 7 ps 10 lc %d tit ''\n", i);//plotデータの読み込み開始
    }
    for(i=0;i<num;i++){
        fprintf(gp,"%f %f\n",x[i],y[i]);
        fprintf(gp,"e\n"); //データ書き込み終了を知らせる．
    }
    fflush(gp); //バッファーに溜まったデーターを吐き出す．これをしないと，real time animation にならない．
}


void move(const int num, const double t, double *x, double *y)
{
    int i;
    double init_phase[num];
    for (i =0;i<num;i++){
        init_phase[i] =  i*2.0*M_PI/num;
    }

    for (i=0;i<num;i++){
        x[i] = cos(t+init_phase[i]);
        y[i] = sin(t+init_phase[i]);
    }
}

計算結果はこんな感じに

HTMLって鍵括弧(<>)を＆lt ＆gtって入力しないと行けないのか・・・プログラムコピペする時めんどくさいな・・・

今回はspleep関数使ってるけど，FPSとかってどうやって制御すれば良いんだ？

Hello world

男は黙ってHello world.

print("Hello world")