void test_fp_utilities( void )
{
#if __STDC_VERSION__ >= 199901L
    printf( "Testing C99 miscellaneous functions...\n" );

    VERIFY( CompDbl( copysign( -2.0, 1.0), 2.0 ) );
    VERIFY( CompDbl( copysign( -2.0, -1.0), -2.0 ) );
    VERIFY( CompDbl( copysign( 2.0, -1.0), -2.0 ) );
    VERIFY( CompDbl( copysign( 2.0, 1.0), 2.0 ) );
    
    VERIFY( CompDbl( fmax( 2.0, 1.0), 2.0 ) );
    VERIFY( CompDbl( fmax( -2.0, -1.0), -1.0 ) );
    VERIFY( CompDbl( fmin( 2.0, 1.0), 1.0 ) );
    VERIFY( CompDbl( fmin( -2.0, -1.0), -2.0 ) );
    
    VERIFY( CompDbl( fma( 2.0, 3.0, 4.0), 10.0 ) );
    VERIFY( CompDbl( fma( 2.0, 3.0, -4.0), 2.0 ) );
    VERIFY( CompDbl( fma( -2.0, 3.0, 4.0), -2.0 ) );
    VERIFY( CompDbl( fma( -2.0, -3.0, 4.0), 10.0 ) );
    
    VERIFY( CompDbl( fdim( 3.0, 2.0), 1.0 ) );
    VERIFY( CompDbl( fdim( 2.0, 3.0), 0.0 ) );
    
    VERIFY( CompDbl( nextafter( 1.0, 2.0), 1.0+1.0E-16 ) );
    VERIFY( CompDbl( nextafter( 1.0, 0.0), 1.0-1.0E-16 ) );
    
    VERIFY( CompDbl( scalbn( 1.0, 3.0), 8.0 ) );
    VERIFY( CompDbl( scalbn( 4.0, 3.0), 32.0 ) );
#endif
}

void test_scalbn()
{
    static_assert((std::is_same<decltype(scalbn((double)0, (int)0)), double>::value), "");
    static_assert((std::is_same<decltype(scalbnf(0, (int)0)), float>::value), "");
    static_assert((std::is_same<decltype(scalbnl(0, (int)0)), long double>::value), "");
    assert(scalbn(1, 1) == 2);
}

fcomplex
cexpf(fcomplex z) {
	fcomplex	ans;
	float		x, y, c, s;
	double		t;
	int		n, ix, iy, hx, hy;

	x = F_RE(z);
	y = F_IM(z);
	hx = THE_WORD(x);
	hy = THE_WORD(y);
	ix = hx & 0x7fffffff;
	iy = hy & 0x7fffffff;
	if (iy == 0) {		/* y = 0 */
		F_RE(ans) = expf(x);
		F_IM(ans) = y;
	} else if (ix == 0x7f800000) {	/* x is +-inf */
		if (hx < 0) {
			if (iy >= 0x7f800000) {
				F_RE(ans) = zero;
				F_IM(ans) = zero;
			} else {
				sincosf(y, &s, &c);
				F_RE(ans) = zero * c;
				F_IM(ans) = zero * s;
			}
		} else {
			if (iy >= 0x7f800000) {
				F_RE(ans) = x;
				F_IM(ans) = y - y;
			} else {
				sincosf(y, &s, &c);
				F_RE(ans) = x * c;
				F_IM(ans) = x * s;
			}
		}
	} else {
		sincosf(y, &s, &c);
		if (ix >= 0x42B171AA) {	/* |x| > 88.722... ~ log(2**128) */
#if defined(__i386) && !defined(__amd64)
			int	rp = __swapRP(fp_extended);
#endif
			t = __k_cexp(x, &n);
			F_RE(ans) = (float)scalbn(t * (double)c, n);
			F_IM(ans) = (float)scalbn(t * (double)s, n);
#if defined(__i386) && !defined(__amd64)
			if (rp != fp_extended)
				(void) __swapRP(rp);
#endif
		} else {
			t = expf(x);
			F_RE(ans) = t * c;
			F_IM(ans) = t * s;
		}
	}
	return (ans);
}

