如何分解这个变态的卡米歇尔强伪素数？

mathematica

2887148238050771212671429597130393991977609459279722700926516024197432\
3037991527331163289831446392259419778031109293496555784189494417409338\
0561511397999942154241693397290542371100275104208013496673175515285922\
6962916775325475044445856101949404200039904432116776619949629539250452\
6987193290703735640322737012784538991261203092448414947289768854060249\
76768122077071687938121709811322297802059565867

利用同余式分解强伪素数的方法见
https://bbs.emath.ac.cn/forum.ph ... 302&fromuid=865

但是这个整数对同余式似乎没办法用，

这个伪素数可以写成
n=(353*m+1)*(313*m+1)*(m+1)

如果不知道这个分解形式，应该如何分解呢？

mathematica

1. Clear["Global`*"];
   
2. n=2887148238050771212671429597130393991977609459279722700926516024197432\
   
3. 3037991527331163289831446392259419778031109293496555784189494417409338\
   
4. 0561511397999942154241693397290542371100275104208013496673175515285922\
   
5. 6962916775325475044445856101949404200039904432116776619949629539250452\
   
6. 6987193290703735640322737012784538991261203092448414947289768854060249\
   
7. 76768122077071687938121709811322297802059565867;
   
8. Do[ans=Solve[(m+1)*(a*m+1)*(b*m+1)==n,{m},Integers];If[ans!={},Print[{a,b}]],{a,1,500},{b,1,a}]
   

复制代码

我只会用这种笨办法来分解！

nyy

两年后，还是我自己来回答这个问题吧，既然是卡米歇尔数，那么假设最小的素数因子是p（这个p通常都很大），
由于是卡米歇尔数，因此任意选择一个与n互质的基a，那么必然满足a^(n-1)=1(mod n)
假设x=a^((n-1)/2)(mod n),那么必然有x^2=1(mod n),
由于此处的n=3 mod4，那么必然存在对于不同的a，必然存在x既不等于1又不等于n-1（且在1到n-1的整数中，至少有3/4是满足这样的基，要不然miller rabin测试就识别不了卡米歇尔数了），
因此只要把a从2到100000，一个一个穷举就可以了。
这个时候，只要计算GCD[x-1,n]，那么就必然能得到n的因子！

代码

1. Clear["Global`*"];
   
2. n = 288714823805077121267142959713039399197760945927972270092651602419\
   
3. 7432303799152733116328983144639225941977803110929349655578418949441740\
   
4. 9338056151139799994215424169339729054237110027510420801349667317551528\
   
5. 5922696291677532547504444585610194940420003990443211677661994962953925\
   
6. 0452698719329070373564032273701278453899126120309244841494728976885406\
   
7. 024976768122077071687938121709811322297802059565867;
   
8. Do[b=PowerMod[a,(n-1)/2,n];If[(b!=1)&&(b!=n-1)&&(Mod[b^2-1,n]==0),Print[a];Break[]],{a,2,2000}]
   
9. 
   
10. a=307;
    
11. x=PowerMod[a,(n-1)/2,n];
    
12. aa=GCD[x-1,n]
    

复制代码

得到
1047509697104598522420442364894558245396251310534812430290126166254072\
4079869634880456766224539126779375883658239075983560088580357347
剩下的再怎么分解，我就不知道了！

nyy

nyy 发表于 2022-9-19 11:54
两年后，还是我自己来回答这个问题吧，既然是卡米歇尔数，那么假设最小的素数因子是p（这个p通常都很大）， ...

经过验证，上面得到的最大的素数因子，假设这个卡米歇尔数是n，上面的素数因子是aa，然后令m=n/aa

由于是卡米歇尔数，
a^(n-1)=1 mod n，那么必然满足a^(n-1)=1 mod m(模换掉了)，计算x=a^((n-1)/2) mod m，
如果找到的x，既不等于1，又不等于m-1，此时根据卡米歇尔数的性质，必然有x^2=1(mod m),
因此此时，只要计算GCD[x-1,m]，那么就能得到m的因子。

现在的问题就是如何找到满足上面条件的x，答案就是穷举法，一个一个测试。得到对应的a=311，

计算得到素数因子
9288117144298564802198256663229369144731633433354002568861458641289650\
529742764072472996680682001931854537068070342161060361829042067

