Calcolo pi greca metodo Brent-Salamin

[irelimax]

Buongiorno a tutti, devo scrivere un codice in C che calcoli il valore di pi greco tramite la seguente procedura iterativa, definendo una funzione esterna:
$a_0=1, b_0=1/ sqrt(2), p_0=1, t_0=1/4$ e iterando $a_{n+1}= (a_n + b_n)/2, b_{n+1}=sqrt(a_n*b_n), p_{n+1}=2*p_n, t_{n+1}= t_n - p_n(a_n-a_{n+1})^2$

Il codice che ho scritto è il seguente ma c'è qualcosa che non va..penso che il problema sia nel ciclo for ma ho provato a farlo in diversi modi e non viene, potreste aiutarmi??:

#include<conio.h>
#include<iostream>
using namespace std;

double brent ()
{

int i,j,N=10000;
	double y[10000], x[10000];
	double a_0=1, b_0= 1/(sqrt(2.0));
	
	double pi;
	double t_0=1/4, p_0=1;
	double t[10000], p[10000];
	

	for(i=1; i<=N; i++)
	{
	x[0]=a_0;
	y[0]=b_0;
	for(j=1;j<=N; j++)
	{
	x[i]=(x[i-1]+y[j-1])/2;
	y[j]=sqrt((x[i-1])*(y[j-1]));
	
	t[0]=t_0;
	p[0]=p_0;
	p[j]=2*p[j-1];
	t[j]=t[j-1]-p[j-1]*((x[i-1]-x[i])*(x[i-1]-x[i]));
	}
}	

	pi=((x[N]+y[N])*(x[N]+y[N]))/4*t[N];

return (pi);

}

int main()
{
	double PI;

	PI=brent();
	cout<<"The value of PI is "<<PI<<endl;
	getch();
	return 0;
}

[apatriarca]

Sono sul cellulare per cui non leggo molto bene, ma che senso hanno due cicli annidati? Credo tu non abbia compreso il metodo..

[irelimax]

All'inizio infatti avevo fatto il codice con un solo ciclo for ma mi dava un errore (Stack around the variable 'p' was corrupted) quindi ho pensato di fare due cicli for...

[hyoukarou]

È abbastanza incasinato quel codice, qualche tips:

[*:3tgi1q83] nel ciclo for j assume il valore N, quindi x[j] è out of range(che è quello che dicevi ti diceva l'os)[/*:m:3tgi1q83]
[*:3tgi1q83] a t_0, 1/4 lo valuta 0, avresti dovuto scrivere 1/4.0[/*:m:3tgi1q83]
[*:3tgi1q83] il ciclo esterno è praticamente inutile(e dannoso), x[0] = a_0 &co. lo dovresti eseguire una sola volta[/*:m:3tgi1q83]
[*:3tgi1q83] sprechi molta memoria creando quegli array, puoi fare il tutto con qualche valore, un po' come si fa con il classico algoritmo per trovare l'n-esimo numero di fibonacci[/*:m:3tgi1q83]
[*:3tgi1q83] a_0, b_0 e quegli altri valori che non modifichi sarebbe bene mettere costanti(anzi, sarebbe meglio non metterli proprio in questo caso)[/*:m:3tgi1q83]
[*:3tgi1q83] ho fatto un po' di correzioni e mi esce che pi vale 0.163991, probabilmente c'è qualche errore in qualche formula[/*:m:3tgi1q83]
[*:3tgi1q83] cerca di essere il più sintetico possibile, mantieni le cose semplici[/*:m:3tgi1q83][/list:u:3tgi1q83]

buona fortuna

[irelimax]

ok l'ho sistemato e anche a me adesso viene 0.164053...mmmm... ho tolto gli array ed è venuto molto più semplice..
comunque ho controllato la formula e anche wikipedia conferma la formula che ho io...dv è l'intoppo?

[hyoukarou]

Hai dimenticato le parentesi al denominatore di \(\pi\).

[irelimax]

ups....ecco infatti ora funziona...grazie mille per averi aiutata ;P