big_float.cpp

Go to the documentation of this file.

/*****************************************************************************
    
    big_float.cpp -- Big Float Numbers implementation.

    This file is a part of the Arageli library.

    Copyright (C) Alexander Pshenichnikov, 2005--2006
    Copyright (C) Nikolay Santalov, 2005--2006
    
    University of Nizhni Novgorod, Russia

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

//#if 0 // WARNING! DISABLED WHILE PORTING!

#include "config.hpp"

#include <sstream>

#include "powerest.hpp"
#include "logarithm.hpp"
#include "big_float.hpp"


#pragma warning(disable : 4244) // here: conversion from long double to size_t,
                                // possible loss of data


namespace Arageli
{


//a-la bigarith 
void big_float_warning ( const char *s )
{
  std::cerr << "Big float warning: " << s << "\n";
}

void big_float_error ( const char *s )
{
  std::cerr << "Big float error: " << s << "\n";
}

void big_float_fatal_error(const char *s)
{
  std::cerr << "Big float fatal error: \n " << s << "\n";
  exit(1);
}


#define CHECK_PREC(p)  if ( p < PREC_MIN || p > PREC_MAX ) \
    {\
        std::cout << p;\
        big_float_warning ( "attempt to use illegal precision. PREC_MIN is used as default precision" );\
        p = PREC_MIN;\
    }\
    //else prec = p;
#define CHECK_MODE(m) if ( m < EXACT || m > TO_INF ) \
 {\
     big_float_warning ( "attempt to set illegal rounding mode. TO_NEAREST is used as default rounding mode" );\
     m = TO_NEAREST;\
 }\
    
//returns |n| mod m
unsigned fracsion ( int n, unsigned m )
{ 
    if ( n < 0 ) 
    return (m-((-n)%m))%m; 
    else return (n%m);
}

int get_exponent_of_double(double d)
{
 
 int *pd = (int*)&d;
 pd += sizeof ( double ) / sizeof ( int ) - 1;
 *pd <<= 1;
 *pd += 0x80000000;//? only for 32 bits 
 *pd >>= ( sizeof ( int ) * 0x8 - 0xB );
  return *pd + 1;
}

_Internal::digit * get_mantissa_of_double (double d) //array of digits representing mantissa of d
{
    _Internal::digit *man;
    _Internal::digit *pint;
    man = new _Internal::digit [ 4 ];
    if (!man)
    {
        std::cerr << "Big arith fatal error: \n " << "the heap overflow" << "\n";
        exit ( 1 );
    }
    pint = ( _Internal::digit * ) &d;
    man [ 0 ] = ((*pint) << 16 ) >> 16;
    man [ 1 ] =  ( *pint ) >> 16;
    pint++;
    man [ 2 ] = ((*pint) << 16 ) >> 16;
    man [ 3 ] =  ((( *pint ) << 12 ) >> 28) + 16;
    
    return man;
}

// TODO: Delete this function! It's already implemented (and better).
//swap  a  and b
//inline void swap( big_int &a, big_int &b )
//{
//  a = a + b;
//  b = a - b;
//  a = a - b;
//}

// Unused now 
//inline long get_number_of_digits ( big_int b )
//{
//  return (b.length() - 1 ) / _Internal::bits_per_digit + 1;
//}


// return sign of big_int 
inline int get_sign ( big_int s )
{
    if ( s > 0 ) return 1;
    else return -1;
}

//íîðìàëèçàöèÿ ñ ðàçëè÷íûìè ðåæèìàìè îêðóãëåíèÿ
/*void normalize ( long prec,  int mode, big_int &s, big_int &e )
{
    int is_digit = 0;
    big_int stemp ( s ); 
    
    if ( s == 0 ) 
    {
         e = 0;
        return;
    }//åñëè s=0, òî ïîêàçàòåëü òîæå çàíóëÿåì
    if ( mode == EXACT) return;//íè÷åãî íå äåëàåì

    while ( get_number_of_digits ( s ) < prec )
    {
        s = s << bits_per_digit;// *(max_digit+1);
        e = e - 1;
    }//Äîïèñàëè íóëè â êîíåö
    if ( get_number_of_digits ( s ) == prec ) return;
    while ( get_number_of_digits ( s ) > ( prec + 1 ) )
    {
        stemp = s >> bits_per_digit;
        if ( !( is_digit ) && (stemp << bits_per_digit) != s ) is_digit = 1;
        e = e + 1;
        s = stemp;
    }
    switch (mode)
    {   
        case TO_NEAREST:
                       s = ( s + BASE * get_sign ( s ) / 2 ) >> bits_per_digit;
                       e = e + 1;
                       if ( get_number_of_digits ( s ) > prec) 
                       {
                           s = s >> bits_per_digit;
                           e = e+1;
                       }
                       break;
                     
    case TO_ZERO :
                       s = s >> bits_per_digit;
                       e = e + 1;
                       break;
    case TO_P_INF  :
                       if ( s < 0 ) s = s >> bits_per_digit;
                       else 
                       {
                            if ( is_digit ) s = ( s >> bits_per_digit) + 1;
                            else 
                            {
                                   stemp = s >> bits_per_digit;
                                   if ( ( stemp << bits_per_digit ) != s)  s = stemp + 1;
                                   else  s = s >> bits_per_digit;
                            }
                        }
                       e = e + 1;

                       if ( get_number_of_digits ( s ) > prec ) 
                       {
                           s = s >> bits_per_digit;
                           e = e + 1; 
                       }
                       break;
      case TO_M_INF :
                       if ( s > 0 ) s = s >> bits_per_digit;
                       else 
                       {
                             if ( is_digit ) s = s >> ( bits_per_digit ) - 1;
                             else 
                             {
                                    stemp = s >> bits_per_digit;
                                    if ( ( stemp << bits_per_digit ) != s) s = stemp - 1;
                                    else s = s >> bits_per_digit;
                              }
                       }
                       e=e+1;

                      if ( get_number_of_digits ( s ) > prec ) 
                      {
                          s = s >> bits_per_digit; 
                          e = e + 1;
                      }
                      break;
    case TO_INF     :
                       //s=s/(max_digit+1);
                       if ( is_digit )  s = ( s >> bits_per_digit) + get_sign ( s ); 
                       else 
                       {
                             stemp = s >> bits_per_digit;
                             if ( (stemp << bits_per_digit ) != s) s = stemp + get_sign ( s );
                             else        s=s/(max_digit+1);
                       }
                       e = e + 1;
                       if ( get_number_of_digits ( s ) > prec )
                       {
                           s = s >> bits_per_digit;
                           e = e + 1; 
                       }
                       break;
                     
    default         : break;

    }//switch(mode)
    return;
}*/