这个是中间大的素数因子p2

所有的素数因子（从小到大排列）是：
{296744956686855105501541746429053327307719917998530433509950755312768\
38753171770199594238596428121188033664754218345562493168782883,

9288117144298564802198256663229369144731633433354002568861458641289650\
529742764072472996680682001931854537068070342161060361829042067,

1047509697104598522420442364894558245396251310534812430290126166254072\
4079869634880456766224539126779375883658239075983560088580357347}
经过判定，上面的三个因子全部是素数因子！

全部代码如下：

1. Clear["Global`*"];
   
2. n = 288714823805077121267142959713039399197760945927972270092651602419\
   
3. 7432303799152733116328983144639225941977803110929349655578418949441740\
   
4. 9338056151139799994215424169339729054237110027510420801349667317551528\
   
5. 5922696291677532547504444585610194940420003990443211677661994962953925\
   
6. 0452698719329070373564032273701278453899126120309244841494728976885406\
   
7. 024976768122077071687938121709811322297802059565867;
   
8. Do[b=PowerMod[a,(n-1)/2,n];If[(b!=1)&&(b!=n-1)&&(Mod[b^2-1,n]==0),Print[a];Break[]],{a,2,2000}]
   
9. 
   
10. a=307;
    
11. x=PowerMod[a,(n-1)/2,n];
    
12. aa=GCD[x-1,n](*此处得到的是p3,最大的素数因子*)
    
13. m=n/aa;
    
14. 
    
15. Do[b=PowerMod[a,(n-1)/2,m];If[(b!=1)&&(b!=m-1)&&(Mod[b^2-1,m]==0),Print[a];Break[]],{a,2,2000}]
    
16. a=311;
    
17. x=PowerMod[a,(n-1)/2,m];
    
18. bb=GCD[x-1,m](*此处得到的是p2,中间的素数因子*)
    
19. 
    
20. {p1,p2,p3}=Sort[{aa,bb,n/aa/bb}](*素数因此从小到大排序!*)
    
21. 
    

复制代码

nyy

总结一下：
先找witness（这个witness最好也满足费马测试），也就是能证明卡米歇尔数是合数的miller rabin的测试基，
对于素数，永远满足x^2=1(mod p)，所以只有±1的平凡解，所以从liar的嘴里得不到因子，只能去witness去找因子。

多找几个witness基，然后不断开平方找因子。

1. (*分解卡米歇尔大合数*)
   
2. Clear["Global`*"];
   
3. n=1267342751027148335067248740406802390769;(*找到一个卡米歇尔数,等待分解的卡米歇尔数,合数*)
   
4. nold=n;(*备份一下*)
   
5. (*miller rabin test,n0 big odd integer,a0 base*)
   
6. MillerRabin[n0_,a0_]:=Module[{n=n0,a=a0,s,m,t1,k},
   
7.     s=0;m=n-1;While[Mod[m,2]==0,m=m/2;s=s+1];
   
8.     t1=PowerMod[a,m,n];
   
9.     If[t1==1,Return[True]];
   
10.     k=0;While[k<s-1&&t1!=n-1,k=k+1;t1=Mod[t1^2,n]];
    
11.     If[t1==n-1,Return[True],Return[False]]
    
12. ](*miller rabin子函数*)
    
13. (*子函数,指数不断除以2,看能否得到因子.参数:基底;待分解的大整数;同余式的模;大整数里面有几个2*)
    
14. GetFactor[a_,nold_,n_,kmax_]:=Module[{x,k},(*m不指定为局部变量,m全局变量*)
    
15.     Do[x=PowerMod[a,(nold-1)/2^k,n];
    
16.       If[(x!=1)&&(x!=n-1)&&Mod[x^2-1,n]==0,m={GCD[x-1,n],GCD[x+1,n],k};Print[m]],
    
17.       {k,1,kmax}](*nold-1最多4个2,所以k的最大值是4*)
    
18. ]
    
19. (*找出n的miller rabin的witness(没通过MR测试的底能提供因子),因为需要这些数来分解,因此多找几个备用*)
    
