自学教程：C++ IMAGPART函数代码示例

intelem_equal(elem_srcptr op1, elem_srcptr op2, const ring_t ring){    switch (ring->type)    {        case TYPE_FMPZ:            return fmpz_equal(op1, op2);        case TYPE_LIMB:            return *((mp_srcptr) op1) == *((mp_srcptr) op2);        case TYPE_POLY:            return elem_poly_equal(op1, op2, ring);        case TYPE_MOD:            return elem_equal(op1, op2, ring->parent);        case TYPE_FRAC:            return elem_equal(NUMER(op1, ring), NUMER(op2, ring), RING_NUMER(ring)) &&                    elem_equal(DENOM(op1, ring), DENOM(op2, ring), RING_DENOM(ring));        case TYPE_COMPLEX:            return elem_equal(REALPART(op1, ring), REALPART(op2, ring), RING_PARENT(ring)) &&                    elem_equal(IMAGPART(op1, ring), IMAGPART(op2, ring), RING_PARENT(ring));        default:            NOT_IMPLEMENTED("equal", ring);    }}

static inline __complex128 mult_c128 (__complex128 x, __complex128 y){  __float128 r1 = REALPART(x), i1 = IMAGPART(x);  __float128 r2 = REALPART(y), i2 = IMAGPART(y);  __complex128 res;  COMPLEX_ASSIGN(res, r1*r2 - i1*i2, i2*r1 + i1*r2);  return res;}

// Careful: the algorithm for the division sucks. A lot.static inline __complex128 div_c128 (__complex128 x, __complex128 y){  __float128 n = hypotq (REALPART (y), IMAGPART (y));  __float128 r1 = REALPART(x), i1 = IMAGPART(x);  __float128 r2 = REALPART(y), i2 = IMAGPART(y);  __complex128 res;  COMPLEX_ASSIGN(res, r1*r2 + i1*i2, i1*r2 - i2*r1);  return res / n;}

voidelem_sub(elem_ptr res, elem_srcptr op1, elem_srcptr op2, const ring_t ring){    switch (ring->type)    {        case TYPE_FMPZ:            fmpz_sub(res, op1, op2);            break;        case TYPE_LIMB:            *((mp_ptr) res) = *((mp_srcptr) op1) - *((mp_srcptr) op2);            break;        case TYPE_POLY:            elem_poly_sub(res, op1, op2, ring);            break;        case TYPE_MOD:            {                switch (RING_PARENT(ring)->type)                {                    case TYPE_LIMB:                        *((mp_ptr) res) = n_submod(*((mp_srcptr) op1), *((mp_srcptr) op2), ring->nmod.n);                        break;                    case TYPE_FMPZ:                        fmpz_sub(res, op1, op2);                        if (fmpz_sgn(res) < 0)                            fmpz_add(res, res, RING_MODULUS(ring));                        break;                    default:                        NOT_IMPLEMENTED("sub (mod)", ring);                }            }            break;        case TYPE_FRAC:            elem_frac_sub(res, op1, op2, ring);            break;        case TYPE_COMPLEX:            elem_sub(REALPART(res, ring), REALPART(op1, ring), REALPART(op2, ring), ring->parent);            elem_sub(IMAGPART(res, ring), IMAGPART(op1, ring), IMAGPART(op2, ring), ring->parent);            break;        default:            NOT_IMPLEMENTED("sub", ring);    }}

doublecimag(double complex z){	const double_complex z1 = { .f = z };	return (IMAGPART(z1));}

floatcimagf(float complex z){	const float_complex z1 = { .f = z };	return (IMAGPART(z1));}

/* Complex argument.  The angle made with the +ve real axis.   Range -pi-pi.  */GFC_REAL_4cargf (GFC_COMPLEX_4 z){  GFC_REAL_4 arg;  return atan2f (IMAGPART (z), REALPART (z));}

long doublecimagl(long double complex z){	const long_double_complex z1 = { .f = z };	return (IMAGPART(z1));}

__complex128ccoshq (__complex128 a){  __float128 r = REALPART (a), i = IMAGPART (a);  __complex128 v;  COMPLEX_ASSIGN (v, coshq (r) * cosq (i),  sinhq (r) * sinq (i));  return v;}