void big_float::normalize_1 ( long prec, long mode )
{
    if ( !s.sign() ) 
    {
        e = 0;
        return;
    }//significant s=0, then exponenta = 0 
    
    int s_len = s.length() - prec; 
        
    if ( mode == EXACT && s.length() <= PREC_MAX ||  (!s_len)  ) return;// do nothing - mantissa has apropriate length 
    if ( s_len < 0  ) 
    {
        s = s << -s_len;
        e = e + s_len;
        return;
    }//write zero at the end and return
    if ( mode == EXACT )
    {
        mode = TO_NEAREST; // s.lenght > PREC_MAX 
        big_float_error ( "can't perfom operarion with EXACT rounding mode" );
    }
    int is_digit = 0; // needs then rounds to +-infinity
    is_digit = ( (s >> s_len - 1) << s_len - 1 != s ); 
    s = s >> s_len - 1;
    e = e + s_len - 1; 

    switch (mode)
    {   
        case TO_NEAREST:
                       s = ( s + s.sign() ) >> 1;
                       e = e + 1;
                       if ( s.length() > prec ) 
                       {
                           s = s >> 1;
                           e = e + 1;
                       }
                       break;
                     
    case TO_ZERO :
                       s = s >> 1;
                       e = e + 1;
                       break;
    case TO_P_INF  :
                       if ( s.sign() < 0 ) s = s >> 1;
                       else 
                       {
                            if ( is_digit ) s = ( s >> 1) + 1;
                            else 
                                s = s + 1 >> 1;
                       }
                       e = e + 1;

                       if ( s.length() > prec ) 
                       {
                           s = s >> 1;
                           e = e + 1; 
                       }
                       break;
      case TO_M_INF :
                       if ( s.sign() > 0  ) s = s >> 1;
                       else 
                       {
                            if ( is_digit ) s = ( s >> 1) - 1;
                            else 
                                s = s - 1 >> 1;
                       }
                       e = e + 1;

                       if ( s.length() > prec ) 
                       {
                           s = s >> 1;
                           e = e + 1; 
                       }
                       break;
    case TO_INF     :
                       //s=s/(max_digit+1);
                       if ( is_digit )  s = ( s >> 1 ) + s.sign();
                       else 
                       {
                             if ( s [ 0 ] )
                                 s = ( s >> 1 ) + s.sign();
                             else    s = s >> 1;
                       }
                       e = e + 1;
                       
                       if ( s.length() > prec )
                       {
                           s = s >> 1;
                           e = e + 1; 
                       }
                       break;
                     
    default         : break;

    }//switch(mode)
    return;
}
//returns exponent
big_int big_float::get_exponent( void ) const
{
    return e;
}

//returns mantissa
big_int big_float::get_significant( void ) const
{
    return s;
}

//default constructor
big_float::big_float(void)// : e(0),s(0)
{
    //e = s = 0;
}
//copy constructor
big_float::big_float ( const big_float &b )
{
    e = b.e;
    s = b.s;
}

//constructor (from string to big_float) 
big_float::big_float ( const char *str )
{
    std::istringstream s ( str );
    s >> *this;
}

