求解b进制的n位自守数
假设有一个n位b进制数：a，
a×a 在b进制下的最后n位等于a，则称a为自守数。
例如当b=10,n=6时，
a=890625、a=109376 都为所求
890625×890625 ＝ 793212890625
109376×109376 ＝ 11963109376
当b=12,n=6时,
自守数有：1B3854 和 A08369
若给出b，n,(2<=b<=36,1<=n<=2000)
如何把相应的自守数找到？
问题点数：0、回复次数：35Top
1 楼Knuthocean（摘天上的星星）回复于 2005-06-26 22:51:15 得分 0
是不是从最后一位开始求，先保证最后一位满足条件，然后往高位求？
mmmcd (超超) :
你先把自己的思路给大家学习一下吧！Top
2 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-06-27 08:11:44 得分 0
定理1：如果 x 为 b 进制 n 位“自守数”，则 (3x^2 - 2x^3) mod (b^2n) 为 b 进制 2n 位“自守数”。
定理2：如果 x 为 b 进制 n 位“自守数”，则 (x^2 - 1)^2 mod (b^2n) 为 b 进制 2n 位“自守数”。
由自守数有：1B3854（b=12,n=6），通过上述定理可分别得：2B21B61B3854、909A05A08369（n=12）
一般来说，b 进制 n 位“自守数”存在“共轭互补现象”：其和 = b^n + 1Top
3 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-06-27 08:15:40 得分 0
更正一下：
应为“其两两∑≡1(mod b^n)”更严谨，如果将“000..0”与“000..1”也看作一对的话。Top
4 楼mmmcd（超超）回复于 2005-06-27 20:45:23 得分 0
如果b是一个素数p的n次方，b=p^n，是不是除了00...00, 00...01，再没别的自守数了。Top
5 楼mmmcd（超超）回复于 2005-06-27 20:58:50 得分 0
有一种思路：
例如求12进制下的6位自守数x，
x*x mod 12^6 = x => x*x-x mod 12^6 = 0
x mod 12^6 = 0或1
对12^6分解，可以得到同余方程组：
x mod 3^6 = 0 或 1
x mod 4^6 = 0 或 1
用孙子定理解之：
M1=4^6=4096 M1'=118 M1*M1'=483328
M2=3^6=729 M2'=3433 M2*M2'=2502657
x=b1*M1*M1'+b2*M2*M2'mod 12^6
b1=0,b2=0 => x=0 mod 12^6
b1=1,b2=1 => x=1 mod 12^6
b1=1,b2=0 => x=483328 mod 12^6 x写成12进制：1B3854
b1=0,b2=1 => x=2502657 mod 12^6 x写成12进制：A08369
当b=10,n=4过程一样：
x mod 2^4 = 0, 1
x mod 5^4 = 0, 1
M1=625 M1'=1 M1*M1'=625
M2=16 M2'=586 M2*M2'=9376
625，9376就是所求
当n是一个比较大的数，特别是大奇数时，求Mi,Mi'有点郁闷。
Top
6 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-06-28 10:48:35 得分 0
楼主上述方法在 n 比较小时可行；n 较大时在则效率迅速下降。
以下是对比数据（b=10，单位：秒）：
位数n 64 128 256 512 1024
程序A 0.000544 0.002002 0.009507 0.053972 0.344025
程序B 0.010537 0.010596 0.010622 0.010751 0.010897
其中“程序A”为：根据楼主上述算法，调用最新的 HugeCalc V5.0.0.1，并作适当优化所得。
其中“程序B”为：http://maths.myrice.com/download/zss.zip (37.1 KB (38,021 字节))
测试环境为：P4 1.7GHz / 256M RAM / WinXP
B的时间复杂度为：O(nlogn)；A的时间复杂度接近为：O(n^2)，主要耗费在计算 Mi' 上Top
7 楼mmmcd（超超）回复于 2005-06-28 20:22:47 得分 0
要是b=10,
x=(((5^2)^2)^2)^2... 就是其中一解
1000...01 - x 是另一解
但是b=12时，如何快速计算？Top
8 楼Albert_Andrew（阿尔伯特.安德鲁）回复于 2005-06-28 20:51:10 得分 0
使用dfs搜索.从n-1位推导n位.
关键是乘法有一点trick.
zju上通过了,0.69秒,不过这个代码由于不是我的.所以我就不帖了.Top
9 楼Albert_Andrew（阿尔伯特.安德鲁）回复于 2005-06-28 20:53:02 得分 0
因子分解的方法是行不通的.
另外,有17种base只有0,1两个sr数.
我的用时0.01,作弊的。呵呵.Top
10 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-06-29 10:25:08 得分 0
四年前，我曾写过一个程序专门研究“自守数问题”，用的是不断平方的算法，其效率最高可达到 O(n^2)；
去年起，我用递归的算法，并调用比较成熟的 HugeCalc 大数算法库，效率优化到达到 O(nlogn)。
下面给出“快速计算 n 位 12 进制自首数”的源程序，算法用的是孙子定理，
该程序已高度优化，但其效率仍然为 O(n^2)，也就是说不适合 n 比较大的情形；
当 n 比较大时，用我前面给出的两个定理才是上策。
由于调用了最新的 HugeCalc V5.x，所以请按如下步骤操作：
1、从 http://maths.myrice.com/software.htm#02 下载 HugeCalc 的完全版；
2、解压缩 HugeCalc.rar;
3、在 \HugeCalc_Dll_Import\ 的同级目录里新建目录，如：\FactorialTail\
4、在新目录新建一个cpp文件，将如下代码复制进去编译即可。
// AutomorphicNumber12.cpp
//
#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include "../HugeCalc_Dll_Import/HugeCalc.h" // 公共接口
#include "../HugeCalc_Dll_Import/HugeInt.h" // 10进制系统
#include "../HugeCalc_Dll_Import/HugeIntX.h" // 16进制系统
#ifndef _UNICODE
#pragma comment( lib, "../HugeCalc_Dll_Import/HugeCalc.lib" )
#else
#pragma comment( lib, "../HugeCalc_Dll_Import/HugeCalcU.lib" )
#endif
// Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:
// Win32 Debug: Debug Multithreaded DLL
// Win32 Release: Multithreaded DLL
int main(int argc, char* argv[])
{
// printf("Hello World!\n");
cout << "Call " << HugeCalc::GetVersion() << endl << endl;
if ( HC_LICENSE_NONE == HugeCalc::GetLicenseLevel() )
{
cout << endl << "警告：您未通过 HugeCalc.dll 的许可认证！" \
<< endl << endl << "解决方案可选下列方案之一：" \
<< endl << " 一、请移动和[或]修改文件名为：../CopyrightByGuoXianqiang/[../]HugeCalc.exe；" \
<< endl << " 二、请在 HugeCalc.chm 中进行注册（一劳永逸）。" \
<< flush << flush;
}
else
{
UINT32 n, n2, k;
// 也可全用 CHugeInt 类
CHugeIntX Pow3;
CHugeIntX AutomorphicNumber12;
cout << endl << "下面将快速计算 n 位 12 进制自首数（当 n == 0 时退出）:" << endl;
while ( TRUE )
{
cout << endl << "n = ";
cin >> n;
// 检验入参
if ( 0 == n )
{
cout << endl;
break;
}
else
{
if ( 1 > n )
{
n = 1;
}
else if ( 2000 < n )
{
n = 2000;
}
}
// 开始计算
HugeCalc::EnableTimer();
HugeCalc::ResetTimer();
n2 = n * 2;
( Pow3 = 3 ).Pow( n );
AutomorphicNumber12 = 1;
for ( k = 0; n2 != k; ++k )
{
if ( AutomorphicNumber12.IsOdd() )
AutomorphicNumber12 += Pow3;
AutomorphicNumber12 /= 2;
}
AutomorphicNumber12 <<= n2;
// 计算结束
cout << "timer: " << HugeCalc::ShowTimer( FALSE, FALSE ) << " s" << endl;
// 输出结果
cout << AutomorphicNumber12.ConvertToStr( 12, NULL, FS_NORMAL, 4, NULL ) << endl;
}
}
system( "pause" );
return 0;
}
运行结果：
n = 64 timer: 0.000074 s
n = 128 timer: 0.000136 s
n = 256 timer: 0.000338 s
n = 512 timer: 0.001226 s
n = 1024 timer: 0.003597 s
n = 2000 timer: 0.013060 sTop
11 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-06-29 10:27:10 得分 0
n = 2000
timer: 0.013060 s
47636048395B836A0904486626532B4B85A39219938A217815B40278B155AB66312A032A32860528508A49A168A1AB3633B44999A14210617897B92056931868B1A01147995327B8A7667278397328B34B0066594059B529088A2723094B0586B85930B27161724B714A78819998343A542008B613688840566350A402934AB666B3032699764321888B0964B17724054035BB969A9942844662A5B70A8096696A1163298A6B416587A258503B579961499A706B1276029823A0821786072311897651A3A394988149749855562B9BB2454A008498419094474625823663971A246BA464189922790696931539602408BBA2234677080338B9B595650B875498480A55AA12526A709806BABA0845977B6145257504566920515BB45B406A32129506B578048A48B348869986776236383185B0258197560766A1806408A254A99B21A9641461A713608A48337A3787B135890904A408122A56A346187B90B7B4373845384B8BB0A168A051955A5A6214049636234B2BB70845B09187A6586B4045B4235BB9802665935708500070B57A4AA59B45B813B7823948A9A3B532059224A092518238B220B603007019554B3B0183B96331760391545B888B6367736A73181B717540888B4A2B9BA8439637199198657A381914930671B4710A86818B0587726911AB4738192736261BA99813303382280263840220689A9212B33124484B1B1586673751B4214B0337A5B4A67299A2409573401357823BA81628B29304479175751803B809154A1930324718109168829B4067A32043221B5B029425086A572158B8217B6155475B8B80BB82BBA2B067419266233205651B632756A65428673A080213B3BA189B760A1346BAA76842B5893B5A2A74472013415BB842A41389314A24113641567300A59B5AA31A579041BB6A840758938745254441250349528B8251025A33B56A15380A933982702923903B69390A16BB150B709012086308A3A159A8744B60149BAB75662875B2BB5564A01494692961364A34746991955778B76426A5098976127955AA3B682639B99A640525261693B1AA6A6559A02909A78A7A51B647522A9787775262137BAA491655B838312B16A31115620B8425A746B0B0607B021B17B39B2643102A8415590B84A64BAB72943763B43175713386241467870230478624786326813990B6A8448895237780086194532990421231555A0899A245031870A450A2530769A0868262B1400167784423090479B6B210013736A3B4B9843B24253867A44A8A03385043254797647B6B71B00244855A5960B39854AA76313BA1A516379014204335A1AA78B113A4B2A1B2B37062787B50274566720B2A1525678B79013599830402B6468A133A81A3A16986B267B3452B21B61B3854
另一个用 12^2000+1 - 它 即可。Top
12 楼mathe（）回复于 2005-06-29 10:54:34 得分 0
你的程序是计算4^(-n) (mod 3^n),不需要这么复杂.
首先
2^(-1)= (1+3^n)/2 (mod 3^n)
然后计算
(1+3^n)/2 关于模3^n的2n次方就可以了.
mul=(1+3^n)/2;
base=1;
m=3^n;
times=2*n;
while(times){
if(is_odd(times)){
base*=mul;
base %=m;
}
mul*=mul;
mul%=m;
}
Top
13 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-06-29 11:09:45 得分 0
其实，这正是我优化的一处地方（当 n 比较小时），
因为，用我的算法可以避免大数的乘法（或平方）及取模，而这通常是比较耗时的！Top
14 楼mathe（）回复于 2005-06-29 11:23:22 得分 0
可能有些道理,不过n大一些时,明显用乘法和取模快,是O(n log^2(n))的复杂堵Top
15 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-06-29 11:53:16 得分 0
下面再优化一些：尽量减少大整数左右移位的次数，
计算 n = 2000，由原来的 0.013060s 减少到 0.008394s！
请将如下源代码替换原来相应之处：
AutomorphicNumber12 = 1;
k = 0;
while ( TRUE )
{
AutomorphicNumber12 += Pow3;
// 注意：AutomorphicNumber12 必须选用 CHugeIntX 类
UINT32 u32BitIndex = 0;
while ( !AutomorphicNumber12.IsBitOne( u32BitIndex ) )
{
++u32BitIndex;
}
k += u32BitIndex;
if ( k < n2 )
{
AutomorphicNumber12 >>= u32BitIndex;
}
else
{
u32BitIndex -= ( k - n2 );
AutomorphicNumber12 <<= n2 - u32BitIndex;
break;
}
}Top
16 楼mmmcd（超超）回复于 2005-07-01 22:02:38 得分 0
如果是30进制(b=30)，该如何计算？Top
17 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-07-02 12:45:29 得分 0
当 b = 30，求 n 位自守数 x，可解如下同余方程组获得之一：
x≡0（mod 2^n），x≡1（mod 15^n）
下面据此给出“快速计算 n 位 30 进制自守数”的源程序，算法用的是孙子定理，
该程序已高度优化，但其效率仍然为 O(n^2)，也就是说不适合 n 比较大的情形；
当 n 比较大时，用我前面给出的两个定理才是上策。
由于调用了最新的 HugeCalc V5.x，所以请按如下步骤操作：
1、从 http://maths.myrice.com/software.htm#02 下载 HugeCalc 的完全版；
2、解压缩 HugeCalc.rar;
3、在 \HugeCalc_Dll_Import\ 的同级目录里新建目录，如：\AutomorphicNumber30\
4、在新目录新建一个AutomorphicNumber30.cpp文件，将如下代码复制进去编译即可。
// AutomorphicNumber30.cpp
//
#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include "../HugeCalc_Dll_Import/HugeCalc.h" // 公共接口
#include "../HugeCalc_Dll_Import/HugeInt.h" // 10进制系统
#include "../HugeCalc_Dll_Import/HugeIntX.h" // 16进制系统
#ifndef _UNICODE
#pragma comment( lib, "../HugeCalc_Dll_Import/HugeCalc.lib" )
#else
#pragma comment( lib, "../HugeCalc_Dll_Import/HugeCalcU.lib" )
#endif
// Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:
// Win32 Debug: Debug Multithreaded DLL
// Win32 Release: Multithreaded DLL
int main(int argc, char* argv[])
{
// printf("Hello World!\n");
cout << "Call " << HugeCalc::GetVersion() << endl << endl;
if ( HC_LICENSE_NONE == HugeCalc::GetLicenseLevel() )
{
cout << endl << "警告：您未通过 HugeCalc.dll 的许可认证！" \
<< endl << endl << "解决方案可选下列方案之一：" \
<< endl << " 一、请移动和[或]修改文件名为：../CopyrightByGuoXianqiang/[../]HugeCalc.exe；" \
<< endl << " 二、请在 HugeCalc.chm 中进行注册（一劳永逸）。" \
<< flush << flush;
}
else
{
UINT32 n, k;
CHugeIntX Pow15;
CHugeIntX AutomorphicNumber30;
cout << endl << "下面将快速计算 n 位 30 进制自守数（1 ≤n≤2000；当 n == 0 时退出）..." << endl;
while ( TRUE )
{
cout << endl << "n = ";
cin >> n;
// 检验入参
if ( 0 == n )
{
cout << endl;
break;
}
else
{
if ( 1 > n )
{
n = 1;
}
else if ( 2000 < n )
{
n = 2000;
}
}
// 开始计算
HugeCalc::EnableTimer();
HugeCalc::ResetTimer();
( Pow15 = 15 ).Pow( n );
AutomorphicNumber30 = 1;
k = 0;
while ( TRUE )
{
AutomorphicNumber30 += Pow15;
UINT32 u32BitIndex = 0;
while ( !AutomorphicNumber30.IsBitOne( u32BitIndex ) )
{
++u32BitIndex;
}
k += u32BitIndex;
if ( k < n )
{
AutomorphicNumber30 >>= u32BitIndex;
}
else
{
u32BitIndex -= ( k - n );
AutomorphicNumber30 <<= n - u32BitIndex;
break;
}
}
// 计算结束
cout << "timer: " << HugeCalc::ShowTimer( FALSE, FALSE ) << " s" << endl;
// 输出结果
cout << AutomorphicNumber30.ConvertToStr( 30, NULL, FS_NORMAL, 4, NULL ) \
<< endl;
}
}
system( "pause" );
return 0;
}Top
18 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-07-02 12:47:50 得分 0
运行结果（赛扬 466MHz / 64MB RAM / Win98Sec）：
n = 256, timer: 0.000818 s
n = 512, timer: 0.002289 s
n = 1024, timer: 0.008472 s
n = 2000, timer: 0.031413 s
L1E0JEH8J54LA52HR6LQMROR15N8TFJE3ISCQA32M2IIF1Q4H7GN6E4NMR8LIAPS5RAP5DN6C7M0IKQPD2C4HI0GH2CT0L0L816CCDM833TGGKT8D470GOI2IP9D5SJDMMS3NCDALIMGFTP50ONLRKT4OFO8JNRDE82F9P5JFGND9EJ1LC6C4440DFQSPNICRR7QJOBA7A6Q7ET8JNS770AS877P866SQBIS2OF2348IJ21OFHMQQTCSGGHRPE9GOFM4TGK7LEF25BQHLT6OQQF4OEOBKKK1883D9PEO3PGTESMC2BTHPBL7N8TQNGDNKMLPJRL1GBSNI3AKTELN78OGT47Q0HHE3MEH3QR6T6RHK60NKLAGQNIPLRGB4BGFD6ADGS59B2BDHPL9O98D4TMQ4SIJBDPIQD034993GMS8880L7FMTC1PMD8CI8I6597NTHIDS59EH0RAG0GGMFHFS95PMK4HGECEE7TJ40HQPFBMDRHEJKHMLONFRK6M0B57G5HONB2SFQ6L0OCH65R882TRQPMFKK3ARKT1G6EKTQ5068CK2O289617J6MGKHLLEPBEKJGAOSP3FG451M7FT97HC33GQQ2AJBHA9B3G242DPJMH1BI1ILP6LCL3APSTEA6EF1JS3ACNRDB9FH3KMA7GD3CE6ETHCKT7J15AE6LLJ4B1N7E4KIRI599LD93L9S8HQKRHFHLH2C4ACOCQJS5DNL2MPN94R418GG0D91MAEJCGSCELMO8KDMIFFSJ0N602OLNBGM1QFHSEKSER7RGQLMMFCH6OHRBG1GOB4Q2DOHDGSF8NSI56FQG59B54IKEK568NE9HHI5BHT0AB4TJ8PJJMNENGHELBJRHGN4HA4LK54F2RTCF8N0TH07GKQ5LPJDJI2CG87S16II3GNHL5M4EL5BDAN2PAKB24SG5I0DCDC0HRJSD5K40QACF5GKESKRRGGCC2FD4G0NG76CP8062173JJ97J53T6TQ3378OQL7HKQS8B5FA426GGESJ0I0JO8A425GI59A1NQBM99978IF02CAEGJ4QB4QHN83PD9ED7S903SJE5GNRPRKD1SF3TOMDQ3C7LQEKNRCL4TFCN8KTLANQBC27HIN63H7KEFDI0FP49QC8K0582OR5JAEA4C14JT21EP277N04C1BNL0QFDIOQDTG15BOFEE2E3BTEL80AHBO90HSITMJQR0N7P55KEPQLM004QKE9K6DT1GDKBHT5JN439LRTMA4CMEL0E46RR1TRAQ4HLJ026PRAQNQBFTE6LHALER5QD0JAGOE66KH4400C36Q0TROMJGR2N7CJPH4GK4MP8E6LDGB5E766A7PT5M8GTGDG1318S8CTA9NJIQ827GSPITN5NRF4HQK0ND4PS146R6HM658KLC818H8ONR9QDQB0ORN5FEJKGN0H9PHJA8RMDFDQABO1PK2QJ7F0BAHOOADBQDMDCC7ASH3EPIBKMRFTK2M03QHSBQQSF7D7ST4AO80N7QD3BES809L9PL7AIFMHEO7A9IM9L974D3MFLJRAOA6Q72QORTQ6HBHHM26DA8CGK8GG73229SMCLP0HSHR9JGTP9TBHRF4B9QAI0SS36R6RP21E7KQNS6G097CJ6AQ0F03DS5IOJQOQ3IEOIQFF1CC6E8RHF5J3C68LTF792I70J21LTK1S7D8JCPP6H6BB1R9G8OR7LI6O42SLB9KIAGNMP087TOOMQSRO367DGB87T87TLKG15O07G0SDF7411N893R6OB3GHQI4MKA6A9TARH0NL4MA81RO392M6Q7I8O7LNQ5C1GCP0F7BSR0M7JKMLP3QIM10F74G57S243LKAPEM3CD9BJ4DN5RCQ0GM7QOHDM25P0R53M46HQRC8TG6GO1EJ87B8KBOBLK8AIT9QTKTBOTE8J9AJ4AB6CN47FQ324085HCGG10RT62BT1Q6IRKS95AGEPN2DSAB2A2NEQEFS3MGTop
19 楼mmmcd（超超）回复于 2005-07-02 15:30:26 得分 0
如果b=15，如何计算？Top
20 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-07-02 18:06:16 得分 0
将 b 分解为 b=p*q（p,q为互素且均大于1的整数），而后通过解同余方程组得到一组解。
该解与p,q相关；即：相同的b，不同的p,q，结果可能不一致。
无论 b 的奇偶性如何，如果用孙子定理，所依据的理论是一致的，
但偶数情形因可将 b 分解为 b=2^r*m 型（m为>1的奇数），从而可用一些“移位”运算优化。
下面据此给出“下面将用孙子定理计算 n 位 15 进制自守数”的源程序，
该程序已初步优化，但其效率仍然为 O(n^2)，也就是说不适合 n 比较大的情形；
当 n 比较大时，用我前面给出的两个定理才是上策。
由于调用了最新的 HugeCalc V5.x，所以请按如下步骤操作：
1、从 http://maths.myrice.com/software.htm#02 下载 HugeCalc 的完全版；
2、解压缩 HugeCalc.rar;
3、在 \HugeCalc_Dll_Import\ 的同级目录里新建目录，如：\AutomorphicNumber15\
4、在新目录新建一个AutomorphicNumber15.cpp文件，将如下代码复制进去编译即可。
// AutomorphicNumber15.cpp
//
#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include "../HugeCalc_Dll_Import/HugeCalc.h" // 公共接口
#include "../HugeCalc_Dll_Import/HugeInt.h" // 10进制系统
#include "../HugeCalc_Dll_Import/HugeIntX.h" // 16进制系统
#ifndef _UNICODE
#pragma comment( lib, "../HugeCalc_Dll_Import/HugeCalc.lib" )
#else
#pragma comment( lib, "../HugeCalc_Dll_Import/HugeCalcU.lib" )
#endif
// Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:
// Win32 Debug: Debug Multithreaded DLL
// Win32 Release: Multithreaded DLL
int main(int argc, char* argv[])
{
// printf("Hello World!\n");
cout << "Call " << HugeCalc::GetVersion() << endl << endl;
if ( HC_LICENSE_NONE == HugeCalc::GetLicenseLevel() )
{
cout << endl << "警告：您未通过 HugeCalc.dll 的许可认证！" \
<< endl << endl << "解决方案可选下列方案之一：" \
<< endl << " 一、请移动和[或]修改文件名为：../CopyrightByGuoXianqiang/[../]HugeCalc.exe；" \
<< endl << " 二、请在 HugeCalc.chm 中进行注册（一劳永逸）。" \
<< flush << flush;
}
else
{
UINT32 k;
BOOL bOut2File = FALSE;
cout << endl << "下面将用孙子定理计算 n 位 15 进制自守数..." << endl;
cout << endl << "请问是否将结果保存在文件中(y or n)？";
k = cin.get();
bOut2File = ( 'y' == k || 'Y' == k );
cout << endl << "请输入位数 n（1 ≤n≤2000；当 n == 0 时退出）:" << endl;
while ( TRUE )
{
UINT32 n;
cout << endl << "n = ";
cin >> n;
// 检验入参
if ( 0 == n )
{
cout << endl;
break;
}
else
{
if ( 1 > n )
{
n = 1;
}
else if ( 2000 < n )
{
n = 2000;
}
}
// 开始计算
HugeCalc::EnableTimer();
HugeCalc::ResetTimer();
// m1 = p^n, m2 = q^n,
// x ≡ 0 ( mod m1 ), x ≡ 1 ( mod m2 )
// x ≡ k*m1 ( mod m1*m2 )，其中 k 满足 k*m1 ≡ 1 ( mod m2 )
// 推论：同样的 b=p*q，不同的分解 p,q 可能得到不同的结果！
// 针对 b=15, 令 p=3, q=5,
CHugeIntX m1( 3 );
CHugeIntX m2( 5 );
m1.Pow( n ); // m1 = 3^n
m2.Pow( n ); // m2 = 5^n
// 下面计算 k = m1^-1 mod m2
// 用数论知识，可推导出：k ≡ 3^[5^(n-1)*4-n] (mod m2)
CHugeIntX phi;
(( phi = m2 ) /= 5 ) <<= 2;
CHugeIntX AutomorphicNumber15;
AutomorphicNumber15.PowMod( 3, phi-=n, m2 ) \
*= m1;
// 计算结束
cout << "timer: " << HugeCalc::ShowTimer( FALSE, FALSE ) << "s" << endl;
if ( !bOut2File )
{
// 输出结果
cout << AutomorphicNumber15.ConvertToStr( 15, NULL, FS_NORMAL, 4, NULL ) \
<< endl;
}
else
{
// 打开文件
ofstream outputFile( "AutomorphicNumber15.txt", ios::ate );
if ( !outputFile )
{
cerr << "File could not be opened" << endl;
exit( 1 );
}
// 输出结果
outputFile << "n = " << n << ", timer: " \
<< HugeCalc::ShowTimer( FALSE, FALSE ) << "s" << endl;
outputFile << \
AutomorphicNumber15.ConvertToStr( 15, NULL, FS_NORMAL, 4, NULL ) \
<< endl << endl;
outputFile.close();
}
}
}
system( "pause" );
return 0;
}
Top
21 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-07-02 18:07:04 得分 0
运行结果（赛扬 466MHz / 64MB RAM / Win98Sec）：
n = 256, timer: 0.189361s
1757771958573C9DAD23C7116B6193AA065D330B4889BD78BD83B60A01A3655B704416E572E410D22A2018D41E3AD14691A8334CE1B40C120B08776B006C2927AE1AEC7EB2D167AD9643D9B28E1B85D00315CB6831B5C2722E88502450191C02DDE82E066354D9DE5C92883813E65CD4CAE7658978C2E9CE8570624D4BDA86
n = 512, timer: 0.864722s
D14D2D39E63B3D916341B84D093B012D9AED88251835939677B1927D06AD9CA062E38A1136897415AE324A235054C9338E9E6CAB280A8A17B8CE6C089898EE2952E599B51A3C00B15A7C61BD6116837604883D73E0462ACB6765A2D8141098ACCD6B564A5C08C06B3A30EE107731365E000CD8D4E3D38E6BB740306DE01EAC9001757771958573C9DAD23C7116B6193AA065D330B4889BD78BD83B60A01A3655B704416E572E410D22A2018D41E3AD14691A8334CE1B40C120B08776B006C2927AE1AEC7EB2D167AD9643D9B28E1B85D00315CB6831B5C2722E88502450191C02DDE82E066354D9DE5C92883813E65CD4CAE7658978C2E9CE8570624D4BDA86
n = 1024, timer: 4.759449s
240700C1723DC9011172195D60D0D20D8AA2047B47CB299E03A85ECD9ED74B4C5954B78DB337414E807C35C383547C6CE17968479E6ABB63579DA310CC2360B0C70DEAAD3C30C3489A231A14525ABA1D8E5E082EB628E3D6A5A9E2602232E67D32737EB6A0487B05D84877A4697710EB5995DEC5C218893757800E5D850833960DCC627BB189284984D9DB29C34E793BC9E2EE9C2ABD763B7D7964C8D978440A41BC17041090493618BCA5445DB086B606571DB12072EC562CE32D57DE225A7B95A9CEB107D5D673C2B84AD55C9406B31E4E29815BD5A2C3C0D0091CDD501B71C259C607E132714C23749B8948215422541DA8A6E390B818D15859A46B85AD680D14D2D39E63B3D916341B84D093B012D9AED88251835939677B1927D06AD9CA062E38A1136897415AE324A235054C9338E9E6CAB280A8A17B8CE6C089898EE2952E599B51A3C00B15A7C61BD6116837604883D73E0462ACB6765A2D8141098ACCD6B564A5C08C06B3A30EE107731365E000CD8D4E3D38E6BB740306DE01EAC9001757771958573C9DAD23C7116B6193AA065D330B4889BD78BD83B60A01A3655B704416E572E410D22A2018D41E3AD14691A8334CE1B40C120B08776B006C2927AE1AEC7EB2D167AD9643D9B28E1B85D00315CB6831B5C2722E88502450191C02DDE82E066354D9DE5C92883813E65CD4CAE7658978C2E9CE8570624D4BDA86
n = 2000, timer: 16.682695s
A573C89A32D937640C3CDC5034466A80B9BE9CC52966827289826B17381971A2679827426208BEACBCCC99176755D63E109A7E21E3D14C2270735D8A2D35BCE8C7AB959750AA465D23DD1399706AB7B83741CE14455499DD8239E69104E269C83E6D7D710C7493832C859405E62E988077D98B54E6678850420AED0C47A56CD8C212CBA31DD477DC97855642B4738D9E4D17DAB7C249B7357B674AB6CD332B07DC9A90A5A4A4ED33D69AC0638BE98A23BD184B04D3AAAE4C61C85614470D0142CC3434AC39EB85BD0CD05CCA97D998AC60C9A0D5C24D5BED96D24920797B333884141DB35A6CDEC9D8E1A6BD900D946787209DE430C59E27E2ADBC5B967CE9555B0C8A5A737D5906667C0728777E0B4212006EE2A205E88C80EB47443E98284D91CDDCD7C89EC604079B2629306B83A5ED620C9638231021B4A7A138E39A29E1A81535041B60E76E44BCB23CBE509DB62BD0A5ABDA5AD6E1DB2CC1132ECEEC2A8907C1624966CC62290C3B85ED74E561A5EB82CB35880E4E6D9B80438C4EBE4D33CAD3229B3C581870235A3414B4DC9D5944247D4AC0B604B9798E333707777C7DC7D2E675EEB94EE187E7603B5367A1227DDB57A9051324B4B356AC6EC2261E7BB22D999B1C6D8B0A25E72523C3D82C03833262618B1D7ABA277DE9A88C0727D966B356A5127E87240700C1723DC9011172195D60D0D20D8AA2047B47CB299E03A85ECD9ED74B4C5954B78DB337414E807C35C383547C6CE17968479E6ABB63579DA310CC2360B0C70DEAAD3C30C3489A231A14525ABA1D8E5E082EB628E3D6A5A9E2602232E67D32737EB6A0487B05D84877A4697710EB5995DEC5C218893757800E5D850833960DCC627BB189284984D9DB29C34E793BC9E2EE9C2ABD763B7D7964C8D978440A41BC17041090493618BCA5445DB086B606571DB12072EC562CE32D57DE225A7B95A9CEB107D5D673C2B84AD55C9406B31E4E29815BD5A2C3C0D0091CDD501B71C259C607E132714C23749B8948215422541DA8A6E390B818D15859A46B85AD680D14D2D39E63B3D916341B84D093B012D9AED88251835939677B1927D06AD9CA062E38A1136897415AE324A235054C9338E9E6CAB280A8A17B8CE6C089898EE2952E599B51A3C00B15A7C61BD6116837604883D73E0462ACB6765A2D8141098ACCD6B564A5C08C06B3A30EE107731365E000CD8D4E3D38E6BB740306DE01EAC9001757771958573C9DAD23C7116B6193AA065D330B4889BD78BD83B60A01A3655B704416E572E410D22A2018D41E3AD14691A8334CE1B40C120B08776B006C2927AE1AEC7EB2D167AD9643D9B28E1B85D00315CB6831B5C2722E88502450191C02DDE82E066354D9DE5C92883813E65CD4CAE7658978C2E9CE8570624D4BDA86Top
22 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-07-02 21:42:28 得分 0
关于 n 位 b 进制“自守数”的个数，我进行了简单的推导，得到如下结论：
定理：如果自然数 b 的质因数个数为 k，则 n 位 b 进制“自守数”的个数有且仅有 2^k 个（包括平凡解 0...00、0...01）
如当 n=8，b=30(=2*3*5) 时，有且仅有如下 8 组解：
1) 0000 0000
2) bq6l b2j6
3) eqef s3mg
4) qtm5 csql
5) 307O h139
6) f3fe 1q7e
7) i3n8 irao
8) 0000 0001
注：HugeCalc 直接支持用小写字符集 [0-9a-z] 输出结果（可避免某些数字与字母的混淆）。Top
23 楼JTZY（GxQ结贴专用马甲）回复于 2005-07-02 22:02:46 得分 0
更正一下（并按个位上的数从小到大排列）：
1) 0000 0000
2) 0000 0001
3) bq6l b2j6
4) 307O h139
5) f3fe 1q7e
6) eqef s3mg
7) qtm5 csql
8) i3n8 iram
注：因受连续3次回帖限制，不得不请出我的“结贴专用马甲”:)Top
24 楼mmmcd（超超）回复于 2005-07-03 10:41:15 得分 0
如果求逆的时候：
x*a mod m = 1
转化为
a*x + m*y = 1
用欧几里德算法求x
时间上会不会比 a乘(m的欧拉函数值)次方 好一点？Top
25 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-07-05 10:42:40 得分 0
下面对“用孙子定理计算 n 位 15 进制自守数”的算法再优化！
计算 3^(-n) mod 5^n 的优化算法：
∵5^n mod 3 = (-1)^n mod 3，令 m = 5^n
∴ 当 n 为奇数时，3^(-n) mod m = ((1+m)/3)^n mod m；
当 n 为偶数时，3^(-n) mod m = ((1+2*m)/3)^n mod m；
请注意，此时上面右式中已全部为普通的整型运算。该算法比用Euler定理要快很多！
n位15进制自守数 | 优化前 | 优化后
----------------+----------+-----------
n = 256, timer: 0.054091s | 0.000339s
n = 512, timer: 0.251948s | 0.000686s
n = 1024, timer: 2.174461s | 0.001965s
n = 2000, timer: 7.916857s | 0.007955s
（P4 1.7Ghz / 256M RAM / WinXP）
附：优化后的源代码为：
// 。。。
// 开始计算
HugeCalc::EnableTimer();
HugeCalc::ResetTimer();
// 注意：为什么全部不用 CHugeIntX 类了？请感兴趣的朋友揣摩一下:)
CHugeInt AutomorphicNumber15;
const CHugeInt m2( CHugeInt( 5 ).Pow( n ) ); // m2 = 5^n
CHugeInt m1( ( 2 - ( n & 1 )) * m2 );
( AutomorphicNumber15.PowMod( ( ++m1 ) /= 3, n, CHugeInt( 10 ).Pow( n ) ) \
%= m2 ) \
*= CHugeInt( 3 ).Pow( n ) ;
// 计算结束
// 。。。
to mmmcd(超超):
“用欧几里德算法求x时间上会不会比 a乘(m的欧拉函数值)次方 好一点？” 验证了吗？Top
26 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-07-07 08:44:22 得分 0
为防止CSDN“私吞”帖子，我已将本贴及本版4051695号贴收录在了我的个人主页中：
http://maths.diy.myrice.com/florilegium/
为保持原滋原味，收录时是从原贴直接按保存得到的，未作任何删改编辑。Top
27 楼mmmcd（超超）回复于 2005-07-31 08:45:03 得分 0
谁把整个求解过程总结一下，使之适用于处理任意整数：b，n ？Top
28 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-08-01 08:03:41 得分 0
to mmmcd(超超):
不知需要的是理论方面的、还是程序方面的总结？（其实现在难度都不大了）
我现在正在开发下一版 HugeCalc，新增或加强了一些关于数论方面的函数，导出函数总数已达400个。新版现正进行各项测试与完善，估计9月初可正式release.Top
29 楼mmmcd（超超）回复于 2005-08-01 22:23:40 得分 0
程序方面的总结Top
30 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-08-02 07:58:10 得分 0
在下一版 HugeCalc V5.1.0.1(预期2005.09.01正式release)中，新增了 G.c.dex（扩展最大公约数）、InvertMod（倒数取模）等数论函数，可更方便对该问题的解决。
/*=============================================================
回复人： mmmcd(超超) ( ) 信誉：124 2005-7-3 10:41:15 得分: 0
如果求逆的时候：
x*a mod m = 1
转化为
a*x + m*y = 1
用欧几里德算法求x
时间上会不会比 a乘(m的欧拉函数值)次方 好一点？
=============================================================*/
经过测试，在某些时候，确实要好一些；但当 x、a 均为 n 次幂时，即：x=x'^n，a=a'^n，则应先求出 x'a'≡1(mod m)，而后计算 x = x'^n mod m 将更快！
如果需要，我可以将任意进制、任意位数的“自守数”的程序抛出来（但必须依赖于特定的大数运算库；好在是 HugeCalc 早就具备快速任意进制转换功能，且不仅仅是局限于2～36间的！）。Top
31 楼mathe（）回复于 2005-09-01 08:02:37 得分 0
对于b进制,计算n位自守数,可以如下:
i)首先计算b进制1位自守数
对于b的任何一种分解b=u*v, 其中(u,v)=1 (如果b有k个不同素因子,总共有2^k种分解方法, u*v 和v*u看成两种分解方法)
我们可以求解方程: x=0 (mod u), x=1 (mod v)
对于这个x,有x^2=x (mod b),即x是一位自守数
ii)对于每个b进制s位自首数x
x=0 (mod u^s) x=1 (mod v^s)
y=(x^2-1)^2 (mod b^(2s))是b进制2s位自守数
而且 y=0 (mod v^(2s)), y=1 (mod u^(2s))
通过这种递归方法,可以计算出2^k个n位自守数.Top
32 楼mathe（）回复于 2005-09-01 08:11:22 得分 0
直接计算的话,就是直接计算
x=0 (mod u^n),
x=1 (mod v^n)
设
x=a*u^n
所以
a=u^(-n) (mod v^n)
我们需要计算a,
可以采用先计算 t=u^(-1) (mod v^n) (欧几里得算法)
然后计算t^n (mod v^n), 需要log(n)次n位数得乘法和求模运算
采用上面一种算法,也是需要将x连续平方和求模log(n)次,但是上面计算过程中开始几步需要整数的运算位数都远远低于n位.
所以如果假设n位整数的乘法和求模需要的时间复杂度都是 O(n*log(n)),那么上面算法需要的时间复杂度还是O(n*log(n)),但是这个算法需要O(n*log^2(n)),在n较大时速度要慢一些,但是总体差别不是很明显,毕竟log(n)比较小Top
33 楼gxqcn（★) HugeCalc ← http://hugecalc.ik8.com/ (☆）回复于 2005-09-01 08:25:23 得分 0
上述方法是通用算法（但步骤ii用y=(3x^2 - 2x^3) (mod b^(2s))也许更快些）。
用递归方法，算法效率是很高的（有赖于乘法的效率）；
但当计算位数n较小时，据我自测直接解同余方程组效率也许更高些。
另：原定于今天（9月1号）release 新版的 HugeCalc，由于一些事情耽误了进度，无法预期完成，非常抱歉！（为了摆脱我家那台“老爷机”（赛扬466，64M RAM）的限制：只能跑Win98，无法上宽带；于是狠狠心，配置了一台新机器，以后在家调试就方便多了，只是这段时间为此分了不少心。。。）Top
34 楼languagec（各有所求）回复于 2005-09-02 23:12:08 得分 0
35 楼mathe（）回复于 2005-09-07 07:02:00 得分 0
呵呵，又有人被迫升级电脑了Top
