Merge pull request #69 from JuliaLang/vs/longdouble

Long double versions of math functions

Author: ViralBShah

Make.inc

@@ -3,8 +3,8 @@
 OS := $(shell uname)
 # Do not forget to bump SOMINOR when changing VERSION,
 # and SOMAJOR when breaking ABI in a backward-incompatible way
-VERSION = 0.4
-SOMAJOR = 1
+VERSION = 0.5
+SOMAJOR = 2
 SOMINOR = 0
 DESTDIR =
 prefix = /usr/local
@@ -64,7 +64,7 @@ endif
 
 ifeq ($(ARCH),x86_64)
 override ARCH := amd64
-endif 
+endif
 
 ifneq (,$(findstring MINGW,$(OS)))
 override OS=WINNT

ld128/e_acoshl.c

@@ -0,0 +1,58 @@
+/* @(#)e_acosh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/* acoshl(x)
+ * Method :
+ *	Based on
+ *		acoshl(x) = logl [ x + sqrtl(x*x-1) ]
+ *	we have
+ *		acoshl(x) := logl(x)+ln2,	if x is large; else
+ *		acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else
+ *		acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1.
+ *
+ * Special cases:
+ *	acoshl(x) is NaN with signal if x<1.
+ *	acoshl(NaN) is NaN without signal.
+ */
+
+#include <math.h>
+
+#include "math_private.h"
+
+static const long double
+one	= 1.0,
+ln2	= 0.6931471805599453094172321214581766L;
+
+long double
+acoshl(long double x)
+{
+	long double t;
+	u_int64_t lx;
+	int64_t hx;
+	GET_LDOUBLE_WORDS64(hx,lx,x);
+	if(hx<0x3fff000000000000LL) {		/* x < 1 */
+	    return (x-x)/(x-x);
+	} else if(hx >=0x4035000000000000LL) {	/* x > 2**54 */
+	    if(hx >=0x7fff000000000000LL) {	/* x is inf of NaN */
+		return x+x;
+	    } else
+		return logl(x)+ln2;	/* acoshl(huge)=logl(2x) */
+	} else if(((hx-0x3fff000000000000LL)|lx)==0) {
+	    return 0.0L;			/* acosh(1) = 0 */
+	} else if (hx > 0x4000000000000000LL) {	/* 2**28 > x > 2 */
+	    t=x*x;
+	    return logl(2.0L*x-one/(x+sqrtl(t-one)));
+	} else {			/* 1<x<2 */
+	    t = x-one;
+	    return log1pl(t+sqrtl(2.0L*t+t*t));
+	}
+}

ld128/e_atanhl.c

@@ -0,0 +1,65 @@
+/* @(#)e_atanh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/* atanhl(x)
+ * Method :
+ *    1.Reduced x to positive by atanh(-x) = -atanh(x)
+ *    2.For x>=0.5
+ *                   1              2x                          x
+ *	atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
+ *                   2             1 - x                      1 - x
+ *
+ * 	For x<0.5
+ *	atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x))
+ *
+ * Special cases:
+ *	atanhl(x) is NaN if |x| > 1 with signal;
+ *	atanhl(NaN) is that NaN with no signal;
+ *	atanhl(+-1) is +-INF with signal.
+ *
+ */
+
+#include <math.h>
+
+#include "math_private.h"
+
+static const long double one = 1.0L, huge = 1e4900L;
+
+static const long double zero = 0.0L;
+
+long double
+atanhl(long double x)
+{
+	long double t;
+	u_int32_t jx, ix;
+	ieee_quad_shape_type u;
+
+	u.value = x;
+	jx = u.parts32.mswhi;
+	ix = jx & 0x7fffffff;
+	u.parts32.mswhi = ix;
+	if (ix >= 0x3fff0000) /* |x| >= 1.0 or infinity or NaN */
+	  {
+	    if (u.value == one)
+	      return x/zero;
+	    else
+	      return (x-x)/(x-x);
+	  }
+	if(ix<0x3fc60000 && (huge+x)>zero) return x;	/* x < 2^-57 */
+
+	if(ix<0x3ffe0000) {		/* x < 0.5 */
+	    t = u.value+u.value;
+	    t = 0.5*log1pl(t+t*u.value/(one-u.value));
+	} else
+	    t = 0.5*log1pl((u.value+u.value)/(one-u.value));
+	if(jx & 0x80000000) return -t; else return t;
+}

ld128/e_coshl.c