//constructor  ( from two big_int to big_float )
big_float::big_float ( const big_int &s, const big_int &e )
{
    this->e = e;
    this->s = s;
    normalize_1 ( ); //Realy needed? 
    //Yes, needed for present implementation. Not needed if each big_float number keeps it's own precision and rounding mode 
}

//constructor (from double to big_float) !WARNING! Only for 32 bits digit.
big_float::big_float( double b )
{
    if ( b == 0.0 ) 
    {
        big_float ();
        return;
    }
    _Internal::digit *man = get_mantissa_of_double ( b );

    if ( b > 0.0 )
    {
        s.free_mem_and_alloc_number( 1, man, 4 );
    }
    else 
    {
        s.free_mem_and_alloc_number( -1, man, 4 );
    }
    e = get_exponent_of_double ( b ) - 52;
    normalize_1 (  );
}
big_float::big_float( int i )
{
    *this = big_float( i, 0 );
}

//destructor
big_float::~big_float(){}

/*
    operator=
*/
big_float & big_float::operator = ( const big_float &b  )  
{
    if ( this == &b  ) return *this;
    s = b.s;
    e = b.e;
    return *this;
}

big_float & big_float::operator = ( const double d  )  
{
    if ( d == 0.0 ) 
    {
        s = 0;
        e = 0;
        return *this;
    }
    _Internal::digit *man = get_mantissa_of_double ( d );

    if ( d > 0.0 )
    {
        s.free_mem_and_alloc_number( 1, man, 4 );
    }
    else 
    {
        s.free_mem_and_alloc_number( -1, man, 4 );
    }
    e = get_exponent_of_double ( d ) - 52;
    
        normalize_1 (  );
    return *this;
}

big_float & big_float::operator = ( const int i )
{
    return big_float::operator=( double(i) );
}

big_float & big_float::operator = ( const char *str )
{
    return *this = big_float::big_float( str );
}

/*
*  compare operations 
*/
int cmp ( const big_float & a, const big_float & b ) 
{
    big_float temp = a - b;
    return temp.get_significant().sign();
}

int operator == ( const big_float & a, const big_float & b )
{
    return  cmp ( a, b ) == 0;
}

int operator != ( const big_float & a, const big_float & b )
{
    return cmp ( a, b ) != 0;    
}

int operator > (const big_float & a, const big_float & b)
{
    return cmp ( a, b ) == 1; 
}

int operator >= (const big_float & a, const big_float & b)
{
    return cmp ( a, b ) != -1;
}

int operator < (const big_float & a, const big_float & b)
{
    return cmp ( a, b ) == -1;
}

int operator <= (const big_float & a, const big_float & b)
{
    return cmp ( a, b ) != 1;
}


//set default global bits precision 
void big_float::set_global_precision ( long p )
{
    CHECK_PREC(p)
    prec = p;
}

//set global rounding mode
void big_float::set_round_mode ( int m )
{
    CHECK_MODE(m)
    mode = m;
}

//unary plus
big_float operator + ( const big_float &a )
{
    return a;// normalizaton needed may be 
}

//unary minus
big_float operator - ( const big_float &a )
{
    big_float b;
    b.s = -a.s;
    b.e = a.e;
    return b;// normalizaton needed may be ???
}

//operator + 
big_float operator + (const big_float & b, const big_float & c)
{
    return add ( b, c, big_float::prec, big_float::mode);
}

//operator -
big_float operator - (const big_float & b, const big_float & c)
{
    return sub( b, c, big_float::prec, big_float::mode);
}

//operator *
big_float operator * (const big_float & b, const big_float & c)
{
    return mul( b, c, big_float::prec, big_float::mode);
}

//operator /
big_float operator / (const big_float & b, const big_float & c)
{
    return div ( b, c, big_float::prec, big_float::mode );
}

// adding ( with default precision and default rounding mode )
big_float add(const big_float & b, const big_float & c )
{
    return add ( b, c, big_float::prec, big_float::mode );
}

// subtraction ( with default precision and default rounding mode )
big_float sub (const big_float & b, const big_float & c )
{
    return sub( b, c, big_float::prec, big_float::mode );
}

//multiplying ( with default precision and default rounding mode )
big_float mul( const big_float & b, const big_float & c )
{
    return mul( b, c, big_float::prec, big_float::mode );
}

// division ( with default precision and default rounding mode )
big_float div ( const big_float &b, const big_float & c )
{
    return div (b, c, big_float::prec, big_float::mode );
}

//adding ( with default rounding mode )
big_float add(const big_float & b, const big_float & c, long prec)
{
    return add (b, c, prec, big_float::mode );
}

//subtraction ( with default rounding mode )
big_float sub(const big_float & b, const big_float & c, long prec)
{
    return sub(b,c,prec,big_float::mode);
}

//multiplying ( with default rounding mode )
big_float mul(const big_float & b, const big_float & c, long prec)
{
    return mul(b,c,prec,big_float::mode);
}

