sqrt() double trouble

**Walter Roberson** · Jun 27 '08, 07:40 PM

Re: sqrt() double trouble

In article <g2ajo7$efb$1@i news.gazeta.pl> , copx <copx@gazeta.pl wrote:

>For those who will certainly complain that discussing the use of CPU
>specific instructions is OT here - I say this question is on topic here
>because my problem is caused by the C standard which only features double
>precision sqrt() and thus does not allow me to exploit the capabilities of
>the target CPUs.

If you have a problem with the C standard and want it changed, then
take the issue up in comp.std.c .

Personally, though, I wouldn't bother doing that: I'd just use
the standard-mandated float square root function,

#include <math.h>

float sqrtf(float x)

You will find this in section 7.12.7.5 of the C99 standard.
--
"The Romans believed that every man had his Genius, and every
woman her Juno." -- Thomas Bulfinch

**copx** · Jun 27 '08, 07:40 PM

Re: sqrt() double trouble

"Walter Roberson" <roberson@ibd.n rc-cnrc.gc.caschri eb im Newsbeitrag
news:g2akk8$ro0 $1@canopus.cc.u manitoba.ca...

In article <g2ajo7$efb$1@i news.gazeta.pl> , copx <copx@gazeta.pl wrote:
>

>>For those who will certainly complain that discussing the use of CPU
>>specific instructions is OT here - I say this question is on topic here
>>because my problem is caused by the C standard which only features double
>>precision sqrt() and thus does not allow me to exploit the capabilities of
>>the target CPUs.

>
If you have a problem with the C standard and want it changed, then
take the issue up in comp.std.c .
>
Personally, though, I wouldn't bother doing that: I'd just use
the standard-mandated float square root function,
>
#include <math.h>
>
float sqrtf(float x)
>
You will find this in section 7.12.7.5 of the C99 standard.

Is this a new function of C99 or is it also part of C89?

**copx** · Jun 27 '08, 07:40 PM

Re: sqrt() double trouble

"copx" <copx@gazeta.pl schrieb im Newsbeitrag
news:g2akru$kp4 $1@inews.gazeta .pl...

>
"Walter Roberson" <roberson@ibd.n rc-cnrc.gc.caschri eb im Newsbeitrag
news:g2akk8$ro0 $1@canopus.cc.u manitoba.ca...

>In article <g2ajo7$efb$1@i news.gazeta.pl> , copx <copx@gazeta.pl wrote:
>>

>>>For those who will certainly complain that discussing the use of CPU
>>>specific instructions is OT here - I say this question is on topic here
>>>because my problem is caused by the C standard which only features double
>>>precision sqrt() and thus does not allow me to exploit the capabilities
>>>of
>>>the target CPUs.

>>
>If you have a problem with the C standard and want it changed, then
>take the issue up in comp.std.c .
>>
>Personally, though, I wouldn't bother doing that: I'd just use
>the standard-mandated float square root function,
>>
>#include <math.h>
>>
> float sqrtf(float x)
>>
>You will find this in section 7.12.7.5 of the C99 standard.

>
Is this a new function of C99 or is it also part of C89?

I have looked this up in 3 different C standard library references and the
result suggests that this is a C99-only thing. Well, at least it seems
someone on the standard comitee was aware of the issue.

**Tim Prince** · Jun 27 '08, 07:40 PM

Re: sqrt() double trouble

copx wrote:

"Walter Roberson" <roberson@ibd.n rc-cnrc.gc.caschri eb im Newsbeitrag
news:g2akk8$ro0 $1@canopus.cc.u manitoba.ca...

>In article <g2ajo7$efb$1@i news.gazeta.pl> , copx <copx@gazeta.pl wrote:
>>

>>For those who will certainly complain that discussing the use of CPU
>>specific instructions is OT here - I say this question is on topic here
>>because my problem is caused by the C standard which only features double
>>precision sqrt() and thus does not allow me to exploit the capabilities of
>>the target CPUs.

>If you have a problem with the C standard and want it changed, then
>take the issue up in comp.std.c .
>>
>Personally, though, I wouldn't bother doing that: I'd just use
>the standard-mandated float square root function,
>>
>#include <math.h>
>>
> float sqrtf(float x)
>>
>You will find this in section 7.12.7.5 of the C99 standard.

>
Is this a new function of C99 or is it also part of C89?
>
>
>

