namespace BigNum
{
  /*
    CalcGCD : ユークリッドの互除法を使った最大公約数の計算

    x, y : 最大公約数を求める二数
  */
  template< class T >
  T CalcGCD( T x, T y )
  {
    if ( x < y ) std::swap( x, y );
    if ( y == 0 ) return( x );

    do {
      T r = x % y;
      x = y;
      y = r;
    } while ( y != 0 );

    return( x );
  }

  /*
    Reduce : 約分(最大公約数で二数を割る)

    x, y : 対象の二数へのポインタ
  */
  template< class T >
  void Reduce( T* x, T* y )
  {
    T gcd = CalcGCD( *x, *y );
    if ( gcd > 1 ) {
      *x /= gcd;
      *y /= gcd;
    }
  }

  /**
     BugNum::Rational : 多倍長有理数

     多倍長整数を分子・分母とする有理数の型
  **/
  class Rational
  {
    bool isPos_; // 符号は正か ?
    Unsigned a_; // 分子
    Unsigned b_; // 分母

  public:

    // コンストラクタ
    Rational( const Unsigned& a = 0, const Unsigned& b = 1, bool isPos = true );

    // 小数表記で文字列を返す
    std::string circulator() const;

    // 無効値判定
    bool isNaN() const { return( b_ == 0 ); }
  };

  /*
    Rational コンストラクタ(分子・分母指定)

    a : 分子
    b : 分母
    isPos : 符号は正か
  */
  Rational::Rational( const Unsigned& a, const Unsigned& b, bool isPos )
  {
    if ( a.isNaN() || b.isNaN() ) {
      b_ = 0;
    } else {
      a_ = a;
      b_ = b;
      isPos_ = isPos;
    }

    Reduce( &a_, &b_ );
  }

  /*
    Rational::circulator : 小数表記に変換した文字列を返す

    循環小数の場合、循環節を[]で示す
  */
  string Rational::circulator() const
  {
    if ( isNaN() ) return( "NaN" );

    Unsigned num( a_ );

    std::stringstream ss;
    ss << ( ( isPos_ ) ? "+" : "-" );

    // 商と余りを求める
    Unsigned r;
    num.div( b_, &r );
    // 整数部まで出力
    ss << num;
    if ( r != 0 ) ss << ".";

    // 小数部の計算
    std::map< Unsigned, std::vector< Unsigned >::size_type > flag; // 余りと位置を記録するmap
    std::vector< Unsigned > numVec; // 小数部の数値を記録するvector
    std::vector< Unsigned >::size_type i = 0;
    while ( r != 0 ) {
      if ( flag.find( r ) != flag.end() ) break;
      flag[r] = i++;
      num = r * 10;
      num.div( b_, &r );
      numVec.push_back( num );
    }

    // 小数部の出力
    for ( i = 0 ; i < numVec.size() ; ++i ) {
      if ( r != 0 && i == ( flag.find( r ) )->second ) ss << "[";
      ss << numVec[i];
    }
    if ( r != 0 ) ss << "]";

    return( ss.str() );
  }
}

namespace BigNum
{
  /*
    LTrim : 左側の空白をスキップした位置を返す

    str : 対象の文字列
    s : チェック開始位置

    戻り値 : 空白をスキップした最初の文字の位置(全て空白なら文字列長と等しい)
  */
  string::size_type LTrim( const string& str, string::size_type s )
  {
    while ( s < str.length() ) {
      if ( ! isspace( str[s] ) ) break;
      ++s;
    }

    return( s );
  }

  /*
    RTrim : 右側の空白をスキップした位置を返す

    str : 対象の文字列
    e : チェック開始位置(チェック対象の末尾の次の位置)

    戻り値 : 空白をスキップした末尾の次の文字の位置(全て空白ならゼロと等しい)
  */
  string::size_type RTrim( const string& str, string::size_type e )
  {
    if ( e > str.length() )
      e = str.length();

    while ( e > 0 ) {
      if ( ! isspace( str[e - 1] ) ) break;
      --e;
    }

    return( e );
  }