20. Select[Range[2,10],Not@MillerRabin[n,#]&](*输出{2, 3, 4, 5, 6, 7, 8, 9, 10}*)
    
21. a=2;(*上面生成的结果,选择2为基*)
    
22. GetFactor[a,nold,n,4]
    
23. PrimeQ[m](*素性判定,看哪个是素因子*)
    
24. p1=m[[2]];(*第一次找到的因子*)
    
25. n=nold/p1;(*将n除以上面找到素因子,这样模就变小了*)
    
26. a=3;(*计算表明2这个基再也给不出因子,因此换成基3*)
    
27. GetFactor[a,nold,n,4]
    
28. PrimeQ[m](*素性判定,看哪个是素因子*)
    

复制代码

nyy

nyy 发表于 2022-9-20 13:39
总结一下：
先找witness（这个witness最好也满足费马测试），也就是能证明卡米歇尔数是合数的miller rabin ...

1. Clear["Global`*"];(*清除所有变量*)
   
2. nn=803837457453639491257079614341942108138837688287558145837488917\
   
3. 522297427376533365218650233616396004545791504202360320876656996\
   
4. 676098728404396540823292873879185086916685732826776177102938969\
   
5. 773947016708230428687109997439976544144845341155872450633409279\
   
6. 022275296229414984230688168540432645753401832978611129896064484\
   
7. 5216191652872597534901
   
8. m=2;n=(nn-1)/2^m
   
9. 
   
10. (*miller rabin子函数*)
    
11. MR[n0_,a0_]:=Module[{n=n0,a=a0,s,m,t1,k},
    
12.     s=0;m=n-1;While[Mod[m,2]==0,m=m/2;s=s+1];
    
13.     t1=PowerMod[a,m,n];
    
14.     If[t1==1,Return[True]];
    
15.     k=0;While[k<s-1&&t1!=n-1,k=k+1;t1=Mod[t1^2,n]];
    
16.     If[t1==n-1,Return[True],Return[False]]
    
17. ]
    
18. Do[If[Not@MR[nn,k],Print[{k}]],{k,1,307}]
    
19. 
    
20. xx=PowerMod[14,n,nn]
    

复制代码

先找到一个witness，也就是14，然后不断开平方找因子

nyy

再来举一个例子！
https://faculty.lynchburg.edu/~nicely/misc/mpzspsp.html

找到
k=11, B=2,3,5,7,11,13,17,19,23,29,31: N=56897193526942024370326972321

这个数是强伪素数

1. Clear["Global`*"];(*清除所有变量*)
   
2. nn=56897193526942024370326972321
   
3. m=5;n=(nn-1)/2^m
   
4. (*miller rabin子函数*)
   
5. MR[n0_,a0_]:=Module[{n=n0,a=a0,s,m,t1,k},
   
6.     s=0;m=n-1;While[Mod[m,2]==0,m=m/2;s=s+1];
   
7.     t1=PowerMod[a,m,n];
   
8.     If[t1==1,Return[True]];
   
9.     k=0;While[k<s-1&&t1!=n-1,k=k+1;t1=Mod[t1^2,n]];
   
10.     If[t1==n-1,Return[True],Return[False]]
    
11. ]
    
12. Do[If[Not@MR[nn,k],Print[{k}]],{k,1,35}]
    
13. 
    
14. xx=PowerMod[14,n,nn]
    
15. 
    
16. {GCD[xx-1,nn],GCD[xx+1,nn]}
    

复制代码

输出结果
Out[28]= 56897193526942024370326972321

Out[29]= 1778037297716938261572717885

\:6B63\:5728\:8BA1\:7B97In[28]:= {14}

\:6B63\:5728\:8BA1\:7B97In[28]:= {22}

\:6B63\:5728\:8BA1\:7B97In[28]:= {35}

Out[32]= 56897193526941611221950985163

Out[33]= {137716125329053, 413148375987157}

nyy

总结一下：
卡米歇尔强伪素数必然满足费马小定理a^(n-1)=1(mod n)
卡米歇尔强伪素数必然能找到通不过miller rabin算法的基（且这种基超过3/4，因此很容易找到）。
因此对于卡米歇尔强伪素数，很容易找到witness基，满足x^2=1 (mod n)，
因此也就容易分解质因数了，因此看起来强大，反而更容易分解！