计算机程序设计 C1 作业记

发表于 2022-09-23 分类于作业

来到妮可大学的第一次 CS 大作业。

作业内容

目标：编写一个程序，随机产生一个 blum 数，并对其进行素性测试。

blum 数是两个大素数 $p$ 、 $q$ 的乘积，其中 $p$ 和 $q$ 均是二进制下不小于五百位的、模 $4$ 余 $3$ 的大素数。

限定使用 C99，只允许使用 stdio、stdlib、time 库。

time 库用于随机数初始化、统计运行时间；stdlib 库用于产生随机数
其他功能函数均要求自己实现

要求：

使用 Miller Rabin 算法来测试素数，在判定素数时设定随机选取 $10$ 个数测试通过。
产生的 blum 数以 $16$ 进制形式在屏幕显示。
使用 $20$ 个随机数对产生的 blum 数进行 Miller Rabin 素性测试，给出 $20$ 次的测试结果。
两分钟出一组结果。

过程

$9$ 月 $6$ 号，也就是第一次上课当天，老师就布置了这个作业，不过是 $10$ 月底交。上网了解了一下，发现其实我们要做的就是生成 RSA 密钥，于是又去学习了一下著名的蒙哥马利算法。

然后开始写，高精度还调了一个小时，花了不少时间。

然后 Miller Rabin 的时候快速幂直接暴力模乘，又调了半个小时，可以正确输出结果了。一看，全过程要花 $90\sim 100\,\text{s}$ 的时间，那看来 $2$ 分钟这个要求还是很低的。

然后就开始写蒙哥马利，写到要睡觉了还是没有写完。

第二天中午继续写，终于写完了，一运行，啥都输出不出来，寄。

晚上接着调，调了一整个晚上可以正确输出了，时间降为 $30\sim 40\,\text{s}$ ，优化还是挺大的。

到了 $9$ 月 $8$ 号，看到群里 zzeh 大佬的程序已经起飞了，于是效仿他将位权改成了 $2^{28}$ ，直接变成 $4\sim 5\,\text{s}$ 出结果。我之前一直认为 long long 会慢，于是没去压这么多位。

之后一直没什么进展，而群里已经有大佬卡到一秒之内了。在 $9$ 月 $16$ 号，我再一次效仿 zzeh 大佬，去学习了一下 Barrett Reduction 算法，发现挺好理解的，代码也很好写（本人写了一篇学习笔记）。试了一下，把蒙哥马利去掉，加上 Barrett Reduction，花了二十分钟时间，一运行，我去，只要用 $0.4\sim 1\,\text{s}$ 了~~（垃圾蒙哥马利毁我青春）~~。

于是就不卡了，都已经这么快了。到了 $9$ 月 $20$ 号第一次上机课，我就把程序拿给助教检查了~~（我才不会告诉你检查我程序的助教是学姐呢）~~，助教想检查我生成的 $p$ 和 $q$ 是不是素数，但是她用的那个网站只支持输入十进制数，而我用了 Barrett Reduction，改成十进制输出比较困难。于是助教就叫我讲一下 Miller Rabin 的算法流程，我讲了一遍，她就说 ok 了，思路正确代码看起来也没有什么问题，给了我满分。

之后找助教要了一下素性测试的网址~~（加到学姐 qq 啦！！）~~，测了一下，我程序输出的 $p,q$ 确实是质数。应该就是没有什么问题了。

但后来另一个助教说要自己写一个 Miller Rabin 来测我们的，之前的成绩无效，那看来还要再等一段时间才能知道成绩。

Updata（ $2022.10.22$ ）：助教用自己写得 Miller Rabin 重新测了一边，没有问题，过了。

~~吐槽一下，C 语言真的难用，可能是 C++ 用太多了。~~

还没完

$10$ 月 $27$ 号，学姐发 qq 过来正式通知我说我过了。但是她又补了一句，说我的代码在 Linux 与 Mac 系统下都要跑十分钟，但她最后记录的我在我自己 Windows 电脑上跑的 $0.8\,\text{s}$ 的结果。

我就感觉很疑惑，即使是系统差异也不应该差这么多啊。学姐说可能是不同系统下随机数生成有差异，这突然提醒了我。

我看到了我代码的第 $135$ 行（是我下面展示的代码唯一加了注释的一行），原来我写的是：

