4
\$\begingroup\$

I wrote a small program for fun to try to prove Pi by taking a certain precision and radius and using it to calculate the area of the circle. My method should be giving me an area that is just slightly larger than the actual area of the circle, approaching the actual area as the precision gets smaller.

What is being calculated

The funny thing is that this area should always be larger than the actual area, but when I set the precision to be less than 1.0e-9L (say 1.0e-10L), I get a result smaller than 3.141592654. Can you look at my code and tell me where I messed up if anywhere?

#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <cmath>
#include <future>

using namespace std;

namespace Settings{
    register constexpr long double precision = 1.0e-10L;
    register constexpr long double radius = 1.0e0L;
}

// Make the math cleaner
inline long double SQUARE(register long double Val) { return Val * Val; }

// Main calculate function
long double get_Area(register long double radius = Settings::radius);

int main()
{
    time_t benchmark = time(0);
    register long double area = get_Area();
    register long double pi = area / SQUARE(Settings::radius);

    cout << "Pi should be around:" << setprecision(50) << pi << endl << endl;
    cout << "It took " << time(0) - benchmark << " to complete the calculation" << endl;
    system("pause");

    return 0;
}

long double get_Area(register long double radius)
{
    register long double area = 0.0e0L, max_Area = SQUARE(radius);

    // Allow for threading
    register auto part = [](register long double start_point, register long double end_point, register long double m_Area){
        register long double area = 0.0L, y = 0.0L;

        for ( register long double x = start_point; x < end_point; x += Settings::precision ){
            y = sqrt(m_Area - SQUARE(x));
            area += y * Settings::precision;
        }

        return area;
    };

    // Divide the problem into eight parts since my pc has eight threads
    future<long double> A1 = async(launch::async, part,                   0.0L, (1.0L / 8.0L) * radius, max_Area);
    future<long double> A2 = async(launch::async, part, (1.0L / 8.0L) * radius, (2.0L / 8.0L) * radius, max_Area);
    future<long double> A3 = async(launch::async, part, (2.0L / 8.0L) * radius, (3.0L / 8.0L) * radius, max_Area);
    future<long double> A4 = async(launch::async, part, (3.0L / 8.0L) * radius, (4.0L / 8.0L) * radius, max_Area);
    future<long double> A5 = async(launch::async, part, (4.0L / 8.0L) * radius, (5.0L / 8.0L) * radius, max_Area);
    future<long double> A6 = async(launch::async, part, (5.0L / 8.0L) * radius, (6.0L / 8.0L) * radius, max_Area);
    future<long double> A7 = async(launch::async, part, (6.0L / 8.0L) * radius, (7.0L / 8.0L) * radius, max_Area);
    future<long double> A8 = async(launch::async, part, (7.0L / 8.0L) * radius,                 radius, max_Area);

    area = A1.get() + A2.get() + A3.get() + A4.get() + A5.get() + A6.get() + A7.get() + A8.get();

    return area * 4.0L;
}
\$\endgroup\$
1
  • \$\begingroup\$ Suggestion: set the precision to 1/4 of the radius. Pick nice numbers for the radius and do the math by hand. Then run your program and see if you get an exact match \$\endgroup\$
    – Kam
    Apr 24, 2016 at 16:06

2 Answers 2

6
\$\begingroup\$

I see a number of things that could help you improve your program.

Don't abuse using namespace std

Putting using namespace std at the top of every program is a bad habit that you'd do well to avoid. Know when to use it and when not to (as when writing include headers). In this particular case, I happen to think it's perfectly appropriate because it's a single short program that and not a header. Some people seem to think it should never be used under any circumstance, but my view is that it can be used as long as it is done responsibly and with full knowledge of the consequences.

Don't use system("pause")

There are two reasons not to use system("cls") or system("pause"). The first is that it is not portable to other operating systems which you may or may not care about now. The second is that it's a security hole, which you absolutely must care about. Specifically, if some program is defined and named cls or pause, your program will execute that program instead of what you intend, and that other program could be anything. First, isolate these into a seperate functions cls() and pause() and then modify your code to call those functions instead of system. Then rewrite the contents of those functions to do what you want using C++. For example:

void pause() {
    getchar();
}

Don't use register

As of C++11, the use of the register keyword is deprecated. In C++17, it will be gone, so it's wise to avoid its use.

Be careful with floating point

Adding a small number to a large number using floating point is risky because in certain cases, as potentially with your for loop within get_Area(). Specifically, with floating point numbers, adding a large number to a small number can lead to loss of precision. In extreme cases, if a is very large and b is very small, we can have a + b = a which is obviously incorrect for non-zero values of b. (See the oft-cited and still excellent What Every Computer Scientist Should Know About Floating-Point Arithmetic, by David Goldberg for a readable technical discussion for why this is so.) I should mention that it will not happen with this particular code with these particular values, but it's a potential problem that we programmers should all understand.

Consider using portable scalability

The choice of 8 parts is just fine, but you might consider using std::thread::hardware_concurrency for a portable way to scale on multiprocessor systems.

Don't repeat yourself

Rather than calculating the low and high values for each slice manually, and naming each std::future by hand, have the computer do it instead:

