1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
#pragma once


#include <iterator><--- Include file:  not found. Please note: Cppcheck does not need standard library headers to get proper results.

#include "snippet/iterations.hpp"
#include "internal/types.hpp"
#include "internal/dev_env.hpp"
#include "numeric/prime_enumerator.hpp"


// Thanks to: https://nyaannyaan.github.io/library/multiplicative-function/divisor-multiple-transform.hpp.html
namespace uni {

namespace divisor_transform {


using size_type = internal::size_t;


template<std::size_t OFFSET, std::ranges::sized_range R>
    requires (!requires { typename R::mapped_type; })
auto zeta(R&& range) noexcept(NO_EXCEPT) {
    const auto n = std::ranges::ssize(range) + OFFSET - 1;
    ITR(p, uni::prime_enumerator(n)) {
        for(size_type k=1; k*p <= n; ++k) range[k*p - OFFSET] += range[k - OFFSET];
    }
}

template<std::size_t OFFSET, std::ranges::sized_range R>
    requires (!requires { typename R::mapped_type; })
auto mobius(R&& range) noexcept(NO_EXCEPT) {
    const auto n = std::ranges::ssize(range) + OFFSET - 1;
    ITR(p, uni::prime_enumerator(n)) {
        for(size_type k=n/p; k>0; --k) range[k*p - OFFSET] -= range[k - OFFSET];
    }
}


template<std::ranges::sized_range R>
auto zeta(R&& range) noexcept(NO_EXCEPT) { return zeta<0>(range); }

template<std::ranges::sized_range R>
auto mobius(R&& range) noexcept(NO_EXCEPT) { return mobius<0>(range); }


template<std::ranges::input_range R>
    requires requires { typename R::mapped_type; }
auto zeta(R&& range) noexcept(NO_EXCEPT) {
    auto begin = std::ranges::begin(range);
    auto end = std::ranges::end(range);
    for(auto itr1 = end; itr1-- != begin; ) {
        for(auto itr2 = begin; itr2 != end; ++itr2) {
            if(itr1->first == itr2->first) break;
            if(itr1->first % itr2->first == 0) itr1->second += itr2->second;
        }
    }
}

template<std::ranges::input_range R>
    requires requires { typename R::mapped_type; }
auto mobius(R&& range) noexcept(NO_EXCEPT) {
    auto begin = std::ranges::begin(range);
    auto end = std::ranges::end(range);
    for(auto itr2 = begin; itr2 != end; ++itr2) {
        for(auto itr1 = end; (itr1--) != begin; ) {
            if(itr1->first == itr2->first) break;
            if(itr1->first % itr2->first == 0) itr1->second -= itr2->second;
        }
    }
}


} // namespace divisor_transform


namespace multiple_transform  {


using size_type = internal::size_t;


template<std::size_t OFFSET, std::ranges::sized_range R>
    requires (!requires { typename R::mapped_type; })
auto zeta(R&& range) noexcept(NO_EXCEPT) {
    const auto n = std::ranges::ssize(range) + OFFSET - 1;
    ITR(p, uni::prime_enumerator(n)) {
        for(size_type k=n/p; k>0; --k) range[k - OFFSET] += range[k*p - OFFSET];
    }
}

template<std::size_t OFFSET, std::ranges::sized_range R>
    requires (!requires { typename R::mapped_type; })
auto mobius(R&& range) noexcept(NO_EXCEPT) {
    const auto n = std::ranges::ssize(range) + OFFSET - 1;
    ITR(p, uni::prime_enumerator(n)) {
        for(size_type k=1; k*p <= n; ++k) range[k - OFFSET] -= range[k*p - OFFSET];
    }
}


template<std::ranges::sized_range R>
auto zeta(R&& range) noexcept(NO_EXCEPT) { return zeta<0>(range); }

template<std::ranges::sized_range R>
auto mobius(R&& range) noexcept(NO_EXCEPT) { return mobius<0>(range); }


template<std::ranges::input_range R>
    requires requires { typename R::mapped_type; }
auto zeta(R&& range) noexcept(NO_EXCEPT) {
    auto begin = std::ranges::begin(range);
    auto end = std::ranges::end(range);
    for(auto itr2 = begin; itr2 != end; ++itr2) {
        for(auto itr1 = end; --itr1 != itr2; ) {
            if(itr1->first % itr2->first == 0) itr2->second += itr1->second;
        }
    }
}

template<std::ranges::input_range R>
    requires requires { typename R::mapped_type; }
auto mobius(R&& range) noexcept(NO_EXCEPT) {
    auto begin = std::ranges::begin(range);
    auto end = std::ranges::end(range);
    for(auto itr2 = end; itr2-- != begin; ) {
        for(auto itr1 = end; --itr1 != itr2; ) {
            if(itr1->first % itr2->first == 0) itr2->second -= itr1->second;
        }
    }
}


} // namespace multiple_transform

} // namespace uni