1	p=rand()/4*rand()/32

我突然想起来只有 Windows 下系统 rand() 生成数的范围是 $0\sim 32767$ ，其他系统下都是 $0\sim 2^{31}-1$ ，而我 $p$ 又是用 int 存的，那不就溢出成负数了！！？？

于是叫学姐帮忙改成了：

p=rand()

但是还是寄，还是要跑 $10$ 分钟。

我突然又想到，在第 $111$ 行（即 Pow 函数的第一行），我直接将 p 赋值给了 a.s[1]，也就是 a 这个高精度数的第一位。但是我设置的位权只有 $2^28$ ，可 rand() 生成的数的最大值达到 $2^{31}-1$ ，那又不对了。

于是我又叫学姐改成：

1	p=(ll)rand()*rand()%B

这下就正常了，平均 $0.7\,\text{s}$ 出一个结果。

~~垃圾 Windows 毁我青春。~~

只能说自己太蠢了，系统差异问题都没有考虑到。

Code

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define ll long long
#define lf double
const ll B=1<<28;

typedef struct{ll s[124],l;}Int;
Int p,q,r,s,t,y,z; int k,c; int n,pr[10005],o[100010];

int Max(int a,int b){return a>b?a:b;}

void Clear(ll *a,int l){int i; for (i=0; i<=l; i++) a[i]=0;}

void Copy(ll *a,ll *b,int l){int i; for (i=0; i<l; i++) a[i]=b[i];}

void Print(Int a){int i; for (printf("0x%llx",a.s[a.l]),i=a.l-1; i; i--) printf("%07llx",a.s[i]);}

int Cmp(Int a,Int b)
{
    if (a.l>b.l) return 1; if (a.l<b.l) return -1; int i;
    for (i=a.l; i; i--)
        if (a.s[i]>b.s[i]) return 1;
        else if (a.s[i]<b.s[i]) return -1;
    return 0;
}

Int Lsh(Int x,int k)
{
    if (k>=x.l) return Clear(x.s,x.l),x.l=1,x;
    return Copy(x.s+1,x.s+k+1,x.l-k),x.l-=k,Clear(x.s+x.l+1,k),x;
}

Int Add(Int a,Int b)
{
    a.l=Max(a.l,b.l); int i;
    for (i=1; i<=a.l; i++) a.s[i]+=b.s[i],a.s[i+1]+=a.s[i]/B,a.s[i]%=B;
    return a.l+=a.s[a.l+1]>0,a;
}

Int add(Int a,int b)
{
    a.s[1]+=b; int i;
    for (i=1; i<=a.l; i++) a.s[i+1]+=a.s[i]/B,a.s[i]%=B;
    return a.l+=a.s[a.l+1]>0,a;
}

Int Sub(Int a,Int b)
{
    a.l=Max(a.l,b.l); int i;
    for (i=1; i<=a.l; i++) if (a.s[i]-=b.s[i],a.s[i]<0) a.s[i]+=B,a.s[i+1]--;
    for ( ; !a.s[a.l]&&a.l>1; a.l--); return a;
}

Int Mul(Int a,Int b)
{
    Int E; Clear(E.s,123),E.l=a.l+b.l-1; int i,j;
    for (i=1; i<=a.l; i++) for (j=1; j<=b.l; j++)
        E.s[i+j-1]+=a.s[i]*b.s[j],E.s[i+j]+=E.s[i+j-1]/B,E.s[i+j-1]%=B;
    return E.l+=E.s[E.l+1]>0,E;
}

Int mul(Int a,int b)
{
    int i,x; for (i=1; i<=a.l; i++) a.s[i]*=b;
    for (i=1; i<=a.l; i++) a.s[i+1]+=a.s[i]/B,a.s[i]%=B;
    for (x=a.s[a.l+1]; x; x/=B) a.s[++a.l]=x%B;
    return a;
}

