/

2010-11-13  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
            Tobias Burnus  <burnus@net-b.de>

        PR fortran/32049
        * Makefile.def: Add libquadmath; build it with language=fortran.
        * configure.ac: Add libquadmath.
        * Makefile.tpl: Handle multiple libs in check-[+language+].
        * Makefile.in: Regenerate.
        * configure: Regenerate.

libquadmath/
2010-11-13  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
            Tobias Burnus  <burnus@net-b.de>

        PR fortran/32049
        Initial implementation and checkin.

gcc/fortran/
2010-11-13  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
            Tobias Burnus  <burnus@net-b.de>

        PR fortran/32049
        * gfortranspec.c (find_spec_file): New function.
        (lang_specific_driver): Try to find .spec file and use it.
        * trans-io.c (iocall): Define
        * IOCALL_X_REAL128/COMPLEX128(,write).
        (gfc_build_io_library_fndecls): Build decl for __float128 I/O.
        (transfer_expr): Call __float128 I/O functions.
        * trans-types.c (gfc_init_kinds): Allow kind-16 belonging
        to __float128.

gcc/testsuite/
2010-11-13  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
            Tobias Burnus  <burnus@net-b.de>

        PR fortran/32049
        * gfortran.dg/quad_1.f90: New.
        * lib/gcc-defs.exp (gcc-set-multilib-library-path): Use also
        compiler arguments.
        * lib/gfortran.exp (gfortran_link_flags): Add libquadmath to
        library search path; call gcc-set-multilib-library-path with
        arguments such that libgfortran.spec is found.
        (gfortran_init): Add path for libgfortran.spec to
GFORTRAN_UNDER_TEST.

libgomp/
2010-11-13  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
            Tobias Burnus  <burnus@net-b.de>

        PR fortran/32049
        * configure.ac: 
        * configure: Regenerate.

libgfortran/
2010-11-13  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
            Tobias Burnus  <burnus@net-b.de>

        PR fortran/32049
        * Makefile.am: Add missing pow_r16_i4.c, add transfer128.c,
        link libquadmath, if used.
        * acinclude.m4 (LIBGFOR_CHECK_FLOAT128): Add.
        * configure.ac: Use it, touch spec file.
        * gfortran.map: Add pow_r16_i4 and
        transfer_(real,complex)128(,write) functions.
        * intrinsics/cshift0.c (cshift0): Handle __float128 type.
        * intrinsics/erfc_scaled_inc.c: Ditto.
        * intrinsics/pack_generic.c (pack): Ditto
        * intrinsics/spread_generic.c (spread): Ditto.
        * intrinsics/unpack_generic.c (unpack1): Ditto.
        * io/read.c (convert_real): Ditto.
        * io/transfer.c: Update comments.
        * io/transfer128.c: New file.
        * io/write_float.def (write_float): Handle __float128 type.
        * libgfortran.h: #include quadmath_weak.h, define __builtin_infq
        and nanq.
        * m4/mtype.m4: Handle __float128 type.
        * runtime/in_pack_generic.c (internal_pack): Ditto.
        * runtime/in_unpack_generic.c (internal_unpack): Ditto.
        * kinds-override.h: New file.
        * libgfortran.spec.in: Ditto.
        * generated/pow_r16_i4.c: Generated.
        * Makefile.in: Regenerate.
        * configure: Regenerate.
        * config.h: Regenerate.
        * bessel_r10.c: Regenerate.
        * bessel_r16.c: Regenerate.
        * bessel_r4.c: Regenerate.
        * bessel_r8.c: Regenerate.
        * exponent_r16.c: Regenerate.
        * fraction_r16.c: Regenerate.
        * nearest_r16.c: Regenerate.
        * norm2_r10.c: Regenerate.
        * norm2_r16.c: Regenerate.
        * norm2_r4.c: Regenerate.
        * norm2_r8.c: Regenerate.
        * rrspacing_r16.c: Regenerate.
        * set_exponent_r16.c: Regenerate.
        * spacing_r16.c: Regenerate.



git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@166825 138bc75d-0d04-0410-961f-82ee72b054a4

1 files changed, 124 insertions, 0 deletions

