On Nov 7, 6:40*pm, pete <(E-Mail Removed)> wrote:

> ImpalerCore wrote:

>

> > On Nov 7, 1:40 pm, James Dow Allen <(E-Mail Removed)> 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.