Int Div(Int a,Int b)
{
    Int E,c; int i,j,l,r,I;
    Clear(E.s,123),E.l=0,Clear(c.s,123),c.l=1;
    for (i=a.l; i; i--)
    {
        c=add(mul(c,B),a.s[i]);
        for (l=j=0,r=B-1; l<=r; )
        {
            if (I=l+r>>1,Cmp(mul(b,I+1),c)<1) l=I+1;
            else j=I,r=I-1;
        }
        if (E.s[++E.l]=j) c=Sub(c,mul(b,j));
    }
    for (i=1; i+i<=E.l; i++) j=E.s[i],E.s[i]=E.s[E.l-i+1],E.s[E.l-i+1]=j;
    for ( ; !E.s[E.l]&&E.l>1; E.l--); return E;
}

Int div2(Int a)
{
    int i; for (i=a.l; i; i--) a.s[i-1]+=a.s[i]%2*B,a.s[i]/=2;
    for ( ; !a.s[a.l]&&a.l>1; a.l--); return a;
}

Int Mod(Int a)
{
    Int E=Sub(a,Mul(Lsh(Mul(a,s),k),t));
    while (Cmp(E,t)>=0) E=Sub(E,t);
    return E;
}

int mod(Int a,int b)
{
    int i,c=0;
    for (i=a.l; i; i--) c=(c*B+a.s[i])%b;
    return c;
}

Int Pow(int x,Int b)
{
    Int a,E; int i,j; Clear(a.s,123),a.l=1,a.s[1]=x,Clear(E.s,123),E.l=E.s[1]=1;
    for (i=1,j=0; i<b.l||b.s[i]>=(1<<j); j=(j+1)%28,i+=!j)
    {
        if (b.s[i]>>j&1) E=Mod(Mul(E,a));
        a=Mod(Mul(a,a));
    }
    return E;
}

int MR(int p)
{
    Int a=Pow(p,y),b=a; int i; c++;
    for (i=0; i<2&&(a.l>1||a.s[1]>1); i++,a=b)
        if (b=Mod(Mul(b,b)),b.l<2&&b.l==1&&Cmp(a,z)) return 0;
	return a.l<2&&a.s[1]==1;
}

void BR(Int x){t=x,k=2*x.l,Clear(s.s,123),s.s[s.l=k+1]=1,s=Div(s,x);}

int Isp(Int x,int t)
{
    BR(y=x),y.s[1]--,z=y,y=div2(y); int i,p;
    for (i=1; i<=t; i++)
    {
        if (p=(ll)rand()*rand()%B/*rand()/4*rand()/32 rand()*/,p&1^1) p++;
        if (!MR(p)) return printf("Fail at test %d.\n",i),0;
    }
    return puts("Success!"),1;
}

Int Ranp()
{
    lf s=clock(); Int E; int i; c=0;
    while (1)
    {
        puts("Generating prime..."); lf t=clock();
        Clear(E.s,123),E.l=rand()%2+18;
        for (i=2; i<E.l; i++) E.s[i]=(ll)rand()*rand()%B;
        E.s[1]=(ll)rand()*rand()%(B/4)*4+3,E.s[E.l]=(ll)rand()*rand()%(B-1)+1;
        for (i=1; i<=n; i++) if (!mod(E,pr[i])){puts("Fail at pretest."); goto NG;}
        if (!Isp(E,10)) goto NG;
        printf("Total time cost %dms.\nAvg time cost %dms.\n",(int)(clock()-s+0.5),(int)((clock()-s)/c+0.5));
        return E; NG: printf("Time cost %dms.\n",(int)(clock()-t+0.5));
    }
}

int main()
{
    lf s=clock(),t; freopen("c1.txt","w",stdout);

    srand(time(0)); int i,j,k;
    for (i=2; i<100005; i++)
    {
        if (!o[i]) pr[++n]=i;
        for (j=1; j<=n&&(k=i*pr[j])<100005; j++) if (o[k]=1,i%pr[j]<1) break;
    }

    puts("prime p:"),p=Ranp(),puts("");
    puts("prime q:"),q=Ranp(),puts("");

    puts("blum:\nMultiplying..."),r=Mul(p,q);
    c=0,t=clock(),Isp(r,20);
    printf("Time cost %dms.\nAvg time cost %dms.\n\n",(int)(clock()-t+0.5),(int)((clock()-t)/c+0.5));

    printf("Total time cost %dms\n\n",(int)(clock()-s+0.5));
    printf("p: "),Print(p),puts("");
    printf("q: "),Print(q),puts("");
    printf("blum: "),Print(r),puts("");

    return 0;
}