double
log1p(double x)
{
	static const double zero=0.0, negone= -1.0, one=1.0,
		      half=1.0/2.0, small=1.0E-20;   /* 1+small == 1 */
	double z,s,t,c;
	int k;

	if (isnan(x))
		return (x);

	if(finite(x)) {
	   if( x > negone ) {

	   /* argument reduction */
	      if(copysign(x,one)<small) return(x);
	      k=logb(one+x); z=scalbn(x,-k); t=scalbn(one,-k);
	      if(z+t >= sqrt2 )
		  { k += 1 ; z *= half; t *= half; }
	      t += negone; x = z + t;
	      c = (t-x)+z ;		/* correction term for x */

 	   /* compute log(1+x)  */
              s = x/(2+x); t = x*x*half;
	      c += (k*ln2lo-c*x);
	      z = c+s*(t+__log__L(s*s));
	      x += (z - t) ;

	      return(k*ln2hi+x);
	   }
	/* end of if (x > negone) */

	    else {
#if defined(__vax__)
		if ( x == negone )
		    return (infnan(-ERANGE));	/* -INF */
		else
		    return (infnan(EDOM));	/* NaN */
#else	/* defined(__vax__) */
		/* x = -1, return -INF with signal */
		if ( x == negone ) return( negone/zero );

		/* negative argument for log, return NaN with signal */
	        else return ( zero / zero );
#endif	/* defined(__vax__) */
	    }
	}
    /* end of if (finite(x)) */

    /* log(-INF) is NaN */
	else if(x<0)
	     return(zero/zero);

    /* log(+INF) is INF */
	else return(x);
}

inline T normalize_value(const T& val, const mpl::true_&)
{
   BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
   BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);

   boost::intmax_t shift = std::numeric_limits<T>::digits - ilogb(val) - 1;
   T result = scalbn(val, shift);
   result = round(result);
   return scalbn(result, -shift);
}

double attribute_hidden __ieee754_scalb(double x, double fn)
{
	if (isnan(x)||isnan(fn)) return x*fn;
	if (!isfinite(fn)) {
	    if(fn>0.0) return x*fn;
	    else       return x/(-fn);
	}
	if (rint(fn)!=fn) return (fn-fn)/(fn-fn);
	if ( fn > 65000.0) return scalbn(x, 65000);
	if (-fn > 65000.0) return scalbn(x,-65000);
	return scalbn(x,(int)fn);
}

double
scalb(double x, double fn)
{
	if (isnan(x)||isnan(fn)) return x*fn;
	if (!finite(fn)) {
	    if(fn>0.0) return x*fn;
	    else       return x/(-fn);
	}
	if (rint(fn)!=fn) return (fn-fn)/(fn-fn);
	if ( fn > 65000.0) return scalbn(x, 65000);
	if (-fn > 65000.0) return scalbn(x,-65000);
	return scalbn(x,(int)fn);
}

/*
 * Convert a floating point number in native format to TMS32032
 * floating point single precision (32 bit) format.  Note that the
 * TMS floating point value is returned as an unsigned integer.
 */
static unsigned int
float2tms32(float x)
{
    unsigned int zero = 0x80000000; /* Zero value is special case */
    int nfracbits = 23;    /* Not including hidden bit / sign bit */
    int signbit = 1 << nfracbits;
    int fracmask = ~((~0)<<nfracbits);
    int iexp;
    int sign;
    int ifrac;
    unsigned int rtn;

    if (x == 0){
        rtn = zero;
    }else{
        iexp = ilogb(x);  /* Binary exponent if 1 <= |fraction| < 2 */
        ifrac = (int)scalbn(x, nfracbits-iexp); /* Frac part as integer */
        if (x<0 && (ifrac & signbit)){
            /* Force top bit of negative fraction to be 0 */
            ifrac <<= 1;
            iexp--;
        }
        sign = x<0 ? signbit : 0;
        rtn = (iexp << (nfracbits+1)) | sign | (ifrac & fracmask);
    }
    return rtn;
}

