/* Pi berechnen mit dem Chudnovsky-Algorithmus
 */

void setup() {
  // put your setup code here, to run once:
  Serial.begin(115200);
  Serial.println("Hello, ESP32!");
}

void loop() {
  double Pi = gauss_legendre(1);
  Serial.printf("Pi(1) = %15.14f \n", Pi);
  delay(1000); // this speeds up the simulation
  Pi = gauss_legendre(2);
  Serial.printf("Pi(2) = %15.14f \n", Pi);
  delay(1000); // this speeds up the simulation
  Pi = gauss_legendre(3);
  Serial.printf("Pi(3) = %20.19f \n", Pi);
  delay(1000); // this speeds up the simulation
  Pi = gauss_legendre(4);
  Serial.printf("Pi(4) = %20.19lf \n", Pi);
  delay(1000); // this speeds up the simulation
}

double chudnovsky(int iterations) {
    double sum = 0.0;
    double k = 0.0;
    double a = 13591409.0;
    double b = 545140134.0;
    double c = 640320.0;
    
    for (int i = 0; i < iterations; i++) {
        double numerator = factorial(6*k) * (a + b*k);
        double denominator = factorial(3*k) * pow(factorial(k), 3) * pow(c, 3*k);
        sum += numerator / denominator * pow(-1, k);
        k += 1.0;
    }
    
    return 1 / (12 * sum);
}

double gauss_legendre(int iterations) {
    double a = 1.0;
    double b = 1.0 / sqrt(2);
    double t = 1.0 / 4.0;
    double p = 1.0;
    
    for (int i = 0; i < iterations; i++) {
        double a_next = (a + b) / 2;
        double b_next = sqrt(a * b);
        double t_next = t - p * pow(a - a_next, 2);
        double p_next = 2 * p;
        
        a = a_next;
        b = b_next;
        t = t_next;
        p = p_next;
    }
    
    return pow(a + b, 2) / (4 * t);
}

unsigned long long factorial(int n) {
    unsigned long long result = 1;
    for (int i = 2; i <= n; i++) {
        result *= i;
    }
    return result;
}