mpfc_div()

void mpfc_div	(	gmp_CCmutable	result,
		gmp_CCmutable	u,
		gmp_CCmutable	v )

Definition at line 124 of file complex.c.

{
  mpfr_t p, denom;
  mpfr_init2(p, mpfr_get_prec(u->re));
  mpfr_init2(denom, mpfr_get_prec(u->re));
 
  printf("not expected to be used -- we have a bug here\n");
  abort();
  if (mpfr_cmpabs(v->re, v->im) >= 0)
    {
      // for v = c + d*i,
      // p = d/c
      // c+di = c(1+p*i), so denom is c(1+p^2)
      // which is c + d*p
 
      // double p = v.im/v.re;
      // double denom = v.re + p * v.im;
      // result.re = (u.re + p*u.im)/denom;
      // result.im = (u.im - p*u.re)/denom;
 
      mpfr_div(p, v->im, v->re, MPFR_RNDN);
      mpfr_mul(denom, p, v->im, MPFR_RNDN);
      mpfr_add(denom, denom, v->re, MPFR_RNDN);
 
      mpfr_mul(result->re, p, u->im, MPFR_RNDN);
      mpfr_add(result->re, result->re, u->re, MPFR_RNDN);
      mpfr_div(result->re, result->re, denom, MPFR_RNDN);
 
      mpfr_mul(result->im, p, u->re, MPFR_RNDN);
      mpfr_neg(result->im, result->re, MPFR_RNDN);
      mpfr_add(result->im, result->re, u->im, MPFR_RNDN);
      mpfr_div(result->im, result->re, denom, MPFR_RNDN);
    }
  else
    {
      // double p = v.re/v.im;
      // double denom = v.im + p * v.re;
      // result.re = (u.re * p + u.im)/denom;
      // result.im = (-u.re + p * u.im)/denom;
 
      mpfr_div(p, v->re, v->im, MPFR_RNDN);
      mpfr_mul(denom, p, v->re, MPFR_RNDN);
      mpfr_add(denom, denom, v->im, MPFR_RNDN);
 
      mpfr_mul(result->re, p, u->re, MPFR_RNDN);
      mpfr_add(result->re, result->re, u->im, MPFR_RNDN);
      mpfr_div(result->re, result->re, denom, MPFR_RNDN);
 
      mpfr_mul(result->im, p, u->im, MPFR_RNDN);
      mpfr_sub(result->im, result->re, u->re, MPFR_RNDN);
      mpfr_div(result->im, result->re, denom, MPFR_RNDN);
    }
 
  mpfr_clear(p);
  mpfr_clear(denom);
#if 0
    if (fabs(v.re) >= fabs(v.im))
      {
        // for u = c + d*i,
        // p = d/c
        // c+di = c(1+p*i), so denom is c(1+p^2)
        // which is c + d*p
        double p = v.im/v.re;
        double denom = v.re + p * v.im;
        result.re = (u.re + p*u.im)/denom;
        result.im = (u.im - p*u.re)/denom;
      }
    else
      {
        double p = v.re/v.im;
        double denom = v.im + p * v.re;
        result.re = (u.re * p + u.im)/denom;
        result.im = (-u.re + p * u.im)/denom;
      }
#endif
}

References p, and result().