Greg wrote:
[color=blue]
>
int2str@gmail.com wrote:[color=green]
>> Mark wrote:[color=darkred]
>> > always wondered how a computer might go about doing this...
>> > i'm taking calculus at university... and we went over newton's method.
>> > i found this interesting...
>> >
>> > ------------------------------
>> > #include <cstdlib>
>> > #include <iostream>
>> > #include <cmath>
>> >
>> > using namespace std;
>> >
>> > double f1(double x);
>> > double f2(double x, double y);
>> >
>> > int main(int argc, char *argv[])
>> > {
>> > double x = 5;
>> > cout << f1(x) << endl;
>> > cout << sqrt(x) << endl;
>> >
>> > system("PAUSE");
>> > return EXIT_SUCCESS;
>> > }
>> >
>> > double f1(double x)
>> > {
>> > return f2(x,x/2);
>> > }
>> >
>> > double f2(double x, double y)
>> > {
>> > double z = y - (y*y-x)/(2*x);
>> > if( y == z ) return z;[/color]
>>
>> I thought it was bad to directly compare floats?[/color]
>
> It would depend on why the floats are being compared. Testing two
> different floating point values for equality is a bad idea since the
> values are compared exactly even though the values may not be
> represented exactly. So two values that would be considered equal by
> the programmer may not always compare as equal in the program.
>
> In this case, however, the comparison is fine. One floating value is
> being compared against itself, or more precisely, the result of one
> iteration is being compared against the result of the previous
> iteration. Since the difference between the values produced at each
> iteration becomes smaller and smaller, at some point the difference
> will have become so minute that the floating point variable can no
> longer detect it.[/color]
This reasoning is based on calculus and does not apply to floating point
numbers of limitedprecision.
[color=blue]
> At that point, the two floating point values will compare equal (since
> the new value appears unchanged from the previous one). Since there is
> no point in any further iterations when that happens, the function
> stops iterating and returns with its result.[/color]
Nope:
#include <iostream>
#include <iomanip>
#include <limits>
#include <cmath>
double my_sqrt( const double & x )
{
double y = 0;
double z = 1;
while ( y != z )
{
z = y;
y = y - (y * y - x) / (2 * x);
std::cout << std::setprecision(17)
<< y << " "
<< z << '\n';
}
return y;
}
double my_sqrt_b ( const double & x )
{
double root = 0;
double quot = x;
while ( std::abs( root - quot ) >
std::numeric_limits<double>::epsilon() ) {
root = ( root + quot ) / 2;
quot = x / root;
std::cout << root << " " << quot << '\n';
}
return ( root );
}
int main ( void ) {
std::cout << my_sqrt( 0.5 ) << '\n';
}
prints on my machine:
0.5 0
0.75 0.5
0.6875 0.75
0.71484375 0.6875
0.7038421630859375 0.71484375
0.70844837254844606 0.7038421630859375
0.7065492759819042 0.70844837254844606
0.7073373965913512 0.7065492759819042
0.70701120397472073 0.7073373965913512
0.70714636142893661 0.70701120397472073
0.70709038494675236 0.70714636142893661
0.70711357246260598 0.70709038494675236
0.70710396810177689 0.70711357246260598
0.70710794639649821 0.70710396810177689
0.70710629853942519 0.70710794639649821
0.70710698110529846 0.70710629853942519
0.70710669837744955 0.70710698110529846
0.70710681548719212 0.70710669837744955
0.70710676697875419 0.70710681548719212
0.70710678707160801 0.70710676697875419
0.70710677874887562 0.70710678707160801
0.70710678219626433 0.70710677874887562
0.70710678076830913 0.70710678219626433
0.70710678135978755 0.70710678076830913
0.7071067811147892 0.70710678135978755
0.7071067812162708 0.7071067811147892
0.70710678117423575 0.7071067812162708
0.70710678119164727 0.70710678117423575
0.70710678118443515 0.70710678119164727
0.70710678118742254 0.70710678118443515
0.70710678118618508 0.70710678118742254
0.70710678118669767 0.70710678118618508
0.70710678118648529 0.70710678118669767
0.70710678118657333 0.70710678118648529
0.7071067811865368 0.70710678118657333
0.70710678118655201 0.7071067811865368
0.70710678118654569 0.70710678118655201
0.70710678118654824 0.70710678118654569
0.70710678118654724 0.70710678118654824
0.70710678118654768 0.70710678118654724
0.70710678118654746 0.70710678118654768
0.70710678118654757 0.70710678118654746
0.70710678118654746 0.70710678118654757
0.70710678118654757 0.70710678118654746
0.70710678118654746 0.70710678118654757
0.70710678118654757 0.70710678118654746
0.70710678118654746 0.70710678118654757
0.70710678118654757 0.70710678118654746
0.70710678118654746 0.70710678118654757
0.70710678118654757 0.70710678118654746
0.70710678118654746 0.70710678118654757
0.70710678118654757 0.70710678118654746
0.70710678118654746 0.70710678118654757
0.70710678118654757 0.70710678118654746
0.70710678118654746 0.70710678118654757
0.70710678118654757 0.70710678118654746
0.70710678118654746 0.70710678118654757
As you can see, you run the risk that y and z enter a cycle. At this point,
their difference does not get any smaller.
I just know enough about numerical analysis to stay away from these
problems. Dealing with floats and doubles is *hard* and best left to
experts who can cope with a total break down of Calculus based intuition
and reasoning.
Best
Kai-Uwe Bux
