main.cpp

#include <iostream>
#include <vector>
#include <map>
#include <random>
#include <chrono>
#include <cmath>
#include <iomanip>
#include <string>

#include "include/DoubleDouble.h"

using namespace doubledouble;
using dd = DoubleDouble;

// ================================================================
// Configuration
// ================================================================
// Available tests: "prod", "compprod", "ddprod", "blockprod", "logprod", "all"
std::string TEST_MODE = "all";   // change or override via CLI
int N = 100000;                  // number of factors
int M = 1;                       // number of experiments
int SEED = 42;                   // random seed
bool VERBOSE = false;            // print extra info
// ================================================================

// ------------------------------------------------
// compute the relative error
// ------------------------------------------------
inline double rel_error(long double ref, double val) {
    return std::fabs((val - ref) / ref);
}

// ------------------------------------------------
// TwoProduct using FMA if available (C++17 fma)
// ------------------------------------------------
inline void TwoProductFMA(double a, double b, double& x, double& y) {
    x = a * b;
    y = std::fma(a, b, -x);
}

// ------------------------------------------------
// Classic product algorithm
// ------------------------------------------------
double Prod(const std::vector<double>& a) {
    double p = a[0];
    for (size_t i = 1; i < a.size(); ++i)
        p = p * a[i];
    return p;
}

// ------------------------------------------------
// Compensated product algorithm
// ------------------------------------------------
// 
// ------------------------------------------------
double CompProd(const std::vector<double>& a) {
    double p = a[0];
    double e = 0.0;
    for (size_t i = 1; i < a.size(); ++i) {
        double pi, delta;
        TwoProductFMA(p, a[i], pi, delta);
        e = std::fma(e, a[i], delta);
        p = pi;
    }
    return std::nextafter(p + e, 0.0);
}

// ------------------------------------------------
// Double Double Product
// ------------------------------------------------
// Use the double double library to compute the product
// ------------------------------------------------
double DDProd(const std::vector<double>& a) {
    dd p = dd(a[0]);
    for (size_t i = 1; i < a.size(); ++i) {
        p *= a[i];
    }
    return p.upper + p.lower;
}

// ------------------------------------------------
// Block‐wise high‐dynamic‐range accumulator
// ------------------------------------------------
// multiply elements in blocks to limit propagation of rounding
// errors and reduce underflow/overflow. Each block is computed in
// extended precision (long double) and combined progressively
// ------------------------------------------------
double BLOCKProd(const std::vector<double>& a, size_t blockSize = 1024) {
    const size_t n = a.size();
    if (n == 0) return 1.0;
    size_t i = 0;

    long double totalMant = 1.0L;
    int64_t totalExp = 0;  // track exponent shifts (base 2)

    while (i < n) {
        size_t end = std::min(i + blockSize, n);
        long double blockMant = 1.0L;
        int64_t blockExp = 0;

        for (size_t j = i; j < end; ++j) {
            long double val = a[j];
            blockMant *= val;
            // Normalize blockMant to avoid overflow/underflow
            if (blockMant > 1e300L) {
                blockMant *= 1.0L / 1e300L;
                blockExp += 300;
            }
            else if (blockMant != 0.0L && std::fabsl(blockMant) < 1e-300L) {
                blockMant *= 1e300L;
                blockExp -= 300;
            }
        }

        // Combine block into total
        totalMant *= blockMant;
        totalExp += blockExp;

        // Normalize totalMant
        if (std::fabsl(totalMant) > 1e300L) {
            totalMant *= 1.0L / 1e300L;
            totalExp += 300;
        }
        else if (totalMant != 0.0L && std::fabsl(totalMant) < 1e-300L) {
            totalMant *= 1e300L;
            totalExp -= 300;
        }

        i = end;
    }

    // Convert back to double: totalMant × 2^{totalExp}
    double result = static_cast<double>(totalMant) * std::ldexp(1.0, totalExp);
    return result;
}

double PairwiseSum(const std::vector<double>& a) {
    if (a.empty()) {
        return 0.0;
    }
    if (a.size() == 1) {
        return a[0];
    }
    size_t mid = a.size() / 2;
    std::vector<double> left(a.begin(), a.begin() + mid);
    std::vector<double> right(a.begin() + mid, a.end());
    return PairwiseSum(left) + PairwiseSum(right);
}