C89 reserved the name for this usage, but made its presence optional. For
some versions of MSVC you had to set an option to make it available. In
C99 there's also <tgmath.h. Apparently, you are among those who don't
consider C99 to be C (do you take K&R 1 as your standard?), but you are a
bit late to be asking for a different way of doing things.

**copx** · Jun 27 '08, 07:40 PM

Re: sqrt() double trouble

"Tim Prince" <tprince@nospam computer.orgsch rieb im Newsbeitrag
news:cT42k.991$ L_.890@flpi150. ffdc.sbc.com...
[snip]

C89 reserved the name for this usage, but made its presence optional. For
some versions of MSVC you had to set an option to make it available. In
C99 there's also <tgmath.h. Apparently, you are among those who don't
consider C99 to be C (do you take K&R 1 as your standard?), but you are a
bit late to be asking for a different way of doing things.

It's a question of portability. C99 is a trainwreck of a standard - at least
on the desktop - because those who matter don't support it. MS owns 90% of
the desktop market, their compiler is the standard for almost all desktop
development - and they do not support C99 nor seem to have the intention to
ever add support for it. Borland (now rather irrelevant I know) does not
support C99 either I think. You are stuck with a hobbiest port of GCC
(MinGW) if you want to develop non-trivial C99 apps on Windows. All the
other alternatives are even worse.
I wish C99 would be supported better so that it could actually replace C89,
but so far fully conforming C99 compilers are exotic beasts while C89
compilers are available almost everywhere. That's why C99 is not a
replacement for C89 - unfortunately.
It seems that the language is considered a dead/legacy support thing by most
of the industry. I mean, it is not like C99 is such a complex standard to
implement, right? MS could probably update their compiler to C99 in a month
or Apple/Redhat/IBM etc. could pay some people to fix GCC's last C99
conformance issues in a similar time frame, but there seems to be not enough
interest.

**Keith Thompson** · Jun 27 '08, 07:40 PM

Re: sqrt() double trouble

"copx" <copx@gazeta.pl writes:

"Walter Roberson" <roberson@ibd.n rc-cnrc.gc.caschri eb im Newsbeitrag
news:g2akk8$ro0 $1@canopus.cc.u manitoba.ca...

>In article <g2ajo7$efb$1@i news.gazeta.pl> , copx <copx@gazeta.pl wrote:
>>

>>>For those who will certainly complain that discussing the use of CPU
>>>specific instructions is OT here - I say this question is on topic here
>>>because my problem is caused by the C standard which only features double
>>>precision sqrt() and thus does not allow me to exploit the capabilities of
>>>the target CPUs.

>>
>If you have a problem with the C standard and want it changed, then
>take the issue up in comp.std.c .
>>
>Personally, though, I wouldn't bother doing that: I'd just use
>the standard-mandated float square root function,
>>
>#include <math.h>
>>
> float sqrtf(float x)
>>
>You will find this in section 7.12.7.5 of the C99 standard.

>
Is this a new function of C99 or is it also part of C89?

C90 only has sqrt(double). C99 adds sqrtf(float) and sqrtl(long
double) (as well as the type-generic sqrt() macro in <tgmath.h>).

If your implementation doesn't provide sqrtf(), you can probably find
an open source implementation of it.

--
Keith Thompson (The_Other_Keit h) kst-u@mib.org <http://www.ghoti.net/~kst>
Nokia
"We must do something. This is something. Therefore, we must do this."
-- Antony Jay and Jonathan Lynn, "Yes Minister"

**jacob navia** · Jun 27 '08, 07:40 PM

Re: sqrt() double trouble

copx wrote:

"copx" <copx@gazeta.pl schrieb im Newsbeitrag
news:g2akru$kp4 $1@inews.gazeta .pl...

>"Walter Roberson" <roberson@ibd.n rc-cnrc.gc.caschri eb im Newsbeitrag
>news:g2akk8$ro 0$1@canopus.cc. umanitoba.ca...

>>In article <g2ajo7$efb$1@i news.gazeta.pl> , copx <copx@gazeta.pl wrote:
>>>
>>>For those who will certainly complain that discussing the use of CPU
>>>specific instructions is OT here - I say this question is on topic here
>>>because my problem is caused by the C standard which only features double
>>>precision sqrt() and thus does not allow me to exploit the capabilities
>>>of
>>>the target CPUs.
>>If you have a problem with the C standard and want it changed, then
>>take the issue up in comp.std.c .
>>>
>>Personally, though, I wouldn't bother doing that: I'd just use
>>the standard-mandated float square root function,
>>>
>>#include <math.h>
>>>
>> float sqrtf(float x)
>>>
>>You will find this in section 7.12.7.5 of the C99 standard.

