/* See aaiafuncs.h */
/* Last edited on 2007-01-04 00:12:32 by stolfi */

#include <aa.h>
#include <ia.h>
#include <flt.h>

#include <aaiafuncs.h>

#include <math.h>

/* AA/IA FUNCTION WITH LARGE OVERESTIMATION */

Float g_fp (Float x)
  {
    ROUND_NEAR;
    return (sqrt(x*x - x + 0.50)/sqrt(x*x + 0.50));
  }

Interval g_ia (Interval x)
  {
    Interval x2  = ia_sqr(x);
    Interval hlf = {Half, Half};
    Interval dif = ia_sub(x2, x);
    Interval sum1 = ia_add(dif, hlf);
    Interval sum2 = ia_add(x2, hlf);
    Interval res = ia_div(ia_sqrt(sum1), ia_sqrt(sum2));
    return (res);
  }

AAP g_aa (AAP x)
  {
    AAP x2  = aa_sqr(x);
    AAP hlf = aa_const(Half, Zero);
    AAP dif = aa_sub(x2, x);
    AAP sum1 = aa_add(dif, hlf);
    AAP sum2 = aa_add(x2, hlf);
    AAP res = aa_div(aa_sqrt(sum1), aa_sqrt(sum2));
    return (res);
  }
  
/* ITERATED AA/IA FUNCTION WITH ERROR EXPLOSION */

Float gg_fp (Float x)
  {
    ROUND_NEAR;
    return g_fp(g_fp(x));
  }

Interval gg_ia (Interval x)
  {
    return (g_ia(g_ia(x)));
  }

AAP gg_aa (AAP x)
  {
    MemP frame = aa_top();
    AAP res = g_aa(g_aa(x));
    return (aa_return(frame, res));
  }

/* AA/IA SQUARE ROOT FUNCTION */

Float sqrt_fp (Float x)
  {
    ROUND_NEAR;
    if (x >= 0) return (sqrt(x));
    else return (Zero);
  }

Interval sqrt_ia (Interval x)
  {
    if (x.hi >= 0) return (ia_sqrt(x));
    else return (ia_const(Zero, Zero));
  }

AAP sqrt_aa (AAP x)
  {
    Interval r = aa_range(x);
    if (r.hi >= 0) return (aa_sqrt(x));
    else return (aa_zero());
  }