| 
注册时间2007-12-27最后登录1970-1-1威望 星金币 枚贡献 分经验 点鲜花 朵魅力 点上传 次下载 次积分48888在线时间 小时 
 | 
 
 发表于 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 | 太强悍了 |  
查看全部评分
 |