Re: long double versions of functions in gcc under Cygwin

"jaysome" <jaysome@xxxxxxxxxxx> wrote in message
news:58ggd29fb960dal5thfnqgj7dagndnpddk@xxxxxxxxxx
On Mon, 7 Aug 2006 21:11:08 -0700, "Dann Corbit" <dcorbit@xxxxxxxxx>
wrote:

"lcw1964" <leslie.wright@xxxxxxxxxxxxx> wrote in message
news:1155007718.733220.249820@xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Dann Corbit wrote:

The qfloat library is by S. Moshier. You can find qfloat along with
the
Cephes collection (which has tons of long double math functions) here:

http://www.moshier.net/#Cephes


Thank you! I should have given Mr. Moshier proper credit. I am also
aware of the link you referred me to, my right now porting those
libraries to GCC is a little beyond my skill set, so being able to
access the qfloat functionality thru Mr. Navia's lcc-win32 "wrapper" is
a good start.

You don't have to know anything. They come with their own makefiles.

At most, you will have to know what kind of machine you are compiling on
(if
it is not a 32 bit platform or has odd endianness or something).

The standard makefile will probably fit your situation.

Just expand this archive:
http://www.moshier.net/qlib.zip
and type "make"

The Cephes functions are even the default math functions used in some
linux
distributions (IIRC).

There are 477 usages of "goto" in this source. That gives me a queasy
feeling.

I'm in the Knuth camp:
http://portal.acm.org/citation.cfm?id=356640&dl=&coll=GUIDE&CFID=15151515&CFTOKEN=6184618


One part of me sees me sitting through a code review and vehemently
rebuking this code after coming across about the 5th goto statement.
The other part of me sees me reviewing the black box test results that
passed and not caring about how this was coded, as long as it was
coded in Standard C. Oh the dichotomy.

I think Moshier's code is beautiful. It was originally written in 1984, and
later reworked so that you could compile it with ANSI style prototypes.

At the top of every file is a detailed explanation of what the code does,
which equations are used to solve the problem, and what the measured errors
were in calibration of the function's useful range.

Here is a sample, complete with goto (which could be removed, of course):

/* psi.c

* Psi (digamma) function
*
*
* SYNOPSIS:
*
* double x, y, psi();
*
* y = psi( x );
*
*
* DESCRIPTION:
*
* d -
* psi(x) = -- ln | (x)
* dx
*
* is the logarithmic derivative of the gamma function.
* For integer x,
* n-1
* -
* psi(n) = -EUL + > 1/k.
* -
* k=1
*
* This formula is used for 0 < n <= 10. If x is negative, it
* is transformed to a positive argument by the reflection
* formula psi(1-x) = psi(x) + pi cot(pi x).
* For general positive x, the argument is made greater than 10
* using the recurrence psi(x+1) = psi(x) + 1/x.
* Then the following asymptotic expansion is applied:
*
* inf. B
* - 2k
* psi(x) = log(x) - 1/2x - > -------
* - 2k
* k=1 2k x
*
* where the B2k are Bernoulli numbers.
*
* ACCURACY:
* Relative error (except absolute when |psi| < 1):
* arithmetic domain # trials peak rms
* DEC 0,30 2500 1.7e-16 2.0e-17
* IEEE 0,30 30000 1.3e-15 1.4e-16
* IEEE -30,0 40000 1.5e-15 2.2e-16
*
* ERROR MESSAGES:
* message condition value returned
* psi singularity x integer <=0 MAXNUM
*/


/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier
*/

#include "mconf.h"

#ifdef UNK
static double A[] =
{
8.33333333333333333333E-2,
-2.10927960927960927961E-2,
7.57575757575757575758E-3,
-4.16666666666666666667E-3,
3.96825396825396825397E-3,
-8.33333333333333333333E-3,
8.33333333333333333333E-2
};
#endif

#ifdef DEC
static unsigned short A[] =
{
0037252, 0125252, 0125252, 0125253,
0136654, 0145314, 0126312, 0146255,
0036370, 0037017, 0101740, 0174076,
0136210, 0104210, 0104210, 0104211,
0036202, 0004040, 0101010, 0020202,
0136410, 0104210, 0104210, 0104211,
0037252, 0125252, 0125252, 0125253
};
#endif

#ifdef IBMPC
static unsigned short A[] =
{
0x5555, 0x5555, 0x5555, 0x3fb5,
0x5996, 0x9599, 0x9959, 0xbf95,
0x1f08, 0xf07c, 0x07c1, 0x3f7f,
0x1111, 0x1111, 0x1111, 0xbf71,
0x0410, 0x1041, 0x4104, 0x3f70,
0x1111, 0x1111, 0x1111, 0xbf81,
0x5555, 0x5555, 0x5555, 0x3fb5
};
#endif

#ifdef MIEEE
static unsigned short A[] =
{
0x3fb5, 0x5555, 0x5555, 0x5555,
0xbf95, 0x9959, 0x9599, 0x5996,
0x3f7f, 0x07c1, 0xf07c, 0x1f08,
0xbf71, 0x1111, 0x1111, 0x1111,
0x3f70, 0x4104, 0x1041, 0x0410,
0xbf81, 0x1111, 0x1111, 0x1111,
0x3fb5, 0x5555, 0x5555, 0x5555
};
#endif

#define EUL 0.57721566490153286061

#ifdef ANSIPROT
extern double floor(double);
extern double log(double);
extern double tan(double);
extern double polevl(double, void *, int);
#else
double floor(), log(), tan(), polevl();
#endif
extern double PI,
MAXNUM;


double psi(x)
double x;
{
double p,
q,
nz,
s,
w,
y,
z;
int i,
n,
negative;

negative = 0;
nz = 0.0;

if (x <= 0.0) {
negative = 1;
q = x;
p = floor(q);
if (p == q) {
mtherr("psi", SING);
return (MAXNUM);
}
/* Remove the zeros of tan(PI x)
* by subtracting the nearest integer from x
*/
nz = q - p;
if (nz != 0.5) {
if (nz > 0.5) {
p += 1.0;
nz = q - p;
}
nz = PI / tan(PI * nz);
} else {
nz = 0.0;
}
x = 1.0 - x;
}
/* check for positive integer up to 10 */
if ((x <= 10.0) && (x == floor(x))) {
y = 0.0;
n = x;
for (i = 1; i < n; i++) {
w = i;
y += 1.0 / w;
}
y -= EUL;
goto done;
}
s = x;
w = 0.0;
while (s < 10.0) {
w += 1.0 / s;
s += 1.0;
}

if (s < 1.0e17) {
z = 1.0 / (s * s);
y = z * polevl(z, A, 6);
} else
y = 0.0;

y = log(s) - (0.5 / s) - y - w;

done:

if (negative) {
y -= nz;
}
return (y);
}

--
Jay



.

Relevant Pages