* [patch,libgfortran] Fix binary128 ERFC_SCALED
@ 2013-11-17 21:41 FX
2013-11-20 20:13 ` Steve Kargl
0 siblings, 1 reply; 8+ messages in thread
From: FX @ 2013-11-17 21:41 UTC (permalink / raw)
To: gfortran; +Cc: gcc-patches
[-- Attachment #1: Type: text/plain, Size: 575 bytes --]
This patch fixes libgfortran’s binary128 [aka real(kind=16)] variant of ERFC_SCALED. The original code, which I had lifted from netlib, gives only 18 significant decimal digits, which is not enough for binary128 (33 decimal digits).
I thus implemented a new variant for binary128. For arguments < 12, it simply calls erfcq() then multiplies by expq(x*x). For larger arguments, it uses a power expansion in 1/x. The new implementation provides answers within to 2 ulp of the correct value.
Regtested on x86_64-apple-darwin13, comes with a testcase. OK to commit?
FX
[-- Attachment #2: erfc_scaled.ChangeLog --]
[-- Type: application/octet-stream, Size: 344 bytes --]
2013-11-17 Francois-Xavier Coudert <fxcoudert@gcc.gnu.org>
PR libfortran/49024
* intrinsics/erfc_scaled.c (erfc_scaled_r16): New function.
* intrinsics/erfc_scaled_inc.c: Do not provide quadruple
precision variant.
2013-11-17 Francois-Xavier Coudert <fxcoudert@gcc.gnu.org>
PR libfortran/49024
* gfortran.dg/erf_3.F90: New file.
[-- Attachment #3: erfc_scaled.diff --]
[-- Type: application/octet-stream, Size: 3174 bytes --]
Index: libgfortran/intrinsics/erfc_scaled.c
===================================================================
--- libgfortran/intrinsics/erfc_scaled.c (revision 204446)
+++ libgfortran/intrinsics/erfc_scaled.c (working copy)
@@ -45,8 +45,55 @@ see the files COPYING3 and COPYING.RUNTI
#include "erfc_scaled_inc.c"
#endif
-#ifdef HAVE_GFC_REAL_16
+#if defined(HAVE_GFC_REAL_16) && defined(GFC_REAL_16_IS_LONG_DOUBLE)
#undef KIND
#define KIND 16
#include "erfc_scaled_inc.c"
#endif
+
+
+/* For quadruple-precision (__float128), netlib's implementation is
+ not accurate enough. We provide another one. */
+
+
+extern GFC_REAL_16 erfc_scaled_r16 (GFC_REAL_16);
+export_proto(erfc_scaled_r16);
+
+
+GFC_REAL_16
+erfc_scaled_r16 (GFC_REAL_16 x)
+{
+ if (x < 12)
+ {
+ /* Compute directly as ERFC_SCALED(x) = ERFC(x) * EXP(X**2).
+ This is not perfect, but much better than netlib. */
+ return erfcq(x) * expq(x*x);
+ }
+ else
+ {
+ /* Calculate ERFC_SCALED(x) using a power series in 1/x:
+ ERFC_SCALED(x) = 1 / (x * sqrt(pi))
+ * (1 + Sum_n (-1)**n * (1 * 3 * 5 * ... * (2n-1))
+ / (2 * x**2)**n)
+ */
+ GFC_REAL_16 sum = 0, oldsum;
+ GFC_REAL_16 inv2x2 = 1 / (2 * x * x);
+ GFC_REAL_16 fac = 1;
+ int n = 1;
+
+ while (n < 200)
+ {
+ fac *= - (2*n - 1) * inv2x2;
+ oldsum = sum;
+ sum += fac;
+
+ if (sum == oldsum)
+ break;
+
+ n++;
+ }
+
+ return (1 + sum) / x * (M_2_SQRTPIq / 2);
+ }
+}
+
Index: libgfortran/intrinsics/erfc_scaled_inc.c
===================================================================
--- libgfortran/intrinsics/erfc_scaled_inc.c (revision 204446)
+++ libgfortran/intrinsics/erfc_scaled_inc.c (working copy)
@@ -48,11 +48,6 @@ see the files COPYING3 and COPYING.RUNTI
# define TRUNC(x) truncl(x)
# endif
-#elif (KIND == 16 && defined(GFC_REAL_16_IS_FLOAT128))
-
-# define EXP(x) expq(x)
-# define TRUNC(x) truncq(x)
-
#else
# error "What exactly is it that you want me to do?"
Index: gcc/testsuite/gfortran.dg/erf_3.F90
===================================================================
--- gcc/testsuite/gfortran.dg/erf_3.F90 (revision 0)
+++ gcc/testsuite/gfortran.dg/erf_3.F90 (working copy)
@@ -0,0 +1,44 @@
+! { dg-do run }
+! { dg-options "-fno-range-check -ffree-line-length-none -O0" }
+! { dg-add-options ieee }
+!
+! Check that simplification functions and runtime library agree on ERF,
+! ERFC and ERFC_SCALED, for quadruple-precision.
+
+program test
+ implicit none
+
+ real(kind=16) :: x16
+
+#define CHECK(a) \
+ x16 = a ; \
+ call check(erf(real(a,kind=16)), erf(x16)) ; \
+ call check(erfc(real(a,kind=16)), erfc(x16)) ; \
+ call check(erfc_scaled(real(a,kind=16)), erfc_scaled(x16))
+
+ CHECK(0.0)
+ CHECK(0.9)
+ CHECK(1.9)
+ CHECK(10.)
+ CHECK(11.)
+ CHECK(12.)
+ CHECK(13.)
+ CHECK(14.)
+ CHECK(49.)
+ CHECK(190.)
+
+ CHECK(-0.0)
+ CHECK(-0.9)
+ CHECK(-1.9)
+ CHECK(-19.)
+ CHECK(-190.)
+
+contains
+
+ subroutine check (a, b)
+ real(kind=16), intent(in) :: a, b
+ print *, abs(a-b) / spacing(a)
+ if (abs(a - b) > 10 * spacing(a)) call abort
+ end subroutine
+
+end program test
^ permalink raw reply [flat|nested] 8+ messages in thread
* Re: [patch,libgfortran] Fix binary128 ERFC_SCALED
2013-11-17 21:41 [patch,libgfortran] Fix binary128 ERFC_SCALED FX
@ 2013-11-20 20:13 ` Steve Kargl
2013-11-20 22:19 ` FX
0 siblings, 1 reply; 8+ messages in thread
From: Steve Kargl @ 2013-11-20 20:13 UTC (permalink / raw)
To: FX; +Cc: gfortran, gcc-patches
On Sun, Nov 17, 2013 at 11:29:50AM +0100, FX wrote:
> This patch fixes libgfortran?s binary128 [aka real(kind=16)] variant of ERFC_SCALED. The original code, which I had lifted from netlib, gives only 18 significant decimal digits, which is not enough for binary128 (33 decimal digits).
>
> I thus implemented a new variant for binary128. For arguments < 12, it simply calls erfcq() then multiplies by expq(x*x). For larger arguments, it uses a power expansion in 1/x. The new implementation provides answers within to 2 ulp of the correct value.
>
> Regtested on x86_64-apple-darwin13, comes with a testcase. OK to commit?
> FX
>
This is probably ok as it fixes a 2**64ish ULP problem.
There is a missed optimization in
+ if (x < 12)
+ {
+ /* Compute directly as ERFC_SCALED(x) = ERFC(x) * EXP(X**2).
+ This is not perfect, but much better than netlib. */
+ return erfcq(x) * expq(x*x);
+ }
If x is less than approximately -8192, then erfc(x)*exp(x*x)
overflows. Something like the following shortcircuits the
2 function calls and their possible argument reduction steps.
volatile GFC_REAL_16 zero = 0;
if (x < 12)
{
if (x < some_overflow_value_to_be_determined)
return (1 / zero);
return erfcq(x) * expq(x*x);
}
PS: There is missing whitespace: x*x should be x * x.
--
Steve
^ permalink raw reply [flat|nested] 8+ messages in thread
* Re: [patch,libgfortran] Fix binary128 ERFC_SCALED
2013-11-20 20:13 ` Steve Kargl
@ 2013-11-20 22:19 ` FX
2013-11-20 23:02 ` Steve Kargl
2013-11-21 10:32 ` Andreas Schwab
0 siblings, 2 replies; 8+ messages in thread
From: FX @ 2013-11-20 22:19 UTC (permalink / raw)
To: Steve Kargl; +Cc: gfortran, gcc-patches
[-- Attachment #1: Type: text/plain, Size: 495 bytes --]
> There is a missed optimization in
>
> + if (x < 12)
> + {
> + /* Compute directly as ERFC_SCALED(x) = ERFC(x) * EXP(X**2).
> + This is not perfect, but much better than netlib. */
> + return erfcq(x) * expq(x*x);
> + }
>
> If x is less than approximately -8192, then erfc(x)*exp(x*x)
> overflows.
Hum, I get roughly -106 where you have -8192, so I’ll not commit immediately with the value I have, and let you check it first.
OK to commit?
FX
[-- Attachment #2: erfc_scaled.diff --]
[-- Type: application/octet-stream, Size: 3274 bytes --]
Index: libgfortran/intrinsics/erfc_scaled.c
===================================================================
--- libgfortran/intrinsics/erfc_scaled.c (revision 205019)
+++ libgfortran/intrinsics/erfc_scaled.c (working copy)
@@ -45,8 +45,59 @@ see the files COPYING3 and COPYING.RUNTI
#include "erfc_scaled_inc.c"
#endif
-#ifdef HAVE_GFC_REAL_16
+#if defined(HAVE_GFC_REAL_16) && defined(GFC_REAL_16_IS_LONG_DOUBLE)
#undef KIND
#define KIND 16
#include "erfc_scaled_inc.c"
#endif
+
+
+/* For quadruple-precision (__float128), netlib's implementation is
+ not accurate enough. We provide another one. */
+
+
+extern GFC_REAL_16 erfc_scaled_r16 (GFC_REAL_16);
+export_proto(erfc_scaled_r16);
+
+
+GFC_REAL_16
+erfc_scaled_r16 (GFC_REAL_16 x)
+{
+ if (x < -106.566990228185312813205074546585730Q)
+ {
+ return __builtin_infq();
+ }
+ if (x < 12)
+ {
+ /* Compute directly as ERFC_SCALED(x) = ERFC(x) * EXP(X**2).
+ This is not perfect, but much better than netlib. */
+ return erfcq(x) * expq(x * x);
+ }
+ else
+ {
+ /* Calculate ERFC_SCALED(x) using a power series in 1/x:
+ ERFC_SCALED(x) = 1 / (x * sqrt(pi))
+ * (1 + Sum_n (-1)**n * (1 * 3 * 5 * ... * (2n-1))
+ / (2 * x**2)**n)
+ */
+ GFC_REAL_16 sum = 0, oldsum;
+ GFC_REAL_16 inv2x2 = 1 / (2 * x * x);
+ GFC_REAL_16 fac = 1;
+ int n = 1;
+
+ while (n < 200)
+ {
+ fac *= - (2*n - 1) * inv2x2;
+ oldsum = sum;
+ sum += fac;
+
+ if (sum == oldsum)
+ break;
+
+ n++;
+ }
+
+ return (1 + sum) / x * (M_2_SQRTPIq / 2);
+ }
+}
+
Index: libgfortran/intrinsics/erfc_scaled_inc.c
===================================================================
--- libgfortran/intrinsics/erfc_scaled_inc.c (revision 205019)
+++ libgfortran/intrinsics/erfc_scaled_inc.c (working copy)
@@ -48,11 +48,6 @@ see the files COPYING3 and COPYING.RUNTI
# define TRUNC(x) truncl(x)
# endif
-#elif (KIND == 16 && defined(GFC_REAL_16_IS_FLOAT128))
-
-# define EXP(x) expq(x)
-# define TRUNC(x) truncq(x)
-
#else
# error "What exactly is it that you want me to do?"
Index: gcc/testsuite/gfortran.dg/erf_3.F90
===================================================================
--- gcc/testsuite/gfortran.dg/erf_3.F90 (revision 0)
+++ gcc/testsuite/gfortran.dg/erf_3.F90 (working copy)
@@ -0,0 +1,44 @@
+! { dg-do run }
+! { dg-options "-fno-range-check -ffree-line-length-none -O0" }
+! { dg-add-options ieee }
+!
+! Check that simplification functions and runtime library agree on ERF,
+! ERFC and ERFC_SCALED, for quadruple-precision.
+
+program test
+ implicit none
+
+ real(kind=16) :: x16
+
+#define CHECK(a) \
+ x16 = a ; \
+ call check(erf(real(a,kind=16)), erf(x16)) ; \
+ call check(erfc(real(a,kind=16)), erfc(x16)) ; \
+ call check(erfc_scaled(real(a,kind=16)), erfc_scaled(x16))
+
+ CHECK(0.0)
+ CHECK(0.9)
+ CHECK(1.9)
+ CHECK(10.)
+ CHECK(11.)
+ CHECK(12.)
+ CHECK(13.)
+ CHECK(14.)
+ CHECK(49.)
+ CHECK(190.)
+
+ CHECK(-0.0)
+ CHECK(-0.9)
+ CHECK(-1.9)
+ CHECK(-19.)
+ CHECK(-190.)
+
+contains
+
+ subroutine check (a, b)
+ real(kind=16), intent(in) :: a, b
+ print *, abs(a-b) / spacing(a)
+ if (abs(a - b) > 10 * spacing(a)) call abort
+ end subroutine
+
+end program test
^ permalink raw reply [flat|nested] 8+ messages in thread
* Re: [patch,libgfortran] Fix binary128 ERFC_SCALED
2013-11-20 22:19 ` FX
@ 2013-11-20 23:02 ` Steve Kargl
2013-11-21 1:00 ` FX
2013-11-21 10:32 ` Andreas Schwab
1 sibling, 1 reply; 8+ messages in thread
From: Steve Kargl @ 2013-11-20 23:02 UTC (permalink / raw)
To: FX; +Cc: gfortran, gcc-patches
On Wed, Nov 20, 2013 at 09:13:14PM +0100, FX wrote:
> > There is a missed optimization in
> >
> > + if (x < 12)
> > + {
> > + /* Compute directly as ERFC_SCALED(x) = ERFC(x) * EXP(X**2).
> > + This is not perfect, but much better than netlib. */
> > + return erfcq(x) * expq(x*x);
> > + }
> >
> > If x is less than approximately -8192, then erfc(x)*exp(x*x)
> > overflows.
>
> Hum, I get roughly -106 where you have -8192, so I?ll not
> commit immediately with the value I have, and let you check it first.
>
> OK to commit?
>
The optimization can come later. Your current version is fine,
because it fixes the ULP issue.
As to the optimization, for large negative x, one has
efrc(x)*exp(x*x) --> 2 * exp(x*x) = 2**emax
exp(x*x) = 2**(emax-1)
x*x = (emax-1)*log(2)
x = - sqrt((emax-1)*log(2)) ! Choosing minus branch
x = -106.56... ! emax=16384 (could be off-by-1, did not check)
So, yeah, you're correct. My suggestion was based on the not
so careful mistake of replacing x*x by x+x and dropping log(2).
That is, I had x+x = -emax --> x = - emax / 2.
--
Steve
^ permalink raw reply [flat|nested] 8+ messages in thread
* Re: [patch,libgfortran] Fix binary128 ERFC_SCALED
2013-11-20 23:02 ` Steve Kargl
@ 2013-11-21 1:00 ` FX
0 siblings, 0 replies; 8+ messages in thread
From: FX @ 2013-11-21 1:00 UTC (permalink / raw)
To: gfortran, gcc-patches; +Cc: Steve Kargl
> So, yeah, you're correct. My suggestion was based on the not
> so careful mistake of replacing x*x by x+x and dropping log(2).
> That is, I had x+x = -emax --> x = - emax / 2.
Committed as rev. 205151, thanks for the review!
FX
^ permalink raw reply [flat|nested] 8+ messages in thread
* Re: [patch,libgfortran] Fix binary128 ERFC_SCALED
2013-11-20 22:19 ` FX
2013-11-20 23:02 ` Steve Kargl
@ 2013-11-21 10:32 ` Andreas Schwab
2013-11-21 10:43 ` FX
1 sibling, 1 reply; 8+ messages in thread
From: Andreas Schwab @ 2013-11-21 10:32 UTC (permalink / raw)
To: FX; +Cc: Steve Kargl, gfortran, gcc-patches
../../../libgfortran/intrinsics/erfc_scaled.c:59:1: error: unknown type name 'GFC_REAL_16'
extern GFC_REAL_16 erfc_scaled_r16 (GFC_REAL_16);
^
../../../libgfortran/intrinsics/erfc_scaled.c:59:1: warning: parameter names (without types) in function declaration [enabled by default]
../../../libgfortran/intrinsics/erfc_scaled.c:63:1: error: unknown type name 'GFC_REAL_16'
GFC_REAL_16
^
../../../libgfortran/intrinsics/erfc_scaled.c:64:18: error: unknown type name 'GFC_REAL_16'
erfc_scaled_r16 (GFC_REAL_16 x)
^
../../../libgfortran/intrinsics/erfc_scaled.c:66:3: error: unsupported non-standard suffix on floating constant
if (x < -106.566990228185312813205074546585730Q)
^
make[3]: *** [erfc_scaled.lo] Error 1
Andreas.
--
Andreas Schwab, SUSE Labs, schwab@suse.de
GPG Key fingerprint = 0196 BAD8 1CE9 1970 F4BE 1748 E4D4 88E3 0EEA B9D7
"And now for something completely different."
^ permalink raw reply [flat|nested] 8+ messages in thread
* Re: [patch,libgfortran] Fix binary128 ERFC_SCALED
2013-11-21 10:32 ` Andreas Schwab
@ 2013-11-21 10:43 ` FX
2013-11-21 11:43 ` Andreas Schwab
0 siblings, 1 reply; 8+ messages in thread
From: FX @ 2013-11-21 10:43 UTC (permalink / raw)
To: gfortran, gcc-patches; +Cc: Steve Kargl, Andreas Schwab
[-- Attachment #1: Type: text/plain, Size: 516 bytes --]
> ../../../libgfortran/intrinsics/erfc_scaled.c:59:1: error: unknown type name ‘GFC_REAL_16'
I’m really sorry about that, I should have tested on a system without binary128, that would have caught it.
Attached patch committed as rev. 205193 after checking it on system both with and without binary128.
Sorry again,
FX
2013-11-20 Francois-Xavier Coudert <fxcoudert@gcc.gnu.org>
PR libfortran/59227
* intrinsics/erfc_scaled.c (erfc_scaled_r16): Don't define if
__float128 is not available.
[-- Attachment #2: fix.diff --]
[-- Type: application/octet-stream, Size: 635 bytes --]
Index: intrinsics/erfc_scaled.c
===================================================================
--- intrinsics/erfc_scaled.c (revision 205151)
+++ intrinsics/erfc_scaled.c (working copy)
@@ -52,14 +52,14 @@ see the files COPYING3 and COPYING.RUNTI
#endif
+#ifdef GFC_REAL_16_IS_FLOAT128
+
/* For quadruple-precision (__float128), netlib's implementation is
not accurate enough. We provide another one. */
-
extern GFC_REAL_16 erfc_scaled_r16 (GFC_REAL_16);
export_proto(erfc_scaled_r16);
-
GFC_REAL_16
erfc_scaled_r16 (GFC_REAL_16 x)
{
@@ -101,3 +101,5 @@ erfc_scaled_r16 (GFC_REAL_16 x)
}
}
+#endif
+
^ permalink raw reply [flat|nested] 8+ messages in thread
* Re: [patch,libgfortran] Fix binary128 ERFC_SCALED
2013-11-21 10:43 ` FX
@ 2013-11-21 11:43 ` Andreas Schwab
0 siblings, 0 replies; 8+ messages in thread
From: Andreas Schwab @ 2013-11-21 11:43 UTC (permalink / raw)
To: FX; +Cc: gfortran, gcc-patches, Steve Kargl
FX <fxcoudert@gmail.com> writes:
> 2013-11-20 Francois-Xavier Coudert <fxcoudert@gcc.gnu.org>
>
> PR libfortran/59227
There is no connection to PR59227
Andreas.
--
Andreas Schwab, SUSE Labs, schwab@suse.de
GPG Key fingerprint = 0196 BAD8 1CE9 1970 F4BE 1748 E4D4 88E3 0EEA B9D7
"And now for something completely different."
^ permalink raw reply [flat|nested] 8+ messages in thread
end of thread, other threads:[~2013-11-21 9:24 UTC | newest]
Thread overview: 8+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2013-11-17 21:41 [patch,libgfortran] Fix binary128 ERFC_SCALED FX
2013-11-20 20:13 ` Steve Kargl
2013-11-20 22:19 ` FX
2013-11-20 23:02 ` Steve Kargl
2013-11-21 1:00 ` FX
2013-11-21 10:32 ` Andreas Schwab
2013-11-21 10:43 ` FX
2013-11-21 11:43 ` Andreas Schwab
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).