gistfile1.cpp

/*
 * main.cpp
 * Date: 2012-10-31
 * Author: Karsten Ahnert (karsten.ahnert@gmx.de)
 */

#include <cstdlib>
#include <cmath>
#include <iostream>
#include <iterator>
#include <utility>
#include <algorithm>
#include <cassert>
#include <vector>
#include <complex>

extern "C" {
#include <quadmath.h>
}

const __float128 zero =strtoflt128 ("0.0", NULL);

namespace std {

    inline __float128 abs( __float128 x )
    {
        return fabsq( x );
    }

    inline __float128 sqrt( __float128 x )
    {
        return sqrtq( x );
    }

    inline __float128 pow( __float128 x , __float128 y )
    {
        return powq( x , y );
    }

    inline __float128 abs( std::complex< __float128 > x )
    {
        return sqrtq( x.real() * x.real() + x.imag() * x.imag() );
    }

    inline std::complex< __float128 > pow( std::complex< __float128> x , __float128 y )
    {
        __float128 r = pow( abs(x) , y );
        __float128 phi = atanq( x.imag() / x.real() );
        return std::complex< __float128 >( r * cosq( y * phi ) , r * sinq( y * phi ) );
    }
}

#include <boost/array.hpp>
#include <boost/range/algorithm.hpp>
#include <boost/range/adaptor/filtered.hpp>
#include <boost/range/numeric.hpp>
#include <boost/numeric/odeint.hpp>



using namespace boost::numeric::odeint;
using namespace std;

typedef __float128  my_float;
typedef std::vector< std::complex < my_float > > state_type;

struct radMod
{
    my_float m_om;
    my_float m_l;

    radMod( my_float om  , my_float l )
        : m_om( om ) , m_l( l ) { }

    void operator()( const state_type &x , state_type &dxdt , my_float r ) const
    {      

        dxdt[0] = x[1];
        dxdt[1] = -(2*(r-1)/(r*(r-2)))*x[1]-((m_om*m_om*r*r/((r-2)*(r-2)))-(m_l*(m_l+1)/(r*(r-2))))*x[0];
    }
};


inline ostream& operator<< (ostream& os, const __float128& f) {

    char* y = new char[1000];
    quadmath_snprintf(y, 1000, "%.30Qg", f) ;
    os.precision(30); 
    os<<y;
    delete[] y;
    return os;
}


struct streaming_observer
{
    std::ostream& m_out;

    streaming_observer( std::ostream &out ) : m_out( out ) { }

    template < class State >
    void operator()( const State &x , my_float t ) const
    {
        m_out << t <<", ("<< x[0].real() <<", "<< x[0].imag()<<"I )"<<endl;
    }
};


struct push_back_state_and_time
{
    std::vector< state_type > & m_states;
    std::vector< my_float >   & m_times;

    push_back_state_and_time( std::vector< state_type > &states , std::vector< my_float > &times )
        : m_states( states ) , m_times( times ) { }

    void operator()( const state_type &x , my_float t )
    {
        m_states.push_back( x );
        m_times.push_back( t );
    }
};





