library

This documentation is automatically generated by online-judge-tools/verification-helper


Project maintained by tko919 Hosted on GitHub Pages — Theme by mattgraham

:heavy_check_mark: Divisor Multiple Transform
(Convolution/divisor.hpp)

使い方

void DivisorTransform::zeta(vector<T>& a): $a’[n]=\sum_{n \bmod d=0} a[d]$ を計算。 mobius(vector<T>& a) は逆変換。
void MultipleTransform::zeta(vector<T>& a): $a’[n]=\sum_{k \bmod n=0} a[k]$ を計算。 mobius(vector<T>& a) は逆変換。

Depends on

Required by

Verified with

Code

#pragma once
#include "Math/sieve.hpp"


namespace DivisorTransform{
    int n;
    vector<int> ps;
    template<typename T>void zeta(vector<T>& a){
        int N=a.size()-1;
        if(n<N){
            ps=sieve(N);
            n=N;
        }
        for(auto& p:ps){
            for(int k=1;k*p<=N;k++)a[k*p]+=a[k];
        }
    }
    template<typename T>void mobius(vector<T>& a){
        int N=a.size()-1;
        if(n<N){
            ps=sieve(N);
            n=N;
        }
        for(auto& p:ps){
            for(int k=N/p;k;k--)a[k*p]-=a[k];
        }
    }
};

namespace MultipleTransform{
    int n;
    vector<int> ps;
    template<typename T>void zeta(vector<T>& a){
        int N=a.size()-1;
        if(n<N){
            ps=sieve(N);
            n=N;
        }
        for(auto& p:ps){
            for(int k=N/p;k;k--)a[k]+=a[k*p];
        }
    }
    template<typename T>void mobius(vector<T>& a){
        int N=a.size()-1;
        if(n<N){
            ps=sieve(N);
            n=N;
        }
        for(auto& p:ps){
            for(int k=1;k*p<=N;k++)a[k]-=a[k*p];
        }
    }
};

/**
 * @brief Divisor Multiple Transform
 * @docs docs/divisor.md
 */
#line 2 "Math/sieve.hpp"

template <int L = 101010101> vector<int> sieve(int N) {
    bitset<L> isp;
    int n, sq = ceil(sqrt(N));
    for (int z = 1; z <= 5; z += 4) {
        for (int y = z; y <= sq; y += 6) {
            for (int x = 1; x <= sq and (n = 4 * x * x + y * y) <= N; ++x) {
                isp[n].flip();
            }
            for (int x = y + 1; x <= sq and (n = 3 * x * x - y * y) <= N;
                 x += 2) {
                isp[n].flip();
            }
        }
    }
    for (int z = 2; z <= 4; z += 2) {
        for (int y = z; y <= sq; y += 6) {
            for (int x = 1; x <= sq and (n = 3 * x * x + y * y) <= N; x += 2) {
                isp[n].flip();
            }
            for (int x = y + 1; x <= sq and (n = 3 * x * x - y * y) <= N;
                 x += 2) {
                isp[n].flip();
            }
        }
    }
    for (int y = 3; y <= sq; y += 6) {
        for (int z = 1; z <= 2; ++z) {
            for (int x = z; x <= sq and (n = 4 * x * x + y * y) <= N; x += 3) {
                isp[n].flip();
            }
        }
    }
    for (int n = 5; n <= sq; ++n)
        if (isp[n]) {
            for (int k = n * n; k <= N; k += n * n) {
                isp[k] = false;
            }
        }
    isp[2] = isp[3] = true;

    vector<int> ret;
    for (int i = 2; i <= N; i++)
        if (isp[i]) {
            ret.push_back(i);
        }
    return ret;
}

/**
 * @brief Prime Sieve
 */
#line 3 "Convolution/divisor.hpp"

namespace DivisorTransform{
    int n;
    vector<int> ps;
    template<typename T>void zeta(vector<T>& a){
        int N=a.size()-1;
        if(n<N){
            ps=sieve(N);
            n=N;
        }
        for(auto& p:ps){
            for(int k=1;k*p<=N;k++)a[k*p]+=a[k];
        }
    }
    template<typename T>void mobius(vector<T>& a){
        int N=a.size()-1;
        if(n<N){
            ps=sieve(N);
            n=N;
        }
        for(auto& p:ps){
            for(int k=N/p;k;k--)a[k*p]-=a[k];
        }
    }
};

namespace MultipleTransform{
    int n;
    vector<int> ps;
    template<typename T>void zeta(vector<T>& a){
        int N=a.size()-1;
        if(n<N){
            ps=sieve(N);
            n=N;
        }
        for(auto& p:ps){
            for(int k=N/p;k;k--)a[k]+=a[k*p];
        }
    }
    template<typename T>void mobius(vector<T>& a){
        int N=a.size()-1;
        if(n<N){
            ps=sieve(N);
            n=N;
        }
        for(auto& p:ps){
            for(int k=1;k*p<=N;k++)a[k]-=a[k*p];
        }
    }
};

/**
 * @brief Divisor Multiple Transform
 * @docs docs/divisor.md
 */
Back to top page