反素数

发表于 2023-06-27 | 分类于算法学习

反素数

模板

反素数定义：对于正整数$n$，任何比$n$小的数正约数个数都小于$n$的正约数个数，则称$n$是反素数。

反素数性质：

性质1：若正整数$n$是反素数，则$n$的质因数是若干个最小的质数的乘积，且它们的指数单调不增。

性质2：不大于$n$的最大反素数是$1,2,\cdots n$中约数最多的数中最小的。

基于反素数的以上性质，我们可以采用dfs搜索查找反素数，分别确定每个质因数的指数，更新最大的约数个数。

long long p[16] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53};//列举前几个素数
long long n,ans,ans_num;
void dfs(int d,long long tmp,long long num,int mi)
{
    if(d >= 16 || tmp > n)
    return ;//超出边界则返回
    if(num > ans_num || (num == ans_num && tmp < ans))
    {
        ans_num = num;
        ans = tmp;
    }//更新约数个数
    for(int i = 1;i <= mi;i++)//枚举第d个质数的指数
    {
        if(tmp * p[d] > n)
        break;
        tmp = tmp * p[d];
        dfs(d + 1,tmp,num * (i + 1),i);
    }
}

例题

CF27E A Number With The Given Amount Of Divisors

题目链接：https://codeforces.com/problemset/problem/27/E

题意

给定一个正整数$n$，求约数个数恰好是$n$的最小正整数。

数据范围

$1\leq n \leq 1000$，所求正整数不大于$10^{18}$。

题解

易知这样的$n$满足反素数的性质，在dfs的时候将约数等于$n$作为更新的条件即可。

代码

#include<iostream>
#define maxn 1e18
using namespace std;
long long p[16] = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53};
long long ans,n;
void dfs(int d,long long tmp,long long num,int mi)
{
    if(d > 16 || tmp > maxn)
    return ;
    if(n == num && tmp < ans)
    ans = tmp;
    for(int i = 1;i <= mi;i++)
    {
        if(p[d] * tmp > maxn)
        break;
        tmp = tmp * p[d];
        dfs(d + 1,tmp,num * (i + 1),i);
    }
}
int main()
{
    cin >> n;
    ans = maxn;
    dfs(0,1,1,63);
    cout << ans;
    system("pause");
    return 0;
}

线性筛

发表于 2023-06-24 | 分类于算法学习

线性筛

模板

线性筛是在埃氏筛的基础上保证每个质数只会被它的最小质因数筛一次的筛法，可以将时间复杂度优化到$O(n)$。

对于数$i$，从小到大枚举$prime$数组（已经筛出的质数数组），当第一次出现$prime[j] \mid i$ 时，可知$prime[j]$是$i\times prime[j]$的最小质因数，而对于$i \times prime[j+1]$，它的最小质因数是$prime[j]$，此时为了保证只被最小质因数筛掉，停止枚举。

对于合数$x$，设它的最大质因数为$x_0$，在枚举到$\frac{x}{x_0}$时一定会被筛掉，则所有的合数都能够被筛掉。

const int N = 100010;
int prime[N],n,now;
bool isprime[N];
void get_prime()
{
    for(int i = 2;i <= n;i++)
    {
        if(!isprime[i])
        now++,prime[now] = i;//如果没有被筛掉，就是质数
        for(int j = 1;j <= now && i * prime[j] <= n;j++)
        {
            isprime[i * prime[j]] = 1;
            if(i % prime[j] == 0)
            break;//保证每个合数只会被它的最小质因数筛掉
        }
    }
}

例题

acwing196 质数距离

题目链接：https://www.acwing.com/problem/content/198/

题意

给定两个正整数$L,U$，在闭区间$[L,U]$中找到距离最近和最远的相邻质数。

数据范围

$1 \leq L<U\leq 2^{31}-1,U-L\leq 10^{6}$

题解

可以考虑求出$L$到$U$区间内的全部质数，但$L$和$U$本身数据范围较大，线性筛无法做到。

考虑使用二次筛法，对于合数$X$，由于$X$必然存在一个不大于$\sqrt{X}$的质数因子，于是先用线性筛筛出所有小于$10^{5}$的质数，分别对这些质数标记它们的倍数中在$[L,U]$的，剩下没有被标记过的就是所有的质数。

最后由于线性查找其中相邻两个质数中差最大和最小的即可，时间复杂度$O(U-L)$

代码

