From owner-freebsd-numerics@FreeBSD.ORG Fri Sep 12 21:15:58 2014 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id 284EFE7C for ; Fri, 12 Sep 2014 21:15:58 +0000 (UTC) Received: from mail-yk0-x22d.google.com (mail-yk0-x22d.google.com [IPv6:2607:f8b0:4002:c07::22d]) (using TLSv1 with cipher ECDHE-RSA-RC4-SHA (128/128 bits)) (Client CN "smtp.gmail.com", Issuer "Google Internet Authority G2" (verified OK)) by mx1.freebsd.org (Postfix) with ESMTPS id DCCC0AA7 for ; Fri, 12 Sep 2014 21:15:57 +0000 (UTC) Received: by mail-yk0-f173.google.com with SMTP id 19so749971ykq.18 for ; Fri, 12 Sep 2014 14:15:57 -0700 (PDT) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=google.com; s=20120113; h=mime-version:from:date:message-id:subject:to:content-type; bh=XRCNJY+dPAKlOTMGGI2H5Vpn/HdzvUkMQS8VVgb5esU=; b=KadNqIgv35EirTotVTL0nNJQy6aeMLsI23GzXGYnw9btvElqAYGW7bkrre1VBQobTL ld6FHn6fcZ3vwZWkFxko80IZGFEkgTYXemCBP+KbLR7764PL+wUWUmHq9J7iD+/0xY5Z E+NSKs/mTUiioiT6KqmsKCv2yl/NwGZfS+uOznw1H1W784Wgmy78iTYDNiQ9UOw/BRGz upSFIXhmOPTkwhE38ZsTuCprzJYgWnW1A0YiTIhukgO86m6UZ6j+OiY93DpmmlsXx90w ctUC15S6AYpvyvp5ktQjarWpK8co9JU5i8xE0rbq60fj9eMPWwBcuW9ZKxog913Zm7bN CmhQ== X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20130820; h=x-gm-message-state:mime-version:from:date:message-id:subject:to :content-type; bh=XRCNJY+dPAKlOTMGGI2H5Vpn/HdzvUkMQS8VVgb5esU=; b=b7TQcakSULNcxGgUgGJ24+FRxtBG6s86coF1RpdJ0xjrv5UxaQQnPJphuX5j87QUs3 fDc2H//QCd/OfcSGRTX5evvJ3yIeLjptGM3aV0X6Hh/ljKNrxprQ+tX+Jt2I61zHb7HZ g7JCvlGCaSougv9GpNY/4Dvl2ec8Ua9UyPhaEO1mIybHcRKTWGfHr6fyNfnds1xPl2nS qq+DWIHg7h8mOVD8NTQNpVTK7z5vybvcy5U3KrcqZp7eMEiPbdABeiYlSAF04BEo66h9 6wh30R/K/bFN/X+4CG2TuVSi0V2WhA3Na0Q6jwhSHm/lIjLal9o+Ljg02+IYVRJnj2jk Ku4w== X-Gm-Message-State: ALoCoQnBongM20T5HM7bS+hcNkcWdjKJzfA3i9m/m2+g/1qhbt7q0l0PFxb4NwYWMfOWHKjjhQZs X-Received: by 10.236.71.227 with SMTP id r63mr12882433yhd.2.1410556557104; Fri, 12 Sep 2014 14:15:57 -0700 (PDT) MIME-Version: 1.0 Received: by 10.170.147.69 with HTTP; Fri, 12 Sep 2014 14:15:37 -0700 (PDT) From: enh Date: Fri, 12 Sep 2014 14:15:37 -0700 Message-ID: Subject: lgamma_r and lgammaf_r return the wrong sign for -0.f To: freebsd-numerics@freebsd.org Content-Type: text/plain; charset=UTF-8 X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.18-1 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: Fri, 12 Sep 2014 21:15:58 -0000 if I pass -0.f to lgammaf_r, the sign returned in *signgamp is 1. this is incorrect --- it should be -1. both lgamma_r and lgammaf_r are affected, but the other special cases in those functions look fine to me. this is fixed in OpenBSD and glibc, but FreeBSD and NetBSD both have the same bug. patch below (whitespace mangled courtesy of gmail). i'd prefer to wait for this to be fixed in FreeBSD and pull down the fix rather than just fix it locally. btw, it looks like you're missing coshl/sinhl/tanhl for ld128 now? (they've been removed from imprecise.c without having ld128 implementations added afaics.) --elliott (Android libc maintainer) Index: lib/msun/src/e_lgammaf_r.c =================================================================== --- lib/msun/src/e_lgammaf_r.c (revision 271480) +++ lib/msun/src/e_lgammaf_r.c (working copy) @@ -140,7 +140,11 @@ *signgamp = 1; ix = hx&0x7fffffff; if(ix>=0x7f800000) return x*x; - if(ix==0) return one/vzero; + if(ix==0) { + if(hx<0) + *signgamp = -1; + return one/vzero; + } if(ix<0x35000000) { /* |x|<2**-21, return -log(|x|) */ if(hx<0) { *signgamp = -1; Index: lib/msun/src/e_lgamma_r.c =================================================================== --- lib/msun/src/e_lgamma_r.c (revision 271480) +++ lib/msun/src/e_lgamma_r.c (working copy) @@ -212,7 +212,11 @@ *signgamp = 1; ix = hx&0x7fffffff; if(ix>=0x7ff00000) return x*x; - if((ix|lx)==0) return one/vzero; + if((ix|lx)==0) { + if(hx<0) + *signgamp = -1; + return one/vzero; + } if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */ if(hx<0) { *signgamp = -1; From owner-freebsd-numerics@FreeBSD.ORG Fri Sep 12 22:05:45 2014 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id 3BC52CDD for ; Fri, 12 Sep 2014 22:05:45 +0000 (UTC) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (Client CN "troutmask.apl.washington.edu", Issuer "troutmask.apl.washington.edu" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id 19006FEB for ; Fri, 12 Sep 2014 22:05:45 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.9/8.14.9) with ESMTP id s8CM5im0097520 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO); Fri, 12 Sep 2014 15:05:44 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.9/8.14.9/Submit) id s8CM5i6d097519; Fri, 12 Sep 2014 15:05:44 -0700 (PDT) (envelope-from sgk) Date: Fri, 12 Sep 2014 15:05:44 -0700 From: Steve Kargl To: enh Subject: Re: lgamma_r and lgammaf_r return the wrong sign for -0.f Message-ID: <20140912220544.GA96714@troutmask.apl.washington.edu> References: MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: User-Agent: Mutt/1.5.23 (2014-03-12) Cc: freebsd-numerics@freebsd.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.18-1 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: Fri, 12 Sep 2014 22:05:45 -0000 On Fri, Sep 12, 2014 at 02:15:37PM -0700, enh via freebsd-numerics wrote: > if I pass -0.f to lgammaf_r, the sign returned in *signgamp is 1. this > is incorrect --- it should be -1. > > both lgamma_r and lgammaf_r are affected, but the other special cases > in those functions look fine to me. > > this is fixed in OpenBSD and glibc, but FreeBSD and NetBSD both have > the same bug. Are you sure FreeBSD has a bug? From n1256.pdf (committee draft of C99) F.9.5.3 The lgamma functions -- lgamma(1) returns +0. -- lgamma(2) returns +0. -- lgamma(x) returns +inf and raises the ``divide-by-zero'' floating-pointr exception for x a negative integer or zero. -- lgamma(-inf) returns +inf. -- lgamma(+inf) returns +inf. See the 3rd bullet. -0 is 0 and -0 is a negative integer. POSIX appears to defer to ISO C. n1570.pdf (committe draft for C11) has (almost?) identical text. > patch below (whitespace mangled courtesy of gmail). i'd prefer to wait > for this to be fixed in FreeBSD and pull down the fix rather than just > fix it locally. I have a bigger patch coming with ld80 and ld128 version of lgammal and lgammal_r. > btw, it looks like you're missing coshl/sinhl/tanhl for ld128 now? > (they've been removed from imprecise.c without having ld128 > implementations added afaics.) The ld80 and ld128 versions are in the same file. troutmask:fvwm:kargl[226] grep LDBL /usr/src/lib/msun/src/e_coshl.c #if LDBL_MAX_EXP != 0x4000 #define BIAS (LDBL_MAX_EXP - 1) #if LDBL_MANT_DIG == 64 #elif LDBL_MANT_DIG == 113 #endif /* LDBL_MANT_DIG == 64 */ GET_LDBL_EXPSIGN(ix,x); if (ix Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id 6BCADF6 for ; Fri, 12 Sep 2014 22:17:24 +0000 (UTC) Received: from mail-yh0-x236.google.com (mail-yh0-x236.google.com [IPv6:2607:f8b0:4002:c01::236]) (using TLSv1 with cipher ECDHE-RSA-RC4-SHA (128/128 bits)) (Client CN "smtp.gmail.com", Issuer "Google Internet Authority G2" (verified OK)) by mx1.freebsd.org (Postfix) with ESMTPS id 2840F175 for ; Fri, 12 Sep 2014 22:17:24 +0000 (UTC) Received: by mail-yh0-f54.google.com with SMTP id z6so700746yhz.41 for ; Fri, 12 Sep 2014 15:17:23 -0700 (PDT) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=google.com; s=20120113; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :cc:content-type; bh=XtaJ6NrXjM9DqybcZTc9CuSaAXjCWQNOQtkmcvHKyz0=; b=El2prp6mhqVruz5t6o7dLGJByix25cU92n0MJ4fnt7WkwlFNI0DutQ6rFCqp5DK2nx Cg3a8cuozXaPshA13q80GVDAiIgWfKrQhhJexog83eyGiDgEIon3lQ/mHbIPSYi6wzD6 RpJ1L4cUrYjuv9Ify1Jwh/+WBvA7KFEYP8Ke3OzKYzjqXM1CxIYb5uwudZBHdpSnyXtT 3mwvij/PQ9IrhF9t6YQsRQcHX/mdaAZVb7B08ZIPwKg2W7Tq152fOtWvLBSPqtuR9xAZ RIlmqbBShsV3dlc+Sx0yLUnFDx/XN39R+Li57dB40qJU4mUU2PC5Ea9iskhTLxCcBJhb mbEQ== X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20130820; h=x-gm-message-state:mime-version:in-reply-to:references:from:date :message-id:subject:to:cc:content-type; bh=XtaJ6NrXjM9DqybcZTc9CuSaAXjCWQNOQtkmcvHKyz0=; b=l0fmagT9NT2GQ90EP2IWALbtqIEGcvKTWNKtccs08QspFSi+mZ17WEAQoYYA5rynmA D9mC890NHI4xEtCQrUfvhDR9MBfHmCep2YMXQUMLMZw27VWRmHTaQlTdq7stHfyRrqty CHLlUm53K4lFNsC9O/jyPcag8bO5FKSOM6ApHLIxmgy412kVEHiDGng1KXcPzoZOiL5q iscl6wWFQKHDMokYJjWU4O/WN2MnpGkuMPV63pc3W2f8jxN/VjWkVS38V2/wqghPr8e0 tKLsyKSUG81cXv4Mwhgzj0W6vl5nrJrXBJmhhXH9NBrSaBiWN1hrlmm8vH+8/4BPSRjZ vT2g== X-Gm-Message-State: ALoCoQkvTajk8K0g4K+ybFa7IzFYcYA2gNcIqS/aLr94lw20Dw5rCka6c8ggt31JwNME6T3pH0DC X-Received: by 10.236.229.138 with SMTP id h10mr3735867yhq.122.1410560243094; Fri, 12 Sep 2014 15:17:23 -0700 (PDT) MIME-Version: 1.0 Received: by 10.170.147.69 with HTTP; Fri, 12 Sep 2014 15:17:02 -0700 (PDT) In-Reply-To: <20140912220544.GA96714@troutmask.apl.washington.edu> References: <20140912220544.GA96714@troutmask.apl.washington.edu> From: enh Date: Fri, 12 Sep 2014 15:17:02 -0700 Message-ID: Subject: Re: lgamma_r and lgammaf_r return the wrong sign for -0.f To: Steve Kargl Content-Type: text/plain; charset=UTF-8 Cc: freebsd-numerics@freebsd.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.18-1 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: Fri, 12 Sep 2014 22:17:24 -0000 On Fri, Sep 12, 2014 at 3:05 PM, Steve Kargl wrote: > On Fri, Sep 12, 2014 at 02:15:37PM -0700, enh via freebsd-numerics wrote: >> if I pass -0.f to lgammaf_r, the sign returned in *signgamp is 1. this >> is incorrect --- it should be -1. >> >> both lgamma_r and lgammaf_r are affected, but the other special cases >> in those functions look fine to me. >> >> this is fixed in OpenBSD and glibc, but FreeBSD and NetBSD both have >> the same bug. > > Are you sure FreeBSD has a bug? From n1256.pdf (committee draft of C99) note that i'm reporting bugs in lgamma_r and lgammaf_r. lgamma and lgammaf don't have an out argument for the sign. > F.9.5.3 The lgamma functions > -- lgamma(1) returns +0. > -- lgamma(2) returns +0. > -- lgamma(x) returns +inf and raises the ``divide-by-zero'' floating-pointr > exception for x a negative integer or zero. > -- lgamma(-inf) returns +inf. > -- lgamma(+inf) returns +inf. > > See the 3rd bullet. -0 is 0 and -0 is a negative integer. > > POSIX appears to defer to ISO C. n1570.pdf (committe draft for C11) > has (almost?) identical text. > >> patch below (whitespace mangled courtesy of gmail). i'd prefer to wait >> for this to be fixed in FreeBSD and pull down the fix rather than just >> fix it locally. > > I have a bigger patch coming with ld80 and ld128 version of lgammal > and lgammal_r. sorry, i haven't even looked at the *l variants. >> btw, it looks like you're missing coshl/sinhl/tanhl for ld128 now? >> (they've been removed from imprecise.c without having ld128 >> implementations added afaics.) > > The ld80 and ld128 versions are in the same file. > > troutmask:fvwm:kargl[226] grep LDBL /usr/src/lib/msun/src/e_coshl.c > #if LDBL_MAX_EXP != 0x4000 > #define BIAS (LDBL_MAX_EXP - 1) > #if LDBL_MANT_DIG == 64 > #elif LDBL_MANT_DIG == 113 > #endif /* LDBL_MANT_DIG == 64 */ > GET_LDBL_EXPSIGN(ix,x); > if (ix #if LDBL_MANT_DIG == 64 > #elif LDBL_MANT_DIG == 113 > > see also e_sinhl.c and s_tanhl.c. d'oh. yes, i need to import those files! --elliott From owner-freebsd-numerics@FreeBSD.ORG Fri Sep 12 22:41:13 2014 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id 86E03EB1 for ; Fri, 12 Sep 2014 22:41:13 +0000 (UTC) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (Client CN "troutmask.apl.washington.edu", Issuer "troutmask.apl.washington.edu" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id 493023D6 for ; Fri, 12 Sep 2014 22:41:13 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.9/8.14.9) with ESMTP id s8CMfCKk097730 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO); Fri, 12 Sep 2014 15:41:12 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.9/8.14.9/Submit) id s8CMfCWm097729; Fri, 12 Sep 2014 15:41:12 -0700 (PDT) (envelope-from sgk) Date: Fri, 12 Sep 2014 15:41:12 -0700 From: Steve Kargl To: enh Subject: Re: lgamma_r and lgammaf_r return the wrong sign for -0.f Message-ID: <20140912224112.GA97578@troutmask.apl.washington.edu> References: <20140912220544.GA96714@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: User-Agent: Mutt/1.5.23 (2014-03-12) Cc: freebsd-numerics@freebsd.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.18-1 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: Fri, 12 Sep 2014 22:41:13 -0000 On Fri, Sep 12, 2014 at 03:17:02PM -0700, enh via freebsd-numerics wrote: > On Fri, Sep 12, 2014 at 3:05 PM, Steve Kargl > > F.9.5.3 The lgamma functions > > -- lgamma(1) returns +0. > > -- lgamma(2) returns +0. > > -- lgamma(x) returns +inf and raises the ``divide-by-zero'' floating-point > > exception for x a negative integer or zero. > > -- lgamma(-inf) returns +inf. > > -- lgamma(+inf) returns +inf. > > > > See the 3rd bullet. -0 is 0 and -0 is a negative integer. Of course, neither POSIX nor ISO C specify lgamma_r() only lgamma. I just spent too much time on the 'divide-by-zero' bug, and conflated that with signgam. > > POSIX appears to defer to ISO C. n1570.pdf (committe draft for C11) > > has (almost?) identical text. > > > >> patch below (whitespace mangled courtesy of gmail). i'd prefer to wait > >> for this to be fixed in FreeBSD and pull down the fix rather than just > >> fix it locally. > > > > I have a bigger patch coming with ld80 and ld128 version of lgammal > > and lgammal_r. > > sorry, i haven't even looked at the *l variants. I only just finished writing the *l variants. It took me a long time to unwind the comments in e_lgamma_r.c, so that I could write the long double versions. -- Steve From owner-freebsd-numerics@FreeBSD.ORG Fri Sep 12 22:52:13 2014 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [8.8.178.115]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id C35BE5C4 for ; Fri, 12 Sep 2014 22:52:13 +0000 (UTC) Received: from mail-yk0-x230.google.com (mail-yk0-x230.google.com [IPv6:2607:f8b0:4002:c07::230]) (using TLSv1 with cipher ECDHE-RSA-RC4-SHA (128/128 bits)) (Client CN "smtp.gmail.com", Issuer "Google Internet Authority G2" (verified OK)) by mx1.freebsd.org (Postfix) with ESMTPS id 4F7A478C for ; Fri, 12 Sep 2014 22:52:13 +0000 (UTC) Received: by mail-yk0-f176.google.com with SMTP id 131so797204ykp.7 for ; Fri, 12 Sep 2014 15:52:12 -0700 (PDT) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=google.com; s=20120113; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :cc:content-type; bh=SBM9TEcw4kr5eGifVq7ae5LfltAR3pgwULdRUe9opAY=; b=bvhsDXjvusuKisBYeQvpCa/MHfy0vQSZVAFeZbYqZG3gG2qWl3KeRpT+pw5dRcGchp 9L23qf1oIuVk9y3XO6UusIcYlO+P2HGR1QbFOVa6+O1cBlMvsyP6/EU05lkABY9D53ty fPiLEO8kFZ6g5PK7ZthfKtxXujh9N6qAX7d6UqPNwNePouboCHvuvNuITzniL5SbZ0se o+qGPQSJ8T32IYhNctKLuxlqqZi3fWu/ba0ooWjV5D31WIsu3MOD/MjJbm7nIpsLzCTu NMgBHp/IthNPX8OghtIgU0BoRcDsu56/7GyzQVRt/3M6qsyUiaIa29kdslzfuqPBs3el lPjA== X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20130820; h=x-gm-message-state:mime-version:in-reply-to:references:from:date :message-id:subject:to:cc:content-type; bh=SBM9TEcw4kr5eGifVq7ae5LfltAR3pgwULdRUe9opAY=; b=I6Ne6dIItEqYl/NKzIA0aW2ItDESSEvvFYw9RPTePB7jb9GSrTX2a3rzDZApuT2tMI D8l+k9wFre50IuV3OGjWnwD31mhxI2Bv4qd8o6K5oLw/nD755ORlJmZr/Cp1Ye1ADURn n/Yv2nio+OT2enh64m6+b0tCg4TCiMvI2tycJqGe/vxPSJFGCm0Wp/uC1a8c9TA7xLas LA154Xlb98tqUNMyTdhhISjBwvojNDYFirCx4/wpUQAt8PHGomGhN6hJPuGenJms1VIo +hss3nZrxDWGLf+y5Cie0STwgFelNr9rKYiUhM6Mt4gLDV9d0/AQMGVaJzKwRe9HM06B 271w== X-Gm-Message-State: ALoCoQmKfCgU1Gm+fPsWymTx7k1l7Wi6NNIIBoEezzKJGWV/7/hjp/fD4lZFiz0KHcaCCNg6npAL X-Received: by 10.236.85.208 with SMTP id u56mr13686400yhe.48.1410562332361; Fri, 12 Sep 2014 15:52:12 -0700 (PDT) MIME-Version: 1.0 Received: by 10.170.147.69 with HTTP; Fri, 12 Sep 2014 15:51:52 -0700 (PDT) In-Reply-To: <20140912224112.GA97578@troutmask.apl.washington.edu> References: <20140912220544.GA96714@troutmask.apl.washington.edu> <20140912224112.GA97578@troutmask.apl.washington.edu> From: enh Date: Fri, 12 Sep 2014 15:51:52 -0700 Message-ID: Subject: Re: lgamma_r and lgammaf_r return the wrong sign for -0.f To: Steve Kargl Content-Type: text/plain; charset=UTF-8 Cc: freebsd-numerics@freebsd.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.18-1 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: Fri, 12 Sep 2014 22:52:14 -0000 On Fri, Sep 12, 2014 at 3:41 PM, Steve Kargl wrote: > On Fri, Sep 12, 2014 at 03:17:02PM -0700, enh via freebsd-numerics wrote: >> On Fri, Sep 12, 2014 at 3:05 PM, Steve Kargl >> > F.9.5.3 The lgamma functions >> > -- lgamma(1) returns +0. >> > -- lgamma(2) returns +0. >> > -- lgamma(x) returns +inf and raises the ``divide-by-zero'' floating-point >> > exception for x a negative integer or zero. >> > -- lgamma(-inf) returns +inf. >> > -- lgamma(+inf) returns +inf. >> > >> > See the 3rd bullet. -0 is 0 and -0 is a negative integer. > > Of course, neither POSIX nor ISO C specify lgamma_r() only lgamma. > I just spent too much time on the 'divide-by-zero' bug, and > conflated that with signgam. i think the reporter's feeling was that *signgamp should be -1 because as x approaches 0 from the negative side, Gamma(x) approaches -Infinity. (though i suspect they noticed this when porting code from glibc.) >> > POSIX appears to defer to ISO C. n1570.pdf (committe draft for C11) >> > has (almost?) identical text. >> > >> >> patch below (whitespace mangled courtesy of gmail). i'd prefer to wait >> >> for this to be fixed in FreeBSD and pull down the fix rather than just >> >> fix it locally. >> > >> > I have a bigger patch coming with ld80 and ld128 version of lgammal >> > and lgammal_r. >> >> sorry, i haven't even looked at the *l variants. > > I only just finished writing the *l variants. It took me > a long time to unwind the comments in e_lgamma_r.c, > so that I could write the long double versions. is there a freebsd-numerics-commits or equivalent i could subscribe to? From owner-freebsd-numerics@FreeBSD.ORG Fri Sep 12 23:37:50 2014 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [8.8.178.115]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id ADB3D33A for ; Fri, 12 Sep 2014 23:37:50 +0000 (UTC) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (Client CN "troutmask.apl.washington.edu", Issuer "troutmask.apl.washington.edu" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id 84AC6B16 for ; Fri, 12 Sep 2014 23:37:50 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.9/8.14.9) with ESMTP id s8CNbnpA098025 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO); Fri, 12 Sep 2014 16:37:49 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.9/8.14.9/Submit) id s8CNbnuL098024; Fri, 12 Sep 2014 16:37:49 -0700 (PDT) (envelope-from sgk) Date: Fri, 12 Sep 2014 16:37:49 -0700 From: Steve Kargl To: enh Subject: Re: lgamma_r and lgammaf_r return the wrong sign for -0.f Message-ID: <20140912233749.GA97910@troutmask.apl.washington.edu> References: <20140912220544.GA96714@troutmask.apl.washington.edu> <20140912224112.GA97578@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: User-Agent: Mutt/1.5.23 (2014-03-12) Cc: freebsd-numerics@freebsd.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.18-1 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: Fri, 12 Sep 2014 23:37:50 -0000 On Fri, Sep 12, 2014 at 03:51:52PM -0700, enh wrote: > On Fri, Sep 12, 2014 at 3:41 PM, Steve Kargl > wrote: > > On Fri, Sep 12, 2014 at 03:17:02PM -0700, enh via freebsd-numerics wrote: > >> On Fri, Sep 12, 2014 at 3:05 PM, Steve Kargl > >> > F.9.5.3 The lgamma functions > >> > -- lgamma(1) returns +0. > >> > -- lgamma(2) returns +0. > >> > -- lgamma(x) returns +inf and raises the ``divide-by-zero'' floating-point > >> > exception for x a negative integer or zero. > >> > -- lgamma(-inf) returns +inf. > >> > -- lgamma(+inf) returns +inf. > >> > > >> > See the 3rd bullet. -0 is 0 and -0 is a negative integer. > > > > Of course, neither POSIX nor ISO C specify lgamma_r() only lgamma. > > I just spent too much time on the 'divide-by-zero' bug, and > > conflated that with signgam. > > i think the reporter's feeling was that *signgamp should be -1 because > as x approaches 0 from the negative side, Gamma(x) approaches > -Infinity. (though i suspect they noticed this when porting code from > glibc.) Mathematically, gamma(x) -> +inf as x -> 0 for x > 0 and gamma(x) -> -inf as x -> 0 for x < 0. Assuming that a processor supports a signed zero, it happens that x = +-0 is the only integral floating point value that allows one to infer a sign for signgam. All other negative integral values, of course, do not have this possible, so signgam is arbitrarily assigned 1. > >> > POSIX appears to defer to ISO C. n1570.pdf (committe draft for C11) > >> > has (almost?) identical text. > >> > > >> >> patch below (whitespace mangled courtesy of gmail). i'd prefer to wait > >> >> for this to be fixed in FreeBSD and pull down the fix rather than just > >> >> fix it locally. > >> > > >> > I have a bigger patch coming with ld80 and ld128 version of lgammal > >> > and lgammal_r. > >> > >> sorry, i haven't even looked at the *l variants. > > > > I only just finished writing the *l variants. It took me > > a long time to unwind the comments in e_lgamma_r.c, > > so that I could write the long double versions. > > is there a freebsd-numerics-commits or equivalent i could subscribe to? No. I tend to post my patches here. Then Bruce Evans and I tend to go off on private email threads with code reviews. A few months later I get around to address issues raised by Bruce, and submit a new patch. This cycle goes on for a year or so, and then I finally commit my patches. There are also a few developers (most MIPS/ARMS) that have been fixing fenv.h issues. They commit their patches without posting to freebsd-numerics. I'll note that my lgamma.diff file is clocking in at 1400 LOC. mailman may eat my post. In the end, it is probably prudent to simply look periodically at FreeBSD's src/lib/msun through the svnweb interface. -- Steve From owner-freebsd-numerics@FreeBSD.ORG Sat Sep 13 06:00:33 2014 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [8.8.178.115]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id 08EE68C2 for ; Sat, 13 Sep 2014 06:00:33 +0000 (UTC) Received: from mail106.syd.optusnet.com.au (mail106.syd.optusnet.com.au [211.29.132.42]) by mx1.freebsd.org (Postfix) with ESMTP id BBC26EEE for ; Sat, 13 Sep 2014 06:00:32 +0000 (UTC) Received: from c122-106-147-133.carlnfd1.nsw.optusnet.com.au (c122-106-147-133.carlnfd1.nsw.optusnet.com.au [122.106.147.133]) by mail106.syd.optusnet.com.au (Postfix) with ESMTPS id 3E2503C3083; Sat, 13 Sep 2014 16:00:23 +1000 (EST) Date: Sat, 13 Sep 2014 16:00:22 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: enh Subject: Re: lgamma_r and lgammaf_r return the wrong sign for -0.f In-Reply-To: Message-ID: <20140913152505.E950@besplex.bde.org> References: <20140912220544.GA96714@troutmask.apl.washington.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.1 cv=AOuw8Gd4 c=1 sm=1 tr=0 a=7NqvjVvQucbO2RlWB8PEog==:117 a=PO7r1zJSAAAA:8 a=srA-gRaI36IA:10 a=lreBrSRemiwA:10 a=kj9zAlcOel0A:10 a=JzwRw_2MAAAA:8 a=UaAATAUXzxPjzBn7DBQA:9 a=CjuIK1q_8ugA:10 Cc: freebsd-numerics@freebsd.org, Steve Kargl X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.18-1 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: Sat, 13 Sep 2014 06:00:33 -0000 On Fri, 12 Sep 2014, enh via freebsd-numerics wrote: > On Fri, Sep 12, 2014 at 3:05 PM, Steve Kargl > wrote: >> On Fri, Sep 12, 2014 at 02:15:37PM -0700, enh via freebsd-numerics wrote: >>> if I pass -0.f to lgammaf_r, the sign returned in *signgamp is 1. this >>> is incorrect --- it should be -1. >>> >>> both lgamma_r and lgammaf_r are affected, but the other special cases >>> in those functions look fine to me. >>> >>> this is fixed in OpenBSD and glibc, but FreeBSD and NetBSD both have >>> the same bug. >> >> Are you sure FreeBSD has a bug? From n1256.pdf (committee draft of C99) C99 doesn't even have signgam. But it has the usual complications for +-0 when the sign of the result can be determined for a real function, so it has tgamma(+-0) = +-Inf. POSIX specifies that signgam has the same sign as gamma(), where gamma is spelled with a greek letter. (Bug; it should say the same sign as tgamma(). The one spelled with the greek letter should mean the mathematical one, and that has a complex pole at 0. Complex poles don't have signs.) I looked at what tgamma() does. It is correct for +-0, but for large finite args and +Inf it returns 1/zero. This is only correct for +Inf. >>> patch below (whitespace mangled courtesy of gmail). i'd prefer to wait >>> for this to be fixed in FreeBSD and pull down the fix rather than just >>> fix it locally. The patch can possibly be improved a little by duplicating the zero check ((ix|lx) == 0) instead of duplicating the sign check (hx < 0). This gives the same number of checks, but 1 fewer setting of signgam. In logl(), the checks are arranged more subtly to try to optimize for the usual case. This works best in ld80 long double precision. The single check (hx <= 0) classifies all negative and 0 x. This saves as much as 1 whole cycle or 1.5% on modern x86. It works best in for ld80 since hx doesn't contain any mantissa bits then. When hx contains mantissa bits, its sign bit is usually checked using the C '<' operator as in (hx < 0) and not as a bit. Sometimes the compiler can combine this with another operator and do it as a bit test. The ((ix|lx) == 0) check may be too special for this to be possible except in 64-bit mode where hx, ix, and lx can be treated as a single 64-bit value (then (hlx = 0) works, where int64_t hlx is hx and lx combined). The optimization is going the other way an converting several bit tests into a single signed comparison. Bruce From owner-freebsd-numerics@FreeBSD.ORG Sat Sep 13 19:36:36 2014 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id 63A95CAE for ; Sat, 13 Sep 2014 19:36:36 +0000 (UTC) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (Client CN "troutmask.apl.washington.edu", Issuer "troutmask.apl.washington.edu" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id 39FB5F71 for ; Sat, 13 Sep 2014 19:36:36 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.9/8.14.9) with ESMTP id s8DJaYv9002358 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO) for ; Sat, 13 Sep 2014 12:36:34 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.9/8.14.9/Submit) id s8DJaYKL002357 for freebsd-numerics@freebsd.org; Sat, 13 Sep 2014 12:36:34 -0700 (PDT) (envelope-from sgk) Date: Sat, 13 Sep 2014 12:36:34 -0700 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: lgammal[_r] patch Message-ID: <20140913193634.GA2323@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline User-Agent: Mutt/1.5.23 (2014-03-12) X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.18-1 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: Sat, 13 Sep 2014 19:36:36 -0000 Following my .sig is a patch that contains the ld80 and ld128 implementations of lgamma[_r](). It also contains an update to lgammaf_r() and lgamma_r(). Specifically, it has * Makefile: . Attach e_lgammal[_r].c to the build. . Create links for lgammal[_r].3. * Symbol.map: . Sort lgammal into its proper place. . Add FBSD_1.4 to contain the new lgammal_r symbol. * ld128/e_lgammal_r.c: . 128-bit implementataion of lgammal_r(). * ld80/e_lgammal_r.c: . Intel 80-bit format implementation of lgammal_r(). * man/lgamma.3: . Document the new functions. * src/e_lgamma.c: . Expose lgammal as a weak reference for platforms where long double is mapped to double. * src/e_lgamma_r.c: . Use integer literal constants instead of real literal constants. Let compiler(s) do the job of conversion to the appropriate type. . Expose lgammal_r as a weak reference for platforms where long double is mapped to double. * src/e_lgammaf_r.c: . Fixed the Cygnus Support conversion of e_lgamma_r.c to float. This includes the generation of new polynomial and rational approximations with fewer terms. For each approximation, include a comment on an estimate of the accuracy over the relevant range. . Use integer literal constants instead of real literal constants. Let compiler(s) do the job of conversion to the appropriate type. This allows the removal of several explicit casts of double values to float. * src/e_lgammal.c: . Wrapper for lgammal() about lgammal_r(). * src/imprecise.c: . Remove the lgamma. * src/math.h: . Add a prototype for lgammal_r(). -- Steve Index: Makefile =================================================================== --- Makefile (revision 271479) +++ Makefile (working copy) @@ -98,6 +98,7 @@ # If long double != double use these; otherwise, we alias the double versions. COMMON_SRCS+= e_acoshl.c e_acosl.c e_asinl.c e_atan2l.c e_atanhl.c \ e_coshl.c e_fmodl.c e_hypotl.c \ + e_lgammal.c e_lgammal_r.c \ e_remainderl.c e_sinhl.c e_sqrtl.c \ invtrig.c k_cosl.c k_sinl.c k_tanl.c \ s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.c \ @@ -188,7 +189,8 @@ ilogb.3 logb.3 ilogb.3 logbf.3 ilogb.3 logbl.3 MLINKS+=j0.3 j1.3 j0.3 jn.3 j0.3 y0.3 j0.3 y1.3 j0.3 y1f.3 j0.3 yn.3 MLINKS+=j0.3 j0f.3 j0.3 j1f.3 j0.3 jnf.3 j0.3 y0f.3 j0.3 ynf.3 -MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 lgamma.3 lgammaf.3 \ +MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 \ + lgamma.3 lgammaf.3 lgamma.3 lgammal.3 \ lgamma.3 tgamma.3 lgamma.3 tgammaf.3 MLINKS+=log.3 log10.3 log.3 log10f.3 log.3 log10l.3 \ log.3 log1p.3 log.3 log1pf.3 log.3 log1pl.3 \ Index: Symbol.map =================================================================== --- Symbol.map (revision 271479) +++ Symbol.map (working copy) @@ -269,6 +269,7 @@ erfl; expl; expm1l; + lgammal; log10l; log1pl; log2l; @@ -276,7 +277,11 @@ sinhl; tanhl; /* Implemented as weak aliases for imprecise versions */ - lgammal; powl; tgammal; }; + +/* First added in 11.0-CURRENT */ +FBSD_1.4 { + lgammal_r; +}; Index: ld128/e_lgammal_r.c =================================================================== --- ld128/e_lgammal_r.c (revision 0) +++ ld128/e_lgammal_r.c (working copy) @@ -0,0 +1,327 @@ +/* @(#)e_lgamma_r.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include +__FBSDID("$FreeBSD$"); +/* + * See e_lgamma_r.c for complete comments. + * + * Converted to long double by Steven G. Kargl. + */ +#include "fpmath.h" +#include "math.h" +#include "math_private.h" + +static const volatile double vzero = 0; + +static const double +zero= 0, +half= 0.5, +one = 1; + +static const long double +pi = 3.14159265358979323846264338327950288e+00L; +/* + * Domain y in [0x1p-119, 0.28], range ~[-1.4065e-36, 1.4065e-36]: + * |(lgamma(2 - y) + y / 2) / y - a(y)| < 2**-119.1 + */ +static const long double +a0 = 7.72156649015328606065120900824024296e-02L, +a1 = 3.22467033424113218236207583323018498e-01L, +a2 = 6.73523010531980951332460538330282217e-02L, +a3 = 2.05808084277845478790009252803463129e-02L, +a4 = 7.38555102867398526627292839296001626e-03L, +a5 = 2.89051033074152328576829509522483468e-03L, +a6 = 1.19275391170326097618357349881842913e-03L, +a7 = 5.09669524743042462515256340206203019e-04L, +a8 = 2.23154758453578096143609255559576017e-04L, +a9 = 9.94575127818397632126978731542755129e-05L, +a10 = 4.49262367375420471287545895027098145e-05L, +a11 = 2.05072127845117995426519671481628849e-05L, +a12 = 9.43948816959096748454087141447939513e-06L, +a13 = 4.37486780697359330303852050718287419e-06L, +a14 = 2.03920783892362558276037363847651809e-06L, +a15 = 9.55191070057967287877923073200324649e-07L, +a16 = 4.48993286185740853170657139487620560e-07L, +a17 = 2.13107543597620911675316728179563522e-07L, +a18 = 9.70745379855304499867546549551023473e-08L, +a19 = 5.61889970390290257926487734695402075e-08L, +a20 = 6.42739653024130071866684358960960951e-09L, +a21 = 3.34491062143649291746195612991870119e-08L, +a22 = -1.57068547394315223934653011440641472e-08L, +a23 = 1.30812825422415841213733487745200632e-08L; +/* + * Domain x in [tc-0.24, tc+0.28], range ~[-6.3201e-37, 6.3201e-37]: + * |(lgamma(x) - tf) - t(x - tc)| < 2**-120.3. + */ +static const long double +tc = 1.46163214496836234126265954232572133e+00L, +tf = -1.21486290535849608095514557177691584e-01L, +tt = 1.57061739945077675484237837992951704e-36L, +t0 = -1.99238329499314692728655623767019240e-36L, +t1 = -6.08453430711711404116887457663281416e-35L, +t2 = 4.83836122723810585213722380854828904e-01L, +t3 = -1.47587722994530702030955093950668275e-01L, +t4 = 6.46249402389127526561003464202671923e-02L, +t5 = -3.27885410884813055008502586863748063e-02L, +t6 = 1.79706751152103942928638276067164935e-02L, +t7 = -1.03142230366363872751602029672767978e-02L, +t8 = 6.10053602051788840313573150785080958e-03L, +t9 = -3.68456960831637325470641021892968954e-03L, +t10 = 2.25976482322181046611440855340968560e-03L, +t11 = -1.40225144590445082933490395950664961e-03L, +t12 = 8.78232634717681264035014878172485575e-04L, +t13 = -5.54194952796682301220684760591403899e-04L, +t14 = 3.51912956837848209220421213975000298e-04L, +t15 = -2.24653443695947456542669289367055542e-04L, +t16 = 1.44070395420840737695611929680511823e-04L, +t17 = -9.27609865550394140067059487518862512e-05L, +t18 = 5.99347334438437081412945428365433073e-05L, +t19 = -3.88458388854572825603964274134801009e-05L, +t20 = 2.52476631610328129217896436186551043e-05L, +t21 = -1.64508584981658692556994212457518536e-05L, +t22 = 1.07434583475987007495523340296173839e-05L, +t23 = -7.03070407519397260929482550448878399e-06L, +t24 = 4.60968590693753579648385629003100469e-06L, +t25 = -3.02765473778832036018438676945512661e-06L, +t26 = 1.99238771545503819972741288511303401e-06L, +t27 = -1.31281299822614084861868817951788579e-06L, +t28 = 8.60844432267399655055574642052370223e-07L, +t29 = -5.64535486432397413273248363550536374e-07L, +t30 = 3.99357783676275660934903139592727737e-07L, +t31 = -2.95849029193433121795495215869311610e-07L, +t32 = 1.37790144435073124976696250804940384e-07L; +/* + * Domain y in [-0.1, 0.232], range ~[-1.4046e-37, 1.4181e-37]: + * |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-122.8 + */ +static const long double +u0 = -7.72156649015328606065120900824024311e-02L, +u1 = 4.24082772271938167430983113242482656e-01L, +u2 = 2.96194003481457101058321977413332171e+00L, +u3 = 6.49503267711258043997790983071543710e+00L, +u4 = 7.40090051288150177152835698948644483e+00L, +u5 = 4.94698036296756044610805900340723464e+00L, +u6 = 2.00194224610796294762469550684947768e+00L, +u7 = 4.82073087750608895996915051568834949e-01L, +u8 = 6.46694052280506568192333848437585427e-02L, +u9 = 4.17685526755100259316625348933108810e-03L, +u10 = 9.06361003550314327144119307810053410e-05L, +v1 = 5.15937098592887275994320496999951947e+00L, +v2 = 1.14068418766251486777604403304717558e+01L, +v3 = 1.41164839437524744055723871839748489e+01L, +v4 = 1.07170702656179582805791063277960532e+01L, +v5 = 5.14448694179047879915042998453632434e+00L, +v6 = 1.55210088094585540637493826431170289e+00L, +v7 = 2.82975732849424562719893657416365673e-01L, +v8 = 2.86424622754753198010525786005443539e-02L, +v9 = 1.35364253570403771005922441442688978e-03L, +v10 = 1.91514173702398375346658943749580666e-05L, +v11 = -3.25364686890242327944584691466034268e-08L; +/* + * Domain x in (2, 3], range ~[-1.3341e-36, 1.3536e-36]: + * |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-120.1 + * with y = x - 2. + */ +static const long double +s0 = -7.72156649015328606065120900824024297e-02L, +s1 = 1.23221687850916448903914170805852253e-01L, +s2 = 5.43673188699937239808255378293820020e-01L, +s3 = 6.31998137119005233383666791176301800e-01L, +s4 = 3.75885340179479850993811501596213763e-01L, +s5 = 1.31572908743275052623410195011261575e-01L, +s6 = 2.82528453299138685507186287149699749e-02L, +s7 = 3.70262021550340817867688714880797019e-03L, +s8 = 2.83374000312371199625774129290973648e-04L, +s9 = 1.15091830239148290758883505582343691e-05L, +s10 = 2.04203474281493971326506384646692446e-07L, +s11 = 9.79544198078992058548607407635645763e-10L, +r1 = 2.58037466655605285937112832039537492e+00L, +r2 = 2.86289413392776399262513849911531180e+00L, +r3 = 1.78691044735267497452847829579514367e+00L, +r4 = 6.89400381446725342846854215600008055e-01L, +r5 = 1.70135865462567955867134197595365343e-01L, +r6 = 2.68794816183964420375498986152766763e-02L, +r7 = 2.64617234244861832870088893332006679e-03L, +r8 = 1.52881761239180800640068128681725702e-04L, +r9 = 4.63264813762296029824851351257638558e-06L, +r10 = 5.89461519146957343083848967333671142e-08L, +r11 = 1.79027678176582527798327441636552968e-10L; +/* + * Domain z in [8, 0x1p70], range ~[-9.8214e-35, 9.8214e-35]: + * |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-113.0 + */ +static const long double +w0 = 4.18938533204672741780329736405617738e-01L, +w1 = 8.33333333333333333333333333332852026e-02L, +w2 = -2.77777777777777777777777727810123528e-03L, +w3 = 7.93650793650793650791708939493907380e-04L, +w4 = -5.95238095238095234390450004444370959e-04L, +w5 = 8.41750841750837633887817658848845695e-04L, +w6 = -1.91752691752396849943172337347259743e-03L, +w7 = 6.41025640880333069429106541459015557e-03L, +w8 = -2.95506530801732133437990433080327074e-02L, +w9 = 1.79644237328444101596766586979576927e-01L, +w10 = -1.39240539108367641920172649259736394e+00L, +w11 = 1.33987701479007233325288857758641761e+01L, +w12 = -1.56363596431084279780966590116006255e+02L, +w13 = 2.14830978044410267201172332952040777e+03L, +w14 = -3.28636067474227378352761516589092334e+04L, +w15 = 5.06201257747865138432663574251462485e+05L, +w16 = -6.79720123352023636706247599728048344e+06L, +w17 = 6.57556601705472106989497289465949255e+07L, +w18 = -3.26229058141181783534257632389415580e+08L; + +static long double +sin_pil(long double x) +{ + volatile long double vz; + long double y,z; + uint64_t lx, n; + uint16_t hx; + + y = -x; + + vz = y+0x1.p112; + z = vz-0x1.p112; + if (z == y) + return zero; + + vz = y+0x1.p110; + EXTRACT_LDBL128_WORDS(hx,lx,n,vz); + z = vz-0x1.p110; + if (z > y) { + z -= 0.25; + n--; + } + n &= 7; + y = y - z + n * 0.25L; + + switch (n) { + case 0: y = __kernel_sinl(pi*y,zero,0); break; + case 1: + case 2: y = __kernel_cosl(pi*(0.5-y),zero); break; + case 3: + case 4: y = __kernel_sinl(pi*(one-y),zero,0); break; + case 5: + case 6: y = -__kernel_cosl(pi*(y-1.5),zero); break; + default: y = __kernel_sinl(pi*(y-2.0),zero,0); break; + } + return -y; +} + + +long double +lgammal_r(long double x, int *signgamp) +{ + long double nadj,p,p1,p2,p3,q,r,t,w,y,z; + uint64_t llx,lx; + int i; + uint16_t hx; + + EXTRACT_LDBL128_WORDS(hx, lx, llx, x); + + if((hx & 0x7fff) == 0x7fff) { /* erfl(nan)=nan */ + i = (hx>>15)<<1; + return (1-i)+one/x; /* erfl(+-inf)=+-1 */ + } + + /* purge off +-inf, NaN, +-0, tiny and negative arguments */ + *signgamp = 1; + if((hx & 0x7fff) == 0x7fff) /* x is +-Inf or NaN */ + return x*x; + if((hx==0||hx==0x8000)&&lx==0) return one/vzero; + + /* purge off tiny and negative arguments */ + if(fabsl(x)<0x1p-119L) { + if(hx&0x8000) { + *signgamp = -1; + return -logl(-x); + } else return -logl(x); + } + if(hx&0x8000) { + if(fabsl(x)>=0x1p112) + return one/vzero; + t = sin_pil(x); + if(t==zero) return one/vzero; + nadj = logl(pi/fabsl(t*x)); + if(t=0.7315998077392578) {y = 1-x; i= 0;} + else if(x>=0.2316399812698364) {y= x-(tc-1); i=1;} + else {y = x; i=2;} + } else { + r = 0; + if(x>=1.7316312789916992) {y=2-x;i=0;} + else if(x>=1.2316322326660156) {y=x-tc;i=1;} + else {y=x-1;i=2;} + } + switch(i) { + case 0: + z = y*y; + p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*(a10+z*(a12+z*(a14+z*(a16+ + z*(a18+z*(a20+z*a22)))))))))); + p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*(a13+z*(a15+ + z*(a17+z*(a19+z*(a21+z*a23))))))))))); + p = y*p1+p2; + r += (p-y/2); break; + case 1: + p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+ + y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+ + y*(t17+y*(t18+y*(t19+y*(t20+y*(t21+y*(t22+y*(t23+ + y*(t24+y*(t25+y*(t26+y*(t27+y*(t28+y*(t29+y*(t30+ + y*(t31+y*t32)))))))))))))))))))))))))))))); + r += (tf + p); break; + case 2: + p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*(u6+y*(u7+ + y*(u8+y*(u9+y*u10)))))))))); + p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*(v6+y*(v7+ + y*(v8+y*(v9+y*(v10+y*v11)))))))))); + r += (-y/2 + p1/p2); + } + } + else if(x<8) { + i = x; + y = x-i; + p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*(s6+y*(s7+y*(s8+ + y*(s9+y*(s10+y*s11))))))))))); + q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*(r6+y*(r7+y*(r8+ + y*(r9+y*(r10+y*r11)))))))))); + r = y/2+p/q; + z = 1; /* lgamma(1+s) = log(s) + lgamma(s) */ + switch(i) { + case 7: z *= (y+6); /* FALLTHRU */ + case 6: z *= (y+5); /* FALLTHRU */ + case 5: z *= (y+4); /* FALLTHRU */ + case 4: z *= (y+3); /* FALLTHRU */ + case 3: z *= (y+2); /* FALLTHRU */ + r += logl(z); break; + } + } else if (x < 0x1p119L) { + t = logl(x); + z = one/x; + y = z*z; + w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*(w6+y*(w7+y*(w8+ + y*(w9+y*(w10+y*(w11+y*(w12+y*(w13+y*(w14+y*(w15+y*(w16+ + y*(w17+y*w18))))))))))))))))); + r = (x-half)*(t-one)+w; + } else + r = x*(logl(x)-1); + if(hx&0x8000) r = nadj - r; + return r; +} Property changes on: ld128/e_lgammal_r.c ___________________________________________________________________ Added: svn:mime-type ## -0,0 +1 ## +text/plain \ No newline at end of property Added: svn:keywords ## -0,0 +1 ## +FreeBSD=%H \ No newline at end of property Added: svn:eol-style ## -0,0 +1 ## +native \ No newline at end of property Index: ld80/e_lgammal_r.c =================================================================== --- ld80/e_lgammal_r.c (revision 0) +++ ld80/e_lgammal_r.c (working copy) @@ -0,0 +1,343 @@ +/* @(#)e_lgamma_r.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include +__FBSDID("$FreeBSD"); +/* + * See s_lgamma_r.c for complete comments. + * + * Converted to long double by Steven G. Kargl. + */ +#include +#ifdef __i386__ +#include +#endif + +#include "fpmath.h" +#include "math.h" +#include "math_private.h" + +static const volatile double vzero = 0; + +static const double +zero= 0, +half= 0.5, +one = 1; + +static const union IEEEl2bits +piu = LD80C(0xc90fdaa22168c235, 1, 3.14159265358979323851e+00L); +#define pi (piu.e) +/* + * Domain y in [0x1p-70, 0.27], range ~[-4.5264e-22, 4.5264e-22]: + * |(lgamma(2 - y) + y / 2) / y - a(y)| < 2**-70.9 + */ +static const union IEEEl2bits +a0u = LD80C(0x9e233f1bed863d26, -4, 7.72156649015328606028e-02L), +a1u = LD80C(0xa51a6625307d3249, -2, 3.22467033424113218889e-01L), +a2u = LD80C(0x89f000d2abafda8c, -4, 6.73523010531979398946e-02L), +a3u = LD80C(0xa8991563eca75f26, -6, 2.05808084277991211934e-02L), +a4u = LD80C(0xf2027e10634ce6b6, -8, 7.38555102796070454026e-03L), +a5u = LD80C(0xbd6eb76dd22187f4, -9, 2.89051035162703932972e-03L), +a6u = LD80C(0x9c562ab05e0458ed, -10, 1.19275351624639999297e-03L), +a7u = LD80C(0x859baed93ee48e46, -11, 5.09674593842117925320e-04L), +a8u = LD80C(0xe9f28a4432949af2, -13, 2.23109648015769155122e-04L), +a9u = LD80C(0xd12ad0d9b93c6bb0, -14, 9.97387167479808509830e-05L), +a10u= LD80C(0xb7522643c78a219b, -15, 4.37071076331030136818e-05L), +a11u= LD80C(0xca024dcdece2cb79, -16, 2.40813493372040143061e-05L), +a12u= LD80C(0xbb90fb6968ebdbf9, -19, 2.79495621083634031729e-06L), +a13u= LD80C(0xba1c9ffeeae07b37, -17, 1.10931287015513924136e-05L); +#define a0 (a0u.e) +#define a1 (a1u.e) +#define a2 (a2u.e) +#define a3 (a3u.e) +#define a4 (a4u.e) +#define a5 (a5u.e) +#define a6 (a6u.e) +#define a7 (a7u.e) +#define a8 (a8u.e) +#define a9 (a9u.e) +#define a10 (a10u.e) +#define a11 (a11u.e) +#define a12 (a12u.e) +#define a13 (a13u.e) +/* + * Domain x in [tc-0.24, tc+0.28], range ~[-6.1205e-22, 6.1205e-22]: + * |(lgamma(x) - tf) - t(x - tc)| < 2**-70.5 + */ +static const union IEEEl2bits +tcu = LD80C(0xbb16c31ab5f1fb71, 0, 1.46163214496836234128e+00L), +tfu = LD80C(0xf8cdcde61c520e0f, -4, -1.21486290535849608093e-01L), +ttu = LD80C(0xd46ee54b27d4de99, -69, -2.81152980996018785880e-21L), +t0u = LD80C(0x80b9406556a62a6b, -68, 3.40728634996055147231e-21L), +t1u = LD80C(0xc7e9c6f6df3f8c39, -67, -1.05833162742737073665e-20L), +t2u = LD80C(0xf7b95e4771c55d51, -2, 4.83836122723810583532e-01L), +t3u = LD80C(0x97213c6e35e119ff, -3, -1.47587722994530691476e-01L), +t4u = LD80C(0x845a14a6a81dc94b, -4, 6.46249402389135358063e-02L), +t5u = LD80C(0x864d46fa89997796, -5, -3.27885410884846056084e-02L), +t6u = LD80C(0x93373cbd00297438, -6, 1.79706751150707171293e-02L), +t7u = LD80C(0xa8fcfca7eddc8d1d, -7, -1.03142230361450732547e-02L), +t8u = LD80C(0xc7e7015ff4bc45af, -8, 6.10053603296546099193e-03L), +t9u = LD80C(0xf178d2247adc5093, -9, -3.68456964904901200152e-03L), +t10u = LD80C(0x94188d58f12e5e9f, -9, 2.25976420273774583089e-03L), +t11u = LD80C(0xb7cbaef14e1406f1, -10, -1.40224943666225639823e-03L), +t12u = LD80C(0xe63a671e6704ea4d, -11, 8.78250640744776944887e-04L), +t13u = LD80C(0x914b6c9cae61783e, -11, -5.54255012657716808811e-04L), +t14u = LD80C(0xb858f5bdb79276fe, -12, 3.51614951536825927370e-04L), +t15u = LD80C(0xea73e744c34b9591, -13, -2.23591563824520112236e-04L), +t16u = LD80C(0x99aeabb0d67ba835, -13, 1.46562869351659194136e-04L), +t17u = LD80C(0xd7c6938325db2024, -14, -1.02889866046435680588e-04L), +t18u = LD80C(0xe24cb1e3b0474775, -15, 5.39540265505221957652e-05L); +#define tc (tcu.e) +#define tf (tfu.e) +#define tt (ttu.e) +#define t0 (t0u.e) +#define t1 (t1u.e) +#define t2 (t2u.e) +#define t3 (t3u.e) +#define t4 (t4u.e) +#define t5 (t5u.e) +#define t6 (t6u.e) +#define t7 (t7u.e) +#define t8 (t8u.e) +#define t9 (t9u.e) +#define t10 (t10u.e) +#define t11 (t11u.e) +#define t12 (t12u.e) +#define t13 (t13u.e) +#define t14 (t14u.e) +#define t15 (t15u.e) +#define t16 (t16u.e) +#define t17 (t17u.e) +#define t18 (t18u.e) +/* + * Domain y in [-0.1, 0.232], range ~[-8.1938e-22, 8.3815e-22]: + * |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-71.2 + */ +static const union IEEEl2bits +u0u = LD80C(0x9e233f1bed863d27, -4, -7.72156649015328606095e-02L), +u1u = LD80C(0x98280ee45e4ddd3d, -1, 5.94361239198682739769e-01L), +u2u = LD80C(0xe330c8ead4130733, 0, 1.77492629495841234275e+00L), +u3u = LD80C(0xd4a213f1a002ec52, 0, 1.66119622514818078064e+00L), +u4u = LD80C(0xa5a9ca6f5bc62163, -1, 6.47122051417476492989e-01L), +u5u = LD80C(0xc980e49cd5b019e6, -4, 9.83903751718671509455e-02L), +u6u = LD80C(0xff636a8bdce7025b, -9, 3.89691687802305743450e-03L), +v1u = LD80C(0xbd109c533a19fbf5, 1, 2.95413883330948556544e+00L), +v2u = LD80C(0xd295cbf96f31f099, 1, 3.29039286955665403176e+00L), +v3u = LD80C(0xdab8bcfee40496cb, 0, 1.70876276441416471410e+00L), +v4u = LD80C(0xd2f2dc3638567e9f, -2, 4.12009126299534668571e-01L), +v5u = LD80C(0xa07d9b0851070f41, -5, 3.91822868305682491442e-02L), +v6u = LD80C(0xe3cd8318f7adb2c4, -11, 8.68998648222144351114e-04L); +#define u0 (u0u.e) +#define u1 (u1u.e) +#define u2 (u2u.e) +#define u3 (u3u.e) +#define u4 (u4u.e) +#define u5 (u5u.e) +#define u6 (u6u.e) +#define v1 (v1u.e) +#define v2 (v2u.e) +#define v3 (v3u.e) +#define v4 (v4u.e) +#define v5 (v5u.e) +#define v6 (v6u.e) +/* + * Domain x in (2, 3], range ~[-3.3648e-22, 3.4416e-22]: + * |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-72.3 + * with y = x - 2. + */ +static const union IEEEl2bits +s0u = LD80C(0x9e233f1bed863d27, -4, -7.72156649015328606095e-02L), +s1u = LD80C(0xd3ff0dcc7fa91f94, -3, 2.07027640921219389860e-01L), +s2u = LD80C(0xb2bb62782478ef31, -2, 3.49085881391362090549e-01L), +s3u = LD80C(0xb49f7438c4611a74, -3, 1.76389518704213357954e-01L), +s4u = LD80C(0x9a957008fa27ecf9, -5, 3.77401710862930008071e-02L), +s5u = LD80C(0xda9b389a6ca7a7ac, -9, 3.33566791452943399399e-03L), +s6u = LD80C(0xbc7a2263faf59c14, -14, 8.98728786745638844395e-05L), +r1u = LD80C(0xbf5cff5b11477d4d, 0, 1.49502555796294337722e+00L), +r2u = LD80C(0xd9aec89de08e3da6, -1, 8.50323236984473285866e-01L), +r3u = LD80C(0xeab7ae5057c443f9, -3, 2.29216312078225806131e-01L), +r4u = LD80C(0xf29707d9bd2b1e37, -6, 2.96130326586640089145e-02L), +r5u = LD80C(0xd376c2f09736c5a3, -10, 1.61334161411590662495e-03L), +r6u = LD80C(0xc985983d0cd34e3d, -16, 2.40232770710953450636e-05L), +r7u = LD80C(0xe5c7a4f7fc2ef13d, -25, -5.34997929289167573510e-08L); +#define s0 (s0u.e) +#define s1 (s1u.e) +#define s2 (s2u.e) +#define s3 (s3u.e) +#define s4 (s4u.e) +#define s5 (s5u.e) +#define s6 (s6u.e) +#define r1 (r1u.e) +#define r2 (r2u.e) +#define r3 (r3u.e) +#define r4 (r4u.e) +#define r5 (r5u.e) +#define r6 (r6u.e) +#define r7 (r7u.e) +/* + * Domain z in [8, 0x1p70], range ~[-3.0235e-22, 3.0563e-22]: + * |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-71.7 + */ +static const union IEEEl2bits +w0u = LD80C(0xd67f1c864beb4a69, -2, 4.18938533204672741776e-01L), +w1u = LD80C(0xaaaaaaaaaaaaaaa1, -4, 8.33333333333333332678e-02L), +w2u = LD80C(0xb60b60b60b5491c9, -9, -2.77777777777760927870e-03L), +w3u = LD80C(0xd00d00cf58aede4c, -11, 7.93650793490637233668e-04L), +w4u = LD80C(0x9c09bf626783d4a5, -11, -5.95238023926039051268e-04L), +w5u = LD80C(0xdca7cadc5baa517b, -11, 8.41733700408000822962e-04L), +w6u = LD80C(0xfb060e361e1ffd07, -10, -1.91515849570245136604e-03L), +w7u = LD80C(0xcbd5101bb58d1f2b, -8, 6.22046743903262649294e-03L), +w8u = LD80C(0xad27a668d32c821b, -6, -2.11370706734662081843e-02L); +#define w0 (w0u.e) +#define w1 (w1u.e) +#define w2 (w2u.e) +#define w3 (w3u.e) +#define w4 (w4u.e) +#define w5 (w5u.e) +#define w6 (w6u.e) +#define w7 (w7u.e) +#define w8 (w8u.e) + +static long double +sin_pil(long double x) +{ + volatile long double vz; + long double y,z; + uint64_t n; + uint16_t hx; + + y = -x; + + vz = y+0x1p63L; + z = vz-0x1p63L; + if (z == y) + return zero; + + vz = y+0x1p61; + EXTRACT_LDBL80_WORDS(hx,n,vz); + z = vz-0x1p61; + if (z > y) { + z -= 0.25; /* adjust to round down */ + n--; + } + n &= 7; /* octant of y mod 2 */ + y = y - z + n * 0.25; /* y mod 2 */ + + switch (n) { + case 0: y = __kernel_sinl(pi*y,zero,0); break; + case 1: + case 2: y = __kernel_cosl(pi*(0.5-y),zero); break; + case 3: + case 4: y = __kernel_sinl(pi*(one-y),zero,0); break; + case 5: + case 6: y = -__kernel_cosl(pi*(y-1.5),zero); break; + default: y = __kernel_sinl(pi*(y-2.0),zero,0); break; + } + return -y; +} + +long double +lgammal_r(long double x, int *signgamp) +{ + long double nadj,p,p1,p2,p3,q,r,t,w,y,z; + uint64_t lx; + int i; + uint16_t hx; + + EXTRACT_LDBL80_WORDS(hx,lx,x); + + /* purge off +-inf, NaN, +-0 */ + *signgamp = 1; + if((hx & 0x7fff) == 0x7fff) /* x is +-Inf or NaN */ + return x*x; + if((hx==0||hx==0x8000)&&lx==0) return one/vzero; + + ENTERI(); + + /* purge off tiny and negative arguments */ + if(fabsl(x)<0x1p-70L) { /* |x|<2**-70, return -log(|x|) */ + if(hx&0x8000) { + *signgamp = -1; + RETURNI(-logl(-x)); + } else RETURNI(-logl(x)); + } + if(hx&0x8000) { + if(fabsl(x)>=0x1p63) /* |x|>=2**(p-1), must be -integer */ + RETURNI(one/vzero); + t = sin_pil(x); + if(t==zero) RETURNI(one/vzero); /* -integer */ + nadj = logl(pi/fabsl(t*x)); + if(t=0.7315998077392578) {y = 1-x; i= 0;} + else if(x>=0.2316399812698364) {y= x-(tc-1); i=1;} + else {y = x; i=2;} + } else { + r = 0; + if(x>=1.7316312789916992) {y=2-x;i=0;} + else if(x>=1.2316322326660156) {y=x-tc;i=1;} + else {y=x-1;i=2;} + } + switch(i) { + case 0: + z = y*y; + p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*(a10+z*a12))))); + p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*a13)))))); + p = y*p1+p2; + r += (p-y/2); break; + case 1: + p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+ + y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+ + y*(t17+y*t18)))))))))))))))); + r += (tf + p); break; + case 2: + p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*u6)))))); + p2 = 1+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*v6))))); + r += (-y/2 + p1/p2); + } + } + else if(x<8) { + i = x; + y = x-i; + p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); + q = 1+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*(r6+y*r7)))))); + r = y/2+p/q; + z = 1; /* lgamma(1+s) = log(s) + lgamma(s) */ + switch(i) { + case 7: z *= (y+6); /* FALLTHRU */ + case 6: z *= (y+5); /* FALLTHRU */ + case 5: z *= (y+4); /* FALLTHRU */ + case 4: z *= (y+3); /* FALLTHRU */ + case 3: z *= (y+2); /* FALLTHRU */ + r += logl(z); break; + } + /* 8.0 <= x < 2**70 */ + } else if (x < 0x1p70L) { + t = logl(x); + z = one/x; + y = z*z; + w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*(w6+y*(w7+y*w8))))))); + r = (x-half)*(t-one)+w; + } else + /* 2**70 <= x <= inf */ + r = x*(logl(x)-1); + if(hx&0x8000) r = nadj - r; + RETURNI(r); +} Property changes on: ld80/e_lgammal_r.c ___________________________________________________________________ Added: svn:mime-type ## -0,0 +1 ## +text/plain \ No newline at end of property Added: svn:keywords ## -0,0 +1 ## +FreeBSD=%H \ No newline at end of property Added: svn:eol-style ## -0,0 +1 ## +native \ No newline at end of property Index: man/lgamma.3 =================================================================== --- man/lgamma.3 (revision 271479) +++ man/lgamma.3 (working copy) @@ -28,7 +28,7 @@ .\" from: @(#)lgamma.3 6.6 (Berkeley) 12/3/92 .\" $FreeBSD$ .\" -.Dd January 14, 2005 +.Dd September 12, 2014 .Dt LGAMMA 3 .Os .Sh NAME @@ -36,6 +36,8 @@ .Nm lgamma_r , .Nm lgammaf , .Nm lgammaf_r , +.Nm lgammal , +.Nm lgammal_r , .Nm gamma , .Nm gamma_r , .Nm gammaf , @@ -58,6 +60,10 @@ .Fn lgammaf "float x" .Ft float .Fn lgammaf_r "float x" "int *signgamp" +.Ft "long double" +.Fn lgammal "long double x" +.Ft "long double" +.Fn lgammal_r "long double x" "int *signgamp" .Ft double .Fn gamma "double x" .Ft double @@ -66,14 +72,15 @@ .Fn gammaf "float x" .Ft float .Fn gammaf_r "float x" "int *signgamp" -.Ft double +.Ft "long double" .Fn tgamma "double x" .Ft float .Fn tgammaf "float x" .Sh DESCRIPTION -.Fn lgamma x +.Fn lgamma x , +.Fn lgammaf x , and -.Fn lgammaf x +.Fn lgammal x .if t \{\ return ln\||\(*G(x)| where .Bd -unfilled -offset indent @@ -87,13 +94,15 @@ .Fa signgam returns the sign of \(*G(x). .Pp -.Fn lgamma_r x signgamp +.Fn lgamma_r x signgamp , +.Fn lgammaf_r x signgamp , and -.Fn lgammaf_r x signgamp +.Fn lgammal_r x signgamp provide the same functionality as -.Fn lgamma x +.Fn lgamma x , +.Fn lgammaf x , and -.Fn lgammaf x +.Fn lgammal x , but the caller must provide an integer to store the sign of \(*G(x). .Pp The @@ -115,6 +124,7 @@ and .Fn lgammaf_r , respectively. + .Sh IDIOSYNCRASIES Do not use the expression .Dq Li signgam\(**exp(lgamma(x)) @@ -139,14 +149,18 @@ will lose up to 10 significant bits. .Sh RETURN VALUES .Fn gamma , +.Fn gammaf , +.Fn gammal , .Fn gamma_r , -.Fn gammaf , .Fn gammaf_r , +.Fn gammal_r , .Fn lgamma , +.Fn lgammaf , +.Fn lgammal , .Fn lgamma_r , -.Fn lgammaf , +.Fn lgammaf_r , and -.Fn lgammaf_r +.Fn lgammal_r return appropriate values unless an argument is out of range. Overflow will occur for sufficiently large positive values, and non-positive integers. @@ -159,6 +173,7 @@ The .Fn lgamma , .Fn lgammaf , +.Fn lgammal , .Fn tgamma , and .Fn tgammaf Index: src/e_lgamma.c =================================================================== --- src/e_lgamma.c (revision 271479) +++ src/e_lgamma.c (working copy) @@ -20,6 +20,7 @@ * * Method: call __ieee754_lgamma_r */ +#include #include "math.h" #include "math_private.h" @@ -31,3 +32,7 @@ { return __ieee754_lgamma_r(x,&signgam); } + +#if (LDBL_MANT_DIG == 53) +__weak_reference(lgamma, lgammal); +#endif Index: src/e_lgamma_r.c =================================================================== --- src/e_lgamma_r.c (revision 271479) +++ src/e_lgamma_r.c (working copy) @@ -82,6 +82,7 @@ * lgamma(-inf) = inf (bug for bug compatible with C99!?) * */ +#include #include "math.h" #include "math_private.h" @@ -250,7 +251,7 @@ p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); p = y*p1+p2; - r += (p-0.5*y); break; + r += (p-y/2); break; case 1: z = y*y; w = z*y; @@ -273,11 +274,11 @@ r = half*y+p/q; z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ switch(i) { - case 7: z *= (y+6.0); /* FALLTHRU */ - case 6: z *= (y+5.0); /* FALLTHRU */ - case 5: z *= (y+4.0); /* FALLTHRU */ - case 4: z *= (y+3.0); /* FALLTHRU */ - case 3: z *= (y+2.0); /* FALLTHRU */ + case 7: z *= (y+6); /* FALLTHRU */ + case 6: z *= (y+5); /* FALLTHRU */ + case 5: z *= (y+4); /* FALLTHRU */ + case 4: z *= (y+3); /* FALLTHRU */ + case 3: z *= (y+2); /* FALLTHRU */ r += __ieee754_log(z); break; } /* 8.0 <= x < 2**58 */ @@ -293,3 +294,8 @@ if(hx<0) r = nadj - r; return r; } + +#if (LDBL_MANT_DIG == 53) +__weak_reference(lgamma_r, lgammal_r); +#endif + Index: src/e_lgammaf_r.c =================================================================== --- src/e_lgammaf_r.c (revision 271479) +++ src/e_lgammaf_r.c (working copy) @@ -1,5 +1,6 @@ /* e_lgammaf_r.c -- float version of e_lgamma_r.c. * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + * Conversion to float fixed By Steven G. Kargl. */ /* @@ -22,72 +23,63 @@ static const volatile float vzero = 0; static const float -zero= 0.0000000000e+00, -half= 5.0000000000e-01, /* 0x3f000000 */ -one = 1.0000000000e+00, /* 0x3f800000 */ +zero= 0, +half= 0.5, +one = 1, pi = 3.1415927410e+00, /* 0x40490fdb */ -a0 = 7.7215664089e-02, /* 0x3d9e233f */ -a1 = 3.2246702909e-01, /* 0x3ea51a66 */ -a2 = 6.7352302372e-02, /* 0x3d89f001 */ -a3 = 2.0580807701e-02, /* 0x3ca89915 */ -a4 = 7.3855509982e-03, /* 0x3bf2027e */ -a5 = 2.8905137442e-03, /* 0x3b3d6ec6 */ -a6 = 1.1927076848e-03, /* 0x3a9c54a1 */ -a7 = 5.1006977446e-04, /* 0x3a05b634 */ -a8 = 2.2086278477e-04, /* 0x39679767 */ -a9 = 1.0801156895e-04, /* 0x38e28445 */ -a10 = 2.5214456400e-05, /* 0x37d383a2 */ -a11 = 4.4864096708e-05, /* 0x383c2c75 */ -tc = 1.4616321325e+00, /* 0x3fbb16c3 */ -tf = -1.2148628384e-01, /* 0xbdf8cdcd */ -/* tt = -(tail of tf) */ -tt = 6.6971006518e-09, /* 0x31e61c52 */ -t0 = 4.8383611441e-01, /* 0x3ef7b95e */ -t1 = -1.4758771658e-01, /* 0xbe17213c */ -t2 = 6.4624942839e-02, /* 0x3d845a15 */ -t3 = -3.2788541168e-02, /* 0xbd064d47 */ -t4 = 1.7970675603e-02, /* 0x3c93373d */ -t5 = -1.0314224288e-02, /* 0xbc28fcfe */ -t6 = 6.1005386524e-03, /* 0x3bc7e707 */ -t7 = -3.6845202558e-03, /* 0xbb7177fe */ -t8 = 2.2596477065e-03, /* 0x3b141699 */ -t9 = -1.4034647029e-03, /* 0xbab7f476 */ -t10 = 8.8108185446e-04, /* 0x3a66f867 */ -t11 = -5.3859531181e-04, /* 0xba0d3085 */ -t12 = 3.1563205994e-04, /* 0x39a57b6b */ -t13 = -3.1275415677e-04, /* 0xb9a3f927 */ -t14 = 3.3552918467e-04, /* 0x39afe9f7 */ -u0 = -7.7215664089e-02, /* 0xbd9e233f */ -u1 = 6.3282704353e-01, /* 0x3f2200f4 */ -u2 = 1.4549225569e+00, /* 0x3fba3ae7 */ -u3 = 9.7771751881e-01, /* 0x3f7a4bb2 */ -u4 = 2.2896373272e-01, /* 0x3e6a7578 */ -u5 = 1.3381091878e-02, /* 0x3c5b3c5e */ -v1 = 2.4559779167e+00, /* 0x401d2ebe */ -v2 = 2.1284897327e+00, /* 0x4008392d */ -v3 = 7.6928514242e-01, /* 0x3f44efdf */ -v4 = 1.0422264785e-01, /* 0x3dd572af */ -v5 = 3.2170924824e-03, /* 0x3b52d5db */ -s0 = -7.7215664089e-02, /* 0xbd9e233f */ -s1 = 2.1498242021e-01, /* 0x3e5c245a */ -s2 = 3.2577878237e-01, /* 0x3ea6cc7a */ -s3 = 1.4635047317e-01, /* 0x3e15dce6 */ -s4 = 2.6642270386e-02, /* 0x3cda40e4 */ -s5 = 1.8402845599e-03, /* 0x3af135b4 */ -s6 = 3.1947532989e-05, /* 0x3805ff67 */ -r1 = 1.3920053244e+00, /* 0x3fb22d3b */ -r2 = 7.2193557024e-01, /* 0x3f38d0c5 */ -r3 = 1.7193385959e-01, /* 0x3e300f6e */ -r4 = 1.8645919859e-02, /* 0x3c98bf54 */ -r5 = 7.7794247773e-04, /* 0x3a4beed6 */ -r6 = 7.3266842264e-06, /* 0x36f5d7bd */ -w0 = 4.1893854737e-01, /* 0x3ed67f1d */ -w1 = 8.3333335817e-02, /* 0x3daaaaab */ -w2 = -2.7777778450e-03, /* 0xbb360b61 */ -w3 = 7.9365057172e-04, /* 0x3a500cfd */ -w4 = -5.9518753551e-04, /* 0xba1c065c */ -w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */ -w6 = -1.6309292987e-03; /* 0xbad5c4e8 */ +/* + * Domain y in [0x1p-27, 0.27], range ~[-3.4599e-10, 3.4590e-10]: + * |(lgamma(2 - y) + 0.5 * y) / y - a(y)| < 2**-31.4 + */ +a0 = 7.72156641e-02, /* 0x3d9e233f */ +a1 = 3.22467119e-01, /* 0x3ea51a69 */ +a2 = 6.73484802e-02, /* 0x3d89ee00 */ +a3 = 2.06395667e-02, /* 0x3ca9144f */ +a4 = 6.98275631e-03, /* 0x3be4cf9b */ +a5 = 4.11768444e-03, /* 0x3b86eda4 */ +/* + * Domain x in [tc-0.24, tc+0.28], range ~[-5.6577e-10, 5.5677e-10]: + * |(lgamma(x) - tf) - t(x - tc)| < 2**-30.8. + */ +tc = 1.46163213e+00, /* 0x3fbb16c3 */ +tf = -1.21486291e-01, /* 0xbdf8cdce */ +t0 = -2.94064460e-11, /* 0xae0154b7 */ +t1 = -2.35939837e-08, /* 0xb2caabb8 */ +t2 = 4.83836412e-01, /* 0x3ef7b968 */ +t3 = -1.47586212e-01, /* 0xbe1720d7 */ +t4 = 6.46013096e-02, /* 0x3d844db1 */ +t5 = -3.28450352e-02, /* 0xbd068884 */ +t6 = 1.86483748e-02, /* 0x3c98c47a */ +t7 = -9.89206228e-03, /* 0xbc221251 */ +/* + * Domain y in [-0.1, 0.232], range ~[-8.4931e-10, 8.7794e-10]: + * |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-31.2 + */ +u0 = -7.72156641e-02, /* 0xbd9e233f */ +u1 = 7.36789703e-01, /* 0x3f3c9e40 */ +u2 = 4.95649040e-01, /* 0x3efdc5b6 */ +v1 = 1.10958421e+00, /* 0x3f8e06db */ +v2 = 2.10598111e-01, /* 0x3e57a708 */ +v3 = -1.02995494e-02, /* 0xbc28bf71 */ +/* + * Domain x in (2, 3], range ~[-5.5189e-11, 5.2317e-11]: + * |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-35.0 + * with y = x - 2. + */ +s0 = -7.72156641e-02, /* 0xbd9e233f */ +s1 = 2.69987404e-01, /* 0x3e8a3bca */ +s2 = 1.42851010e-01, /* 0x3e124789 */ +s3 = 1.19389519e-02, /* 0x3c439b98 */ +r1 = 6.79650068e-01, /* 0x3f2dfd8c */ +r2 = 1.16058730e-01, /* 0x3dedb033 */ +r3 = 3.75673687e-03, /* 0x3b763396 */ +/* + * Domain z in [8, 0x1p24], range ~[-1.2640e-09, 1.2640e-09]: + * |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-29.6. + */ +w0 = 4.18938547e-01, /* 0x3ed67f1d */ +w1 = 8.33332464e-02, /* 0x3daaaa9f */ +w2 = -2.76129087e-03; /* 0xbb34f6c6 */ static float sin_pif(float x) @@ -168,7 +160,7 @@ else {y = x; i=2;} } else { r = zero; - if(ix>=0x3fdda618) {y=(float)2.0-x;i=0;} /* [1.7316,2] */ + if(ix>=0x3fdda618) {y=2-x;i=0;} /* [1.7316,2] */ else if(ix>=0x3F9da620) {y=x-tc;i=1;} /* [1.23,1.73] */ else {y=x-one;i=2;} } @@ -175,48 +167,43 @@ switch(i) { case 0: z = y*y; - p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); - p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); + p1 = a0+z*(a2+z*a4); + p2 = z*(a1+z*(a3+z*a5)); p = y*p1+p2; - r += (p-(float)0.5*y); break; + r += (p-y/2); break; case 1: - z = y*y; - w = z*y; - p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */ - p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); - p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); - p = z*p1-(tt-w*(p2+y*p3)); + p = t0+y*t1+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*t7))))); r += (tf + p); break; case 2: - p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); - p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); - r += (-(float)0.5*y + p1/p2); + p1 = y*(u0+y*(u1+y*u2)); + p2 = one+y*(v1+y*(v2+y*v3)); + r += (p1/p2-y/2); } } else if(ix<0x41000000) { /* x < 8.0 */ - i = (int)x; - y = x-(float)i; - p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); - q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); - r = half*y+p/q; + i = x; + y = x-i; + p = y*(s0+y*(s1+y*(s2+y*s3))); + q = one+y*(r1+y*(r2+y*r3)); + r = y/2+p/q; z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ switch(i) { - case 7: z *= (y+(float)6.0); /* FALLTHRU */ - case 6: z *= (y+(float)5.0); /* FALLTHRU */ - case 5: z *= (y+(float)4.0); /* FALLTHRU */ - case 4: z *= (y+(float)3.0); /* FALLTHRU */ - case 3: z *= (y+(float)2.0); /* FALLTHRU */ + case 7: z *= (y+6); /* FALLTHRU */ + case 6: z *= (y+5); /* FALLTHRU */ + case 5: z *= (y+4); /* FALLTHRU */ + case 4: z *= (y+3); /* FALLTHRU */ + case 3: z *= (y+2); /* FALLTHRU */ r += __ieee754_logf(z); break; } - /* 8.0 <= x < 2**58 */ - } else if (ix < 0x5c800000) { + /* 8.0 <= x < 2**24 */ + } else if (ix < 0x4b800000) { t = __ieee754_logf(x); z = one/x; y = z*z; - w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); + w = w0+z*(w1+y*w2); r = (x-half)*(t-one)+w; } else - /* 2**58 <= x <= inf */ + /* 2**24 <= x <= inf */ r = x*(__ieee754_logf(x)-one); if(hx<0) r = nadj - r; return r; Index: src/e_lgammal.c =================================================================== --- src/e_lgammal.c (revision 0) +++ src/e_lgammal.c (working copy) @@ -0,0 +1,27 @@ + +/* @(#)e_lgamma.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +#include +__FBSDID("$FreeBSD$"); + +#include "math.h" +#include "math_private.h" + +extern int signgam; + +long double +lgammal(long double x) +{ + return lgammal_r(x,&signgam); +} Property changes on: src/e_lgammal.c ___________________________________________________________________ Added: svn:mime-type ## -0,0 +1 ## +text/plain \ No newline at end of property Added: svn:keywords ## -0,0 +1 ## +FreeBSD=%H \ No newline at end of property Added: svn:eol-style ## -0,0 +1 ## +native \ No newline at end of property Index: src/imprecise.c =================================================================== --- src/imprecise.c (revision 271479) +++ src/imprecise.c (working copy) @@ -60,5 +60,4 @@ long double imprecise_ ## f ## l(long double v) { return f(v); }\ DECLARE_WEAK(f ## l) -DECLARE_IMPRECISE(lgamma); DECLARE_IMPRECISE(tgamma); Index: src/math.h =================================================================== --- src/math.h (revision 271479) +++ src/math.h (working copy) @@ -496,8 +496,12 @@ long double tanl(long double); long double tgammal(long double); long double truncl(long double); +#endif /* __ISO_C_VISIBLE >= 1999 */ -#endif /* __ISO_C_VISIBLE >= 1999 */ +#if __BSD_VISIBLE +long double lgammal_r(long double, int *); +#endif /* __BSD_VISIBLE */ + __END_DECLS #endif /* !_MATH_H_ */ From owner-freebsd-numerics@FreeBSD.ORG Sat Sep 13 19:51:54 2014 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [8.8.178.115]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id 06DFB391 for ; Sat, 13 Sep 2014 19:51:54 +0000 (UTC) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (Client CN "troutmask.apl.washington.edu", Issuer "troutmask.apl.washington.edu" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id D78F5184 for ; Sat, 13 Sep 2014 19:51:53 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.9/8.14.9) with ESMTP id s8DJprNR002429 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO) for ; Sat, 13 Sep 2014 12:51:53 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.9/8.14.9/Submit) id s8DJprNV002428 for freebsd-numerics@freebsd.org; Sat, 13 Sep 2014 12:51:53 -0700 (PDT) (envelope-from sgk) Date: Sat, 13 Sep 2014 12:51:53 -0700 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: Re: lgammal[_r] patch Message-ID: <20140913195153.GB2323@troutmask.apl.washington.edu> References: <20140913193634.GA2323@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20140913193634.GA2323@troutmask.apl.washington.edu> User-Agent: Mutt/1.5.23 (2014-03-12) X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.18-1 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: Sat, 13 Sep 2014 19:51:54 -0000 On Sat, Sep 13, 2014 at 12:36:34PM -0700, Steve Kargl wrote: > Following my .sig is a patch that contains the ld80 and ld128 > implementations of lgamma[_r](). It also contains an update > to lgammaf_r() and lgamma_r(). Specifically, it has > Yes, a follow-up to my own post. The polynomial approximation for the "t" coefficients has been changed from whatever occurs in lgamma_r() to one that actually works in lgamma[fl]_r(). To explain, I tried to match the approximation in lgamma_r() in the other functions in tc-0.23 <= x <= tc+0.27, but my attempts were fraught with instabilities. Specifically, I tried approximations of the form 1. f(x) = (lgamma(x) - lgamma(tc)) / (x - tc)**2 = t(x - tc) 2. f(x) = (lgamma(x) - lgamma(tc)) = (x - tc)**2 * t(x - tc) 3. f(x) = (lgamma(x) - lgamma(tc)) = t(x - tc) For 1., there were stability issue for x -> tc. For 2., the special Remes algorithm I developed would loose an extrema. I speculate, but haven't pursued, that multiple zeros occurring in the rhs at x = tc is the root of the problem. For 3., I had sensitive issues near the endpoints of the domain. So, I adjusted the domain to [tc-0.24, tc+0.28]. This yielded the coefficients found in the code. -- steve