  /*
    Str2Num : 文字列→整数値変換

    str : 文字列
    e : 読み込み開始位置・終了位置 + 1
    num : 変換結果(絶対値)へのポインタ
    isPos : 変換結果(符号)へのポインタ

    (+/-)[0-9] 形式で、空白文字はチェックしない
    isPos が NULL ポインタの場合は 符号は取得しない

    戻り値 : 文字列が数値変換できなかった場合は false を返す
  */
  bool Str2Num( const string& str, string::size_type s, string::size_type e, Unsigned* num, bool* isPos )
  {
    if ( s >= e ) return( false );

    // 符号のチェック
    if ( isPos != 0 ) {
      if ( str[s] == '-' || str[s] == '+' ) {
	if ( str[s] == '-' ) *isPos = false;
	++s;
	if ( s >= e ) return( false );
      } else {
        *isPos = true;
      }
    }

    *num = Unsigned( str.substr( s, e - s ) );

    return( ! num->isNaN() );
  }

  /*
    Rational::toDecimal : 小数表記の文字列から有理数を構築

    const string& s : 文字列
  */
  void Rational::toDecimal( const string& s )
  {
    b_ = 0; // NaN にしておく

    string::size_type sp = LTrim( s, 0 );          // 開始位置
    string::size_type ep = RTrim( s, s.length() ); // 末尾の次の位置

    string::size_type dp = s.find( '.' ); // 小数点
    string::size_type rdp;                // 循環節の開始位置

    // 整数表記の場合
    if ( dp == string::npos || dp + 1 == ep ) {
      if ( dp == string::npos ) dp = ep;
      if ( ! Str2Num( s, sp, dp, &a_, &isPos_ ) ) return;
      b_ = 1;
      return;
    }

    // 有限小数を取得
    if ( ! Str2Num( s, sp, dp, &a_, &isPos_ ) ) return;
    for ( rdp = dp + 1 ; rdp < ep ; ++rdp ) {
      if ( s[rdp] == '[' ) break;
      if ( ! isdigit( s[rdp] ) ) return;
    }
    if ( rdp > dp + 1 ) {
      Unsigned fdA; // 有限小数部
      if ( ! Str2Num( s, dp + 1, rdp, &fdA, 0 ) ) return;
      b_ = 10;
      b_.pow( rdp - dp - 1 );
      a_ = a_ * b_ + fdA;
    }

    // 循環節を取得
    if ( rdp < ep ) {
      Unsigned rdA; // 循環小数部
      if ( s[ep - 1] != ']' ) return;
      if ( ! Str2Num( s, rdp + 1, ep - 1, &rdA, 0 ) ) return;
      if ( b_ == 0 ) b_ = 1;

      Unsigned rdB( 10 );
      rdB = ( rdB.pow( ep - rdp - 2 ) - 1 ) * b_;
      Rational rd( rdA, rdB, true );
      if ( isPos_ )
	operator+=( rd );
      else
	operator-=( rd );
    }

    Reduce( &a_, &b_ );
  }
}

namespace BigNum
{
  /*
    比較演算子の多重定義

    num : 比較対象の有理数

    大小比較は以下のように行う

    a/b < num.a/num.b → a*num.b < num.a*b
    a/b > num.a/num.b → a*num.b > num.a*b

    少なくとも一方が無効値なら常に false を返す
  */
  bool Rational::operator==( const Rational& num ) const
  {
    if ( isNaN() || num.isNaN() ) return( false );
    return( a_ == num.a_ && b_ == num.b_ && isPos_ == num.isPos_ );
  }
  bool Rational::operator!=( const Rational& num ) const
  {
    if ( isNaN() || num.isNaN() ) return( false );
    return( ! operator==( num ) );
  }
  bool Rational::operator<( const Rational& num ) const
  {
    if ( isNaN() || num.isNaN() ) return( false );

    if ( isPos_ && ( ! num.isPos_ ) ) return( false );
    if ( ( ! isPos_ ) && num.isPos_ ) return( true );

    return( ( isPos_ ) ? a_ * num.b_ < num.a_ * b_ : a_ * num.b_ > num.a_ * b_ );
  }
  bool Rational::operator>( const Rational& num ) const
  {
    if ( isNaN() || num.isNaN() ) return( false );

    if ( isPos_ && ( ! num.isPos_ ) ) return( true );
    if ( ( ! isPos_ ) && num.isPos_ ) return( false );

    return( ( isPos_ ) ? a_ * num.b_ > num.a_ * b_ : a_ * num.b_ < num.a_ * b_ );
  }
  bool Rational::operator<=( const Rational& num ) const
  {
    if ( isNaN() || num.isNaN() ) return( false );
    return( ! operator>( num ) );
  }
  bool Rational::operator>=( const Rational& num ) const
  {
    if ( isNaN() || num.isNaN() ) return( false );
    return( ! operator<( num ) );
  }

