// f107_o.h    -    float_107 type for C++
//
//    Copyright (C) 1998-2013 Robert P. Munafo.
//
//    This program is free software; you can redistribute it and/or
//    modify it under the terms of the GNU General Public License as
//    published by the Free Software Foundation; either version 2 of the
//    License, or (at your option) any later version.
//
//    This program is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
//    General Public License (immediately following this text) for more
//    details.
//
//    To receive a copy of the GNU General Public License, write to the
//    Free Software Foundation, Inc., 59 Temple Place, Suite 330,
//    Boston, MA 02111-1307 USA
//
// This is a quadruple-precision float type built upon IEEE 64-bit doubles.
//
// Usage instructions and other information are in f107_o.cpp


#ifndef _f107_o_
#define _f107_o_

/* The following sets the default value of F107_FMA for cases where the
   makefile did not specify it on the compile line, or in a project header,
   etc.

   F107_FMA should be set to 1 if you are building for a target that correctly
   implements the Fused Multiply-Add instructions. This includes all PowerPC
   targets, PowerPC emulators that CORRECTLY implement FMA, and the forthcoming
   Sandy Bridge architecture from Intel. Most emulators, including Sheepshaver
   and Rosetta, do not correctly emulate FMA.

   Note that "#pragma fp_contract on" is not implemented by the gcc compiler
   and generates an appropriate warning. The #ifndef __GNUC__ here is to avoid
   getting the warnings. When compiling with gcc (e.g. in XCode on MacOS X, or
   anything under Linux, etc.) in addition to specifying F107_FMA in your
   makefile or compile line, you should also tell the compiler whether to
   generate FMA's using the -mfused-madd or -mno-fused-madd options.
 */
#ifdef __POWERPC__
# ifndef __GNUC__
#   ifndef F107_FMA
#     define F107_FMA 1
#   endif
# endif
#else
# ifdef __FMA__
#   ifdef __AVX2__
#     include <immintrin.h>
#     define F107_FMA 1
#   endif
# endif
#endif

/* Now make sure F107_FMA has some numeric value */
#ifndef F107_FMA
# define F107_FMA 0
#endif

#ifndef F107_APPROX
# define F107_APPROX 1
#endif

// "pf" stands for "primitive float"; it is the floating point type upon
// which the library is built.
typedef double f107_pf;
typedef unsigned char f107_u8;
typedef short f107_s16;
typedef unsigned int f107_u32;

#define F107_WORD_0(x) (*(   (f107_u32 *)&(x) ))
#define F107_WORD_1(x) (*(1+((f107_u32 *)&(x))))

extern const f107_pf K107_p27_1;
extern int f107_need_init;

class f107_o;

void f107_fatal_error(const char * s, f107_s16 err, f107_s16 param);
void f107_o_init(void);
int f107_isinf(f107_pf x);
int f107_isinf_raw(f107_pf x);
int f107_isnan(f107_pf x);
int f107_isnan_raw(f107_pf x);
void f107_stats(void);
void print_f107_a(f107_o a, int cr);

void div_f107(f107_o a, f107_o b, f107_o &res);
f107_o square ( const f107_o z );
f107_o fabs ( const f107_o z ); 
f107_o trunc ( const f107_o z );
f107_o floor(const f107_o x);

void pr_raw_double_107(double d, int sexp, int hm, int mant);
void pr_mant_comp(double mc);
void cfa(f107_o x, f107_s16 order, f107_o *n, f107_o *d);
f107_o f107_intpow(f107_o x, int p);
f107_o sqrt(f107_o x);
f107_o f107_sincos(f107_o xq, f107_o zq, f107_o uq, double i);
f107_o sin(f107_o x);
f107_o cos(f107_o x);
f107_pf f107_pow2_int(int x);
f107_o restrict_tanh(f107_o x);
f107_o exp(f107_o x);
f107_o log(f107_o target);
f107_o sinh(f107_o x);
f107_o tanh(f107_o x);

