Hatena::ブログ(Diary)

arupaka-_-arupakaの日記

2013-11-19

cygwinのlapackで連立方程式をとく。

cygwinlapack連立方程式をとく。

まず、準備として、setup.exeからlapackっぽいのをいれておく。

とりあえずのプログラム

http://auewe.hatenablog.com/entry/2013/05/04/075717を参考にさせていただきました。

http://www.jsces.org/Issue/Journal/pdf/nakata-0411.pdf


ポイントは行列は横にみる。 N×M行列のAijは

a[(i-1)+M*(j-1)]

ただし,ladのときは違う。


簡単な例:

a(1,1)=1

a(1,2)=4

a(2,1)=5

a(2,2)=6

b(1)=3

b(2)=4

#include <stdio.h>
#include<stdlib.h>
#define N 2
int main(void)
{
  long i,j,n=N,nrhs=1,ipiv[N],lda=N,ldb=N,info;
   double* A;
   double* b;
   A=(double*)malloc(sizeof(double)*N*N);
   b=(double*)malloc(sizeof(double)*N);

   i=1;j=1;
   A[(i-1)+N*(j-1)]=1;
   i=1;j=2;
   A[(i-1)+N*(j-1)]=4;
   i=2;j=1;
   A[(i-1)+N*(j-1)]=5;
   i=2;j=2;
   A[(i-1)+N*(j-1)]=6;

    i=1;
    b[i-1]=3;
    i=2;
    b[i-1]=4;



  dgesv_(&n,&nrhs,A,&lda,ipiv,b,&ldb,&info);
  for(i=0;i<N;i++)
    printf("%lf\n",b[i]);
  return 0;
}

コンパイル

gcc -O3 solve.c -llapack

ランダムな例:

サンプルプログラム

$ cat lapack_test.c

#include <stdio.h>
#include<stdlib.h>
#define N 2000   // 2000x2000行列
int main(void)
{
  long i,j,n=N,nrhs=1,ipiv[N],lda=N,ldb=N,info;
   double* A;
   double* b;
   A=(double*)malloc(sizeof(double)*N*N);

   b=(double*)malloc(sizeof(double)*N);

 
  for(i=1;i<N*N;i++){
        A[i]=rand();
  }

  for(i=1;i<N;i++){

        b[i]=rand();
  }


  dgesv_(&n,&nrhs,A,&lda,ipiv,b,&ldb,&info);
  for(i=0;i<N;i++)
    printf("%lf\n",b[i]);
  return 0;
}

コンパイル

gcc -O3 lapack_test.c -llapack

実行

./a.exe

結果

0.346819

  • 0.144207

0.165927

  • 1.834181
  • 0.967015
  • 0.040578

1.918667

  • 0.092941
  • 2.201499

...

2013-10-22

windows上でのRとCの連携2。 RからC言語のgslのライブラリを関数(プログラム)を使用する。 cygwin環境編。

  • RとCの連携の基本は、

Windows上でのRとCの連携

http://d.hatena.ne.jp/arupaka-_-arupaka/20131020/1382252302 を参考に。

  • 具体的な手順

(1)まず、準備。mingwを利用して、cygwin上でcygwinに依存しないgslを作る。

(Rはcygwinライブラリを参照できないので、windowsdllを利用する必要がある)

http://d.hatena.ne.jp/arupaka-_-arupaka/20131022/1382412374 

(2)C言語のソースを用意.

http://d.hatena.ne.jp/arupaka-_-arupaka/20131020/1382252302 を参考に。

に、void型、値の受け渡しは、すべてポイント型。


例)平均mu 標準偏差sigma の正規分布乱数をnumber個発生する。

gsl_test2.c

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

void c_rnorm(int* number,double* mean,double* sigma,double* out){
        const gsl_rng_type *T;
        int i;
        gsl_rng *r;
        double v;
        T = gsl_rng_default;
        r = gsl_rng_alloc(T);
        gsl_rng_set(r,2);
        for(i=0;i<*number;i++){

                out[i]=gsl_ran_gaussian(r,*sigma)+*mean;
      
        }

}

(2)コンパイル

>|sh| 

x86_64-w64-mingw32-gcc gsl_test2.c -I/cygdrive/c/gsl/include -L/cygdrive/c/gsl/lib -lgsl -o gsl_test2.o -c