//division ( with default rounding mode )
big_float div(const big_float & b, const big_float & c, long prec)
{
    return div(b,c,prec,big_float::mode);
}

//adding
big_float add ( const big_float & b, const big_float & c, long prec, int mode )
{

  CHECK_PREC(prec)//
  CHECK_MODE(mode)//âîîáùå ãîâîðÿ, íàêëàäíî // Êõå... âñòàâòå â ARAGELI_ASSERT_0 // Òàê è ñäåëàåì
  big_int b_e(b.e), c_e(c.e), b_s(b.s), c_s(c.s);

  big_float temp;
  
  // some checks   
  if ( !c_s.length() ) // c == 0, so b + c == b  
  {
    temp = b;
    temp.normalize_1 ( prec, mode );
    return temp;
  }
  if ( !b_s.length() ) // b == 0, so b + c == c  
  {
    temp = c;
    temp.normalize_1 ( prec, mode );
    return temp;
  }
  
  // followed if ( ... )... needed for adding no normalized numbers make equal length
  long delta = b_s.length() - c_s.length();
  
  if ( delta > 0 )
  {
    c_s = c_s << delta;
    c_e = c_e - delta;
  }
  else 
  {
    b_s = b_s << -delta;
    b_e = b_e + delta;
  }
  
  if ( b_e < c_e )  
  { 
      swap ( b_e, c_e );
      swap ( b_s, c_s );
  }
     
  big_int ed = b_e - c_e;//ëèáî òðè ðàçà âû÷èñëÿòü, ëèáî îäèí ðàç ñîõðàíèòü
  
  if ( ed > prec  && mode != EXACT/*+1 */) // ìîæåò ñëó÷èòüñÿ òàê ÷òî b_e - c_e î÷åíü áîëüøîå è mode == EXACT !!! 
    //  óñëîâèå âåðíîå, òîëüêî åñëè mode == TO_NAREST
  {
        temp.s = b_s;//åñëè mode=EXACT, òîãäà, âîîáùå ãîâîðÿ, ýòîé ïðîâåðêè íå äîëæíî áûòü
        temp.e = b_e;
        if ( mode != TO_NEAREST )
        {
            temp.s = ( temp.s << 1 ) + c_s.sign();
            temp.e = temp.e - 1;
        }
  }
  else
  {
      if ( ed.length() > _Internal::bits_per_digit /* PREC_MAX - b_s.length()*/ ) //ed > PREC_MAX and mode == EXACT
            big_float_fatal_error( "can't perfom adding with EXACT rounding mode" ); 
      
      // WARNING! I have commented '/*.to_digit()*/' at the following
      // source line because I think it's needless. (Sergey Lyalin)
      b_s = b_s << ed/*.to_digit ()*/; //ìîæåì ïîòåðÿòü çíà÷åíèå, åñëè mode == EXCACT è  b_e - c_e î÷åíü áîëüøîå
       
       temp.s = c_s + b_s;
       temp.e = c_e;
  }
  temp.normalize_1 ( prec, mode );
  return temp;
}

//subtraction
big_float sub(const big_float & b, const big_float & c, long prec, int mode)
{
 big_float temp(c);
 temp.s = -temp.s;
 temp = add ( b, temp, prec, mode);
 return temp;
}

//myltiplying
big_float mul( const big_float & b, const big_float & c, long prec, int mode )
{
  CHECK_PREC(prec)
  CHECK_MODE(mode)
  big_float temp;
  temp.s = b.get_significant() * c.get_significant();
  temp.e = b.get_exponent() + c.get_exponent();
  temp.normalize_1 ( prec, mode );
  return temp;
}