void f107_to_digits(f107_o x, char *s, int * expn, int * sign, int precision);
void f107_spf(char *s1, char sign_flag, int precision, f107_o x);
void f107_spfg(char *s1, int length, char sign_flag, int precision, f107_o x);
int f107_strlen(char * s);
void f107_strncpy(char * to, char * fr, int n);
void f107_snprinf_int(char * to, int len, int x);
f107_o f107_sscan(char * s);



class f107_o
{
  private:
  
  public:
    f107_pf  c0;
    f107_pf  c1;


    void print_f107_a(f107_o a, int cr);

    void div_f107(f107_o numer, f107_o denom, f107_o &res);

    // constructors

    // initialize
    //  input: ()
    // output: f107 a = 0
    //    ops: 2
    f107_o ( ) : c0(0), c1(0) { }

    // convert_12
    //  input: f53 x
    // output: f107 x = x
    //    ops: 2
    f107_o ( f107_pf x ) : c0(x), c1(0.0) { }

    // convert_22
    //  input: f107 x
    // output: f107 x = x
    //    ops: 2
    f107_o ( const f107_o & x ): c0(x.c0), c1(x.c1) { }

    //                                            f107_o x = "1.23";
    f107_o ( const char * s ) : c0(0), c1(0)
    {
      // printf("Constructing f107 from string '%s'\n", s);
      *this = f107_sscan((char *) s);      
    }

    // Conversion operators (demotion)
    operator short ( )
    {
      short  result;

      result = (short) c0;
      return result;
    }

    operator int ( )
    {
      long  result;

      // assumes that int has no greater precision than f107_pf
      result = (long) c0;
      return result;
    }

    operator long ( )
    {
      long  result;

      // long might have greater precision than f107_pf (e.g. Intel Core
      // 2 Duo, long is 64 bits)
      result = ((long) c0) + ((long) c1);
      return result;
    }

    operator float ( )
    {
      float  result;

      result = c0;
      return result;
    }

    operator double ( )
    {
      double  result;

      result = c0;
      return result;
    }

    operator long double ( )
    {
      long double  result;

      result =  ((long double) c0) + ((long double) c1);
      return result;
    }

    // assignment

    // operator =
    //  input: f53 a
    // output: f107 a = a
    //    ops: 2
    f107_o & operator = ( f107_pf x)
    {
      c0 = x;
      c1 = 0;
      
      return *this;
    }

    // operator =
    //  input: f107 a
    // output: a
    //    ops: 2
    f107_o & operator = ( const f107_o x)
    {
      c0 = x.c0;
      c1 = x.c1;

      return *this;
    }

    // operator =
    //  input: f53 a
    // output: f107 a = a
    //    ops: 2
    f107_o & operator = ( const char * s)
    {
      *this = f107_sscan((char *) s);
      
      return *this;
    }


    // Very simple calculations and building blocks for addition

    // native multiply by 2

    // f107_native_2x
    //
    // scale2_22
    //  input: f107 x
    // output: x * 2.0
    //    ops: 2
    static f107_o f107_native_2x ( const f107_o x )
    {
      f107_o res;

      res.c0 = x.c0 * 2.0;
      res.c1 = x.c1 * 2.0;

      return(res);
    }

    // f107_native_half
    //
    // scale2_22
    //  input: f107 x
    // output: x * 0.5
    //    ops: 2
    static f107_o f107_native_half ( const f107_o x )
    {
      f107_o res;

      res.c0 = x.c0 * 0.5;
      res.c1 = x.c1 * 0.5;

      return(res);
    }

    // f107_native_scale
    //
    // scale2_22
    //  input: f107 x
    // output: x * 0.5
    //    ops: 2
    static f107_o f107_native_scale ( const f107_o x, f107_pf sc )
    {
      f107_o res;

      res.c0 = x.c0 * sc;
      res.c1 = x.c1 * sc;

      return(res);
    }