>Is this a new function of C99 or is it also part of C89?

>
I have looked this up in 3 different C standard library references and the
result suggests that this is a C99-only thing. Well, at least it seems
someone on the standard comitee was aware of the issue.
>
>
>

/* sqrtf.c
*
* Square root
*
*
*
* SYNOPSIS:
*
* float x, y, sqrtf();
*
* y = sqrtf( x );
*
*
*
* DESCRIPTION:
*
* Returns the square root of x.
*
* Range reduction involves isolating the power of two of the
* argument and using a polynomial approximation to obtain
* a rough value for the square root. Then Heron's iteration
* is used three times to converge to an accurate value.
*
*
*
* ACCURACY:
*
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,1.e38 100000 8.7e-8 2.9e-8
*
*
* ERROR MESSAGES:
*
* message condition value returned
* sqrtf domain x < 0 0.0
*
*/

/*
Cephes Math Library Release 2.2: June, 1992
Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

/* Single precision square root
* test interval: [sqrt(2)/2, sqrt(2)]
* trials: 30000
* peak relative error: 8.8e-8
* rms relative error: 3.3e-8
*
* test interval: [0.01, 100.0]
* trials: 50000
* peak relative error: 8.7e-8
* rms relative error: 3.3e-8
*
* Copyright (C) 1989 by Stephen L. Moshier. All rights reserved.
*/
#include "mconf.h"

float frexpf( float, int * );
float ldexpf( float, int );

float sqrtf( float xx )#else
float frexpf(), ldexpf();

float sqrtf(float xx)
{
float f, x, y;
int e;

f = xx;
if( f <= 0.0 ) {
if( f < 0.0 )
//mtherr( "sqrtf", DOMAIN );
return( 0.0 );
}

x = frexpf( f, &e ); /* f = x * 2**e, 0.5 <= x < 1.0 */
/* If power of 2 is odd, double x and decrement the power of 2. */
if( e & 1 ) {
x = x + x;
e -= 1;
}

e >>= 1; /* The power of 2 of the square root. */

if( x 1.41421356237 ) {
/* x is between sqrt(2) and 2. */
x = x - 2.0;
y =
((((( -9.8843065718E-4 * x
+ 7.9479950957E-4) * x
- 3.5890535377E-3) * x
+ 1.1028809744E-2) * x
- 4.4195203560E-2) * x
+ 3.5355338194E-1) * x
+ 1.41421356237E0 ;
goto sqdon;
}

if( x 0.707106781187 ) {
/* x is between sqrt(2)/2 and sqrt(2). */
x = x - 1.0;
y =
((((( 1.35199291026E-2 * x
- 2.26657767832E-2) * x
+ 2.78720776889E-2) * x
- 3.89582788321E-2) * x
+ 6.24811144548E-2) * x
- 1.25001503933E-1) * x * x
+ 0.5 * x
+ 1.0;
goto sqdon;
}

/* x is between 0.5 and sqrt(2)/2. */
x = x - 0.5;
y =
((((( -3.9495006054E-1 * x
+ 5.1743034569E-1) * x
- 4.3214437330E-1) * x
+ 3.5310730460E-1) * x
- 3.5354581892E-1) * x
+ 7.0710676017E-1) * x
+ 7.07106781187E-1;

sqdon:
y = ldexpf( y, e ); /* y = y * 2**e */
return( y);
}

--
jacob navia
jacob at jacob point remcomp point fr
logiciels/informatique

403 Forbidden

http://www.cs.virginia.edu/~lcc-win32

**Keith Thompson** · Jun 27 '08, 07:40 PM

Re: sqrt() double trouble

Tim Prince <tprince@nospam computer.orgwri tes:

copx wrote:

>"Walter Roberson" <roberson@ibd.n rc-cnrc.gc.caschri eb im
>Newsbeitrag news:g2akk8$ro0 $1@canopus.cc.u manitoba.ca...

>>In article <g2ajo7$efb$1@i news.gazeta.pl> , copx <copx@gazeta.pl wrote:
>>>
>>>For those who will certainly complain that discussing the use of CPU
>>>specific instructions is OT here - I say this question is on topic here
>>>because my problem is caused by the C standard which only features double
>>>precision sqrt() and thus does not allow me to exploit the capabilities of
>>>the target CPUs.
>>If you have a problem with the C standard and want it changed, then
>>take the issue up in comp.std.c .
>>>
>>Personally, though, I wouldn't bother doing that: I'd just use
>>the standard-mandated float square root function,
>>>
>>#include <math.h>
>>>
>> float sqrtf(float x)
>>>
>>You will find this in section 7.12.7.5 of the C99 standard.

