sqrtl(), etc. for Cygwin

Discussion in 'C Programming' started by James Dow Allen, Nov 7, 2012.

  1. I've been using Cygwin lately, and am largely satisfied.
    It has gcc, but lacks sqrtl(), atan2l(), etc.

    What's the recommended remedy?
    I Googled "source sqrtl" and got to codeforge.com
    but only AFTER passing one painful Gotcha test
    was I told I'd need another password.

    I'd bother with making a free(?) account if I were sure
    codeforge.com is the best piece of cheese in the maze,
    but it was just a random Google hit and
    it seemed wiser to ask the experts here.
    (I downloaded Cygwin just a few months ago --
    why doesn't its gcc have sqrtl?)

    James
    James Dow Allen, Nov 7, 2012
    #1
    1. Advertising

  2. James Dow Allen

    Eric Sosman Guest

    On 11/7/2012 1:40 PM, James Dow Allen wrote:
    > I've been using Cygwin lately, and am largely satisfied.
    > It has gcc, but lacks sqrtl(), atan2l(), etc.
    >
    > What's the recommended remedy?
    > I Googled "source sqrtl" and got to codeforge.com
    > but only AFTER passing one painful Gotcha test
    > was I told I'd need another password.
    >
    > I'd bother with making a free(?) account if I were sure
    > codeforge.com is the best piece of cheese in the maze,
    > but it was just a random Google hit and
    > it seemed wiser to ask the experts here.
    > (I downloaded Cygwin just a few months ago --
    > why doesn't its gcc have sqrtl?)


    A web search for "Cygwin long double" suggests that this
    lack is a long-standing issue, possibly with roots in licensing
    land. Maybe if you grovel through the search results more
    thoroughly than I did, you'll find something more helpful ...

    --
    Eric Sosman
    d
    Eric Sosman, Nov 7, 2012
    #2
    1. Advertising

  3. James Dow Allen <> writes:
    > I've been using Cygwin lately, and am largely satisfied.
    > It has gcc, but lacks sqrtl(), atan2l(), etc.
    >
    > What's the recommended remedy?
    > I Googled "source sqrtl" and got to codeforge.com
    > but only AFTER passing one painful Gotcha test
    > was I told I'd need another password.
    >
    > I'd bother with making a free(?) account if I were sure
    > codeforge.com is the best piece of cheese in the maze,
    > but it was just a random Google hit and
    > it seemed wiser to ask the experts here.
    > (I downloaded Cygwin just a few months ago --
    > why doesn't its gcc have sqrtl?)


    gcc doesn't have sqrtl() on any system. It's part of the runtime
    library, not the compiler. That's not just a quibble; most or all
    Linux-based systems use the GNU C library, but Cygwin apparently uses
    something else.

    I don't know *why* Cygwin's runtime library doesn't provide sqrtl(), but
    here's part of the Cygwin man page for sqrt():

    PORTABILITY
    `sqrt' is ANSI C. `sqrtf' is an extension.

    Both sqrtf and sqrtl were added to the ISO C standard in 1999, so the
    man page is badly out of date (though apparently no more so than the
    runtime library itself).

    I know there are problems with the representation of long double; gcc
    uses one size, and Microsoft uses another. This causes problems with
    the MinGW port of gcc. The Cygwin problem could be related to that.

    Of course you can always use sqrt() instead of sqrtl(), at the cost of
    some loss of precision (but if you didn't care about that you probably
    wouldn't be using long double in the first place).

    --
    Keith Thompson (The_Other_Keith) <http://www.ghoti.net/~kst>
    Will write code for food.
    "We must do something. This is something. Therefore, we must do this."
    -- Antony Jay and Jonathan Lynn, "Yes Minister"
    Keith Thompson, Nov 7, 2012
    #3
  4. On Wed, 07 Nov 2012 12:15:58 -0800, Keith Thompson wrote:

    > James Dow Allen <> writes:
    >> I've been using Cygwin lately, and am largely satisfied.
    >> It has gcc, but lacks sqrtl(), atan2l(), etc.
    >>
    >> What's the recommended remedy?
    >> I Googled "source sqrtl" and got to codeforge.com
    >> but only AFTER passing one painful Gotcha test
    >> was I told I'd need another password.
    >>
    >> I'd bother with making a free(?) account if I were sure
    >> codeforge.com is the best piece of cheese in the maze,
    >> but it was just a random Google hit and
    >> it seemed wiser to ask the experts here.
    >> (I downloaded Cygwin just a few months ago --
    >> why doesn't its gcc have sqrtl?)

    >
    >
    > Of course you can always use sqrt() instead of sqrtl(), at the cost of
    > some loss of precision (but if you didn't care about that you probably
    > wouldn't be using long double in the first place).


    As long as one is disinteresting in x = inf, NaN, and possibly -0. You
    can do

    #include <math.h>
    long double
    foo_sqrtl(long double x)
    {
    long double r;
    r = sqrt((double)x);
    return ((r + (x / r))/2);
    }

    The tricky part is making the above compliant with IEEE 754.

    --
    steve
    Steven G. Kargl, Nov 7, 2012
    #4
  5. Steven G. Kargl <> wrote:

    (snip)

    > As long as one is disinteresting in x = inf, NaN, and possibly -0. You
    > can do


    > #include <math.h>
    > long double
    > foo_sqrtl(long double x)
    > {
    > long double r;
    > r = sqrt((double)x);
    > return ((r + (x / r))/2);
    > }


    I might have done two rounds to be sure.

    Also, the last one should be

    return(x/r+(r-x/r)/2);

    if you might be on a non-base2 system, such as
    IBM hex floating point or decimal floating point.

    > The tricky part is making the above compliant with IEEE 754.


    -- glen
    glen herrmannsfeldt, Nov 7, 2012
    #5
  6. James Dow Allen

    ImpalerCore Guest

    On Nov 7, 1:40 pm, James Dow Allen <> wrote:
    > I've been using Cygwin lately, and am largely satisfied.
    > It has gcc, but lacks sqrtl(), atan2l(), etc.
    >
    > What's the recommended remedy?
    > I Googled "source sqrtl" and got to codeforge.com
    > but only AFTER passing one painful Gotcha test
    > was I told I'd need another password.
    >
    > I'd bother with making a free(?) account if I were sure
    > codeforge.com is the best piece of cheese in the maze,
    > but it was just a random Google hit and
    > it seemed wiser to ask the experts here.
    > (I downloaded Cygwin just a few months ago --
    > why doesn't its gcc have sqrtl?)


    \code
    #include <stdio.h>
    #define c_sizeof_array(arr) (sizeof (arr) / sizeof ((arr)[0]))

    const double tolerance = 0.0000001;
    const int max_iterations = 10;

    double heron_sqrt( double x );

    int main( void )
    {
    size_t idx;
    char tmp_str[40];

    double input[] = {
    0, 1, 10, 100, 1000, 10000, 0.1, 0.25, -1
    };

    for ( idx = 0; idx < c_sizeof_array (input); ++idx )
    {
    snprintf( tmp_str, sizeof (tmp_str), "sqrt( %.2f )", input[idx] );
    printf( "%-20s = %.10f\n", tmp_str, heron_sqrt( input[idx] ) );
    }

    return 0;
    }

    double heron_sqrt( double x )
    {
    /*
    * The algorithm used to calculate the square root 'x' of number 'S'
    * is based on Heron's method.
    *
    * Let x(0) be an arbitrary positive number.
    * Let x(n+1) be the average of x(n) and 'S / x(n)'
    * Repeat until the desired accuracy is achieved.
    *
    * The closer x(0) is to the square root value, the fewer number
    * of iterations are needed to converge to the accuracy required.
    *
    * The global variable 'max_iterations' allows the person to
    * limit the accuracy of the result.
    */

    double sqrt_x = 0;
    double guess = 0;
    double error = 0;
    int itr = 0;

    if ( x == 0.0 ) {
    sqrt_x = 0;
    }
    else
    {
    sqrt_x = x * 0.5;

    do
    {
    guess = sqrt_x;
    sqrt_x = 0.5 * ( guess + x / guess );
    error = sqrt_x * sqrt_x - x;
    error = error < 0 ? -error : error;
    ++itr;
    } while ( error > tolerance && itr < max_iterations );
    }

    return sqrt_x;
    }
    \endcode

    There are some issues, like how you want to deal with negative
    numbers, but you should be able to patch the algorithm up to create a
    workable substitute for 'sqrtl' by replacing the type with 'long
    double'.

    Best regards,
    John D.
    ImpalerCore, Nov 7, 2012
    #6
  7. On Wed, 07 Nov 2012 20:52:33 +0000, glen herrmannsfeldt wrote:

    > Steven G. Kargl <> wrote:
    >
    > (snip)
    >
    >> As long as one is disinteresting in x = inf, NaN, and possibly -0. You
    >> can do

    >
    >> #include <math.h>
    >> long double
    >> foo_sqrtl(long double x)
    >> {
    >> long double r;
    >> r = sqrt((double)x);
    >> return ((r + (x / r))/2);
    >> }

    >
    > I might have done two rounds to be sure.


    To be sure of what? OP is using cygwin. sqrt()
    will be a 53-bit double precision result, and OP
    wants long double sqrtl() which is a 64-bit
    floating point type. One round would give
    106-bit.

    > Also, the last one should be
    >
    > return(x/r+(r-x/r)/2);


    I beg to differ.

    #include <stdio.h>
    #include <math.h>
    long double foo_sqrtl(long double x)
    {
    long double r;
    r = sqrt((double)x);
    return ((r + (x / r))/2);
    }

    int main(void)
    {
    long double x;
    x = M_PI / M_LN2;
    printf("%.18Le\n%.18Le\n", sqrtl(x), foo_sqrtl(x));
    return (0);
    }

    % cc -o foo foo.c -lm && ./foo
    2.128934038862452440e+00
    2.128934038862452440e+00

    > if you might be on a non-base2 system, such as
    > IBM hex floating point or decimal floating point.


    OP is using cygwin.

    --
    steve
    Steven G. Kargl, Nov 7, 2012
    #7
  8. James Dow Allen, Nov 7, 2012
    #8
  9. Steven G. Kargl <> wrote:
    >> Steven G. Kargl <> wrote:


    (snip)
    >>> #include <math.h>
    >>> long double
    >>> foo_sqrtl(long double x)
    >>> {
    >>> long double r;
    >>> r = sqrt((double)x);
    >>> return ((r + (x / r))/2);
    >>> }


    >> I might have done two rounds to be sure.


    > To be sure of what? OP is using cygwin. sqrt()
    > will be a 53-bit double precision result, and OP
    > wants long double sqrtl() which is a 64-bit
    > floating point type. One round would give
    > 106-bit.


    I thought that there were some systems with software 128 bit
    floating point for long double. Yes that should be fine to
    get from 53 to 64 bits.

    >> Also, the last one should be


    >> return(x/r+(r-x/r)/2);


    > I beg to differ.


    (snip)

    >> if you might be on a non-base2 system, such as
    >> IBM hex floating point or decimal floating point.


    > OP is using cygwin.


    Well, much of GNU libc in general was mentioned in the thread.

    There seem to be very few hardware implementations of 128 bit
    floating point, other than IBM (back to the 360/85 and continuing
    to present day processors) and VAX H-float (microcode on some systems).

    Single precision sqrt for S/360:

    DE FR0,BUFF GIVE TWO PASSES OF NEWTON-RAPHSON
    AU FR0,BUFF ITERATION
    HER FR0,FR0
    DER FR2,FR0 (X/Y1+Y1)/2 = (Y1-X/Y1)/2+X/Y1 TO GUARD
    AU FR0,ROUND LAST DIGIT-. ADD ROUNDING FUDGE
    SER FR0,FR2
    HER FR0,FR0
    AER FR0,FR2

    -- glen
    glen herrmannsfeldt, Nov 7, 2012
    #9
  10. Keith Thompson wrote:

    [snip]

    > gcc doesn't have sqrtl() on any system. It's part of the runtime
    > library, not the compiler. That's not just a quibble; most or all
    > Linux-based systems use the GNU C library, but Cygwin apparently uses
    > something else.

    [...]

    Just for my own information and for future reference. Are functions like
    this sqrtl and other GNU C function OT here? Isn't this group restricted to
    ANSI and ISO C ?

    Bill
    Bill Cunningham, Nov 7, 2012
    #10
  11. James Dow Allen

    Fred K Guest

    On Wednesday, November 7, 2012 2:29:23 PM UTC-8, Bill Cunningham wrote:
    > Keith Thompson wrote: [snip] > gcc doesn't have sqrtl() on any system. It's part of the runtime > library, not the compiler. That's not just a quibble; most or all > Linux-based systems use the GNU C library, but Cygwin apparently uses > something else. [...] Just for my own information and for future reference. Are functions like this sqrtl and other GNU C function OT here? Isn't this group restricted to ANSI and ISO C ? Bill


    Read Keith Thompson's first reply:
    "Both sqrtf and sqrtl were added to the ISO C standard in 1999"
    Fred K, Nov 7, 2012
    #11
  12. On Nov 7, 4:29 pm, "Bill Cunningham" <> wrote:
    > Keith Thompson wrote:
    >
    > [snip]
    >
    > > gcc doesn't have sqrtl() on any system.  It's part of the runtime
    > > library, not the compiler.  That's not just a quibble; most or all
    > > Linux-based systems use the GNU C library, but Cygwin apparently uses
    > > something else.

    >
    > [...]
    >
    > Just for my own information and for future reference. Are functions like
    > this sqrtl and other GNU C function OT here? Isn't this group restricted to
    > ANSI and ISO C ?
    >
    > Bill


    Couple of lines below the point where you snipped his post, Keith
    states "Both sqrtf and sqrtl were added to the ISO C standard in
    1999".

    - Anand
    Anand Hariharan, Nov 7, 2012
    #12
  13. Fred K wrote:
    > On Wednesday, November 7, 2012 2:29:23 PM UTC-8, Bill Cunningham
    > wrote:
    >> Keith Thompson wrote: [snip] > gcc doesn't have sqrtl() on any
    >> system. It's part of the runtime > library, not the compiler. That's
    >> not just a quibble; most or all > Linux-based systems use the GNU C
    >> library, but Cygwin apparently uses > something else. [...] Just for
    >> my own information and for future reference. Are functions like this
    >> sqrtl and other GNU C function OT here? Isn't this group restricted
    >> to ANSI and ISO C ? Bill

    >
    > Read Keith Thompson's first reply:
    > "Both sqrtf and sqrtl were added to the ISO C standard in 1999"


    Oh I see I'm sorry. Missed that. So it's on topic since it's ISO C. But
    GNU extensions are not?

    Bill
    Bill Cunningham, Nov 7, 2012
    #13
  14. "Steven G. Kargl" <> writes:
    > On Wed, 07 Nov 2012 12:15:58 -0800, Keith Thompson wrote:

    [...]
    >> Of course you can always use sqrt() instead of sqrtl(), at the cost of
    >> some loss of precision (but if you didn't care about that you probably
    >> wouldn't be using long double in the first place).

    >
    > As long as one is disinteresting in x = inf, NaN, and possibly -0. You
    > can do
    >
    > #include <math.h>
    > long double
    > foo_sqrtl(long double x)
    > {
    > long double r;
    > r = sqrt((double)x);
    > return ((r + (x / r))/2);
    > }
    >
    > The tricky part is making the above compliant with IEEE 754.


    Cool.

    Experiment shows that it works correctly for 1.0 and 2.0. Therefore, by
    mathematical induction, it works perfectly for all argument values.

    :cool:}

    (Incidentally, the cast is unnecessary but harmless.)

    --
    Keith Thompson (The_Other_Keith) <http://www.ghoti.net/~kst>
    Will write code for food.
    "We must do something. This is something. Therefore, we must do this."
    -- Antony Jay and Jonathan Lynn, "Yes Minister"
    Keith Thompson, Nov 7, 2012
    #14
  15. On Nov 8, 6:40 am, pete <> wrote:

    > I have some portable freestanding implementations of
    >  double fs_log2(double x);
    >  [et cetera]
    >    http://www.mindspring.com/~pfilandr/C/fs_math/fs_math.h
    >    http://www.mindspring.com/~pfilandr/C/fs_math/fs_math.c


    Thank you! I've downloaded these and will start using them as
    needed. I particularly like that they're "freestanding"
    and in just two files. (Often I've given up trying to use
    on-line code, rather than face the annoyance of tracking down
    prerequisites.)

    James
    James Dow Allen, Nov 8, 2012
    #15
  16. Keith Thompson <> wrote:

    (snip)
    >> #include <math.h>
    >> long double
    >> foo_sqrtl(long double x)
    >> {
    >> long double r;
    >> r = sqrt((double)x);
    >> return ((r + (x / r))/2);
    >> }


    (snip)
    > Experiment shows that it works correctly for 1.0 and 2.0. Therefore, by
    > mathematical induction, it works perfectly for all argument values.


    > :cool:}


    > (Incidentally, the cast is unnecessary but harmless.)


    It reminds me of the narrowing casts required by Java.
    (Except I think Java doesn't have long double yet, but when it does...)

    -- glen
    glen herrmannsfeldt, Nov 8, 2012
    #16
  17. James Dow Allen

    ImpalerCore Guest

    On Nov 7, 6:40 pm, pete <> wrote:
    > ImpalerCore wrote:
    >
    > > On Nov 7, 1:40 pm, James Dow Allen <> wrote:
    > > > I've been using Cygwin lately, and am largely satisfied.
    > > > It has gcc, but lacks sqrtl(), atan2l(), etc.

    >
    > > > What's the recommended remedy?
    > > > I Googled "source sqrtl" and got to codeforge.com
    > > > but only AFTER passing one painful Gotcha test
    > > > was I told I'd need another password.

    >
    > > > I'd bother with making a free(?) account if I were sure
    > > > codeforge.com is the best piece of cheese in the maze,
    > > > but it was just a random Google hit and
    > > > it seemed wiser to ask the experts here.
    > > > (I downloaded Cygwin just a few months ago --
    > > > why doesn't its gcc have sqrtl?)

    >
    > > \code
    > > #include <stdio.h>
    > > #define c_sizeof_array(arr) (sizeof (arr) / sizeof ((arr)[0]))

    >
    > > const double tolerance = 0.0000001;
    > > const int max_iterations = 10;

    >
    > > double heron_sqrt( double x );

    >
    > > int main( void )
    > > {
    > >   size_t idx;
    > >   char tmp_str[40];

    >
    > >   double input[] = {
    > >     0, 1, 10, 100, 1000, 10000, 0.1, 0.25, -1
    > >   };

    >
    > >   for ( idx = 0; idx < c_sizeof_array (input); ++idx )
    > >   {
    > >     snprintf( tmp_str, sizeof (tmp_str), "sqrt( %.2f )", input[idx]);
    > >     printf( "%-20s = %.10f\n", tmp_str, heron_sqrt( input[idx] ) );
    > >   }

    >
    > >   return 0;
    > > }

    >
    > > double heron_sqrt( double x )
    > > {
    > >   /*
    > >    * The algorithm used to calculate the square root 'x' of number 'S'
    > >    * is based on Heron's method.
    > >    *
    > >    * Let x(0) be an arbitrary positive number.
    > >    * Let x(n+1) be the average of x(n) and 'S / x(n)'
    > >    * Repeat until the desired accuracy is achieved.
    > >    *
    > >    * The closer x(0) is to the square root value, the fewer number
    > >    * of iterations are needed to converge to the accuracy required.
    > >    *
    > >    * The global variable 'max_iterations' allows the person to
    > >    * limit the accuracy of the result.
    > >    */

    >
    > >   double sqrt_x = 0;
    > >   double guess = 0;
    > >   double error = 0;
    > >   int itr = 0;

    >
    > >   if ( x == 0.0 ) {
    > >     sqrt_x = 0;
    > >   }
    > >   else
    > >   {
    > >     sqrt_x = x * 0.5;

    >
    > >     do
    > >     {
    > >       guess = sqrt_x;
    > >       sqrt_x = 0.5 * ( guess + x / guess );
    > >       error = sqrt_x * sqrt_x - x;
    > >       error = error < 0 ? -error : error;
    > >       ++itr;
    > >     } while ( error > tolerance && itr < max_iterations );
    > >   }

    >
    > >   return sqrt_x;
    > > }
    > > \endcode

    >
    > > There are some issues, like how you want to deal with negative
    > > numbers, but you should be able to patch the algorithm up to create a
    > > workable substitute for 'sqrtl' by replacing the type with 'long
    > > double'.

    >
    > I prefer this implementation of Hero's method
    > for finding square roots:
    >
    > double
    > square_root(double x)
    > {
    >     double root = x;
    >
    >     if (root > 0) {
    >         double next = root / 2 + 0.5;
    >
    >         do {
    >             root = next;
    >             next = (x / root + root) / 2;
    >         } while (root > next);
    >     }
    >     return root;
    >
    > }
    >
    > (square_root(x)) returns the value of the largest double
    > which is less than or equal to the square root of (x).


    Thanks for the alternate implementation. The version I gave was from
    an example from when I was mucking around contract programming, but
    with the contracts stripped out.

    \code
    #include <stdio.h>
    #include <clover/macros.h>
    #include <clover/contract.h>

    const double tolerance = 0.0000001;
    const int max_iterations = 10;

    double heron_sqrt( double x );

    int main( void )
    {
    size_t idx;
    char tmp_str[40];

    double input[] = {
    0, 1, 10, 100, 1000, 10000, 0.1, 0.25, -1
    };

    for ( idx = 0; idx < c_sizeof_array (input); ++idx )
    {
    snprintf( tmp_str, sizeof (tmp_str), "sqrt( %.2f )", input[idx] );
    printf( "%-20s = %.10f\n", tmp_str, heron_sqrt( input[idx] ) );
    }

    return 0;
    }

    double heron_sqrt( double x )
    {
    /*
    * The algorithm used to calculate the square root 'x' of number 'S'
    * is based on Heron's method.
    *
    * Let x(0) be an arbitrary positive number.
    * Let x(n+1) be the average of x(n) and 'S / x(n)'
    * Repeat until the desired accuracy is achieved.
    *
    * The closer x(0) is to the square root value, the fewer number
    * of iterations are needed to converge to the accuracy required.
    *
    * The global variable 'max_iterations' allows the person to
    * limit the accuracy of the result to cause postcondition
    * violations.
    */

    double sqrt_x = 0;
    double guess = 0;
    double error = 0;
    int itr = 0;
    C_POST( double x_sq; )

    C_REQUIRE( x >= 0 );

    if ( x == 0.0 ) {
    sqrt_x = 0;
    }
    else
    {
    sqrt_x = x * 0.5;

    do
    {
    guess = sqrt_x;
    sqrt_x = 0.5 * ( guess + x / guess );
    error = sqrt_x * sqrt_x - x;
    error = error < 0 ? -error : error;
    ++itr;
    } while ( error > tolerance && itr < max_iterations );
    }

    C_POST( x_sq = sqrt_x * sqrt_x; )
    C_ENSURE( x_sq - tolerance <= x && x <= x_sq + tolerance );

    return sqrt_x;
    }
    \endcode

    C_REQUIRE and C_ENSURE are basically glorified asserts that separate
    pre and post conditions (which are either elided or evaluated
    depending on the level of contract evaluation), and C_POST is a macro
    to hide or declare the variables and statements when evaluating
    postconditions. I used 'max_iterations' as a driver to test
    postcondition evaluation.

    > I have some portable freestanding implementations of
    >  double fs_log2(double x);
    >  double fs_exp2(double x);
    >  long double fs_powl(long double x, long double y);
    >  long double fs_sqrtl(long double x);
    >  long double fs_logl(long double x);
    >  long double fs_expl(long double x);
    >  long double fs_cosl(long double x);
    >  long double fs_fmodl(long double x, long double y);
    > declared here:
    >    http://www.mindspring.com/~pfilandr/C/fs_math/fs_math.h
    > and defined here:
    >    http://www.mindspring.com/~pfilandr/C/fs_math/fs_math.c


    I've visited your site several times in the past and appreciate your
    generosity of sharing your code. I hope to do the same eventually
    once I get my code to a respectable state.

    > In terms of accuracy, my square root function is as good as any.
    > The pow, log, and exp functions are not quite as good.
    >
    > I like the performance of my cosine function.
    > Given:
    >     #include <math.h>
    >     #define DEGREES_TO_RADIANS(A)   ((A) * atan(1) / 45)
    > on the implementation that I am using,
    > using my fs_cos,
    > the value of (fs_cos(DEGREES_TO_RADIANS(-3540)) - 0.5))
    > is 9.992007e-016.
    > The value obtatined by using the supplied standard library
    > cos function in (cos(DEGREES_TO_RADIANS(-3540)) - 0.5)
    > is -1.060133e-015.
    > Mathematically, the value should be zero.


    One of these days I need to take a deeper look at your math and sort
    function implementations.

    Best regards,
    John D.
    ImpalerCore, Nov 8, 2012
    #17
  18. James Dow Allen

    ImpalerCore Guest

    On Nov 7, 6:40 pm, Keith Thompson <> wrote:
    > "Steven G. Kargl" <> writes:
    > > On Wed, 07 Nov 2012 12:15:58 -0800, Keith Thompson wrote:

    > [...]
    > >> Of course you can always use sqrt() instead of sqrtl(), at the cost of
    > >> some loss of precision (but if you didn't care about that you probably
    > >> wouldn't be using long double in the first place).

    >
    > > As long as one is disinteresting in x = inf, NaN, and possibly -0.  You
    > > can do

    >
    > > #include <math.h>
    > > long double
    > > foo_sqrtl(long double x)
    > > {
    > >    long double r;
    > >    r = sqrt((double)x);
    > >    return ((r + (x / r))/2);
    > > }

    >
    > > The tricky part is making the above compliant with IEEE 754.

    >
    > Cool.
    >
    > Experiment shows that it works correctly for 1.0 and 2.0.  Therefore, by
    > mathematical induction, it works perfectly for all argument values.
    >
    > :cool:}
    >
    > (Incidentally, the cast is unnecessary but harmless.)


    I believe that the cast is there to make the loss of precision more
    explicit, e.g. for documentation rather than semantic purposes, and to
    avoid potential slaps on the wrist from the compiler. I tend to favor
    that casting style as well for any type conversion that loses
    precision for statements that do not include a declaration of the
    variable type.

    \code snippet
    double x;
    ....
    int y = x; /* ok, type explicit in statement */
    \endcode

    \code snippet
    double a;
    int b;

    ....

    b = (int)a;
    \endcode

    Best regards,
    John D.
    ImpalerCore, Nov 8, 2012
    #18
  19. James Dow Allen

    Noob Guest

    Noob, Nov 8, 2012
    #19
  20. James Dow Allen

    Philip Lantz Guest

    Steven G. Kargl wrote:
    > glen herrmannsfeldt wrote:
    > > Steven G. Kargl wrote:
    > >> #include <math.h>
    > >> long double
    > >> foo_sqrtl(long double x)
    > >> {
    > >> long double r;
    > >> r = sqrt((double)x);
    > >> return ((r + (x / r))/2);
    > >> }

    >
    > > ... the last one should be
    > >
    > > return(x/r+(r-x/r)/2);

    >
    > I beg to differ.


    Mathematically, ((r + x/r)/2) and (x/r + (r - x/r)/2) are equivalent. Is
    there a numerical reason to prefer one over the other? (Obviously, if
    there are no such concerns, one might as well use the simpler one.)
    Philip Lantz, Nov 8, 2012
    #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. Replies:
    0
    Views:
    404
  2. Replies:
    0
    Views:
    410
  3. WELCOME to

    , Sep 3, 2003, in forum: Python
    Replies:
    1
    Views:
    496
    Robin Becker
    Sep 5, 2003
  4. Replies:
    0
    Views:
    395
  5. Kevin Walzer

    Re: PIL (etc etc etc) on OS X

    Kevin Walzer, Aug 1, 2008, in forum: Python
    Replies:
    4
    Views:
    356
    Fredrik Lundh
    Aug 13, 2008
Loading...

Share This Page