mpfc_abs()

void mpfc_abs	(	gmp_RRmutable	result,
		gmp_CCmutable	c )

Definition at line 225 of file complex.c.

{
  mpfr_t a, b;
 
  mpfr_init2(a, mpfr_get_prec(c->re));
  mpfr_init2(b, mpfr_get_prec(c->re));
  mpfr_abs(a, c->re, MPFR_RNDN);
  mpfr_abs(b, c->im, MPFR_RNDN);
  if (mpfr_zero_p(a))
    mpfr_set(result, b, MPFR_RNDN);
  else if (mpfr_zero_p(b))
    mpfr_set(result, a, MPFR_RNDN);
  else if (mpfr_greater_p(a, b))
    {
      // double d = b/a;
      // But use b for d, as it is not needed later.
      // result = a * ::sqrt(1.0 + d*d);
      mpfr_div(b, b, a, MPFR_RNDN);
      mpfr_sqr(b, b, MPFR_RNDN);
      mpfr_add_si(b, b, 1, MPFR_RNDN);
      mpfr_sqrt(b, b, MPFR_RNDN);
      mpfr_mul(result, a, b, MPFR_RNDN);
    }
  else
    {
      // double d = a/b;
      // But use a for d, as it is not needed later.
      // result = b * ::sqrt(1.0 + d*d);
      mpfr_div(a, a, b, MPFR_RNDN);
      mpfr_sqr(a, a, MPFR_RNDN);
      mpfr_add_si(a, a, 1, MPFR_RNDN);
      mpfr_sqrt(a, a, MPFR_RNDN);
      mpfr_mul(result, b, a, MPFR_RNDN);
    }
  mpfr_clear(a);
  mpfr_clear(b);
 
#if 0
  static void abs(double &result, elem c)
  {
    double a = fabs(c.re);
    double b = fabs(c.im);
    if (a == 0.0) result = b;
    else if (b == 0.0) result = a;
    else if (a > b)
      {
        double d = b/a;
        result = a * ::sqrt(1.0 + d*d);
      }
    else
      {
        double d = a/b;
        result = b * ::sqrt(1.0 + d*d);
      }
  }
#endif
}

References abs, and result().