//íàñòîÿùåå äåëåíèå
/*big_float div(const big_float & b, const big_float & c, long prec, int mode)
{
  //Çàäàíû 2 ÷èñëà, íàäî ïîëó÷èòü b/c ñ òî÷íîñòüþ prec
  if (c.s==0) return big_float();//ÍÀÄÎ !!!!!!!!!!!!!!!áåñêîíå÷íîñòü ñ íóæíûì çíàêîì
  big_float x,y,z;
    x.s=1;
        int q=get_number_of_digits(c.s);
        int i=0;
    x.e=-c.e-q; //x.e=x.e-q;
        normalize(prec,mode,x.s,x.e);
        y=x;
        normalize(prec,mode,y.s,y.e);
        z=y-x;
        normalize(prec,mode,z.s,z.e);
        //ïîëó÷èì ñíà÷àëà 1/c
        std::cout<<(-c.e-q)<<"\n"<<x<<"\t\t"<<y<<"\n";
        long p=prec;
        while (/*(z.e>(-p-2-b.e-get_number_of_digits(b.s))&&(i<100))//)
    {
                i++;
        y=x;
                x=2*y-c*y*y;
                z=x-y;
                std::cout<<"x="<<x<<"\t\t"<<"y="<<y<<"\t\t";
                std::cout<<"z.e="<<z.e<<std::endl;
                //prec=prec*2;
                //getch();
        }
        prec=p;
        if (get_sign(x.s)!=get_sign(c.s)) x=-x;
        z=b*x; normalize(prec,mode,z.s,z.e);
    return (z);
} */
//division
big_float div ( const big_float & b, const big_float & c, long prec, int mode)
{
  //std::cout << "1: " << b << "\t 2: " << c << "  prec" << prec << std::endl; 
  if ( c.s==0 ) return big_float();//ÍÀÄÎ !!!!!!!!!!!!!!!áåñêîíå÷íîñòü ñ íóæíûì çíàêîì
  big_float x;
  big_int temp = pow2<big_int> ( prec << 1 );
  
  //std::cout << "temp = " << temp << std::endl;
  //for ( long i=0; i<prec; i++, temp=temp*(max_digit+1),temp=temp*(max_digit+1));
  //std::cout << "c.s = " << c.s << std::endl;
  x.s = temp / c.s;
  x.e = -c.e - prec - prec;
  //x.out( std::cout, 'd');
  x.normalize_1 (prec,mode );
  //std::cout << "  / : "<< x * b << std::endl;
  return mul( x, b, prec, mode );
}
big_float divnu ( const big_float & b, const big_float & c, long prec, int mode)
{
    //Ïîëó÷èì ñíà÷àëà 1/ñ, à ïîòîì óìíîæèì íà b
    //äëÿ ýòîãî, ÷òîáû ñêîðîñòü ñõîäèìîñòè áûëà õîðîøåé ñ äîëæíî áûòü ìåæäó 1/2 è 1
    short k = 0;
    big_float x;
    big_float temp = 1;
    temp.out ( std::cout, 'd');
    x.s = c.s;
    while ( x.s.length() % _Internal::bits_per_digit > 0 ) 
    {
        x.s = x.s << 1;
        k++;
    }
    x.e = -( int )x.s.length();// / bits_per_digit;
    for ( int i=0; i < prec + 1; i++ )
    {
        temp = mul (2,temp,prec,mode) - mul ( mul ( temp,temp,prec,mode ),x,prec,mode );
    }
    for (int i = 0; i < k; i++)
    {
        temp = temp * 2;
    }
    temp.e = temp.e + x.e - c.e;
    temp = mul ( b, temp, prec, mode );
    return temp;
}

big_float div_i ( const big_int & c )
{
    std::size_t l  = c.length();

    std::size_t step = 1;
    big_float res;
        
    res.s = 1;
    
    for ( step = 1 ; step < 32; step <<= 1 )
    {
        res.s = res.s * ( pow2<big_int> ( l * step + 1  ) - c * res.s ); 
    }
    //res.e =- l * ( step );
    std::cout << 10000 * res.s / pow2<big_int> ( 500)<< std::endl;
    //std::cout << res.s << std::endl << res.e << std::endl;
    return res;
}


/*std::ostream & operator << (std::ostream & s, const big_float & x)//ïðîñòàÿ âåðñèÿ
{
 big_int temp=1;
 big_int i=0;

 if ( x.get_exponent() >= 0 )//åñëè ÷èñëî öåëîå, òî è âûâîäèì êàê öåëîå
 {
  for( ; i < x.get_exponent () ; i = i + 1, temp = temp >>  bits_per_digit);
  s << ( x.get_significant( ) * temp ) << std::endl;
 }
 else//äåëåíèå ñòîëáèêîì
 {
  big_int num,denum=1;
  num=x.get_significant()*get_sign(x.s);
  for(i=0;i<-x.get_exponent();i=i+1,denum=denum*(max_digit+1));
  temp=num/denum;
  if(x.get_significant()<0) std::cout<<'-';
  s<<temp<<',';
  num=(num-temp*denum)*10;
  for(i=0;i<PREC;i=i+1)
  {
   temp=num/denum;
   s << temp;
   num=(num-temp*denum)*10;

  }
 }
 return s;
}*/

/*
* stream operators 
*/

void reduce_final_zeros ( std::basic_string<char> &str )
{
    std::size_t last_zero = str.length() - 1;  
    
    while ( str[ last_zero ] == '0' )
        last_zero--;
    str.erase(last_zero + 1);
    if ( str[last_zero] == '.' )
        str.erase(last_zero--);
    if ( str[last_zero] == '.' )
        str.erase(last_zero);
}