int main( int argc , char **argv )
{
    __float128 x1 = 8.0 ;
    __float128 x2 = 9.0;
    __float128 x3 = -9.0;
    __float128 x4 = pow( x1 , x2 ); 
    __float128 x5 = abs( x3 );

    cout<< x1 <<endl;
    cout<< x2 <<endl;
    cout<< x3 <<endl;
    cout<< x4 <<endl;
    cout<< x5 <<endl;




    state_type x(2);

    my_float re0 = strtoflt128 ("-0.00008944230755601224204687038354994353820468", NULL);
    my_float im0 = strtoflt128 ("0.00004472229441850588228136889483397204368247", NULL);
    my_float re1 = strtoflt128 ("-4.464175354293244250869336196695966076150E-6 ", NULL);
    my_float im1 = strtoflt128 ("-8.950483248390306670770345406051469584488E-6", NULL);

    x[0] = complex< my_float >( re0 ,im0 );
    x[1] = complex< my_float >( re1 ,im1 );

    __float128 x6 = abs( x[0] );
    complex< my_float > x7 = pow( x[0] , 1 );
    cout << x6 << endl;
    cout << x7.real() << ", " << x7.imag() << endl;


    const my_float dt =strtoflt128 ("-0.001", NULL);
    const my_float start =strtoflt128 ("10000.0", NULL);
    const my_float end =strtoflt128 ("9990.0", NULL);
    const my_float omega =strtoflt128 ("2.0", NULL);
    const my_float ell =strtoflt128 ("1.0", NULL);


    // //[***************basic RK4 with float128 with streaming***************
    // typedef runge_kutta4< state_type , my_float > stepper_type;
    // integrate_const( stepper_type() , radMod(omega, ell ) , x , start , end, dt  , streaming_observer(cout) );
    // //]

    //[***************basic RK4 with float128 with pushback***************
    // std:: vector<state_type> x_vec;   
    // std:: vector<my_float> times;    
    // typedef runge_kutta4< state_type , my_float > stepper_type;
    // size_t steps=integrate_const( stepper_type() , radMod(omega, ell ) , x , start , end, dt  , push_back_state_and_time( x_vec , times ) );
    // for( size_t i=0; i<=steps; i++ )
    // {
    //     cout << times[i] << '\t' <<" ("<< x_vec[i][0].real() << ", " << x_vec[i][0].imag()<<"I) " << endl;     
    // }
    //]

//[*************testing iteration of these objects****************************
//std::vector<std::complex < my_float > > myFloatVector;
//std::vector<std::complex < my_float > >::iterator myFloatIterator;    
//myFloatVector.push_back(complex<my_float>(re0,im0));
//myFloatVector.push_back(complex<my_float>(re1,im1));
// myFloatVector.push_back(complex<my_float>(start,end));
//for(myFloatIterator = myFloatVector.begin(); 
//   myFloatIterator != myFloatVector.end();
//   myFloatIterator++)
//  {
// cout<<"( "<<(*myFloatIterator).real()<<" , "<<(*myFloatIterator).imag()<<" I)"<<endl;
//   }
//]

//[*************basic RK4 with accumulator(doesn't work because it can't match accumulate call************
/*
  runge_kutta4< state_type, my_float > stepper;
  my_float res = std::accumulate( make_const_step_iterator_begin( stepper , radMod(omega , ell) , x , start , end , dt ) ,
  make_const_step_iterator_end( stepper , radMod(omega , ell) , x ) ,
  zero ,
  []( my_float sum , const state_type &x ) {
  return sum + x[0].real(); } );    
  cout << res << endl;
*/

    //]

    //[**************adaptive stepper float128*********************************
    my_float abs_err=strtoflt128("1.0E-15", NULL), rel_err=strtoflt128("1.0E-10", NULL);
    my_float  a_x=strtoflt128("1.0", NULL) ,a_dxdt = strtoflt128("1.0", NULL);  
    typedef runge_kutta_dopri5< state_type, my_float > dopri5_type;
    typedef controlled_runge_kutta< dopri5_type > controlled_dopri5_type;
    typedef dense_output_runge_kutta< controlled_dopri5_type > dense_output_dopri5_type;
    
    dense_output_dopri5_type dopri5( controlled_dopri5_type( default_error_checker< my_float >( abs_err , rel_err , a_x , a_dxdt )  ) );

    std::for_each( make_adaptive_time_iterator_begin(dopri5 , radMod(omega , ell) , x , start , end , dt) ,
                   make_adaptive_time_iterator_end(dopri5 , radMod(omega , ell) , x ) , []( const std::pair< state_type&, my_float > &x ) {
                       std::cout << x.second << ", " << x.first[0].real() << "\n"; }
        );
    //]

//[***************Boost::const_step_iterator_accumulate_range(exactly same issues it seems-no matching call)***********
    /*
      runge_kutta4< state_type > stepper;

      my_float res = boost::accumulate( make_const_step_range( stepper , radMod(omega , ell) , x , start , end , dt ) , zero ,
      []( my_float sum , const state_type &x ) {
      return sum + x[0]; } );
      cout << res << endl;
    */
    //]



    return 0;
}