#include <iostream>
#include <cmath>

using namespace std;

// Definición de la función f(x, y, z) que queremos integrar
double f(double x, double y, double z) {
    return sin(x) * cos(y) * exp(z);
}

// Función para calcular la aproximación numérica de la integral triple
double calcularIntegralTriple(double x0, double x1, double y0, double y1, double z0, double z1, int n) {
    double hx = (x1 - x0) / n;
    double hy = (y1 - y0) / n;
    double hz = (z1 - z0) / n;

    double sum = 0.0;

    for (int i = 0; i <= n; ++i) {
        double x = x0 + i * hx;
        for (int j = 0; j <= n; ++j) {
            double y = y0 + j * hy;
            for (int k = 0; k <= n; ++k) {
                double z = z0 + k * hz;

                // Calcula el valor de f(x, y, z) en el punto (x, y, z)
                double funcValue = f(x, y, z);

                // Añade la contribución al sumatorio
                if (i == 0 || i == n) {
                    if (j == 0 || j == n) {
                        if (k == 0 || k == n) {
                            sum += funcValue;
                        } else {
                            sum += 2 * funcValue;
                        }
                    } else {
                        if (k == 0 || k == n) {
                            sum += 2 * funcValue;
                        } else {
                            sum += 4 * funcValue;
                        }
                    }
                } else {
                    if (j == 0 || j == n) {
                        if (k == 0 || k == n) {
                            sum += 2 * funcValue;
                        } else {
                            sum += 4 * funcValue;
                        }
                    } else {
                        if (k == 0 || k == n) {
                            sum += 4 * funcValue;
                        } else {
                            sum += 8 * funcValue;
                        }
                    }
                }
            }
        }
    }

    // Multiplica por hx, hy, hz y 1/8 para obtener el valor aproximado de la integral
    double volume = hx * hy * hz;
    double integral = sum * volume / 8.0;

    return integral;
}

int main() {
    // Definir los límites de integración y el número de divisiones por dimensión (n)
    double x0 = 0.0, x1 = 1.0;
    double y0 = 0.0, y1 = 1.0;
    double z0 = 0.0, z1 = 1.0;
    int n = 100;  // Número de divisiones por dimensión

    // Calcular la aproximación numérica de la integral triple
    double resultado = calcularIntegralTriple(x0, x1, y0, y1, z0, z1, n);

    // Multiplicar por 0.125 como indica la fórmula dada
    resultado *= 0.125;

    // Mostrar el resultado
    cout << "El resultado de la integral triple aproximada es: " << resultado << endl;

    return 0;
}