找回密码
 欢迎注册
楼主: mathe

[擂台] 平方数数字和

[复制链接]
发表于 2008-7-10 11:21:54 | 显示全部楼层
看mathe的方法,对于很多的x,应该都可以构造01向量来达到。
那么,接下来的困难得多的问题是,对于小于N的自然数x=a的大概有多少呢?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-7-10 11:44:52 | 显示全部楼层
这个问题比较难。可靠的方法很难想到。
不过我们可以假设平方数各位数都是均匀分布的,给出一个近似统计数据来(计算公式估计还是有点麻烦)。
而更加精确的模型可以对平方数的首位和末位的分布各个更加精确一些的统计模型,比如末位应该是$1/5$的概率取1,4,9,6,然后$1/10$的概率取0和5
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-7-10 13:10:22 | 显示全部楼层
刚才把 36# 算法再次优化,
完全消除了除法及64位的运算,
算出 X=130,仅需 11.243 s;遍历18位以内的平方数,仅需 126.074 s

大家看看还有哪里可以优化的?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2008-7-10 14:07:46 | 显示全部楼层
函数SumOfDigits可以提速的,比如我们事先计算1~9999所有数字的sumofdigit,并且保存
然后对于每个整数,每四位的对应和再累加,可以减少一些乘除次数.
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-7-10 14:16:21 | 显示全部楼层
类似想法我在 38# 曾提到过。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-7-10 19:40:28 | 显示全部楼层
我想

折半算法更好吧
由于这个求数位和是并行操作
所以考虑如果采取精细的汇编
假设考虑到32位十进制数字
则首先除以10^16
分开,分开的部分
分别再除以10^8
则只需要
一次除以10^16
二次除以10^8
四次除以10^4
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-7-10 19:41:58 | 显示全部楼层


更黑的方法在下面
如果一开始就存储
1-9999内的乘法
共5*10^7个结果
则可以考虑整体的数字
采取10000进制
呵呵
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-7-10 19:42:51 | 显示全部楼层

但上述方法只能对计算超过2^64以上的
数字的平方和才有好的效果
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-7-10 20:11:07 | 显示全部楼层

更加优化的算法

