分类目录归档：研究

走别人没有走过的路。

AKS素性检测算法的C语言实现

按道理说这个算法提出之后，应该是有很多语言的实现版本。但是对于多项式的方幂运算，C语言版本总感觉缺少点什么。以前在其他地方见过，但多项式方幂的实现算法复杂度较高。因此自己也实现一个版本。但目前遇到一个问题，就是64bit的两个整数相乘，会溢出。虽然是在模n（也是64bit的）下运算的，但为了避免溢出，好好的乘法非要改用加法实现。希望在这方面可以有改进的地方( Mul64Mod )。

备注：2022年11月1日修订。根据文献[3]的提示，GNU G++17可以使用 unsigned __int128 数据类型，虽然我的编译器（ g++.exe (MinGW-W64 x86_64-ucrt-posix-seh, built by Brecht Sanders) 12.2.0 ）提示使用的数据类型应该为 __uint128_t 。但起码目前模乘法的速度问题算是基本解决了。另外，我将当前版本保存了一份在网上：https://onlinegdb.com/dFtGoQgL6 或者 https://gist.github.com/tpu01yzx/172c65d3003bb3a09a941e69bc6c370b

#include <stdio.h>
#include <stdint.h>
#include <inttypes.h>
#include <memory.h>
#include <math.h>

#define N 100
#define MAX_FACTORS 32
#define MAXR 320

typedef unsigned long long int uint64;
typedef unsigned int uint32;
typedef unsigned char uint8;

uint64 gcd(uint64 a, uint64 b) {
    uint64 c = a;
    while(b) {
        c = a % b;
        a = b;
        b = c;
    }
    return a;
}

uint8 BitCount(uint64 n) {
    uint8 ans = 0;
    while(n) {
        n>>=1;
        ans++;
    }
    return ans;
}

uint32 SquareRoot(uint64 x) {
    if(x < (1ULL<<32)) {
        return (uint32)sqrt((uint32)x);
    }
    uint8 logx = BitCount(x);    
    uint32 l = pow(2.0, (double)(logx-1) / 2);
    uint32 r = pow(2.0, (double)(logx) / 2);
    uint32 m = l;
    uint64 m2 = x;

    while(l <= r) {        
        m = (l + r ) / 2;        
        m2 = m * m;
        
        if(m2 == x) return m;
        if(m2 < x) {
            l = m + 1;
        } else {
            r = m - 1;
        }
    }
    return m;
}

uint64 Power(uint32 a, uint8 k) {
    if(k == 0) return 1;

    uint64 ans = 1;
    uint64 a2 = a;
    while(k > 1) {
        if(k & 0x01) {
            ans *= a2;
        }
        a2 = a2 * a2;
        k>>=1;
    }
    ans *= a2;
    return ans;
}

uint32 PowerMod(uint32 a, uint8 k, uint32 mod) {
    if(k == 0) return 1;

    uint64 ans = 1;
    uint64 a2 = a % mod;    
    while(k > 1) {
        if(k & 0x01) {
            ans = (ans * a2) % mod;            
        }
        a2 = (a2 * a2) % mod;        
        k>>=1;
    }
    ans = (ans * a2) % mod; 
    return (uint32)ans;
}

uint64 Mul64Mod(uint64 a, uint64 b, uint64 mod) {
    return (__uint128_t)a * b % mod;
}

void MulPoly(uint64 *p1, uint64 *p2, uint64 n, uint32 r) {
    uint64 ans[MAXR];
    int i, j;

    memset(ans, 0, sizeof(uint64) * r);
    for(i = 0; i < r; i++) {
        for(j = 0; j <= i; j++) {
            ans[i] += Mul64Mod(p1[j], p2[i - j], n);
            if(ans[i] >= n) ans[i] -= n;
        }
        for(j = i + 1; j < r; j++) {
            ans[i] += Mul64Mod(p1[j], p2[r + i - j], n);
            if(ans[i] >= n) ans[i] -= n;
        }        
    }    
    memcpy(p1, ans, sizeof(uint64) * r);
}

void PowerPoly(uint64 *coff, uint64 n, uint32 r) {
    if(n == 0) {
        memset(coff, 0, sizeof(uint64) * r);
        coff[0] = 1ULL;
        return;
    }
    uint64 n0 = n;
    uint64 ans[MAXR];
    uint64 coff2[MAXR];
    memset(ans, 0, sizeof(uint64) * r);    ans[0] = 1ULL; 
    memcpy(coff2, coff, sizeof(uint64) * r);

    while(n > 1) {
        if(n & 0x01) {
            MulPoly(ans, coff2, n0, r);
        }
        MulPoly(coff2, coff2, n0, r);
        n>>=1;
    }
    MulPoly(ans, coff2, n0, r);
    
    memcpy(coff, ans, sizeof(uint64) * r);
}

