long double versions of functions in gcc under Cygwin

Discussion in 'C Programming' started by lcw1964, Aug 8, 2006.

  1. lcw1964

    lcw1964 Guest

    Greetings, all,

    I am trying to port a little bit of math code to gcc, that in the
    original version used the long double version of several functions (in
    particular, atanl, fabsl, and expl).

    I get a complie-time "unidentified reference" error to the expl()
    calls, but gcc seems to digest atanl and fabsl just fine. Changing expl
    to exp cures the compile time problem, but I get at best double
    precision in the final results. I am assuming that the use of exp() vs.
    expl() is the weak link.

    The GCC documentation seems to imply that expl() is supported, but I
    have no idea where to find it or how to link it in properly. For that
    matter, I can't seem to find prototypes in math.h for fabsl or atanl,
    and they don't make gcc cough at all.

    I hope this tenderfoot can find some direction, or I may resort to
    singing the praises of the egregiously un-portable lcc-win32 with its
    impressive 100+ digit precision qfloat library ;)

    cheers,

    Les

    p.s. I am trying to keep this simple, so if there is a solution within
    the main gcc offerings without me having to turn to the GSL, I would
    like to try that first.
    lcw1964, Aug 8, 2006
    #1
    1. Advertising

  2. "lcw1964" <> writes:
    > I am trying to port a little bit of math code to gcc, that in the
    > original version used the long double version of several functions (in
    > particular, atanl, fabsl, and expl).
    >
    > I get a complie-time "unidentified reference" error to the expl()
    > calls, but gcc seems to digest atanl and fabsl just fine. Changing expl
    > to exp cures the compile time problem, but I get at best double
    > precision in the final results. I am assuming that the use of exp() vs.
    > expl() is the weak link.
    >
    > The GCC documentation seems to imply that expl() is supported, but I
    > have no idea where to find it or how to link it in properly. For that
    > matter, I can't seem to find prototypes in math.h for fabsl or atanl,
    > and they don't make gcc cough at all.
    >
    > I hope this tenderfoot can find some direction, or I may resort to
    > singing the praises of the egregiously un-portable lcc-win32 with its
    > impressive 100+ digit precision qfloat library ;)
    >
    > cheers,
    >
    > Les
    >
    > p.s. I am trying to keep this simple, so if there is a solution within
    > the main gcc offerings without me having to turn to the GSL, I would
    > like to try that first.


    gcc is a compiler, not a complete C implementation. (Actually gcc is
    a collection of compilers, but for our purposes here we can consider
    only the C compiler.) The math functions are implemented by the
    runtime library, not by the compiler.

    In some implementations, the compiler and the runtime library are
    provided together. gcc, however, generally uses whatever runtime
    library is provided by the underlying operating system. On some
    systems, the C runtime library happens to be one that, like gcc, is
    also provided by the GNU project. I suspect you're using a system
    where that isn't the case.

    I suggest you ask in a newsgroup that deals with your operating system
    (probably MS Windows given your mention of lcc-win32 as an
    alternative).

    (I don't know whether lcc-win32 provides its own C runtime library;
    check the web site or ask in comp.compilers.lcc if you want more
    information.)

    --
    Keith Thompson (The_Other_Keith) <http://www.ghoti.net/~kst>
    San Diego Supercomputer Center <*> <http://users.sdsc.edu/~kst>
    We must do something. This is something. Therefore, we must do this.
    Keith Thompson, Aug 8, 2006
    #2
    1. Advertising

  3. lcw1964

    lcw1964 Guest

    Keith Thompson wrote:
    >
    > In some implementations, the compiler and the runtime library are
    > provided together. gcc, however, generally uses whatever runtime
    > library is provided by the underlying operating system. On some
    > systems, the C runtime library happens to be one that, like gcc, is
    > also provided by the GNU project. I suspect you're using a system
    > where that isn't the case.
    >



    I think you are right, and I was afraid it had something to do with
    that.

    I am learning quickly that if I wish to make the most out of higher
    precision math programming, it behooves me to expand my horizons and
    develop some facility with a high or even arbitrary precision package,
    though I must admit that some of the source code I have contemplated
    looks daunting indeed for this beginner.

    In the meantime, lcc-win32 seems a reasonable, though admittedly
    non-portable option. Messr. Navia's implementation of the Cephes qfloat
    library seems pretty robust and at least in my barely ept hands
    produces prodigious results with surprising little change to code. For
    example, the very code that I am having trouble getting full long
    double results with gcc/Cygwin with lcc-win32/qfloat routinely gives
    102-105 digit accuracy most of the time, and at least 100 digits all of
    the time. It is simply a matter of including qfloat.h, change various
    commands to their qfloat equivalents (atoq, expq, atanq, etc.),
    appending a q to floating point constants, and properly formatting the
    output strings for printf or whatever. For my limited personal
    purposes, it is about the best option I have hit upon so far, though I
    do admit that cross-platform and cross-compiler compatibility would be
    much more vital if my interests were less parochial.

    Thanks for the feedback, though I must admit it leaves me with 3.1 gigs
    of Cygwin on my hard drive that I don't know quite what to do with ;)

    Les
    lcw1964, Aug 8, 2006
    #3
  4. lcw1964

    Dann Corbit Guest

    "lcw1964" <> wrote in message
    news:...
    > Greetings, all,
    >
    > I am trying to port a little bit of math code to gcc, that in the
    > original version used the long double version of several functions (in
    > particular, atanl, fabsl, and expl).
    >
    > I get a complie-time "unidentified reference" error to the expl()
    > calls, but gcc seems to digest atanl and fabsl just fine. Changing expl
    > to exp cures the compile time problem, but I get at best double
    > precision in the final results. I am assuming that the use of exp() vs.
    > expl() is the weak link.
    >
    > The GCC documentation seems to imply that expl() is supported, but I
    > have no idea where to find it or how to link it in properly. For that
    > matter, I can't seem to find prototypes in math.h for fabsl or atanl,
    > and they don't make gcc cough at all.
    >
    > I hope this tenderfoot can find some direction, or I may resort to
    > singing the praises of the egregiously un-portable lcc-win32 with its
    > impressive 100+ digit precision qfloat library ;)


    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

    > cheers,
    >
    > Les
    >
    > p.s. I am trying to keep this simple, so if there is a solution within
    > the main gcc offerings without me having to turn to the GSL, I would
    > like to try that first.
    >
    Dann Corbit, Aug 8, 2006
    #4
  5. lcw1964

    lcw1964 Guest

    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.

    Les
    lcw1964, Aug 8, 2006
    #5
  6. lcw1964

    Dann Corbit Guest

    "lcw1964" <> wrote in message
    news:...
    >
    > 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).
    Dann Corbit, Aug 8, 2006
    #6
  7. lcw1964

    jaysome Guest

    On Mon, 7 Aug 2006 21:11:08 -0700, "Dann Corbit" <>
    wrote:

    >"lcw1964" <> wrote in message
    >news:...
    >>
    >> 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.

    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.

    --
    Jay
    jaysome, Aug 8, 2006
    #7
  8. lcw1964

    jacob navia Guest

    jaysome a écrit :
    > There are 477 usages of "goto" in this source. That gives me a queasy
    > feeling.
    >


    I rewrote all the basic functions in 386 and AMD64 assembly.
    The speed gain is considerable, and the gotos are even worst:

    Who hasn't written a

    jmp label

    in assembly?

    Seriously, the code is well written, and if you look at the
    dates in there you will se code from eighties. And it still runs,
    twenty years later.

    I would like to see what code you have written in 20 years, even
    if it doesn't use gotos.

    Stephen Moshier has written a very good package.


    > One part of me sees me sitting through a code review and vehemently
    > rebuking this code after coming across about the 5th goto statement.


    This is just dogmatic. gotos arre part of C. And they are used in
    the Cephes library in a reasonable and very clear way.

    > 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.
    >
    jacob navia, Aug 8, 2006
    #8
  9. lcw1964

    pete Guest

    lcw1964 wrote:

    > expl()


    #include <float.h>

    long double fs_expl(long double x);
    long double fs_logl(long double x);
    long double fs_sqrtl(long double x);

    long double fs_expl(long double x)
    {
    long unsigned n, square;
    long double b, e;
    static long double x_max, x_min;

    if (1 > x_max) {
    x_max = fs_logl(LDBL_MAX);
    x_min = fs_logl(LDBL_MIN);
    }
    if (x_max >= x && x >= x_min) {
    for (square = 0; x > 1; x /= 2) {
    ++square;
    }
    while (-1 > x) {
    ++square;
    x /= 2;
    }
    e = b = n = 1;
    do {
    b /= n++;
    b *= x;
    e += b;
    b /= n++;
    b *= x;
    e += b;
    } while (b > LDBL_EPSILON / 4);
    while (square-- != 0) {
    e *= e;
    }
    } else {
    e = x > 0 ? LDBL_MAX : 0;
    }
    return e;
    }

    long double fs_logl(long double x)
    {
    long int n;
    long double a, b, c, epsilon;
    static long double A, B, C;

    if (LDBL_MAX >= x && x > 0) {
    if (1 > A) {
    A = fs_sqrtl(2);
    B = A / 2;
    C = fs_logl(A);
    }
    for (n = 0; x > A; x /= 2) {
    ++n;
    }
    while (B > x) {
    --n;
    x *= 2;
    }
    a = (x - 1) / (x + 1);
    x = C * n + a;
    c = a * a;
    n = 1;
    epsilon = LDBL_EPSILON * x;
    if (0 > a) {
    if (epsilon > 0) {
    epsilon = -epsilon;
    }
    do {
    n += 2;
    a *= c;
    b = a / n;
    x += b;
    } while (epsilon > b);
    } else {
    if (0 > epsilon) {
    epsilon = -epsilon;
    }
    do {
    n += 2;
    a *= c;
    b = a / n;
    x += b;
    } while (b > epsilon);
    }
    x *= 2;
    } else {
    x = -LDBL_MAX;
    }
    return x;
    }

    long double fs_sqrtl(long double x)
    {
    long int n;
    long double a, b;

    if (LDBL_MAX >= x && x > 0) {
    for (n = 0; x > 2; x /= 4) {
    ++n;
    }
    while (0.5 > x) {
    --n;
    x *= 4;
    }
    a = x;
    b = (1 + x) / 2;
    do {
    x = b;
    b = (a / x + x) / 2;
    } while (x > b);
    while (n > 0) {
    x *= 2;
    --n;
    }
    while (0 > n) {
    x /= 2;
    ++n;
    }
    } else {
    if (x != 0) {
    x = LDBL_MAX;
    }
    }
    return x;
    }

    --
    pete
    pete, Aug 8, 2006
    #9
  10. lcw1964

    jacob navia Guest

    pete a écrit :
    > lcw1964 wrote:
    >
    >
    >>expl()

    >
    >
    > #include <float.h>
    >
    > long double fs_expl(long double x);
    > long double fs_logl(long double x);
    > long double fs_sqrtl(long double x);
    >
    > long double fs_expl(long double x)
    > {
    > long unsigned n, square;
    > long double b, e;
    > static long double x_max, x_min;
    >
    > if (1 > x_max) {
    > x_max = fs_logl(LDBL_MAX);
    > x_min = fs_logl(LDBL_MIN);
    > }
    > if (x_max >= x && x >= x_min) {
    > for (square = 0; x > 1; x /= 2) {
    > ++square;
    > }
    > while (-1 > x) {
    > ++square;
    > x /= 2;
    > }
    > e = b = n = 1;
    > do {
    > b /= n++;
    > b *= x;
    > e += b;
    > b /= n++;
    > b *= x;
    > e += b;
    > } while (b > LDBL_EPSILON / 4);
    > while (square-- != 0) {
    > e *= e;
    > }
    > } else {
    > e = x > 0 ? LDBL_MAX : 0;
    > }
    > return e;
    > }
    >
    > long double fs_logl(long double x)
    > {
    > long int n;
    > long double a, b, c, epsilon;
    > static long double A, B, C;
    >
    > if (LDBL_MAX >= x && x > 0) {
    > if (1 > A) {
    > A = fs_sqrtl(2);
    > B = A / 2;
    > C = fs_logl(A);
    > }
    > for (n = 0; x > A; x /= 2) {
    > ++n;
    > }
    > while (B > x) {
    > --n;
    > x *= 2;
    > }
    > a = (x - 1) / (x + 1);
    > x = C * n + a;
    > c = a * a;
    > n = 1;
    > epsilon = LDBL_EPSILON * x;
    > if (0 > a) {
    > if (epsilon > 0) {
    > epsilon = -epsilon;
    > }
    > do {
    > n += 2;
    > a *= c;
    > b = a / n;
    > x += b;
    > } while (epsilon > b);
    > } else {
    > if (0 > epsilon) {
    > epsilon = -epsilon;
    > }
    > do {
    > n += 2;
    > a *= c;
    > b = a / n;
    > x += b;
    > } while (b > epsilon);
    > }
    > x *= 2;
    > } else {
    > x = -LDBL_MAX;
    > }
    > return x;
    > }
    >
    > long double fs_sqrtl(long double x)
    > {
    > long int n;
    > long double a, b;
    >
    > if (LDBL_MAX >= x && x > 0) {
    > for (n = 0; x > 2; x /= 4) {
    > ++n;
    > }
    > while (0.5 > x) {
    > --n;
    > x *= 4;
    > }
    > a = x;
    > b = (1 + x) / 2;
    > do {
    > x = b;
    > b = (a / x + x) / 2;
    > } while (x > b);
    > while (n > 0) {
    > x *= 2;
    > --n;
    > }
    > while (0 > n) {
    > x /= 2;
    > ++n;
    > }
    > } else {
    > if (x != 0) {
    > x = LDBL_MAX;
    > }
    > }
    > return x;
    > }
    >


    I find this code well DOCUMENTED isn't it?

    The source of the code (who wrote it originally), the algorithms
    used are well explained, the places in the code where you have to
    watch for accuracy are pointed out, a nice package.

    fs_sqrt is approximately 20 times slower
    than the library function.
    jacob navia, Aug 8, 2006
    #10
  11. lcw1964

    pete Guest

    jacob navia wrote:
    >
    > pete a écrit :


    > > long double fs_sqrtl(long double x);


    > fs_sqrt is approximately 20 times slower
    > than the library function.


    OP didn't say if sqrtl was available.

    --
    pete
    pete, Aug 8, 2006
    #11
  12. lcw1964

    lcw1964 Guest

    jacob navia wrote:
    >
    > I find this code well DOCUMENTED isn't it?
    >
    > The source of the code (who wrote it originally), the algorithms
    > used are well explained, the places in the code where you have to
    > watch for accuracy are pointed out, a nice package.
    >
    > fs_sqrt is approximately 20 times slower
    > than the library function.


    i appreciate this code being shared. I don't want to get get caught in
    the middle of barbed repartee, but I have to admit some editorial
    comments would have been helpful, especially for a neophyte like me.

    i must admit i am a little perplexed that my initial question, which i
    thought was straightforward, was not answered--namely, when I use gcc
    under Cygwin and try to call long double expl(long double), a function
    that works perfectly well when I compile under Borland C++ or
    lcc-win32, I get an unknown identifier error, whereas fabsl and sqrtl
    work fine. Forgive my naivete, but I thought that expl() was a standard
    C function, and if I am including the wrong files or should be linking
    in something else, I would like to know.

    if i have already been given my answer please forgive me for missing
    the point and I will try to reread what has been offered and change my
    ways.

    many thanks,

    Les
    lcw1964, Aug 8, 2006
    #12
  13. In article <>,
    lcw1964 <> wrote:
    ....
    >i must admit i am a little perplexed that my initial question, which i
    >thought was straightforward, was not answered--namely, when I use gcc
    >under Cygwin and try to call long double expl(long double), a function
    >that works perfectly well when I compile under Borland C++ or
    >lcc-win32, I get an unknown identifier error, whereas fabsl and sqrtl


    Anything that involves brand names is OT here.
    Kenny McCormack, Aug 8, 2006
    #13
  14. lcw1964

    Richard Guest

    (Kenny McCormack) writes:

    > In article <>,
    > lcw1964 <> wrote:
    > ...
    >>i must admit i am a little perplexed that my initial question, which i
    >>thought was straightforward, was not answered--namely, when I use gcc
    >>under Cygwin and try to call long double expl(long double), a function
    >>that works perfectly well when I compile under Borland C++ or
    >>lcc-win32, I get an unknown identifier error, whereas fabsl and sqrtl

    >
    > Anything that involves brand names is OT here.
    >


    Even if a "standards" solution solves his issues? Frequently using the
    correct way fixes issues with the "platorm special" way - its why many
    people ask C questions here I would have thought.
    Richard, Aug 8, 2006
    #14
  15. Troll Alert: Kenny McCormack

    The only way to deal with trolls is to limit your reaction to reminding
    others not to respond to trolls.

    Information on trolls: http://members.aol.com/intwg/trolls.htm

    --

    Frederick Gotham
    Frederick Gotham, Aug 8, 2006
    #15
  16. lcw1964

    jacob navia Guest

    lcw1964 wrote:

    > i must admit i am a little perplexed that my initial question, which i
    > thought was straightforward, was not answered--namely, when I use gcc
    > under Cygwin and try to call long double expl(long double), a function
    > that works perfectly well when I compile under Borland C++ or
    > lcc-win32, I get an unknown identifier error, whereas fabsl and sqrtl
    > work fine. Forgive my naivete, but I thought that expl() was a standard
    > C function, and if I am including the wrong files or should be linking
    > in something else, I would like to know.
    >
    > if i have already been given my answer please forgive me for missing
    > the point and I will try to reread what has been offered and change my
    > ways.
    >
    > many thanks,
    >
    > Les
    >


    Obviously this is a bug. You could report it to them, maybe they
    are interested in knowing about it. There must be some mailing
    list in the cygwin docs. I was subscribed ages ago.

    Under linux:
    [root@gateway tmp]# cat texpl.c
    #include <stdio.h>
    #include <math.h>
    int main(void)
    {
    long double n = expl(1.0L);
    printf("%Lg\n",n);
    return 0;
    }
    [root@gateway tmp]# gcc texpl.c -lm
    [root@gateway tmp]# ./a.out
    2.71828
    [root@gateway tmp]#

    this works, so it must be a bug in the cygwin environment/library.
    jacob navia, Aug 8, 2006
    #16
  17. lcw1964

    Dann Corbit Guest

    "jaysome" <> wrote in message
    news:...
    > On Mon, 7 Aug 2006 21:11:08 -0700, "Dann Corbit" <>
    > wrote:
    >
    >>"lcw1964" <> wrote in message
    >>news:...
    >>>
    >>> 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
    >
    Dann Corbit, Aug 8, 2006
    #17
  18. lcw1964

    lcw1964 Guest

    jacob navia wrote:
    > lcw1964 wrote:
    >
    > > i must admit i am a little perplexed that my initial question, which i
    > > thought was straightforward, was not answered--namely, when I use gcc
    > > under Cygwin and try to call long double expl(long double), a function
    > > that works perfectly well when I compile under Borland C++ or
    > > lcc-win32, I get an unknown identifier error, whereas fabsl and sqrtl
    > > work fine. Forgive my naivete, but I thought that expl() was a standard
    > > C function, and if I am including the wrong files or should be linking
    > > in something else, I would like to know.
    > >
    > > if i have already been given my answer please forgive me for missing
    > > the point and I will try to reread what has been offered and change my
    > > ways.
    > >
    > > many thanks,
    > >
    > > Les
    > >

    >
    > Obviously this is a bug. You could report it to them, maybe they
    > are interested in knowing about it. There must be some mailing
    > list in the cygwin docs. I was subscribed ages ago.
    >
    > Under linux:
    > [root@gateway tmp]# cat texpl.c
    > #include <stdio.h>
    > #include <math.h>
    > int main(void)
    > {
    > long double n = expl(1.0L);
    > printf("%Lg\n",n);
    > return 0;
    > }
    > [root@gateway tmp]# gcc texpl.c -lm
    > [root@gateway tmp]# ./a.out
    > 2.71828
    > [root@gateway tmp]#
    >
    > this works, so it must be a bug in the cygwin environment/library.



    Messr. Navia, that worked for me.

    The key is the addition of the -lm parameter to the command line, which
    I did not have before.

    I suspect that I have made an embarassing beginner's error and have
    stirred up a lot of hubbub unnecessarily!

    I have no idea what those three little characters mean (-lm), but
    something tells me that before I post my next newbie question I do a
    little more homework first.

    Les
    lcw1964, Aug 8, 2006
    #18
  19. "lcw1964" <> writes:
    [...]
    > i must admit i am a little perplexed that my initial question, which i
    > thought was straightforward, was not answered--namely, when I use gcc
    > under Cygwin and try to call long double expl(long double), a function
    > that works perfectly well when I compile under Borland C++ or
    > lcc-win32, I get an unknown identifier error, whereas fabsl and sqrtl
    > work fine. Forgive my naivete, but I thought that expl() was a standard
    > C function, and if I am including the wrong files or should be linking
    > in something else, I would like to know.


    expl() is defined by the C99 standard, but not by the older C90
    standard. Since C99 is not yet implemented as widely as C90, a given
    implementation might not support expl(). (C90 has exp() which takes
    and returns a double; C99 adds expf() and expl().)

    You might find that exp() is good enough. If not, you'll have to find
    another solution.

    For more information, try the Cygwin mailing list; see www.cygwin.com
    for links.

    --
    Keith Thompson (The_Other_Keith) <http://www.ghoti.net/~kst>
    San Diego Supercomputer Center <*> <http://users.sdsc.edu/~kst>
    We must do something. This is something. Therefore, we must do this.
    Keith Thompson, Aug 8, 2006
    #19
  20. lcw1964

    lcw1964 Guest

    lcw1964 wrote:
    > jacob navia wrote:
    > > lcw1964 wrote:
    > >
    > > > i must admit i am a little perplexed that my initial question, which i
    > > > thought was straightforward, was not answered--namely, when I use gcc
    > > > under Cygwin and try to call long double expl(long double), a function
    > > > that works perfectly well when I compile under Borland C++ or
    > > > lcc-win32, I get an unknown identifier error, whereas fabsl and sqrtl
    > > > work fine. Forgive my naivete, but I thought that expl() was a standard
    > > > C function, and if I am including the wrong files or should be linking
    > > > in something else, I would like to know.
    > > >
    > > > if i have already been given my answer please forgive me for missing
    > > > the point and I will try to reread what has been offered and change my
    > > > ways.
    > > >
    > > > many thanks,
    > > >
    > > > Les
    > > >

    > >
    > > Obviously this is a bug. You could report it to them, maybe they
    > > are interested in knowing about it. There must be some mailing
    > > list in the cygwin docs. I was subscribed ages ago.
    > >
    > > Under linux:
    > > [root@gateway tmp]# cat texpl.c
    > > #include <stdio.h>
    > > #include <math.h>
    > > int main(void)
    > > {
    > > long double n = expl(1.0L);
    > > printf("%Lg\n",n);
    > > return 0;
    > > }
    > > [root@gateway tmp]# gcc texpl.c -lm
    > > [root@gateway tmp]# ./a.out
    > > 2.71828
    > > [root@gateway tmp]#
    > >
    > > this works, so it must be a bug in the cygwin environment/library.

    >
    >
    > Messr. Navia, that worked for me.
    >
    > The key is the addition of the -lm parameter to the command line, which
    > I did not have before.
    >
    > I suspect that I have made an embarassing beginner's error and have
    > stirred up a lot of hubbub unnecessarily!
    >
    > I have no idea what those three little characters mean (-lm), but
    > something tells me that before I post my next newbie question I do a
    > little more homework first.
    >


    Actually I spoke too soon!

    The following variant of M. Navia's example generates the error too:

    #include <stdio.h>
    #include <math.h>
    int main(void)
    { long double x;
    x = 1.0L;
    long double n = expl(x);
    printf("%.20Lg\n",n);
    return 0;
    }

    So too does this variant, where the argument is other than 1.0L:

    #include <stdio.h>
    #include <math.h>
    int main(void)
    {
    long double n = expl(3.789L);
    printf("%Lg\n",n);
    return 0;
    }

    Now I am no C god, and I present myself to you as a lowly plebeian as a
    superstitious emperor presents himself to the Oracle of Delphi, but
    something tells me that the compiler recognizes the string "expl(1.0L)"
    as the constant 2.7182818284590452354, NOT as a call to the long
    double version of exp. Try to put anything else between those
    parentheses except 1.0L, 1.0, 1--a long double variable or some other
    long double constant other than unity--and the compiler regurgitates it
    back contemptuously.

    It could be a bug in my platform or setup, but if M. Navia or someone
    else could experiment with this in gcc under Linux or any other
    platform, I would be much obliged. This is driving me bonkers now!!!!!

    Les
    lcw1964, Aug 8, 2006
    #20
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Sydex
    Replies:
    12
    Views:
    6,453
    Victor Bazarov
    Feb 17, 2005
  2. ing. Jan Chládek

    wstring/wcout in GCC under DJGPP/Cygwin

    ing. Jan Chládek, Apr 24, 2005, in forum: C++
    Replies:
    7
    Views:
    5,370
    ing. Jan Chládek
    Apr 26, 2005
  3. Daniel Rudy

    unsigned long long int to long double

    Daniel Rudy, Sep 19, 2005, in forum: C Programming
    Replies:
    5
    Views:
    1,178
    Peter Shaggy Haywood
    Sep 20, 2005
  4. lcw1964

    long double in gcc implementations

    lcw1964, Jul 29, 2006, in forum: C Programming
    Replies:
    67
    Views:
    1,323
  5. Francesco S. Carta
    Replies:
    6
    Views:
    675
    Francesco S. Carta
    Sep 20, 2009
Loading...

Share This Page