Quelle cxsc.C Sprache: C

/****************************************************************************
**
*W  cxsc.C                       GAP source                 Laurent Bartholdi
**
*Y  Copyright (C) 2008-2012 Laurent Bartholdi
**
**  This file contains the implementation of the CXSC float package.
*/
#include "floatconfig.h"

#undef PACKAGE
#undef PACKAGE_BUGREPORT
#undef PACKAGE_NAME
#undef PACKAGE_STRING
#undef PACKAGE_TARNAME
#undef PACKAGE_URL
#undef PACKAGE_VERSION

#include "floattypes.h"

#define ERROR_CXSC(gap_name,obj)          \
  ErrorQuit(#gap_name ": argument must be a CXSC float, not a %s",    \
     (Int)TNAM_OBJ(obj),0)

#ifdef __cplusplus
static inline bool HAS_FILTER(Obj obj, Obj filter)
{
  return DoFilter(filter,obj) == True;
  return IS_DATOBJ(obj) && DoFilter(filter,obj) == True;
}
#endif
#define IS_RP(obj) HAS_FILTER(obj,IS_CXSC_RP)
#define TEST_IS_RP(gap_name,obj)    \
  if (!IS_RP(obj))      \
    ErrorQuit(#gap_name ": expected a real, not a %s",  \
        (Int)TNAM_OBJ(obj),0)

#define IS_CP(obj) HAS_FILTER(obj,IS_CXSC_CP)
#define TEST_IS_CP(gap_name,obj)    \
  if (!IS_CP(obj))      \
    ErrorQuit(#gap_name ": expected a complex, not a %s", \
        (Int)TNAM_OBJ(obj),0)

#define IS_RI(obj) HAS_FILTER(obj,IS_CXSC_RI)
#define TEST_IS_RI(gap_name,obj)           \
  if (!IS_RI(obj))      \
    ErrorQuit(#gap_name ": expected an interval, not a %s", \
        (Int)TNAM_OBJ(obj),0)

#define IS_CI(obj) HAS_FILTER(obj,IS_CXSC_CI)
#define TEST_IS_CI(gap_name,obj)           \
  if (!IS_CI(obj))             \
    ErrorQuit(#gap_name ": expected a complex interval, not a %s",\
        (Int)TNAM_OBJ(obj),0)

/****************************************************************
* cxsc data are stored as follows:
* +--------------------+----------+
* | TYPE_CXSC_RI       | interval  |
* +--------------------+----------+
****************************************************************/
#define RP_OBJ(obj) (*(cxsc::real *) (ADDR_OBJ(obj)+1))
#define RI_OBJ(obj) (*(cxsc::interval *) (ADDR_OBJ(obj)+1))
#define CP_OBJ(obj) (*(cxsc::complex *) (ADDR_OBJ(obj)+1))
#define CI_OBJ(obj) (*(cxsc::cinterval *) (ADDR_OBJ(obj)+1))

#include "except.hpp"
#include "real.hpp"
#include "complex.hpp"
#include "interval.hpp"
#include "cinterval.hpp"

#include "cxsc_poly.h"

#include "cpoly.hpp"
#include "cipoly.hpp"
#include "cpzero.hpp"
#include "rpoly.hpp"
#include "rpeval.hpp"

// this function is missing from CXSC 2.5.1
cxsc::complex cxsc::sqr(const cxsc::complex&z) throw() { return z*z; }

// this function is missing from CXSC 2.5.1
bool Disjoint(cxsc::cinterval &a, cxsc::cinterval &b) {
  return Disjoint(Re(a),Re(b)) || Disjoint(Im(a),Im(b)); }

Obj TYPE_CXSC_RP, TYPE_CXSC_CP, TYPE_CXSC_RI, TYPE_CXSC_CI;
Obj IS_CXSC_RP, IS_CXSC_CP, IS_CXSC_RI, IS_CXSC_CI;
Obj FAMILY_CXSC;
Obj GAPLog2Int; // sometimes it's exported as FuncLog2Int, sometimes not

/****************************************************************
* creators
****************************************************************/
static inline Obj NEW_RP (void)
{
  return NEW_DATOBJ(sizeof(cxsc::real), TYPE_CXSC_RP);
}
static inline Obj NEW_CP (void)
{
  return NEW_DATOBJ(sizeof(cxsc::complex), TYPE_CXSC_CP);
}
static inline Obj NEW_RI (void)
{
  return NEW_DATOBJ(sizeof(cxsc::interval), TYPE_CXSC_RI);
}
static inline Obj NEW_CI (void)
{
  return NEW_DATOBJ(sizeof(cxsc::cinterval), TYPE_CXSC_CI);
}

static inline Obj OBJ_RP (cxsc::real i)
{
  Obj f = NEW_RP();
  RP_OBJ(f) = i;
  return f;
}
static inline Obj OBJ_CP (cxsc::complex i)
{
  Obj f = NEW_CP();
  CP_OBJ(f) = i;
  return f;
}
static inline Obj OBJ_RI (cxsc::interval i)
{
  Obj f = NEW_RI();
  RI_OBJ(f) = i;
  return f;
}
static inline Obj OBJ_CI (cxsc::cinterval i)
{
  Obj f = NEW_CI();
  CI_OBJ(f) = i;
  return f;
}

static inline cxsc::real RP_GET(cxsc::real x) { return x; }
static inline cxsc::complex CP_GET(cxsc::real x) { return _complex(x); }
static inline cxsc::interval RI_GET(cxsc::real x) { return _interval(x); }
static inline cxsc::cinterval CI_GET(cxsc::real x) { return _cinterval(x); }

