求解b进制的n位自守数

[mmmcd]

假设有一个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)
如何把相应的自守数找到？

[mmmcd]

如果是30进制(b=30)，该如何计算？

[mmmcd]

如果b=15，如何计算？

[gxqcn]

将 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;
}

[gxqcn]

运行结果（赛扬 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
A573C89A32D937640C3CDC5034466A80B9BE9CC52966827289826B17381971A2679827426208BEACBCCC99176755D63E109A7E21E3D14C2270735D8A2D35BCE8C7AB959750AA465D23DD1399706AB7B83741CE14455499DD8239E69104E269C83E6D7D710C7493832C859405E62E988077D98B54E6678850420AED0C47A56CD8C212CBA31DD477DC97855642B4738D9E4D17DAB7C249B7357B674AB6CD332B07DC9A90A5A4A4ED33D69AC0638BE98A23BD184B04D3AAAE4C61C85614470D0142CC3434AC39EB85BD0CD05CCA97D998AC60C9A0D5C24D5BED96D24920797B333884141DB35A6CDEC9D8E1A6BD900D946787209DE430C59E27E2ADBC5B967CE9555B0C8A5A737D5906667C0728777E0B4212006EE2A205E88C80EB47443E98284D91CDDCD7C89EC604079B2629306B83A5ED620C9638231021B4A7A138E39A29E1A81535041B60E76E44BCB23CBE509DB62BD0A5ABDA5AD6E1DB2CC1132ECEEC2A8907C1624966CC62290C3B85ED74E561A5EB82CB35880E4E6D9B80438C4EBE4D33CAD3229B3C581870235A3414B4DC9D5944247D4AC0B604B9798E333707777C7DC7D2E675EEB94EE187E7603B5367A1227DDB57A9051324B4B356AC6EC2261E7BB22D999B1C6D8B0A25E72523C3D82C03833262618B1D7ABA277DE9A88C0727D966B356A5127E87240700C1723DC9011172195D60D0D20D8AA2047B47CB299E03A85ECD9ED74B4C5954B78DB337414E807C35C383547C6CE17968479E6ABB63579DA310CC2360B0C70DEAAD3C30C3489A231A14525ABA1D8E5E082EB628E3D6A5A9E2602232E67D32737EB6A0487B05D84877A4697710EB5995DEC5C218893757800E5D850833960DCC627BB189284984D9DB29C34E793BC9E2EE9C2ABD763B7D7964C8D978440A41BC17041090493618BCA5445DB086B606571DB12072EC562CE32D57DE225A7B95A9CEB107D5D673C2B84AD55C9406B31E4E29815BD5A2C3C0D0091CDD501B71C259C607E132714C23749B8948215422541DA8A6E390B818D15859A46B85AD680D14D2D39E63B3D916341B84D093B012D9AED88251835939677B1927D06AD9CA062E38A1136897415AE324A235054C9338E9E6CAB280A8A17B8CE6C089898EE2952E599B51A3C00B15A7C61BD6116837604883D73E0462ACB6765A2D8141098ACCD6B564A5C08C06B3A30EE107731365E000CD8D4E3D38E6BB740306DE01EAC9001757771958573C9DAD23C7116B6193AA065D330B4889BD78BD83B60A01A3655B704416E572E410D22A2018D41E3AD14691A8334CE1B40C120B08776B006C2927AE1AEC7EB2D167AD9643D9B28E1B85D00315CB6831B5C2722E88502450191C02DDE82E066354D9DE5C92883813E65CD4CAE7658978C2E9CE8570624D4BDA86

[gxqcn]

关于 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] 输出结果（可避免某些数字与字母的混淆）。

[JTZY]

更正一下（并按个位上的数从小到大排列）：

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次回帖限制，不得不请出我的“结贴专用马甲”:)

[gxqcn]

当 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;
}

[gxqcn]

