int
elem_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;
}

void
elem_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);
    }
}

double
cimag(double complex z)
{
	const double_complex z1 = { .f = z };

	return (IMAGPART(z1));
}

float
cimagf(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_4
cargf (GFC_COMPLEX_4 z)
{
  GFC_REAL_4 arg;

  return atan2f (IMAGPART (z), REALPART (z));
}

long double
cimagl(long double complex z)
{
	const long_double_complex z1 = { .f = z };

	return (IMAGPART(z1));
}

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

__complex128
ctanhq (__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 complex
cexpl (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 complex
cexp (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 complex
cexpf (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 complex
ccoshf (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;
}

__complex128
cexpq (__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 complex
ccoshl (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 complex
ccosh (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_4
ccoshf (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_4
cexpf (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 complex
ctanf (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;
}