    /* normalize_22 is called QUICK_TWO_SUM in the QD paper (Hida et al. 2000).
       Normalize 2 f107_pf's, given the greater and lesser. (If it is
           unknown which is greater, use add_112)
    */
    // normalize_22
    //  input: f53 a; f53 b
    // output: (s,e) = a + b
    //    ops: 3
    static f107_pf normalize_22(f107_pf a, f107_pf b, f107_pf &e)
    {
      f107_pf s = a+b;
      e = b - (s - a);
      return s;
    }

    /* split_53 is called SPLIT in the QD paper (Hida et al. 2000).
           It divides a f107_pf into two parts with at most 27 bits of
           precision. We treat f107_pf as if it has 54 bits of precision; see
           the note in f107_o.cpp */
    // split_f53
    //  input: f53 a
    // output: hi and lo components, each with 27 bits, such that hi + lo = a
    //    ops: 4
    static f107_pf split_f53(f107_pf a, f107_pf &lo)
    {
      f107_pf t = K107_p27_1*a;
      f107_pf hi = t-(t-a);
      lo=a-hi;
      return hi;
    }

/* On PowerPC, specify that we want to generate FMA (fused
   multiply-add) instructions */
#ifdef __POWERPC__
# ifndef __GNUC__
#   if F107_FMA
#     pragma fp_contract on
#   endif
# endif
#else
# if F107_FMA
    /* NOTE: If this was defined it is because we were compiled
       with the options -mavx2 -mfma, which will only work on
       Intel Haskell and later or comparable AMD processors. */
/* #error "F107_FMA is only implemented for PowerPC" */
# endif
#endif

/* f107_fmsub(a, b, c) returns a*b-c using an FMA operation */
#ifdef __POWERPC__
# if __GNUC__
        /* The inline assembly is from <ppc_intrinsics.h> */
        static inline f107_pf 
        f107_fmsub (f107_pf  a, f107_pf c, f107_pf b) {
          f107_pf result;
          __asm__ ("fmsub %0, %1, %2, %3" : "=f"(result) : "f"(a), "f"(c), "f"(b));
          return result;
        }
# else
        /* Metrowerks generates the instruction "properly" */
        static inline f107_pf 
        f107_fmsub (f107_pf  a, f107_pf c, f107_pf b) { return (a * c) - b; }
# endif
#else
        static inline f107_pf 
        f107_fmsub (f107_pf  a, f107_pf c, f107_pf b) {
# if F107_FMA
          __m128d aa, bb, cc, rr;
          aa = _mm_set_sd(a); /* aa[0]=a, aa[1]=undefined */
          bb = _mm_set_sd(b); /* bb[0]=b, bb[1]=undefined */
          cc = _mm_set_sd(c); /* cc[0]=c, cc[1]=undefined */
          rr = _mm_fmsub_pd(aa,cc,bb);
          return _mm_cvtsd_f64(rr);
# else
          return (a * c) - b;
# endif
        }
#endif

    // mul_112
    //  input: f53 a; f53 b
    // output: (r0, r1) = a * b
    //    ops: 3
    static f107_pf mul_112(f107_pf a, f107_pf b, f107_pf &e)
    {
#if F107_FMA
      f107_pf p;
      p = a * b;
          e = f107_fmsub(a,b,p);
#else
      f107_pf alo, blo, p, ahi, bhi;
      p = a*b;
      ahi = f107_o::split_f53(a, alo);
      bhi = f107_o::split_f53(b, blo);
      e = ((ahi*bhi - p) + ahi*blo + bhi*alo) + alo*blo;
#endif

      return p;
    }

    // square_12
    //  input: f53 a
    // output: (r0, r1) = a * a
    //    ops: 3
    static f107_pf square_12(f107_pf a, f107_pf &e)
    {
#if F107_FMA
      f107_pf p;
          p = a * a;
          e = f107_fmsub(a,a,p);
#else
      f107_pf alo;
      f107_pf p = a*a;
      f107_pf ahi = f107_o::split_f53(a, alo);
      e=((ahi*ahi - p) + 2.0*ahi*alo) + alo*alo;
#endif

      return p;
    }

/* We're done using the FMA instructions; tell the GCC compiler to stop
   generating them */
#ifdef __POWERPC__
# ifndef __GNUC__
#   if F107_FMA
#     pragma fp_contract reset
#   endif
# endif
#endif


