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

[擂台] [转载]整数分解为2的幂的和

[复制链接]
发表于 2010-4-1 09:48:27 | 显示全部楼层
找到一参考链接,希望对大家有帮助:http://www.research.att.com/~njas/sequences/A002577
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2010-4-1 10:00:01 | 显示全部楼层
参考结果:
(09:58) gp > fun(10^10)
%42 = 54049175382734938696472608138751495789943756108317583013021460100933300805
7412691704996494093758954301613502661738339294616171428
(09:58) gp > fun(9999999999)
%43 = 54049175081216462692786814372794682219839080412323933520822604994459632936
5195818680963578081187969378596524625340118423702393914
(09:59) gp > fun(5000000000)
%44 = 30151847600368579376595681357010467569599364949219885510647366786922168730
24032916012570984923016978036398220870913777514
(09:59) gp > fun(10^20)
%45 = 43042563299313952734013069432214177548042434846399288362760684201324989381
19428155153121889280159376379011783329890393930499395722101736097081625686224245
66366498980394146270463052257252966304363501172774521050318507344828008882298390
82379847955523132491831150650058309077548396707331549148234717834318442932294597
34780853201507119441850242992445445817195949632660066541654548725984655426816444
32289838453044740887422260486644915610243351008926612883770007880770769723370148
31687732527977588782742511987239798196320278659104637075388654703281904856603609
51657838864444
(09:59) gp > fun(10^30)
%46 = 39606084812969213925551982340151089705142237904741895617156498724063177667
21320341876595581006723343069015458259478849604917703362220601821556551674309071
10828619568537552685942227204940134428356222759829765518170238870314147611278236
85494523721691763316758488454598459961114530156082623316213258082188390899358233
72366252515528858296790857652494807925453554349845754320856168299780608051609926
13912228388330678580301180800512316075280178526424996161715000541127476693494567
90382411541215886828108211479735297373035902774857439924701078778055954751999050
40968676452305585927594245021545034369759677922652619318017171782845499012871911
87028434441725326473690698895110662621263236239885101370883988957059225135578740
69358636456792234062778343807699071369813058497552167267109805166395202370147096
14296017491335856619518297428998486688551358778330293627349883576997523892906181
84644230104826242306366421162126227084760048563582961567808404817513925828159827
62468612278282564430097288659530822485252701065924658725956424537058031728240707
01735948861891791382577589823526647095294846813850862730940205503140119846350116
59079838548998674775996761572814016168120026321896415137607685922915924382951526
85397702061554260601110450776855230852618256512050198991855176738551767164723643
528912674547290921860946750646161732695717302956498566820
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2010-4-1 10:06:36 | 显示全部楼层
对应Pari/Gp代码
  1. bitcount(n)=
  2. {
  3.     local(r);
  4.     r=if(n<=0,0,1+bitcount(n>>1));
  5.     r
  6. }

  7. bitmap(n)=
  8. {
  9.     local(v,c);
  10.     v=vector(0);
  11.     c=0;if(n%2==1,n=n-1);
  12.     while(n>0, if(n%2!=0,v=concat(v,c));c=c+1;n=floor(n/2));
  13.     v
  14. }
  15. mapm(b)=
  16. {
  17.    local(T);
  18.    T=matrix(b-1,b-1);
  19.    T[1,1]=1;
  20.    for(u=2,b-1,
  21.       T[u,u]=1;
  22.       for(v=1,u-1,
  23.          h=T[u-1,v];
  24.          for(w=1,v-1,
  25.             h=h+T[u-1,w]*T[u-1-w,v-w]
  26.          );
  27.          T[u,v]=h+T[u-1,v];
  28.       )
  29.     );
  30.     T
  31. }

  32. fun1(n)=
  33. {
  34.     local(b,T,h,R,m);
  35.     b=bitcount(n);
  36.     T=mapm(b);
  37.     m=bitmap(n);
  38.     R=T[m[1],];
  39.     for(u=2,length(m),
  40.        for(v=0,m[u]-1,
  41.           h=T[m[u],m[u]-v];
  42.           for(w=1,m[u]-v-1,
  43.              h=h+R[w]*T[m[u]-w,m[u]-v-w]
  44.           );
  45.           R[m[u]-v]=R[m[u]-v]+h
  46.        )        
  47.     );
  48.     h=0;
  49.     for(u=1,b-1,h=h+R[u]);
  50.     h+1
  51. }

  52. fun(n)=
  53. {
  54.     if(n<=1,1,fun1(n))
  55. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2010-4-1 10:21:52 | 显示全部楼层
关于这个问题我现在找到两个时间复杂度为$O(log^3(n))$的算法。
上面贴出的这个方法实现起来比较简单些,但是消耗的内存要大一些。

评分

参与人数 1鲜花 +8 收起 理由
wayne + 8

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2010-4-1 10:38:19 | 显示全部楼层
这种方法直接使用这个模型,而没有使用前面的递推式.
对于整数n,我们可以将它分解成二进制形式
$n=2^{i_1}+2^{i_2}+...+2^{i_t};0<=i_1<i_2<...<i_t$
这个分解对应于n的一个划分,我们称为标准划分
我们将一条长度为n的线段按上面方法划分成t段,长度依次为$2^{i_k}$,其中左边的子线段长度总是小于右边的。
同样对于n的任意一个划分,我们也同样可以将长度为n的线段划分成若干个长度为2的幂的子线段,而其左边子线段的长度总是不大于右边的子线段。
对于n的一个任意划分,如果我们将它的图案同标准划分的图案比较,我们会发现,所有标准划分的分界必然是任意划分的边界,如图:
g1.GIF
如果我们已经知道对于整数$2^t$,将其划分成2的幂,其中最大幂为$2^h$的不同划分数目为T(t,h)
那么对于整数n,我们可以通过动态规划来计算
记其中将$2^{i_1}+2^{i_2}+...+2^{i_k}$划分成2的幂,而且使用的最大幂为$2^h$的划分数目为R(k,h)
那么我们可以得到
$R(k+1,h)=sum_{u=0}^{i_k}R(k,u)*T(i_k-u,h-u)$
然后最后累加计算$\sum_hR(t,h)$就可以了
15#就是对应的Pari/Gp代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2010-4-1 10:40:19 | 显示全部楼层
效率很高啊

(10:39) gp > fun(10^50)
%24 = 88863001401273174864626472666793341626417026891733807347542089314849881492
05801173987246954080765467213743791365013997329999025139675552015911965213784053
79141842866987988178682899264337190104350090364254052063533650381484017787209412
93223311849395737463936753873458960430637734530634238590902843090810094273528506
83639450768730950848617344562622338345412650523186091529046901354820421416097666
53280597368964391921158792458600892166405397818698377972913269407949665415258011
72746616481058052372111233384089597714828642051266768443878139916009676594260401
99957430620974686395559469290473472260360523193558808822879960049512632931608206
97684197273647351236972001893520131038080958663345065889852897707962779139628399
16154809120127766656205976553198254345261560905684202909384094842669214209202592
95150949611207820621288675657661390446369007847539582928573989132025106669988250
56239044305899297508682701337788038772991681078121842729446779488912581108202828
00232767097517412554440678220050169461427776383981443800281539495885368046288134
22182975617805624647592897892747530215994731583772035156379580509117249276181734
10681063571780970329818845281282176042609061438250222462602230013860273036065974
99258848282158583346452784226914339480942171280409545970986667536119745475644133
27457666155265650205054184019866995981423602516743389305472986434972692938322337
86278616371572686162392304278129728048865745385853941415456446394140899479904335
60791293197836097890696406916644842409328482374566170528301373040669281623244963
62667935189923887177348339482956576377806811899569967550967812604369716634894472
32092318949503471483020878426127285498549441988436086967787967849223419687514754
44085822687972678299547955707012123110870254946909086039183672236350532249321328
46700849459167493506906011115032109749251336742743428509373074089425684369744412
57906165833050332394136197406880235958042728504828042204028304249930048908215003
42699269519756561727221981701152724294801645378633161952899565207056189890147114
13995740966464575976769740720933626283493031768269395926959495714438942701385487
92818681917302347723325710605549361104122340075217474359725637805262472615963639
04327480211717831959882429333096193081553094722734957220826418165866622607749391
23819416520235056141177141859875325007448606523150498012656869042037220308877277
30733980529933686438419380978586990015115934550158375515353646497416476233161306
08595131799145388011136515340497233625518043229000833046579029251080361520353402
42534161052586878921973858740809952004749226442519128939524712095878325932100707
37993350541100974601931138512557159028882168756839210615586122086587715245429419
26935715208577694412028516948889826582868577669843095661416750761837147414371532
40438873534736633332367084849612979242935741194148313349901637126323915510678956
03510270675925545118778709686882529969065079762628756901521868282350719088905352
45955489764736818650303445819550694491787718506397991703506902253539660292751578
89520714934530561477632120006643679500686589299493560361542015664880653898965262
03989048305115429495145032655196191753124185695782083170118201656393039552676351
10169549909210371986120260430236693478294767463898259146659530793298519915592561
91357480695167831765549423591775604149094168893322306060522247581894915279555389
68103544979470559893115616729615929307272335879925884729843365474946434832308592
22297423181105788228796806607900143982449806593639698132668658504979421119799376
19228099922417406999357339933520864015535671723826209881138001009302917193402822
91141705297460605200666754722272538671850945940022403020544338679607743198649641
80217099925171794779646162406975761477790110112834689184402095945401464487216147
20352927745864406933839509624623482898277664762931622413390526026497560351038295
13619613110693665893847369743836080882710470301235633610650969673873553834726102
54556
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2010-4-1 11:12:12 | 显示全部楼层
另外一种方法我认为更加有效一些,因为使用的内存会更小,可以用于计算更加大的范围的结果,但是实现起来更加麻烦一些,而且我不能像Fans一样,
可以直接将编程工作交给“秘书”,这里就纸上谈兵一把:
我们记$g(n)=f(2n),g(0)=1$,那么我们可以知道
$g(2n)=2g(0)+2g(1)+...+2g(n-1)+g(n)$
$g(2n+1)=2g(0)+2g(1)+...+2g(n-1)+2g(n)$
我们记$n_h=n>>h$,也就是$n_h$是n右移h位,特别的$n_0=n$.而设n的二进制表示为b位,所以$n_b=0$
所以我们知道,如果$n_0$是奇数,那么必然有
$g(n)=2g(0)+2g(1)+...+2g(n_1-1)+2g(n_1)$
而如果$n_0$是偶数,那么必然有
$g(n)=2g(0)+2g(1)+...+2g(n_1-1)+1g(n_1)$
更加一般的,我们假设
$g(n)=u_{n,h}(0)g(0)+u_{n,h}(1)g(1)+...+u_{n,h}(n_{h+1}-1)g(n_{h+1}-1)+c_{n,h}g(n_{h+1})$
其中$u_{n,h}(x)$是一个h次多项式,
我们将右边再次使用上面的基本递推公式得到
$g(n)=$
$u_{n,h}(0)*1*g(0)+$
$u_{n,h}(1)*2*g(0)+$
$u_{n,h}(2)*2*g(0)+u_{n,h}(2)*1*g(1)+$
$...$
$u_{n,h}(n_{h+1}-1)*2*g(0)+....$
$c_{n,h}*2*g(0)+...$
我们现在查看右边$g(0),g(1),...,g(n_{h+2})$的系数可以知道
它们的通常形式是$g(t)$的系数为$u_{n,h}(2t)*1+2*sum_{u=2t+1}^{n_{h+1}-1}u_{n,h}(u)+2*c_{n,h}$
而只有$n_{h+1}$是偶数时,最后一项会例外,它变成$c_{n,h}$
而上面通常形式中第一项是t的h次多项式,那个求和式可以变成t的h+1次多项式,所以总体会变成一个h+1次多项式
也就是说,同前面类似,除了最后一项系数可能例外,其余的系数可以用一个h+1次多项式统一描绘。
所以如果我们一直通过这种方法递推下去,找到$u_{n,b-1}$和$c_{n,b-1}$,那么就可以计算出$g(n)$了。
但是这里还有一个问题,就是里面那个多项式累加求和如何计算的问题,我们当然可以直接使用$sum_{x=t}^s x^h$的求和公式,但是比较复杂
而且会出现分数表示。一个更加好的办法是记
$v_0(x)=1,v_1(x)=x/{1!},v_2(x)={x(x-1)}/{2!},...,v_t(x)={x(x-1)...(x-t+1)}/{t!},...$
然后我们需要事先计算表格M(a,b)使得
$v_s(2x)=sum_{t=0}^s M(s,t)v_t(x)$
这个表格的计算也比较简单,首先,$v_0(2x)=v_0(x)$
其次
$v_{s+1}(2x)=v_s(2x)*{2x-s}/{s+1}=sum_{t=0}^s M(s,t)v_t(x){2x-2t+2t-s}/{s+1}=sum_{t=0}^s ({2(t+1)}/{s+1}M(s,t)v_{t+1}(x)+{2t-s}/{s+1}M(s,t)v_t(x))$
我们可以递推出整个表格M,只是每次计算时会牵涉到一个分母s+1,但是我们知道合并结果以后必然是整数,所以可以计算过程都乘上s+1,最后除去s+1就可以了。
先我们可以设$u_{n,h}(x)=sum_{t=0}^h U_{n,t}*v_t(x)$
我们现在主要需要计算$sum_{u=2t+1}^{n_{h+1}-1}u_{n,h}(u)$
由于$sum_{u=2t+1}^av_s(u)=v_{s+1}(a+1)-v_{s+1}(2t+1)$
所以合并以后每个$g(t)$的系数形式为$sum_s a_s*v_s(2t)+b_s*v_s(2t+1)+c$
这时我们就需要表格M(a,b)将$v_s(2t)$转化为$v_s(t)$的线性组合。而$v_s(2t+1)=v_{s-1}(2t)+v_s(2t)$,于是它也可以通过同一个表格转化
由此我们可以得到$u_{n,h+1}$
这个算法整体时间复杂度也是$O(log^3(n))$此大数乘法运算,内存复杂度看上去同前面算法相同,也是$O(log^2(n))$,但是前面算法中数组T保存的是$2^t$的划分数,是很大的整数,
而这里的表格M中保存的数相对都比较小(不会超过n),所以花费的内存要小很多
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2010-4-1 13:16:35 | 显示全部楼层
15#忘了介绍T(u,h)的计算方法了。
如果h==u,那么显然T(u,h)=1,如果h<u,我们将$2^u$长度的线段等分成两段,类似15#中n的划分方案,对于任意一个方案,这两段都会包含完整的子线段。于是我们根据前面一个子线段中最大线段长度,可以将划分数目分解为:$T(u,h)=\sum_{t=0}^hT(u-1,t)*T(u-1-t,h-t)$
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2010-4-1 13:36:20 | 显示全部楼层
找到一参考链接,希望对大家有帮助:

准确的来说,是这个链接:

http://www.research.att.com/~njas/sequences/A000123
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2010-4-1 13:45:03 | 显示全部楼层
我试图找更贴切的,但没有找到。
还是LS的厉害!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2024-5-6 14:45 , Processed in 0.049385 second(s), 18 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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