double
cosh(double x) {
	double t, w;

	w = fabs(x);
	if (!finite(w))
		return (w * w);
	if (w < 0.3465) {
		t = expm1(w);
		w = 1.0 + t;
		if (w != 1.0)
			w = 1.0 + (t * t) / (w + w);
		return (w);
	} else if (w < 22.0) {
		t = exp(w);
		return (0.5 * (t + 1.0 / t));
	} else if (w <= lnovft) {
		return (0.5 * exp(w));
	} else {
		w = (w - 1024 * ln2hi) - 1024 * ln2lo;
		if (w >= ln2)
			return (_SVID_libm_err(x, x, 5));
		else
			return (scalbn(exp(w), 1023));
	}
}

T ulp_imp(const T& val, const mpl::false_&, const Policy& pol)
{
   BOOST_STATIC_ASSERT(std::numeric_limits<T>::is_specialized);
   BOOST_STATIC_ASSERT(std::numeric_limits<T>::radix != 2);
   BOOST_MATH_STD_USING
   int expon;
   static const char* function = "ulp<%1%>(%1%)";

   int fpclass = (boost::math::fpclassify)(val);

   if(fpclass == (int)FP_NAN)
   {
      return policies::raise_domain_error<T>(
         function,
         "Argument must be finite, but got %1%", val, pol);
   }
   else if((fpclass == (int)FP_INFINITE) || (fabs(val) >= tools::max_value<T>()))
   {
      return (val < 0 ? -1 : 1) * policies::raise_overflow_error<T>(function, 0, pol);
   }
   else if(fpclass == FP_ZERO)
      return detail::get_smallest_value<T>();
   //
   // This code is almost the same as that for float_next, except for negative integers,
   // where we preserve the relation ulp(x) == ulp(-x) as does Java:
   //
   expon = 1 + ilogb(fabs(val));
   T diff = scalbn(T(1), expon - std::numeric_limits<T>::digits);
   if(diff == 0)
      diff = detail::get_smallest_value<T>();
   return diff;
}

vmod_hash_backend(const struct vrt_ctx *ctx, struct vmod_directors_hash *rr,
    const char *arg, ...)
{
	struct SHA256Context sha_ctx;
	va_list ap;
	const char *p;
	unsigned char sha256[SHA256_LEN];
	VCL_BACKEND be;
	double r;

	CHECK_OBJ_NOTNULL(ctx, VRT_CTX_MAGIC);

	CHECK_OBJ_NOTNULL(rr, VMOD_DIRECTORS_HASH_MAGIC);
	SHA256_Init(&sha_ctx);
	va_start(ap, arg);
	p = arg;
	while (p != vrt_magic_string_end) {
		SHA256_Update(&sha_ctx, arg, strlen(arg));
		p = va_arg(ap, const char *);
	}
	va_end(ap);
	SHA256_Final(sha256, &sha_ctx);

	r = vbe32dec(sha256);
	r = scalbn(r, -32);
	assert(r >= 0 && r <= 1.0);
	be = vdir_pick_be(rr->vd, r, rr->nloops);
	return (be);
}

	double ldexp(double value, int exp)
{
	if(!finite(value)||value==0.0) return value;
	value = scalbn(value,exp);
	if(!finite(value)||value==0.0) errno = ERANGE;
	return value;
}

double ldexp(double value, int exp)
{
    if(!finite(value)||value==0.0) return value;
    value = scalbn(value,exp);
    if(!finite(value)||value==0.0) libm_errno = 34;
    return value;
}

int main(void)
{
	#pragma STDC FENV_ACCESS ON
	double y;
	float d;
	int e, i, err = 0;
	struct di_d *p;

	for (i = 0; i < sizeof t/sizeof *t; i++) {
		p = t + i;

		if (p->r < 0)
			continue;
		fesetround(p->r);
		feclearexcept(FE_ALL_EXCEPT);
		y = scalbn(p->x, p->i);
		e = fetestexcept(INEXACT|INVALID|DIVBYZERO|UNDERFLOW|OVERFLOW);

		if (!checkexceptall(e, p->e, p->r)) {
			printf("%s:%d: bad fp exception: %s scalbn(%a, %lld)=%a, want %s",
				p->file, p->line, rstr(p->r), p->x, p->i, p->y, estr(p->e));
			printf(" got %s\n", estr(e));
			err++;
		}
		d = ulperr(y, p->y, p->dy);
		if (!checkcr(y, p->y, p->r)) {
			printf("%s:%d: %s scalbn(%a, %lld) want %a got %a, ulperr %.3f = %a + %a\n",
				p->file, p->line, rstr(p->r), p->x, p->i, p->y, y, d, d-p->dy, p->dy);
			err++;
		}
	}
	return !!err;
}