>Is this a new function of C99 or is it also part of C89?
>>

C89 reserved the name for this usage, but made its presence optional.

Interesting. You're right (C90 7.13.4), but I hadn't realized that.

I wonder why the committee reserved the names but didn't require the
functions to be implemented.

For some versions of MSVC you had to set an option to make it
available. In C99 there's also <tgmath.h. Apparently, you are
among those who don't consider C99 to be C (do you take K&R 1 as your
standard?), but you are a bit late to be asking for a different way of
doing things.

How does asking whether something is C99-specific imply anything more
than wanting to know whether it's C99-specific?

--
Keith Thompson (The_Other_Keit h) kst-u@mib.org <http://www.ghoti.net/~kst>
Nokia
"We must do something. This is something. Therefore, we must do this."
-- Antony Jay and Jonathan Lynn, "Yes Minister"

**Martin Ambuhl** · Jun 27 '08, 07:40 PM

Re: sqrt() double trouble

copx wrote:

I am writing a program which uses sqrt() a lot. A historically very slow
function but I read on CPUs with SSE support this is actually fast. Problem:
C's sqrt() (unlike C++ sqrt()) is defined to work on doubles only, and early
versions of SSE did not support double precision. The CPU requirement
"something with at least first generation SEE support" (everything from
P3/Athlon XP and up) is acceptable for me. However, requiring a CPU which
supports double precision SSE is too much.
Is there any work around for this problem? (except for switching to C++)
>
For those who will certainly complain that discussing the use of CPU
specific instructions is OT here - I say this question is on topic here
because my problem is caused by the C standard which only features double
precision sqrt() and thus does not allow me to exploit the capabilities of
the target CPUs.

In no way does C++'s alleged support for long double arguments to sqrt()
have anything to do with your problem. C++'s long double sqrt() cannot
magically turn your CPU into one with SSE support for long double
sqrt(). All C++ can do for you is substitute a call to some
sqrt_for_long_d ouble() for sqrt() behind your back. The same
functionality is available in C if you get a version of the math library
supporting C99's sqrtl() function, or any of the packages of math
routines for C with long double functionality, such as the Cephes one
<http://www.netlib.org/cephes/ldouble.tgz>.

**santosh** · Jun 27 '08, 07:40 PM

Re: sqrt() double trouble

copx wrote:

>
"Tim Prince" <tprince@nospam computer.orgsch rieb im Newsbeitrag
news:cT42k.991$ L_.890@flpi150. ffdc.sbc.com...
[snip]

>C89 reserved the name for this usage, but made its presence optional.
> For
>some versions of MSVC you had to set an option to make it available.
>In
>C99 there's also <tgmath.h. Apparently, you are among those who
>don't consider C99 to be C (do you take K&R 1 as your standard?), but
>you are a bit late to be asking for a different way of doing things.

>
It's a question of portability. C99 is a trainwreck of a standard - at
least on the desktop - because those who matter don't support it. MS
owns 90% of the desktop market, their compiler is the standard for
almost all desktop development - and they do not support C99 nor seem
to have the intention to ever add support for it. Borland (now rather
irrelevant I know) does not support C99 either I think. You are stuck
with a hobbiest port of GCC (MinGW) if you want to develop non-trivial
C99 apps on Windows. All the other alternatives are even worse.

What about Intel's compiler? It's a commercial product and seems to
support the vast majority of C99 features, even though their
documentation tells you otherwise.

<snip>

**Dan** · Jun 27 '08, 07:40 PM

Re: sqrt() double trouble

"copx" <copx@gazeta.pl wrote in message
news:g2ajo7$efb $1@inews.gazeta .pl...

>I am writing a program which uses sqrt() a lot. A historically very slow
>function but I read on CPUs with SSE support this is actually fast.
>Problem: C's sqrt() (unlike C++ sqrt()) is defined to work on doubles only,
>and early versions of SSE did not support double precision. The CPU
>requirement "something with at least first generation SEE support"
>(everything from P3/Athlon XP and up) is acceptable for me. However,
>requiring a CPU which supports double precision SSE is too much.
Is there any work around for this problem? (except for switching to C++)
>
For those who will certainly complain that discussing the use of CPU
specific instructions is OT here - I say this question is on topic here
because my problem is caused by the C standard which only features double
precision sqrt() and thus does not allow me to exploit the capabilities of
the target CPUs.