    // unary negate

    // neg_2
    //  input: f107 a
    // output: f107 x = -a
    //    ops: 2
    friend f107_o operator - ( const f107_o lhs )
    {
      f107_o  res;

      res.c0 = -(lhs.c0);
      res.c1 = -(lhs.c1);

      return res;
    }


    // operators << and >> are used for "quick" comparisons that
    // look at just the most significant part of the f107.

    // operator << (crude magnitude comparison)
    //  input: f107 a; f107 b
    // output: int (a.c0 < b.c0)
    //    ops: %%     len: 3
    friend int operator << ( const f107_o lhs, const f107_o rhs )
    {
      return(lhs.c0 < rhs.c0);
    }

    // operator << (crude magnitude comparison)
    //  input: f107 a; f53 b
    // output: int (a.c0 < b)
    //    ops: %%      len: 3
    friend int operator << ( const f107_o lhs, f107_pf rhs )
    {
      return(lhs.c0 < rhs);
    }

    // operator << (crude magnitude comparison)
    //  input: f53 a; f107 b
    // output: int (a < b.c0)
    //    ops: %%     len: 3
    friend int operator << ( f107_pf lhs, const f107_o rhs )
    {
      return (lhs < rhs.c0);
    }


    // operator >> (crude magnitude comparison)
    //  input: f107 a; f107 b
    // output: int (a.c0 > b.c0)
    //    ops: %%
    friend int operator >> ( const f107_o lhs, const f107_o rhs )
    {
      return (lhs.c0 > rhs.c0);
    }

    // operator >> (crude magnitude comparison)
    //  input: f53 a; f107 b
    // output: int (a > b.c0)
    //    ops: %%
    friend int operator >> ( f107_pf lhs, const f107_o rhs )
    {
      return (lhs > rhs.c0);
    }

    // operator >> (crude magnitude comparison)
    //  input: f53 a; f107 b
    // output: int (a.c0 > b)
    //    ops: %%
    friend int operator >> ( const f107_o lhs, f107_pf rhs )
    {
      return (lhs.c0 > rhs);
    }


    /* add_112 is called TWO_SUM in the QD paper (Hida et al. 2000).
       Add 2 f107_pf's, returning the 2 components of a normalized dd_real
    */
    // add_112
    //  input: f53 a; f53 b
    // output: (s,e) = a + b
    //    ops: 6
    static f107_pf add_112(f107_pf a, f107_pf b, f107_pf &e)
    {
      f107_pf s = a + b;
      f107_pf v = s - a;
      e = (a-(s-v))+(b-v);
      return s;
    }

    /* sub_112 is called TWO_DIFF in the QD paper (Hida et al. 2000).
       Subtract 2 f107_pf's, returning the 2 components of a normalized dd_real
    */
    // sub_112
    //  input: f53 a; f53 b
    // output: (s,e) = a - b
    //    ops: 6
    static f107_pf sub_112(f107_pf a, f107_pf b, f107_pf &e)
    {
      f107_pf s = a - b;
      f107_pf t = s - a;
      e = (a-(s-t))-(b+t);
      return s;
    }


    // Here are the real comparison operators; these compare both parts
    // of the f107 when necessary.

    // operator <
    //  input: f107 a; f53 b
    // output: int (a < b)
    //    ops: %%    len: 7
    friend int operator < ( const f107_o lhs, f107_pf rhs )
    {
      return ((lhs.c0 < rhs) || ((lhs.c0 == rhs) && (lhs.c1 < 0)));
    }

    // operator <
    //  input: f53 a; f107 b
    // output: int (a < b)
    //    ops: %%     len: 7
    friend int operator < ( f107_pf lhs, const f107_o rhs )
    {
      return ((lhs < rhs.c0) || ((lhs == rhs.c0) && (rhs.c1 > 0)));
    }

