From patchwork Fri Dec 11 12:11:50 2015 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Siddhesh Poyarekar X-Patchwork-Id: 58275 Delivered-To: patch@linaro.org Received: by 10.112.157.166 with SMTP id wn6csp90426lbb; Fri, 11 Dec 2015 04:13:01 -0800 (PST) X-Received: by 10.98.65.76 with SMTP id o73mr14730871pfa.119.1449835981261; Fri, 11 Dec 2015 04:13:01 -0800 (PST) Return-Path: Received: from sourceware.org (server1.sourceware.org. [209.132.180.131]) by mx.google.com with ESMTPS id 12si960838pfb.199.2015.12.11.04.13.00 for (version=TLS1_2 cipher=ECDHE-RSA-AES128-GCM-SHA256 bits=128/128); Fri, 11 Dec 2015 04:13:01 -0800 (PST) Received-SPF: pass (google.com: domain of libc-alpha-return-65497-patch=linaro.org@sourceware.org designates 209.132.180.131 as permitted sender) client-ip=209.132.180.131; Authentication-Results: mx.google.com; spf=pass (google.com: domain of libc-alpha-return-65497-patch=linaro.org@sourceware.org designates 209.132.180.131 as permitted sender) smtp.mailfrom=libc-alpha-return-65497-patch=linaro.org@sourceware.org; dkim=pass header.i=@sourceware.org DomainKey-Signature: a=rsa-sha1; c=nofws; d=sourceware.org; h=list-id :list-unsubscribe:list-subscribe:list-archive:list-post :list-help:sender:from:to:cc:subject:date:message-id:in-reply-to :references; q=dns; s=default; b=jxtLWqenT2WqMGQiWyCtnKe+ifGCPGR ebEmbTozyaR2iLTBbTgeGa7OziBlRspfVtmzpRH6YxTQ9x3ZsKnNukSOa/QsoBdR AO8qEVMhqGhzxrrl9dr3xtmF/FhAXVQgHGrstq9ngtlVphknfbd3LHwC942xPThw q5PQuh68d6DQ= DKIM-Signature: v=1; a=rsa-sha1; c=relaxed; d=sourceware.org; h=list-id :list-unsubscribe:list-subscribe:list-archive:list-post :list-help:sender:from:to:cc:subject:date:message-id:in-reply-to :references; s=default; bh=ZX78VtdU5DkWu5WiV21SOBIVULk=; b=pdTkZ KWOi4kKwJNRiXEkrMAcSgnVWZW4AXr0INCQZiPCgeswiUpBpIWpU802ZyqrCKXYN ROop96IbW1zukR6pJ2QS1QbFqE+p7gfNJwRvdgHv8XgSj3nW0F/FX1n3Ef97Zums 9DXc7C8vmyKsyJFOwhC+SBUy7yQTeh5iBhgkyI= Received: (qmail 57417 invoked by alias); 11 Dec 2015 12:12:43 -0000 Mailing-List: contact libc-alpha-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Unsubscribe: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: libc-alpha-owner@sourceware.org Delivered-To: mailing list libc-alpha@sourceware.org Received: (qmail 57113 invoked by uid 89); 11 Dec 2015 12:12:42 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-2.1 required=5.0 tests=AWL, BAYES_00, KAM_ASCII_DIVIDERS, RCVD_IN_DNSWL_LOW, SPF_PASS autolearn=no version=3.3.2 X-HELO: mail-pa0-f44.google.com X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20130820; h=x-gm-message-state:from:to:cc:subject:date:message-id:in-reply-to :references; bh=GJCeSYnoNsOgKK5OAn0QbjNCFF7RPrmY3HAzRAI4nZ4=; b=h6Qlwz8OtVlNiXD4vgB/tq3oPoXwRiqwxcQ8yOAZzIqoFfNlnD0UbSxaUkktebb72c 7NxS6yKE+ROqxM0yBDcXD3t/ENEArPk6JLFmmEYTuCyRyd253Hb35q+j0p2tPc9t3rmW uHkPAiWwS3/0ePjKqyPhmOYVuPjh+zEvPM1h5YfLKbrwFJcj3i5+A13aYDnOMpOemgPt RDjVT6yjV/2Qw8q5/+wc0ik6zk758g6NNgGcNmtfXROAt1RaQLN806CskBAqrKfLdy3C NTS4RWIAfFIHDKEQYbg3nY8ljqFcdaGEZ2AZYZUdzAhUU90l2WsRW+WpWiVmTeT9PM5f Q4fQ== X-Gm-Message-State: ALoCoQmaw7IFupM6boMp2YKiGevgP5VpJbOnmzgsE6HpjFa5mmgWz4GhW3D/LRPBAb4lC/LdD0H2/MtbaA5Q7Wp5Cz4rlvln0A== X-Received: by 10.66.102.4 with SMTP id fk4mr24733527pab.85.1449835954347; Fri, 11 Dec 2015 04:12:34 -0800 (PST) From: Siddhesh Poyarekar To: libc-alpha@sourceware.org Cc: Siddhesh Poyarekar Subject: [PATCH 3/3] Consolidate sincos computation for 2.426265 < |x| < 105414350 Date: Fri, 11 Dec 2015 17:41:50 +0530 Message-Id: <1449835910-4651-4-git-send-email-siddhesh.poyarekar@linaro.org> In-Reply-To: <1449835910-4651-1-git-send-email-siddhesh.poyarekar@linaro.org> References: <1449835910-4651-1-git-send-email-siddhesh.poyarekar@linaro.org> Like the previous change, exploit the fact that computation for sin and cos is identical except that it is apart by a quadrant. Siddhesh * sysdeps/ieee754/dbl-64/s_sin.c [IN_SINCOS]: Avoid defining sloww, sloww1, sloww2 and csloww2 when included in s_sincos.c. * sysdeps/ieee754/dbl-64/s_sincos.c (__sin_local, __cos_local): Move common code... (do_sincos_1): ... into a new function. (sincos_sloww, sincos_sloww1, sincos_sloww2): New functions. (__sincos): Call do_sincos_1. --- sysdeps/ieee754/dbl-64/s_sin.c | 8 +- sysdeps/ieee754/dbl-64/s_sincos.c | 343 ++++++++++++++++++++++---------------- 2 files changed, 208 insertions(+), 143 deletions(-) -- 2.5.0 diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c index a8c5d4a..bc23e84 100644 --- a/sysdeps/ieee754/dbl-64/s_sin.c +++ b/sysdeps/ieee754/dbl-64/s_sin.c @@ -135,9 +135,11 @@ double __mpcos (double x, double dx, bool reduce_range); static double slow (double x); static double slow1 (double x); static double slow2 (double x); +#ifndef IN_SINCOS static double sloww (double x, double dx, double orig); static double sloww1 (double x, double dx, double orig, int m); static double sloww2 (double x, double dx, double orig, int n); +#endif static double bsloww (double x, double dx, double orig, int n); static double bsloww1 (double x, double dx, double orig, int n); static double bsloww2 (double x, double dx, double orig, int n); @@ -145,7 +147,9 @@ int __branred (double x, double *a, double *aa); static double cslow2 (double x); static double csloww (double x, double dx, double orig); static double csloww1 (double x, double dx, double orig, int m); +#ifndef IN_SINCOS static double csloww2 (double x, double dx, double orig, int n); +#endif /* Given a number partitioned into U and X such that U is an index into the sin/cos table, this macro computes the cosine of the number by combining @@ -818,6 +822,7 @@ slow2 (double x) return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); } +#ifndef IN_SINCOS /***************************************************************************/ /* Routine compute sin(x+dx) (Double-Length number) where x is small enough*/ /* to use Taylor series around zero and (x+dx) */ @@ -946,6 +951,7 @@ sloww2 (double x, double dx, double orig, int n) return __mpsin (orig, 0, true); } +#endif /***************************************************************************/ /* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */ @@ -1175,6 +1181,7 @@ csloww1 (double x, double dx, double orig, int m) } +#ifndef IN_SINCOS /***************************************************************************/ /* Routine compute sin(x+dx) (Double-Length number) where x in second or */ /* fourth quarter of unit circle.Routine receive also the original value */ @@ -1207,7 +1214,6 @@ csloww2 (double x, double dx, double orig, int n) return __mpcos (orig, 0, true); } -#ifndef IN_SINCOS # ifndef __cos weak_alias (__cos, cos) # ifdef NO_LONG_DOUBLE diff --git a/sysdeps/ieee754/dbl-64/s_sincos.c b/sysdeps/ieee754/dbl-64/s_sincos.c index 24dc43c..4d72d2b 100644 --- a/sysdeps/ieee754/dbl-64/s_sincos.c +++ b/sysdeps/ieee754/dbl-64/s_sincos.c @@ -33,9 +33,8 @@ static double __always_inline __sin_local (double x, int4 k) { - double xx, res, t, cor, y, s, c, sn, ssn, cs, ccs, xn, a, da, eps; - mynumber u, v; - int4 m, n; + double xx, res, t, cor, y, s, c, sn, ssn, cs, ccs; + mynumber u; double retval = 0; if (k < 0x3e500000) /* if x->0 =>sin(x)=x */ @@ -74,7 +73,7 @@ __sin_local (double x, int4 k) } /* else if (k < 0x3feb6000) */ /*----------------------- 0.855469 <|x|<2.426265 ----------------------*/ - else if (k < 0x400368fd) + else { y = (x > 0) ? hp0 - x : hp0 + x; @@ -92,72 +91,6 @@ __sin_local (double x, int4 k) retval = (res == res + 1.020 * cor) ? ((x > 0) ? res : -res) : slow2 (x); } /* else if (k < 0x400368fd) */ -/*-------------------------- 2.426265<|x|< 105414350 ----------------------*/ - else - { - t = (x * hpinv + toint); - xn = t - toint; - v.x = t; - y = (x - xn * mp1) - xn * mp2; - n = v.i[LOW_HALF] & 3; - da = xn * mp3; - a = y - da; - da = (y - a) - da; - eps = fabs (x) * 1.2e-30; - - switch (n) - { /* quarter of unit circle */ - case 0: - case 2: - xx = a * a; - if (n) - { - a = -a; - da = -da; - } - if (xx < 0.01588) - { - /* Taylor series. */ - res = TAYLOR_SIN (xx, a, da, cor); - cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; - retval = (res == res + cor) ? res : sloww (a, da, x); - } - else - { - if (a > 0) - m = 1; - else - { - m = 0; - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big); - res = do_sin (u, y, da, &cor); - cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; - retval = ((res == res + cor) ? ((m) ? res : -res) - : sloww1 (a, da, x, m)); - } - break; - - case 1: - case 3: - if (a < 0) - { - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big) + da; - res = do_cos (u, y, &cor); - cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; - retval = ((res == res + cor) ? ((n & 2) ? -res : res) - : sloww2 (a, da, x, n)); - break; - } - } /* else if (k < 0x419921FB ) */ - return retval; } @@ -168,9 +101,9 @@ static double __always_inline __cos_local (double x, int4 k) { - double y, xx, res, t, cor, xn, a, da, eps; - mynumber u, v; - int4 m, n; + double y, xx, res, cor, a, da; + mynumber u; + int m; double retval = 0; @@ -187,7 +120,7 @@ __cos_local (double x, int4 k) retval = (res == res + 1.020 * cor) ? res : cslow2 (x); } /* else if (k < 0x3feb6000) */ - else if (k < 0x400368fd) + else { /* 0.855469 <|x|<2.426265 */ ; y = hp0 - fabs (x); a = y + hp1; @@ -221,73 +154,6 @@ __cos_local (double x, int4 k) } /* else if (k < 0x400368fd) */ - - else if (k < 0x419921FB) - { /* 2.426265<|x|< 105414350 */ - t = (x * hpinv + toint); - xn = t - toint; - v.x = t; - y = (x - xn * mp1) - xn * mp2; - n = v.i[LOW_HALF] & 3; - da = xn * mp3; - a = y - da; - da = (y - a) - da; - eps = fabs (x) * 1.2e-30; - - switch (n) - { - case 1: - case 3: - xx = a * a; - if (n == 1) - { - a = -a; - da = -da; - } - if (xx < 0.01588) - { - res = TAYLOR_SIN (xx, a, da, cor); - cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; - retval = (res == res + cor) ? res : csloww (a, da, x); - } - else - { - if (a > 0) - { - m = 1; - } - else - { - m = 0; - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big); - res = do_sin (u, y, da, &cor); - cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; - retval = ((res == res + cor) ? ((m) ? res : -res) - : csloww1 (a, da, x, m)); - } - break; - - case 0: - case 2: - if (a < 0) - { - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big) + da; - res = do_cos (u, y, &cor); - cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; - retval = ((res == res + cor) ? ((n) ? -res : res) - : csloww2 (a, da, x, n)); - break; - } - } /* else if (k < 0x419921FB ) */ - return retval; } @@ -409,6 +275,182 @@ reduce_and_compute_sincos (double x, double *sinxp, double *cosxp) *cosxp = cosx; } +/* Like sloww and csloww, except that it accepts the quadrant and decides to + call mpsin or mpcos based on that. */ +static double +sincos_sloww (double x, double dx, double orig, int4 n) +{ + double y, t, res, cor, w[2], a, da, xn; + res = TAYLOR_SLOW (x, dx, cor); + + if (cor > 0) + cor = 1.0005 * cor + fabs (orig) * 3.1e-30; + else + cor = 1.0005 * cor - fabs (orig) * 3.1e-30; + + if (res == res + cor) + return res; + + (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); + if (w[1] > 0) + cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; + else + cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; + + if (w[0] == w[0] + cor) + return (x > 0) ? w[0] : -w[0]; + + t = (orig * hpinv + toint); + xn = t - toint; + y = (orig - xn * mp1) - xn * mp2; + da = xn * pp3; + t = y - da; + da = (y - t) - da; + y = xn * pp4; + a = t - y; + da = ((t - a) - y) + da; + + if (n & 2) + { + a = -a; + da = -da; + } + (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); + if (w[1] > 0) + cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; + else + cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; + + if (w[0] == w[0] + cor) + return (a > 0) ? w[0] : -w[0]; + + return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); +} + +/* Like sloww1 and csloww1, except that it accepts the quadrant and decides to + call mpsin or mpcos based on that. */ +static double +sincos_sloww1 (double x, double dx, double orig, int n) +{ + mynumber u; + double w[2], y, cor, res, t; + + t = fabs (x); + if (x < 0) + dx = -dx; + + u.x = big + t; + y = t - (u.x - big); + res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); + + if (res == res + cor) + return (x > 0) ? res : -res; + + __dubsin (t, dx, w); + + if (w[1] > 0) + cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); + else + cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); + + if (w[0] == w[0] + cor) + return (x > 0) ? w[0] : -w[0]; + + return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); +} + +/* Like sloww2 and csloww2, except that it accepts the quadrant and decides to + call mpsin or mpcos based on that. */ +static double +sincos_sloww2 (double x, double dx, double orig, int n) +{ + mynumber u; + double w[2], y, cor, res; + + u.x = big + x; + y = x - (u.x - big); + res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); + + if (res == res + cor) + return (n & 2) ? -res : res; + + __docos (x, dx, w); + + if (w[1] > 0) + cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); + else + cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); + + if (w[0] == w[0] + cor) + return (n & 2) ? -w[0] : w[0]; + + return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); +} + +/* Compute sin (A + DA). cos can be computed by shifting the quadrant N + clockwise. */ +static double +__always_inline +do_sincos_1 (double a, double da, double x, int4 n) +{ + double xx, retval, res, cor, y; + mynumber u; + int m; + double eps = fabs (x) * 1.2e-30; + + switch (n) + { /* quarter of unit circle */ + case 2: + a = -a; + da = -da; + case 0: + xx = a * a; + if (xx < 0.01588) + { + /* Taylor series. */ + res = TAYLOR_SIN (xx, a, da, cor); + cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; + retval = (res == res + cor) ? res : sincos_sloww (a, da, x, n); + } + else + { + double b = a, db = da; + if (a > 0) + m = 1; + else + { + m = 0; + b = -b; + db = -db; + } + u.x = big + b; + y = b - (u.x - big); + res = do_sin (u, y, db, &cor); + cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; + retval = ((res == res + cor) ? ((m) ? res : -res) + : sincos_sloww1 (a, da, x, n)); + } + break; + + case 1: + case 3: + if (a < 0) + { + a = -a; + da = -da; + } + u.x = big + a; + y = a - (u.x - big) + da; + res = do_cos (u, y, &cor); + cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; + retval = ((res == res + cor) ? ((n & 2) ? -res : res) + : sincos_sloww2 (a, da, x, n)); + break; + } + + return retval; +} + void __sincos (double x, double *sinx, double *cosx) { @@ -419,12 +461,29 @@ __sincos (double x, double *sinx, double *cosx) u.x = x; k = 0x7fffffff & u.i[HIGH_HALF]; - if (k < 0x419921FB) + + if (k < 0x400368fd) { *sinx = __sin_local (x, k); *cosx = __cos_local (x, k); return; } + if (k < 0x419921FB) + { + double t = (x * hpinv + toint); + double xn = t - toint; + u.x = t; + double y = (x - xn * mp1) - xn * mp2; + int4 n = u.i[LOW_HALF] & 3; + double da = xn * mp3; + double a = y - da; + da = (y - a) - da; + + *sinx = do_sincos_1 (a, da, x, n); + *cosx = do_sincos_1 (a, da, x, (n + 1) & 3); + + return; + } if (k < 0x42F00000) { double t = (x * hpinv + toint);