__complex128ctanhq (__complex128 a){  __float128 rt = tanhq (REALPART (a)), it = tanq (IMAGPART (a));  __complex128 n, d;  COMPLEX_ASSIGN (n, rt, it);  COMPLEX_ASSIGN (d, 1, rt * it);  return C128_DIV(n,d);}

long double complexcexpl (long double complex z){  long double a, b;  long double complex v;  a = REALPART (z);  b = IMAGPART (z);  COMPLEX_ASSIGN (v, cosl (b), sinl (b));  return expl (a) * v;}

double complexcexp (double complex z){  double a, b;  double complex v;  a = REALPART (z);  b = IMAGPART (z);  COMPLEX_ASSIGN (v, cos (b), sin (b));  return exp (a) * v;}

float complexcexpf (float complex z){  float a, b;  float complex v;  a = REALPART (z);  b = IMAGPART (z);  COMPLEX_ASSIGN (v, cosf (b), sinf (b));  return expf (a) * v;}

float complexccoshf (float complex a){  float r, i;  float complex v;  r = REALPART (a);  i = IMAGPART (a);  COMPLEX_ASSIGN (v, coshf (r) * cosf (i), - (sinhf (r) * sinf (i)));  return v;}

__complex128cexpq (__complex128 z){  __float128 a, b;  __complex128 v;  a = REALPART (z);  b = IMAGPART (z);  COMPLEX_ASSIGN (v, cosq (b), sinq (b));  return expq (a) * v;}

long double complexccoshl (long double complex a){  long double r, i;  long double complex v;  r = REALPART (a);  i = IMAGPART (a);  COMPLEX_ASSIGN (v, coshl (r) * cosl (i), - (sinhl (r) * sinl (i)));  return v;}

double complexccosh (double complex a){  double r, i;  double complex v;  r = REALPART (a);  i = IMAGPART (a);  COMPLEX_ASSIGN (v, cosh (r) * cos (i), - (sinh (r) * sin (i)));  return v;}

/* cosh(z) = cosh(a)cos(b) - isinh(a)sin(b)  */GFC_COMPLEX_4ccoshf (GFC_COMPLEX_4 a){  GFC_REAL_4 r;  GFC_REAL_4 i;  GFC_COMPLEX_4 v;  r = REALPART (a);  i = IMAGPART (a);  COMPLEX_ASSIGN (v, coshf (r) * cosf (i), - (sinhf (r) * sinf (i)));  return v;}

/* exp(z) = exp(a)*(cos(b) + isin(b))  */GFC_COMPLEX_4cexpf (GFC_COMPLEX_4 z){  GFC_REAL_4 a;  GFC_REAL_4 b;  GFC_COMPLEX_4 v;  a = REALPART (z);  b = IMAGPART (z);  COMPLEX_ASSIGN (v, cosf (b), sinf (b));  return expf (a) * v;}

float complexctanf (float complex a){  float rt, it;  float complex n, d;  rt = tanf (REALPART (a));  it = tanhf (IMAGPART (a));  COMPLEX_ASSIGN (n, rt, it);  COMPLEX_ASSIGN (d, 1, - (rt * it));  return n / d;}

double complexctan (double complex a){  double rt, it;  double complex n, d;  rt = tan (REALPART (a));  it = tanh (IMAGPART (a));  COMPLEX_ASSIGN (n, rt, it);  COMPLEX_ASSIGN (d, 1, - (rt * it));  return n / d;}

long double complexctanl (long double complex a){  long double rt, it;  long double complex n, d;  rt = tanl (REALPART (a));  it = tanhl (IMAGPART (a));  COMPLEX_ASSIGN (n, rt, it);  COMPLEX_ASSIGN (d, 1, - (rt * it));  return n / d;}