  /*
    Rational::calcNumerator : 分子の加算・減算処理 ( operator+=,-=用 )

    num : 演算対象の有理数
    cmpSign : 二数の符号比較結果 ( += の場合異なる場合が true,-= の場合同じ場合が true )
  */
  void Rational::calcNumerator( const Rational& num, bool cmpSign )
  {
    Unsigned gcd = CalcGCD( b_, num.b_ );
    a_ *= num.b_ / gcd;
    Unsigned numA = num.a_ * b_ / gcd;
    if ( cmpSign ) {
      // 減算処理
      if ( a_ < numA ) {
	a_ = numA - a_;
	isPos_ = ! isPos_;
      } else {
	a_ -= numA;
      }
    } else {
      // 加算処理
      a_ += numA;
    }
    b_ *= num.b_ / gcd;

    Reduce( &a_, &b_ );
  }

  /*
    operator+=(-=)(*=)(/=) : 四則演算子の多重定義

    const Rational& num : 演算対象の有理数
  */
  Rational& Rational::operator+=( const Rational& num )
  {
    if ( this == &num ) {
      a_ *= 2;
      Reduce( &a_, &b_ );
      return( *this );
    }

    calcNumerator( num, isPos_ != num.isPos_ );

    return( *this );
  }

  Rational& Rational::operator-=( const Rational& num )
  {
    if ( this == &num ) {
      a_ = 0;
      return( *this );
    }

    calcNumerator( num, isPos_ == num.isPos_ );

    return( *this );
  }

  Rational& Rational::operator*=( const Rational& num )
  {
    a_ *= num.a_;
    b_ *= num.b_;
    isPos_ = ( isPos_ == num.isPos_ );

    Reduce( &a_, &b_ );

    return( *this );
  }

  Rational& Rational::operator/=( const Rational& num )
  {
    if ( this == &num ) {
      a_ = b_ = 1;
      return( *this );
    }
    if ( isNaN() || num.isNaN() ) {
      b_ = 0;
      return( *this );
    }

    a_ *= num.b_;
    b_ *= num.a_;
    isPos_ = ( isPos_ == num.isPos_ );

    Reduce( &a_, &b_ );

    return( *this );
  }
}

/*
  CF_sqrt : 連分数を用いて平方根を求める

  x : 対象の数 (符号なし整数型 UI を想定)

  戻り値 : x の平方根 ( 浮動小数点型 Res を想定)
*/
template< class Res, class UI >
Res CF_sqrt( UI x )
{
  // C = [x]を求める
  UI c;
  for ( c = 1 ; c * c <= x ; ++c );
  --c;
  if ( x == c * c ) return( c ); // ちょうど割り切れた

  Res r( x - c * c ); // [x]とC^2の誤差
  Res c2( 2 * c );    // 2C
  Res num( r / c2 );  // [x] - C^2 / 2C

  // 有効桁数の限界まで計算を繰り返す
  for (;;) {
    Res d( r / ( num + c2 ) );
    if ( std::abs( d - num ) < std::numeric_limits< Res >::epsilon() ) break;
    num = d;
  }

  num += c;

  return( num );
}

/*
  CF_goldnum : 連分数を用いて黄金比 φ を求める

  init : 演算精度 (演算回数 符号なし整数型 UI を想定)

  戻り値 : 黄金比 (浮動小数点型 Res を想定)
*/
template< class Res >
Res CF_goldnum()
{
  Res ans( 1 );

  for ( ; ; ) {
    Res buff( 1 / ( ans + 1 ) );
    if ( std::abs( buff - ans ) <= std::numeric_limits< Res >::epsilon() )
      break;
    ans = buff;
  }

  ans += 1;

  return( ans );
}

/*
  CF_napier : 連分数を用いてネイピア数 e を求める

  init : 演算精度(init*2を初期値として連分数を計算 符号なし整数型 UI を想定)

  戻り値 : ネイピア数 (浮動小数点型 Res を想定)
*/
template< class Res, class UI >
Res CF_napier( UI init )
{
  Res ans = Res();
  for ( init *= 2 ; init >= 2 ; init -= 2 ) {
    ans = 1 / ( ans + 1 );
    ans = 1 / ( ans + init );
    ans = 1 / ( ans + 1 );
  }

  ans += 2;

  return( ans );
}

/*
  CF_pi : 連分数を用いて円周率 π を求める

  init : 演算精度(initを初期値として連分数を計算 符号なし整数型 UI を想定)

  戻り値 : ネイピア数 (浮動小数点型 Res を想定)
*/
template< class Res, class UI >
Res CF_pi( UI init )
{
  Res ans( 2 * init + 1 );
  for ( UI i = init ; i > 0 ; --i )
    ans = ( 2 * i - 1 ) + ( i * i ) / ans;

  ans = 4 / ans;

  return( ans );
}