||<

x86_64-w64-mingw32-gcc -shared gsl_test2.o -o gsl_test2.dll -I/cygdrive/c/gsl/include -L/cygdrive/c/gsl/lib -lgsl

(3)Rのラッパー関数を作る。

c_norm<-function(n,mean1,sd1){

	.C("c_rnorm",arg1=as.integer(n),arg2=as.double(mean1),arg3=as.double(sd1),arg4=as.double(rep(0,n)))$arg4
}

(4)実行

dyn.load("gsl_test2.dll")
c_norm(1000,5,3)
||< 


(5)結果
>|c|
 1.09068286  6.92593977  3.04377685 -0.30245483  2.24621373
  6.95256996  2.27045627  4.57722179  7.34193082  8.47567727 ...

cygwinからWindows64bit上で動くgslライブラリを作る。 mingwとcygwinの連携

(1)まず、gslのソースをもってくる

http://ftp.jaist.ac.jp/pub/GNU/gsl/

(2)次にコンパイル

./configure --host=x86_64-w64-mingw32 --prefix=/cygdrive/c/gsl CC=x86_64-w64-mingw32-gcc.exe CXX=i686-pc-mingw32-g++.exe

make

make install

configure の意味は

configure --help

を参考に。

コンパイラCC=x86_64-w64-mingw32-gcc.exe をしておくこと(mingw64bitのコンパイラ)、

環境を x86_64-w64-mingw32 (mingw環境)

としておくことがポイント。

mingwがない場合はそれっぽいものをsetup.exeから落としてくることができる。

参考サイト

mingw な ruby2.0 を cygwin 上でビルドする

http://qiita.com/hidachinoiro/items/19d2083ecc1bb72e3ccd

http://www.bookshelf.jp/texi/autoconf-ja/autoconf-ja_13.html

configure host bild の違い

http://d.hatena.ne.jp/maminus/20100129/1264781242

(3)PATHの設定。Windowsの(システム)パスにc:\gsl\bin を加えておく。

   コントロールパネル->システム環境設定->環境変数

注意:WindowsにはユーザーPATHとシステムPATHがあるので注意が必要。

おそらくこれでOK.

(4)テスト

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

void main(){
 const gsl_rng_type *T;
 int i;
 gsl_rng *r;
 double sigma=1;
 double mean=0;
 int number=100;
 double v;
 //gsl_rng_setup();
 T = gsl_rng_default;
 r = gsl_rng_alloc(T);
 gsl_rng_set(r,2);
 for(i=0;i<number;i++){
 v=gsl_ran_gaussian(r,1.0);
 printf("%lf\n",v);
}

}

コンパイル (cygwin上)

x86_64-w64-mingw32-gcc gsl_test.c -I/cygdrive/c/gsl/include -L/cygdrive/c/gsl/lib -lgsl -o a2.out

実行 (cygwin上)

./a2.out

結果

  • 1.303106

0.641980

  • 0.652074
  • 1.767485
  • 0.917929

0.650857

  • 0.909848
  • 0.140926

0.780644

1.158559

  

   

2013-10-19

CをRから使う。

http://www.okada.jp.org/RWiki/?R%A4%AB%A4%E9%C2%BE%B8%C0%B8%EC%CD%F8%CD%D1を参考に

(1)Mingw を導入

(2) PATHを通す. システム→システムの詳細設定→環境変数

PATH ...

(3)Rtoolsを導入 http://essrc.hyogo-u.ac.jp/cran/ をダウンロード

(4)ダウンロードをしたのコピーしてそのフォルダにPATHを通す。

C:\Program Files\R\R-3.0.2\bin\x64とか(元から入っている場合も)。

gslでBスプラインでデータの平滑化

gslでBスプラインでデータの平滑化.

Bスプイラインは自然スプラインより計算が安定。

Bスプライン(拡張スプライン)同じところにノットを重ねることで

不連続を表現できる、3次だと、1回重ねると2回微分、2回重ねると1回微分、

3回重ねると関数そのものが不連続になる。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>

/* number of data points to fit */
#define N        1000

/* number of fit coefficients */
#define NCOEFFS 7