/* sqrt(z).  Algorithm pulled from glibc.  */GFC_COMPLEX_4csqrtf (GFC_COMPLEX_4 z){  GFC_REAL_4 re;  GFC_REAL_4 im;  GFC_COMPLEX_4 v;  re = REALPART (z);  im = IMAGPART (z);  if (im == 0.0)    {      if (re < 0.0)        {          COMPLEX_ASSIGN (v, 0.0, copysignf (sqrtf (-re), im));        }      else        {          COMPLEX_ASSIGN (v, fabsf (sqrt (re)),                          copysignf (0.0, im));        }    }  else if (re == 0.0)    {      GFC_REAL_4 r;      r = sqrtf (0.5 * fabs (im));      COMPLEX_ASSIGN (v, copysignf (r, im), r);    }  else    {      GFC_REAL_4 d, r, s;      d = hypotf (re, im);      /* Use the identity   2  Re res  Im res = Im x         to avoid cancellation error in  d +/- Re x.  */      if (re > 0)        {          r = sqrtf (0.5 * d + 0.5 * re);          s = (0.5 * im) / r;        }      else        {          s = sqrtf (0.5 * d - 0.5 * re);          r = fabsf ((0.5 * im) / s);        }      COMPLEX_ASSIGN (v, r, copysignf (s, im));    }  return v;}

voidelem_set(elem_ptr res, elem_srcptr src, const ring_t ring){    if (res != src)    {        switch (ring->type)        {            case TYPE_FMPZ:                fmpz_set(res, src);                break;            case TYPE_LIMB:                *((mp_ptr) res) = *((mp_srcptr) src);                break;            case TYPE_MOD:                elem_set(res, src, ring->parent);                break;            case TYPE_POLY:                elem_poly_set(res, src, ring);                break;            case TYPE_FRAC:                elem_set(NUMER(res, ring), NUMER(src, ring), RING_NUMER(ring));                elem_set(DENOM(res, ring), DENOM(src, ring), RING_DENOM(ring));                break;            case TYPE_COMPLEX:                elem_set(REALPART(res, ring), REALPART(src, ring), RING_PARENT(ring));                elem_set(IMAGPART(res, ring), IMAGPART(src, ring), RING_PARENT(ring));                break;            default:                NOT_IMPLEMENTED("set", ring);        }    }}

voidelem_set_si(elem_ptr elem, long v, const ring_t ring){    switch (ring->type)    {        case TYPE_FMPZ:            fmpz_set_si(elem, v);            break;        case TYPE_LIMB:            *((mp_ptr) elem) = v;            break;        case TYPE_POLY:            elem_poly_set_si(elem, v, ring);            break;        case TYPE_MOD:            {                switch (RING_PARENT(ring)->type)                {                    case TYPE_FMPZ:                        fmpz_set_si(elem, v);                        fmpz_mod(elem, elem, RING_MODULUS(ring));                        break;                    case TYPE_LIMB:                        *((mp_ptr) elem) = nmod_set_si(v, ring->nmod);                        break;                    default:                        NOT_IMPLEMENTED("set_si (mod)", ring);                }            }            break;        case TYPE_FRAC:            elem_set_si(NUMER(elem, ring), v, RING_NUMER(ring));            elem_one(DENOM(elem, ring), RING_DENOM(ring));            break;        case TYPE_COMPLEX:            elem_set_si(REALPART(elem, ring), v, ring->parent);            elem_zero(IMAGPART(elem, ring), ring->parent);            break;        default:            NOT_IMPLEMENTED("set_si", ring);    }}

float complexcsqrtf (float complex z){  float re, im;  float complex v;  re = REALPART (z);  im = IMAGPART (z);  if (im == 0)    {      if (re < 0)        {          COMPLEX_ASSIGN (v, 0, copysignf (sqrtf (-re), im));        }      else        {          COMPLEX_ASSIGN (v, fabsf (sqrtf (re)), copysignf (0, im));        }    }  else if (re == 0)    {      float r;      r = sqrtf (0.5 * fabsf (im));      COMPLEX_ASSIGN (v, r, copysignf (r, im));    }  else    {      float d, r, s;      d = hypotf (re, im);      /* Use the identity   2  Re res  Im res = Im x         to avoid cancellation error in  d +/- Re x.  */      if (re > 0)        {          r = sqrtf (0.5 * d + 0.5 * re);          s = (0.5 * im) / r;        }      else        {          s = sqrtf (0.5 * d - 0.5 * re);          r = fabsf ((0.5 * im) / s);        }      COMPLEX_ASSIGN (v, r, copysignf (s, im));    }  return v;}