    // operator >
    //  input: f53 a; f107 b
    // output: int (a > b)
    //    ops: %%     len: 7
    friend int operator > ( f107_pf lhs, const f107_o rhs )
    {
      return ((lhs > rhs.c0) || ((lhs == rhs.c0) && (rhs.c1 < 0)));
    }

    // operator >
    //  input: f107 a; f53 b
    // output: int (a > b)
    //    ops: %%     len: 7
    friend int operator > ( const f107_o lhs, f107_pf rhs )
    {
      return ((lhs.c0 > rhs) || ((lhs.c0 == rhs) && (lhs.c1 > 0)));
    }

    /* normalize_32 is used by the strict versions of add and mul */
    // normalize_22
    //  input: f53 x0; f53 x1; f53 x2
    // output: (s0, s1) ~= x0 + x1 + x2
    //    ops: 6    length: 12
    static f107_pf normalize_32(f107_pf x0, f107_pf x1, f107_pf x2, f107_pf &s1)
    {
      f107_pf s, e, t2, s0;

      s = f107_o::normalize_22(x1, x2,   t2);   // 3
      s0 = f107_o::normalize_22(x0, s,   e);    // 3

      if (e != 0) {
        s1 = e + t2;
      } else {
        s0 = f107_o::normalize_22(s0, t2,   s1);
      }
      
      return s0;
    }

    // operator *
    //  input: f53 a; f107 b
    // output: f107 x = a * b
    //    ops: 8
    friend f107_o operator * ( f107_pf a, const f107_o b )
    {
      f107_o  res;
      f107_pf t0, t1, b1;
#if F107_APPROX

      t0 = f107_o::mul_112(a, b.c0,   b1);                // 3
      t1 = b1 + a*b.c1;                           // 2
      res.c0 = f107_o::normalize_22(t0, t1,   res.c1);    // 3
#else
      f107_pf d1, d2, u2, t2;

      t0 = f107_o::mul_112(a, b.c0,   b1);

      d1 = f107_o::mul_112(a, b.c1,   d2);

      t1 = f107_o::add_112(b1, d1,   u2);
      t2 = d2 + u2;

      res.c0 = f107_o::normalize_32(t0, t1, t2,   res.c1);
#endif

      return res;
    }

    // operator *
    //  input: f107 a; f53 b
    // output: f107 x = a * b
    //    ops: 8
    friend f107_o operator * ( const f107_o a, f107_pf b )
    {
      f107_o  res;
      f107_pf t0, t1, b1;
#if F107_APPROX

      t0 = f107_o::mul_112(a.c0, b,   b1);
      t1 = b1 + a.c1*b;
      res.c0 = f107_o::normalize_22(t0, t1,   res.c1);
#else
      f107_pf e1, e2, u2, t2;

      t0 = f107_o::mul_112(a.c0, b,   b1);

      e1 = f107_o::mul_112(a.c1, b,   e2);

      t1 = f107_o::add_112(b1, e1,   u2);
      t2 = e2 + u2;

      res.c0 = f107_o::normalize_32(t0, t1, t2, res.c1);
#endif

      return res;
    }

    // operator *
    //  input: f107 a; int b
    // output: f107 x = a * b
    //    ops: 8
    friend f107_o operator * ( const f107_o a, int i )
    {
      f107_pf b;
      f107_o  res;
      f107_pf t0, t1, b1;

      b = i;
#if F107_APPROX
      t0 = f107_o::mul_112(a.c0, b,   b1);
      t1 = b1 + a.c1*b;
      res.c0 = f107_o::normalize_22(t0, t1,   res.c1);
#else
      f107_pf e1, e2, u2, t2;

      t0 = f107_o::mul_112(a.c0, b,   b1);

      e1 = f107_o::mul_112(a.c1, b,   e2);

      t1 = f107_o::add_112(b1, e1,   u2);
      t2 = e2 + u2;

      res.c0 = f107_o::normalize_32(t0, t1, t2, res.c1);
#endif

      return res;
    }