/* nbreak = ncoeffs + 2 - k = ncoeffs - 2 since k = 4 */
#define NBREAK   (NCOEFFS - 2)

void
spline(double* x_vector,double* y_vector, double* w_vector, int ndata,double* knots_vector,int nknots)
{
  const size_t n = (const size_t) ndata;
  const size_t ncoeffs = (const size_t) nknots+2;
  const size_t nbreak = (const size_t) nknots;
  size_t i, j;
  gsl_bspline_workspace *bw;
  gsl_vector *B;
  double dy;
  gsl_vector *c, *w;
  gsl_vector *x, *y;
  gsl_matrix *X, *cov;
  gsl_multifit_linear_workspace *mw;
  double chisq, Rsq, dof, tss;
 char* out_original_filename="out_original_data.dat";
 char* out_spline_filename="out_spline_data.dat";
  FILE* out_original;
  FILE* out_spline;

 // gsl_rng_env_setup();
//  r = gsl_rng_alloc(gsl_rng_default);

  /* allocate a cubic bspline workspace (k = 4) */
  bw = gsl_bspline_alloc(4, nbreak);
  B = gsl_vector_alloc(ncoeffs);

  x = gsl_vector_alloc(n);
  y = gsl_vector_alloc(n);
  X = gsl_matrix_alloc(n, ncoeffs);
  c = gsl_vector_alloc(ncoeffs);
  w = gsl_vector_alloc(n);
  cov = gsl_matrix_alloc(ncoeffs, ncoeffs);
  mw = gsl_multifit_linear_alloc(n, ncoeffs);

 out_original=fopen(out_original_filename,"w");
 out_spline=fopen(out_spline_filename,"w");

// double *knots_vector;
// knots_vector=(double*)malloc(sizeof(double)*ncoeffs+2);

// double *x_vector;
// x_vector=(double*)malloc(sizeof(double)*n);

// double *y_vector;
// y_vector=(double*)malloc(sizeof(double)*n);

 //double *w_vector;
 // w_vector=(double*)malloc(sizeof(double)*n);


  gsl_vector *knots_vector2;
//input the data
/*
 for(i=0;i<n;i++){


      double sigma;
      double xi = i;
      double yi = 0.01*(i+1);

      sigma = 0.1 * yi;
      dy = gsl_ran_gaussian(r, sigma);
      yi += dy;
      y_vector[i]=yi;
      x_vector[i]=xi;
      w_vector[i]=1.0/sigma*sigma;
 }
*/
//Set the knots
/*
 for(i=0;i<ncoeffs-2;i++){
     double dx=(x_vector[n-1]+0.0)/((ncoeffs-2)-1.0);
     knots_vector[i]=i*dx;
     printf("knots %lf %lf\n",x_vector[n-1],knots_vector[i]);
 }
*/


  printf("#m=0,S=0\n");
  /* this is the data to be fitted */
  for (i = 0; i < n; ++i)
    {
      double wi=y_vector[i];
      double xi = x_vector[i];
      double yi = y_vector[i];

     // sigma = 0.1 * yi;
      //dy = gsl_ran_gaussian(r, sigma);
      //yi += dy;

      gsl_vector_set(x, i, xi);
      gsl_vector_set(y, i, yi);
      gsl_vector_set(w, i, wi);

      fprintf(out_original,"%f %f\n", xi, yi);
    }
    knots_vector2=gsl_vector_alloc(ncoeffs-2);
    // gsl_vector_set(knots_vector2,0,min_x);
    for(i=0;i<ncoeffs-2;i++){
       gsl_vector_set(knots_vector2,i,knots_vector[i]);

     printf("%d\n",i);
     }
     // gsl_vector_set(knots_vector2,i,knots_vector[ncoeffs+1]);
  /* use uniform breakpoints on [0, 15] */
//  gsl_bspline_knots_uniform(0.0, gsl_vector_max(x), bw);
 printf("nbreak%d\n",bw->nbreak);
  gsl_bspline_knots(knots_vector2, bw);
  /* construct the fit matrix X */
  for (i = 0; i < n; ++i)
    {
      double xi = gsl_vector_get(x, i);

      /* compute B_j(xi) for all j */
      gsl_bspline_eval(xi, B, bw);

      /* fill in row i of X */
      for (j = 0; j < ncoeffs; ++j)
        {
          double Bj = gsl_vector_get(B, j);
          gsl_matrix_set(X, i, j, Bj);
        }
    }

  /* do the fit */
  gsl_multifit_wlinear(X, w, y, c, cov, &chisq, mw);

  dof = n - ncoeffs;
  tss = gsl_stats_wtss(w->data, 1, y->data, 1, y->size);
  Rsq = 1.0 - chisq / tss;

  fprintf(stderr, "chisq/dof = %e, Rsq = %f\n",
                   chisq / dof, Rsq);

  /* output the smoothed curve */
  {
    double xi, yi, yerr;

    printf("#m=1,S=0\n");
    double delta=1.0;
    for (xi = 0.0; xi < gsl_vector_max(x); xi += delta)
      {
        gsl_bspline_eval(xi, B, bw);
        gsl_multifit_linear_est(B, c, cov, &yi, &yerr);
        fprintf(out_spline,"%f %f\n", xi, yi);
      }
  }

//  gsl_rng_free(r);
  gsl_bspline_free(bw);
  gsl_vector_free(B);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_matrix_free(X);
  gsl_vector_free(c);
  gsl_vector_free(w);
  gsl_matrix_free(cov);
  gsl_multifit_linear_free(mw);
// free(knots_vector);
// free(x_vector);
// free(y_vector);
// free(w_vector);

 fclose(out_original);
 fclose(out_spline);

//  return 0;
}