下面的代码可以遍历全部的64bit内的完全平方数,且效率非常高:
  1. #include <stdlib.h>
  2. #include <stdio.h>

  3. #include "../../../HugeCalc_API/CppAPI/Include/HugeCalc.h"    // 公共接口

  4. #pragma message( "automatic link to ../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )
  5. #pragma comment( lib, "../../../HugeCalc_API/CppAPI/Lib/HugeCalc.lib" )


  6. #define TEN_POW2    100UL
  7. #define TEN_POW4    10000UL

  8. int main( void )
  9. {
  10.     const UINT32 cycleDelta[ 4 ] = { 1, 3, 3, 2  };    // 0, 1, 4, 7
  11.     UINT32 table[ TEN_POW4 ];
  12.     UINT32 mark[ 9*4*5 + 1 ];
  13.     UINT32 value[ 5 ], delta[ 5 ];
  14.     UINT32 * p = table;

  15.     UINT32 i, j, k, m, s;
  16.     UINT32 i2, i3, i4, s2, s3;

  17.     HugeCalc::EnableTimer();
  18.     HugeCalc::ResetTimer();

  19.     for ( i4 = 0; i4 < 10; ++i4 )
  20.     {
  21.         for ( i3 = 0; i3 < 10; ++i3 )
  22.         {
  23.             s3 = i4 + i3;
  24.             for ( i2 = 0; i2 < 10; ++i2 )
  25.             {
  26.                 s2 = s3 + i2;
  27.                 for ( i = 0; i < 10; ++i )
  28.                 {
  29.                     *p++ = s2 + i;
  30.                 }
  31.             }
  32.         }
  33.     }

  34.     memset( mark, 0, sizeof( mark ));
  35.     memset( value, 0, sizeof( value ));
  36.     memset( delta, 0, sizeof( delta ));

  37.     k = 1;
  38.     j = TEN_POW2;
  39.     value[ 0 ] = 0;
  40.     delta[ 0 ] = 1;
  41.     m = 1;
  42.     for ( ; ; )
  43.     {
  44.         s = 0;
  45.         for ( i = 0; i < m; ++i )
  46.         {
  47.             value[ i ] += delta[ i ];
  48.             if ( value[ i ] >= TEN_POW4 )
  49.             {
  50.                 value[ i ] -= TEN_POW4;
  51.                 ++value[ i + 1 ];
  52.             }

  53.             s += table[ value[ i ] ];
  54.         }
  55.         if ( 0 == mark[ s ] )
  56.         {
  57.             mark[ s ] = k;
  58.         }


  59.         if ( j == ++k )
  60.         {
  61.             if ( 0 == k )
  62.             {
  63.                 break;
  64.             }

  65.             j *= ( 5 == ++m ) ? 0 : TEN_POW2;
  66.         }

  67.         delta[ 0 ] += 2;

  68.         i = 0;
  69.         while ( delta[ i ] >= TEN_POW4 )
  70.         {
  71.             delta[ i ] -= TEN_POW4;
  72.             ++delta[ ++i ];
  73.         }
  74.     }

  75.     HugeCalc::EnableTimer( FALSE );


  76.     s = 0;
  77.     i = 0;
  78.     while ( s < 9*4*5 )
  79.     {
  80.         s += cycleDelta[ i ];
  81.         if ( 4 == ++i )
  82.         {
  83.             i = 0;
  84.         }

  85.         if ( 0 != mark[ s ] )
  86.         {
  87.             if ( mark[ s ] <= 0xFFFF )
  88.             {
  89.                 printf( "X = %u\tS = %u\n", s, mark[ s ] * mark[ s ] );
  90.             }
  91.             else
  92.             {
  93.                 printf( "X = %u\tS = %I64u\n", s, UInt32x32To64( mark[ s ], mark[ s ] ));
  94.             }
  95.         }
  96.         else
  97.         {
  98.             printf( "X = %u\n", s );
  99.         }
  100.     }

  101.     printf( HugeCalc::GetTimerStr( FT_DOT06SEC_s ) );
  102.     printf( "\n" );

  103.     return 0;
  104. }
复制代码
备注:代码中,仅调用了 HugeCalc 的高精度计时部分。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-7-10 20:11:44 | 显示全部楼层

运行结果

X = 1   S = 1
X = 4   S = 4
X = 7   S = 16
X = 9   S = 9
X = 10  S = 64
X = 13  S = 49
X = 16  S = 169
X = 18  S = 576
X = 19  S = 289
X = 22  S = 1849
X = 25  S = 4489
X = 27  S = 3969
X = 28  S = 17956
X = 31  S = 6889
X = 34  S = 27889
X = 36  S = 69696
X = 37  S = 98596
X = 40  S = 97969
X = 43  S = 499849
X = 45  S = 1887876
X = 46  S = 698896
X = 49  S = 2778889
X = 52  S = 4999696
X = 54  S = 9696996
X = 55  S = 19998784
X = 58  S = 46689889
X = 61  S = 66699889
X = 63  S = 79869969
X = 64  S = 277788889
X = 67  S = 478996996
X = 70  S = 876988996
X = 72  S = 3679999569
X = 73  S = 1749999889
X = 76  S = 5599977889
X = 79  S = 7998976969
X = 81  S = 8998988769
X = 82  S = 17999978896
X = 85  S = 36799899889
X = 88  S = 88998998929
X = 90  S = 297889998849
X = 91  S = 299879997769
X = 94  S = 897977978689
X = 97  S = 975979998889
X = 99  S = 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长整型,不含除法,连乘法都极少。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-3-29 23:12 , Processed in 0.058479 second(s), 15 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表