    // square_22
    //  input: f107 a
    // output: f107 x = a * a
    //    ops: 9
    static f107_o f107_native_square(f107_o a)
    {
      f107_o res;
      f107_pf t0, t1, u1, v1;
#if F107_APPROX

      t0 = f107_o::square_12(a.c0,   u1);               // 3
      v1 = a.c0 * a.c1;                                 // 1
      t1 = (2.0 * v1) + u1;                             // 2
      res.c0 = f107_o::normalize_22(t0, t1,   res.c1);  // 3
#else
      f107_pf t2, v2, w2;

      t0 = square_12(a.c0,   u1);

      v1 = f107_o::mul_112(a.c0, a.c1,   v2);

      v1 = 2.0 * v1;
      v2 = 2.0 * v2;

      t1 = f107_o::add_112(u1, v1,   w2);

      t2 = v2 + w2;

      res.c0 = f107_o::normalize_32(t0, t1, t2,   res.c1);
#endif

      return res;
    }

    /* add_1112 is the THREE-SUM that returns only two components.
       Adds 3 f107_pf's, returning the 2 components of a normalized dd_real
    */
    // add_1112
    //  input: f53 a; f53 b; f52 c
    // output: (r0, r1) = a + b + c
    //    ops: 13
    static f107_pf add_1112(f107_pf x, f107_pf y, f107_pf z, f107_pf &r1)
    {
      f107_pf u, v, w, r0;

      u = f107_o::add_112(x, y,   v);    // 6
      r0 = f107_o::add_112(u, z,   w);   // 6
      r1 = v + w;                        // 1
      return r0;
    }

    // operator *
    //  input: f107 a; f107 b
    // output: f107 x = a * b
    //    ops: 10
    friend f107_o operator * ( const f107_o a, const f107_o b )
    {
      f107_o  res;
      f107_pf t0, t1, b1;
#if F107_APPROX

      t0 = f107_o::mul_112(a.c0, b.c0,   b1);           // 3
      t1 = a.c0*b.c1 + b1;                              // 2
      t1 = a.c1*b.c0 + t1;                              // 2
      res.c0 = f107_o::normalize_22(t0, t1,   res.c1);  // 3
#else
      f107_pf d1, d2, e1, e2, g2, u2, t2;

      t0 = f107_o::mul_112(a.c0, b.c0,   b1);

      d1 = f107_o::mul_112(a.c0, b.c1,   d2);
      e1 = f107_o::mul_112(a.c1, b.c0,   e2);

      g2 = a.c1 * b.c1;

      t1 = f107_o::add_1112(b1, d1, e1,   u2);
      t2 = d2 + e2 + g2 + u2;

      res.c0 = f107_o::normalize_32(t0, t1, t2,   res.c1);
#endif

      return res;
    }

    friend f107_o operator * ( const f107_o lhs, const char * rhs );

    friend f107_o operator * ( const char * lhs, const f107_o rhs );

    // operator +
    //  input: f107 a; f53 b
    // output: f107 x = a + b
    //    ops: 10
    friend f107_o operator + ( const f107_o lhs, f107_pf rhs )
    {
      f107_o  res;
      f107_pf t0, t1, u;
#if F107_APPROX

      t0 = f107_o::add_112(lhs.c0, rhs,   u);             // 6
      t1 = u + lhs.c1;                                    // 1
      res.c0 = f107_o::normalize_22(t0, t1,   res.c1);    // 3
#else
      f107_pf t2;

      t0 = f107_o::add_112(lhs.c0, rhs,   u);
      t1 = f107_o::add_112(lhs.c1, u,   t2);
      res.c0 = f107_o::normalize_32(t0, t1, t2,   res.c1);
#endif
  
      return res;
    }

