Bisezione vs Netwon

#include <iostream>
#include <cmath>
#include <iomanip>
#include <stdexcept>
#include <fstream>
#include <vector>

using namespace std;

/*******************************************************************************
* CONSTANTS
*******************************************************************************/
const int MAX_ITERATION = 100;
const double DERIVATE_PRECISION = 0.001;

/*******************************************************************************
* Functions to be analyzed
*******************************************************************************/

// This function if is 0 is the solutioon of sqrt(2) - n = 0
double f(double x){
    return x*x - 2;
}

// not used in this file but can be assigned to func in main
double f2(double x){
    return x*x*x*x - x*x;
}

/*******************************************************************************
* utils
*******************************************************************************/

// function that returns the "sign" of the value -1, 0 for 0 or, 1 for positive
template <typename T>
int sgn(T val) {
    return (T(0) < val) - (val < T(0));
}

// contains the result of the function algorithms, it contains the iteration
// count and the resulting number to calculate the error
struct AlgResult{
    int itr_count;
    double result;
};

// type of function pointer f(x)
typedef double (*fptr_t) (double);


// calculate the error given the expected result
// 0 as expected result is a pain tho
double calc_error(double expected_result, double experimental_result){
    double err = (expected_result - experimental_result) / expected_result;
    err = err * 100;
    return err;
}

// calculate the derivate of the function for newton method
double prime_derivative(fptr_t fbase, double n){
    double h = DERIVATE_PRECISION;
    return (fbase(n + h) - fbase(n))/ h;
}

/*******************************************************************************
* Bisection method
*******************************************************************************/

// calculates how many iterations to reach the given tolerance
AlgResult bisect(double lower, double upper, double tolerance, fptr_t fx){
    int itr_count = 0;
    AlgResult bs_res;

    double half = (upper + lower) / 2.0;
    double res;

    do{
        res = fx(half);
        // save the result
        bs_res.result = half;
        bs_res.itr_count = itr_count;

        // proceed with the algorithm
        if( sgn<double>(res) == sgn<double>(f(lower)) ){
            lower = half;
        }
        else{
            upper = half;
        }

        half = (lower + upper) / 2;

        ++itr_count;

        if(itr_count >= MAX_ITERATION){
            throw runtime_error(
                "Function bisect: Max iteration reached without convergence");
        }

    } while(abs(res) > tolerance);

    return bs_res;
}

/*******************************************************************************
* Newthon method
*******************************************************************************/

AlgResult newton(double x0, fptr_t fx, double tolerance){
    int itr_count = 0;
    AlgResult bs_res;

    // condition to calculate next iteration
    bool cond;

    do{

        double x1 = x0 - fx(x0) / prime_derivative(fx, x0);

        // save the results
        bs_res.itr_count = itr_count;
        bs_res.result = x1;

        // calculate if is necessary to go for an other round
        cond = abs(x1 - x0) <= tolerance * abs(x1);

        if(itr_count >= MAX_ITERATION){
            throw runtime_error(
                "Function newton: Max iteration reached without convergence");
        }

        // assign the new x0
        x0 = x1;
        itr_count += 1;

    }while(!cond);

    return bs_res;
}

/*******************************************************************************
* MAIN PROGRAM
*******************************************************************************/

enum AlgType{
    NEWTON = 0,
    BISECT = 1
};

int main()
{
    //------------------ user inputs --------------------------------------

    // define the search interval or the x0 for the newton method
    double lower = 0, upper = 2.0, x0 = 2.0;

    // chose the algorithm (AlgType.NEWTON or AlgType.BISECT)
    AlgType algtype = AlgType::NEWTON;

    // function to be analyzed
    fptr_t func = &f;

    // define the initial tolerance
    // the tolerance step later are calculated as 1/(power of 10)
    double init_tol = 0.1;

    // define the dataset size
    int data_size = 10;

    // define the theoretical value the algorithm should find
    double tvalue = sqrt(2);

    // variable saved to file
    string filename = "data.dat";
    vector<double> x(data_size);
    vector<double> y(data_size);
    vector<int> itr(data_size);

    //------------------ algorithm starts --------------------------------------
    cout << "Algoritm starts" << endl;
    string alg_name =
        (algtype == AlgType::NEWTON)? "Newton method" : "Bisect method";
    cout << "Method chosen: " << alg_name << endl;
    AlgResult br;
    for(int i = 0; i < data_size; i++){
            double tolerance = init_tol / pow(10, i);

            cout << "------------------------------" << endl;
            cout << "tolerance: " << tolerance << endl;

            if(algtype == AlgType::NEWTON){
                br = newton(x0, f, tolerance);
            }
            else{
                br = bisect(lower, upper, tolerance, func);
            }

            cout << "iterations: " << br.itr_count
                 << " result: " << br.result << endl;
            cout << "error: " << calc_error(tvalue, br.result) << "%" << endl;

            // save the steps in the arrays
            x[i] = tolerance;
            y[i] = br.result;
            itr[i] = br.itr_count;
    }

    cout << "---- algorithme finished -----" << endl;

    // save file
    ofstream file;
    file.open(filename, ios::out | ios::binary);
    if(file.is_open()){
        // since under the vector the data are in a contigous memory
        // &first_element gives the address to the memory block is pointed too
        file.write(reinterpret_cast<char*>(&x[0]),   data_size*sizeof(double));
        file.write(reinterpret_cast<char*>(&y[0]),   data_size*sizeof(double));
        file.write(reinterpret_cast<char*>(&itr[0]), data_size*sizeof(int));
        file.close();

        cout << "Data saved in file: " << filename << endl;
    }
    else{
        throw runtime_error("File problem: can't write");
    }

    return 0;
}