@@ -0,0 +1,105 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/* coshl(x)
+ * Method :
+ * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
+ *      1. Replace x by |x| (coshl(x) = coshl(-x)).
+ *      2.
+ *                                                      [ exp(x) - 1 ]^2
+ *          0        <= x <= ln2/2  :  coshl(x) := 1 + -------------------
+ *                                                         2*exp(x)
+ *
+ *                                                 exp(x) +  1/exp(x)
+ *          ln2/2    <= x <= 22     :  coshl(x) := -------------------
+ *                                                         2
+ *          22       <= x <= lnovft :  coshl(x) := expl(x)/2
+ *          lnovft   <= x <= ln2ovft:  coshl(x) := expl(x/2)/2 * expl(x/2)
+ *          ln2ovft  <  x           :  coshl(x) := huge*huge (overflow)
+ *
+ * Special cases:
+ *      coshl(x) is |x| if x is +INF, -INF, or NaN.
+ *      only coshl(0)=1 is exact for finite x.
+ */
+
+#include <math.h>
+
+#include "math_private.h"
+
+static const long double one = 1.0, half = 0.5, huge = 1.0e4900L,
+ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
+
+long double
+coshl(long double x)
+{
+  long double t, w;
+  int32_t ex;
+  ieee_quad_shape_type u;
+
+  u.value = x;
+  ex = u.parts32.mswhi & 0x7fffffff;
+
+  /* Absolute value of x.  */
+  u.parts32.mswhi = ex;
+
+  /* x is INF or NaN */
+  if (ex >= 0x7fff0000)
+    return x * x;
+
+  /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
+  if (ex < 0x3ffd62e4) /* 0.3465728759765625 */
+    {
+      t = expm1l (u.value);
+      w = one + t;
+      if (ex < 0x3fb80000) /* |x| < 2^-116 */
+	return w;		/* cosh(tiny) = 1 */
+
+      return one + (t * t) / (w + w);
+    }
+
+  /* |x| in [0.5*ln2,40], return (exp(|x|)+1/exp(|x|)/2; */
+  if (ex < 0x40044000)
+    {
+      t = expl (u.value);
+      return half * t + half / t;
+    }
+
+  /* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
+  if (ex <= 0x400c62e3) /* 11356.375 */
+    return half * expl (u.value);
+
+  /* |x| in [log(maxdouble), overflowthresold] */
+  if (u.value <= ovf_thresh)
+    {
+      w = expl (half * u.value);
+      t = half * w;
+      return t * w;
+    }
+
+  /* |x| > overflowthresold, cosh(x) overflow */
+  return huge * huge;
+}

ld128/e_expl.c

@@ -0,0 +1,145 @@
+/*	$OpenBSD: e_expl.c,v 1.3 2013/11/12 20:35:18 martynas Exp $	*/
+
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+/*							expl.c
+ *
+ *	Exponential function, 128-bit long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, expl();
+ *
+ * y = expl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns e (2.71828...) raised to the x power.
+ *
+ * Range reduction is accomplished by separating the argument
+ * into an integer k and fraction f such that
+ *
+ *     x    k  f
+ *    e  = 2  e.
+ *
+ * A Pade' form of degree 2/3 is used to approximate exp(f) - 1
+ * in the basic range [-0.5 ln 2, 0.5 ln 2].
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      +-MAXLOG    100,000     2.6e-34     8.6e-35
+ *
+ *
+ * Error amplification in the exponential function can be
+ * a serious matter.  The error propagation involves
+ * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
+ * which shows that a 1 lsb error in representing X produces
+ * a relative error of X times 1 lsb in the function.
+ * While the routine gives an accurate result for arguments
+ * that are exactly represented by a long double precision
+ * computer number, the result contains amplified roundoff
+ * error for large arguments not exactly represented.
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * exp underflow    x < MINLOG         0.0
+ * exp overflow     x > MAXLOG         MAXNUM
+ *
+ */
+
+/*	Exponential function	*/
+
+#include <float.h>
+#include <math.h>
+
+#include "math_private.h"
+
+/* Pade' coefficients for exp(x) - 1
+   Theoretical peak relative error = 2.2e-37,
+   relative peak error spread = 9.2e-38
+ */
+static long double P[5] = {
+ 3.279723985560247033712687707263393506266E-10L,
+ 6.141506007208645008909088812338454698548E-7L,
+ 2.708775201978218837374512615596512792224E-4L,
+ 3.508710990737834361215404761139478627390E-2L,
+ 9.999999999999999999999999999999999998502E-1L
+};
+static long double Q[6] = {
+ 2.980756652081995192255342779918052538681E-12L,
+ 1.771372078166251484503904874657985291164E-8L,
+ 1.504792651814944826817779302637284053660E-5L,
+ 3.611828913847589925056132680618007270344E-3L,
+ 2.368408864814233538909747618894558968880E-1L,
+ 2.000000000000000000000000000000000000150E0L
+};
+/* C1 + C2 = ln 2 */
+static const long double C1 = -6.93145751953125E-1L;
+static const long double C2 = -1.428606820309417232121458176568075500134E-6L;
+
+static const long double LOG2EL = 1.442695040888963407359924681001892137426646L;
+static const long double MAXLOGL = 1.1356523406294143949491931077970764891253E4L;
+static const long double MINLOGL = -1.143276959615573793352782661133116431383730e4L;
+static const long double huge = 0x1p10000L;
+#if 0 /* XXX Prevent gcc from erroneously constant folding this. */
+static const long double twom10000 = 0x1p-10000L;
+#else
+static volatile long double twom10000 = 0x1p-10000L;
+#endif
+
+long double
+expl(long double x)
+{
+long double px, xx;
+int n;
+
+if( x > MAXLOGL)
+	return (huge*huge);		/* overflow */
+
+if( x < MINLOGL )
+	return (twom10000*twom10000);	/* underflow */
+
+/* Express e**x = e**g 2**n
+ *   = e**g e**( n loge(2) )
+ *   = e**( g + n loge(2) )
+ */
+px = floorl( LOG2EL * x + 0.5L ); /* floor() truncates toward -infinity. */
+n = px;
+x += px * C1;
+x += px * C2;
+/* rational approximation for exponential
+ * of the fractional part:
+ * e**x =  1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
+ */
+xx = x * x;
+px = x * __polevll( xx, P, 4 );
+xx = __polevll( xx, Q, 5 );
+x =  px/( xx - px );
+x = 1.0L + x + x;
+
+x = ldexpl( x, n );
+return(x);
+}