std::istream & operator >> ( std::istream &is , big_float &fnum )
{
    char simbol;
    big_int dm, de, bm,be;
    long add = 0;
    
    std::stringstream bufferm; //integer part of mantissa
    std::stringstream buffere; //exponenta
        
    do 
    {
        is.get( simbol );
    } while ( isspace ( simbol ) && !is.fail() );
    
    if ( simbol == '-' || simbol == '+' ) // signum may be 
    {
        bufferm << simbol;
        is.get ( simbol );
    }

    if ( isdigit ( simbol ) ) // mantissa
    {
        bufferm << simbol;
        is.get ( simbol );
        while ( isdigit ( simbol ) && !is.fail() )
        {
            bufferm << simbol;
            is.get ( simbol );
        }
    }
    
    if ( simbol == '.' )//decimal point
    {
        simbol = is.get ();
        while ( isdigit (simbol) && !is.fail() )
        {
            bufferm << simbol;
            simbol = is.get();
            add++;
        }
    }
    
    if ( simbol == 'e' || simbol == 'E')//exponenta
    {
        is.get ( simbol );
        if ( simbol == '-' || simbol == '+') // signum may be 
        {
            buffere << simbol;
            is.get ( simbol );
        }
        while ( isdigit (simbol)  && !is.fail() )
        {
            buffere << simbol;
            is.get ( simbol );
        }
    }
    
    if ( bufferm.str().empty() || bufferm.str() == "+" || bufferm.str() == "-")
    {
        fnum = big_float ();//zero 
    }
    else 
    {
        bufferm >> dm;
        if ( !bufferm.str().empty() && buffere.str() != "+" && buffere.str() != "-" )
            buffere >> de;
        de = de - add;
        do_bin_convert (dm, de, big_float :: prec /*- 1 bits_per_digit - 1*/, bm, be);
        //if ( be.get_sign() != -1 )
        //{
        //  bm = bm << (be - (be / bits_per_digit) * bits_per_digit); 
        //  be =  be  / bits_per_digit;
        //}
        //else 
        //{
        //  bm = bm << bits_per_digit + be - (be / bits_per_digit) * bits_per_digit; 
        //  be = be / bits_per_digit - 1;
        //}
            //std::cout << "bm = " << bm.length() << "    be = "  << be.length() << std::endl; 
        fnum = big_float ( bm, be );
        //std::cout << fnum.get_significant()<< std::endl;
        //std::cout << fnum.get_exponent()<< std::endl;
    }
    
    is.putback ( simbol );
    return is;
}

std::ostream & operator << ( std::ostream &os , const big_float &fnum )
{
    long flags = os.flags( );

    if ( flags & std::ios::scientific )
        fnum.out ( os, 'e' );
    else if ( flags & std::ios :: fixed )
        fnum.out ( os, 'f');
    else 
        fnum.out(os, 'a');
    return os;
}

void big_float :: in ( std::istream & ins, char c )
{
    switch  ( c )
    {
        case 'f': 
            ins >> *this;
            break;      
        case 'd':
            ins >> s >> e;
            normalize_1 ();
            break;
        default:
            big_float_warning ( "unrecognized input format: " );
            std::cerr << c << '\n';
            ins >> *this;
    }
}