int main(void){
        int i;
        int n=1000;
        int nknots=9;
        double dx,dy;

        //Initilaization of random variable
        gsl_rng *r;
        gsl_rng_env_setup();
        r = gsl_rng_alloc(gsl_rng_default);


        double *x_vector;
        x_vector=(double*)malloc(sizeof(double)*n);

        double *y_vector;
        y_vector=(double*)malloc(sizeof(double)*n);

        double *w_vector;
        w_vector=(double*)malloc(sizeof(double)*n);


        double *knots_vector;
        knots_vector=(double*)malloc(sizeof(double)*nknots);

        //Set the data
        for(i=0;i<n;i++){


             double sigma;
             double xi = i;
             double yi = 0.01*(i+1);

             sigma = 0.1 * yi;
             dy = gsl_ran_gaussian(r, sigma);
             yi += dy;
             y_vector[i]=yi;
             x_vector[i]=xi;
             w_vector[i]=1.0/sigma*sigma;
         }

        //Set the knots
         for(i=0;i<nknots;i++){
             double dx=(x_vector[n-1]+0.0)/((nknots)-1.0);
             knots_vector[i]=i*dx;
             printf("knots %lf %lf\n",x_vector[n-1],knots_vector[i]);
         }

        spline(x_vector,y_vector,w_vector,n,knots_vector,nknots);

         free(knots_vector);
         free(x_vector);
         free(y_vector);
         free(w_vector);


        return 0;


}

2013-10-14

gslで最適化

例:GSLで最適化

http://www.gnu.org/software/gsl/manual/html_node/Multimin-Examples.html#Multimin-Examples

http://www.gnu.org/software/gsl/manual/html_node/Providing-a-function-to-minimize.html#Providing-a-function-to-minimize

を参考に。。


微分がわかる場合。

#include<stdio.h>
#include<gsl/gsl_multimin.h>

double
my_f (const gsl_vector *v, void *params)
{
  double x, y;
  double *p = (double *)params;

  x = gsl_vector_get(v, 0);
  y = gsl_vector_get(v, 1);

  return p[2] * (x - p[0]) * (x - p[0]) +
           p[3] * (y - p[1]) * (y - p[1]) + p[4];
}

/* The gradient of f, df = (df/dx, df/dy). */
void
my_df (const gsl_vector *v, void *params,
       gsl_vector *df)
{
  double x, y;
  double *p = (double *)params;

  x = gsl_vector_get(v, 0);
  y = gsl_vector_get(v, 1);

  gsl_vector_set(df, 0, 2.0 * p[2] * (x - p[0]));
  gsl_vector_set(df, 1, 2.0 * p[3] * (y - p[1]));
}