int CheckPoly(uint64 *p, uint64 a, uint64 n, uint32 r) {
    uint64 n0 = n % r;
    uint32 i;
    if(p[0] != a) return 0;    
    for(i = 1; i < n0; i++) {
        if(p[i] != 0) return 0;
    }
    if(p[n0] != 1ULL) return 0;
    for(i = n0 + 1; i < r; i++) {
        if(p[i] != 0) return 0;
    }
    return 1;
}

uint32 PerfectRoot(uint8 a, uint64 n, uint8 logn) {    
    if(a == 1) return n;
    uint32 l = pow(2.0, (double)(logn-1) / a);;
    uint32 r = pow(2.0, ((double)(logn) / a));
    uint32 m;
    uint64 mp;
    while(l <= r) {
        m = (l + r) / 2;
        mp = Power(m, a);
        if(mp == n) return m;
        if(mp < n) {
            l = m + 1;
        } else {
            r = m - 1;
        }
    }
    return 0;
}

int IsPower(uint64 n) {
    uint8 i, j;
    uint8 cnt = BitCount(n);
    for(i = 2; i < cnt; i++) {
        if(PerfectRoot(i, n, cnt)) return 1;
    }
    return 0;
}

uint8 SmallFactors(uint32 r, uint32 *factors, uint32 *exponents) {
  uint32 i;
  uint32 sqrtr = SquareRoot(r);
  uint8 p = 0;
  for(i = 2; i <= sqrtr; i++) {
      if(r % i == 0) {
          factors[p] = i;
          exponents[p] = 0;
          while(r % i == 0) {
              r /= i;
              exponents[p]++;
          }
          p++;
          sqrtr = SquareRoot(r);
      }
  }
  if(r > 1) {
      factors[p] = r;
      exponents[p] = 1;
      p++;
  }
  return p;
}

uint32 SmallOrder(uint32 n, uint32 r) {
  uint32 i;
  uint32 factors[MAX_FACTORS];
  uint32 exponents[MAX_FACTORS];
  uint32 p;    
  uint32 ans;
  
  ans = 1;
  p = SmallFactors(r, factors, exponents);
  for(i = 0; i < p; i++) {
    if(exponents[i] > 1) {
      ans *= Power(factors[i], exponents[i] - 1);    
    }
    ans *= factors[i] - 1;
  }
  
  p = SmallFactors(ans, factors, exponents);
  for(i = 0; i < p; i++) {
    while(ans % factors[i] == 0) {
      if(PowerMod(n, ans, r) == 1) {
        ans /= factors[i];
      } else {
        break;
      }
    }
    if(ans % factors[i] == 0) {
      ans *= factors[i];
    }
  }
  return ans;
}


uint32 FindR(uint64 n) {
    uint32 r, k;    
    uint8 logn = BitCount(n);
    uint32 maxr = Power(logn, 5);
    uint32 maxk = Power(logn, 2);
    if(maxr < 3) maxr = 3;    
    for(r = 2; r <= maxr; r++) {
        if(gcd(n, r) > 1) continue;   
        k = SmallOrder(n % r, r);
        if(k > maxk) break;
    }
    return r;
}


uint32 SmallPhi(uint32 r) {
    uint32 ans;
    uint32 i;    
    uint32 factors[MAX_FACTORS];
    uint32 exponents[MAX_FACTORS];
    uint32 p;    
    p = SmallFactors(r, factors, exponents);
    ans = 1;
    for(i = 0; i < p; i++) {
        if(exponents[i] > 1) {
          ans *= Power(factors[i], exponents[i] - 1);    
        }
        ans *= factors[i] - 1;
    }
    return ans;
}

int IsPrimeAKS(uint64 n) {
    uint64 i = 0;
    uint32 r = 0;
    uint64 t = 0;
    uint32 logn = 0;
    uint64 maxa = 0;
    uint64 poly[MAXR];

    if(IsPower(n)) return 0;    

    r = FindR(n);

    for(i = 2; i <= r; i++) {
        t = gcd(n, i);
        if(t > 1 && t < n) return 0;
    }

    if(n <= r) return 1;

    logn = BitCount(n);
    maxa = ((uint64)logn) * SquareRoot(SmallPhi(r));
    if(maxa >= n) maxa = n - 1;


    for(i = 1; i <= maxa; i++) {      
        memset(poly, 0, sizeof(uint64) * r);
        poly[0] = i; poly[1] = 1;                
        PowerPoly(poly, n, r);
        if(!CheckPoly(poly, i, n, r)) return 0;
    }

    return 1;
}


