那么,接下来的困难得多的问题是,对于小于N的自然数x=a的大概有多少呢? 这个问题比较难。可靠的方法很难想到。
不过我们可以假设平方数各位数都是均匀分布的,给出一个近似统计数据来(计算公式估计还是有点麻烦)。
而更加精确的模型可以对平方数的首位和末位的分布各个更加精确一些的统计模型,比如末位应该是$1/5$的概率取1,4,9,6,然后$1/10$的概率取0和5 刚才把 36# 算法再次优化,
完全消除了除法及64位的运算,
算出 X=130,仅需 11.243 s;遍历18位以内的平方数,仅需 126.074 s
大家看看还有哪里可以优化的? 函数SumOfDigits可以提速的,比如我们事先计算1~9999所有数字的sumofdigit,并且保存
然后对于每个整数,每四位的对应和再累加,可以减少一些乘除次数. 类似想法我在 38# 曾提到过。:handshake :lol 我想
折半算法更好吧
由于这个求数位和是并行操作
所以考虑如果采取精细的汇编
假设考虑到32位十进制数字
则首先除以10^16
分开,分开的部分
分别再除以10^8
则只需要
一次除以10^16
二次除以10^8
四次除以10^4 :lol
更黑的方法在下面
如果一开始就存储
1-9999内的乘法
共5*10^7个结果
则可以考虑整体的数字
采取10000进制
呵呵 :lol
但上述方法只能对计算超过2^64以上的
数字的平方和才有好的效果
更加优化的算法
下面的代码可以遍历全部的64bit内的完全平方数,且效率非常高:#include <stdlib.h>#include <stdio.h>
#include "../../../HugeCalc_API/CppAPI/Include/HugeCalc.h" // 公共接口
#pragma message( "automatic link to ../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
#pragma comment( lib, "../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
#define TEN_POW2 100UL
#define TEN_POW4 10000UL
int main( void )
{
const UINT32 cycleDelta[ 4 ] = { 1, 3, 3, 2}; // 0, 1, 4, 7
UINT32 table[ TEN_POW4 ];
UINT32 mark[ 9*4*5 + 1 ];
UINT32 value[ 5 ], delta[ 5 ];
UINT32 * p = table;
UINT32 i, j, k, m, s;
UINT32 i2, i3, i4, s2, s3;
HugeCalc::EnableTimer();
HugeCalc::ResetTimer();
for ( i4 = 0; i4 < 10; ++i4 )
{
for ( i3 = 0; i3 < 10; ++i3 )
{
s3 = i4 + i3;
for ( i2 = 0; i2 < 10; ++i2 )
{
s2 = s3 + i2;
for ( i = 0; i < 10; ++i )
{
*p++ = s2 + i;
}
}
}
}
memset( mark, 0, sizeof( mark ));
memset( value, 0, sizeof( value ));
memset( delta, 0, sizeof( delta ));
k = 1;
j = TEN_POW2;
value[ 0 ] = 0;
delta[ 0 ] = 1;
m = 1;
for ( ; ; )
{
s = 0;
for ( i = 0; i < m; ++i )
{
value[ i ] += delta[ i ];
if ( value[ i ] >= TEN_POW4 )
{
value[ i ] -= TEN_POW4;
++value[ i + 1 ];
}
s += table[ value[ i ] ];
}
if ( 0 == mark[ s ] )
{
mark[ s ] = k;
}
if ( j == ++k )
{
if ( 0 == k )
{
break;
}
j *= ( 5 == ++m ) ? 0 : TEN_POW2;
}
delta[ 0 ] += 2;
i = 0;
while ( delta[ i ] >= TEN_POW4 )
{
delta[ i ] -= TEN_POW4;
++delta[ ++i ];
}
}
HugeCalc::EnableTimer( FALSE );
s = 0;
i = 0;
while ( s < 9*4*5 )
{
s += cycleDelta[ i ];
if ( 4 == ++i )
{
i = 0;
}
if ( 0 != mark[ s ] )
{
if ( mark[ s ] <= 0xFFFF )
{
printf( "X = %u\tS = %u\n", s, mark[ s ] * mark[ s ] );
}
else
{
printf( "X = %u\tS = %I64u\n", s, UInt32x32To64( mark[ s ], mark[ s ] ));
}
}
else
{
printf( "X = %u\n", s );
}
}
printf( HugeCalc::GetTimerStr( FT_DOT06SEC_s ) );
printf( "\n" );
return 0;
}备注:代码中,仅调用了 HugeCalc 的高精度计时部分。
运行结果
X = 1 S = 1X = 4 S = 4
X = 7 S = 16
X = 9 S = 9
X = 10S = 64
X = 13S = 49
X = 16S = 169
X = 18S = 576
X = 19S = 289
X = 22S = 1849
X = 25S = 4489
X = 27S = 3969
X = 28S = 17956
X = 31S = 6889
X = 34S = 27889
X = 36S = 69696
X = 37S = 98596
X = 40S = 97969
X = 43S = 499849
X = 45S = 1887876
X = 46S = 698896
X = 49S = 2778889
X = 52S = 4999696
X = 54S = 9696996
X = 55S = 19998784
X = 58S = 46689889
X = 61S = 66699889
X = 63S = 79869969
X = 64S = 277788889
X = 67S = 478996996
X = 70S = 876988996
X = 72S = 3679999569
X = 73S = 1749999889
X = 76S = 5599977889
X = 79S = 7998976969
X = 81S = 8998988769
X = 82S = 17999978896
X = 85S = 36799899889
X = 88S = 88998998929
X = 90S = 297889998849
X = 91S = 299879997769
X = 94S = 897977978689
X = 97S = 975979998889
X = 99S = 3957779999889
X = 100 S = 2699997789889
X = 103 S = 9879498789889
X = 106 S = 9998768898889
X = 108 S = 29998985899689
X = 109 S = 85986989688889
X = 112 S = 97888999968769
X = 115 S = 386999898769969
X = 117 S = 429998989997889
X = 118 S = 578889999977689
X = 121 S = 898999897988929
X = 124 S = 1959999889996996
X = 126 S = 6788999798879769
X = 127 S = 3699998989898689
X = 130 S = 9895699989899689
X = 133 S = 38896878989988889
X = 135 S = 67699789959899889
X = 136 S = 38999699989995889
X = 139 S = 188997899869998769
X = 142 S = 279869897899999969
X = 144 S = 498999778899898896
X = 145 S = 1877896979979898969
X = 148 S = 989879999979599689
X = 151 S = 6979497898999879969
X = 153 S = 5899989587897999889
X = 154 S = 8899988895999696889
X = 157
X = 160
X = 162
X = 163
X = 166
X = 169
X = 171
X = 172
X = 175
X = 178
X = 180
103.824245 s
Press any key to continue
相对于 37#,不仅把 X=145 这个空洞给填补了,还多了更大的三组数据。
发过帖后,才发现本算法已实现了 无心人 在 47# 中提到的 10000 进制进行优化。
注意,我的算法中完全没有64bit长整型,不含除法,连乘法都极少。