From owner-freebsd-numerics@FreeBSD.ORG Mon Jun 17 23:29:50 2013 Return-Path: Delivered-To: freebsd-numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) by hub.freebsd.org (Postfix) with ESMTP id A9A7BA74; Mon, 17 Jun 2013 23:29:50 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail12.syd.optusnet.com.au (mail12.syd.optusnet.com.au [211.29.132.193]) by mx1.freebsd.org (Postfix) with ESMTP id 4597017B6; Mon, 17 Jun 2013 23:29:49 +0000 (UTC) Received: from c122-106-156-23.carlnfd1.nsw.optusnet.com.au (c122-106-156-23.carlnfd1.nsw.optusnet.com.au [122.106.156.23]) by mail12.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id r5HNTfDW017625 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Tue, 18 Jun 2013 09:29:42 +1000 Date: Tue, 18 Jun 2013 09:29:41 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: David Schultz Subject: Re: Implementation for coshl. In-Reply-To: <20130610055834.GA68643@zim.MIT.EDU> Message-ID: <20130618091306.A787@besplex.bde.org> References: <20130610003645.GA16444@troutmask.apl.washington.edu> <20130610110740.V24058@besplex.bde.org> <20130610055834.GA68643@zim.MIT.EDU> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.0 cv=Q6eKePKa c=1 sm=1 a=LM0AswAWfpYA:10 a=kj9zAlcOel0A:10 a=PO7r1zJSAAAA:8 a=JzwRw_2MAAAA:8 a=MtpAjnDH-vAA:10 a=LZ7Z3jF9hK5WwFTye5kA:9 a=CjuIK1q_8ugA:10 a=ebeQFi2P/qHVC0Yw9JDJ4g==:117 Cc: freebsd-numerics@FreeBSD.org, Steve Kargl X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 17 Jun 2013 23:29:50 -0000 [Reply to my latest saved mail in this thread. Mostly not in response to it.] On Sun, 9 Jun 2013, David Schultz wrote: > On Mon, Jun 10, 2013, Bruce Evans wrote: >> On Sun, 9 Jun 2013, Steve Kargl wrote: >>> The former would require a refactoring of >>> s_expl.c into a kernel __kernel_expl(x, hi, lo). I have no plans on >>> pursuing this at the this time. >> >> But you need this soon for __ldexp_exp() and __ldexp_cexp(), which are >> needed for hyperbolic functions (the large args already fixed for float >> and double precision) and for cexp() and trig and hyperbolic complex >> functions. It is much easier to implement these using a kernel. I do >> this only for float precision. > > That's fine for now. We need a better k_exp, too. I believe > Bruce was working on one. The wrappers turned out to be fairly easy. Especially __ldexp_exp*(() and __ldexp_cexp*(), given the main wrapper. Are you going to finish coshl() and sinhl()? I'm getting closer to hacking on them. I'm now in the middle of doing this for cosh(): diff -u2 e_cosh.c~ e_cosh.c % --- e_cosh.c~ Thu May 30 18:14:16 2013 % +++ e_cosh.c Mon Jun 17 23:54:12 2013 % @@ -39,10 +39,23 @@ % #include "math_private.h" % % -static const double one = 1.0, half=0.5, huge = 1.0e300; % +static const double one = 1.0, half = 0.5; % +static const volatile double huge = 1.0e300, tiny = 1.0e-300; Part of fixing old bugs. % +static const double % +/* % + * Domain [-1, 1], range ~[-2.2746e-18, 2.2748e-18]: % + * |cosh(x) - c(x)| < 2**-58.6 % + */ % +C2 = 5.0000000000000000e-1, /* 0x10000000000000.0p-53 */ % +C4 = 4.1666666666665422e-2, /* 0x155555555554a2.0p-57 */ % +C6 = 1.3888888889057476e-3, /* 0x16c16c16c29bca.0p-62 */ % +C8 = 2.4801587216389553e-5, /* 0x1a01a01881edc9.0p-68 */ % +C10 = 2.7557340350974160e-7, /* 0x127e50a5570530.0p-74 */ % +C12 = 2.0873996906899188e-9, /* 0x11ee3d8f036513.0p-81 */ % +C14 = 1.1653016291730710e-11; /* 0x19a010a270ce5a.0p-89 */ New method for small x. % % double % __ieee754_cosh(double x) % { % - double t,w; % + double_t t,x2; % int32_t ix; % double -> double_t prepares for using the kernel and avoids discarding accuracy on i386. % @@ -54,13 +67,14 @@ % if(ix>=0x7ff00000) return x*x; % % - /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ % - if(ix<0x3fd62e43) { % - t = expm1(fabs(x)); % - w = one+t; % - if (ix<0x3c800000) return w; /* cosh(tiny) = 1 */ % - return one+(t*t)/(w+w); % + /* |x| in [0,1), return 1 or c(x) */ Fix comment to match code. Update it for new code. % + if(ix<0x3ff00000) { % + if (ix<0x3c800000) % + return one+tiny; /* cosh(tiny) = 1(+) with inexact */ Fix minor bugs in setting inexact % + x2 = x*x; % + return x2*x2*(x2*x2*x2*(x2*(x2*C14 + C12) + C10) + % + x2*(x2*C8 + C6) + C4) + x2*C2 + 1; More accuracy for |x| < 1. % } % % - /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */ % + /* |x| in [1,22], return (exp(|x|)+1/exp(|x|)/2; */ % if (ix < 0x40360000) { % t = __ieee754_exp(fabs(x)); I just need to change this and later cases to use the kernel and its wrappers. In float precision, this is: % - /* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */ % - if (ix < 0x41100000) { % - t = __ieee754_expf(fabsf(x)); % - return half*t+half/t; % + /* |x| in [1, 12), return accurate exp(|x|)/2+1/exp(|x|)/2 */ % + if (ix<0x41400000) { % + k_hexpf(fabsf(x), &hi, &lo); % + return lo + 0.25F/(hi + lo) + hi; % } Bruce