mnozenie odporne na przepelnienie

#include <iostream>
#include <vector>
#include <algorithm>
#include <random>
#include <cmath>
#include <chrono>


using namespace std;

class stoper
{
    chrono::time_point<chrono::high_resolution_clock> a,b;
public:
    void start(){a = chrono::high_resolution_clock::now();}
    void stop() {b = chrono::high_resolution_clock::now();}
    double value()
    {
        chrono::duration<double> elapsed_seconds = b-a;
        return elapsed_seconds.count();
    }
};

template <typename F>
void test(F fun)
{
    stoper st;
    st.start();
    fun();
    st.stop();
    cout<<st.value()<< " s"<<endl;
}


template <typename iter_t >
typename iter_t::value_type stable_mul_1 ( const iter_t first, const iter_t last ){
    using T = typename iter_t::value_type;
    constexpr T one = T(1);

    auto big=first;
    auto small = first;
    T result =  one;

    while ( abs(*small) >= one && small!=last ) ++small;
    while ( abs(*big)   < one && big  !=last ) ++big;

    while ((big!= last) && (small !=last))    {

        while( abs(result) >= one && small!=last  ) {
            //cout<<"mul "<< *small<<endl;
            result *= *(small++);

            while ( abs(*small) >= one && small!=last ) ++small;
        }
        while( abs(result) < one && big!=last  ) {
          //  cout<<"mul "<< *big<<endl;
            result *= *(big++);

            while ( abs(*big)   < one && big  !=last ) ++big;
        }
    }

    //ogon

    while( small!=last  ) {
        //cout<<"mul "<< *small<<endl;
        result *= *(small++);
        while ( abs(*small) >= one && small!=last ) ++small;

    }

    while( big!=last  )
    {
//        cout<<"mul "<< *big<<endl;
        result *= *(big++);
        while ( abs(*big)   < one && big  !=last ) ++big;

    }
    return result;
}

template <typename iter_t >
typename iter_t::value_type stable_mul_partition (const iter_t first, const iter_t last ){
    using T = typename iter_t::value_type;
    constexpr T one = T(1);

    auto mid = partition(first, last, [](T x) {return abs(x)>one;} );

    auto small = mid;
    auto big = first;

    T result =  one;

    while ((big!= mid ) && (small !=last))  {
        if (abs( result) >= one)
            result*= *(small++);
        else
            result*= *(big++);
    }
    while (big!= mid )
        result*= *(big++);
    while (small != last )
        result*= *(small++);
   return result;
}

template <typename iter_t >
typename iter_t::value_type stable_mul_frexp ( iter_t first, const iter_t last ){
    using T = typename iter_t::value_type;
    constexpr T one = T(1);

    T result =  one;
    int_fast64_t ex=0;
    int tex;

    for (; first !=last;++first )    {
        result*= *first;
        result=frexp(result,&tex);
        ex+=tex;
    }
   return ldexp(result,ex);
}


int main()
{

    vector<double> tab;

    tab.resize(1000000);

    random_device rd;

    mt19937 gen(rd());


    lognormal_distribution<double> dist(0.0,50.0);//liczby z roskładu lognormalnego z ogromnego zasięgu
    generate(begin(tab),end(tab),[&]{return dist(gen);});
    for (auto it=tab.rbegin(); it!= prev(tab.rend()); ++it) {
        *it= *it/(*(next(it)));                         // właściwe liczby to ilorazy, gwarantuje to, że wynik mnozenia jest w zakresie
    }
    shuffle(begin(tab),end(tab),gen); // wymieszajmy kolejność, wtedy zwykłe mnożenie najprawdopodobniej wpadnie na +-inf


    cout<<"probka liczb"<<endl;
    for (int i=0;i<20;i+=2)    {
       cout<<tab[i]<<endl;
    }
    cout<<"koniec probki liczb"<<endl;


    cout<<"mnozenie standardowe "<<accumulate(begin(tab),end(tab),1.0,std::multiplies<double>() )<<endl;


   double logacu=0.0;

   for (auto it=begin(tab); it!= end(tab);++it)
       logacu+=log(*it);
   cout<<"dodawanie logarytmów "<<exp(logacu)<<endl;

    test( [&] {cout<<"podwojny przelot     "<< stable_mul_1(begin(tab),end(tab)) <<" ";});

    test( [&] {cout<<"frexp                "<< stable_mul_frexp(begin(tab),end(tab)) <<" ";});

    test( [&] {cout<<"partition            "<< stable_mul_partition(begin(tab),end(tab)) <<" ";});

    return 0;
}


/*

probka liczb
1.7583e+33
1.93469e+44
4.25265e-75
1.59875e+27
2.43612e-19
4.97215e+39
2.14228e-19
2.37518e+29
3.94031e-62
3.15255e+81
koniec probki liczb
mnozenie standardowe inf
dodawanie logarytmów 4.11519e-05
podwojny przelot     4.11519e-05 0.0829092 s
frexp                4.11519e-05 0.0293514 s
partition            4.11519e-05 0.080956 s
**************************
probka liczb
1.63243e+13
90.8683
7.364e-11
6.92679e+10
5.5909e+83
1.15735e-06
2.54212e+56
4.15332e-49
193442
1.75265e-41
koniec probki liczb
mnozenie standardowe inf
dodawanie logarytmów 107.279
podwojny przelot     107.279 0.0825575 s
frexp                107.279 0.0291725 s
partition            107.279 0.0868694 s

*****************************
probka liczb
7.86603e+19
431585
3.34872e-56
9.67111e+29
1.69078e-08
1.24735e-26
2.94383e+21
1.00652e+52
1.4806e+16
3.64962e-07
koniec probki liczb
mnozenie standardowe inf
dodawanie logarytmów 4.06236e-07
podwojny przelot     4.06236e-07 0.0828723 s
frexp                4.06236e-07 0.029212 s
partition            4.06236e-07 0.0832593 s



*/