Quelle mpc.c Sprache: C

/****************************************************************************
**
*W  mpc.c                       GAP source                  Laurent Bartholdi
**
*Y  Copyright (C) 2008-2012 Laurent Bartholdi
**
**  This file contains the functions for the float package.
**  complex floats are implemented using the MPC 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 <stdio.h>

#include "floattypes.h"
#include <mpc.h>
#include <mp_poly.h>

/****************************************************************
* mpc's are stored as follows:
* +----------+-----------------------------------------+---------------------+
* | TYPE_MPC |               __mpc_struct              |    __mp_limb_t[]    |
* |          | __mpfr_struct real          imag        | limbr ... limbi ... |
* |          | prec exp sign mant   prec exp sign mant |                     |
* +----------+-----------------------------------------+---------------------+
*                               \___________________\____^         ^
*                                                    \____________/
* it is assumed that the real and imaginary mpfr's are allocated with the
* same precision
****************************************************************/
#define MPC_OBJ(obj) ((mpc_ptr) (ADDR_OBJ(obj)+1))
#define REMANTISSA_MPC(p) ((mp_limb_t *) (p+1))
#define IMMANTISSA_MPC(p) (REMANTISSA_MPC(p)+(mpc_get_prec(p)+GMP_NUMB_BITS-1)/GMP_NUMB_BITS)

static inline mpc_ptr GET_MPC(Obj obj) {
  mpc_ptr p = MPC_OBJ(obj);
  mpfr_custom_move (p->re, REMANTISSA_MPC(p));
  mpfr_custom_move (p->im, IMMANTISSA_MPC(p));
  return p;
}

/****************************************************************************
**
*F  allocate new object, initialize to NaN
**
*/
Obj TYPE_MPC;

static inline Obj NEW_MPC( mp_prec_t prec )
{
  Obj f = NEW_DATOBJ(sizeof(__mpc_struct)+2*mpfr_custom_get_size(prec), TYPE_MPC);
  mpc_ptr p = MPC_OBJ(f);
  mpfr_custom_init_set(p->re, MPFR_NAN_KIND, 0, prec, REMANTISSA_MPC(p));
  mpfr_custom_init_set(p->im, MPFR_NAN_KIND, 0, prec, IMMANTISSA_MPC(p));
  return f;
}

/****************************************************************************
**
*F Func1_MPC( <float> ) . . . . . . . . . . . . . . . .1-argument functions
**
*/
#define Func1_MPC(name,mpc_name)   \
  static Obj name(Obj self, Obj f)   \
  {              \
    mp_prec_t prec = mpc_get_prec(GET_MPC(f));         \
    Obj g = NEW_MPC(prec);    \
    mpc_name (MPC_OBJ(g), GET_MPC(f), MPC_RNDNN); \
    return g;      \
  }
#define Func1_MPFRMPC(name,mpc_name)   \
  static Obj name(Obj self, Obj f)   \
  {              \
    mp_prec_t prec = mpc_get_prec(GET_MPC(f));         \
    Obj g = NEW_MPFR(prec);    \
    mpc_name (MPFR_OBJ(g), GET_MPC(f), MPC_RNDNN); \
    return g;      \
  }
