#include "specialAlgorithm.hpp"
#include <iomanip>

arma::vec* special_algorithm(
    double sub_sig,
    double main_sig,
    double sup_sig,
    arma::vec* g_vec
)
{
    int n = g_vec->n_elem;
    arma::vec diag = arma::vec(n);

    for (int i = 1; i < n; i++) {
        // Calculate values for main diagonal based on indices
        diag(i-1) = (double)(i+1) / i;
        (*g_vec)(i) += (*g_vec)(i-1) / diag(i-1);
    }
    // The last element in main diagonal has value (i+1)/i = (n+1)/n
    (*g_vec)(n-1) /= (double)(n+1) / (n);

    for (int i = n-2; i >= 0; i--) {
        (*g_vec)(i) = ((*g_vec)(i) + (*g_vec)(i+1))/ diag(i);
    }

    return g_vec;
}

void special_algorithm_main() 
{
    arma::vec g_vec, *v_vec;
    std::ofstream ofile;
    int steps;
    double sub_sig, main_sig, sup_sig, step_size;

    for (int i = 0; i < 6; i++) {
        steps = std::pow(10, i+1);
        step_size = 1./(double) steps;
        g_vec.resize(steps - 1);

        v_vec = special_algorithm(sub_sig, main_sig, sup_sig, &g_vec);

        ofile.open("output/special_algorithm_" + std::to_string(steps) + ".txt");

        for (int j=0; j < v_vec->n_elem; j++) {
            ofile << std::setprecision(4) << std::scientific << step_size*(j+1) << ","
                << std::setprecision(4) << std::scientific << (*v_vec)(j) << std::endl;
        }
        ofile.close();
    }
}