Does anyone have the code for the routine handy? Is it that hard to change
it to float?

**Chris Dollin** · Jun 27 '08, 07:40 PM

Re: sqrt() double trouble

copx wrote:

I am writing a program which uses sqrt() a lot. A historically very slow
function but I read on CPUs with SSE support this is actually fast. Problem:
C's sqrt() (unlike C++ sqrt()) is defined to work on doubles only, and early
versions of SSE did not support double precision. The CPU requirement
"something with at least first generation SEE support" (everything from
P3/Athlon XP and up) is acceptable for me. However, requiring a CPU which
supports double precision SSE is too much.
Is there any work around for this problem? (except for switching to C++)
>
For those who will certainly complain that discussing the use of CPU
specific instructions is OT here - I say this question is on topic here
because my problem is caused by the C standard which only features double
precision sqrt() and thus does not allow me to exploit the capabilities of
the target CPUs.

Doesn't it?

Suppose the compiler noticed the expression `(float) sqrt( aFloat )` (or
moral equivalents) and implemented it with a machine-specific sqrt-float
instruction or millicode. Would that violate the Standard?

If not, then the Standard isn't preventing you from exploiting the
capabilities of the target CPUs -- it's the implementations .

(I'm assuming you have independent measurements that show that sqrt-speed
really is an important performance bottleneck for your application.)

--
"The letter was not unproductive." /Mansfield Park/

Hewlett-Packard Limited registered no:
registered office: Cain Road, Bracknell, Berks RG12 1HN 690597 England

**Ben Bacarisse** · Jun 27 '08, 07:40 PM

Re: sqrt() double trouble

jacob navia <jacob@nospam.c omwrites:

<snip discussion>

/* sqrtf.c

<snip>

* SYNOPSIS:
*
* float x, y, sqrtf();
*
* y = sqrtf( x );

<snip>

*/
>
/*
Cephes Math Library Release 2.2: June, 1992
Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

Modified by you I think.

<snip>

float sqrtf(float xx)
{
float f, x, y;
int e;
>
f = xx;
if( f <= 0.0 ) {
if( f < 0.0 )
//mtherr( "sqrtf", DOMAIN );
return( 0.0 );
}

Your C99 comment introduces a bug. I hope this is not the code used
in your library. I leave the rest so you can see it assumes that f !=
0.0f from here on down.

x = frexpf( f, &e ); /* f = x * 2**e, 0.5 <= x < 1.0 */
/* If power of 2 is odd, double x and decrement the power of 2. */
if( e & 1 ) {
x = x + x;
e -= 1;
}
>
e >>= 1; /* The power of 2 of the square root. */
>
if( x 1.41421356237 ) {
/* x is between sqrt(2) and 2. */
x = x - 2.0;
y =
((((( -9.8843065718E-4 * x
+ 7.9479950957E-4) * x
- 3.5890535377E-3) * x
+ 1.1028809744E-2) * x
- 4.4195203560E-2) * x
+ 3.5355338194E-1) * x
+ 1.41421356237E0 ;
goto sqdon;
}
>
if( x 0.707106781187 ) {
/* x is between sqrt(2)/2 and sqrt(2). */
x = x - 1.0;
y =
((((( 1.35199291026E-2 * x
- 2.26657767832E-2) * x
+ 2.78720776889E-2) * x
- 3.89582788321E-2) * x
+ 6.24811144548E-2) * x
- 1.25001503933E-1) * x * x
+ 0.5 * x
+ 1.0;
goto sqdon;
}
>
/* x is between 0.5 and sqrt(2)/2. */
x = x - 0.5;
y =
((((( -3.9495006054E-1 * x
+ 5.1743034569E-1) * x
- 4.3214437330E-1) * x
+ 3.5310730460E-1) * x
- 3.5354581892E-1) * x
+ 7.0710676017E-1) * x
+ 7.07106781187E-1;
>
sqdon:
y = ldexpf( y, e ); /* y = y * 2**e */
return( y);
}

--
Ben.

**Bartc** · Jun 27 '08, 07:40 PM

Re: sqrt() double trouble

"jacob navia" <jacob@nospam.c omwrote in message
news:g2aoto$3eh $1@aioe.org...

<snip lots of code>

If floating point hardware is available, then using the standard sqrt()
function, even with conversion from and to float, is going to be faster.

--
Bartc

sqrt() double trouble

sqrt() double trouble

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment