\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] ) );
}
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.
}