int main()
{  
    uint64 n;
    while(scanf("%"SCNu64, &n) != EOF) {
        printf("%s\n", IsPrimeAKS(n) ? "Yes" : "No");
    }
    return 0;
}

参考文献：

两种纯文本终端行下"录屏"的方法

发表评论

方法一: script & scriptreplay
https://aws.amazon.com/cn/premiumsupport/knowledge-center/record-linux-terminal-or-ssh-session/
记录开始: script -a -t timingfile.txt typescript.txt
记录结束: exit
记录重放: scriptreplay --timing=timingfile.txt typescript.txt

方法二: ttyrec & ttyrec2video
https://github.com/jwodder/ttyrec2video
安装: sudo apt-get install ttyrec ffmpeg python3-pip
记录开始: ttyrec record.log
记录结束: exit
记录重放: ttyplay record.log

记录转视频需要用到 ttyrec2video
安装: sudo python3 -m pip install git+https://github.com/jwodder/ttyrec2video.git
修复Bug: sudo pip3 install imageio==2.4.1
转换: ttyrec2video record.log record.mp4

注：pip默认的源在国外很慢, 可以设置使用清华大学的源
sudo pip3 config set global.index-url https://pypi.tuna.tsinghua.edu.cn/simple

一个合并bib文件的小工具

发表评论

因为当写过两篇文章之后，偶尔来个汇报之类的，就需要将之前文章的参考文献汇集在一起处理。然而这些文献本身可能是大量重复的，如何将多个bib文件中的重复文献去掉，仅保留不重复的部分，并最终输出到一个合并之后的bib文件中，是写这个小工具的一个背景。

原来的bib文件已经有自己的citation key，而且自己也习惯于这样的生成方式。虽然有一些文献管理软件也有此功能，但有如下弊端：

1）基于title来去重，而不是citation key去重；
2）bib文件中的其他非文献信息会被清除掉；
3）导出后原来的citation key无法保持原样

考虑到以上因素，所以就写了这个小工具：https://github.com/tpu01yzx/BibtexParser

这个工具的用法很简单，主要是几个参数：

1） -O --output 设置合并后输出到的文件名，缺省是标准输出设备。
2） --onlyregular 是否只输出Reguler（类似article,book,misc之类表示具体的文献记录，而非辅助信息如comment, string）的记录，默认为否。
3） --outputplain 是否输出非@开头的记录，这类块大部分是注释或者用于格式控制的，默认为否。

当然了，仅有这个工具，用起来还是有点麻烦，所以最后送上一个Bat批处理文件，把这个批处理文件，bib合并小工具，以及要合并的所有bib文件放在同一个目录。执行这个批处理文件，就不用每次都去找命令行，对于大部分习惯了Windows的用户来说应该是件好事。

@echo off
SetLocal EnableDelayedExpansion

set allbib=_all.bib
set exec=bib_combiner.exe

echo @comment{this file is generated by BibtexParser.exe} > %allbib%
set list=
for  %%i in (*.bib) do ( 
	if not "%%i" == "%allbib%" (		
		set list=!list! "%%i"
	)
)
%exec% "%allbib%" %list% -O "%allbib%" -R
%exec% --help
set exec=
set allbib=
set list=
pause

给伸手党准备好了，Windows下的打包（包含上面提到的一个工具和一个批处理文件：bib_combiner.exe和一个run.bat），点击这里下载。

一种二元序列的线性复杂度算法的实现

发表评论

#include <stdio.h>
#include <memory.h>

#define MAXL (1000)

typedef unsigned char uchar;

void prints(const uchar *s, const int n)
{
    int i;
    for(i=0;i<n;i++) {printf("%u",s[i]);}
    printf("\n");
}

//implement the Berlekamp–Massey algorithm for binary sequences.
//Please see https://en.wikipedia.org/wiki/Berlekamp%E2%80%93Massey_algorithm#The_algorithm_for_the_binary_field
int bm(const uchar*s, const int n)
{
    int i,j,lc,m,d;
    uchar b[MAXL];
    uchar c[MAXL];
    uchar t[MAXL];
    
    b[0]=1;c[0]=1;
    for(i=1;i<n;i++){b[i]=0;c[i]=0;};
  
    lc=0;m=-1;
    
    for(i=0;i<n;i++) {
        d=s[i];
        for(j=1;j<=lc;j++) { d ^=s[i-j] & c[j]; }
        if(d!=0){
            memcpy(t,c,n);
            for(j=0;j<n-i+m;j++) { c[i-m+j]^=b[j]; }
            
            if(2*lc<=i) {
                lc=i+1-lc;
                m=i;
                memcpy(b,t,n);
            }
        }
    }
  
    return lc;
}