void big_float :: out ( std::ostream & os, char c  ) const 
{
    using namespace _Internal;
    
    if ( !s.sign() ) // always show 0.0 if number is zero 
    {
        os << "0.0";
        return;
    }
    
    std::stringstream str;            // for temporary output 
    big_int dm, de;                   // mantissa and exponenta in decimal representation   
    int  position;                    // decimal point position 
    std::basic_string<char> temp_str; // for temporary output

    if ( os.precision() > D_PREC_MAX )
        os.precision( D_PREC_MAX ); 
    std::ios_base::fmtflags prev_flags = os.flags ();   // ostream flags
    
    //don't take a lot of trouble over signums
    if ( s.sign() == -1 ) 
        os << "-";
            
    else if ( prev_flags & std::ios::showpos )
        os << "+";

    str.setf( std::ios::dec | std::ios::showpos ); // always out signum see below  
    
    switch ( c )
    {
        case 'f': //floating point mode e.g 123.456 
                position = s.length() + e.to_native_int<int>(); // Check this conversion good way to improve this function 
                
                //approximately evaluate floating point position in decimal representation  
                if ( position >= 0 )
                    position = 1 + ( int )((position + 1) * log (2.0l) / log ( 10.0l )); 
                else 
                    position = (int)( (long double) position * log (2.0l) / log ( 10.0l ) );
                
                if ( position + os.precision() > 0 ) //there is significant digits (probably) 
                {
                    // convert to decimal and out in string stream  
                    do_dec_convert(dm, de, os.precision() + position, s, e);
                    str << dm;
                    temp_str = str.str();// save dm in non-const string variable for manipulation 
                    temp_str.erase(0,1); // erase signum ( it have been printed already )  

                    if ( position > 0 )
                    {
                        //insert decimal point in right position 
                        temp_str.insert( position = temp_str.length() + de.to_native_int<int>(), "." );
                    }
                    else //there is leading zeros 
                    {
                        os << "0.";
                        if ( position ) //one zero is produced already 
                        {
                            os.width( -position );
                            os.fill ( '0' );
                            os << "0"; 
                        }
                        if (temp_str.length() <= os.precision() + position ) 
                            position--, os << "0"; 
                    }

                    //erase excess digits and final zeros //without rounding!!  
                    temp_str.erase( position + os.precision() ); 
                    reduce_final_zeros ( temp_str );
                    os << temp_str;
                }
                else // no significant digits on this precision (only zeros) 
                    os << "0.0";
                            
                break;

        case 'e':// exponential mode  e.g 1.23e456
                
                do_dec_convert ( dm, de, os.precision(), s, e );
                str << dm;
                temp_str = str.str();
                temp_str.erase(0,1); // erase signum ( it have been printed already )
                temp_str.insert(1,".");
                reduce_final_zeros ( temp_str ); 

                //out 
                if ( prev_flags & std::ios::uppercase )
                    os << temp_str << 'E' << de + os.precision();
                else 
                    os << temp_str << 'e' << de + os.precision();   

                break; 

        case 'a': // auto format ( if neither scientific notation and not fixed notation are used )
                do_dec_convert ( dm, de, os.precision(), s, e );
                str << dm;
                temp_str = str.str();
                temp_str.erase(0,1);// erase signum ( it have been printed already )

                if ( de.sign() < 0 && de.to_native_int<int>() > -4 - os.precision () ) // show number in fixed format 
                {
                    position = temp_str.length()+ de.to_native_int<int>(); //determine floating point position  
                    if ( position <= 0 ) //there is leading zeros
                    {
                        os << "0.";
                        if ( position ) //one zero is produced already 
                        {
                            os.width( -position );
                            os.fill ( '0' );
                            os << "0"; 
                        }
                    }
                    else if ( position != os.precision() ) 
                        //insert decimal point in right position 
                        temp_str.insert ( position, "." );
                    else 
                        temp_str.insert ( position, ".." ); //will be erased in followed

                    //erase excess digits and //without rounding!!  
                    temp_str.erase( temp_str.length()-1 );
                    reduce_final_zeros ( temp_str );
                    os << temp_str;
                }
                else // use scientific format
                {
                    temp_str.insert( 1,"." );
                    //erase excess digits //without rounding!!  
                    temp_str.erase( temp_str.length() - 1 ); 
                    reduce_final_zeros ( temp_str );

                    //out 
                    if ( prev_flags & std::ios::uppercase )
                        os << temp_str << 'E' << de + os.precision();
                    else 
                        os << temp_str << 'e' << de + os.precision(); 
                }
                
            break;

        case 'b'://binary ( zero and one ). out significant and exponent separately 
            if ( s.sign() == -1 )
                os << '-';
            else if ( prev_flags & std::ios::showpos )
                os << '+';
            for ( int counter = s.length () - 1; counter >= 0 ; counter-- )
                os << s [ counter ]; 

            os << '\t';

            if ( e.sign() == -1 )
                os << '-';
            else if ( prev_flags & std::ios::showpos )
                os << '+';

            for ( int counter = e.length () - 1; counter >= 0 ; counter-- )
                os << e [ counter ]; 
            break;

        case 'd': // out significant and exponent separately, in decimal mode   
            os.setf( std::ios :: dec, std::ios::hex );
            os << s << '\t' << e;
            break;

        case 'h': // out significant and exponent separately, in hex mode 
            os.setf( std::ios::hex, std::ios::dec );
            os << s << '\t' << e;
            break;

        case 'o': // out significant and exponent separately, in oct mode 
            os.setf( std::ios :: oct, std::ios::dec | std::ios::hex );
            os << s << '\t' << e;
            break;

        default:
            std::cerr << "Warning!: unrecognized output format:" << c << std::endl;
            os.setf( std::ios::dec);
            os << s << '\t' << e;
    }
    //restore flags 
    os.setf( prev_flags, os.flags() );

}
//returns a floating-point value representing the smallest integer that is greater than or equal to a
big_float ceil  (const big_float & a)
{
    big_float res ( a );
    big_int e = a.e;
    big_int temp;

    if ( e >=  0 )
        return res;
    if ( e <= -big_int(a.s.length()) ) 
    {
        if ( a.s > 0 )
            return 1.0;
        else return 0.0;
    }
    temp = ( res.s >> e) << e;
    if ( res.s != temp && res.s > 0 ) 
        res.s = temp + pow2<big_int>( e/*.to_digit ()*/ );  
    else res.s = temp;
    return res;
}

// returns a floating-point value representing the largest integer that is less than or equal to a
big_float floor (const big_float & a)
{
    big_float res ( a );
    big_int e = a.e;
    big_int temp;

    if ( e >=  0 )
        return res;
    if ( e <= -big_int(a.s.length()) ) 
    {
        if ( a.s > 0 )
            return 1.0;
        else return 0.0;
    }
    temp = ( res.s >> e) << e;
    if ( res.s != temp && res.s < 0 ) 
        res.s = temp - pow2<big_int>( e );  
    else res.s = temp;
    return res;

}

// Returns the fractal part of the number
big_float frac  (const big_float & a)
{
    if ( a.e >= 0 )
        return 0.0;
    
    big_float res ( a );
    if ( res.e <= - big_int(res.s.length()) )
        return res;
    
    res.s = res.s - ((res.s >> res.e) << res.e);
    return res;
}

