#include "mylib.h"
#include "mayer_primes.h"

#include <cassert>          // assert
#include <vector>
#include <iostream>
#include <cmath>
#include <tuple>
#include <string>
#include <algorithm>
#include <fstream>
#include <stdexcept>
#include <random>
#include <list>
using namespace std;

// -------------- Task A --------------

tuple<double, double, double> means0(double a, double b, double c){
    double arith = (a+b+c) / 3;
    double geo = pow((a*b*c), 1.0/3);
    double harm = 3 / ((1.0/a) + (1.0/b) + (1.0/c));
    return make_tuple(arith, geo, harm);
}

tuple<double, double, double> means(const vector<double>& v){
    size_t n = v.size();
    double sum = 0;
    double logsum = 0;
    double invsum = 0;

    for (size_t i = 0; i<n; ++i){
        sum += v[i];
        logsum += log(v[i]);
        invsum += 1.0/v[i];
    }

    double arith = sum / static_cast<double>(n);
    double geo = exp(1.0/static_cast<double>(n) * logsum);
    double harm = static_cast<double>(n) / invsum;
    return make_tuple(arith, geo, harm);
}

// -------------- Task B --------------

void fill_vector(istream& istr, vector<double>& v)
{
    double d=0;
    while ( istr >> d) v.push_back(d); // Einlesen
    if (!istr.eof())
    { // Fehlerbehandlung
        cout << " Error handling \n";
        if ( istr.bad() )  throw runtime_error("Schwerer Fehler in istr");
        if ( istr.fail() )   // Versuch des Aufraeumens
        {
            cout << " Failed in reading all data.\n";
            istr.clear();
        }
    }
    v.shrink_to_fit();                 // C++11
    return;
}


void read_vector_from_file(const string& file_name, vector<double>& v)
{
    ifstream fin(file_name);           // Oeffne das File im ASCII-Modus
    if( fin.is_open() )                // File gefunden:
    {
       v.clear();                      // Vektor leeren
       fill_vector(fin, v);
    }
    else                               // File nicht gefunden:
    {
        cout << "\nFile " << file_name << " has not been found.\n\n" ;
        assert( fin.is_open() && "File not found." );       // exeption handling for the poor programmer
    }

 return;
}

void write_vector_to_file(const string& file_name, const vector<double>& v)
{
    ofstream fout(file_name);          // Oeffne das File im ASCII-Modus
    if( fout.is_open() )
    {
       for (size_t k=0; k<v.size(); ++k)
       {
          fout << v.at(k) << endl;
       }
    }
    else
    {
        cout << "\nFile " << file_name << " has not been opened.\n\n" ;
        assert( fout.is_open() && "File not opened."  );         // exeption handling for the poor programmer
    }

 return;
}

// -------------- Task C --------------

long int sum_of_spec(int n)
{
    long int sum = 0;
    for (long int i=1; i<n+1; i++){
        if (i % 3 == 0 ||  i % 5 == 0){
            sum += i;
        }
    }
    return sum;
}

double formula(int n)
{
    double div_by_3 = floor( n / 3.0 );
    double div_by_5 = floor( n / 5.0 );
    double div_by_15 = floor( n / 15.0 );

    double S_3 = 3.0 * (div_by_3*((div_by_3+1)/2.0));
    double S_5 = 5.0 * (div_by_5*((div_by_5+1)/2.0));
    double S_15 = 15.0 * (div_by_15*((div_by_15+1)/2.0));

    double sum = S_3 + S_5 - S_15;

    return sum;
}

// -------------- Task E --------------

void insert_into_vector(vector<int>& vec, int n) {
    random_device rd;                       // random device
    mt19937 gen(rd());                      // seed
    uniform_int_distribution<> dist(1, n);  // define range

    for (int i = 0; i < n; ++i) {
        int rand_num = dist(gen);
        auto pos = lower_bound(vec.begin(), vec.end(), rand_num);
        vec.insert(pos, rand_num);
    }
    assert(is_sorted(vec.begin(), vec.end()));
}

void insert_into_list(list<int>& lst, int n) {
    random_device rd;
    mt19937 gen(rd());
    uniform_int_distribution<> dist(1, n);

    for (int i = 0; i < n; ++i) {
        int rand_num = dist(gen);
        auto pos = lower_bound(lst.begin(), lst.end(), rand_num);
        lst.insert(pos, rand_num);
    }
    assert(is_sorted(lst.begin(), lst.end()));
}