int main()
{
    int n;
    int lc;
    char c;
    int i,j;
    uchar s[MAXL];
    uchar ans_s[MAXL];
    int ans_lc;
    int ans_j;
    
    //if you would like to read the sequence from file, please uncomment the follows line.
    //freopen ("s1.txt","r",stdin);
    //skip the first line
    while((c=getchar())!='\n');
    //read the sequence
    n=0;
    while(n<MAXL){
        c=getchar();
        //the line should be ended with '\n' or ']'.
        if(c=='\n' || c==']') break;
        
        //only '0' and '1' are captured.
        if(c>='0' && c<='1') {
            //convert '0' and '1' to 0 and 1.
            s[n++]=c-'0';
        }
    }
    prints(s,n);
    printf("There are %d bits.\n", n);
    
    lc = bm(s,n);
    printf("The original LC is %d\n", lc);
    
    //make sure that current ans_lc is big enough.
    ans_lc=n; ans_j=-1;
    for(i=0;i<n;i++) {
        s[i] ^=1;
        
        lc=bm(s,n);
        if(lc<ans_lc) {
            //save the current optimal answer;
            ans_lc=lc; ans_j=i;
            memcpy(ans_s,s,n);
            printf("Update ans with lc=%d when  s[%u]:%u->%u.\n", 
            ans_lc, ans_j, s[ans_j]^1, ans_s[ans_j]);
        }
        
        s[i] ^=1;
    }
    
    printf("The minimal 1+lc is %d when s[%u]:%u->%u.\n", 
        1+ans_lc, ans_j, s[ans_j], ans_s[ans_j]);
    prints(ans_s,n);
    
    
    return 0;
}

这是基于C语言实现的wiki上二元序列的线性复杂度，如果需要对应的极小多项式可以从最后的c变量数组中直接转换过来，这并不是什么有难度的事情。

同伴矩阵的方幂

发表评论

最近在研究一个矩阵的方幂，这个矩阵如下

事实上这个矩阵和下面这个多项式存在一一对应的关系,

使用Mathematica计算$A^k$, 代码如下：
<< FiniteFields` Assuming[Element[p, Primes] && Element[n, Integers] && n > 0, F = GF[p, n]; A = {{0, 0, -a[0]}, {1, 0, -a[1]}, {0, 1, -a[2]}}; Print[MatrixForm[ Simplify[MatrixPower[A, k], Element[a[0], F] && Element[a[1], F] && Element[a[2], F] && Element[k, Integers] && k > 1]]]; ];
计算结果为：

显然在化简方面，Mathematica并没有充分利用$a_0$,$a_1$,$a_2$是有限域上的元素的特点，所以这个表达式还是挺复杂的，需要自己进一步做化简. 最终的目标是为了计算$det(A^k-I)$.

突然灵机一动, 改动了一下Mathematica的代码如下:
<< FiniteFields` Assuming[Element[p, Primes] && Element[n, Integers] && n > 0, F = GF[p, n]; A = {{0, 0, -a[0]}, {1, 0, -a[1]}, {0, 1, -a[2]}}; Print[MatrixForm[ Simplify[Det[MatrixPower[A, k] - IdentityMatrix[3]], Element[a[0], F] && Element[a[1], F] && Element[a[2], F] && Element[k, Integers] && k > 1]]]; ];
结果为:

其中$x_1$,$x_2$,$x_3$是方程$f(x)$的三个根.

On Some Papers (5)

发表评论

On Generalized Zero-Difference Balanced Functions
对ZDB函数进行了推广得到G-ZDB函数，同样可以用于构造CCC码和DSS结构。基本思路还是做cyclotomic cosets的划分，不过这种划分是不平衡的，大概是这样

$1,1,\ldots,1,m,m,\ldots,m$

的样子，然后就数一下有多少个1和多少个m。

On Some Papers (4)

发表评论

已经几天没写了，始终有些过意不去，既然现在记得就再写点吧。

A note on the behaviour of the NFS in the medium prime case: Smoothness of Norms
作者就NFS算法的理论算法复杂度是实际算法运行时间进行了试验，并对其中多项式的选择给出了一些观察的结果，这可以指导我们实际操作的时候如何设定参数的大小。其中在筛阶段的光滑概率似乎还没有很好的计算方法，多项式的的预选择就是先估计这个光滑概率。
继续阅读 →

没有姓名的鬼

A Ghost Without Any Names

分类目录归档：研究

AKS素性检测算法的C语言实现

参考文献：

两种纯文本终端行下"录屏"的方法

一个合并bib文件的小工具

一种二元序列的线性复杂度算法的实现

同伴矩阵的方幂

On Some Papers (5)

On Some Papers (4)