/****************************************************************
* 1-argument functions
****************************************************************/
#define Func1_CXSC(gap_name,cxsc_name)    \
  static Obj gap_name##_RP(Obj self, Obj f)   \
  {        \
    TEST_IS_RP(gap_name##_RP,f);    \
    if (IsQuietNaN(RP_OBJ(f))) { return f; }   \
    return OBJ_RP(cxsc_name(RP_OBJ(f)));   \
  }        \
  static Obj gap_name##_CP(Obj self, Obj f)   \
  {        \
    TEST_IS_CP(gap_name##_CP,f);    \
    if (IsQuietNaN(Re(CP_OBJ(f)))) { return f; }  \
    return OBJ_CP(cxsc_name(CP_OBJ(f)));   \
  }        \
  static Obj gap_name##_RI(Obj self, Obj f)   \
  {        \
    TEST_IS_RI(gap_name##_RI,f);    \
    if (IsQuietNaN(Inf(RI_OBJ(f)))) { return f; }  \
    return OBJ_RI(cxsc_name(RI_OBJ(f)));   \
  }        \
  static Obj gap_name##_CI(Obj self, Obj f)   \
  {        \
    TEST_IS_CI(gap_name##_CI,f);    \
    if (IsQuietNaN(Inf(Re(CI_OBJ(f))))) { return f; }  \
    return OBJ_CI(cxsc_name(CI_OBJ(f)));   \
  }

#define Func1a_CXSC(gap_name,result,cxsc_name)   \
  static Obj gap_name##_CP(Obj self, Obj f)   \
  {        \
    TEST_IS_CP(gap_name##_CP,f);    \
    if (IsQuietNaN(Re(CP_OBJ(f)))) { return f; }  \
    return OBJ_##result##P(cxsc_name(CP_OBJ(f)));  \
  }        \
  static Obj gap_name##_CI(Obj self, Obj f)   \
  {        \
    TEST_IS_CI(gap_name##_CI,f);    \
    if (IsQuietNaN(Inf(Re(CI_OBJ(f))))) { return f; }  \
    return OBJ_##result##I(cxsc_name(CI_OBJ(f)));  \
  }

#define Func1b_CXSC(gap_name,cxsc_name)    \
  static Obj gap_name##_RI(Obj self, Obj f)   \
  {        \
    TEST_IS_RI(gap_name##_RI,f);    \
    if (IsQuietNaN(Inf(RI_OBJ(f)))) { return f; }  \
    return OBJ_RP(cxsc_name(RI_OBJ(f)));   \
  }        \
  static Obj gap_name##_CI(Obj self, Obj f)   \
  {        \
    TEST_IS_CI(gap_name##_CI,f);    \
    if (IsQuietNaN(Inf(Re(CI_OBJ(f))))) { return f; }  \
    return OBJ_CP(cxsc_name(CI_OBJ(f)));   \
  }

#define Inc1_CXSC_arg(name,arg)     \
  { #name, 1, arg, (ObjFunc) name, "src/cxsc_float.c:" #name }
#define Inc1_CXSC(name)     \
  Inc1_CXSC_arg(name##_RP,"cxsc::rp"),   \
    Inc1_CXSC_arg(name##_CP,"cxsc::cp"),  \
    Inc1_CXSC_arg(name##_RI,"cxsc::ri"),  \
    Inc1_CXSC_arg(name##_CI,"cxsc::ci")
#define Inc1a_CXSC(name)    \
  Inc1_CXSC_arg(name##_CP,"cxsc::cp"),   \
    Inc1_CXSC_arg(name##_CI,"cxsc::ci")
#define Inc1b_CXSC(name)    \
    Inc1_CXSC_arg(name##_RI,"cxsc::ri"),  \
    Inc1_CXSC_arg(name##_CI,"cxsc::ci")

static Obj CXSC_INT (Obj self, Obj f)
{
  TEST_IS_INTOBJ(CXSC_INT,f);
  return OBJ_RP(Double(INT_INTOBJ(f)));
}

#define VAL_MACFLOAT(obj) (*(Double *)ADDR_OBJ(obj))
#define IS_MACFLOAT(obj) (TNUM_OBJ(obj) == T_MACFLOAT)

static Obj CXSC_IEEE754 (Obj self, Obj f)
{
  if (!IS_MACFLOAT(f)) {
    ErrorMayQuit("CXSC_IEEE754: object must be a float, not a %s",
                       (Int)TNAM_OBJ(f),0);
  }
  return OBJ_RP(VAL_MACFLOAT(f));
}

static Obj CP_CXSC_RP_RP (Obj self, Obj f, Obj g)
{
  TEST_IS_RP(CP_CXSC_RP_RP,f);
  TEST_IS_RP(CP_CXSC_RP_RP,g);
  return OBJ_CP(cxsc::complex(RP_OBJ(f),RP_OBJ(g)));
}
static Obj CP_CXSC_RP (Obj self, Obj f)
{
  TEST_IS_RP(CP_CXSC_RP,f);
  return OBJ_CP(cxsc::complex(RP_OBJ(f)));
}
static Obj RI_CXSC_RP (Obj self, Obj f)
{
  TEST_IS_RP(RI_CXSC_RP,f);
  return OBJ_RI(cxsc::interval(RP_OBJ(f)));
}
static Obj RI_CXSC_RP_RP (Obj self, Obj f, Obj g)
{
  TEST_IS_RP(RI_CXSC_RP_RP,f);
  TEST_IS_RP(RI_CXSC_RP_RP,g);
  return OBJ_RI(cxsc::interval(RP_OBJ(f),RP_OBJ(g)));
}
static Obj CI_CXSC_CP (Obj self, Obj f)
{
  TEST_IS_CP(CI_CXSC_CP,f);
  return OBJ_CI(cxsc::cinterval(CP_OBJ(f)));
}
static Obj CI_CXSC_RI_RI (Obj self, Obj f, Obj g)
{
  TEST_IS_RI(CI_CXSC_RI_RI,f);
  TEST_IS_RI(CI_CXSC_RI_RI,g);
  return OBJ_CI(cxsc::cinterval(RI_OBJ(f),RI_OBJ(g)));
}
static Obj CI_CXSC_CP_CP (Obj self, Obj f, Obj g)
{
  TEST_IS_CP(CI_CXSC_CP_CP,f);
  TEST_IS_CP(CI_CXSC_CP_CP,g);
  return OBJ_CI(cxsc::cinterval(CP_OBJ(f),CP_OBJ(g)));
}

static Obj CXSC_NEWCONSTANT (Obj self, Obj f)
{
  TEST_IS_INTOBJ(CXSC_NEWCONSTANT,f);
  switch (INT_INTOBJ(f)) {
  case 0: return OBJ_RP(cxsc::MinReal);
  case 1: return OBJ_RP(cxsc::minreal);
  case 2: return OBJ_RP(cxsc::MaxReal);
  case 3: return OBJ_RP(cxsc::Infinity);
  case 4: return OBJ_RP(cxsc::SignalingNaN);
  case 5: return OBJ_RP(cxsc::QuietNaN);
  case 6: return OBJ_RP(cxsc::Pi_real);        // Pi
  case 7: return OBJ_RP(cxsc::Pi2_real);       // 2*Pi
  case 8: return OBJ_RP(cxsc::Pi3_real);       // 3*Pi
  case 9: return OBJ_RP(cxsc::Pid2_real);      // Pi/2
  case 10: return OBJ_RP(cxsc::Pid3_real);      // Pi/3
  case 11: return OBJ_RP(cxsc::Pid4_real);      // Pi/4
  case 12: return OBJ_RP(cxsc::Pir_real);       // 1/Pi
  case 13: return OBJ_RP(cxsc::Pi2r_real);      // 1/(2*Pi)
  case 14: return OBJ_RP(cxsc::Pip2_real);      // Pi^2
  case 15: return OBJ_RP(cxsc::SqrtPi_real);    // sqrt(Pi)
  case 16: return OBJ_RP(cxsc::Sqrt2Pi_real);   // sqrt(2Pi)
  case 17: return OBJ_RP(cxsc::SqrtPir_real);   // 1/sqrt(Pi)
  case 18: return OBJ_RP(cxsc::Sqrt2Pir_real);  // 1/sqrt(2Pi)
  case 19: return OBJ_RP(cxsc::Sqrt2_real);     // sqrt(2)
  case 20: return OBJ_RP(cxsc::Sqrt2r_real);    // 1/sqrt(2)
  case 21: return OBJ_RP(cxsc::Sqrt3_real);     // sqrt(3)
  case 22: return OBJ_RP(cxsc::Sqrt3d2_real);   // sqrt(3)/2
  case 23: return OBJ_RP(cxsc::Sqrt3r_real);    // 1/sqrt(3)
  case 24: return OBJ_RP(cxsc::Ln2_real);       // ln(2)
  case 25: return OBJ_RP(cxsc::Ln2r_real);      // 1/ln(2)
  case 26: return OBJ_RP(cxsc::Ln10_real);      // ln(10)
  case 27: return OBJ_RP(cxsc::Ln10r_real);     // 1/ln(10)
  case 28: return OBJ_RP(cxsc::LnPi_real);      // ln(Pi)
  case 29: return OBJ_RP(cxsc::Ln2Pi_real);     // ln(2Pi)
  case 30: return OBJ_RP(cxsc::E_real);         // e
  case 31: return OBJ_RP(cxsc::Er_real);        // 1/e
  case 32: return OBJ_RP(cxsc::Ep2_real);       // e^2
  case 33: return OBJ_RP(cxsc::Ep2r_real);      // 1/e^2
  case 34: return OBJ_RP(cxsc::EpPi_real);      // e^(Pi)
  case 35: return OBJ_RP(cxsc::Ep2Pi_real);     // e^(2Pi)
  case 36: return OBJ_RP(cxsc::EpPid2_real);    // e^(Pi/2)
  case 37: return OBJ_RP(cxsc::EpPid4_real);    // e^(Pi/4)

  case 100: return OBJ_RI(cxsc::Pi_interval);        // Pi
  case 101: return OBJ_RI(cxsc::Pi2_interval);       // 2*Pi
  case 102: return OBJ_RI(cxsc::Pi3_interval);       // 3*Pi
  case 103: return OBJ_RI(cxsc::Pid2_interval);      // Pi/2
  case 104: return OBJ_RI(cxsc::Pid3_interval);      // Pi/3
  case 105: return OBJ_RI(cxsc::Pid4_interval);      // Pi/4
  case 106: return OBJ_RI(cxsc::Pir_interval);       // 1/Pi
  case 107: return OBJ_RI(cxsc::Pi2r_interval);      // 1/(2*Pi)
  case 108: return OBJ_RI(cxsc::Pip2_interval);      // Pi^2
  case 109: return OBJ_RI(cxsc::SqrtPi_interval);    // sqrt(Pi)
  case 110: return OBJ_RI(cxsc::Sqrt2Pi_interval);   // sqrt(2Pi)
  case 111: return OBJ_RI(cxsc::SqrtPir_interval);   // 1/sqrt(Pi)
  case 112: return OBJ_RI(cxsc::Sqrt2Pir_interval);  // 1/sqrt(2Pi)
  case 113: return OBJ_RI(cxsc::Sqrt2_interval);     // sqrt(2)
  case 114: return OBJ_RI(cxsc::Sqrt2r_interval);    // 1/sqrt(2)
  case 115: return OBJ_RI(cxsc::Sqrt3_interval);     // sqrt(3)
  case 116: return OBJ_RI(cxsc::Sqrt3d2_interval);   // sqrt(3)/2
  case 117: return OBJ_RI(cxsc::Sqrt3r_interval);    // 1/sqrt(3)
  case 118: return OBJ_RI(cxsc::Ln2_interval);       // ln(2)
  case 119: return OBJ_RI(cxsc::Ln2r_interval);      // 1/ln(2)
  case 120: return OBJ_RI(cxsc::Ln10_interval);      // ln(10)
  case 121: return OBJ_RI(cxsc::Ln10r_interval);     // 1/ln(10)
  case 122: return OBJ_RI(cxsc::LnPi_interval);      // ln(Pi)
  case 123: return OBJ_RI(cxsc::Ln2Pi_interval);     // ln(2Pi)
  case 124: return OBJ_RI(cxsc::E_interval);         // e
  case 125: return OBJ_RI(cxsc::Er_interval);        // 1/e
  case 126: return OBJ_RI(cxsc::Ep2_interval);       // e^2
  case 127: return OBJ_RI(cxsc::Ep2r_interval);      // 1/e^2
  case 128: return OBJ_RI(cxsc::EpPi_interval);      // e^(Pi)
  case 129: return OBJ_RI(cxsc::Ep2Pi_interval);     // e^(2Pi)
  case 130: return OBJ_RI(cxsc::EpPid2_interval);    // e^(Pi/2)
  case 131: return OBJ_RI(cxsc::EpPid4_interval);    // e^(Pi/4)
  }
  return Fail;
}

static Obj INT_CXSC (Obj self, Obj f)
{
  TEST_IS_RP(INT_CXSC,f);
  int n = 0;
  int sign = 1;
  cxsc::real r = RP_OBJ(f);
  if (r < 0.0)
    r = -r, sign = -1;
  for (Int i = 1L << NR_SMALL_INT_BITS; i; i >>= 1)
    if (r >= Double(i))
      r = r-Double(i), n = n+i;
  if (r >= 1.0)
    return Fail;
  else
    return INTOBJ_INT(sign*n);
}

static Obj SIGN_CXSC_RP (Obj self, Obj f)
{
  TEST_IS_RP(SIGN_CXSC_RP,f);
  return INTOBJ_INT(sign(RP_OBJ(f)));
}

static Obj SIGN_CXSC_RI (Obj self, Obj f)
{
  TEST_IS_RI(SIGN_CXSC_RI,f);
  if (Inf(RI_OBJ(f))>0.0)
    return INTOBJ_INT(1);
  if (Sup(RI_OBJ(f))<0.0)
    return INTOBJ_INT(-1);
  if (RI_OBJ(f)==0.0)
    return INTOBJ_INT(0);
  return Fail;
}

#define Func1c_CXSC(gap_name,cxsc_name)    \
  static Obj gap_name##_RP(Obj self, Obj f)   \
  {        \
    TEST_IS_RP(gap_name##_RP,f);    \
    return cxsc_name(RP_OBJ(f)) ? True : False;   \
  }        \
  static Obj gap_name##_CP(Obj self, Obj f)   \
  {        \
    TEST_IS_CP(gap_name##_CP,f);    \
    return cxsc_name(CP_OBJ(f)) ? True : False;   \
  }        \
  static Obj gap_name##_RI(Obj self, Obj f)   \
  {        \
    TEST_IS_RI(gap_name##_RI,f);    \
    return cxsc_name(RI_OBJ(f)) ? True : False;   \
  }        \
  static Obj gap_name##_CI(Obj self, Obj f)   \
  {        \
    TEST_IS_CI(gap_name##_CI,f);    \
    return cxsc_name(CI_OBJ(f)) ? True : False;   \
  }

bool IsQuietNaN(cxsc::complex &x) {
  return IsQuietNaN(Re(x)) || IsQuietNaN(Im(x)); }
bool IsQuietNaN(cxsc::interval &x) {
  return IsQuietNaN(Inf(x)) || IsQuietNaN(Sup(x)); }
bool IsQuietNaN(cxsc::cinterval &x) {
  return IsQuietNaN(Re(x)) || IsQuietNaN(Im(x)); }
Func1c_CXSC(ISNAN_CXSC,IsQuietNaN)
bool IsInfinity(cxsc::complex &x) {
  return IsInfinity(Re(x)) || IsInfinity(Im(x)); }
bool IsInfinity(cxsc::interval &x) {
  return IsInfinity(Inf(x)) || IsInfinity(Sup(x)); }
bool IsInfinity(cxsc::cinterval &x) {
  return IsInfinity(Re(x)) || IsInfinity(Im(x)); }
Func1c_CXSC(ISXINF_CXSC,IsInfinity)
bool IsPInfinity(cxsc::real &x) { return IsInfinity(x) && x > 0.0; }
bool IsPInfinity(cxsc::interval &x) { return IsInfinity(x) && x > 0.0; }
bool IsPInfinity(cxsc::complex &x) { return IsInfinity(x); }
bool IsPInfinity(cxsc::cinterval &x) { return IsInfinity(x); }
Func1c_CXSC(ISPINF_CXSC,IsPInfinity)
bool IsNInfinity(cxsc::real &x) { return IsInfinity(x) && x < 0.0; }
bool IsNInfinity(cxsc::interval &x) { return IsInfinity(x) && x < 0.0; }
bool IsNInfinity(cxsc::complex &x) { return IsInfinity(x); }
bool IsNInfinity(cxsc::cinterval &x) { return IsInfinity(x); }
Func1c_CXSC(ISNINF_CXSC,IsNInfinity)
bool IsZero(cxsc::real &x) { return x == 0.0; }
bool IsZero(cxsc::interval &x) { return x == 0.0; }
bool IsZero(cxsc::complex &x) { return x == 0.0; }
bool IsZero(cxsc::cinterval &x) { return x == 0.0; }
Func1c_CXSC(ISZERO_CXSC,IsZero)
bool IsOne(cxsc::real &x) { return x == 1.0; }
bool IsOne(cxsc::interval &x) { return x == 1.0; }
bool IsOne(cxsc::complex &x) { return x == 1.0; }
bool IsOne(cxsc::cinterval &x) { return x == 1.0; }
Func1c_CXSC(ISONE_CXSC,IsOne)
bool IsNumber(cxsc::real &x) { return !IsInfinity(x) && !IsQuietNaN(x); }
bool IsNumber(cxsc::interval &x) { return !IsInfinity(x) && !IsQuietNaN(x); }
bool IsNumber(cxsc::complex &x) { return !IsInfinity(x) && !IsQuietNaN(x); }
bool IsNumber(cxsc::cinterval &x) { return !IsInfinity(x) && !IsQuietNaN(x); }
Func1c_CXSC(ISNUMBER_CXSC,IsNumber)

static Obj STRING_CXSC (Obj self, Obj f, Obj precision, Obj digits)
{
  std::string s;
  TEST_IS_INTOBJ(STRING_CXSC,precision);
  TEST_IS_INTOBJ(STRING_CXSC,digits);
  s << SetPrecision(INT_INTOBJ(precision),INT_INTOBJ(digits)) << Variable;

  if (IS_RP(f))
    s << RP_OBJ(f);
  else if (IS_CP(f))
    s << CP_OBJ(f);
  else if (IS_RI(f))
    s << RI_OBJ(f);
  else if (IS_CI(f))
    s << CI_OBJ(f);
  else ERROR_CXSC(STRING_CXSC,f);

  Obj str = NEW_STRING(s.length());
  s.copy (CSTR_STRING(str), s.length());

  return str;
}

cxsc::interval abs2(const cinterval& x)
{
  interval a=cxsc::abs(Re(x)), b=cxsc::abs(Im(x)), r;
  int exa=cxsc::expo(Sup(a)), exb=cxsc::expo(Sup(b)), ex;
    if (exb > exa)
    {  // Permutation of a,b:
        r = a;  a = b;  b = r;
        ex = exa;  exa = exb;  exb = ex;
    }
    ex = 511 - exa;
    cxsc::times2pown(a,ex);
    cxsc::times2pown(b,ex);
    r = a*a + b*b;
    times2pown(r,-ex);
    return r;
}

//cxsc::real abs(const complex& x) // why the %^(&^*( does the library not contain it?
//{
//  return sqrtx2y2(Re(x),Im(x));
//}

Func1a_CXSC(REAL_CXSC,R,Re);
Func1a_CXSC(IMAG_CXSC,R,Im);
Func1a_CXSC(ABS_CXSC,R,abs);
Func1a_CXSC(NORM_CXSC,R,abs2);
Func1a_CXSC(CONJ_CXSC,C,conj);

static Obj ISEMPTY_CXSC_RI (Obj self, Obj f)
{
  TEST_IS_RI(ISEMPTY_CXSC_RI,f);
  return IsEmpty(RI_OBJ(f)) ? True : False;
}

static Obj ISEMPTY_CXSC_CI (Obj self, Obj f)
{
  TEST_IS_CI(ISEMPTY_CXSC_RI,f);
  return IsEmpty(Re(CI_OBJ(f))) || IsEmpty(Im(CI_OBJ(f))) ? True : False;
}

cxsc::complex RelDiam(cxsc::cinterval z)
{
  if (z == 0.0)
    return cxsc::complex(0.0);
  return diam(z) / Sup(abs(z));
}

Func1b_CXSC(INF_CXSC,Inf);
Func1b_CXSC(SUP_CXSC,Sup);
Func1b_CXSC(MID_CXSC,mid);
Func1b_CXSC(DIAM_CXSC,diam);
Func1b_CXSC(DIAM_REL_CXSC,RelDiam);

static Obj RP_CXSC_STRING (Obj self, Obj str)
{
  TEST_IS_STRING(RP_CXSC_STRING,str);

  std::string s = CSTR_STRING(str);
  Obj f = NEW_RP();
  s >> RP_OBJ(f);
  return f;
}

static Obj CP_CXSC_STRING (Obj self, Obj str)
{
  TEST_IS_STRING(CP_CXSC_STRING,str);

  Obj f = NEW_CP();
  char *s = CSTR_STRING(str);
  if (s[0] == '(') {
    std::string(s) >> CP_OBJ(f);
    return f;
  }

  int sign = 1;
  bool inreal = true, empty = true;
  double newf;

  for (char *p = s;;) {
    switch (*p) {
    case '-':
    case '+':
    case 0:
      if (!empty) { /* drop the last read float */
if (inreal)
   Re(CP_OBJ(f)) += newf;
else
   Im(CP_OBJ(f)) += newf;
empty = inreal = true;
sign = 1;
      }
      if (!*p)
return f;
      if (*p == '-')
sign = -sign;
    case '*': p++; break;
    case 'i':
    case 'I': if (inreal) {
inreal = false;
if (empty) /* accept 'i' as '1*i' */
   Im(CP_OBJ(f)) = sign;
      } else return Fail;
      p++; break;
    default:
      {
int bytesread;
sscanf (p, "%lf%n", &newf, &bytesread);
if (bytesread == 0 && inreal)
   return Fail; /* no valid characters read */
if (sign == -1)
   newf = -newf;
empty = false;
p += bytesread;
      }
    }
  }
}

static Obj RI_CXSC_STRING (Obj self, Obj str)
{
  TEST_IS_STRING(RI_CXSC_STRING,str);

  std::string s = CSTR_STRING(str);
  Obj f = NEW_RI();
  if (s[0] == '[')
    s >> RI_OBJ(f);
  else {
    real l, r;
    std::string t = CSTR_STRING(str);
    s >> RndDown >> l;
    t >> RndUp >> r;
    RI_OBJ(f) = interval(l,r);
  }
  return f;
}

static Obj CI_CXSC_STRING (Obj self, Obj str)
{
  TEST_IS_STRING(CI_CXSC_STRING,str);

  std::string s = CSTR_STRING(str);
  Obj f = NEW_CI();
  if (s[0] == '[')
    s >> CI_OBJ(f);
  else if (s[0] == '(') {
    cxsc::complex l, r;
    std::string t = CSTR_STRING(str);
    s >> RndDown >> l;
    t >> RndUp >> r;
    CI_OBJ(f) = cinterval(l,r);
  } else {
    real l, r;
    std::string t = CSTR_STRING(str);
    char last = s[s.length()-1];
    s >> RndDown >> l;
    t >> RndUp >> r;
    if (last == 'i' || last == 'I')
      CI_OBJ(f) = cinterval(cxsc::complex(0.0,l),cxsc::complex(0.0,r));
    else
      CI_OBJ(f) = cinterval(cxsc::complex(l),cxsc::complex(r));
  }

  return f;
}

Func1_CXSC(INV_CXSC,1.0 /);
Func1_CXSC(AINV_CXSC,-);
Func1_CXSC(COS_CXSC,cxsc::cos);
Func1_CXSC(SIN_CXSC,cxsc::sin);
Func1_CXSC(TAN_CXSC,cxsc::tan);
Func1_CXSC(COT_CXSC,cxsc::cot);
Func1_CXSC(COSH_CXSC,cxsc::cosh);
Func1_CXSC(SINH_CXSC,cxsc::sinh);
Func1_CXSC(TANH_CXSC,cxsc::tanh);
Func1_CXSC(COTH_CXSC,cxsc::coth);
Func1_CXSC(ACOS_CXSC,cxsc::acos);
Func1_CXSC(ASIN_CXSC,cxsc::asin);
Func1_CXSC(ATAN_CXSC,cxsc::atan);
Func1_CXSC(ACOT_CXSC,cxsc::acot);
Func1_CXSC(ACOSH_CXSC,cxsc::acosh);
Func1_CXSC(ASINH_CXSC,cxsc::asinh);
Func1_CXSC(ATANH_CXSC,cxsc::atanh);
Func1_CXSC(ACOTH_CXSC,cxsc::acoth);
Func1_CXSC(SQR_CXSC,cxsc::sqr);
Func1_CXSC(SQRT_CXSC,cxsc::sqrt);
Func1_CXSC(EXP_CXSC,cxsc::exp);
Func1_CXSC(EXPM1_CXSC,cxsc::expm1);
Func1_CXSC(LOG_CXSC,cxsc::ln);
Func1_CXSC(LOG1P_CXSC,cxsc::lnp1);
Func1_CXSC(LOG2_CXSC,cxsc::log2);
Func1_CXSC(LOG10_CXSC,cxsc::log10);

static Obj ABS_CXSC_RP (Obj self, Obj f)
{
  TEST_IS_RP(ABS_CXSC_RP,f);
  return OBJ_RP(cxsc::abs(RP_OBJ(f)));
}
static Obj ABS_CXSC_RI (Obj self, Obj f)
{
  TEST_IS_RI(ABS_CXSC_RI,f);
  return OBJ_RI(cxsc::abs(RI_OBJ(f)));
}

#define MAXDEGREE 1000
#pragma GCC diagnostic ignored "-Wuninitialized"
static Obj ROOTPOLY_CXSC(Obj self, Obj gapcoeffs, Obj gapintervals)
{
  Int degree = LEN_PLIST(gapcoeffs)-1, numroots;
  bool intervals = false, real = true, complex = false;

  CPolynomial poly(degree);
  cxsc::complex coeffs[MAXDEGREE+1], roots[MAXDEGREE];

  if (degree > MAXDEGREE)
    return Fail;

  for (int i = 0; i <= degree; i++) {
    Obj c = ELM_PLIST(gapcoeffs,i+1);
    cxsc::complex z;
    if (IS_RP(c))
      z = RP_OBJ(c);
    else if (IS_CP(c))
      z = CP_OBJ(c), complex = true, real &= (Im(z) == 0.0);
    else if (IS_RI(c))
      z = cxsc::mid(RI_OBJ(c)), intervals = true;
    else if (IS_CI(c))
      z = cxsc::mid(CI_OBJ(c)), complex = intervals = true,
real &= (Im(CI_OBJ(c)) == 0.0);
    else ERROR_CXSC(ROOTPOLY_CXSC,c);
    poly[i] = coeffs[degree-i] = z;
  }

  numroots = cpoly_CXSC (degree, coeffs, roots, 53);

  if (numroots == -1)
    return Fail;

  Obj result = NEW_PLIST(T_PLIST, degree);
  SET_LEN_PLIST(result, degree);
  if (intervals) {
    cxsc::cinterval iroots[MAXDEGREE];

    for (int i = 0; i < numroots; i++) {
      CIPolynomial rp(degree);
      int error;
      CPolyZero(poly,roots[i],rp,iroots[i],error);
      if (error) {
iroots[i] = roots[i];
//fix! call InfoDecision(InfoFloat,1) and then maybe InfoDoPrint(...)
printf("#W CPOLYZERO failed to find enclosure for root %d; returning approximate root\n",i+1);
      }
    }

    // now, if the polynomial was real, force all isolated roots to be real
    if (real)
      for (int i = 0; i < numroots; i++)
if (in(0.0,Im(iroots[i]))) {
   bool lone = true;
   for (int j = 0; j < numroots; j++)
     if (j != i && !Disjoint(iroots[i],iroots[j])) {
       lone = false; break;
     }
   if (lone)
     iroots[i] = cinterval(Re(iroots[i]));
}

    for (int i = 1; i <= numroots; i++) {
      Obj gapz;
      if (!complex && Im(iroots[i-i]) == 0.0)
gapz = OBJ_RI(Re(iroots[i-1]));
      else
gapz = OBJ_CI(iroots[i-1]);
      SET_ELM_PLIST(result,i,gapz);
    }
  } else {
    // if the polynomial is real, and there doesn't exist a root at distance
    // <= 3*imaginarypart, force the imaginary part to be 0.
    if (real)
      for (int i = 0; i < numroots; i++) {
cxsc::real r = 10.0*sqr(Im(roots[i]));
bool lone = true;
for (int j = 0; j < numroots; j++)
   if (j != i && abs2(roots[i]-roots[j]) <= r) {
     lone = false; break;
   }
if (lone)
   roots[i] = cxsc::complex(Re(roots[i]));
      }

    for (int i = 1; i <= numroots; i++) {
      Obj gapz;
      if (!complex && Im(roots[i-1])==0.0)
gapz = OBJ_RP(Re(roots[i-1]));
      else
gapz = OBJ_CP(roots[i-1]);
      SET_ELM_PLIST(result,i,gapz);
    }
  }
  // some roots we missed
  for (int i = numroots+1; i <= degree; i++)
    SET_ELM_PLIST(result,i,Fail);

  return result;
}

static Obj EVALPOLY_CXSC(Obj self, Obj gapcoeffs, Obj gapt)
{
  TEST_IS_RP(EVALPOLY_CXSC,gapt);

  int degree = LEN_PLIST(gapcoeffs)-1, k, err;
  RPolynomial polyr(degree), polyi(degree);
  bool complex = false;

  for (int i = 0; i <= degree; i++) {
    Obj c = ELM_PLIST(gapcoeffs,i+1);
    if (IS_RP(c))
      polyr[i] = RP_OBJ(c);
    else if (IS_CP(c))
      polyr[i] = Re(CP_OBJ(c)), polyi[i] = Im(CP_OBJ(c)), complex = true;
    else
      ERROR_CXSC(EVALPOLY_CXSC,c);
  }

  interval intz[2];
  real z[2];

  RPolyEval (polyr, RP_OBJ(gapt), z[0], intz[0], k, err);
  if (err)
    return Fail;

  if (complex) {
    RPolyEval (polyi, RP_OBJ(gapt), z[1], intz[1], k, err);
    if (err)
      return Fail;
  }

  Obj list = NEW_PLIST(T_PLIST,2);
  SET_LEN_PLIST(list,2);
  Obj gapz, gapintz;
  if (complex)
    gapz = OBJ_CP(cxsc::complex(z[0],z[1])), gapintz = OBJ_CI(cxsc::cinterval(intz[0],intz[1]));
  else
    gapz = OBJ_RP(z[0]), gapintz = OBJ_RI(intz[0]);

  SET_ELM_PLIST(list,1,gapz);
  SET_ELM_PLIST(list,2,gapintz);

  return list;
}

/****************************************************************
* 2-argument functions
****************************************************************/
#define Inc2_CXSC_arg(name,arg)     \
  { #name, 2, arg, (ObjFunc) name, "src/cxsc_float.c:" #name }

#define Inc2_CXSC_X(name,argf)   \
  Inc2_CXSC_arg(name##_RP,argf ",csxc::rp"), \
    Inc2_CXSC_arg(name##_CP,argf ",cxsc::cp"), \
    Inc2_CXSC_arg(name##_RI,argf ",cxsc::ri"), \
    Inc2_CXSC_arg(name##_CI,argf ",cxsc::ci")

#define Inc2_CXSC(name)    \
  Inc2_CXSC_X(name##_RP,"cxsc::rp"),  \
    Inc2_CXSC_X(name##_CP,"cxsc::cp"),  \
    Inc2_CXSC_X(name##_RI,"cxsc::ri"),  \
    Inc2_CXSC_X(name##_CI,"cxsc::ci")

#define Func2_CXSC_X_X(gap_name,c,i,getf,getg,oper) \
  static Obj gap_name(Obj self, Obj f, Obj g)  \
  {       \
    return OBJ_##c##i(oper(getf(f),getg(g)));  \
  }

#define Func2_CXSC_X(gap_name,c,i,get,oper)  \
  Func2_CXSC_X_X(gap_name##_RP,c,i,get,RP_OBJ,oper); \
  Func2_CXSC_X_X(gap_name##_CP,C,i,get,CP_OBJ,oper); \
  Func2_CXSC_X_X(gap_name##_RI,c,I,get,RI_OBJ,oper); \
  Func2_CXSC_X_X(gap_name##_CI,C,I,get,CI_OBJ,oper); \

#define Func2_CXSC(gap_name,oper)   \
  Func2_CXSC_X(gap_name##_RP,R,P,RP_OBJ,oper);  \
  Func2_CXSC_X(gap_name##_CP,C,P,CP_OBJ,oper);  \
  Func2_CXSC_X(gap_name##_RI,R,I,RI_OBJ,oper);  \
  Func2_CXSC_X(gap_name##_CI,C,I,CI_OBJ,oper);

cxsc::complex pow(cxsc::real &a, cxsc::complex &b) { return pow(cxsc::complex(a),b); }
cxsc::interval pow(cxsc::real &a, cxsc::interval &b) { return pow(cxsc::interval(a),b); }
cxsc::cinterval pow(cxsc::real &a, cxsc::cinterval &b) { return pow(cxsc::cinterval(a),b); }
cxsc::cinterval pow(cxsc::complex &a, cxsc::interval &b) { return pow(cxsc::cinterval(a),cxsc::cinterval(b)); }
cxsc::cinterval pow(cxsc::complex &a, cxsc::cinterval &b) { return pow(cxsc::cinterval(a),b); }
cxsc::interval pow(cxsc::interval &a, cxsc::real &b) { return pow(a,cxsc::interval(b)); }
cxsc::cinterval pow(cxsc::interval &a, cxsc::complex &b) { return pow(cxsc::cinterval(a),cxsc::cinterval(b)); }
cxsc::cinterval pow(cxsc::interval &a, cxsc::cinterval &b) { return pow(cxsc::cinterval(a),b); }
cxsc::cinterval pow(cxsc::cinterval &a, cxsc::real &b) { return pow(a,cxsc::cinterval(b)); }
cxsc::cinterval pow(cxsc::cinterval &a, cxsc::complex &b) { return pow(a,cxsc::cinterval(b)); }

Func2_CXSC(SUM_CXSC,operator +);
Func2_CXSC(DIFF_CXSC,operator -);
Func2_CXSC(PROD_CXSC,operator *);
Func2_CXSC(QUO_CXSC,operator /);
Func2_CXSC(POW_CXSC,pow);

#define Func2a_CXSC_X_X(gap_name,c,i,getf,getg,oper) \
  static Obj gap_name(Obj self, Obj f, Obj g)  \
  {       \
    return OBJ_##c##I(oper(getf(f),getg(g)));  \
  }

#define Func2a_CXSC_X(gap_name,c,i,get,oper)  \
  Func2a_CXSC_X_X(gap_name##_RP,c,i,get,RP_OBJ,oper); \
  Func2a_CXSC_X_X(gap_name##_CP,C,i,get,CP_OBJ,oper); \
  Func2a_CXSC_X_X(gap_name##_RI,c,I,get,RI_OBJ,oper); \
  Func2a_CXSC_X_X(gap_name##_CI,C,I,get,CI_OBJ,oper); \

#define Func2a_CXSC(gap_name,oper)   \
  Func2a_CXSC_X(gap_name##_RP,R,P,RP_OBJ,oper);  \
  Func2a_CXSC_X(gap_name##_CP,C,P,CP_OBJ,oper);  \
  Func2a_CXSC_X(gap_name##_RI,R,I,RI_OBJ,oper);  \
  Func2a_CXSC_X(gap_name##_CI,C,I,CI_OBJ,oper);

cxsc::interval operator &(cxsc::real &a, cxsc::real &b) {
  return cxsc::interval(cxsc::QuietNaN); }
cxsc::cinterval operator &(cxsc::real &a, cxsc::complex &b) {
  return cxsc::cinterval(cxsc::QuietNaN); }
cxsc::cinterval operator &(cxsc::complex &a, cxsc::real &b) {
  return cxsc::cinterval(cxsc::QuietNaN); }
cxsc::cinterval operator &(cxsc::complex &a, cxsc::complex &b) {
  return cxsc::cinterval(cxsc::QuietNaN); }

Func2a_CXSC(OR_CXSC,operator |);
Func2a_CXSC(AND_CXSC,operator &);

#define Func2b_CXSC_X_X(gap_name,c,i,getf,getg,oper) \
  static Obj gap_name(Obj self, Obj f, Obj g)  \
  {       \
    return getf(f) oper getg(g) ? True : False;  \
  }

#define Func2b_CXSC_X(gap_name,c,i,get,oper)  \
  Func2b_CXSC_X_X(gap_name##_RP,c,i,get,RP_OBJ,oper); \
  Func2b_CXSC_X_X(gap_name##_CP,C,i,get,CP_OBJ,oper); \
  Func2b_CXSC_X_X(gap_name##_RI,c,I,get,RI_OBJ,oper); \
  Func2b_CXSC_X_X(gap_name##_CI,C,I,get,CI_OBJ,oper); \

#define Func2b_CXSC(gap_name,oper)   \
  Func2b_CXSC_X(gap_name##_RP,R,P,RP_OBJ,oper);  \
  Func2b_CXSC_X(gap_name##_CP,C,P,CP_OBJ,oper);  \
  Func2b_CXSC_X(gap_name##_RI,R,I,RI_OBJ,oper);  \
  Func2b_CXSC_X(gap_name##_CI,C,I,CI_OBJ,oper);

bool operator == (cxsc::complex &a, cxsc::interval &b) {
  return cxsc::cinterval(a) == cxsc::cinterval(b); }
bool operator == (cxsc::interval &a, cxsc::complex &b) {
  return cxsc::cinterval(a) == cxsc::cinterval(b); }
bool operator < (cxsc::complex &a, cxsc::real &b) { return false; }
bool operator < (cxsc::complex &a, cxsc::complex &b) { return false; }
bool operator < (cxsc::real &a, cxsc::complex &b) { return false; }
bool operator < (cxsc::interval &a, cxsc::complex &b) { return false; }
bool operator < (cxsc::complex &a, cxsc::interval &b) { return false; }

Func2b_CXSC(EQ_CXSC,==);
Func2b_CXSC(LT_CXSC,<);

static Obj POWER_CXSC_RP (Obj self, Obj f, Obj g)
{
  TEST_IS_INTOBJ(POWER_CXSC_RP,g);
  TEST_IS_RP(POWER_CXSC_RP,f);
  return OBJ_RP(power(RP_OBJ(f), INT_INTOBJ(g)));
}

static Obj POWER_CXSC_CP (Obj self, Obj f, Obj g)
{
  TEST_IS_INTOBJ(POWER_CXSC_CP,g);
  TEST_IS_CP(POWER_CXSC_CP,f);
  return OBJ_CP(power(CP_OBJ(f), INT_INTOBJ(g)));
}

static Obj POWER_CXSC_RI (Obj self, Obj f, Obj g)
{
  TEST_IS_INTOBJ(POWER_CXSC_RI,g);
  TEST_IS_RI(POWER_CXSC_RI,f);
  return OBJ_RI(power(RI_OBJ(f), INT_INTOBJ(g)));
}

static Obj POWER_CXSC_CI (Obj self, Obj f, Obj g)
{
  TEST_IS_INTOBJ(POWER_CXSC_CI,g);
  TEST_IS_CI(POWER_CXSC_CI,f);
  return OBJ_CI(power(CI_OBJ(f), INT_INTOBJ(g)));
}

static Obj ROOT_CXSC_RP (Obj self, Obj f, Obj g)
{
  TEST_IS_INTOBJ(ROOT_CXSC_RP,g);
  TEST_IS_RP(ROOT_CXSC_RP,f);
  return OBJ_RP(sqrt(RP_OBJ(f), INT_INTOBJ(g)));
}

static Obj ROOT_CXSC_CP (Obj self, Obj f, Obj g)
{
  TEST_IS_INTOBJ(ROOT_CXSC_CP,g);
  TEST_IS_CP(ROOT_CXSC_CP,f);
  return OBJ_CP(sqrt(CP_OBJ(f), INT_INTOBJ(g)));
}

static Obj ROOT_CXSC_RI (Obj self, Obj f, Obj g)
{
  TEST_IS_INTOBJ(ROOT_CXSC_RI,g);
  TEST_IS_RI(ROOT_CXSC_RI,f);
  return OBJ_RI(sqrt(RI_OBJ(f), INT_INTOBJ(g)));
}

static Obj ROOT_CXSC_CI (Obj self, Obj f, Obj g)
{
  TEST_IS_INTOBJ(ROOT_CXSC_CI,g);
  TEST_IS_CI(ROOT_CXSC_CI,f);
  return OBJ_CI(sqrt(CI_OBJ(f), INT_INTOBJ(g)));
}

static inline real ldexp (real f, int s)
{
  real g = f; times2pown(g, s); return g;
}

static inline cxsc::complex ldexp (cxsc::complex f, int s)
{
  return cxsc::complex(ldexp(Re(f),s),ldexp(Im(f),s));
}
static inline interval ldexp (interval f, int s)
{
  return interval(ldexp(Inf(f),s),ldexp(Sup(f),s));
}
static inline cinterval ldexp (cinterval f, int s)
{
  return cinterval(ldexp(Re(f),s),ldexp(Re(f),s));
}

static Obj LDEXP_CXSC_RP (Obj self, Obj f, Obj i)
{
  TEST_IS_INTOBJ(LDEXP_CXSC_RP,i);
  TEST_IS_RP(LDEXP_CXSC_RP,f);
  return OBJ_RP(ldexp(RP_OBJ(f), INT_INTOBJ(i)));
}

static Obj LDEXP_CXSC_CP (Obj self, Obj f, Obj i)
{
  TEST_IS_INTOBJ(LDEXP_CXSC_CP,i);
  TEST_IS_CP(LDEXP_CXSC_CP,f);
  return OBJ_CP(ldexp(CP_OBJ(f),INT_INTOBJ(i)));
}

static Obj LDEXP_CXSC_RI (Obj self, Obj f, Obj i)
{
  TEST_IS_INTOBJ(LDEXP_CXSC_RI,i);
  TEST_IS_RI(LDEXP_CXSC_RI,f);
  return OBJ_RI(ldexp(RI_OBJ(f),INT_INTOBJ(i)));
}

static Obj LDEXP_CXSC_CI (Obj self, Obj f, Obj i)
{
  TEST_IS_INTOBJ(LDEXP_CXSC_CI,i);
  TEST_IS_CI(LDEXP_CXSC_CI,f);
  return OBJ_CI(ldexp(CI_OBJ(f),INT_INTOBJ(i)));
}

static Obj FREXP_CXSC_RP (Obj self, Obj f)
{
  TEST_IS_RP(FREXP_CXSC_RP,f);
  Obj l = NEW_PLIST(T_PLIST,2);
  SET_ELM_PLIST(l,1,OBJ_RP(mant(RP_OBJ(f))));
  SET_ELM_PLIST(l,2,INTOBJ_INT(expo(RP_OBJ(f))));
  SET_LEN_PLIST(l,2);
  return l;
}

static Obj FREXP_CXSC_CP (Obj self, Obj f)
{
  TEST_IS_CP(FREXP_CXSC_CP,f);
  Obj l = NEW_PLIST(T_PLIST,2);
  int e0 = expo(Re(CP_OBJ(f))), e1 = expo(Im(CP_OBJ(f))), e = (e0 > e1 ? e0 : e1);
  SET_ELM_PLIST(l,1,OBJ_CP(ldexp(CP_OBJ(f),-e)));
  SET_ELM_PLIST(l,2,INTOBJ_INT(e));
  SET_LEN_PLIST(l,2);
  return l;
}

static Obj FREXP_CXSC_RI (Obj self, Obj f)
{
  TEST_IS_RI(FREXP_CXSC_RI,f);
  Obj l = NEW_PLIST(T_PLIST,2);
  int e0 = expo(Inf(RI_OBJ(f))), e1 = expo(Sup(RI_OBJ(f))), e = (e0 > e1 ? e0 : e1);
  SET_ELM_PLIST(l,1,OBJ_RI(ldexp(RI_OBJ(f),-e)));
  SET_ELM_PLIST(l,2,INTOBJ_INT(e));
  SET_LEN_PLIST(l,2);
  return l;
}

static Obj FREXP_CXSC_CI (Obj self, Obj f)
{
  TEST_IS_CI(FREXP_CXSC_CI,f);
  Obj l = NEW_PLIST(T_PLIST,2);
  int e00 = expo(Inf(Re(CI_OBJ(f)))), e01 = expo(Sup(Re(CI_OBJ(f)))), e0 = (e00 > e01 ? e00 : e01);
  int e10 = expo(Inf(Im(CI_OBJ(f)))), e11 = expo(Sup(Im(CI_OBJ(f)))), e1 = (e10 > e11 ? e10 : e11);
  int e = (e0 > e1 ? e0 : e1);
  SET_ELM_PLIST(l,1,OBJ_CI(ldexp(CI_OBJ(f),-e)));
  SET_ELM_PLIST(l,2,INTOBJ_INT(e));
  SET_LEN_PLIST(l,2);
  return l;
}

static void put_real (Obj list, int pos, cxsc::real f)
{
  SET_ELM_PLIST(list,pos,INTOBJ_INT(0));
  if (f == 0.0) {
    if (1.0/f > 0.0)
      SET_ELM_PLIST(list,pos+1,INTOBJ_INT(0));
    else
      SET_ELM_PLIST(list,pos+1,INTOBJ_INT(1));
    return;
  }
  if (IsInfinity(f)) {
    if (f > 0.0)
      SET_ELM_PLIST(list,pos+1,INTOBJ_INT(2));
    else
      SET_ELM_PLIST(list,pos+1,INTOBJ_INT(3));
    return;
  }
  if (IsQuietNaN(f)) {
    SET_ELM_PLIST(list,pos+1,INTOBJ_INT(4));
    return;
  }

  cxsc::real m = mant(f);
  cxsc::times2pown(m,26);
  int m0 = _double(m);
  Obj gapm = INTOBJ_INT(m0);
  m -= m0;
  cxsc::times2pown(m,27);
  gapm = SumInt(ProdInt(gapm,INTOBJ_INT(1 << 27)),INTOBJ_INT(_double(m)));

  while (!INT_INTOBJ(RemInt(gapm,INTOBJ_INT(2))))
    gapm = QuoInt(gapm,INTOBJ_INT(2));

  SET_ELM_PLIST(list,pos,gapm);
  SET_ELM_PLIST(list,pos+1,INTOBJ_INT(expo(f)));
}

static Obj EXTREPOFOBJ_CXSC_RP(Obj self, Obj f)
{
  TEST_IS_RP(EXTREPOBJOBJ_CXSC_RP,f);
  Obj l = NEW_PLIST(T_PLIST,2);
  SET_LEN_PLIST(l,2);
  put_real (l,1,RP_OBJ(f));
  return l;
}

static Obj EXTREPOFOBJ_CXSC_RI(Obj self, Obj f)
{
  TEST_IS_RI(EXTREPOBJOBJ_CXSC_RI,f);
  Obj l = NEW_PLIST(T_PLIST,4);
  SET_LEN_PLIST(l,4);
  put_real (l,1,Inf(RI_OBJ(f)));
  put_real (l,3,Sup(RI_OBJ(f)));
  return l;
}

static Obj EXTREPOFOBJ_CXSC_CP(Obj self, Obj f)
{
  TEST_IS_CP(EXTREPOBJOBJ_CXSC_CP,f);
  Obj l = NEW_PLIST(T_PLIST,4);
  SET_LEN_PLIST(l,4);
  put_real (l,1,Re(CP_OBJ(f)));
  put_real (l,3,Im(CP_OBJ(f)));
  return l;
}

static Obj EXTREPOFOBJ_CXSC_CI(Obj self, Obj f)
{
  TEST_IS_CI(EXTREPOBJOBJ_CXSC_CI,f);
  Obj l = NEW_PLIST(T_PLIST,8);
  SET_LEN_PLIST(l,8);
  put_real (l,1,InfRe(CI_OBJ(f)));
  put_real (l,3,SupRe(CI_OBJ(f)));
  put_real (l,5,InfIm(CI_OBJ(f)));
  put_real (l,7,SupIm(CI_OBJ(f)));
  return l;
}

static cxsc::real get_real (Obj l, int pos)
{
  if (LEN_PLIST(l) < pos+1)
    ErrorQuit("OBJBYEXTREP: length of argument must be at least %d", pos+1,0);

  Obj mant = ELM_PLIST(l,pos), expobj = ELM_PLIST(l,pos+1);

  if (!IS_INTOBJ(expobj) || !(IS_INTOBJ(mant) || TNUM_OBJ(mant)==T_INTPOS || TNUM_OBJ(mant)==T_INTNEG))
    ErrorQuit("OBJBYEXTREP: argument must be a list of integers", 0,0);

  int exp = INT_INTOBJ(expobj);

  if (mant == INTOBJ_INT(0))
    switch (exp) {
    case 0: return 0.0;
    case 1: return 1.0 / (-1.0 / 0.0);
    case 2: return 1.0 / 0.0;
    case 3: return -1.0 / 0.0;
    case 4: return cxsc::QuietNaN;
    }

  cxsc::real m = Double(INT_INTOBJ(RemInt(mant,INTOBJ_INT(1 << 27))));
  cxsc::times2pown(m,-27);
  m += Double(INT_INTOBJ(QuoInt(mant,INTOBJ_INT(1 << 27))));
  cxsc::times2pown(m,exp+27-INT_INTOBJ(CALL_1ARGS(GAPLog2Int,mant)));
  return m;
}

static Obj OBJBYEXTREP_CXSC_RP(Obj self, Obj l)
{
  return OBJ_RP(get_real(l,1));
}

static cxsc::interval get_interval(Obj l, int pos)
{
  return interval(get_real(l,pos),get_real(l,pos+2));
}

static Obj OBJBYEXTREP_CXSC_RI(Obj self, Obj l)
{
  return OBJ_RI(get_interval(l,1));
}

static Obj OBJBYEXTREP_CXSC_CP(Obj self, Obj l)
{
  return OBJ_CP(cxsc::complex(get_real(l,1),get_real(l,3)));
}

static Obj OBJBYEXTREP_CXSC_CI(Obj self, Obj l)
{
  return OBJ_CI(cxsc::cinterval(get_interval(l,1),get_interval(l,5)));
}

static Obj BLOW_CXSC_RI (Obj self, Obj f, Obj g)
{
  TEST_IS_RP(BLOW_CXSC_RI,g);
  TEST_IS_RI(BLOW_CXSC_RI,f);
  return OBJ_RI(Blow(RI_OBJ(f), RP_OBJ(g)));
}

static Obj BLOW_CXSC_CI (Obj self, Obj f, Obj g)
{
  TEST_IS_RP(BLOW_CXSC_CI,g);
  TEST_IS_CI(BLOW_CXSC_CI,f);
  return OBJ_CI(Blow(CI_OBJ(f), RP_OBJ(g)));
}

static Obj DISJOINT_CXSC_RI_RI (Obj self, Obj f, Obj g)
{
  TEST_IS_RI(DISJOINT_CXSC_RI_RI,f);
  TEST_IS_RI(DISJOINT_CXSC_RI_RI,g);
  return Disjoint(RI_OBJ(f),RI_OBJ(g)) ? True : False;
}

static Obj DISJOINT_CXSC_CI_CI (Obj self, Obj f, Obj g)
{
  TEST_IS_CI(DISJOINT_CXSC_CI_CI,f);
  TEST_IS_CI(DISJOINT_CXSC_CI_CI,g);
  return Disjoint(CI_OBJ(f),CI_OBJ(g)) ? True : False;
}

static Obj IN_CXSC_RP_RI (Obj self, Obj f, Obj g)
{
  TEST_IS_RP(IN_CXSC_RP_RI,f);
  TEST_IS_RI(IN_CXSC_RP_RI,g);
  return in(RP_OBJ(f),RI_OBJ(g)) ? True : False;
}

static Obj IN_CXSC_RI_RI (Obj self, Obj f, Obj g)
{
  TEST_IS_RI(IN_CXSC_RI_RI,f);
  TEST_IS_RI(IN_CXSC_RI_RI,g);
  return in(RI_OBJ(f),RI_OBJ(g)) ? True : False;
}

bool in (cxsc::complex &a, cxsc::cinterval &b) {
  return in(Re(a),Re(b)) && in(Im(a),Im(b));
}

static Obj IN_CXSC_RP_CI (Obj self, Obj f, Obj g)
{
  TEST_IS_RP(IN_CXSC_RP_CI,f);
  TEST_IS_CI(IN_CXSC_RP_CI,g);
  return in(_cinterval(RP_OBJ(f)),CI_OBJ(g)) ? True : False;
}

static Obj IN_CXSC_CP_CI (Obj self, Obj f, Obj g)
{
  TEST_IS_CP(IN_CXSC_CP_CI,f);
  TEST_IS_CI(IN_CXSC_CP_CI,g);
  return in(_cinterval(CP_OBJ(f)),CI_OBJ(g)) ? True : False;
}

static Obj IN_CXSC_RI_CI (Obj self, Obj f, Obj g)
{
  TEST_IS_RI(IN_CXSC_RI_CI,f);
  TEST_IS_CI(IN_CXSC_RI_CI,g);
  return in(_cinterval(RI_OBJ(f)),CI_OBJ(g)) ? True : False;
}

static Obj IN_CXSC_CI_CI (Obj self, Obj f, Obj g)
{
  TEST_IS_CI(IN_CXSC_CI_CI,f);
  TEST_IS_CI(IN_CXSC_CI_CI,g);
  return in(CI_OBJ(f),CI_OBJ(g)) ? True : False;
}

static Obj ATAN2_CXSC_RP_RP (Obj self, Obj f, Obj g)
{
  TEST_IS_RP(ATAN2_CXSC_RP_RP,f);
  TEST_IS_RP(ATAN2_CXSC_RP_RP,g);
  return OBJ_RP(std::atan2(_double(RP_OBJ(f)),_double(RP_OBJ(g))));
}

static Obj ATAN2_CXSC_CP (Obj self, Obj f)
{
  TEST_IS_CP(ATAN2_CXSC_CP,f);
  cxsc::complex z = CP_OBJ(f);
  return OBJ_RP(std::atan2(_double(Im(z)),_double(Re(z))));
}

static Obj ATAN2_CXSC_CI (Obj self, Obj f)
{
  TEST_IS_CI(ATAN2_CXSC_CI,f);
  return OBJ_RI(Im(cxsc::ln(CI_OBJ(f))));
}

static Obj HYPOT_CXSC_RP2 (Obj self, Obj f, Obj g)
{
  TEST_IS_RP(HYPOT_CXSC_RP2,f);
  TEST_IS_RP(HYPOT_CXSC_RP2,g);
  return OBJ_RP(cxsc::sqrtx2y2(RP_OBJ(f),RP_OBJ(g)));
}

/****************************************************************
* export functions
****************************************************************/
static StructGVarFunc GVarFuncs [] = {
  Inc1_CXSC_arg(CXSC_INT,"int"),
  Inc1_CXSC_arg(CXSC_IEEE754,"float"),
  Inc1_CXSC_arg(CXSC_NEWCONSTANT,"int"),
  Inc1_CXSC_arg(RP_CXSC_STRING,"string"),
  Inc1_CXSC_arg(RI_CXSC_STRING,"string"),
  Inc1_CXSC_arg(RI_CXSC_RP,"cxsc::rp"),
  Inc2_CXSC_arg(RI_CXSC_RP_RP,"cxsc::rp,cxsc::rp"),
  Inc1_CXSC_arg(CP_CXSC_STRING,"string"),
  Inc2_CXSC_arg(CP_CXSC_RP_RP,"cxsc::rp,cxsc::rp"),
  Inc1_CXSC_arg(CP_CXSC_RP,"cxsc::rp"),
  Inc1_CXSC_arg(CI_CXSC_STRING,"string"),
  Inc2_CXSC_arg(CI_CXSC_RI_RI,"cxsc::ri,cxsc::ri"),
  Inc1_CXSC_arg(CI_CXSC_CP,"cxsc::cp"),
  Inc2_CXSC_arg(CI_CXSC_CP_CP,"cxsc::cp,cxsc::cp"),
  Inc1_CXSC_arg(INT_CXSC,"cxsc::rp"),

  Inc1a_CXSC(REAL_CXSC),
  Inc1a_CXSC(IMAG_CXSC),
  Inc1a_CXSC(NORM_CXSC),
  Inc1a_CXSC(CONJ_CXSC),

  Inc1b_CXSC(DIAM_CXSC),
  Inc1b_CXSC(DIAM_REL_CXSC),
  Inc1b_CXSC(ISEMPTY_CXSC),
  Inc1b_CXSC(MID_CXSC),
  Inc1b_CXSC(INF_CXSC),
  Inc1b_CXSC(SUP_CXSC),

  Inc1_CXSC(AINV_CXSC),
  Inc1_CXSC(INV_CXSC),
  Inc1_CXSC(SIN_CXSC),
  Inc1_CXSC(COS_CXSC),
  Inc1_CXSC(TAN_CXSC),
  Inc1_CXSC(COT_CXSC),
  Inc1_CXSC(SINH_CXSC),
  Inc1_CXSC(COSH_CXSC),
  Inc1_CXSC(TANH_CXSC),
  Inc1_CXSC(COTH_CXSC),
  Inc1_CXSC(ASIN_CXSC),
  Inc1_CXSC(ACOS_CXSC),
  Inc1_CXSC(ATAN_CXSC),
  Inc1_CXSC(ACOT_CXSC),
  Inc1_CXSC(ASINH_CXSC),
  Inc1_CXSC(ACOSH_CXSC),
  Inc1_CXSC(ATANH_CXSC),
  Inc1_CXSC(ACOTH_CXSC),
  Inc1_CXSC(SQR_CXSC),
  Inc1_CXSC(SQRT_CXSC),
  Inc1_CXSC(EXP_CXSC),
  Inc1_CXSC(EXPM1_CXSC),
  Inc1_CXSC(FREXP_CXSC),
  Inc1_CXSC(EXTREPOFOBJ_CXSC),
  Inc1_CXSC(OBJBYEXTREP_CXSC),
  Inc1_CXSC(LOG_CXSC),
  Inc1_CXSC(LOG1P_CXSC),
  Inc1_CXSC(LOG2_CXSC),
  Inc1_CXSC(LOG10_CXSC),
  Inc1_CXSC(ABS_CXSC),
  Inc1_CXSC(ISNAN_CXSC),
  Inc1_CXSC(ISXINF_CXSC),
  Inc1_CXSC(ISPINF_CXSC),
  Inc1_CXSC(ISNINF_CXSC),
  Inc1_CXSC(ISZERO_CXSC),
  Inc1_CXSC(ISONE_CXSC),
  Inc1_CXSC(ISNUMBER_CXSC),
  Inc1_CXSC_arg(ATAN2_CXSC_CP,"cxsc::cp"),
  Inc1_CXSC_arg(ATAN2_CXSC_CI,"cxsc::ci"),

  Inc2_CXSC_arg(HYPOT_CXSC_RP2,"x, y"),
  Inc1_CXSC_arg(SIGN_CXSC_RP,"cxsc::rp"),
  Inc1_CXSC_arg(SIGN_CXSC_RI,"cxsc::ri"),
  Inc2_CXSC(SUM_CXSC),
  Inc2_CXSC(DIFF_CXSC),
  Inc2_CXSC(PROD_CXSC),
  Inc2_CXSC(QUO_CXSC),
  Inc2_CXSC(POW_CXSC),
  Inc2_CXSC(OR_CXSC),
  Inc2_CXSC(AND_CXSC),
  Inc2_CXSC(EQ_CXSC),
  Inc2_CXSC(LT_CXSC),
  Inc2_CXSC_arg(ATAN2_CXSC_RP_RP,"cxsc::rp,cxsc::rp"),
  Inc2_CXSC_arg(POWER_CXSC_RP,"cxsc::rp,int"),
  Inc2_CXSC_arg(POWER_CXSC_CP,"cxsc::cp,int"),
  Inc2_CXSC_arg(POWER_CXSC_RI,"cxsc::ri,int"),
  Inc2_CXSC_arg(POWER_CXSC_CI,"cxsc::ci,int"),
  Inc2_CXSC_arg(ROOT_CXSC_RP,"cxsc::rp,int"),
  Inc2_CXSC_arg(ROOT_CXSC_CP,"cxsc::cp,int"),
  Inc2_CXSC_arg(ROOT_CXSC_RI,"cxsc::ri,int"),
  Inc2_CXSC_arg(ROOT_CXSC_CI,"cxsc::ci,int"),
  Inc2_CXSC_arg(LDEXP_CXSC_RP,"cxsc::rp,int"),
  Inc2_CXSC_arg(LDEXP_CXSC_CP,"cxsc::cp,int"),
  Inc2_CXSC_arg(LDEXP_CXSC_RI,"cxsc::ri,int"),
  Inc2_CXSC_arg(LDEXP_CXSC_CI,"cxsc::ci,int"),
  Inc2_CXSC_arg(BLOW_CXSC_RI,"cxsc::ri,cxsc::rp"),
  Inc2_CXSC_arg(BLOW_CXSC_CI,"cxsc::ci,cxsc::rp"),
  Inc2_CXSC_arg(DISJOINT_CXSC_RI_RI,"cxsc::ri,cxsc::ri"),
  Inc2_CXSC_arg(DISJOINT_CXSC_CI_CI,"ccxsc::ri,ccxsc::ri"),
  Inc2_CXSC_arg(IN_CXSC_RP_RI,"cxsc::rp,cxsc::ri"),
  Inc2_CXSC_arg(IN_CXSC_RI_RI,"cxsc::ri,cxsc::ri"),
  Inc2_CXSC_arg(IN_CXSC_RP_CI,"cxsc::rp,ccxsc::ci"),
  Inc2_CXSC_arg(IN_CXSC_CP_CI,"cxsc::cp,ccxsc::ci"),
  Inc2_CXSC_arg(IN_CXSC_RI_CI,"ccxsc::ri,ccxsc::ci"),
  Inc2_CXSC_arg(IN_CXSC_CI_CI,"ccxsc::ci,ccxsc::ci"),
  Inc2_CXSC_arg(ROOTPOLY_CXSC,"cxsc::c*[],bool"),
  Inc2_CXSC_arg(EVALPOLY_CXSC,"cxsc::rp[],cxsc::rp"),
  { "STRING_CXSC", 3, "cxsc:**,int,int", (ObjFunc) STRING_CXSC, "src/cxsc_float.c:STRING_CXSC" },
  {0}
};

void cxsc_terminate (void)
{
  ErrorQuit("cxsc::terminate: I'll be back!", 0,0);
}

extern "C" {
int InitCXSCKernel (void)
{
  InitHdlrFuncsFromTable (GVarFuncs);

  ImportGVarFromLibrary ("TYPE_CXSC_RP", &TYPE_CXSC_RP);
  ImportGVarFromLibrary ("TYPE_CXSC_CP", &TYPE_CXSC_CP);
  ImportGVarFromLibrary ("TYPE_CXSC_RI", &TYPE_CXSC_RI);
  ImportGVarFromLibrary ("TYPE_CXSC_CI", &TYPE_CXSC_CI);

  ImportGVarFromLibrary ("IsCXSCReal", &IS_CXSC_RP);
  ImportGVarFromLibrary ("IsCXSCComplex", &IS_CXSC_CP);
  ImportGVarFromLibrary ("IsCXSCInterval", &IS_CXSC_RI);
  ImportGVarFromLibrary ("IsCXSCBox", &IS_CXSC_CI);

  ImportGVarFromLibrary ("CXSCFloatsFamily", &FAMILY_CXSC);

  ImportFuncFromLibrary ("Log2Int", &GAPLog2Int);

  set_terminate (cxsc_terminate);
  return 0;
}

int InitCXSCLibrary (void)
{
  InitGVarFuncsFromTable (GVarFuncs);
  return 0;
}
}

/****************************************************************************
**
*E  cxsc.c. . . . . . . . . . . . . . . . . . . . . . . . . . . . . ends here
*/

quality97%

¤ Dauer der Verarbeitung: 0.13 Sekunden (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.