Divisor Multiple Transform
(Convolution/divisor.hpp)

• View this file on GitHub
• Last update: 2023-01-16 20:41:46+09:00
• Include: #include "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

• Prime Sieve (Math/sieve.hpp)

Verified with

• Verify/LC_gcd_convolution.test.cpp
• Verify/LC_lcm_convolution.test.cpp

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=50101010>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