//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
//  Program fracdeco.cc
//  Symplectic integrator のための fractal decomposition の
//  分解係数の定数配列を作成するプログラム
//  Copyright (C) by 渡辺尚貴 1998年 8月14日
//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=

/********* このプログラムが行うfractal分解の定義 *********

非可換演算子 Dk, Du の指数積の m次精度のfractal分解 S_{m}(Dk+Du) を
以下の再帰定義として計算します。

S_{m}(Dk+Du) =  S_{m-2}( s_{m-1}*(Dk+Du) )
              * S_{m-2}( s_{m-1}*(Dk+Du) )
              * S_{m-2}( (1-4s_{m-1})*(Dk+Du) )
              * S_{m-2}( s_{m-1}*(Dk+Du) )
              * S_{m-2}( s_{m-1}*(Dk+Du) )

ここで、  s_{m-1} = 1.0/(4.0-pow(4.0,1.0/(m-1))) です。

最後に  S_{2}(Dk+Du) = exp(0.5*Dk) * exp(Du) * exp(0.5*Dk) と定義します。

この再帰定義により、指数内の Dk, Duの係数が
Dk, Du, Dk, Dk, Du, Dk, Dk, Du, Dk, Dk, Du, Dk, ...
の順番で求まりますが、連続する Dk は一つの Dkに係数をまとめます。
また、この分解係数を用いるプログラムでの利便性のため、このプログラム
が出力する分解係数の先頭には係数 0 の Duの因子を置きます。

このプログラムは分解係数をC++の定数配列の定義の書式で出力します。
配列 symp_k[] には Dk の分解係数をその作用順序で出力します。
配列 symp_u[] には Du の分解係数をその作用順序で出力します。
整数 symp_r には配列の要素数を出力します。
整数 symp_m には分割の次数を出力します。

この係数を利用する際には Duから先に作用してください。
利用の際に便利なようにするため、このプログラムは symp_u[0] を 0 に
します。こうすると、Du, Dkを同じ回数だけ作用することになります。

分解係数などの出力は標準出力 stdout に対して行いますが
それ以外のメッセージは標準エラー出力 stderr に対して行いますので
この分解係数をファイルに redirect することもできます。

***************************************************************/


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

// 演算子の種類を表す列挙型の宣言
enum Operator{ Dk, Du };

// 分解係数を格納する配列のためのポインタの宣言と定義
static double *symp_k, *symp_u;

// その配列のための添字の宣言と定義
static int symp_rk, symp_ru;


//---- 次数の入力を受け付ける関数
int Input( void )
{
  int m;

  while(1){  // 正常な値が入力されるまで繰り返すループ
    rewind(stdin);
    fprintf( stderr, "Fractal分割の次数を入力して下さい (2,4,6,8) :");
    if( scanf("%d", &m ) != 1 ){
      puts("入力エラー");
      continue;
    }
    if( m%2 ){
      fprintf( stderr, "奇数次分割はもうひとつ高次の偶数次分割と同じです。\n");
      continue;
    }
    if( !( 2<=m && m<=8 ) ){
      fprintf( stderr, "分割数が異常です。\n");
      continue;
    }

    break;
  };

  return( m );
}


//---- 分解係数を格納するための配列を用意する関数
// 引数 m には分解次数を指定する
void Init( int m )
{
    // この fractal分解では Du, Dk の因子がこれだけある。
  int num = (int)pow(5, m/2-1) + 1;
  symp_k = new double [ num ]; // 配列のためのメモリ確保
  symp_u = new double [ num ];

  symp_rk = 0; // 添字を先頭に合わせる
  symp_ru = 0;
  symp_u[symp_ru++] = 0.0; // Du の最初の因子の係数は 0 にする。
}

//---- 計算された分解係数を格納する関数
// 引数 op には演算子の種類 (Du or Dk)を指定する
// 引数 s  にはその分解係数を指定する
void Add( Operator op, double s )
{
    // 前のAdd()関数で登録された係数の演算子の種類を記憶する変数
  static Operator prev_op = Du;  // 既に Du の係数 0 の因子がある。

  if( prev_op == op ){ // 前回の演算子と同じなら
    switch( op ){
      case Dk: symp_k[symp_rk-1] += s; break; // 今の添字の前の要素
      case Du: symp_u[symp_ru-1] += s; break; // に加算する
    }
  }else{               // 前回の演算子と異なるなら
    prev_op = op;
    switch( op ){
      case Dk: symp_k[symp_rk++] =  s; break; // 今の添字の要素に代入
      case Du: symp_u[symp_ru++] =  s; break; // して、添字を更新
    }
  }
}

//---- 分解を再帰的に行う関数
void FracDeco( int m, double s )
{
  if( m == 2 ){ // 2次まで分解したら係数を配列に登録する。
    Add( Dk, s*0.5 );
    Add( Du, s );
    Add( Dk, s*0.5 );
  }else{
      // 恒等分解のための係数
    const double s_m = 1.0/(4.0-pow(4.0,1.0/(m-1)));

      // fractal分解の実行
    FracDeco( m-2, s*s_m );
    FracDeco( m-2, s*s_m );
    FracDeco( m-2, s*(1.0-4.0*s_m) );
    FracDeco( m-2, s*s_m );
    FracDeco( m-2, s*s_m );
  }
}

//---- 分解係数をC++プログラム用に出力する関数
void Print( int m )
{
  int i;

  printf("\n\n");

  printf("// %d次 fractal分解の分解係数\n\n", m );

    // 分解の次数の出力
  printf("const double symp_m    = %d;\n", m );
    // 分解係数の要素数の出力
  printf("const double symp_r    = %d;\n\n", symp_rk );

    // Dkの分解係数の出力
  printf("const double symp_k[%d] = {\n", symp_rk );
  for( i=0 ; i<symp_rk ; i++ ){
    printf("\t%+15.15lf,\n", symp_k[i] );
  }
  printf("};\n\n");

    // Duの分解係数の出力
  printf("const double symp_u[%d] = {\n", symp_ru );
  for( i=0 ; i<symp_ru ; i++ ){
    printf("\t%+15.15lf,\n", symp_u[i] );
  }
  printf("};\n");

  fprintf( stderr, "分解係数が求まりました。\n");
}


int main( void )
{
  int m;

  m = Input();           // 次数のキー読み込み
  Init( m );             // 配列の初期化
  FracDeco( m, 1.0 );    // 分割の実行
  Print( m );            // 結果の出力
}