// returns an big-integer value representing the smallest integer that is greater than or equal to a
big_int iceil (const big_float & a)
{
    std::size_t l = a.s.length();
    if ( !l ) return l;
    
    big_int temp ( a.s );
    if ( a.e > 0 )
    {
        if ( a.e + a.s.length() > SHRT_MAX/*??*/  )
            big_float_fatal_error ( "Can't convert from big_float to big_int ( exponent too large )..." );
        return a.s <<a.e;   
    }
    if ( a.e <= -big_int(a.s.length()) )
        return ( a.s > 0 ) ? 1 : 0;
    temp = temp >> a.e;
    if ( temp << a.e != a.s && temp > 0)
        return temp + 1;  
    return temp;
}

// returns a big-integer  value representing the largest integer that is less than or equal to a
big_int ifloor(const big_float & a)
{
    std::size_t l = a.s.length();
    if ( !l ) return l;
    
    big_int temp ( a.s );
    if ( a.e > 0 )
    {
        if ( a.e + a.s.length() > SHRT_MAX/*??*/  )
            big_float_fatal_error ( "Can't convert from big_float to big_int ( exponent too large )..." );
        return a.s <<a.e;   
    }
    if ( a.e <= -big_int(a.s.length()) )
            return ( a.s > 0 ) ? 0 : -1;
    temp = temp >> a.e;
    if ( temp << a.e != a.s && temp > 0)
        return temp;  
    return temp - 1;
}
//sqrt 
big_float fsqrt( const big_float & b )
{
    return fsqrt( b, big_float :: prec, big_float :: mode);
}

//user must call fsqrt( const big_float & b ) function instead two forward
//functions ( because EXACT mode not supported yet )  
big_float fsqrt(const big_float & bf, long prec)
{
    return fsqrt( bf, prec, big_float :: mode);
}

big_float fsqrt(const big_float & bf, long prec, int mode)
{
    big_float res ( bf );
    int l;

    if ( res.e [ 0 ] )
    {   
        res.s = res.s << 1;
        res.e = res.e - 1;
    }
    
    l = (res.s.length() - 2 * prec - 1);
    l = (l / 2) * 2;

    if ( l > 0 )
    {   
        res.s = res.s >> l;
        res.e = res.e + l;
    }
    else 
    {
        res.s = res.s << (-l) + 2;
        res.e = res.e + l - 2;
    }
    res.s = sqrt ( res.s );
    res.e = res.e >> 1;
    res.normalize_1 ( prec, mode );//error may be if mode == TO_INF
    
    return res;
}
//Newton fsqrt
big_float nfsqrt ( const big_float & bf, long prec, int mode )
{
    big_float res;
    //std::cout << bf << std::endl;
    res.e = ((( bf.get_exponent() ) + bf.get_significant().length() - 1) >> 1) + 1;
    res.s = 1;
    //res.normalize_1();
    std::cout << "x0 "<<res << std::endl; 
    //res.out ( std::cout, 'd' ) ;
    //std::cout << std::endl;
    std::size_t n = log ( (long double)prec + 1 ) / log ( 2.0l ) + 1;
    
    for ( std::size_t counter = 1; counter <= 2 * n; counter ++ )
    {
        std::cout << div( bf, res, counter * 2, EXACT ) << std::endl;
        res = add ( res,  div( bf, res, counter * 2, EXACT ), counter * 2, TO_NEAREST );  
        std::cout << res << std::endl;
        res.e = res.e - 1;
        std::cout << counter << " iter \t";
        std::cout << res << std::endl;
        //res.out( std::cout, 'd' );
        //std::cout << std::endl;
    }
    res.normalize_1 ( prec, mode );//error may be if mode == TO_INF
    
    return res;
}
//Random numbers 
big_float frandom  ( long prec )
{
    CHECK_PREC(prec)//
    return big_float ( big_int::random_with_length_or_less( prec ), -prec );
}
//Random numbers 
big_float frandom  ( )
{
    return big_float ( big_int::random_with_length_or_less( big_float::prec ), -big_float::prec );
}

big_float frandom1 ( long bits, const big_int & exp)
{
    CHECK_PREC(bits)//
    return big_float (  big_int::random_with_length_or_less( bits ) + pow2<big_int>( bits - 1 ) , exp );
}

/*not used now*/
int big_float::is_NaN( )
{
  return 0;
}     
int big_float::is_Zero( )
{
 return ( s==0 );
}

int big_float::is_pZero( )
{
 return ( ( s ==0 ) && (get_sign( s ) == 1) );
}

int big_float::is_mZero()
{
 return ((s==0)&&(get_sign(s)==-1));
}

int big_float::is_Special()
{
  return ((is_NaN())||(is_Inf())||(is_Zero()));
}     
   

int big_float::is_Inf()
{
  return 0;
}

long big_float::prec; //bits precision 
long big_float::mode; //rounding mode


}

//#endif

---