运行结果（赛扬 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
L1E0JEH8J54LA52HR6LQMROR15N8TFJE3ISCQA32M2IIF1Q4H7GN6E4NMR8LIAPS5RAP5DN6C7M0IKQPD2C4HI0GH2CT0L0L816CCDM833TGGKT8D470GOI2IP9D5SJDMMS3NCDALIMGFTP50ONLRKT4OFO8JNRDE82F9P5JFGND9EJ1LC6C4440DFQSPNICRR7QJOBA7A6Q7ET8JNS770AS877P866SQBIS2OF2348IJ21OFHMQQTCSGGHRPE9GOFM4TGK7LEF25BQHLT6OQQF4OEOBKKK1883D9PEO3PGTESMC2BTHPBL7N8TQNGDNKMLPJRL1GBSNI3AKTELN78OGT47Q0HHE3MEH3QR6T6RHK60NKLAGQNIPLRGB4BGFD6ADGS59B2BDHPL9O98D4TMQ4SIJBDPIQD034993GMS8880L7FMTC1PMD8CI8I6597NTHIDS59EH0RAG0GGMFHFS95PMK4HGECEE7TJ40HQPFBMDRHEJKHMLONFRK6M0B57G5HONB2SFQ6L0OCH65R882TRQPMFKK3ARKT1G6EKTQ5068CK2O289617J6MGKHLLEPBEKJGAOSP3FG451M7FT97HC33GQQ2AJBHA9B3G242DPJMH1BI1ILP6LCL3APSTEA6EF1JS3ACNRDB9FH3KMA7GD3CE6ETHCKT7J15AE6LLJ4B1N7E4KIRI599LD93L9S8HQKRHFHLH2C4ACOCQJS5DNL2MPN94R418GG0D91MAEJCGSCELMO8KDMIFFSJ0N602OLNBGM1QFHSEKSER7RGQLMMFCH6OHRBG1GOB4Q2DOHDGSF8NSI56FQG59B54IKEK568NE9HHI5BHT0AB4TJ8PJJMNENGHELBJRHGN4HA4LK54F2RTCF8N0TH07GKQ5LPJDJI2CG87S16II3GNHL5M4EL5BDAN2PAKB24SG5I0DCDC0HRJSD5K40QACF5GKESKRRGGCC2FD4G0NG76CP8062173JJ97J53T6TQ3378OQL7HKQS8B5FA426GGESJ0I0JO8A425GI59A1NQBM99978IF02CAEGJ4QB4QHN83PD9ED7S903SJE5GNRPRKD1SF3TOMDQ3C7LQEKNRCL4TFCN8KTLANQBC27HIN63H7KEFDI0FP49QC8K0582OR5JAEA4C14JT21EP277N04C1BNL0QFDIOQDTG15BOFEE2E3BTEL80AHBO90HSITMJQR0N7P55KEPQLM004QKE9K6DT1GDKBHT5JN439LRTMA4CMEL0E46RR1TRAQ4HLJ026PRAQNQBFTE6LHALER5QD0JAGOE66KH4400C36Q0TROMJGR2N7CJPH4GK4MP8E6LDGB5E766A7PT5M8GTGDG1318S8CTA9NJIQ827GSPITN5NRF4HQK0ND4PS146R6HM658KLC818H8ONR9QDQB0ORN5FEJKGN0H9PHJA8RMDFDQABO1PK2QJ7F0BAHOOADBQDMDCC7ASH3EPIBKMRFTK2M03QHSBQQSF7D7ST4AO80N7QD3BES809L9PL7AIFMHEO7A9IM9L974D3MFLJRAOA6Q72QORTQ6HBHHM26DA8CGK8GG73229SMCLP0HSHR9JGTP9TBHRF4B9QAI0SS36R6RP21E7KQNS6G097CJ6AQ0F03DS5IOJQOQ3IEOIQFF1CC6E8RHF5J3C68LTF792I70J21LTK1S7D8JCPP6H6BB1R9G8OR7LI6O42SLB9KIAGNMP087TOOMQSRO367DGB87T87TLKG15O07G0SDF7411N893R6OB3GHQI4MKA6A9TARH0NL4MA81RO392M6Q7I8O7LNQ5C1GCP0F7BSR0M7JKMLP3QIM10F74G57S243LKAPEM3CD9BJ4DN5RCQ0GM7QOHDM25P0R53M46HQRC8TG6GO1EJ87B8KBOBLK8AIT9QTKTBOTE8J9AJ4AB6CN47FQ324085HCGG10RT62BT1Q6IRKS95AGEPN2DSAB2A2NEQEFS3MG

[mmmcd]

如果求逆的时候：
x*a mod m = 1
转化为
a*x + m*y = 1

用欧几里德算法求x
时间上会不会比 a乘(m的欧拉函数值)次方 好一点？

[gxqcn]

下面对“用孙子定理计算 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的欧拉函数值)次方 好一点？” 验证了吗？

[gxqcn]

为防止CSDN“私吞”帖子，我已将本贴及本版4051695号贴收录在了我的个人主页中：
http://maths.diy.myrice.com/florilegium/
为保持原滋原味，收录时是从原贴直接按保存得到的，未作任何删改编辑。

[mmmcd]

谁把整个求解过程总结一下，使之适用于处理任意整数：b，n ？

[gxqcn]

to mmmcd(超超):
不知需要的是理论方面的、还是程序方面的总结？（其实现在难度都不大了）

我现在正在开发下一版 HugeCalc，新增或加强了一些关于数论方面的函数，导出函数总数已达400个。新版现正进行各项测试与完善，估计9月初可正式release.

[mmmcd]

程序方面的总结

[gxqcn]

在下一版 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间的！）。

[mathe]

对于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位自守数.

[mathe]

直接计算的话,就是直接计算
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)比较小

[gxqcn]

上述方法是通用算法（但步骤ii用y=(3x^2 - 2x^3) (mod b^(2s))也许更快些）。

用递归方法，算法效率是很高的（有赖于乘法的效率）；
但当计算位数n较小时，据我自测直接解同余方程组效率也许更高些。

另：原定于今天（9月1号）release 新版的 HugeCalc，由于一些事情耽误了进度，无法预期完成，非常抱歉！（为了摆脱我家那台“老爷机”（赛扬466，64M RAM）的限制：只能跑Win98，无法上宽带；于是狠狠心，配置了一台新机器，以后在家调试就方便多了，只是这段时间为此分了不少心。。。）

[languagec]

[mathe]

呵呵，又有人被迫升级电脑了