// Excerpts from:
// doubledouble.cc
// Excerpted and modified by Steve Vavasis, 2/00 for use in QMG 2.0.
// doubledouble package is originally written by Keith Briggs.  
// These excerpts are covered by the
// QMG license, with permission from Keith Briggs.  
// The entire doubledouble package is available on line at
//   http://www-epidem.plantsci.cam.ac.uk/~kbriggs/doubledouble.html
// and is covered by the Gnu Public License.


// KMB 98 Jan 19
// Keith Briggs.   Last revised 98 Feb 09


/*
Copyright (C) 1997 Keith Martin Briggs

Except where otherwise indicated,
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 for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/

#include "doubledouble2.h"


// member constants
const double doubledouble::Split = 134217729.0L; // 2^27+1, for IEEE double


// Square  (faster than x*x)
doubledouble sqr(const doubledouble& x) {
  double hx,tx,C,c;
  C=doubledouble::Split*x.h(); hx=C-x.h(); hx=C-hx; tx=x.h()-hx;
  C=x.h()*x.h();
  c=((((hx*hx-C)+2.0*hx*tx))+tx*tx)+2.0*x.h()*x.l();
  hx=C+c;
  doubledouble r=doubledouble(hx,c+(C-hx));
  return r;
}

// Added by Steve Vavasis
#ifdef _MSC_VER
#pragma optimize("g", off)
#endif

// square root
doubledouble sqrt(const doubledouble& y) {
  double c,p,q,hx,tx,u,uu,cc,hi=y.h();
  if (hi<0.0) {
    doubledouble::Domain_Error d;
    throw d;
  }
  if (0.0==hi) return y;
  c=sqrt(hi);
  p=doubledouble::Split*c; hx=c-p; hx+=p; tx=c-hx;
  p=hx*hx;
  q=2.0*hx*tx;
  u=p+q;
  uu=(p-u)+q+tx*tx;
  cc=(((y.h()-u)-uu)+y.l())/(c+c);
  u=c+cc;
  doubledouble r=doubledouble(u,cc+(c-u));
  return r;
}

#ifdef _MSC_VER
#pragma optimize("g", on)
#define x86
#endif



doubledouble::FloatControlType doubledouble::set_float_control_to_53_bits() {

#ifdef x86
  volatile FloatControlType old_cw;
#ifdef _MSC_VER
  __asm fnstcw old_cw 
    ;
#else
#ifdef __GNUC__
  asm volatile ("fnstcw %0":"=m" (old_cw)); 
#else
  ERROR_MESSAGE("COMPILER UNKNOWN -- NOT KNOWN HOW TO INSTANTIATE FNSTCW ASSEMBLER INSTRUCTION FOR X86");
#endif
#endif
  volatile FloatControlType new_cw = (old_cw & ~0x300) | 0x200;
  
#ifdef _MSC_VER
  __asm fldcw new_cw
    ;
#else
#ifdef __GNUC__
  asm volatile ("fldcw %0": :"m" (new_cw));
#else
  ERROR_MESSAGE("COMPILER UNKNOWN -- NOT KNOWN HOW TO INSTANTIATE FLDCW ASSEMBLER INSTRUCTION FOR X86");
#endif
#endif
  
  return old_cw;
#else  // not X86
  return 0;
#endif
}

void doubledouble::restore_float_control(FloatControlType w) {
#ifdef x86
#ifdef _MSC_VER
  __asm fldcw w
    ;
#else
#ifdef __GNUC__
  asm volatile ("fldcw %0": :"m" (w));
#else
  ERROR_MESSAGE("COMPILER UNKNOWN -- NOT KNOWN HOW TO INSTANTIATE FLDCW ASSEMBLER INSTRUCTION FOR X86");
#endif
#endif
#endif
}