mpfc_sqrt()

void mpfc_sqrt	(	gmp_CCmutable	result,
		gmp_CC	a )

Definition at line 282 of file complex.c.

{
  // The idea is: write a = a1 + i * a2
  // first make it numerically more stable by dividing by the larger
  //   answer will be multiplied by larger afterwards, if needed
  // To take the square root of 1+di, or -1+di
  //   it is enough to solve (if sqrt is e+fi), e^2+f^2 = sqrt(1+d^2),
  //   and e^2-f^2 = 1 (or -1).
 
  if (mpfr_zero_p(a->re) && mpfr_zero_p(a->im))
    {
      mpfr_set_si(result->re, 0, MPFR_RNDN);
      mpfr_set_si(result->im, 0, MPFR_RNDN);
      return;
    }
 
  mpfr_t b, c, d, d2;
 
  unsigned long p = mpfr_get_prec(a->re);
  mpfr_init2(b, p);
  mpfr_init2(c, p);
  mpfr_init2(d, p);
  mpfr_init2(d2, p);
 
  // b = fabs(a.re);
  // c = fabs(a.im);
  mpfr_abs(b, a->re, MPFR_RNDN);
  mpfr_abs(c, a->im, MPFR_RNDN);
  if (mpfr_greater_p(b, c))
    {
      // d = c/b;
      // e = ::sqrt(b) * ::sqrt(0.5 * (1.0 + ::sqrt(1.0 + d*d)));
      // but use c as d, and also as e
      mpfr_div(d, c, b, MPFR_RNDN);
      mpfr_sqr(d, d, MPFR_RNDN);
      mpfr_add_si(d, d, 1, MPFR_RNDN);
      mpfr_sqrt(d, d, MPFR_RNDN);
      mpfr_add_si(d, d, 1, MPFR_RNDN);
      mpfr_div_2ui(d, d, 1, MPFR_RNDN);
      mpfr_sqrt(d, d, MPFR_RNDN);
      mpfr_sqrt(b, b, MPFR_RNDN);
      mpfr_mul(d, d, b, MPFR_RNDN); /* this is now e ! */
    }
  else
    {
      // d = b/c;
      // e = ::sqrt(c) * ::sqrt(0.5 * (d + ::sqrt(1.0 + d*d)));
      mpfr_div(d, b, c, MPFR_RNDN);
      mpfr_sqr(d2, d, MPFR_RNDN);
      mpfr_add_si(d2, d2, 1, MPFR_RNDN);
      mpfr_sqrt(d2, d2, MPFR_RNDN);
      mpfr_add(d, d2, d, MPFR_RNDN);
      mpfr_div_2ui(d, d, 1, MPFR_RNDN);
      mpfr_sqrt(d, d, MPFR_RNDN);
      mpfr_sqrt(c, c, MPFR_RNDN);
      mpfr_mul(d, d, c, MPFR_RNDN); /* this is now e ! */
    }
  if (mpfr_sgn(a->re) >= 0)
    {
      // result.re = e;
      // result.im = a.im/(2.0 * e);
      mpfr_set(result->re, d, MPFR_RNDN);
      mpfr_mul_2ui(d, d, 1, MPFR_RNDN);
      mpfr_div(result->im, a->im, d, MPFR_RNDN);
    }
  else
    {
      // result.im = (a.im >= 0.0 ? e : -e);
      // result.re = a.re/(2.0*result.im);
      mpfr_set(result->im, d, MPFR_RNDN);
      if (mpfr_sgn(a->im) < 0) mpfr_neg(result->im, result->im, MPFR_RNDN);
      mpfr_mul_2ui(d, result->im, 1, MPFR_RNDN);
      mpfr_div(result->re, a->im, d, MPFR_RNDN);
    }
 
  mpfr_clear(b);
  mpfr_clear(c);
  mpfr_clear(d);
  mpfr_clear(d2);
 
#if 0
  static void sqrt(elem &result, elem a)
  {
    double b,c,d,e;
    if (a.re == 0.0 && a.im == 0.0)
      {
        result.re = 0.0;
        result.im = 0.0;
      }
    else
      {
        b = fabs(a.re);
        c = fabs(a.im);
        if (b > c)
          {
            d = c/b;
            e = ::sqrt(b) * ::sqrt(0.5 * (1.0 + ::sqrt(1.0 + d*d)));
          }
        else
          {
            d = b/c;
            e = ::sqrt(c) * ::sqrt(0.5 * (d + ::sqrt(1.0 + d*d)));
          }
        // Now be careful with the signs:
        if (a.re >= 0.0)
          {
            result.re = e;
            result.im = a.im/(2.0 * e);
          }
        else
          {
            result.im = (a.im >= 0.0 ? e : -e);
            result.re = a.re/(2.0*result.im);
          }
      }
  }
#endif
}

References p, and result().