diff --git a/libquadmath/math/hypotq.c b/libquadmath/math/hypotq.c
new file mode 100644
index 00000000000..e4c2d1bc214
--- /dev/null
+++ b/libquadmath/math/hypotq.c
@@ -0,0 +1,124 @@
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/* From e_hypotl.c -- long double version of e_hypot.c.
+ * Conversion to long double by Jakub Jelinek, jakub@redhat.com.
+ * Conversion to __float128 by FX Coudert, fxcoudert@gcc.gnu.org.
+ */
+
+/* hypotq(x,y)
+ *
+ * Method :
+ *	If (assume round-to-nearest) z=x*x+y*y
+ *	has error less than sqrtl(2)/2 ulp, than
+ *	sqrtl(z) has error less than 1 ulp (exercise).
+ *
+ *	So, compute sqrtl(x*x+y*y) with some care as
+ *	follows to get the error below 1 ulp:
+ *
+ *	Assume x>y>0;
+ *	(if possible, set rounding to round-to-nearest)
+ *	1. if x > 2y  use
+ *		x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
+ *	where x1 = x with lower 64 bits cleared, x2 = x-x1; else
+ *	2. if x <= 2y use
+ *		t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
+ *	where t1 = 2x with lower 64 bits cleared, t2 = 2x-t1,
+ *	y1= y with lower 64 bits chopped, y2 = y-y1.
+ *
+ *	NOTE: scaling may be necessary if some argument is too
+ *	      large or too tiny
+ *
+ * Special cases:
+ *	hypotq(x,y) is INF if x or y is +INF or -INF; else
+ *	hypotq(x,y) is NAN if x or y is NAN.
+ *
+ * Accuracy:
+ * 	hypotq(x,y) returns sqrtl(x^2+y^2) with error less
+ * 	than 1 ulps (units in the last place)
+ */
+
+#include "quadmath-imp.h"
+
+__float128
+hypotq (__float128 x, __float128 y)
+{
+  __float128 a, b, t1, t2, y1, y2, w;
+  int64_t j, k, ha, hb;
+
+  GET_FLT128_MSW64(ha,x);
+  ha &= 0x7fffffffffffffffLL;
+  GET_FLT128_MSW64(hb,y);
+  hb &= 0x7fffffffffffffffLL;
+  if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
+  SET_FLT128_MSW64(a,ha);	/* a <- |a| */
+  SET_FLT128_MSW64(b,hb);	/* b <- |b| */
+  if((ha-hb)>0x78000000000000LL) {return a+b;} /* x/y > 2**120 */
+  k=0;
+  if(ha > 0x5f3f000000000000LL) {	/* a>2**8000 */
+     if(ha >= 0x7fff000000000000LL) {	/* Inf or NaN */
+         uint64_t low;
+         w = a+b;			/* for sNaN */
+         GET_FLT128_LSW64(low,a);
+         if(((ha&0xffffffffffffLL)|low)==0) w = a;
+         GET_FLT128_LSW64(low,b);
+         if(((hb^0x7fff000000000000LL)|low)==0) w = b;
+         return w;
+     }
+     /* scale a and b by 2**-9600 */
+     ha -= 0x2580000000000000LL;
+     hb -= 0x2580000000000000LL;	k += 9600;
+     SET_FLT128_MSW64(a,ha);
+     SET_FLT128_MSW64(b,hb);
+  }
+  if(hb < 0x20bf000000000000LL) {	/* b < 2**-8000 */
+      if(hb <= 0x0000ffffffffffffLL) {	/* subnormal b or 0 */
+          uint64_t low;
+  	GET_FLT128_LSW64(low,b);
+  	if((hb|low)==0) return a;
+  	t1=0;
+  	SET_FLT128_MSW64(t1,0x7ffd000000000000LL); /* t1=2^16382 */
+  	b *= t1;
+  	a *= t1;
+  	k -= 16382;
+      } else {		/* scale a and b by 2^9600 */
+          ha += 0x2580000000000000LL; 	/* a *= 2^9600 */
+  	hb += 0x2580000000000000LL;	/* b *= 2^9600 */
+  	k -= 9600;
+  	SET_FLT128_MSW64(a,ha);
+  	SET_FLT128_MSW64(b,hb);
+      }
+  }
+    /* medium size a and b */
+  w = a-b;
+  if (w>b) {
+      t1 = 0;
+      SET_FLT128_MSW64(t1,ha);
+      t2 = a-t1;
+      w  = sqrtq(t1*t1-(b*(-b)-t2*(a+t1)));
+  } else {
+      a  = a+a;
+      y1 = 0;
+      SET_FLT128_MSW64(y1,hb);
+      y2 = b - y1;
+      t1 = 0;
+      SET_FLT128_MSW64(t1,ha+0x0001000000000000LL);
+      t2 = a - t1;
+      w  = sqrtq(t1*y1-(w*(-w)-(t1*y2+t2*b)));
+  }
+  if(k!=0) {
+      uint64_t high;
+      t1 = 1.0L;
+      GET_FLT128_MSW64(high,t1);
+      SET_FLT128_MSW64(t1,high+(k<<48));
+      return t1*w;
+  } else return w;
+}