// -------------- Task F --------------

int single_goldbach(int k) {
    const vector<int> primes = get_primes(k);
    int count = 0;

    for (size_t i = 0; i < primes.size(); i++) {
        for (size_t j = i; j < primes.size(); j++) {
            if (primes[i] + primes[j] == k) {
                count++;
            }
        }
    }

    return count;
}

vector<int> count_goldbach(int n) {
    const vector<int> primes = get_primes(n);
    vector<int> counts(n+1);

    for (size_t i = 1; i < primes.size(); i++) {
        for (size_t j = i; j < primes.size(); j++) {
            int sum = primes[i] + primes[j];
            if (sum <= n) {
                counts[sum]++;
            }
        }
    }

    return counts;
}


void print_decomps(int k) {
    const vector<int> primes = get_primes(k);
    cout << "\nDecompositions for k = " << k << ": ";

    for (size_t i = 0; i < primes.size(); i++) {
        for (size_t j = i; j < primes.size(); j++) {
            if (primes[i] + primes[j] == k) {
                cout << primes[i] << " + " << primes[j] << ", ";
            }
        }
    }
    cout << endl;
}

// -------------- Task G --------------


double sigmoid(double x)
{
    return 1.0 / (1.0 + std::exp(-x));
}

double DenseMatrix::sigmoid(double x)
{
    return 1.0 / (1.0 + exp(-x));
}

DenseMatrix::DenseMatrix(size_t n, size_t m) : rows(n), cols(m), matrix(n*m)
{
    size_t nm = max(n, m);
    for (size_t i = 0; i < rows; i++)
    {
        for (size_t j = 0; j < cols; j++)
        {
            double x_i = 10.0*static_cast<double>(i)/(nm-1) - 5.0;
            double x_j = 10.0*static_cast<double>(j)/(nm-1) - 5.0;
            matrix[i*cols + j] = sigmoid(x_i) * sigmoid(x_j);
        }
    }
}

vector<double> DenseMatrix::Mult(const vector<double>& vec) const
{
    assert(vec.size() == cols);
    vector<double> result(rows, 0.0);

    for (size_t i = 0; i < rows; ++i)
    {
        for (size_t j = 0; j < cols; ++j)
        {
            result[i] += matrix[i*cols + j] * vec[j];
        }
    }
    return result;
}

vector<double> DenseMatrix::MultT(const vector<double>& vec) const
{
    assert(vec.size() == rows);
    vector<double> result(cols, 0.0);

    for (size_t i = 0; i < cols; ++i)
    {
        for (size_t j = 0; j < rows; ++j)
        {
            result[i] += matrix[j*cols + i] * vec[j];
        }
    }
    return result;
}

void DenseMatrix::print() const {
    size_t count(0);
    cout.precision(5);
    for (double val : matrix) {
        printf("%.6f ", val);
        
        count++;
        if(count%cols == 0){cout << endl;}
    }
}

// #################################################

DenseMatrix2::DenseMatrix2(vector<double> const u, vector<double> const v) : rows(u.size()), cols(v.size()), u_(u), v_(v) {}

vector<double> DenseMatrix2::Mult(const vector<double>& vec) const
{
    assert(vec.size() == cols);
    vector<double> result(rows, 0.0);

    double scalar(0.0);
    for (size_t i = 0; i < vec.size(); i++)
    {
        scalar += v_[i] * vec[i];
    }
    for (size_t i = 0; i < u_.size(); i++)
    {
        result[i] = scalar*u_[i];
    }

    return result;
}

vector<double> DenseMatrix2::MultT(const vector<double>& vec) const
{
    assert(vec.size() == rows);
    vector<double> result(cols, 0.0);

    double scalar(0.0);
    for (size_t i = 0; i < u_.size(); i++)
    {
        scalar += u_[i] * vec[i];
    }
    for (size_t i = 0; i < v_.size(); i++)
    {
        result[i] = scalar*v_[i];
    }

    return result;
}


void task_a();
void task_b();
void task_c();
void task_d();
void task_e();
void task_f();
void task_g();