#include<iostream>
#include<vector>
#include<cstring>
using namespace std;
const int N = 100010;
long long prime[N],l,r,now;
bool isprime[N],isans[N * 10];
vector<long long> v;
void get_prime()
{
    for(int i = 2;i <= N;i++)
    {
        if(!isprime[i])
        now++,prime[now] = i;
        for(int j = 1;j <= now && i * prime[j] <= N;j++)
        {
            isprime[i * prime[j]] = 1;
            if(i % prime[j] == 0)
            break;
        }
    }
}
void get_ans()
{
    memset(isans,0,sizeof(isans));
    for(int i = 1;i <= now;i++)
    {
        long long ll = l / prime[i] * prime[i];
        while(ll <= r)
        {
            if(ll >= l && ll != prime[i])
            isans[ll - l] = 1;
            ll += prime[i];
        }
    }
}
int main()
{
    get_prime();
    while(cin >> l >> r)
    {
        if(l == 1)
        l = 2;
        get_ans();
        v.clear();
        for(int i = 0;i <= r - l;i++)
        {
            if(!isans[i])
            v.push_back(i + l);
        }
        if(v.size() <= 1)
        printf("There are no adjacent primes.\n");
        else
        {
            long long minn = 1e6,maxn = 0,l1,r1,l2,r2;
            for(int i = 1;i < v.size();i++)
            {
                if(v[i] - v[i - 1] < minn)
                l1 = v[i - 1],r1 = v[i],minn = v[i] - v[i - 1];
                if(v[i] - v[i - 1] > maxn)
                l2 = v[i - 1],r2 = v[i],maxn = v[i] - v[i - 1];
            }
            printf("%lld,%lld are closest, %lld,%lld are most distant.\n",l1,r1,l2,r2);
        }
    }
    system("pause");
    return 0;
}

FFT

发表于 2023-03-22 | 分类于算法学习

前置知识

多项式的表示：一般来说，多项式的表示方法有点表示和系数表示两种。

系数表示法：一个$n$次多项式可以由$n+1$个系数唯一确定，也即$f(x)=\sum_{i=0}^{n}a_ix^{i}$。

点表示法：一个$n$次多项式可以由$n+1$个各不相同的点唯一确定，设这$n+1$个点的坐标为$(x_0,f(x_0)),(x_1,f(x_1)),\cdots,(x_n,f(x_n))$，代入表达式可得方程组： \[ \begin{cases} a_0+a_1x_0+\cdots+a_nx_0^{n}=f(x_0)\\ a_0+a_1x_1+\cdots+a_nx_1^{n}=f(x_1)\\ \cdots\\ a_0+a_1x_n+\cdots+a_nx_n^{n}=f(x_n)\end{cases}\] 其系数行列式的值由范德蒙行列式可求得： \[\begin{vmatrix}1 & x_0 & \dots &x_0^{n}\\ 1 & x_1 & \dots & x_1^{n} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & \dots & x_n^{n} \end{vmatrix} = \prod_{0\leq i\leq j \leq n}(x_i-x_j)\neq 0\] 可知多项式的系数被唯一确定。

复数的表示与运算：一个复数可以表示为$z=a+bi(a,b\in \mathbb{R})$，对应复平面上一点$P(a,b)$。其模长$|z|=\sqrt {a^{2}+b^{2}}$,当$z\neq0$时，定义向量$\overrightarrow{OP}$与$x$轴正向夹角为辐角$\theta$，则复数也可表示为$z=|z|(\cos\theta+i\sin\theta)$。

复数的加（减）法：设$z_1=a_1+ib_1$，$z_2=a_2+ib_2$，则$z_1\pm z_2=(a_1\pm a_2)+i(b_1 \pm b_2)$。

复数的乘法：设$z_1=a_1+ib_1$，$z_2=a_2+ib_2$，则$z_1z_2=(a_1a_2-b_1b_2)+i(a_1b_2+a_2b_1)$。复数的乘法也可以理解为模长相乘，辐角相加。

单位根：考虑$x^{n}=1$在复数域上的$n$个解，在单位圆上将圆从$(1,0)$开始等分成n份，这$n$个点均满足上述方程，将它们对应的复数记为$\omega_n^{i}(1 \leq i \leq n-1)$，且有$\omega_n^{i}=\cos{\frac{2\pi i}{n}}+i\sin{\frac{2\pi i}{n}}$。

单位根的性质：1）$\omega_{2n}^{2k}=\omega_n^{k}$；

2）$\omega_n^{k+\frac{n}{2}}=-\omega_n^{k}$；

3）$\omega_n^{0}=\omega_n^{n}=1$；

FFT（快速傅里叶变换）

FFT是用来解决两个多项式相乘的一种算法，FFT大致思路如下：先将两个多项式写成点表示法，则所求多项式的点值等于两个多项式点值的乘积，再将其重新转化为系数表示法。

将系数表示法转化为点表示法

取$n$为2的次幂，设$A(x)=a_0+a_1x+\cdots+a_{n-1}x^{n-1}$，将下标按奇偶性分开，则

$A(x)=(a_0+a_2x^{2}+\cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^{3}+\dots+a_{n-1}x^{n-1})$

$=(a_0+a_2x^{2}+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^{2}+\cdots+a_{n-1}x^{n-2})$

令$A_1(x)=a_0+a_2x+\cdots+a_{n-2}x^{\frac{n}{2}-1}$，$A_2(x)=a_1+a_3x+\cdots+a_{n-1}x^{\frac{n}{2}-1}$

则$A(x)=A_1(x^{2})+xA_2(x^{2})$，将$\omega_n^{k}(k<\frac{n}{2})$代入得：