double scalbln(double x, long n) {
    if (n > INT_MAX)
        n = INT_MAX;
    else if (n < INT_MIN)
        n = INT_MIN;
    return scalbn(x, n);
}

double
exp2(double x) {
	int	ix, hx, k;
	double	t;

	ix = ((int *)&x)[HIWORD];
	hx = ix & ~0x80000000;

	if (hx >= 0x4090e000) {	/* |x| >= 1080 or x is nan */
		if (hx >= 0x7ff00000) {	/* x is inf or nan */
			if (ix == 0xfff00000 && ((int *)&x)[LOWORD] == 0)
				return (zero);
			return (x * x);
		}
		t = (ix < 0)? tiny : huge;
		return (t * t);
	}

	if (hx < 0x3fe00000) {	/* |x| < 0.5 */
		if (hx < 0x3c000000)
			return (one + x);
		return (exp(ln2 * x));
	}

	k = (int)x;
	if (x != (double)k)
		k = (int)((ix < 0)? x - half : x + half);
	return (scalbn(exp(ln2 * (x - (double)k)), k));
}

/*
 * Compute the hypotenuse of a right triangle, avoiding intermediate
 * overflow or underflow.
 *
 * (This example ignores the case of one argument having
 * great magnitude and the other small, causing both overflow
 * and underflow!)
 */
double hypotenuse(double sidea, double sideb)
{
#pragma STDC FENV_ACCESS ON
    double sum, scale, ascaled, bscaled, invscale;
    fenv_t fpenv;
    int fpeflags;

    if ( signbit(sidea))  sidea = fabs(sidea);
    if ( signbit(sideb))  sideb = fabs(sideb);

    feholdexcept(&fpenv);        // Save previous environment,
                                 // clear exceptions,
                                 // switch to nonstop processing.
    invscale = 1.0;
    sum = sidea * sidea + sideb * sideb;    // First try whether a^2 + b^2
                                            // causes any exceptions.

    fpeflags = fetestexcept(FE_UNDERFLOW | FE_OVERFLOW);    // Did it?
    if (fpeflags & FE_OVERFLOW && sidea > 1.0 && sideb > 1.0)
    {
        /* a^2 + b^2 caused an overflow. Scale the triangle down. */
        feclearexcept(FE_OVERFLOW);
        scale = scalbn( 1.0, (DBL_MIN_EXP / 2));

        invscale = 1.0 / scale;
        ascaled = scale * sidea;
        bscaled = scale * sideb;
        sum = ascaled * ascaled + bscaled * bscaled;
    }
    else if (fpeflags & FE_UNDERFLOW && sidea < 1.0 && sideb < 1.0)
    {
        /* a^2 + b^2 caused an underflow. Scale the triangle up. */
        feclearexcept(FE_UNDERFLOW);
        scale = scalbn( 1.0, (DBL_MAX_EXP / 2));

        invscale = 1.0 / scale;
        ascaled = scale * sidea;
        bscaled = scale * sideb;
        sum = ascaled * ascaled + bscaled * bscaled;
    }

    feupdateenv(&fpenv);     // restore the caller's environment, and
                             // raise any new exceptions

    /* c = (1/scale) * sqrt((a * scale)^2 + (b * scale)^2): */
    return invscale * sqrt(sum);
}

long double scalbnl(long double x, int exp)
{
#if defined(__arm__) || defined(_ARM_)
    return scalbn(x, exp);
#else
#error Not supported on your platform yet
#endif
}

double
ldexp(double value, int exp0)
{
	if(!finite(value)||value==0.0) return value;
	value = scalbn(value,exp0);
	/*if(!finite(value)||value==0.0) errno = POK_ERRNO_ERANGE;*/
	return value;
}

void
Math_scalbn(void *fp)
{
	F_Math_scalbn *f;

	f = fp;

	*f->ret = scalbn(f->x, f->n);
}