const unsigned slices = thread::hardware_concurrency();
vector<future<long double>> A(slices);
const long double increment = radius / slices;
long double lo = 0;
long double hi = increment;
for (unsigned i=0; i<slices; ++i) {
    A[i] =  async(launch::async, part, lo, hi, max_Area);
    lo += increment;
    hi += increment;
}
return std::accumulate(A.begin(), A.end(), 0.0L, 
    [](long double accum, std::future<long double> &f){ 
        return accum + f.get(); 
    }) * 4;

Extend the accuracy

Since you're using a simple Riemann sum, on the quarter circle you've chosen the sum will indeed always be greater than the actual value. However, we can very simply take advantage of the symmetry of the problem and have that error mostly cancel out by summing the entire top half of the circle. That is, instead of $$\pi = 4\int_0^{r} \sqrt{r^2 - x^2}\delta x$$ we can instead compute $$\pi = 2\int_{-r}^{r} \sqrt{r^2 - x^2}\delta x$$.

Omit return 0

When a C++ program reaches the end of main the compiler will automatically generate code to return 0, so there is no reason to put return 0; explicitly at the end of main.

\$\endgroup\$
5
  • \$\begingroup\$ KK, aside from the "using namespace std" I made all of those changes, especially the one for the portable multi-threading. Thank you, I love that, that is amazing! Also, I checked and my compiler gives me crazy precision on a long double will give me a precision of 50 digits, which is way more than enough for what I'm doing. It takes 11 seconds for me to run the program with a precision of 1.0e-10 and the execution time goes up an order of magnitude when the precision goes down an order of magnitude. Any idea why I'm getting less than 3.141592654? \$\endgroup\$
    – user103722
    Apr 24, 2016 at 21:56
  • \$\begingroup\$ For the latter question, no, I don't know why it's yielding a smaller number. That's not what I observe here (64-bit Linux, gcc 5.3.1). I get 3.1415926545520452511754683211364636008511297404766082763671875 \$\endgroup\$
    – Edward
    Apr 24, 2016 at 22:05
  • \$\begingroup\$ Might be on my end then....Using gcc4.9.2 on windows targeting 32-bit, maybe it's a glitch in the compiler? \$\endgroup\$
    – user103722
    Apr 24, 2016 at 23:13
  • \$\begingroup\$ By the way, you can portably find out how much precision is available like this: const unsigned digits = numeric_limits<long double>::digits; and then use setprecision(digits). On my machine, digits = 64. \$\endgroup\$
    – Edward
    Apr 24, 2016 at 23:28
  • \$\begingroup\$ Oops, I mean ::digits10. \$\endgroup\$
    – Edward
    Apr 25, 2016 at 2:01
1
\$\begingroup\$

Instead of a circle consider just a linear function f(x)=1−x. This slopes down to the right just as the circle, in order to make the sum “guaranteed” larger than the actual area. But it's simpler.

And instead of 64-bit or 80-bit floating point consider just a 3-digit decimal floating point. The integral of f(x)=1−x is F(x)=xx2/2+C, and evaluated from 0 to 1 gives the expected area F(1)-F(0) = 1/2 − 0 = 0.5. So the sum “should” be slightly larger.

But something's wring – no, wrong!

#include <assert.h>
#include <iostream>
#include <iomanip>
#include <math.h>       // floor
using namespace std;

auto with_3_digits( const double x )
    -> double
{
    assert( 0 <= x );
    if( x == 0 ) { return 0; }
    double scaling = 1;
    double r;
    if( x < 1 )
    {
        for( r = x; r < 1; r *= 10 ) { scaling *= 10; }
    }
    else
    {
        for( r = x; r > 10; r /= 10 ) { scaling /= 10; }
    }
    return floor( 100*r )/(100*scaling);
}

auto main()
    -> int
{
    cout << setprecision( 7 ) << fixed;
    for( int delta_scale = 1; delta_scale < 1000; delta_scale *= 10 )
    {
        const double delta = with_3_digits( (1.0/3.0)/delta_scale );
        double sum = 0;
        double x;
        for( x = 0; x < 1; x = with_3_digits( x + delta ) )
        {
            sum = with_3_digits( sum + x*delta );
        }
        cout << "delta = " << delta << ", final x = " << x << ", sum = " << sum << endl;
    }
}
delta = 0.3330000, final x = 1.3300000, sum = 0.6630000
delta = 0.0333000, final x = 1.0200000, sum = 0.5030000
delta = 0.0033300, final x = 1.0000000, sum = 0.4540000

And it's worse: if one attempts delta value 0.000333, then the program just hangs. It gets nowhere. Because then the delta is too small to make any change once x gets up to 0.1.

So, this is your problem: an accumulated error of limited precision floating point addition.

And one possible solution is to instead divide and conquer, by trying to structure the sums so that roughly equal magnitude numbers are added.

But a better solution is to find a series for pi that converges much faster.


Regarding the coding, register is no-no, threads are premature optimization (you need to measure whether they help: an 8-time speedup of very little can be less than overhead in setting up threads), and system( "pause" ) is both unnecessary, annoying and non-portable.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.