/*
  Fibonacci : フィボナッチ数列と隣り合う項の比 Fn/Fn-1 を求める

  テンプレート引数 D は比の型を表す。

  cnt : 求める項の数 (符号なし整数型 UI を想定)
*/
template< class D, class UI >
void Fibonacci( UI cnt )
{
  UI f0( 1 );
  UI f1( 1 );

  std::cout << "F1 = " << f0 << std::endl;
  std::cout << "F2 = " << f1 << " : F2/F1 = " << D( f1 ) / D( f0 ) << std::endl;

  for ( UI i( 2 ) ; i < cnt ; ++i ) {
    UI f2( f0 + f1 );
    std::cout << "F" << i + 1 << " = " << f2 << " : F" << i + 1 << "/" << "F" << i << " = " << D( f2 ) / D( f1 ) << std::endl;
    f0 = f1;
    f1 = f2;
  }
}

/*
  TS_napier : テイラー展開を用いてネイピア数 eを求める

  戻り値 : ネイピア数
*/
template< class Res >
Res TS_napier()
{
  Res fact( 1 );
  Res d( 2 ); // n=1まで計算

  for ( Res i = 2 ; fact >= numeric_limits< Res >::epsilon() ; i += 1 ) {
    fact /= i;
    d += fact;
  }

  return( d );
}

/*
  Gregory_pi : グレゴリーの公式を用いて円周率πを求める

  cnt : 演算回数(精度)

  戻り値 : 円周率
*/
template< class Res, class UI >
Res Gregory_pi( UI cnt )
{
  Res d = 1; // n=1まで計算

  for ( UI i = 1 ; i < cnt ; ++i )
    d += ( ( i % 2 != 0 ) ? -1 : 1 ) / Res( i * 2 + 1 );

  return( d * 4 );
}

/*
  Machin_pi : マチンの公式を用いて円周率πを求める

  cnt : 演算回数(精度)

  戻り値 : 円周率
*/
template< class Res, class UI >
Res Machin_pi( UI cnt )
{
  const Res d1( 1.0 / 5.0 );
  const Res d2( 1.0 / 239.0 );

  Res d( d1 * Res( 4 ) - d2 ); // n=1まで計算

  for ( UI i = 1 ; i < cnt ; ++i ) {
    d += std::pow( d1, Res( i * 2 + 1 ) ) * Res( 4 ) * ( ( i % 2 != 0 ) ? -1 : 1 ) / Res( i * 2 + 1 );
    d -= pow( d2, Res( i * 2 + 1 ) ) * ( ( i % 2 != 0 ) ? -1 : 1 ) / Res( i * 2 + 1 );
  }

  return( d * Res( 4 ) );
}

/*
  Zeta_pi : ゼータ関数を用いて円周率πを求める

  cnt : 演算回数(精度)

  戻り値 : 円周率
*/
template< class Res, class UI >
Res Zeta_pi( UI cnt )
{
  Res d = 0;

  for ( UI i = 1 ; i < cnt ; ++i )
    d += 1 / std::pow( Res( i ), Res( 2 ) );

  return( std::pow( d * 6, 0.5 ) );
}

/*
  Wallis_pi : ウォリス積を用いて円周率πを求める

  cnt : 演算回数(精度)

  戻り値 : 円周率
*/
template< class Res, class UI >
Res Wallis_pi( UI cnt )
{
  Res d = 1;

  for ( UI i = 1 ; i < cnt ; ++i ) {
    d *= std::pow( Res( i * 2 ), Res( 2 ) );
    d /= i * 2 + 1;
    d /= i * 2 - 1;
  }

  return( d * 2 );
}

/*
  Newton_PowerRoot : ニュートン - ラフソン法を用いてべき乗根を求める

  c : べき乗根を求める数値
  n : べき数

  戻り値 : べき乗根
*/
template< class Res, class UI >
Res Newton_PowerRoot( UI c, UI n )
{
  Res x( c );

  for ( ; ; ) {
    Res fx( std::pow( x, n ) - c );      // f(x)の値
    Res dfx( std::pow( x, n - 1 ) * n ); // df(x)/dxの値
    if ( std::abs( fx / dfx ) <= numeric_limits< double >::epsilon() )
      break;
    x = x - fx / dfx;
  }

  return( x );
}