    // operator +
    //  input: f53 a; f107 b
    // output: f107 x = a + b
    //    ops: 10
    friend f107_o operator + ( f107_pf lhs, const f107_o rhs )
    {
      f107_o  res;
      f107_pf t0, t1, u;
#if F107_APPROX

      t0 = f107_o::add_112(lhs, rhs.c0,   u);             // 6
      t1 = u + rhs.c1;                                    // 1
      res.c0 = f107_o::normalize_22(t0, t1,   res.c1);    // 3
#else
      f107_pf t2;

      t0 = f107_o::add_112(lhs, rhs.c0,   u);
      t1 = f107_o::add_112(u, rhs.c1,   t2);
      res.c0 = f107_o::normalize_32(t0, t1, t2,   res.c1);
#endif
  
      return res;
    }

    // operator -
    //  input: f107 a; f53 b
    // output: f107 x = a - b
    //    ops: 10
    friend f107_o operator - ( const f107_o lhs, f107_pf rhs )
    {
      f107_o  res;
      f107_pf t0, t1, u;
#if F107_APPROX

      t0 = f107_o::sub_112(lhs.c0, rhs,   u);

      t1 = u + lhs.c1;

      res.c0 = f107_o::normalize_22(t0, t1,   res.c1);
#else
      f107_pf t2;

      t0 = f107_o::sub_112(lhs.c0, rhs,   u);

      t1 = f107_o::add_112(lhs.c1, u,   t2);

      res.c0 = f107_o::normalize_32(t0, t1, t2,   res.c1);
#endif

      return res;
    }

    // operator -
    //  input: f53 a; f107 b
    // output: f107 x = a - b
    //    ops: 10
    friend f107_o operator - ( f107_pf lhs, const f107_o rhs )
    {
      f107_o  res;
      f107_pf t0, t1, u;
#if F107_APPROX

      t0 = f107_o::sub_112(lhs, rhs.c0,   u);

      t1 = u - rhs.c1;

      res.c0 = f107_o::normalize_22(t0, t1,   res.c1);
#else
      f107_pf t2;

      t0 = f107_o::sub_112(lhs, rhs.c0,   u);

      t1 = f107_o::sub_112(u, rhs.c1,   t2);

      res.c0 = f107_o::normalize_32(t0, t1, t2,   res.c1);
#endif

      return res;
    }

    // operator +
    //  input: f107 a; f107 b
    // output: f107 x = a * b
    //    ops: 11
    friend f107_o operator + ( const f107_o lhs, const f107_o rhs )
    {
      f107_o  res;
      f107_pf t0, t1, u;
#if F107_APPROX

      t0 = f107_o::add_112(lhs.c0, rhs.c0,   u);       // 6
      t1 = (u + lhs.c1)+(rhs.c1);                      // 2
      res.c0 = f107_o::normalize_22(t0, t1,   res.c1); // 3
#else
      f107_pf t2, v, w, x;

      t0 = f107_o::add_112(lhs.c0, rhs.c0,   u);
      v = f107_o::add_112(lhs.c1, rhs.c1,   w);
      t1 = f107_o::add_112(u, v,   x);
      t2 = w + x;
      res.c0 = f107_o::normalize_32(t0, t1, t2,   res.c1);
#endif
  
      return res;
    }

    friend f107_o operator + ( const f107_o lhs, const char * rhs );

    friend f107_o operator + ( const char * lhs, const f107_o rhs );

    // operator -
    //  input: f107 a; f107 b
    // output: f107 x = a - b
    //    ops: 11
    friend f107_o operator - ( const f107_o lhs, const f107_o rhs )
    {
      f107_o  res;
      f107_pf t0, t1, u;
#if F107_APPROX

      t0 = f107_o::sub_112(lhs.c0, rhs.c0,   u);

      t1 = (u+lhs.c1)-(rhs.c1);

      res.c0 = f107_o::normalize_22(t0, t1,   res.c1);
#else
      f107_pf t2, v, w, x;

      t0 = f107_o::sub_112(lhs.c0, rhs.c0,   u);
      v = f107_o::sub_112(lhs.c1, rhs.c1,   w);

      t1 = f107_o::add_112(u, v,   x);
      t2 = w + x;

      res.c0 = f107_o::normalize_32(t0, t1, t2,   res.c1);
#endif

      return res;
    }