$A(\omega_n^{k})=A_1(\omega_n^{2k})+xA_2(\omega_n^{2k})=A_1(\omega_\frac{n}{2}^{k})+\omega_n^{k}A_2(\omega_\frac{n}{2}^{k})$

$A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+xA_2(\omega_n^{2k+n})=A_1(\omega_\frac{n}{2}^{k})-\omega_n^{k}A_2(\omega_\frac{n}{2}^{k})$

可以发现两部分只有后半部分的符号不一样，则若已知$A_1(\omega_\frac{n}{2}^{k})$与$A_2(\omega_\frac{n}{2}^{k})$可以求得$A(\omega_n^{k})$和$A(\omega_n^{k+\frac{n}{2}})$，从而我们可以使用递归分治的方式求得$A(\omega_n^{k})(0\leq k\leq n-1)$，时间复杂度$O(nlogn)$。

将点表示法重新转化为系数表示法

已知$n+1$个点的坐标$(x_0,A(x_0)),(x_1,A(x_1)),\cdots,(x_{n-1},A(x_{n-1}))$，其中$x_i=\omega_n^{i}(0\leq i\leq n-1)$，先证明：\[ a_k=\frac{1}{n}\sum_{i=0}^{n-1}A(\omega_n^{i})(\omega_n^{-k})^i\]证明：

\[\sum_{i=0}^{n-1}A(\omega_n^{i})(\omega_n^{-k})^i\] \[=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i\] \[=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^{j-k})^i)\] \[=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i)\tag{1}\]

记$S(x)=1+x+x^2+\cdots+x^{n-1}$，考虑$S(\omega_n^{k})$：

1）若$k=0$，$S(\omega_n^{k})=n$；

2）若$k\neq 0$，$S(\omega_n^{k})=\omega_n^0+\omega_n^k+\cdots+\omega_n^{(n-1)k}$，$\omega_n^{k}S(\omega_n^{k})=\omega_n^k+\omega_n^{2k}+\cdots+\omega_n^{nk}$，两式相减可得：$(1-\omega_n^k)S(\omega_n^k)=0$，则$S(\omega_n^{k})=0$。

从而$(1)=na_k$，问题得证。

令$B(x)=\sum_{i=0}^{n-1}A(\omega_n^{i})x^i$，则$a_k=\frac{1}{n}B(\omega_n^{-k})$，问题又转化为了从系数表示法转化为点表示法。

迭代优化

注意到最终的位置为其初始位置二进制翻转后得到的，则可以将最终的位置处理出来后，向上迭代处理就可以得到最初的序列。

代码

#include<iostream>
#include<cmath>
using namespace std;
const double pi = acos(-1);
struct Complex
{
    double x,y;
    Complex operator+ (const Complex& t)const
    {
        return {x + t.x,y + t.y};
    }
    Complex operator- (const Complex& t)const
    {
        return {x - t.x,y - t.y};
    }
    Complex operator* (const Complex& t)const
    {
        return {x * t.x - y * t.y,x * t.y + y * t.x};
    }
}a[300010],b[300010];
int n,m,bit,tot,pos[300010];
void fft(Complex a[],int mode)
{
    for(int i = 0;i < tot;i++)
    {
        if(i < pos[i])
        swap(a[i],a[pos[i]]);
    }
    for(int mid = 1;mid < tot;mid = mid * 2)
    {
        auto wi = Complex({cos(pi / mid),mode * sin(pi / mid)});
        for(int i = 0;i < tot;i = i + mid * 2)
        {
            auto w0 = Complex({1,0});
            for(int j = 0;j < mid;j++)
            {
                auto tx = a[i + j],ty = w0 * a[i + j + mid];
                a[i + j] = tx + ty,a[i + j + mid] = tx - ty;
                w0 = w0 * wi;
            }
        }
    }
}
int main()
{
    scanf("%d%d",&n,&m);
    for(int i = 0;i <= n;i++)
    scanf("%lf",&a[i].x);
    for(int i = 0;i <= m;i++)
    scanf("%lf",&b[i].x);
    while((1 << bit) < n + m + 1)
    bit++;
    tot = 1 << bit;
    for(int i = 0;i < tot;i++)
    pos[i] = (pos[i >> 1] >> 1) | ((i & 1) << (bit - 1));
    fft(a,1),fft(b,1);
    for(int i = 0;i < tot;i++)
    a[i] = a[i] * b[i];
    fft(a,-1);
    for(int i = 0;i <= n + m;i++)
    printf("%d ",(int)(a[i].x / tot + 0.5));
    system("pause");
    return 0;
}

Hello World

发表于 2023-03-19

Welcome to Hexo! This is your very first post. Check documentation for more info. If you get any problems when using Hexo, you can find the answer in troubleshooting or you can ask me on GitHub.

Quick Start

Create a new post

1	$ hexo new "My New Post"

More info: Writing

Run server

1	$ hexo server

More info: Server

Generate static files

1	$ hexo generate

More info: Generating

Deploy to remote sites

1	$ hexo deploy

More info: Deployment