#define Inc1_MPC_arg(name,arg)  \
  { #name, 1, arg, (ObjFunc)(void *)(name), "src/mpc.c:" #name }
#define Inc1_MPC(name) Inc1_MPC_arg(name,"complex")

Func1_MPC(PROJ_MPC,mpc_proj);
Func1_MPC(EXP_MPC,mpc_exp);
Func1_MPC(LOG_MPC,mpc_log);
Func1_MPC(CONJ_MPC,mpc_conj);
Func1_MPC(AINV_MPC,mpc_neg);
Func1_MPC(SQRT_MPC,mpc_sqrt);
Func1_MPC(SQR_MPC,mpc_sqr);

Func1_MPC(SIN_MPC,mpc_sin);
Func1_MPC(COS_MPC,mpc_cos);
Func1_MPC(TAN_MPC,mpc_tan);
Func1_MPC(SINH_MPC,mpc_sinh);
Func1_MPC(COSH_MPC,mpc_cosh);
Func1_MPC(TANH_MPC,mpc_tanh);
Func1_MPC(ASIN_MPC,mpc_asin);
Func1_MPC(ACOS_MPC,mpc_acos);
Func1_MPC(ATAN_MPC,mpc_atan);
Func1_MPC(ASINH_MPC,mpc_asinh);
Func1_MPC(ACOSH_MPC,mpc_acosh);
Func1_MPC(ATANH_MPC,mpc_atanh);

Func1_MPFRMPC(ABS_MPC,mpc_abs);
Func1_MPFRMPC(NORM_MPC,mpc_norm);
Func1_MPFRMPC(REAL_MPC,mpc_real);
Func1_MPFRMPC(IMAG_MPC,mpc_imag);
Func1_MPFRMPC(ARG_MPC,mpc_arg);

int mpc_nan_p(mpc_ptr c) { return mpfr_nan_p(c->re) || mpfr_nan_p(c->im); }
int mpc_inf_p(mpc_ptr c) { return mpfr_inf_p(c->re) || mpfr_inf_p(c->im); }
int mpc_zero_p(mpc_ptr c) { return mpfr_zero_p(c->re) && mpfr_zero_p(c->im); }
int mpc_number_p(mpc_ptr c) { return mpfr_number_p(c->re) && mpfr_number_p(c->im); }

static Obj ISNAN_MPC(Obj self, Obj f)
{
  return mpc_nan_p(GET_MPC(f)) ? True : False;
}

static Obj ISINF_MPC(Obj self, Obj f)
{
  return mpc_inf_p(GET_MPC(f)) ? True : False;
}

static Obj ISZERO_MPC(Obj self, Obj f)
{
  return mpc_zero_p(GET_MPC(f)) ? True : False;
}

static Obj ISNUMBER_MPC(Obj self, Obj f)
{
  return mpc_number_p(GET_MPC(f)) ? True : False;
}

static Obj FREXP_MPC(Obj self, Obj f)
{
  mp_prec_t prec = mpc_get_prec(GET_MPC(f));
  Obj g = NEW_MPC(prec);
  mpc_set (MPC_OBJ(g), GET_MPC(f), MPC_RNDNN);
  mp_exp_t ere = mpfr_get_exp(MPC_OBJ(f)->re),
    eim = mpfr_get_exp(MPC_OBJ(f)->im);
  mp_exp_t e = ere > eim ? ere : eim;
  mpfr_set_exp(MPC_OBJ(g)->re, ere-e);
  mpfr_set_exp(MPC_OBJ(g)->im, eim-e);
  Obj l = NEW_PLIST(T_PLIST,2);
  SET_ELM_PLIST(l,1,g);
  SET_ELM_PLIST(l,2, ObjInt_Int(e));
  SET_LEN_PLIST(l,2);
  return l;
}

static Obj LDEXP_MPC(Obj self, Obj f, Obj exp)
{
  mp_exp_t e = Int_ObjInt(exp);
  mp_prec_t prec = mpc_get_prec(GET_MPC(f));
  Obj g = NEW_MPC(prec);
  mpfr_mul_2si (MPC_OBJ(g)->re, GET_MPC(f)->re, e, GMP_RNDN);
  mpfr_mul_2si (MPC_OBJ(g)->im, MPC_OBJ(f)->im, e, GMP_RNDN);
  return g;
}

static Obj EXTREPOFOBJ_MPC(Obj self, Obj f)
{
  mp_prec_t prec = mpc_get_prec(GET_MPC(f));
  int i;
  Obj l = NEW_PLIST(T_PLIST,4);
  SET_LEN_PLIST(l,4);
  Obj g = NEW_MPFR(prec);

  if (mpc_zero_p(GET_MPC(f))) {
    SET_ELM_PLIST(l,1, INTOBJ_INT(0));
    SET_ELM_PLIST(l,2, INTOBJ_INT(0));
    return l;
  }
  if (!mpc_number_p(GET_MPC(f))) {
    SET_ELM_PLIST(l,1, INTOBJ_INT(0));
    if (mpc_nan_p(MPC_OBJ(f)))
      SET_ELM_PLIST(l,2, INTOBJ_INT(4));
    else if (mpc_inf_p(MPC_OBJ(f)))
      SET_ELM_PLIST(l,2, INTOBJ_INT(2));
    return l;
  }

  mpz_t z;
  mpz_init2 (z, prec);

  for (i = 0; i < 2; i++) {
    if (i == 0)
      mpfr_set (MPFR_OBJ(g), GET_MPC(f)->re, GMP_RNDN);
    else
      mpfr_set (MPFR_OBJ(g), GET_MPC(f)->im, GMP_RNDN);

    mp_exp_t e = mpfr_get_exp(MPFR_OBJ(g));
    mpfr_set_exp(MPFR_OBJ(g), prec);
    mpfr_get_z (z, MPFR_OBJ(g), GMP_RNDZ);
    Obj ig = INT_mpz(z);

    SET_ELM_PLIST(l,2*i+1, ig);
    SET_ELM_PLIST(l,2*i+2, ObjInt_Int(e));
  }
  mpz_clear(z);

  return l;
}

#pragma GCC diagnostic ignored "-Wuninitialized"
static Obj OBJBYEXTREP_MPC(Obj self, Obj list)
{
  int i;
  mp_prec_t prec = 0;

  if (LEN_LIST(list) != 4) {
    ErrorMayQuit("OBJBYEXTREP_MPC: object must be a list of length 4, not a %s",
         (Int)TNAM_OBJ(list),0);
  }

  for (i = 0; i < 4; i += 2) {
    Obj m = ELM_PLIST(list,i+1);
    mp_prec_t newprec;
    if (IS_INTOBJ(m))
      newprec = 8*sizeof(long);
    else if (TNUM_OBJ(m) == T_INTPOS || TNUM_OBJ(m) == T_INTNEG)
      newprec = 8*sizeof(mp_limb_t)*SIZE_INT(m);
    else
      ErrorQuit("OBJBYEXTREP_MPC: invalid argument %d", i+1,0);
    if (newprec > prec)
      prec = newprec;
  }
  Obj f = NEW_MPC(prec);

  for (i = 0; i < 4; i++) {
    Obj m = ELM_PLIST(list,i+1);
    Int iarg;
    mpz_ptr zarg;
    int argtype;

    if (IS_INTOBJ(m))
      iarg = INT_INTOBJ(m), argtype = 0;
    else
      zarg = mpz_MPZ(MPZ_LONGINT(m)), argtype = 1;

    if (i & 1 && argtype) /* exponent, must be small int */
      iarg = mpz_get_si(zarg), argtype = 0;

    mpfr_ptr reim;
    if (i < 2)
      reim = GET_MPC(f)->re;
    else
      reim = GET_MPC(f)->im;

    if (i & 1)
      mpfr_set_exp(reim, iarg);
    else if (argtype)
      mpfr_set_z(reim, zarg, GMP_RNDN);
    else {
      if (iarg == 0) { /* special case */
switch (INT_INTOBJ(ELM_PLIST(list,i+2))) {
case 0: /* 0 */
case 1: /* -0 */
   mpfr_set_si(reim, 0, GMP_RNDN); break;
case 2: /* inf */
case 3: /* -inf */
   mpfr_set_inf (reim, 1); break;
case 4: /* nan */
case 5: /* -nan */
   mpfr_set_nan (reim); break;
default:
   ErrorQuit("OBJBYEXTREP_MPC: invalid argument [%d,%d]",
      iarg, INT_INTOBJ(ELM_PLIST(list,i+2)));
}
i++; /* skip "exponent" */
continue;
      }
      mpfr_set_si(reim, iarg, GMP_RNDN);
    }
  }
  return f;
}

static Obj ZERO_MPC(Obj self, Obj f)
{
  mp_prec_t prec = mpc_get_prec(GET_MPC(f));
  Obj g = NEW_MPC(prec);
  mpc_set_ui(MPC_OBJ(g), 0, MPC_RNDNN);
  return g;
}

static Obj ONE_MPC(Obj self, Obj f)
{
  mp_prec_t prec = mpc_get_prec(GET_MPC(f));
  Obj g = NEW_MPC(prec);
  mpc_set_ui (MPC_OBJ(g), 1, MPC_RNDNN);
  return g;
}

static Obj MPC_MAKENAN(Obj self, Obj prec)
{
  TEST_IS_INTOBJ("MPC_MAKENAN",prec);

  Obj g = NEW_MPC(INT_INTOBJ(prec));
  mpfr_set_nan (MPC_OBJ(g)->re);
  mpfr_set_nan (MPC_OBJ(g)->im);
  return g;
}

static Obj MPC_MAKEINFINITY(Obj self, Obj prec)
{
  TEST_IS_INTOBJ("MPC_MAKEINFINITY",prec);

  int p = INT_INTOBJ(prec);
  Obj g = NEW_MPC(p < 0 ? -p : p);
  mpfr_set_inf (MPC_OBJ(g)->re, p);
  mpfr_set_inf (MPC_OBJ(g)->im, p);
  return g;
}

static Obj INV_MPC(Obj self, Obj f)
{
  Obj g = NEW_MPC(mpc_get_prec(GET_MPC(f)));
  mpc_ui_div (MPC_OBJ(g), 1, GET_MPC(f), MPC_RNDNN);
  return g;
}

static Obj MPC_INT(Obj self, Obj i)
{
  Obj g;
  if (IS_INTOBJ(i)) {
    g = NEW_MPC(8*sizeof(long));
    mpc_set_si(MPC_OBJ(g), INT_INTOBJ(i), MPC_RNDNN);
  } else {
    Obj f = MPZ_LONGINT(i);

    g = NEW_MPC(8*sizeof(mp_limb_t)*SIZE_INT(i));
    mpfr_set_z(MPC_OBJ(g)->re, mpz_MPZ(f), GMP_RNDN);
    mpfr_set_ui(MPC_OBJ(g)->im, 0, GMP_RNDN);
  }
  return g;
}

static Obj MPC_INTPREC(Obj self, Obj i, Obj prec)
{
  Obj g;
  TEST_IS_INTOBJ("MPC_INTPREC",prec);
  if (IS_INTOBJ(i)) {
    g = NEW_MPC(INT_INTOBJ(prec));
    mpc_set_si(MPC_OBJ(g), INT_INTOBJ(i), MPC_RNDNN);
  } else {
    Obj f = MPZ_LONGINT(i);
    g = NEW_MPC(INT_INTOBJ(prec));
    mpfr_set_z(MPC_OBJ(g)->re, mpz_MPZ(f), GMP_RNDN);
    mpfr_set_ui(MPC_OBJ(g)->im, 0, GMP_RNDN);
  }
  return g;
}

static Obj PREC_MPC(Obj self, Obj f)
{
  return INTOBJ_INT(mpc_get_prec(GET_MPC(f)));
}

static Obj MPC_MPCPREC(Obj self, Obj f, Obj prec)
{
  TEST_IS_INTOBJ("MPC_MPCPREC",prec);

  Obj g = NEW_MPC(INT_INTOBJ(prec));
  mpc_set (MPC_OBJ(g), GET_MPC(f), MPC_RNDNN);
  return g;
}

static Obj STRING_MPC(Obj self, Obj f, Obj digits)
{
  mp_prec_t prec = mpc_get_prec(GET_MPC(f));
  Obj str = NEW_STRING(2*(prec*302/1000+10)+3);
  int slen = 0, smid, n;

  TEST_IS_INTOBJ("STRING_MPC",digits);
  n = INT_INTOBJ(digits);
  if (n == 1) n = 2;

  char *c = CSTR_STRING(str);
  slen += PRINT_MPFR(c+slen, 0, n, GET_MPC(f)->re, GMP_RNDN);
  c[slen++] = '+';
  smid = slen;
  slen += PRINT_MPFR(c+slen, 0, n, GET_MPC(f)->im, GMP_RNDN);
  if (c[smid] == '-') {
    while (smid <= slen) { c[smid-1] = c[smid]; smid++; }
    slen--;
  }
  c[slen++] = 'i';
  c[slen] = 0;
  SET_LEN_STRING(str, slen);
  SHRINK_STRING(str);

  return str;
}

static Obj VIEWSTRING_MPC(Obj self, Obj f, Obj digits)
{
  mp_prec_t prec = mpc_get_prec(GET_MPC(f));
  Obj str = NEW_STRING(2*(prec*302/1000+10)+3);
  int slen = 0, n;

  TEST_IS_INTOBJ("VIEWSTRING_MPC",digits);
  n = INT_INTOBJ(digits);
  if (n == 1) n = 2;

  char *c = CSTR_STRING(str);
  if (mpc_inf_p(GET_MPC(f))) {
    strcat(c+slen, CSTR_STRING(FLOAT_INFINITY_STRING));
    slen += GET_LEN_STRING(FLOAT_INFINITY_STRING);
  } else if (mpc_nan_p(GET_MPC(f))) {
    c[slen++] = 'n';
    c[slen++] = 'a';
    c[slen++] = 'n';
  } else {
    slen += PRINT_MPFR(c+slen, 0, n, GET_MPC(f)->re, GMP_RNDN);
    Obj im = NEW_MPFR(prec);
    // Reassign in case GC occurred
    c = CSTR_STRING(str);
    mpfr_add(MPFR_OBJ(im), GET_MPC(f)->re, GET_MPC(f)->im, GMP_RNDN);
    mpfr_sub(MPFR_OBJ(im), MPFR_OBJ(im), GET_MPC(f)->re, GMP_RNDN); /* round off small im */
    if (!mpfr_zero_p(MPFR_OBJ(im))) {
      if (mpfr_sgn(MPFR_OBJ(im)) < 0)
c[slen++] = '-';
      else
c[slen++] = '+';
      mpfr_abs (MPFR_OBJ(im), GET_MPC(f)->im, GMP_RNDN);
      slen += PRINT_MPFR(c+slen, 0, n, MPFR_OBJ(im), GMP_RNDN);
      strcat(c+slen, CSTR_STRING(FLOAT_I_STRING));
      slen += GET_LEN_STRING(FLOAT_I_STRING);
    }
  }
  CSTR_STRING(str)[slen] = 0;
  SET_LEN_STRING(str, slen);
  SHRINK_STRING(str);

  return str;
}

static Obj MPC_STRING(Obj self, Obj s, Obj prec)
{
  TEST_IS_STRING(MPC_STRING, s);
  TEST_IS_INTOBJ("MPC_STRING",prec);
  int n = INT_INTOBJ(prec);
  if (n == 0)
    n = GET_LEN_STRING(s)*1000 / 301;

  Obj newg = NEW_MPFR(INT_INTOBJ(prec));
  Obj g = NEW_MPC(INT_INTOBJ(prec));
  // Force newg internal pointers to be updated
  GET_MPFR(newg);

  int sign = 1;
  mpc_set_ui(MPC_OBJ(g), 0, MPC_RNDNN);
  mpfr_ptr f = MPC_OBJ(g)->re;
  char *p = CSTR_STRING(s), *newp;

  for (;;) {
    switch (*p) {
    case '-':
    case '+':
    case 0:
      if (!mpfr_nan_p(MPFR_OBJ(newg))) { /* drop the last read float */
mpfr_add (f, f, MPFR_OBJ(newg), GMP_RNDN);
mpfr_set_nan (MPFR_OBJ(newg));
f = GET_MPC(g)->re;
sign = 1;
      }
      if (!*p)
return g;
      if (*p == '-')
sign = -sign;
    case '*': p++; break;
    case 'i':
    case 'I': if (f == GET_MPC(g)->re) {
f = GET_MPC(g)->im;
if (mpfr_nan_p(MPFR_OBJ(newg)))
   mpfr_set_si (MPFR_OBJ(newg), sign, GMP_RNDN); /* accept 'i' as '1*i' */
      } else return Fail;
      p++; break;
    default:
      mpfr_strtofr(MPFR_OBJ(newg), p, &newp, 10, GMP_RNDN);
      if (newp == p && f != GET_MPC(g)->im)
return Fail; /* no valid characters read */
      if (sign == -1)
mpfr_neg(MPFR_OBJ(newg), MPFR_OBJ(newg), GMP_RNDN);
      p = newp;
    }
  }
  return g;
}

static Obj MPC_MPFR(Obj self, Obj f)
{
  Obj g = NEW_MPC (mpfr_get_prec(GET_MPFR(f)));
  mpfr_set (MPC_OBJ(g)->re, GET_MPFR(f), GMP_RNDN);
  mpfr_set_ui (MPC_OBJ(g)->im, 0, GMP_RNDN);

  return g;
}

/****************************************************************************
**
*F Func2_MPC( <float>, <float> ) . . . . . . . . . . . 2-argument functions
**
*/
#define Func2_MPC(name,mpc_name)     \
  static Obj name##_MPC(Obj self, Obj fl, Obj fr)   \
  {         \
    mp_prec_t precl = mpc_get_prec(GET_MPC(fl)),   \
      precr = mpc_get_prec(GET_MPC(fr));    \
         \
    Obj g = NEW_MPC(precl > precr ? precl : precr);   \
    mpc_name (MPC_OBJ(g), GET_MPC(fl), GET_MPC(fr), MPC_RNDNN);  \
    return g;        \
  }
#define Func2_MPC_MPFR(name,mpc_name)     \
  static Obj name##_MPC_MPFR(Obj self, Obj fl, Obj fr)   \
  {         \
    mp_prec_t precl = mpc_get_prec(GET_MPC(fl)),   \
      precr = mpfr_get_prec(GET_MPFR(fr));    \
         \
    Obj g = NEW_MPC(precl > precr ? precl : precr);   \
    mpc_name (MPC_OBJ(g), GET_MPC(fl), GET_MPFR(fr), MPC_RNDNN); \
    return g;        \
  }
#define Func2_MPFR_MPC(name,mpc_name)     \
  static Obj name##_MPFR_MPC(Obj self, Obj fl, Obj fr)   \
  {         \
    mp_prec_t precl = mpfr_get_prec(GET_MPFR(fl)),   \
      precr = mpc_get_prec(GET_MPC(fr));    \
         \
    Obj g = NEW_MPC(precl > precr ? precl : precr);   \
    mpc_name (MPC_OBJ(g), GET_MPFR(fl), GET_MPC(fr), MPC_RNDNN); \
    return g;        \
  }
#define Inc2_MPC_arg(name,arg)   \
  { #name, 2, arg, (ObjFunc)(void *)(name), "src/mpc.c:" #name }
#define Inc2_MPC(name) Inc2_MPC_arg(name##_MPC,"complex, complex"), \
    Inc2_MPC_arg(name##_MPC_MPFR,"complex, real"),   \
    Inc2_MPC_arg(name##_MPFR_MPC,"real, complex")   \

Func2_MPC(SUM,mpc_add);
Func2_MPC(DIFF,mpc_sub);
Func2_MPC(PROD,mpc_mul);
Func2_MPC(QUO,mpc_div);
Func2_MPC(POW,mpc_pow);
Func2_MPC_MPFR(SUM,mpc_add_fr);
Func2_MPC_MPFR(DIFF,mpc_sub_fr);
Func2_MPC_MPFR(PROD,mpc_mul_fr);
Func2_MPC_MPFR(QUO,mpc_div_fr);
Func2_MPC_MPFR(POW,mpc_pow_fr);
#define mpc_fr_add(a,b,c,d) mpc_add_fr(a,c,b,d)
Func2_MPFR_MPC(SUM,mpc_fr_add);
Func2_MPFR_MPC(DIFF,mpc_fr_sub);
#define mpc_fr_mul(a,b,c,d) mpc_mul_fr(a,c,b,d)
Func2_MPFR_MPC(PROD,mpc_fr_mul);
Func2_MPFR_MPC(QUO,mpc_fr_div);

static Obj POW_MPFR_MPC(Obj self, Obj fl, Obj fr)
{
  mp_prec_t precl = mpfr_get_prec(GET_MPFR(fl)),
    precr = mpc_get_prec(GET_MPC(fr));

  Obj h = NEW_MPC (precl);
  mpfr_set (MPC_OBJ(h)->re, GET_MPFR(fl), GMP_RNDN);
  mpfr_set_ui (MPC_OBJ(h)->im, 0, GMP_RNDN);

  Obj g = NEW_MPC(precl > precr ? precl : precr);
  mpc_pow (MPC_OBJ(g), GET_MPC(h), GET_MPC(fr), MPC_RNDNN);
  return g;
}

static Obj LQUO_MPC(Obj self, Obj fl, Obj fr)
{
  return QUO_MPC(self, fr, fl);
}

static Obj LQUO_MPC_MPFR(Obj self, Obj fl, Obj fr)
{
  return QUO_MPFR_MPC(self, fr, fl);
}

static Obj LQUO_MPFR_MPC(Obj self, Obj fl, Obj fr)
{
  return QUO_MPC_MPFR(self, fr, fl);
}

static Obj EQ_MPC (Obj self, Obj fl, Obj fr)
{
  return mpc_cmp(GET_MPC(fl),GET_MPC(fr)) == 0 ? True : False;
}

static Obj EQ_MPC_MPFR (Obj self, Obj fl, Obj fr)
{
  return mpfr_cmp(GET_MPC(fl)->re,GET_MPFR(fr)) == 0 && mpfr_zero_p(MPC_OBJ(fl)->im) ? True : False;
}

static Obj EQ_MPFR_MPC (Obj self, Obj fl, Obj fr)
{
  return mpfr_cmp(GET_MPFR(fl),GET_MPC(fr)->re) == 0 && mpfr_zero_p(MPC_OBJ(fr)->im) ? True : False;
}

static Obj LT_MPC (Obj self, Obj fl, Obj fr)
{
  int c = mpc_cmp(GET_MPC(fl),GET_MPC(fr)), cre = MPC_INEX_RE(c);
  return (cre < 0 || (cre == 0 && MPC_INEX_IM(c) < 0)) ? True : False;
}

static Obj LT_MPC_MPFR (Obj self, Obj fl, Obj fr)
{
  int c = mpfr_cmp(GET_MPC(fl)->re, GET_MPFR(fr));
  if (!c) c = mpfr_cmp_si(GET_MPC(fl)->im, 0);
  return c < 0 ? True : False;
}

static Obj LT_MPFR_MPC (Obj self, Obj fl, Obj fr)
{
  int c = mpfr_cmp(GET_MPFR(fl), GET_MPC(fr)->re);
  if (!c) c = -mpfr_cmp_si(GET_MPC(fr)->im, 0);
  return c < 0 ? True : False;
}

static Obj MPC_2MPFR (Obj self, Obj fl, Obj fr)
{
  mp_prec_t precl = mpfr_get_prec(GET_MPFR(fl)),
    precr = mpfr_get_prec(GET_MPFR(fr));

  Obj g = NEW_MPC(precl > precr ? precl : precr);
  mpfr_set (MPC_OBJ(g)->re, GET_MPFR(fl), GMP_RNDN);
  mpfr_set (MPC_OBJ(g)->im, GET_MPFR(fr), GMP_RNDN);

  return g;
}

static Obj ROOTPOLY_MPC (Obj self, Obj coeffs, Obj precision)
{
  Obj result;
  Int i, numroots, degree = LEN_PLIST(coeffs)-1;
  mpc_t op[degree+1], zero[degree];
  mp_prec_t prec = INT_INTOBJ(precision);

  if (degree < 1) {
    result = NEW_PLIST(T_PLIST,0);
    SET_LEN_PLIST(result,0);
    return result;
  }

  for (i = 0; i <= degree; i++) {
    mpc_init2(op[degree-i],mpc_get_prec(GET_MPC(ELM_PLIST(coeffs,i+1))));
    mpc_set(op[degree-i],GET_MPC(ELM_PLIST(coeffs,i+1)),MPC_RNDNN);
    if (!mpc_number_p(op[degree-i]))
      return Fail;
  }

  for (i = 0; i < degree; i++)
    mpc_init2(zero[i],prec);

  /* !!! first call routine on low-precision IEEE numbers; if got n roots, refine them using accelerated Newton's method (x - ff'/(f'f' - ff")), or (x-jf/f' with j such that f(new pt) in minimal. if fails, call the high-precision method */

  numroots = cpoly_MPC (degree, op, zero, (int)prec);

  for (i = 0; i <= degree; i++)
    mpc_clear(op[degree-i]);

  if (numroots == -1)
    result = Fail;
  else {
    result = NEW_PLIST(T_PLIST,numroots);
    SET_LEN_PLIST(result,numroots);
    for (i = 1; i <= numroots; i++) {
      Obj t = NEW_MPC(mpc_get_prec(zero[i-1]));
      mpc_set (MPC_OBJ(t), zero[i-1], MPC_RNDNN);
      SET_ELM_PLIST(result,i, t);
    }
  }

  for (i = 0; i < degree; i++)
    mpc_clear(zero[i]);

  return result;
}

/****************************************************************
* export functions
****************************************************************/
static StructGVarFunc GVarFuncs [] = {
  Inc1_MPC(AINV_MPC),
  Inc1_MPC(ABS_MPC),
  Inc1_MPC(INV_MPC),
  Inc1_MPC(EXP_MPC),
  Inc1_MPC(LOG_MPC),
  Inc1_MPC(SQRT_MPC),
  Inc1_MPC(SQR_MPC),
  Inc1_MPC(CONJ_MPC),
  Inc1_MPC(PROJ_MPC),
  Inc1_MPC(SIN_MPC),
  Inc1_MPC(COS_MPC),
  Inc1_MPC(TAN_MPC),
  Inc1_MPC(ASIN_MPC),
  Inc1_MPC(ACOS_MPC),
  Inc1_MPC(ATAN_MPC),
  Inc1_MPC(SINH_MPC),
  Inc1_MPC(COSH_MPC),
  Inc1_MPC(TANH_MPC),
  Inc1_MPC(ASINH_MPC),
  Inc1_MPC(ACOSH_MPC),
  Inc1_MPC(ATANH_MPC),

  Inc1_MPC(ZERO_MPC),
  Inc1_MPC(ONE_MPC),
  Inc1_MPC(ISZERO_MPC),
  Inc1_MPC(ISNUMBER_MPC),
  Inc1_MPC(ISNAN_MPC),
  Inc1_MPC(ISINF_MPC),
  Inc1_MPC_arg(MPC_MAKENAN,"int"),
  Inc1_MPC_arg(MPC_MAKEINFINITY,"int"),
  Inc1_MPC_arg(MPC_INT,"int"),
  Inc1_MPC(PREC_MPC),
  Inc1_MPC(REAL_MPC),
  Inc1_MPC(IMAG_MPC),
  Inc1_MPC(NORM_MPC),
  Inc1_MPC(ARG_MPC),
  Inc1_MPC(MPC_MPFR),

  Inc1_MPC(FREXP_MPC),
  Inc2_MPC_arg(LDEXP_MPC,"complex, int"),
  Inc1_MPC(EXTREPOFOBJ_MPC),
  Inc1_MPC_arg(OBJBYEXTREP_MPC,"list"),

  Inc2_MPC(SUM),
  Inc2_MPC(DIFF),
  Inc2_MPC(PROD),
  Inc2_MPC(QUO),
  Inc2_MPC(LQUO),
  Inc2_MPC(POW),
  Inc2_MPC(EQ),
  Inc2_MPC(LT),
  Inc2_MPC_arg(MPC_STRING,"string, int"),
  Inc2_MPC_arg(STRING_MPC,"complex, int"),
  Inc2_MPC_arg(VIEWSTRING_MPC,"complex, int"),
  Inc2_MPC_arg(MPC_INTPREC,"int, int"),
  Inc2_MPC_arg(MPC_MPCPREC,"complex, int"),
  Inc2_MPC_arg(MPC_2MPFR,"complex, complex"),

  Inc2_MPC_arg(ROOTPOLY_MPC,"coeffs, prec"),

  {0}
};

int InitMPCKernel (void)
{
  InitHdlrFuncsFromTable (GVarFuncs);

  ImportGVarFromLibrary ("TYPE_MPC", &TYPE_MPC);
  return 0;
}

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

/****************************************************************************
**
*E  mpc.c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . .ends here
*/

quality91%

¤ Dauer der Verarbeitung: 0.2 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.