long double complexcsqrtl (long double complex z){  long double re, im;  long double complex v;  re = REALPART (z);  im = IMAGPART (z);  if (im == 0)    {      if (re < 0)        {          COMPLEX_ASSIGN (v, 0, copysignl (sqrtl (-re), im));        }      else        {          COMPLEX_ASSIGN (v, fabsl (sqrtl (re)), copysignl (0, im));        }    }  else if (re == 0)    {      long double r;      r = sqrtl (0.5 * fabsl (im));      COMPLEX_ASSIGN (v, copysignl (r, im), r);    }  else    {      long double d, r, s;      d = hypotl (re, im);      /* Use the identity   2  Re res  Im res = Im x         to avoid cancellation error in  d +/- Re x.  */      if (re > 0)        {          r = sqrtl (0.5 * d + 0.5 * re);          s = (0.5 * im) / r;        }      else        {          s = sqrtl (0.5 * d - 0.5 * re);          r = fabsl ((0.5 * im) / s);        }      COMPLEX_ASSIGN (v, r, copysignl (s, im));    }  return v;}

/* tanh(z) = (tanh(a) + itan(b)) / (1 - itanh(a)tan(b))  */GFC_COMPLEX_4ctanhf (GFC_COMPLEX_4 a){  GFC_REAL_4 rt;  GFC_REAL_4 it;  GFC_COMPLEX_4 n;  GFC_COMPLEX_4 d;  rt = tanhf (REALPART (a));  it = tanf (IMAGPART (a));  COMPLEX_ASSIGN (n, rt, it);  COMPLEX_ASSIGN (d, 1, - (rt * it));  return n / d;}

/* Square root algorithm from glibc.  */__complex128csqrtq (__complex128 z){  __float128 re = REALPART(z), im = IMAGPART(z);  __complex128 v;  if (im == 0)  {    if (re < 0)    {      COMPLEX_ASSIGN (v, 0, copysignq (sqrtq (-re), im));    }    else    {      COMPLEX_ASSIGN (v, fabsq (sqrtq (re)), copysignq (0, im));    }  }  else if (re == 0)  {    __float128 r = sqrtq (0.5 * fabsq (im));    COMPLEX_ASSIGN (v, r, copysignq (r, im));  }  else  {    __float128 d = hypotq (re, im);    __float128 r, s;    /* Use the identity   2  Re res  Im res = Im x	to avoid cancellation error in  d +/- Re x.  */    if (re > 0)      r = sqrtq (0.5 * d + 0.5 * re), s = (0.5 * im) / r;    else      s = sqrtq (0.5 * d - 0.5 * re), r = fabsq ((0.5 * im) / s);    COMPLEX_ASSIGN (v, r, copysignq (s, im));  }  return v;}

/* Both parameters will already have been converted to the result type.  */GFC_COMPLEX_4__dot_product_c4 (gfc_array_c4 * a, gfc_array_c4 * b){  GFC_COMPLEX_4 *pa;  GFC_COMPLEX_4 *pb;  GFC_COMPLEX_4 res;  GFC_COMPLEX_4 conjga;  index_type count;  index_type astride;  index_type bstride;  assert (GFC_DESCRIPTOR_RANK (a) == 1          && GFC_DESCRIPTOR_RANK (b) == 1);  if (a->dim[0].stride == 0)    a->dim[0].stride = 1;  if (b->dim[0].stride == 0)    b->dim[0].stride = 1;  astride = a->dim[0].stride;  bstride = b->dim[0].stride;  count = a->dim[0].ubound + 1 - a->dim[0].lbound;  res = 0;  pa = a->data;  pb = b->data;  while (count--)    {      COMPLEX_ASSIGN(conjga, REALPART (*pa), -IMAGPART (*pa));      res += conjga * *pb;      pa += astride;      pb += bstride;    }  return res;}