// ------------------------------------------------
// Logarithmic product algorithm
// use log and exp property to compute a sum instead
// ------------------------------------------------
double LogProd(const std::vector<double>& a) {
    int sign = 1;
    std::vector<double> log_values;
    log_values.reserve(a.size());

    for (double x_i : a) {
        if (x_i == 0.0) {
            return 0.0;
        }
        if (x_i < 0) {
            sign = -sign;
        }
        log_values.push_back(std::log(std::fabs(x_i)));
    }

    double log_sum = PairwiseSum(log_values);
    return sign * std::exp(log_sum);
}

// ------------------------------------------------
// Reference computed in long double (≈ 80–128 bits)
// ------------------------------------------------
long double reference_prod(const std::vector<double>& a) {
    long double p = static_cast<long double>(a[0]);
    for (size_t i = 1; i < a.size(); ++i)
        p *= static_cast<long double>(a[i]);
    return p;
}

// ------------------------------------------------
// Random vector generator
// ------------------------------------------------
std::vector<double> generate_vector(int n, int seed = 42) {
    std::mt19937_64 gen(seed);
    std::uniform_real_distribution<double> dist(0.5, 1.5);
    std::vector<double> a(n);
    for (auto& x : a) x = dist(gen);
    return a;
}

// ------------------------------------------------
// A struct to hold test results
// ------------------------------------------------
struct TestResult {
    std::string name;
    double total_time = 0.0;
    double total_rel_error = 0.0;
    int runs = 0;

    void record(double time, double rel_error) {
        total_time += time;
        total_rel_error += rel_error;
        runs++;
    }

    void print_avg() const {
        if (runs == 0) return;
        double avg_time = total_time / runs;
        double avg_rel_error = total_rel_error / runs;
        std::cout << std::left << std::setw(10) << name
                  << " | avg_time = " << std::setw(10) << avg_time << " s"
                  << " | avg_rel_error = " << std::scientific << avg_rel_error << "\n";
    }
};

// ------------------------------------------------
// CLI Parsing (optional overrides)
// ------------------------------------------------
void parse_args(int argc, char** argv) {
    for (int i = 1; i < argc; ++i) {
        std::string arg = argv[i];
        if (arg == "--n" && i + 1 < argc) N = std::stoi(argv[++i]);
        else if (arg == "--m" && i + 1 < argc) M = std::stoi(argv[++i]);
        else if (arg == "--mode" && i + 1 < argc) TEST_MODE = argv[++i];
        else if (arg == "--seed" && i + 1 < argc) SEED = std::stoi(argv[++i]);
        else if (arg == "--verbose") VERBOSE = true;
    }
}

// ------------------------------------------------
// Main
// ------------------------------------------------
int main(int argc, char** argv) {
    parse_args(argc, argv);
    std::cout << "Running test mode: " << TEST_MODE
              << " with n = " << N << " over m = " << M << " experiments.\n";

    std::map<std::string, TestResult> results;
    if (TEST_MODE == "prod" || TEST_MODE == "all") results["Prod"].name = "Prod";
    if (TEST_MODE == "compprod" || TEST_MODE == "all") results["CompProd"].name = "CompProd";
    if (TEST_MODE == "ddprod" || TEST_MODE == "all") results["DDProd"].name = "DDProd";
    if (TEST_MODE == "blockprod" || TEST_MODE == "all") results["BLOCKProd"].name = "BLOCKProd";
    if (TEST_MODE == "logprod" || TEST_MODE == "all") results["LogProd"].name = "LogProd";

    std::mt19937_64 exp_seed_gen(SEED);

    for (int i = 0; i < M; ++i) {
        int current_seed = exp_seed_gen();
        auto data = generate_vector(N, current_seed);
        long double ref = reference_prod(data);

        auto run_and_record = [&](const std::string& name, double (*func)(const std::vector<double>&)) {
            if (results.find(name) == results.end()) return;

            auto start = std::chrono::high_resolution_clock::now();
            double result = func(data);
            auto end = std::chrono::high_resolution_clock::now();
            double t = std::chrono::duration<double>(end - start).count();
            double err = rel_error(ref, result);
            results[name].record(t, err);
        };

        run_and_record("Prod", Prod);
        run_and_record("CompProd", CompProd);
        run_and_record("DDProd", DDProd);
        run_and_record("BLOCKProd", [](const std::vector<double>& v){ return BLOCKProd(v, 1024); });
        run_and_record("LogProd", LogProd);
    }

    std::cout << "\n--- Average Results over " << M << " experiments ; n = " << N << " ---\n";
    for (const auto& pair : results) {
        pair.second.print_avg();
    }

    return 0;
}