From patchwork Fri Mar 26 15:34:57 2021 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Paul Zimmermann X-Patchwork-Id: 1458858 Return-Path: X-Original-To: incoming@patchwork.ozlabs.org Delivered-To: patchwork-incoming@bilbo.ozlabs.org Authentication-Results: ozlabs.org; spf=pass (sender SPF authorized) smtp.mailfrom=sourceware.org (client-ip=2620:52:3:1:0:246e:9693:128c; helo=sourceware.org; envelope-from=libc-alpha-bounces@sourceware.org; receiver=) Received: from sourceware.org (server2.sourceware.org [IPv6:2620:52:3:1:0:246e:9693:128c]) (using TLSv1.3 with cipher TLS_AES_256_GCM_SHA384 (256/256 bits) key-exchange X25519 server-signature RSA-PSS (4096 bits) server-digest SHA256) (No client certificate requested) by ozlabs.org (Postfix) with ESMTPS id 4F6Qyr1hxBz9sRf for ; Sat, 27 Mar 2021 02:35:15 +1100 (AEDT) Received: from server2.sourceware.org (localhost [IPv6:::1]) by sourceware.org (Postfix) with ESMTP id 6106D385702C; Fri, 26 Mar 2021 15:35:12 +0000 (GMT) X-Original-To: libc-alpha@sourceware.org Delivered-To: libc-alpha@sourceware.org Received: from mail3-relais-sop.national.inria.fr (mail3-relais-sop.national.inria.fr [192.134.164.104]) by sourceware.org (Postfix) with ESMTPS id 3C4D83857C77 for ; Fri, 26 Mar 2021 15:35:02 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.3.2 sourceware.org 3C4D83857C77 Authentication-Results: sourceware.org; dmarc=none (p=none dis=none) header.from=inria.fr Authentication-Results: sourceware.org; spf=pass smtp.mailfrom=Paul.Zimmermann@loria.fr IronPort-HdrOrdr: A9a23:BZyULa/gJJD1wZzoqNRuk+HPdb1zdoIgy1knxilNYDZSddGVkN3roeQD2XbP5Ao5dVMBvZS7OKeGSW7B7pId2+IsFJqrQQWOggSVBa5v4YboyzfjXw3Sn9Qy6Y5OaK57YeeQMXFejdzmpCm0euxB/PCj8Lq0wdvY1WtnVwtwa6pthj0JaTqzN01tWU17AoAkH4CX/cpNq36Yf24LB/7LZEUte8jmi5n1mIn9YRgAbiRXjTWmqT+z8rb1H1y5034lIldy6I44/WLIngD/7KnLiYDd9jbm22Xe74tblbLao7N+LfeBk80Ebgjr4zzYArhJYaGIvzw+vYiUhWoCrd+kmX0dFvU2wWnYOkuypQLp0wnmzSxG0Q6a9XaoxV/qocn0X1sBerJ8rLMcVhPY7kom+Pt11LtC32yf3qAndi/orWDY79jMVxYvvE+1rX0++NRj9UB3YM8FaLQUrYF3xjIpLKs9 X-IronPort-AV: E=Sophos;i="5.81,280,1610406000"; d="scan'208";a="376975055" Received: from tomate.loria.fr ([152.81.10.51]) by mail3-relais-sop.national.inria.fr with ESMTP/TLS/DHE-RSA-AES256-GCM-SHA384; 26 Mar 2021 16:35:00 +0100 Received: from zimmerma by tomate.loria.fr with local (Exim 4.94) (envelope-from ) id 1lPoUF-007q0B-Ku; Fri, 26 Mar 2021 16:34:59 +0100 From: Paul Zimmermann To: libc-alpha@sourceware.org Subject: [PATCH] Improve the accuracy of tgamma (BZ #26983) Date: Fri, 26 Mar 2021 16:34:57 +0100 Message-Id: <20210326153457.1868146-1-Paul.Zimmermann@inria.fr> X-Mailer: git-send-email 2.30.2 MIME-Version: 1.0 X-Spam-Status: No, score=-10.0 required=5.0 tests=BAYES_00, GIT_PATCH_0, HEADER_FROM_DIFFERENT_DOMAINS, KAM_DMARC_STATUS, KAM_LOTSOFHASH, KAM_SHORT, RCVD_IN_MSPIKE_H3, RCVD_IN_MSPIKE_WL, SPF_HELO_NONE, SPF_PASS, TXREP autolearn=ham autolearn_force=no version=3.4.2 X-Spam-Checker-Version: SpamAssassin 3.4.2 (2018-09-13) on server2.sourceware.org X-BeenThere: libc-alpha@sourceware.org X-Mailman-Version: 2.1.29 Precedence: list List-Id: Libc-alpha mailing list List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Errors-To: libc-alpha-bounces@sourceware.org Sender: "Libc-alpha" With this patch, the maximal known error for tgamma is now reduced to 9 ulps for dbl-64, for all rounding modes. Since exhaustive testing is not possible for dbl-64, it might be that there are still cases with an error larger than 9 ulps, but all known cases are fixed (intensive tests were done to find cases with large errors). Tested on x86_64 and powerpc. --- math/auto-libm-test-in | 6 +- math/auto-libm-test-out-tgamma | 276 ++++++++++++++++++++++++ math/mul_split.h | 74 +++++++ sysdeps/ieee754/dbl-64/e_gamma_r.c | 39 +++- sysdeps/ieee754/dbl-64/gamma_product.c | 4 + sysdeps/ieee754/ldbl-96/gamma_product.c | 5 + sysdeps/x86_64/fpu/libm-test-ulps | 6 +- 7 files changed, 395 insertions(+), 15 deletions(-) diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in index 4edaaa8ee1..af4af15cf4 100644 --- a/math/auto-libm-test-in +++ b/math/auto-libm-test-in @@ -8245,8 +8245,12 @@ tgamma -0x6.ce9158p+0 tgamma -0xd.cbf53d0e7d06p+0 # the next value generates larger error bounds on x86_64 (binary32) tgamma -0x3.0aa534p+0 -# the next value generates larger error bounds on x86_64 (binary64) +# the next values generate large error bounds on x86_64 (binary64) tgamma -0x1.62b8c36c7180bp+4 +tgamma -0x1.62c4d519e8677p+3 +tgamma -0x1.c033cc426752fp+2 +tgamma -0x1.62cfd0d34ade2p+3 +tgamma -0x1.8814da6eb7dbp+5 y0 0.125 y0 0.75 diff --git a/math/auto-libm-test-out-tgamma b/math/auto-libm-test-out-tgamma index 032ad3a7d7..20a9916afa 100644 --- a/math/auto-libm-test-out-tgamma +++ b/math/auto-libm-test-out-tgamma @@ -27397,3 +27397,279 @@ tgamma -0x1.62b8c36c7180bp+4 = tgamma tonearest ibm128 -0x1.62b8c36c7180bp+4 : -0xf.3fb73d327d6a801abb3039b104p-72 : inexact-ok = tgamma towardzero ibm128 -0x1.62b8c36c7180bp+4 : -0xf.3fb73d327d6a801abb3039b1p-72 : inexact-ok = tgamma upward ibm128 -0x1.62b8c36c7180bp+4 : -0xf.3fb73d327d6a801abb3039b1p-72 : inexact-ok +tgamma -0x1.62c4d519e8677p+3 += tgamma downward binary32 -0xb.1626ap+0 : 0x3.fac968p-24 : inexact-ok += tgamma tonearest binary32 -0xb.1626ap+0 : 0x3.fac96cp-24 : inexact-ok += tgamma towardzero binary32 -0xb.1626ap+0 : 0x3.fac968p-24 : inexact-ok += tgamma upward binary32 -0xb.1626ap+0 : 0x3.fac96cp-24 : inexact-ok += tgamma downward binary64 -0xb.1626ap+0 : 0x3.fac96a100e60ep-24 : inexact-ok += tgamma tonearest binary64 -0xb.1626ap+0 : 0x3.fac96a100e61p-24 : inexact-ok += tgamma towardzero binary64 -0xb.1626ap+0 : 0x3.fac96a100e60ep-24 : inexact-ok += tgamma upward binary64 -0xb.1626ap+0 : 0x3.fac96a100e61p-24 : inexact-ok += tgamma downward intel96 -0xb.1626ap+0 : 0x3.fac96a100e60f1f8p-24 : inexact-ok += tgamma tonearest intel96 -0xb.1626ap+0 : 0x3.fac96a100e60f1fcp-24 : inexact-ok += tgamma towardzero intel96 -0xb.1626ap+0 : 0x3.fac96a100e60f1f8p-24 : inexact-ok += tgamma upward intel96 -0xb.1626ap+0 : 0x3.fac96a100e60f1fcp-24 : inexact-ok += tgamma downward m68k96 -0xb.1626ap+0 : 0x3.fac96a100e60f1f8p-24 : inexact-ok += tgamma tonearest m68k96 -0xb.1626ap+0 : 0x3.fac96a100e60f1fcp-24 : inexact-ok += tgamma towardzero m68k96 -0xb.1626ap+0 : 0x3.fac96a100e60f1f8p-24 : inexact-ok += tgamma upward m68k96 -0xb.1626ap+0 : 0x3.fac96a100e60f1fcp-24 : inexact-ok += tgamma downward binary128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcdfb6p-24 : inexact-ok += tgamma tonearest binary128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcdfb8p-24 : inexact-ok += tgamma towardzero binary128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcdfb6p-24 : inexact-ok += tgamma upward binary128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcdfb8p-24 : inexact-ok += tgamma downward ibm128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcdfp-24 : inexact-ok += tgamma tonearest ibm128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcep-24 : inexact-ok += tgamma towardzero ibm128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcdfp-24 : inexact-ok += tgamma upward ibm128 -0xb.1626ap+0 : 0x3.fac96a100e60f1fa454b8cdcep-24 : inexact-ok += tgamma downward binary32 -0xb.1626bp+0 : 0x3.fac6p-24 : inexact-ok += tgamma tonearest binary32 -0xb.1626bp+0 : 0x3.fac6p-24 : inexact-ok += tgamma towardzero binary32 -0xb.1626bp+0 : 0x3.fac6p-24 : inexact-ok += tgamma upward binary32 -0xb.1626bp+0 : 0x3.fac604p-24 : inexact-ok += tgamma downward binary64 -0xb.1626bp+0 : 0x3.fac600633d74p-24 : inexact-ok += tgamma tonearest binary64 -0xb.1626bp+0 : 0x3.fac600633d742p-24 : inexact-ok += tgamma towardzero binary64 -0xb.1626bp+0 : 0x3.fac600633d74p-24 : inexact-ok += tgamma upward binary64 -0xb.1626bp+0 : 0x3.fac600633d742p-24 : inexact-ok += tgamma downward intel96 -0xb.1626bp+0 : 0x3.fac600633d7414fp-24 : inexact-ok += tgamma tonearest intel96 -0xb.1626bp+0 : 0x3.fac600633d7414f4p-24 : inexact-ok += tgamma towardzero intel96 -0xb.1626bp+0 : 0x3.fac600633d7414fp-24 : inexact-ok += tgamma upward intel96 -0xb.1626bp+0 : 0x3.fac600633d7414f4p-24 : inexact-ok += tgamma downward m68k96 -0xb.1626bp+0 : 0x3.fac600633d7414fp-24 : inexact-ok += tgamma tonearest m68k96 -0xb.1626bp+0 : 0x3.fac600633d7414f4p-24 : inexact-ok += tgamma towardzero m68k96 -0xb.1626bp+0 : 0x3.fac600633d7414fp-24 : inexact-ok += tgamma upward m68k96 -0xb.1626bp+0 : 0x3.fac600633d7414f4p-24 : inexact-ok += tgamma downward binary128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d97364607cp-24 : inexact-ok += tgamma tonearest binary128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d97364607cp-24 : inexact-ok += tgamma towardzero binary128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d97364607cp-24 : inexact-ok += tgamma upward binary128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d97364607ep-24 : inexact-ok += tgamma downward ibm128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d973646p-24 : inexact-ok += tgamma tonearest ibm128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d973646p-24 : inexact-ok += tgamma towardzero ibm128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d973646p-24 : inexact-ok += tgamma upward ibm128 -0xb.1626bp+0 : 0x3.fac600633d7414f378d9736461p-24 : inexact-ok += tgamma downward binary64 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7ep-24 : inexact-ok += tgamma tonearest binary64 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7ep-24 : inexact-ok += tgamma towardzero binary64 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7ep-24 : inexact-ok += tgamma upward binary64 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d8p-24 : inexact-ok += tgamma downward intel96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e69p-24 : inexact-ok += tgamma tonearest intel96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e694p-24 : inexact-ok += tgamma towardzero intel96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e69p-24 : inexact-ok += tgamma upward intel96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e694p-24 : inexact-ok += tgamma downward m68k96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e69p-24 : inexact-ok += tgamma tonearest m68k96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e694p-24 : inexact-ok += tgamma towardzero m68k96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e69p-24 : inexact-ok += tgamma upward m68k96 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e694p-24 : inexact-ok += tgamma downward binary128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b139429ecp-24 : inexact-ok += tgamma tonearest binary128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b139429ecp-24 : inexact-ok += tgamma towardzero binary128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b139429ecp-24 : inexact-ok += tgamma upward binary128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b139429eep-24 : inexact-ok += tgamma downward ibm128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b139429p-24 : inexact-ok += tgamma tonearest ibm128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b13942ap-24 : inexact-ok += tgamma towardzero ibm128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b139429p-24 : inexact-ok += tgamma upward ibm128 -0xb.1626a8cf433b8p+0 : 0x3.fac7890382d7e693c33b13942ap-24 : inexact-ok +tgamma -0x1.c033cc426752fp+2 += tgamma downward binary32 -0x7.00cf3p+0 : 0xf.f7014p-8 : inexact-ok += tgamma tonearest binary32 -0x7.00cf3p+0 : 0xf.f7014p-8 : inexact-ok += tgamma towardzero binary32 -0x7.00cf3p+0 : 0xf.f7014p-8 : inexact-ok += tgamma upward binary32 -0x7.00cf3p+0 : 0xf.f7015p-8 : inexact-ok += tgamma downward binary64 -0x7.00cf3p+0 : 0xf.f7014642d9848p-8 : inexact-ok += tgamma tonearest binary64 -0x7.00cf3p+0 : 0xf.f7014642d9848p-8 : inexact-ok += tgamma towardzero binary64 -0x7.00cf3p+0 : 0xf.f7014642d9848p-8 : inexact-ok += tgamma upward binary64 -0x7.00cf3p+0 : 0xf.f7014642d985p-8 : inexact-ok += tgamma downward intel96 -0x7.00cf3p+0 : 0xf.f7014642d9849fcp-8 : inexact-ok += tgamma tonearest intel96 -0x7.00cf3p+0 : 0xf.f7014642d9849fcp-8 : inexact-ok += tgamma towardzero intel96 -0x7.00cf3p+0 : 0xf.f7014642d9849fcp-8 : inexact-ok += tgamma upward intel96 -0x7.00cf3p+0 : 0xf.f7014642d9849fdp-8 : inexact-ok += tgamma downward m68k96 -0x7.00cf3p+0 : 0xf.f7014642d9849fcp-8 : inexact-ok += tgamma tonearest m68k96 -0x7.00cf3p+0 : 0xf.f7014642d9849fcp-8 : inexact-ok += tgamma towardzero m68k96 -0x7.00cf3p+0 : 0xf.f7014642d9849fcp-8 : inexact-ok += tgamma upward m68k96 -0x7.00cf3p+0 : 0xf.f7014642d9849fdp-8 : inexact-ok += tgamma downward binary128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea2ae8p-8 : inexact-ok += tgamma tonearest binary128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea2afp-8 : inexact-ok += tgamma towardzero binary128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea2ae8p-8 : inexact-ok += tgamma upward binary128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea2afp-8 : inexact-ok += tgamma downward ibm128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea28p-8 : inexact-ok += tgamma tonearest ibm128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea2cp-8 : inexact-ok += tgamma towardzero ibm128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea28p-8 : inexact-ok += tgamma upward ibm128 -0x7.00cf3p+0 : 0xf.f7014642d9849fc3f8a3bdea2cp-8 : inexact-ok += tgamma downward binary32 -0x7.00cf38p+0 : 0xf.f6627p-8 : inexact-ok += tgamma tonearest binary32 -0x7.00cf38p+0 : 0xf.f6628p-8 : inexact-ok += tgamma towardzero binary32 -0x7.00cf38p+0 : 0xf.f6627p-8 : inexact-ok += tgamma upward binary32 -0x7.00cf38p+0 : 0xf.f6628p-8 : inexact-ok += tgamma downward binary64 -0x7.00cf38p+0 : 0xf.f6627d248f8ep-8 : inexact-ok += tgamma tonearest binary64 -0x7.00cf38p+0 : 0xf.f6627d248f8ep-8 : inexact-ok += tgamma towardzero binary64 -0x7.00cf38p+0 : 0xf.f6627d248f8ep-8 : inexact-ok += tgamma upward binary64 -0x7.00cf38p+0 : 0xf.f6627d248f8e8p-8 : inexact-ok += tgamma downward intel96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25dp-8 : inexact-ok += tgamma tonearest intel96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ep-8 : inexact-ok += tgamma towardzero intel96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25dp-8 : inexact-ok += tgamma upward intel96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ep-8 : inexact-ok += tgamma downward m68k96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25dp-8 : inexact-ok += tgamma tonearest m68k96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ep-8 : inexact-ok += tgamma towardzero m68k96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25dp-8 : inexact-ok += tgamma upward m68k96 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ep-8 : inexact-ok += tgamma downward binary128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81585fcbp-8 : inexact-ok += tgamma tonearest binary128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81585fcb8p-8 : inexact-ok += tgamma towardzero binary128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81585fcbp-8 : inexact-ok += tgamma upward binary128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81585fcb8p-8 : inexact-ok += tgamma downward ibm128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81585fcp-8 : inexact-ok += tgamma tonearest ibm128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81585fcp-8 : inexact-ok += tgamma towardzero ibm128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81585fcp-8 : inexact-ok += tgamma upward ibm128 -0x7.00cf38p+0 : 0xf.f6627d248f8e25ddceb81586p-8 : inexact-ok += tgamma downward binary64 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad9999808p-8 : inexact-ok += tgamma tonearest binary64 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999981p-8 : inexact-ok += tgamma towardzero binary64 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad9999808p-8 : inexact-ok += tgamma upward binary64 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999981p-8 : inexact-ok += tgamma downward intel96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4ap-8 : inexact-ok += tgamma tonearest intel96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4ap-8 : inexact-ok += tgamma towardzero intel96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4ap-8 : inexact-ok += tgamma upward intel96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4bp-8 : inexact-ok += tgamma downward m68k96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4ap-8 : inexact-ok += tgamma tonearest m68k96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4ap-8 : inexact-ok += tgamma towardzero m68k96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4ap-8 : inexact-ok += tgamma upward m68k96 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4bp-8 : inexact-ok += tgamma downward binary128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608814ecp-8 : inexact-ok += tgamma tonearest binary128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608814ecp-8 : inexact-ok += tgamma towardzero binary128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608814ecp-8 : inexact-ok += tgamma upward binary128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608814ec8p-8 : inexact-ok += tgamma downward ibm128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608814cp-8 : inexact-ok += tgamma tonearest ibm128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608815p-8 : inexact-ok += tgamma towardzero ibm128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608814cp-8 : inexact-ok += tgamma upward ibm128 -0x7.00cf31099d4bcp+0 : 0xf.f6ecad999980e4a4868608815p-8 : inexact-ok +tgamma -0x1.62cfd0d34ade2p+3 += tgamma downward binary32 -0xb.167e8p+0 : 0x3.e85604p-24 : inexact-ok += tgamma tonearest binary32 -0xb.167e8p+0 : 0x3.e85608p-24 : inexact-ok += tgamma towardzero binary32 -0xb.167e8p+0 : 0x3.e85604p-24 : inexact-ok += tgamma upward binary32 -0xb.167e8p+0 : 0x3.e85608p-24 : inexact-ok += tgamma downward binary64 -0xb.167e8p+0 : 0x3.e85606c693e9ap-24 : inexact-ok += tgamma tonearest binary64 -0xb.167e8p+0 : 0x3.e85606c693e9cp-24 : inexact-ok += tgamma towardzero binary64 -0xb.167e8p+0 : 0x3.e85606c693e9ap-24 : inexact-ok += tgamma upward binary64 -0xb.167e8p+0 : 0x3.e85606c693e9cp-24 : inexact-ok += tgamma downward intel96 -0xb.167e8p+0 : 0x3.e85606c693e9bf24p-24 : inexact-ok += tgamma tonearest intel96 -0xb.167e8p+0 : 0x3.e85606c693e9bf24p-24 : inexact-ok += tgamma towardzero intel96 -0xb.167e8p+0 : 0x3.e85606c693e9bf24p-24 : inexact-ok += tgamma upward intel96 -0xb.167e8p+0 : 0x3.e85606c693e9bf28p-24 : inexact-ok += tgamma downward m68k96 -0xb.167e8p+0 : 0x3.e85606c693e9bf24p-24 : inexact-ok += tgamma tonearest m68k96 -0xb.167e8p+0 : 0x3.e85606c693e9bf24p-24 : inexact-ok += tgamma towardzero m68k96 -0xb.167e8p+0 : 0x3.e85606c693e9bf24p-24 : inexact-ok += tgamma upward m68k96 -0xb.167e8p+0 : 0x3.e85606c693e9bf28p-24 : inexact-ok += tgamma downward binary128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8c3ap-24 : inexact-ok += tgamma tonearest binary128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8c3cp-24 : inexact-ok += tgamma towardzero binary128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8c3ap-24 : inexact-ok += tgamma upward binary128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8c3cp-24 : inexact-ok += tgamma downward ibm128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8cp-24 : inexact-ok += tgamma tonearest ibm128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8cp-24 : inexact-ok += tgamma towardzero ibm128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8cp-24 : inexact-ok += tgamma upward ibm128 -0xb.167e8p+0 : 0x3.e85606c693e9bf25ff07132c8dp-24 : inexact-ok += tgamma downward binary32 -0xb.167e9p+0 : 0x3.e852b8p-24 : inexact-ok += tgamma tonearest binary32 -0xb.167e9p+0 : 0x3.e852b8p-24 : inexact-ok += tgamma towardzero binary32 -0xb.167e9p+0 : 0x3.e852b8p-24 : inexact-ok += tgamma upward binary32 -0xb.167e9p+0 : 0x3.e852bcp-24 : inexact-ok += tgamma downward binary64 -0xb.167e9p+0 : 0x3.e852b838ddf22p-24 : inexact-ok += tgamma tonearest binary64 -0xb.167e9p+0 : 0x3.e852b838ddf24p-24 : inexact-ok += tgamma towardzero binary64 -0xb.167e9p+0 : 0x3.e852b838ddf22p-24 : inexact-ok += tgamma upward binary64 -0xb.167e9p+0 : 0x3.e852b838ddf24p-24 : inexact-ok += tgamma downward intel96 -0xb.167e9p+0 : 0x3.e852b838ddf2348cp-24 : inexact-ok += tgamma tonearest intel96 -0xb.167e9p+0 : 0x3.e852b838ddf2349p-24 : inexact-ok += tgamma towardzero intel96 -0xb.167e9p+0 : 0x3.e852b838ddf2348cp-24 : inexact-ok += tgamma upward intel96 -0xb.167e9p+0 : 0x3.e852b838ddf2349p-24 : inexact-ok += tgamma downward m68k96 -0xb.167e9p+0 : 0x3.e852b838ddf2348cp-24 : inexact-ok += tgamma tonearest m68k96 -0xb.167e9p+0 : 0x3.e852b838ddf2349p-24 : inexact-ok += tgamma towardzero m68k96 -0xb.167e9p+0 : 0x3.e852b838ddf2348cp-24 : inexact-ok += tgamma upward m68k96 -0xb.167e9p+0 : 0x3.e852b838ddf2349p-24 : inexact-ok += tgamma downward binary128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd745531ep-24 : inexact-ok += tgamma tonearest binary128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd745531ep-24 : inexact-ok += tgamma towardzero binary128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd745531ep-24 : inexact-ok += tgamma upward binary128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd745532p-24 : inexact-ok += tgamma downward ibm128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd74553p-24 : inexact-ok += tgamma tonearest ibm128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd74553p-24 : inexact-ok += tgamma towardzero ibm128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd74553p-24 : inexact-ok += tgamma upward ibm128 -0xb.167e9p+0 : 0x3.e852b838ddf2348f185dd74554p-24 : inexact-ok += tgamma downward binary64 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa6p-24 : inexact-ok += tgamma tonearest binary64 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa8p-24 : inexact-ok += tgamma towardzero binary64 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa6p-24 : inexact-ok += tgamma upward binary64 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa8p-24 : inexact-ok += tgamma downward intel96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84p-24 : inexact-ok += tgamma tonearest intel96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84p-24 : inexact-ok += tgamma towardzero intel96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84p-24 : inexact-ok += tgamma upward intel96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e88p-24 : inexact-ok += tgamma downward m68k96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84p-24 : inexact-ok += tgamma tonearest m68k96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84p-24 : inexact-ok += tgamma towardzero m68k96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84p-24 : inexact-ok += tgamma upward m68k96 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e88p-24 : inexact-ok += tgamma downward binary128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e7829836ap-24 : inexact-ok += tgamma tonearest binary128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e7829836ap-24 : inexact-ok += tgamma towardzero binary128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e7829836ap-24 : inexact-ok += tgamma upward binary128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e7829836cp-24 : inexact-ok += tgamma downward ibm128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e782983p-24 : inexact-ok += tgamma tonearest ibm128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e782983p-24 : inexact-ok += tgamma towardzero ibm128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e782983p-24 : inexact-ok += tgamma upward ibm128 -0xb.167e869a56f1p+0 : 0x3.e854a96acdfa7e84ba5e782984p-24 : inexact-ok +tgamma -0x1.8814da6eb7dbp+5 += tgamma downward binary32 -0x3.1029b4p+4 : 0x0p+0 : inexact-ok underflow errno-erange += tgamma tonearest binary32 -0x3.1029b4p+4 : 0x0p+0 : inexact-ok underflow errno-erange += tgamma towardzero binary32 -0x3.1029b4p+4 : 0x0p+0 : inexact-ok underflow errno-erange += tgamma upward binary32 -0x3.1029b4p+4 : 0x8p-152 : inexact-ok underflow errno-erange-ok += tgamma downward binary64 -0x3.1029b4p+4 : 0x3.fd9063a22653ep-204 : inexact-ok += tgamma tonearest binary64 -0x3.1029b4p+4 : 0x3.fd9063a22653ep-204 : inexact-ok += tgamma towardzero binary64 -0x3.1029b4p+4 : 0x3.fd9063a22653ep-204 : inexact-ok += tgamma upward binary64 -0x3.1029b4p+4 : 0x3.fd9063a22654p-204 : inexact-ok += tgamma downward intel96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e4p-204 : inexact-ok += tgamma tonearest intel96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e4p-204 : inexact-ok += tgamma towardzero intel96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e4p-204 : inexact-ok += tgamma upward intel96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e8p-204 : inexact-ok += tgamma downward m68k96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e4p-204 : inexact-ok += tgamma tonearest m68k96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e4p-204 : inexact-ok += tgamma towardzero m68k96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e4p-204 : inexact-ok += tgamma upward m68k96 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e8p-204 : inexact-ok += tgamma downward binary128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a5bcp-204 : inexact-ok += tgamma tonearest binary128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a5bcp-204 : inexact-ok += tgamma towardzero binary128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a5bcp-204 : inexact-ok += tgamma upward binary128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a5bep-204 : inexact-ok += tgamma downward ibm128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a5p-204 : inexact-ok += tgamma tonearest ibm128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a6p-204 : inexact-ok += tgamma towardzero ibm128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a5p-204 : inexact-ok += tgamma upward ibm128 -0x3.1029b4p+4 : 0x3.fd9063a22653e8e5a73bb2d5a6p-204 : inexact-ok += tgamma downward binary32 -0x3.1029b8p+4 : 0x0p+0 : inexact-ok underflow errno-erange += tgamma tonearest binary32 -0x3.1029b8p+4 : 0x0p+0 : inexact-ok underflow errno-erange += tgamma towardzero binary32 -0x3.1029b8p+4 : 0x0p+0 : inexact-ok underflow errno-erange += tgamma upward binary32 -0x3.1029b8p+4 : 0x8p-152 : inexact-ok underflow errno-erange-ok += tgamma downward binary64 -0x3.1029b8p+4 : 0x3.fd2a955e20f56p-204 : inexact-ok += tgamma tonearest binary64 -0x3.1029b8p+4 : 0x3.fd2a955e20f58p-204 : inexact-ok += tgamma towardzero binary64 -0x3.1029b8p+4 : 0x3.fd2a955e20f56p-204 : inexact-ok += tgamma upward binary64 -0x3.1029b8p+4 : 0x3.fd2a955e20f58p-204 : inexact-ok += tgamma downward intel96 -0x3.1029b8p+4 : 0x3.fd2a955e20f57628p-204 : inexact-ok += tgamma tonearest intel96 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762cp-204 : inexact-ok += tgamma towardzero intel96 -0x3.1029b8p+4 : 0x3.fd2a955e20f57628p-204 : inexact-ok += tgamma upward intel96 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762cp-204 : inexact-ok += tgamma downward m68k96 -0x3.1029b8p+4 : 0x3.fd2a955e20f57628p-204 : inexact-ok += tgamma tonearest m68k96 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762cp-204 : inexact-ok += tgamma towardzero m68k96 -0x3.1029b8p+4 : 0x3.fd2a955e20f57628p-204 : inexact-ok += tgamma upward m68k96 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762cp-204 : inexact-ok += tgamma downward binary128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f01716p-204 : inexact-ok += tgamma tonearest binary128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f01718p-204 : inexact-ok += tgamma towardzero binary128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f01716p-204 : inexact-ok += tgamma upward binary128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f01718p-204 : inexact-ok += tgamma downward ibm128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f017p-204 : inexact-ok += tgamma tonearest ibm128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f017p-204 : inexact-ok += tgamma towardzero ibm128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f017p-204 : inexact-ok += tgamma upward ibm128 -0x3.1029b8p+4 : 0x3.fd2a955e20f5762a464af7f018p-204 : inexact-ok += tgamma downward binary64 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5cp-204 : inexact-ok += tgamma tonearest binary64 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5cp-204 : inexact-ok += tgamma towardzero binary64 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5cp-204 : inexact-ok += tgamma upward binary64 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5ep-204 : inexact-ok += tgamma downward intel96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d4p-204 : inexact-ok += tgamma tonearest intel96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d8p-204 : inexact-ok += tgamma towardzero intel96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d4p-204 : inexact-ok += tgamma upward intel96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d8p-204 : inexact-ok += tgamma downward m68k96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d4p-204 : inexact-ok += tgamma tonearest m68k96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d8p-204 : inexact-ok += tgamma towardzero m68k96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d4p-204 : inexact-ok += tgamma upward m68k96 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d8p-204 : inexact-ok += tgamma downward binary128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d82ep-204 : inexact-ok += tgamma tonearest binary128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d82ep-204 : inexact-ok += tgamma towardzero binary128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d82ep-204 : inexact-ok += tgamma upward binary128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d83p-204 : inexact-ok += tgamma downward ibm128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d8p-204 : inexact-ok += tgamma tonearest ibm128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d8p-204 : inexact-ok += tgamma towardzero ibm128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d8p-204 : inexact-ok += tgamma upward ibm128 -0x3.1029b4dd6fb6p+4 : 0x3.fd7a5e1abda5c6d62103dbf6d9p-204 : inexact-ok diff --git a/math/mul_split.h b/math/mul_split.h index 8cc60ec630..24138c8c39 100644 --- a/math/mul_split.h +++ b/math/mul_split.h @@ -47,4 +47,78 @@ mul_split (double *hi, double *lo, double x, double y) #endif } +/* Define WANT_MUL_EXPANSION to use mul_expansion(). + Define WANT_DIV_EXPANSION to use div_expansion(). */ + +#if defined(WANT_MUL_EXPANSION) || defined(WANT_DIV_EXPANSION) +/* Add a + b exactly, such that *hi + *lo = a + b. + Assumes |a| >= |b| and rounding to nearest. */ +static void +fast_two_sum (double *hi, double *lo, double a, double b) +{ + double e; + + *hi = a + b; + e = *hi - a; /* exact */ + *lo = b - e; /* exact */ + /* Now *hi + *lo = a + b exactly. */ +} +#endif + +#ifdef WANT_MUL_EXPANSION +/* Multiplication of two floating-point expansions: *hi + *lo is an + approximation of (h1+l1)*(h2+l2), assuming |l1| <= 1/2*ulp(h1) + and |l2| <= 1/2*ulp(h2) and rounding to nearest. */ +static void +mul_expansion (double *hi, double *lo, double h1, double l1, + double h2, double l2) +{ + double r, e; + + mul_split (hi, lo, h1, h2); + r = h1 * l2 + h2 * l1; + /* Now add r to (hi,lo) using fast two-sum, where we know |r| < |hi|. */ + fast_two_sum (hi, &e, *hi, r); + *lo -= e; +} +#endif + +#ifdef WANT_DIV_EXPANSION +/* Calculate X / Y and store the approximate result in *HI + *LO. It is + assumed that Y is not zero, that no overflow nor underflow occurs, and + rounding is to nearest. */ +static void +div_split (double *hi, double *lo, double x, double y) +{ + double a, b; + + *hi = x / y; + mul_split (&a, &b, *hi, y); + /* a + b = hi*y, which should be near x. */ + a = x - a; /* huge cancellation */ + a = a - b; + /* Now x ~ hi*y + a thus x/y ~ hi + a/y. */ + *lo = a / y; +} + +/* Division of two floating-point expansions: *hi + *lo is an + approximation of (h1+l1)/(h2+l2), assuming |l1| <= 1/2*ulp(h1) + and |l2| <= 1/2*ulp(h2), h2+l2 is not zero, and rounding to nearest. */ +static void +div_expansion (double *hi, double *lo, double h1, double l1, + double h2, double l2) +{ + double r, e; + + div_split (hi, lo, h1, h2); + /* (h1+l1)/(h2+l2) ~ h1/h2 + (l1*h2 - l2*h1)/h2^2 */ + r = (l1 * h2 - l2 * h1) / (h2 * h2); + /* Now add r to (hi,lo) using fast two-sum, where we know |r| < |hi|. */ + fast_two_sum (hi, &e, *hi, r); + *lo += e; + /* Renormalize since |lo| might be larger than 0.5 ulp(hi). */ + fast_two_sum (hi, lo, *hi, *lo); +} +#endif + #endif /* _MUL_SPLIT_H */ diff --git a/sysdeps/ieee754/dbl-64/e_gamma_r.c b/sysdeps/ieee754/dbl-64/e_gamma_r.c index fe69920c76..7c8dc800cf 100644 --- a/sysdeps/ieee754/dbl-64/e_gamma_r.c +++ b/sysdeps/ieee754/dbl-64/e_gamma_r.c @@ -24,6 +24,9 @@ #include #include #include +#define WANT_MUL_EXPANSION +#define WANT_DIV_EXPANSION +#include /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's approximation to gamma function. */ @@ -88,7 +91,6 @@ gamma_positive (double x, int *exp2_adj) Compute gamma (X_ADJ + X_EPS) using Stirling's approximation, starting by computing pow (X_ADJ, X_ADJ) with a power of 2 factored out. */ - double exp_adj = -eps; double x_adj_int = round (x_adj); double x_adj_frac = x_adj - x_adj_int; int x_adj_log2; @@ -99,18 +101,22 @@ gamma_positive (double x, int *exp2_adj) x_adj_mant *= 2.0; } *exp2_adj = x_adj_log2 * (int) x_adj_int; - double ret = (__ieee754_pow (x_adj_mant, x_adj) - * __ieee754_exp2 (x_adj_log2 * x_adj_frac) - * __ieee754_exp (-x_adj) - * sqrt (2 * M_PI / x_adj) - / prod); - exp_adj += x_eps * __ieee754_log (x_adj); + double h1, l1, h2, l2; + mul_split (&h1, &l1, __ieee754_pow (x_adj_mant, x_adj), + __ieee754_exp2 (x_adj_log2 * x_adj_frac)); + mul_split (&h2, &l2, __ieee754_exp (-x_adj), sqrt (2 * M_PI / x_adj)); + mul_expansion (&h1, &l1, h1, l1, h2, l2); + /* Divide by prod * (1 + eps). */ + div_expansion (&h1, &l1, h1, l1, prod, prod * eps); + double exp_adj = x_eps * __ieee754_log (x_adj); double bsum = gamma_coeff[NCOEFF - 1]; double x_adj2 = x_adj * x_adj; for (size_t i = 1; i <= NCOEFF - 1; i++) bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i]; exp_adj += bsum / x_adj; - return ret + ret * __expm1 (exp_adj); + /* Now return (h1+l1) * exp(exp_adj), where exp_adj is small. */ + l1 += h1 * __expm1 (exp_adj); + return h1 + l1; } } @@ -188,9 +194,20 @@ __ieee754_gamma_r (double x, int *signgamp) ? __sin (M_PI * frac) : __cos (M_PI * (0.5 - frac))); int exp2_adj; - double tret = M_PI / (-x * sinpix - * gamma_positive (-x, &exp2_adj)); - ret = __scalbn (tret, -exp2_adj); + double h1, l1, h2, l2; + h2 = gamma_positive (-x, &exp2_adj); + mul_split (&h1, &l1, sinpix, h2); + /* sinpix*gamma_positive(.) = h1 + l1 */ + mul_split (&h2, &l2, h1, x); + /* h1*x = h2 + l2 */ + /* (h1 + l1) * x = h1*x + l1*x = h2 + l2 + l1*x */ + l2 += l1 * x; + /* x*sinpix*gamma_positive(.) ~ h2 + l2 */ + h1 = 0x3.243f6a8885a3p+0; /* binary64 approximation of Pi */ + l1 = 0x8.d313198a2e038p-56; /* |h1+l1-Pi| < 3e-33 */ + /* Now we divide h1 + l1 by h2 + l2. */ + div_expansion (&h1, &l1, h1, l1, h2, l2); + ret = __scalbn (-h1, -exp2_adj); math_check_force_underflow_nonneg (ret); } } diff --git a/sysdeps/ieee754/dbl-64/gamma_product.c b/sysdeps/ieee754/dbl-64/gamma_product.c index 643df11b53..2e3b5148dc 100644 --- a/sysdeps/ieee754/dbl-64/gamma_product.c +++ b/sysdeps/ieee754/dbl-64/gamma_product.c @@ -16,6 +16,10 @@ License along with the GNU C Library; if not, see . */ +/* Note: this function is only used on an architecture that doesn't use + ldbl-96 (i.e., anything other than x86_64, i386, ia64 or m68k). + On other architectures the one in ldbl-96 is used. */ + #include #include #include diff --git a/sysdeps/ieee754/ldbl-96/gamma_product.c b/sysdeps/ieee754/ldbl-96/gamma_product.c index 9e00e71452..f881649801 100644 --- a/sysdeps/ieee754/ldbl-96/gamma_product.c +++ b/sysdeps/ieee754/ldbl-96/gamma_product.c @@ -16,6 +16,11 @@ License along with the GNU C Library; if not, see . */ +/* Note: this function is used on an architecture that uses + ldbl-96 (i.e., x86_64, i386, ia64 or m68k), and in particular + it is called from dbl-64/e_gamma_r.c. + On other architectures the one in dbl-64 is used. */ + #include #include #include diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps index bd1fa63702..e0a362a347 100644 --- a/sysdeps/x86_64/fpu/libm-test-ulps +++ b/sysdeps/x86_64/fpu/libm-test-ulps @@ -1734,16 +1734,16 @@ float128: 4 ldouble: 5 Function: "tgamma_downward": -double: 8 +double: 9 float: 7 float128: 5 -ldouble: 5 +ldouble: 6 Function: "tgamma_towardzero": double: 9 float: 7 float128: 5 -ldouble: 5 +ldouble: 6 Function: "tgamma_upward": double: 9