/* Compute both f and df together. */
void
my_fdf (const gsl_vector *x, void *params,
        double *f, gsl_vector *df)
{
  *f = my_f(x, params);
  my_df(x, params, df);
}


int main (void)
{
  size_t iter = 0;
  int status;

  const gsl_multimin_fdfminimizer_type *T;
  gsl_multimin_fdfminimizer *s;

  /* Position of the minimum (1,2), scale factors
     10,20, height 30. */
  double par[5] = { 1.0, 2.0, 10.0, 20.0, 30.0 };

  gsl_vector *x;
  gsl_multimin_function_fdf my_func;

  my_func.n = 2;
  my_func.f = my_f;
  my_func.df = my_df;
  my_func.fdf = my_fdf;
  my_func.params = par;

  /* Starting point, x = (5,7) */
  x = gsl_vector_alloc (2);
  gsl_vector_set (x, 0, 5.0);
  gsl_vector_set (x, 1, 7.0);

  T = gsl_multimin_fdfminimizer_conjugate_fr;
  s = gsl_multimin_fdfminimizer_alloc (T, 2);

  gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-4);

  do
    {
      iter++;
      status = gsl_multimin_fdfminimizer_iterate (s);

      if (status)
        break;

      status = gsl_multimin_test_gradient (s->gradient, 1e-3);

      if (status == GSL_SUCCESS)
        printf ("Minimum found at:\n");

      printf ("%5d %.5f %.5f %10.5f\n", iter,
              gsl_vector_get (s->x, 0),
              gsl_vector_get (s->x, 1),
              s->f);

    }
  while (status == GSL_CONTINUE && iter < 100);

  gsl_multimin_fdfminimizer_free (s);
  gsl_vector_free (x);

  return 0;
}

微分を指定しない場合

#include<stdio.h>
#include<gsl/gsl_multimin.h>

double
my_f (const gsl_vector *v, void *params)
{
  double x, y;
  double *p = (double *)params;

  x = gsl_vector_get(v, 0);
  y = gsl_vector_get(v, 1);

  return p[2] * (x - p[0]) * (x - p[0]) +
           p[3] * (y - p[1]) * (y - p[1]) + p[4];
}


/* Compute both f and df together. */
void
my_fdf (const gsl_vector *x, void *params,
        double *f, gsl_vector *df)
{
  *f = my_f(x, params);
  my_df(x, params, df);
}

int
main(void)
{
  double par[5] = {1.0, 2.0, 10.0, 20.0, 30.0};

  const gsl_multimin_fminimizer_type *T =
    gsl_multimin_fminimizer_nmsimplex2;
  gsl_multimin_fminimizer *s = NULL;
  gsl_vector *ss, *x;
  gsl_multimin_function minex_func;

  size_t iter = 0;
  int status;
  double size;

  /* Starting point */
  x = gsl_vector_alloc (2);
  gsl_vector_set (x, 0, 5.0);
  gsl_vector_set (x, 1, 7.0);

  /* Set initial step sizes to 1 */
  ss = gsl_vector_alloc (2);
  gsl_vector_set_all (ss, 1.0);

  /* Initialize method and iterate */
  minex_func.n = 2;
  minex_func.f = my_f;
  minex_func.params = par;

  s = gsl_multimin_fminimizer_alloc (T, 2);
  gsl_multimin_fminimizer_set (s, &minex_func, x, ss);

  do
    {
      iter++;
      status = gsl_multimin_fminimizer_iterate(s);

      if (status)
        break;

      size = gsl_multimin_fminimizer_size (s);
      status = gsl_multimin_test_size (size, 1e-2);

      if (status == GSL_SUCCESS)
        {
          printf ("converged to minimum at\n");
        }

      printf ("%5d %10.3e %10.3e f() = %7.3f size = %.3f\n",
              iter,
              gsl_vector_get (s->x, 0),
              gsl_vector_get (s->x, 1),
              s->fval, size);
    }
  while (status == GSL_CONTINUE && iter < 100);

  gsl_vector_free(x);
  gsl_vector_free(ss);
  gsl_multimin_fminimizer_free (s);

  return status;
}

2013-10-13

gslで線形方程式、連立一次方程式をとく。

gslで線形方程式、連立一次方程式をとく。

#include<stdio.h>
#include<gsl/gsl_linalg.h>

