- 注册时间
- 2007-12-27
- 最后登录
- 1970-1-1
- 威望
- 星
- 金币
- 枚
- 贡献
- 分
- 经验
- 点
- 鲜花
- 朵
- 魅力
- 点
- 上传
- 次
- 下载
- 次
- 积分
- 48163
- 在线时间
- 小时
|
发表于 2025-9-18 13:26:03
|
显示全部楼层
使用上面算法刚刚弄出了一个pair/gp代码,
但是代码还有点小bug,也就是如果输入数据正好是底数的幂或者比幂小1,再高精度也算不出。
比如算5^k以24开头,由于24=5*5-1,计算会报错。
调用格式为
findappr(base, target, digits, MOD)
其中MOD需要使用一个很大的整数,然后后面计算都会采用模MOD计算,MOD太小会导致计算精度不够报错。比如可以选择MOD=10^100,
而代码中会有一小部分前置浮点计算,需要有足够精度,所以使用前需要设置足够高浮点精度,比如MOD使用10^100,那么计算精度要至少略高于100位,比如可以使用
\p 150
表示使用150度浮点计算精度。
这个函数会返回最小的整数k,使得base^k的前digits位和target前digits位相同。
- testinrange(s,s2,U,V,M)=
- {
- if(U<V,
- if(s>=U+1&&s<=s2&&s2<=V,1,
- if(s2<=U&&(s<=s2||s>=V+1)||s>=V+1&&(s<=s2||s2<=U),0,-1)),
- if(s>=U+1&&(s2<=V||s2>=s)||s2<=V&&(s<=s2||s>=U+1),1,
- if(s<=U&&(s2>=U||s2<s) || s2>=V&&(s<=V||s>=U),0,-1))
- )
- }
- dist(s,s2,C,M)=
- {
- local(d,s3,e);
- s3=(s+s2)/2; if(s2<s,s3+=M/2);
- e=(s2-s)/2;if(s2<s,e+=M/2);
- d=abs(s3-C); if(d>=M/2,d=M-d);
- d+e
- }
- sgndist(s,s2,C,M)=
- {
- local(d,s3,e);
- s3=(s+s2)/2; if(s2<s,s3+=M/2);
- e=(s2-s)/2;if(s2<s,e+=M/2);
- d=s3-C; if(d>=M/2,d-=M);
- if(d<=-M/2,d+=M);
- if(d<0,d-e,d+e)
- }
- findk2(U,M,t,s)=
- {
- local(lu,ld,ru,rd,mu,md,tL,tD,r1,r2,pass);
- lu=0;ld=1;ru=1;rd=1;
- mu=lu+ru;md=ld+rd;
- tL=if(s<0,M-t,0);
- tD=if(s<0,M,t);
- while(1,
- r1=(md*U)%M;r2=(md*(U+1))%M;
- pass=testinrange(r1,r2,tL,tD,M);
- if(pass!=0,break());
- if(mu*M<md*U,ld=md;lu=mu, if(mu*M>=(U+1)*md,rd=md;ru=mu,pass=-1;break()));
- mu=lu+ru;md=rd+ld;
- );
- if(pass<0, md=0);
- md
- }
- findmink(U,V,R,M)=
- {
- local(t,C,s,s2,pass,d,t2,dt,dt2);
- C=(V+R)/2;
- t=2;pass=0;
- while(1,
- s=(t*U)%M;
- s2=(t*(U+1))%M;
- pass=testinrange(s,s2,V,R,M);
- if(pass!=0,break());
- d=dist(s,s2,C,M);
- if(d<M/8-1,break());
- t=t+1
- );
- dt2=0;
- while(pass==0,
- s=(t*U)%M;s2=(t*(U+1))%M;
- pass=testinrange(s,s2,V,R,M);
- if(pass!=0,break());
- d=sgndist(s,s2, C, M);
- t2=if(d>0, findk2(U,M,2*d,-1), findk2(U,M,-2*d,1));
- if(t2==0, print("Out of precision");pass=-1;break());
- dt=if(d>0, M-(t2*U)%M, (t2*(U+1))%M);
- dt2=if(d>0, floor((d-R+C)/dt),floor((-d+V-C)/dt));
- if(dt2<=0,dt2=1);
- t=t+dt2*t2;
- );
- if(pass==-1, print("Precision not enough, need larger M");0,
- t)
- }
- findappr(u,v,n,M)=
- {
- local(lu, lv, rv, t, U, V,R);
- lu=log(u)/log(10);
- lu-=floor(lu);
- lv=log(v)/log(10);
- t=floor(lv);
- v=floor(v*10^(n-1-t));
- lv=log(v)/log(10)-(n-1);
- rv = log(v+1)/log(10)-(n-1);
- U=floor(lu*M);
- V=floor(lv*M);
- R=floor(rv*M);
- findmink(U,V,R,M)
- }
复制代码
比如使用
\p 1000
for(u=5,100, s=u!;r=findappr(18, s, length(digits(s)), 10^500);print(u,",",r))
得到
5,63
6,1680
7,3642
8,14387
9,1696738
10,3998016
11,10947946
12,825514266
13,1821688762
14,197680918419
15,38981198865221
16,15023912390981
17,258918396056895
18,18553374137359214
19,149897101165929134
20,5337998514529961988
21,7957523787372878900
22,48885793990344288895
23,44142110555196810262850
24,578605659839921996371529
25,112783699177135451463259207
26,520754961861071375905926338
27,327864967267938420194303002269
28,137464934434774447528980538258
29,32519371460891860088818610758697
30,360763921805840529622482428966766
31,20399045306454851654826214649915206
32,1395130806326082837353193042296705650
33,1064638931316749760463561942911005277
34,1181657051954195954103216667655460729949
35,2190509075652518134534062095972756779758
36,354905161799969822808667538624939986888082
37,6221940206327830814751851889041152222871268
38,533838105191500155427303371577710575524399974
39,19047773280283685209020700170999048776339446246
40,2206338681742557865698306229856813923814833133780
41,740459273205894939640204120318661175921839689068662
42,2617068225891253087772769203038826284017857620091974
43,45093668414994481240471655185424120348236541757053456
44,2228074893681735129913146374325417488510975394162385176
45,2483565816826088747586892429548905862844976212722004781487
46,9543446793427689131150448190840824011346487795502921412743
47,787612393061307296052807988771781403034104422792132872300041
48,5740567549470619851258056631738062761148307204320654912761887
49,34661582247726859084814384362756574521636589109545298480753594
50,52170379843798184322766971281459512181129341689265295073044042359
51,4022057962132918240344054205288062180490393260439457462217680028009
52,44653284601058096160984785999423552771560818435411682769197307997155
53,4828189017124606116539571494688283230887687456761079401045806529144021
54,307002273867246981574968668480440961117680960416092740936067165010197831
55,5414075147786470809279420576334489995464619956425226222060492681823827039
56,1442230095340767886589041773607378383564517369148963196857585087185910109473
57,17797521219970685899750161058380736023235356579527335706359792369317311260919
58,9758968532745756119153517117502849160860317723107621325115970095657211568476106
59,712958597667451819156065551902607864662439399091027474715518998088714619149761763
60,53910589585359494086614544509902881978604850627952502573908038041823774354388121278
61,443000650692405870948190021958944226140448417448771731286615886823922168344460403277
62,40032166125381085209094078302970040195332847588092497067795628486128022027580773840292
63,10530369284983441260729679487822890757545907873895261889246609734155294752356063392740246
64,249180871957654731718506091741365827336362356020266210130615408917069371945497725773013753
65,2042290384357301242673044785192649442664290240096571703146569007144820692516202745983442280
66,2327404037988383368865409837980472081872947936330541753109541588627147009684035106777623023766
67,68431695519980734068154069027582980371920117753223643528209729141198797714291437412835087830146
68,8214151422728267235745545424161455011130292332926308834995015364539861497641947052703917405765741
69,104248553808740001097835222540744404410923075914212457887416832860788613954661928129826580233667847
70,3585198975382389799054693641162767056212371186035342945050902915173455818422622424811275360783777298
71,1323152352166533317560312955388886229240778344827991287499119514257177992404751625369852500473404454395
72,88417931980305887898400316774467547613640058968253357225342854789284884116263458248301676068786528229413
73,13683012506857002212063306337117729545358176772429135768448366663823858089110873102842753069715427318736983
74,32137917832999726935326841685387144045071975797938686043874716308993136939254916126935386913179954530476223
75,4374659841635175330367741933866944510404383372238198821287467838939007070456473005899169519368181050050405301
76,5937662002552857686021927085332125322275663144632588109976232839073687531689054702432083906040236294096287453903
77,121172222860532316962162865143200877053113843413904030935208908544249884906641401642807911390820851213799918525072
78,97248778733239690723159158376615079517366686059000553893695592886013663955971468090684553890225038683323948839874998
79,4447269259571380209794996557692844198979354573319469845771393290590741415452907585780548760727903483026858700625323868
80,86665995599721937049816666151044977776315566481720036246611387998534277779062465520374008650134141963813665582167815067
81,58193062395149795845715912553195151247374474129885240511229425784189608011302826161510000814779493822767226779817936804653
82,656671109657534638806325175521880315227564629407223615210670215571475571815420583126313830713602093341059329487705437076219
83,36516856169848310286682743527492075523725783974453817378841108134843472098969341824268600916978744744122235542476362654287506
84,16927382189046029646896532074071547135274179991067829685826632106167638388263959227986201025616157390295952501376464613095479036
85,83764139052087301285798285544474938687593865968665947524867349070835071785891002127483285590845886089114530731285830828886178393
86,43576396000362044816138671117964851379733989262897646483436911694862704815097637167007961948514526460555223255173819427762592755455
87,3644835550555795678900612873812501434739557986916738342160324377985875798499507575125487894814976430545670221702563898015753891586667
88,90040612337868602431764199610995212153297867828968775808974277258670057798174376336607697505391816032421862353063751147217192953346215
89,12969180943355325774828350298494384593201167512243444167341160800780118461064934898112137406102370900524502179852016284255629786399380174
90,1584845169218190463155746987199256235048768451788098873719084968425327000085600317670867584360829758481252247214136321411675625064218185562
91,379742512898420158507821638817141315039483409261291152717973380264550730464221141677937860532350260486439371044620795379171325600392659442236
92,18164310105427713000149354461060225889020899426423458119794029912199162382266408572018701593505471246065596486531277634541018937603824971821015
93,3193393236342603624613023695254161839651553017437888175593112866096963625918653075278791509089150343455792608903261211335307391867499904017428965
94,92475329450142826051366437267821699714166409997074849944607827301578152944194472408133787793662642267061643416042525503081123185304486141185815593
95,40408540662145922816967444482183318237414821791273956884749966647107557500172082406257513942943817554149412334630610164270812162191938440940554385066
96,885692151722805825889899968877889263519281541533369364842904839127053587180604799143305523695409238624115475056399001232903764930443537644980667343169
97,115700783207673269482459568194675841659797139016133115365966016270119881404391438034065380949535235789104793143281984816893981935044463991124654754156929
98,90330590414660647065397551980237659174887334302608803887479927424258646831323795701785878192011878120573402871099721290620342184510448284672735686959037
99,11337826467067076076629219480315647406573250353843966131702069365700797672831705805416672607396602145085179910242885910812066433782653931259297912219077423427
100,205434530482046557360596524990087979830518401961820077315640155576171813279107446192065776079720231337742869534959018105639656814524453726449505716744816269258
比如还可以
\p 1000
findappr(Pi,1/9,100,10^500)
得到
\(\pi^{3570227448978643820097344439027815950348875034718098681378739266626694283136290986408344492598865460}\)的前100位都是1 |
评分
-
参与人数 2 | 威望 +16 |
金币 +16 |
贡献 +16 |
经验 +16 |
鲜花 +16 |
收起
理由
|
王守恩
| + 8 |
+ 8 |
+ 8 |
+ 8 |
+ 8 |
一骑绝尘! |
northwolves
| + 8 |
+ 8 |
+ 8 |
+ 8 |
+ 8 |
太强悍了 |
查看全部评分
|