/* sqrt and mpz_sqrt+mpz_get_d don't round the same

   */


#include <float.h>


#ifndef DBL_MANT_DIG
#define DBL_MANT_DIG  53  /* by default assume IEEE double */
#endif


SCM scm_sqrt (SCM x);
SCM_DEFINE (scm_sqrt, "sqrt", 1, 0, 0,
            (SCM x),
	    "Return the square root of @var{x}.\n"
	    "\n"
	    "If @var{x} is exact and a perfect square then the return is\n"
	    "exact, otherwise the return is inexact (either real or\n"
	    "complex).")
#define FUNC_NAME s_scm_sqrt
{
  double root_d;
  SCM    root;
  int    neg;

  if (SCM_INUMP (x))
    {
      long x_l = SCM_INUM (x);
      neg = (x_l < 0);
      x_l = (neg ? -x_l : x_l);

      /* inums which fit the mantissa of a double are handled with sqrt().
         On 32-bit systems inums always fit, but on 64-bit systems we need
         to check.  */
#if SCM_I_FIXNUM_BIT > DBL_MANT_DIG
      if ((unsigned long) x_l < (1UL << DBL_MANT_DIG))
#endif
        {
          root_d = sqrt ((double) x_l);
          if (! neg && floor (root_d) == root_d)
            {
              /* perfect square, return as inum */
              long root_l = (long) root_d;
              return SCM_MAKINUM (neg ? -root_l : root_l);
            }
          goto inexact;
        }

      /* inums which don't fit the mantissa of a double are handled with the
         bignum code, to ensure full accuracy in the result  */
      root = scm_i_mkbig ();
      mpz_set_si (SCM_I_BIG_MPZ (root), x_l);
      x = root;
      goto big;
    }
  else if (SCM_BIGP (x))
    {
      int  scale;

      root = scm_i_mkbig ();
    big:
      if (mpz_perfect_square_p (SCM_I_BIG_MPZ (x)))
        {
          /* Perfect square, return as exact.
             Enhancement: mpz_perfect_square_p already calculates the root
             when proving squareness, it'd be nice to have some sort of
             "square root if perfect square" function.  */
          mpz_sqrt (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (x));
          scm_remember_upto_here_1 (x);
          return scm_i_normbig (root);
        }

      neg = (mpz_sgn (SCM_I_BIG_MPZ (x)) < 0);
      mpz_abs (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (x));
      scm_remember_upto_here_1 (x);

      /* We need DBL_MANT_DIG bits of root, which means twice that of
         operand.  If we've got less, then pad with zeros, effectively
         generating fractional bits for the root.  If we've got more, then
         truncate so as to save work in mpz_sqrt.  */
      scale = ((long) mpz_sizeinbase (SCM_I_BIG_MPZ (root), 2)
               - 1 - 2*DBL_MANT_DIG) / 2;
      if (scale > 0)
        mpz_tdiv_q_2exp (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (root),
                         2 * scale);
      else
        mpz_mul_2exp (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (root),
                      - 2 * scale);

      mpz_sqrt (SCM_I_BIG_MPZ (root), SCM_I_BIG_MPZ (root));
      root_d = ldexp (mpz_get_d (SCM_I_BIG_MPZ (root)), scale);
      scm_remember_upto_here_1 (root);
      goto inexact;
    }
  else if (SCM_REALP (x))
    {
      double x_d = SCM_REAL_VALUE (x);
      neg = (x_d < 0.0);
      root_d = sqrt (fabs (x_d));

    inexact:
      if (neg)
        return scm_make_complex (0.0, root_d);
      else
        return scm_make_real (root_d);
    }
  else if (SCM_COMPLEXP (x))
    {
      /* sqrt(x) = (make-polar (sqrt (magnitude x)) (/ (angle x) 2))
         ENHANCE-ME: Use csqrt() when available.  */
      double re = SCM_COMPLEX_REAL (x);
      double im = SCM_COMPLEX_IMAG (x);
      double mag = sqrt (hypot (re, im));
      double ang = atan2 (re, im) / 2.0;
      double s, c;
#if HAVE_SINCOS
      sincos (ang, &s, &c);
#else
      s = sin (ang);
      c = cos (ang);
#endif
      return scm_make_complex (mag * c, mag * s);
    }
  else
    SCM_WRONG_TYPE_ARG (SCM_ARG1, x);
}
#undef FUNC_NAME