void solve(double** m,double* v,int len){
        using namespace std;
        int s;
        double *mb;
        int k=0;
        mb=(double*)malloc(sizeof(double)*len*len);
        for(int i=0;i<len;i++){
                for(int j=0;j<len;j++){
                      
                        mb[k]=m[i][j];
                        k++;

                }

        }


        gsl_matrix_view m2=gsl_matrix_view_array(mb,len,len);
        gsl_vector_view v2=gsl_vector_view_array(v,len);
        gsl_vector *x=gsl_vector_alloc(len);
        gsl_permutation *p =gsl_permutation_alloc(len);
        gsl_linalg_LU_decomp(&m2.matrix,p,&s);
        gsl_linalg_LU_solve(&m2.matrix,p,&v2.vector,x);
        printf("x=\n");

        gsl_vector_fprintf(stdout,x,"%g");
        gsl_permutation_free(p);
        gsl_vector_free(x);

//      return *v;

}


int main(){

        double** m;
        double* v;
        double* x;
        int len=3;


        m=(double**) malloc(sizeof(double*)*(len));
        v=(double*) malloc(sizeof(double)*(len));
        x=(double*) malloc(sizeof(double)*(len));

        for(int i=0;i<len;i++){
                m[i]=(double*)malloc(sizeof(double)*(len));
        }

        m[0][0]=3;
        m[0][1]=2;
        m[0][2]=5;
        m[1][0]=7;
        m[1][1]=4;
        m[1][2]=1;
        m[2][0]=9;
        m[2][1]=2;
        m[2][2]=6;

        v[0]=2;
        v[1]=4;
        v[2]=7;

        solve(m,v,3);


}

x=

0.818182

マイナス0.454545

0.0909091

帰り値あり

#include<stdio.h>
#include<gsl/gsl_linalg.h>

void solve(double** m,double* v,int len,double* output){
        using namespace std;
        int s;
        double *mb;
        int k=0;
        mb=(double*)malloc(sizeof(double)*len*len);
        for(int i=0;i<len;i++){
                for(int j=0;j<len;j++){
//                      printf("%lf\n",m[i][j]);
                        mb[k]=m[i][j];
                        k++;

                }

        }


        gsl_matrix_view m2=gsl_matrix_view_array(mb,len,len);
        gsl_vector_view v2=gsl_vector_view_array(v,len);
        gsl_vector *x=gsl_vector_alloc(len);
        gsl_permutation *p =gsl_permutation_alloc(len);
        gsl_linalg_LU_decomp(&m2.matrix,p,&s);
        gsl_linalg_LU_solve(&m2.matrix,p,&v2.vector,x);
        printf("x=\n");

        gsl_vector_fprintf(stdout,x,"%g");
        gsl_permutation_free(p);
        //gsl_vector_free(x);
        for(int i=0;i<len;i++){
                output[i]=gsl_vector_get(x,i);

        };
//      return *v;

}


int main(){

        double** m;
        double* v;
        double* x;
        double* output;
        int time_length=3;


        m=(double**) malloc(sizeof(double*)*(time_length));
        v=(double*) malloc(sizeof(double)*(time_length));
        x=(double*) malloc(sizeof(double)*(time_length));

        output=(double*) malloc(sizeof(double)*(time_length));
        for(int i=0;i<time_length;i++){
                m[i]=(double*)malloc(sizeof(double)*(time_length));
        }

        m[0][0]=3;
        m[0][1]=2;
        m[0][2]=5;
        m[1][0]=7;
        m[1][1]=4;
        m[1][2]=1;
        m[2][0]=9;
        m[2][1]=2;
        m[2][2]=6;


/*
        m[0][0]=1;
        m[0][1]=0;
        m[0][2]=0;
        m[1][0]=0;
        m[1][1]=1;
        m[1][2]=0;
        m[2][0]=0;
        m[2][1]=0;
        m[2][2]=1;
*/
        v[0]=2;
        v[1]=4;
        v[2]=7;

        solve(m,v,3,output);
        printf("answer\n");
        for(int i;i<time_length;i++){
                printf("%lf\n",output[i]);

        }

}

出力

x=

0.818182

  • 0.454545

0.0909091

answer

0.818182

  • 0.454545

0.090909