    friend f107_o operator - ( const f107_o lhs, const char * rhs );

    friend f107_o operator - ( const char * lhs, const f107_o rhs );



    // cmp_f107
    //  input: f107 a; f107 b
    // output: int (a <=> b)
    //    ops: %%    len: 12
    // 20080311: changed "static" to "friend"
    // 20120112: changed it back: it should be static, which works so long
    //    as the caller invokes it as "f107_o::cmp_f107(1,b)"
    static int cmp_f107(f107_o a, f107_o b)
    {
      int  result;
      
      if (a.c0 > b.c0) {
        result = 1;
      } else if (a.c0 < b.c0) {
        result = -1;
      } else {
        if (a.c1 > b.c1) {
          result = 1;
        } else if (a.c1 < b.c1) {
          result = -1;
        } else {
          result = 0;
        }
      }
      
      return(result);
    }

    // operator <
    //  input: f107 a; f107 b
    // output: int (a < b)
    //    ops: %%     len: 14
    friend int operator < ( const f107_o lhs, const f107_o rhs )
    {
      return(cmp_f107(lhs, rhs) < 0);
    }

    // operator >
    //  input: f107 a; f107 b
    // output: int (a > b)
    //    ops: %%     len: 14
    friend int operator > ( const f107_o lhs, const f107_o rhs )
    {
      return (cmp_f107(lhs, rhs) > 0);
    }



    friend f107_o operator / ( const f107_o, const f107_o );
    friend f107_o operator / ( f107_pf, const f107_o );
    friend f107_o operator / ( const f107_o, f107_pf );

    f107_o & operator += ( f107_o );
    f107_o & operator += ( f107_pf );
    f107_o & operator -= ( f107_o );
    f107_o & operator -= ( f107_pf );
    f107_o & operator *= ( f107_o );
    f107_o & operator *= ( f107_pf );
    f107_o & operator /= ( f107_o );
    f107_o & operator /= ( f107_pf );

    friend int operator >= ( const f107_o, const f107_o );
    friend int operator >= ( const f107_o, f107_pf );
    friend int operator >= ( f107_pf, const f107_o );
    friend int operator <= ( const f107_o, const f107_o );
    friend int operator <= ( const f107_o, f107_pf );
    friend int operator <= ( f107_pf, const f107_o );
    friend int operator != ( const f107_o, const f107_o );
    friend int operator != ( f107_pf, const f107_o );
    friend int operator != ( const f107_o, f107_pf );
    friend int operator == ( const f107_o, const f107_o );
    friend int operator == ( f107_pf, const f107_o );
    friend int operator == ( const f107_o, f107_pf );

    f107_o square ( const f107_o z );
    f107_o trunc ( const f107_o z );
    f107_o fabs ( const f107_o z );
    f107_o floor ( const f107_o x);

    f107_o sqrt ( const f107_o x );
    f107_o f107_sincos(f107_o xq, f107_o zq, f107_o uq, double i);
    f107_o sin ( const f107_o x );
    f107_o cos ( const f107_o x );

    f107_o restrict_tanh(f107_o x);
    f107_o exp(f107_o x);
    f107_o sinh(f107_o x);
    f107_o tanh(f107_o x);
    f107_o log(f107_o x);

//    friend void f107_to_digits(f107_o x, char *s, int *expn, int *sign,
//                  int precision);

    f107_pf add_P4_2(f107_pf a, f107_pf b, f107_pf c, f107_pf d, f107_pf & r1);

    // 20080311: commented out the following 2
    // f107_pf mul_112(f107_pf a, f107_pf b, f107_pf &e);
    // f107_pf square_12(f107_pf a, f107_pf &e);
};

#endif

/* end of f107_o.h */