Błąd liczenia dla wysokich N

0

Cześć.
W poniższym kodzie dla N=1000 obliczanie pi metodą mnożenia w wersji połączonej działa poprawnie, ale dla N=100000 wynik tego obliczenia pojawia się jako -0.(0). Gdzie tkwi błąd? Gdzie szukać rozwiązania?
Poniżej kod programu.
Pozdrawiam i czekam na pomoc.

#include <stdio.h>
#include <stdlib.h>
#include <sys/types.h>
#include <omp.h>
#include <unistd.h>
#include <time.h>
#include <math.h>
#define N 100000

void pi1()
{
    long double pi = (double)4/3;
    long double tmp2;
    long long int i = 2;
    long int limit = N;
#pragma omp parallel for default(none) shared(limit) private(i, tmp2) reduction(*:pi)
//#pragma omp parallel for shared(limit, tmp2) private(i) reduction(*:pi)
    for(i=2; i<limit; ++i)
    {
        tmp2 = (double)(4*i*i)/(4*i*i-1);
        pi = pi * tmp2;
    }
    printf("liczba pi (mnozenie)  = %.30Lf\n", 2.0*pi);
}

void pi2()
{
    long double pi = 0;
    long double tmp1;
    long int i = 0;
    long int limit = N;
#pragma omp parallel for default(none) shared(limit) private(i, tmp1) reduction(+:pi)
    for(i = 0; i<limit; ++i)
    {
        if(i % 2 == 0)
        {
            tmp1 = (double)(1)/(2*i+1);
        }
        else
        {
            tmp1 = (double)(-1)/(2*i+1);
        }
        pi = pi + tmp1;
    }
    printf("liczba pi (dodawanie) = %.30Lf\n", 4.0*pi);
}

void pi3()
{
    long double pi_mnoz = 1;
    long double pi_dod = 1;
    long double tmp1, tmp2;
    long int i = 0;
    long int limit = N;
#pragma omp parallel for default(none) shared(limit) private(i, tmp1, tmp2) reduction(+:pi_dod) reduction(*:pi_mnoz)
    for(i = 1; i<limit; ++i)
    {
        if(i % 2 == 0)
        {
            tmp1 = (double)(1)/(2*i+1);
        }
        else
        {
            tmp1 = (double)(-1)/(2*i+1);
        }

	if (i < 1)
	{
		printf(":D\n");
	}
	else
	{
		tmp2 = (double)(4*i*i)/(4*i*i-1);
	}

        pi_mnoz = pi_mnoz * tmp2;
        pi_dod = pi_dod + tmp1;
    }
    printf("\nwersja polaczona:\n");
    printf("liczba pi (dodawanie) = %.30Lf\n", 4.0*pi_dod);
    printf("liczba pi (mnozenie)  = %.30Lf\n", 2.0*pi_mnoz);
}


int main()
{
    clock_t t1 = time(NULL), t2, t3;
    omp_set_num_threads(8);
    printf("\nN = %d\n", N);

    printf("\noddzielnie:\n");
    pi1();
    pi2();
    t2 = time(NULL);
    printf("=> Liczone oddzielnie: %f\n", difftime(t2, t1));

    pi3();
    t3 = time(NULL);
    printf("=> Liczone razem: %f.\n",difftime(t3, t2));

    return 0;
}
0

Nie sprawdzałem kodu, ale czy masz na uwadze to że double ma ograniczoną precyzje? To są tylko 64 bity...

1

Robisz tak potworny bałagan w kodzie ...
http://ideone.com/gHZYwq

1 użytkowników online, w tym zalogowanych: 0, gości: 1