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

[求助] 求最小正整数 k,使得 5^k 起始9位数为12!=479001600

[复制链接]
 楼主| 发表于 4 天前 | 显示全部楼层
这个最终提交了前15项,今天早上刚正式发布:

A386547
a(n) is the least integer k such that 5^k begins with n!.
+30
0
0, 0, 2, 4, 12, 116, 2080, 6017, 149704, 2436360, 10819405, 19295517, 1664457026, 55459990176, 223913486556, 700785493886
(list; graph; refs; listen; history; text; internal format)
OFFSET
0,3
LINKS
Table of n, a(n) for n=0..15.
FORMULA
a(n) = A018862(n!).
EXAMPLE
a(4) = 12 because 5^12 = 244140625 is the smallest power of 5 beginning with 4! = 24.
MATHEMATICA
a[n_] := Module[{target = IntegerDigits[n!], k = 0}, While[UnsameQ[Take[IntegerDigits[5^k], Length@target], target], k++]; k];
Table[a[n], {n, 0, 8}]
CROSSREFS
Cf. A000142, A018862, A374922, A374923, A387304.
KEYWORD
nonn,base,more,new
AUTHOR
Zhining Yang, Jul 25 2025
STATUS
approved
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 4 天前 | 显示全部楼层
{2,3,5,7}^k start with n!:

A374923, A374922, A386547,A387304.

点评

可以有——k ! start with {2,3,5,7}^n 。  发表于 3 天前
可以有——k ! start with n!: {2,3,5,7}^n 。  发表于 3 天前
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 3 天前 | 显示全部楼层
7个代码放一起。比较一下。——班门弄斧。
              
AbsoluteTiming[Module[{t = IntegerDigits[8!], k = 0}, While[Take[IntegerDigits[5^k], UpTo[Length[t]]] != t, k++]; k]]——{355.19, 149704}。——这是10#的公式。

AbsoluteTiming[Module[{t = StringJoin[ToString /@IntegerDigits[8!]], k = 0, u = StringLength[ToString[8!]]}, While[True, s = ToString[5^k]; If[StringLength[ s ] >= u && StringTake[s, u] == t, Break[]]; k++;]; k]]{293.302, 149704}。

AbsoluteTiming[a = 8!; b = Ceiling[N[Log[10, a]]]; k = 0; While[c = Ceiling[N[Log[10, 5^k]]]; a != (5^k - Mod[5^k, 10^(c - b)])/10^(c - b), k++]; k]——{157.942, 149704}。

AbsoluteTiming[a = 8!; b = IntegerLength[a]; k = 0; While[c = IntegerLength[5^k]; d = Quotient[5^k, 10^(c - b)]; d != a, k++]; k]{96.0528, 149704}。

AbsoluteTiming[a = 8!; b = IntegerLength[a]; k = 0; While[True, c = 5^k; d = IntegerLength[c]; If[d >= b && Quotient[c, 10^(d - b)] == a, Break[]]; k++;]; k]{60.369, 149704}。

AbsoluteTiming[a = 8!; b = IntegerLength[a]; k = 0; While[c = 5^k; d = IntegerLength[c]; ! (d >= b && Quotient[c, 10^(d - b)] == a), k++]; k]{60.2769, 149704}。

AbsoluteTiming[a = 8!; b = IntegerLength[a]; k = 0; c = 1;  While[True, d = IntegerLength[c]; If[d >= b, t = Quotient[c, 10^(d - b)]; If[t == a, Break[]]]; k++; c *= 5;]; k]{24.0573, 149704}。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 3 天前 | 显示全部楼层
对数运算不一定快。再譬如——A059449——a(n)! begins with 2^n.——1, 2, 8, 14, 89, 164, 18, 266, 304, 292, 2522, 557, 26438, 27045, 57134, 21936。

AbsoluteTiming[Do[a = 2^n; b = Floor[N[Log[10, a]] + 1]; k = 1; While[c = Floor[N[Log[10, k!]] + 1]; a != (k! - Mod[k!, 10^(c - b)])/10^(c - b), k++]; Print[k], {n, 0, 15}]]——A059449——自带的公式。
{419.582, {1, 2, 8, 14, 89, 164, 18, 266, 304, 292, 2522, 557, 26438, 27045, 57134, 21936}}

AbsoluteTiming[Table[a = 2^n; b = IntegerLength[a]; k = 1; f = 1; While[True, f = f*k; c = IntegerLength[f]; If[c >= b, d = Quotient[f, 10^(c - b)]; If[d == a, Break[]]];  k++;]; k, {n, 0, 15}]]
{38.2066, {1, 2, 8, 14, 89, 164, 18, 266, 304, 292, 2522, 557, 26438, 27045, 57134, 21936}}

我这里出来就比A059449还要多了——{1, 2, 8, 14, 89, 164, 18, 266, 304, 292, 2522, 557, 26438, 27045, 57134, 21936, 224636, 529309}
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 前天 11:46 | 显示全部楼层
任意给定一个连分数
\([a_0;a_1,a_2,...]\)
我们定义\(\frac{p_n}{q_n}=[a_0;a_1,a_2,...,a_n]\), 也就是连分数前n+1个数截断以后代表的既约分数。
链接中给出结论\(\frac{p_{n-1}x+q_n}{q_{n-1}x+q_n}=[a_0;a_1,a_2,...,a_n+x]\).
本题中,我们的目标搜索一个整数k使得5^k前面h位为给定整数X, 那么实际上我们的目标就是寻找整数k使得
\(k\lg(5)\)的小数部分在\(\lg(X)\)和\(\lg(X+1)\)的小数部分之间。我们记\(LOW=\{\lg(X)\}, UPP=\{lg(X+1)\}\)
假设\(m=\lfloor k\lg(5)\rfloor, t=\{ k\log(5)\}\), 于是\(LOW\le t \lt UPP\), 而且\(\lg(5)=\frac{m+t}k\)
我们现在的目标是搜索\(\frac m k\)的连分数表示,并且假设它的连分数表示为\([a_0;a_1,a_2,...]\),具体多少项我们不清楚,但是我们知道,前面少数几项会和\(\lg(5)\)的相同,显然第一项应该是0,也就是\(a_0=0\).
假设我们已经确定了前面h项,也就是确定了\([a_0;a_1,a_2,...,a_h]\),现在我们设\(\lg(5)=[a_0;a_1,a_2,...,a_h+x*]\),
而\(\frac m k=[a_0;a_1,a_2,...,a_h+x]\), 对于给定的\([a_0;a_1,a_2,...,a_h]\),我们可以非常容易计算出对应的\(x*\),
由于我们知道\(\frac {LOW} k \le \lg(5) -\frac m k \lt \frac{UPP}k\),  假设\(x=\frac uv\), 我们根据链接公式有\(\frac m k = \frac{ p_{h-1} u + p_h v}{q_{h-1}u+q_h v}\), 或者说\(m=s(p_{h-1} u + p_h v), k = s(q_{h-1}u+q_h v)\),其中\(s\)是m,k的公约数。
由此可以得到不等式
\(\frac{LOW}s \le \frac{(p_hq_{h-1}-q_np_{h-1})(u-x^*v)}{q_{h-1}x^*+q_h}\lt \frac{UPP}s\)
由于\(p_hq_{h-1}-q_hp_{h-1}=\pm 1\), 由此我们很容易得到一个关于\(v(x-x^*)=u-x^*v\)的带参数s的上下界a/s,b/s。而我们的目标是优先寻找最小的v*s, 对应最小的k. 于是我们对应在搜索\(\frac uv\)的连分数形式时.
由于我们要求\(x-x^*\)在集合\((\frac a{vs}, \frac b{vs})=\cup_{k=1}^{\infty}(\frac ak,\frac bk)\), 我们可以放大到区间(0,b) 或者(-a,0) 等(根据a,b的符号),
于是我们可以估计x的范围在\((x^*,b+x^*)\)或者\((x^*-a,x^*)\)等之间,由此我们可以确定x进行连分数展开是第一项的所有可能范围, 因为第一项为k意味x在\([\frac1{k+1},\frac1k]\).
如果其中第一项的可能性是唯一的,那就比较好,我们等于确定了\(a_{h+1}\)的唯一取值,不然我们还需要枚举多种可能的\(a_{h+1}\),从中找出对应\(k = s(q_{h-1}u+q_h v)\)最小的解,并且用当前已经找到的最优值限制搜索深度,不过感觉代码会非常复杂
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 前天 15:48 | 显示全部楼层
我们现在举个实际的例子,以n=6为例子,X=720, 于是LOW=\lg(720)-2=0.85733249643126846023127249068370969871, UPP=\lg(721)-2=0.85793526471942903588328745317287561787.
取\(a_0=0\), 于是\(x^*=\lg(5)=0.69897000433601880478626110527550697323\)
于是要求第一步的{x-x^*}在集合\(\cup_{k=1}^{\infty} [\frac {LOW}k, \frac{UPP}k)\)
分别对应x在区间[1.5563025007672872650175335959592166719,1.5569052690554478406695485584483825911), [1.1276362525516530349018973506173618226,1.1279376366957333227279048318619447822),
[0.98474750314644162486335193550341020613, 0.98494842590916181674735692299979884585)
而后面那些区间倍[0.69897000433601880478626110527550697323,  0.98494842590916181674735692299979884585) 统一包含 (\(vs>=3\))。
0.98494842590916181674735692299979884585)
由于x必须在[0,1),所以前面两部分数据淘汰,我们只余下区间[0.69897000433601880478626110527550697323,  0.98494842590916181674735692299979884585),这个区间在[1/2,1)之间,由此得出x连分数展开式第一项只能为1 ,得到\(a_1=1\)。 对应\(p_0=0,q_0=1,p_1=1,q_1=1\) ,当然如果这时x=1/1不落在上面任何区间\(\cup_{k=1}^{\infty} [\frac {LOW}k, \frac{UPP}k)\)内,需要继续搜索。

由此得出下一轮\(x^*=\frac1{\lg(5)}-1=0.43067655807339305067010656876396563208,p_1q_0-q_1p_0=1\)
得到\(\frac{LOW (x^*+1)}{vs} \le x-x^* \lt \frac{UPP(x^*+1)}{vs}\), 我只知道这时至少要求\(x \in (x^*,1)\), 对应\(a_2=1或2\)
而\(a_2=1\)时,区间大于0.5,对应\(vs\le 17\), 同样1/1落在任意区间内,需要继续搜索, 而且搜索的限制式分母不超过17. 可以直接穷举或继续搜索。
   如果穷举,可以对每个可能的vs,计算x的范围,然后利用x=us/vs, 代表这个区间范围乘上vs后必须包含整数us才可以。
对应区间范围(分别对应vs=3,4,...,17)
[0.83953172644631194787322185412064422214,0.83981918193291906904943474572887593511)
[0.73731793435308222357244303278147457463,0.73753352596803756445460270148764835935)
[0.67598965909714438899197573997797278612,0.67616213238910866169770347494291181390)
[0.63510414225985249927166421144230492711,0.63524787000315605985977065724642078359)
[0.60590020166178686375715597677397074210,0.60602339687033277283267578746321290480)
[0.58399724621323763712127480077272010335,0.58410504202071530756235463512580699571)
[0.56696161419769934973781166388285849543,0.56705743269323505679654929441893573309)
[0.55333310858526871983104115437096920910,0.55341934523125085618390502185343872299)
[0.54218251308418911354368346477033252028,0.54226091003508196477355970793621389654)
[0.53289035016662277497088539010313527959,0.53296221403827455526493861300519320783)
[0.52502775077483587310159471153858376825,0.52509408665636059337302845575586800970)
[0.51828837986758995721363127276896818709,0.51834997747186291175139117811358926844)
[0.51244759174797683011072962583530135009,0.51250508284529825434597220415694769268)
[0.50733690214331534389569068476834286771,0.50739080004705417911623060194488631389)
[0.50282747013920226782359750147396773621,0.50287819757801528920763507058130274438)
[0.49881908613554620020395911632341206376,0.49886699538331405373332793159145068258)
对应每个区间乘上vs后取值为
[1.6572420631921497422794524248340014023,1.6581044296519711058080910996586965412)
[2.0879186212655427929495589935979670344,2.0887809877253641564781976684226621733)
[2.5185951793389358436196655623619326664,2.5194575457987572071483042371866278053)
[2.9492717374123288942897721311258982985,2.9501341038721502578184108059505934374)
[3.3799482954857219449598786998898639306,3.3808106619455433084885173747145590695)
[3.8106248535591149956299852686538295627,3.8114872200189363591586239434785247016)
[4.2413014116325080463000918374177951947,4.2421637780923294098287305122424903336)
[4.6719779697059010969701984061817608268,4.6728403361657224604988370810064559657)
[5.1026545277792941476403049749457264589,5.1035168942391155111689436497704215978)
[5.5333310858526871983104115437096920910,5.5341934523125085618390502185343872299)
[5.9640076439260802489805181124736577231,5.9648700103859016125091567872983528620)
[6.3946842019994732996506246812376233551,6.3955465684592946631792633560623184940)
[6.8253607600728663503207312500015889872,6.8262231265326877138493699248262841261)
[7.2560373181462594009908378187655546193,7.2568996846060807645194764935902497582)
[7.6867138762196524516609443875295202514,7.6875762426794738151895830623542153903)
[8.1173904342930455023310509562934858834,8.1182528007528668658596896311181810223)
[8.5480669923664385530011575250574515155,8.5489293588262599165297961998821466544)
[8.9787435504398316036712640938214171476,8.9796059168996529671999027686461122865)
可以看出,每段都不包含整数点,所以全部淘汰

由此只能\(a_2=2\).
对应\(vs \ge 18, x^*=0.3219280948873623478703194294893901758, p_2=2,q_2=3, p_2q_1-q_2p_1=-1\),

并且得到不等式\(-\frac{UPP(x^*+3)}{sv}\lt x-x^* \le -\frac{LOW(x^*+3)}{vs}\),得到\(sv\ge 9\)时,才可能\(x \in [0,1)\)
对应区间
sv=9, (0.0052615105022403708967412881912270047138,0.0054839941590345684217462674442520823159], 但是这些值小于1/9=1/(vs)<1/v不符合要求
sv=10, (0.037128404231867346366603583648765891669,0.036928168940752568594099102321043321826], 值小于1/10=1/(vs)<1/v不符合要求,
sv=11, (0.062837253117717093982846404790893035826,0.063019285200548710139668660543368099319], 只需要1/11=1/(vs)<1/v不符合要求.
sv=12, (0.084428156598520865140135823515767797493,0.084595019341116513283889557955536605695],
...
总体上,我们会得到0.084428156598520865140135823515767797493<x<x^*=0.3219280948873623478703194294893901758, 由此得到\(3\le a_3 \le 11\)

至此,我们有[0,1,2,3],[0,1,2,4],[0,1,2,5],...,[0,1,2,11]. 共11个候选。
我们应该优先搜索靠近x^*的结果,也就是\(a_3=3\),也就是递增顺序搜索.

在\(a_3=3\)时,计算有\(x^*=0.10628371950538987600239849178357475597, p_3=7,q_3=10,p_3*q_2-q_3*p_2=1\),
由此得到不等式

\(\frac{LOW3}{vs}=\frac{LOW (3x^*+10)}{vs} \le x-x^* \lt \frac{UPP(3x^*+10)}{vs}=\frac{UPP3}{vs}\)

由于\(x^*\)比较小,我们可以得到在\(vs\ge 10\)时,可以得到\(x^*<x<0.99157434953358557384915516687296267908\)
类似上面方案,对于一些较小的n,我们可以先枚举区间\([LOW3+vs*x^*,UPP3+vs*x^*)\)是否包含整数点,然后对于余下范围,在对a_4进行分类。
实际上对于vs=199,可以有区间[29.997146605605939738505051036767673284,30.003366481854542302944866615825255670)
也就是我们得到可以有解us=30,vs=199, 即对应我们的目标解  m/k=[0;1,2,3+30/199]=1453/2080.
当然如果我们这里穷举范围小一些,比如知道vs=100. 那么余下区间在(0.10628371950538987600239849178357475597, 0.19393625713194390549217638040628643152), 两者倒数5.1563333993790198460217886837373475432到9.4087787353856001368325082430578137867,
对应需要继续搜索空间\(a_4=5,6,7,8,9\),优先搜索靠近x^*部分,也就是\(a_4=9\),也就是递减顺序搜索。

搜索过程中可能会提前找到一些其它不是最优的解,比如\(a_4=9\)可以找到x=3/211,对应m/k=13736/19653, 那么后面继续搜索中,我们就可以用\(vs q_h\lt us q_{h-1}+vs q_h\lt k\)来做剪枝,得出vs的一个搜索上界了

评分

参与人数 1威望 +8 金币 +8 贡献 +8 经验 +8 鲜花 +8 收起 理由
northwolves + 8 + 8 + 8 + 8 + 8 很给力!

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 昨天 15:28 | 显示全部楼层
现在代码比较粗糙,验证了N=10~15和wolves相同,另外得出
N=16
700785493882
不过N=17还是出不来,可能需要调整
  1. #include <gmpxx.h>
  2. #include <iostream>
  3. #ifndef N
  4. #define N 10
  5. #endif

  6. #define PREC 3000
  7. #if N==10

  8. mpf_class lgNfac("0.5597630328767937511747612399600555908932271418594001136856622817018191866832707151595484631470708116273899442969485014350620900263856693662453264103027331671177926665919099403914095284940514926910190038154483975520500476678019753243237638601809284996882257656525282203638995000349747345868231323035409926787747523076550826236151137832855580203281034749161818795132939859815213483271758017382846906260428003724227082312295982498627429482231923562145809314578884266948456955425575151223348285543101296610017182622467743146323189318357758827168396998081811369036806837543875120253218217950902097453860936609562706516413101610406306032824675179298025117697243828023871529437211881090000134242485533909791004821407786459813112778974700758110724783487012193394826852717625256134653300428714888746662799516304918143451336441941097801424292472369570765450275724395583098555691235046595365459031065068100081118324693217085246233632885174190388990055960853540334624097213749793093543555133809798468884685371002",PREC,10);
  9. mpf_class lgNa1fac("0.5597631525566940111241786939544799467484747745028933440989937595315976904390877869792913678149940658639563780233135458496085838538867596497996492452805059675370107543531774119982840893409927092665680560488344871416125071286159968924372725815610297539717991318596071753747897880133927606701531700223278208492917279643346559011621406447304101952660661034232112959696393925223494547623725925612695410807375205592293008292619168888826853181722537658449214837495411314701009829451354075244431048323597459612753405388901524365372634129476015584173732335666816138143502093177725484063131657110986664492470835170934602414831746072663353972364086434857470203486954791550406639889613725556775900748837146559905640327739532770141278823149469329865834263494577972489809710643439244913986354857263803099261158738371578551210062211447791872798922459473384067091969993242986110657858150491948394933385676154884871479026599459821640883508715006047530987996752936228350555461675745864089624124185159496815690306991175",PREC,10);
  10. #elif N==11

  11. mpf_class lgNfac("0.6011557180350187919249612112030798325999293323258532082822013003815722189165200627308432495328019906308647572891646830292346685461006768020139157535245259462375327971604393772289963575065535432058725893980165601104069556492804409649896091901048435268183468516139102955794532845373104557119770980883941136350522388102716123605334334721762009454980187580147625189581659453737407766135561994238024117948052508909395040720986700699315605512216152988290765447700321069271183775129315138343581146521187869763012701992114603629057098215496963652371643167939769558826294545711990686297356125280073813118322374502262192478038358106054408992069992524498290961758463823386793038645203065070281652829853264745796653636762147112027067721889709257247793249150031156803422021657763024705770074706212228207303978134185592810519711003502902887860714752436345679124113611700267478347458848102172260955864466364938900866680671575769051738544618993319757509965601023809298421679203894949619411650720597411737319377178737",PREC,10);
  12. mpf_class lgNa1fac("0.6011557289150110874811105665359089961351038515110119290833738457560862941418995834309068620055602929887426537237758034453323420148189637580570031477337836227346736961530354098923189190060786776492012362675660636536945314560720992161843647745963546887933405141944445688919585518488047865731600066042934275808724434625860811765383811618580575006175860471455385013224198556226610842704968915924678216820189585670336575341409564344101554593553718544506737398309441171767597644610848493364940539215076587624149856776564022367256352590147850856678877266892641607737067343647250403213469082768394819162353238322682767882768484131385259529601549006461994799922473170917046062330591416251922230150182695950291598959772833504102404332488207039863642345210268128082566220942356697122472791387788258585104758155683388806470084827751740176974282744567818328201485331836623850107980872645435377681964838956335412982189903091611948770251479347352801790957312671246923454819763780407847338604077991237194465003548780",PREC,10);
  13. #endif

  14. #include "powpref.h"

  15. mpf_class lg5("0.6989700043360188047862611052755069732318101185378914586895725388728918107255754905130727478818138279593155228085690046209052321188664764940003076662953044249354970357458065973381802656883970564988160971018214173828455604681380709536461153004797606891503875374595997366874053785211541526817173273160176738034572064923686824516490728610350530822142310819492099924004519121845402854149680351223738775077091708818090485010028283801395223234999321794820874426713713316579995970794901629154277751045057024378502927553402913863103907780905172387856085034717648321735076851959722537567558366884612617406961169606193667838697609481194178680843145383070946984948680730146215115812816799342464305316070282578679890941031091494143753590127816031233514601437648387226973610721217391501633189696915685844391860563823254511433365754618762660675775304056509397879554957031725393115214538843152315893562020499534030082254342459135981535920543470455658922591706000254599262782983198051109445143089305996245883100365842",PREC,10);
  16. mpf_class one("1.0",PREC,10);
  17. #define LOW lgNfac
  18. #define UPP lgNa1fac
  19. #define M   lg5
  20. #define ONE one
  21. mpz_class bestGK=-1;
  22. #define CHECK_COUNT 5000

  23. bool getKrange(int sgn, const mpf_class& L, const mpf_class& U, const mpf_class& xp, const mpz_class& q1, int& uk_infty, mpz_class& LK, mpz_class& UK, mpf_class& LX, mpf_class& UX, int& retz, int &checkl, int&checku)
  24. {
  25.   retz=0;
  26.   checkl=checku=0;
  27.   if(sgn>0){
  28.     if(xp>=ONE){
  29.       return false;
  30.     }else if(0<xp&&xp<ONE){
  31.       uk_infty=1;
  32.       LK=(mpz_class)ceil(L/(ONE-xp));
  33.       if(bestGK>0){
  34.         mpz_class UB=bestGK/q1;
  35.         UK=UB;
  36.         uk_infty=0;
  37.         if(UK<LK)return false;
  38.       }
  39.       if(uk_infty){
  40.          LX=xp;
  41.       }else{
  42.         checku=1;
  43.         if(UK<=CHECK_COUNT){
  44.           LX=xp+L;
  45.         }else{
  46.           LX=xp+L/(UK-CHECK_COUNT);
  47.         }
  48.       }
  49.       UX=xp+U/LK;
  50.       if(UX>=ONE)UX=ONE;
  51.     }else{
  52.       uk_infty=0;
  53.       LK=(mpz_class)ceil(L/(ONE-xp));
  54.       UK=(mpz_class)floor(-U/xp);
  55.       if(bestGK>0){
  56.         mpz_class UB=bestGK/q1;
  57.         if(UK>UB)UK=UB;
  58.       }
  59.       if(UK<LK)return false;
  60.       checku=1;
  61.       if(UK<=CHECK_COUNT){
  62.         LX=xp+L;
  63.       }else{
  64.         LX=xp+L/(UK-CHECK_COUNT);
  65.       }
  66.       UX=xp+U/LK;
  67.       if(UX>=ONE)UX=ONE;
  68.     }
  69.   }else{
  70.     if(xp<=0){
  71.       return false;
  72.     }else if(xp<=ONE){
  73.       uk_infty=1;
  74.       LK=(mpz_class)ceil(L/xp);
  75.       if(bestGK>0){
  76.         mpz_class UB=bestGK/q1;
  77.         if(UK>UB)UK=UB;
  78.         uk_infty=0;
  79.         if(UK<LK)return false;
  80.       }
  81.       checkl=1;
  82.       LX=xp-U/(LK+CHECK_COUNT);
  83.       if(LX<=0){retz=1;LX=0;}
  84.       UX=ONE;
  85.     }else{
  86.       uk_infty=0;
  87.       LK=(mpz_class)ceil(L/xp);
  88.       UK=(mpz_class)floor(U/(xp-ONE));
  89.       if(bestGK>0){
  90.         mpz_class UB=bestGK/q1;
  91.         if(UK>UB)UK=UB;
  92.       }
  93.       if(UK<LK)return false;
  94.       checkl=1;
  95.       LX=xp-U/(LK+CHECK_COUNT);if(LX<=0){retz=1;LX=0;}
  96.       UX=xp-L/UK;if(UX>ONE)UX=ONE;
  97.     }
  98.   }
  99.   return true;
  100. }

  101. bool check_range(const mpz_class& f, const mpz_class& t, const mpf_class& L, const mpf_class& U, const mpf_class& xp, mpz_class& RU, mpz_class& RD)
  102. {
  103.   mpz_class i;
  104.   for(i=f;i<=t;i++){
  105.     mpf_class a=L+i*xp;
  106.     mpf_class b=U+i*xp;
  107.     mpz_class ai=(mpz_class)floor(a);
  108.     mpz_class bi=(mpz_class)floor(b);
  109.     if(bi>ai){
  110.       RU=bi;RD=i;
  111.       return true;
  112.     }
  113.   }
  114.   return false;
  115. }

  116. #ifdef DUMP
  117. int level;
  118. #define UPLEVEL (level++)
  119. #define DOWNLEVEL (level--)
  120. #else
  121. #define UPLEVEL
  122. #define DOWNLEVEL
  123. #endif

  124. void search(const mpz_class& p0, const mpz_class& q0, const mpz_class& p1, const mpz_class& q1)
  125. {
  126.   mpf_class xp=(q1*M-p1)/(p0-q0*M);
  127.   mpf_class mul = q0*xp+q1;
  128.   int sgn = mul<0?-1:1;
  129.   if(mul<0)mul=-mul;
  130.   mpf_class L = LOW*mul;
  131.   mpf_class R = UPP*mul;
  132.   sgn*=(p1*q0>q1*p0)?1:-1;
  133.   mpz_class LK,UK;
  134.   mpf_class LX,UX;
  135.   int uk_infty,retz;
  136.   int checkl,checku;
  137.   UPLEVEL;
  138.   if(!getKrange(sgn, L, R, xp, q1, uk_infty, LK, UK, LX, UX, retz,checkl,checku)){
  139.     DOWNLEVEL;
  140.     return;
  141.   }
  142.   if(retz){//solution p1/q1 found
  143.         if(bestGK<0||q1<bestGK){
  144.           std::cout<<q1<<std::endl;
  145.           bestGK=q1;
  146.           DOWNLEVEL;
  147.           return;
  148.         }
  149.   }

  150.   if(checkl){
  151.      mpz_class end=LK+CHECK_COUNT-1;
  152.      if(end>UK)end=UK;
  153.      mpz_class RU,RT;
  154.      if(check_range(LK,end,L,R,xp,RU, RT)){//f
  155.         mpz_class nv = q0*RU+q1*RT;
  156.         if(bestGK<0||nv<bestGK){
  157.           std::cout<<nv<<std::endl;
  158.           bestGK=nv;
  159.         }
  160.      }
  161.    }

  162.   if(checku){
  163.      mpz_class i;
  164.      i=UK-CHECK_COUNT+1;
  165.      if(i<LK)i=LK;
  166.      mpz_class RU,RT;
  167.      if(check_range(i,UK,L,R,xp, RU,RT)){
  168.        mpz_class nv=q0*RU+q1*RT;
  169.        if(bestGK<0||nv<bestGK){
  170.          std::cout<<nv<<std::endl;
  171.          bestGK=nv;
  172.        }
  173.      }
  174.   }

  175.   mpz_class lowN=(mpz_class)floor(ONE/UX);
  176.   mpz_class uppN=(mpz_class)floor(ONE/LX);
  177.   mpz_class i;
  178.   for(i=lowN;i<=uppN;i++){
  179.     mpz_class p2=p0+p1*i;
  180.     mpz_class q2=q0+q1*i;
  181. #ifdef DUMP
  182.     int j;
  183.     for(j=0;j<level;j++)std::cout<<" ";
  184.     std::cout<<"a["<<level<<"]="<<i<<"("<<p1<<"/"<<q1<<","<<p2<<"/"<<q2<<std::endl;
  185. #endif
  186.    search(p1,q1,p2,q2);
  187.   }
  188.   DOWNLEVEL;
  189. }

  190. int main()
  191. {
  192.   int sgn=1;
  193.   mpf_set_default_prec(1000);
  194.   mpf_class xp=M;
  195.   gmp_printf("%F\n",xp.get_mpf_t());
  196.   mpf_class L=LOW, R=UPP;
  197.   int uk_infty, retz;
  198.   mpz_class LK, UK;
  199.   mpf_class LX, UX;
  200.   mpz_class p0=0,q0=1;
  201.   int checkl, checku;
  202.   if(!getKrange(sgn, L, R, xp,q0, uk_infty, LK,UK,LX, UX, retz, checkl,checku)){
  203.     std::cerr<<"Fail to search\n";
  204.     return -1;
  205.   }
  206.   if(retz){//finished search;
  207.     std::cout<<"Best result is 0\n";
  208.     return 0;
  209.   }
  210.   if(checkl){
  211.     mpz_class end=LK+CHECK_COUNT-1;
  212.     if(end>UK)end=UK;
  213.     mpz_class RU,RT;
  214.     if(check_range(LK,end,L,R,xp,RU,RT)){
  215.        bestGK=RT;
  216.         std::cout<<RT<<"\n";
  217.     }
  218.   }
  219.   if(checku){
  220.     mpz_class i;
  221.     i=UK-CHECK_COUNT+1;
  222.     if(i<LK)i=LK;
  223.     mpz_class RU,RT;
  224.     if(check_range(i,UK,L,R,xp, RU,RT)){
  225.       bestGK=RT;
  226.         std::cout<<RT<<"\n";
  227.     }
  228.   }
  229.   if(LX>UX){
  230.     std::cerr<<"No more result found\n";
  231.     return 0;
  232.   }
  233.   mpz_class lowN = (mpz_class)floor(ONE/UX);
  234.   mpz_class uppN = (mpz_class)floor(ONE/LX);
  235.   mpz_class i;
  236.   for(i=lowN;i<=uppN;i++){
  237.     mpz_class p1=1,q1=i;
  238. #ifdef DUMP
  239.     int j;
  240.     for(j=0;j<level;j++)std::cout<<" ";
  241.     std::cout<<"a["<<level<<"]="<<i<<"("<<p0<<"/"<<q0<<","<<p1<<"/"<<q1<<std::endl;
  242. #endif
  243.     search(p0,q0,p1,q1);
  244.   }
  245.   std::cout<<"Final:\n";
  246.   std::cout<<bestGK<<std::endl;
  247.   return 0;
  248. }
复制代码



其中关于n!的对数需要预先生成,余下还有
  1. #if N==12
  2. mpf_class lgNfac("0.6803369640826436196474669039071811953364379594407661557329218629410177502490302047476661094167906362226984044402787196809871734505929511653036441270054465062997816775219975400424210304252362293868929522632025468139481715736939373161976391771397157661302940978205316753930551774578136672494354819261141816530211011335645432798139172032595543253883622742642755307298771021473452612088261212387962089334998612641501912012204058881683074851657156647916621921994363817284244174316701159488500174899528364781858809056425191121900834044324754090127971080683636926939254610528691481657309813108224952087925376077837271091234194790213078376317448275052701281411155197416905596464150410852663991726956129991523211257106728418535113520894302464790250535720070604513252263923426077953864715297299629423756143965637616379007835657569444013740292555904497976092012750035013403066089988399251639380922348040010012122776503121670386816966351670977851930116144761250071175677204258385623610810095627839276800963807105",PREC,10);
  3. mpf_class lgNa1fac("0.6803369649893096546876814048373598328245022358903928781241989742599995152346372807883471852196552739121086679909949902749676414037490586059231310550052182140387745464307316838169181157153022960532871632038954027218255469671485488073514798897617437628319218103778731493965531968031650207430659888328737210530831928479341616957925101713271860008890809196333009454077639639678561989942195160612907961670842690102890907027950893342761507315517808598601635609675084140045368344431495754994321350790241459081333749593261304549698012355701062256907757535171221070932624307231911112578772229272992144216524127062511971852297052018681744702201280988372046974527260006801708129250382057245447795196411412776782663336787480517213517179535861213189289270599153771108368196351806148459045494689828277104207718913328588429396888057347895606196081381698830724948978184179346555610483685938201072717182269934420642913521898436366218532694895462136363413746520586280851371515256802502025293815831783781727769197744055",PREC,10);
  4. #elif N==13
  5. mpf_class lgNfac("0.7942803163894803888539720618495096261661671478278344275341128379385708418792244455553135636756902958956763599270078632122241334948045293212197431043372213409681048401454297936164050914886767191440908058762362853297273841248845225199166798403503242234488552499760168862369184264599811537920643754530708428243674723221659157958907942063874831797041517235833759154898822836553315567708789490609788937174688469221537640990098172146690046484108633015500757248604683264968822041563526417105027699157972128612084836223247554916009819823951172792195633852103632528182860795383256448306683953583656006669959761429972174554193955027624928210172583406066399859243531247299504256210605698479171883590822764172126630131020335609010838342422227025213891806762276031862079987428566511919765610679879300981074281679936298698762297640159437616545768491062025258364841826005276663438626250855107658638803837023329229911942945651771293207044729103828937203444152300795636555084647695239253346788575725032629111994583114",PREC,10);
  6. mpf_class lgNa1fac("0.7942803164592239300781125139000298425719773553920629682350976489960174982539659257810485818954522893061583865610868299638345533069276893733754917948458808975891601606427394419711087482328044129088542551613312713757659867080084540192347634819450806351165850187539436548073248227962466216265432309884745862943314625930488233607898025522786213665236075999115079277879292891141429602421176505358525620416569209421747285928630434765859227701304574493481945740902072583448306897705866118617865341879523677083722034449748906566604141889210391355106562453640360898504321403960933100690869141338183809470178139202724982065381088068439757814562498465180086193619084538806793950819722608218769647528067632592149829712789862759004476985324860244112555480385068167495512301317060172657443619718131398857458443969055643404800662242379059430936160797819377985337206327082288869573527129327280319479859640866995118390421916708225160757120378031938624758579122307179626965505877570213139195665127377886241679605136091",PREC,10);
  7. #elif N==14
  8. mpf_class lgNfac("0.9404083520677184147799272151666388464179294256139083753481752527839324710557281211083523436413650412777851471937932041775605075466216717523052330218303589930304108780536872907395751312428062976145800669254764839190758866906182995784699778509524816984381294184232603772841942192494805841839019904576785420208951558595432273874367544822344684338281878876538757707149759059888141104720709581721667632656056931563065783107439870511661928903578576554680517385866256590125049487854852596285866618052540372011208807928264206501646342484282252704770464095057029459562596817781160774417412319568115696024525182918653563101730881457814020663477999433933942743483967388764280206933577798412249746296685989362971731647956522360360507154995106859454243639656185272736665038919887380557140779653497527328752836178963359033961038201978082932424423837389458642304856125966699883774149308183928562716917313892542775759121285307843254519229053156019971745122820203809082374840910397135001758276941495132202421032343117",PREC,10);
  9. mpf_class lgNa1fac("0.9404083520727000962963086808645228506874973980191117832101219348881457560631825969886829646106702645917585608108412462617830232656822999062582750057438361098076968721946949500827236762847969911880918938213892000250991307195446399634000284071079692232685893450808267629796983290752542758686866430479231704227760307302187326290397582297865323794289678532826631680953135584919656202608590719986817887140724035650525968184888093236740958443450329848157738383070569038107208336343108673658296854691524457938781753487380040988804804819255654668407033293603015445401104318807438228500714554882337593446878449208720500991673242825479444760832453018684575188051398663906111276061943034627347187378510262466817614861488760615015595593636419489487205858484865681468534698201467774210890572677456114948612256937836378970473989336748450504256182207861464411227440697863285415206209504260955473803625217450263472118880084932705914043368119219480813843676367901751270992250474316594548778114227538684861012963979395",PREC,10);
  10. #elif N==15
  11. mpf_class lgNfac("0.1164996111233996568612162236972611288498684083424956988676134319620534345649647346643934471707951707475653627706142546920287088077133755975958843941971928278991508656526652455675406012266801532920487210961267227711537840194460087906163537394266360052012392770086809671600122477334462537665123562434516314492356376599422056616644567964229810603612246497510187596880427993160394513122449853542821929274278161749444109428742080098215067948017545598768997140301439287878097799426943504893618899566051940165563697652783535583801801655825560306095047111953841792880787438477029182480041108050104687215011693312809645231019546585555226090254891333701194017982700953233038112190975542222025001141999702344737962091394378149299813724383148160697244730539169882127416113349202608310134989335331986078380583701885146136749260119900912040331327553027139132911420135396621987838424065010464888822661257066215977261980144231152534205412416245045042933048481942013652917187860354724334292865584443548479551920069012",PREC,10);
  12. mpf_class lgNa1fac("0.1164996111237317689623100992113834504723136688953161078635694523693258266031272670572748526316480690623874211211927894103675976828508456624609529067265096585870387726872702471557224814727524206484716676273200439352802288745856957566535592948941722818932064346707327175460008719831039260378908650558342662246014999462846184669321391875642316104276799689080893821336848409313084535282044716405568936061099935893852005512936221053398183010083096346560034372308927313775619950414150090888972161358997982088746062041480412700789292869716309974094574983520477160058400503860388916750492687254794033203400345780330563023577232382091031924005500408950276825900898476847345363818295730901838859448133289386609340068813813630328588850480934575517252796919247269679990499084606175492033387175729555952440745065070446229887051937238221830059095710274548993529767282809373861984209688880499749209130814369143071598792785785399475670524595925616376930025528006569028832406254526692788439290190658052451790937796519",PREC,10);
  13. #elif N==16
  14. mpf_class lgNfac("0.3206195937793244377161718025952332359226279341909298641093232764704861916626627726121024556435398589103032715363382362084077803322474696215946537290159751281571627226694388562148195384730919272967843326888410532397715421468937249760318925375075932485996891271702820204103907336488296430396430469793809362354068116904674758550681653522827687315043003219541787900862351505778783096523728448647866828965911326477082169388628944892634175008020258419485499433446586021558113916247336988276507895385823842651551987439171880131386170532204870754670706973083248505940480030638139032209807640511654217587167014888034973876229108660778511366882309801417406078187978032648177648939708344852167779877718572030018398327270012172724799363871884035763186324788576333219521670464333042303602230547669242700813141446592128091015797101425861397628226336801101541393200307269720415377565909637855625248413175068079856932962774394608608061730242363222407242681657940995255866055927562519896512293227219563496019518605643",PREC,10);
  15. mpf_class lgNa1fac("0.3206195937793451947224901772554389700371849027443560565872429807304242121260830424979386911309252170665727683807762909376035842124320358851993141753211164055175864128527032528996071854348603134570819153767421127829479620612947552539199830303063034525359984626351319502260637235571729504046777433614726585184984552631110242498462631039121103464897835612489582474201635173092591573306128391604237023339837382188072404486236174010566077100402967396806458862692919723930062272021911354461075180456728647827186333470457648414052931454276301325951771417706829111530781933481920636634392351683266822552124813071798691107265043679516729181590897464003892718934055372480319290223135686154094730265164421195598289421457293665620712252996663986883362296914542311734593381921836667650032246922242138788235418973528616216282798423553514391900466042709620459858662120869480735220787961522463786925451896824636770805076843702200401675132712393727018575613267330423476482543393456888350109487759703113060000316047627",PREC,10);
  16. #elif N==17
  17. mpf_class lgNfac("0.5510685151575983662563416969235702659301953126159762614896917588151802539197809584077871257420049984157890567861372012205825140753625842956485420090675252517189242454396160295891249920787080301037034952668473738987461627800888385701731521370152207784600999660728449693954666597391252681072101274886622259937955273442042337563821712929161513238746005833851132593910405857457270204404108582850578655076282754482915853575425933560042336436979028216845340662973472463920513795825950953385855746070719515440745868762865315158241262398386417192709614113558595389460677174650873275915383738589246205303335915258287039472850530367119243346333853347981903489573622222739877861720254714898210507062121633290484215249702421495025432485171124389245785702638303088228916524005515758632651973761985407388081503087438443644519470230903018545759164760099823769261472438531403189514806945414927424986206538332548636706007335854679314794985057160076460612342608526707104367974252941077509599005671198271659597680847322",PREC,10);
  18. mpf_class lgNa1fac("0.5510685151575995872567133660486917495803643928118332270076706769798979036027312131221448356361168829642734335272291105376888177841332538661145665277172339569338308734552740419561444426139083848220937733906099800713143717781531720357997525064017534666891239751607362728534227060943030249612485685768763425556037602755941474941653542152331426660114000937927945792774679996260878096042199067846454542970089145187072428525775090059534747944521032702651627813702799121992043508241301343100990505681206225401410616516828509560584294943674190816159053496770183899973772906699936340807573347813190021439569888598278120014838372065540584097789180572530926368399706041302924025347164795845553931546458809162725834773804499135858726966919527182118595973309247686431104142464761572005006162075237286791399988357786744489056499818078934145418782694253886789001319109746758353842874799102113632149356363874644145315009537351241347979405800895562361598091499985868523771090975088895348382804174019555267298482986737",PREC,10);
  19. #elif N==18
  20. mpf_class lgNfac("0.8063410202609044360601363981582939110986429224594765324598505005527467487615277139806510891554237734774029195140722883868032202409465625042289217549152796466499131134001521472305151349812645719621925123026835994543652760333300441335275580125242473245351570708648869390848865811435941492272788771121553834401048714678961429498903573381880173302979808617317692601358186358465568458758548072821670923081450088401401051187945501328962944820859000920559676891702697909686622506585427883138527055962571578613946865551695231733240457398757515231790525043199275090891827860163442479337966214899386335463425427228018200215335032179019118157358200998303565083725047161238748324796599804490368100804538211517977003313484858590472784853563758898031244080091327145329497840700494039633740824034907267354161416441812254315795816802892388779545646279205656556835935586296071218297712842538543129517008363181292949464962026323889929559590152926759626220418814002352447664319203262102846332753689177115476210155201585",PREC,10);
  21. mpf_class lgNa1fac("0.8063410202609045038934903797764462735174675805069685685804515133250058403716327035570427130955833556988666835830436027631130189996934664551784935763103169754809996132192333595440233328805382249575726542139348884589650801121769157595972661666683934757977803939536480208350053760000001843362048578849457328688196761728920740175616819043814464154649927736236486882658116880477735517155724375974398317721462744365383336410983910842211967082240027994360091361300335138448463461377630924505255045425073447100237495138444508449716416040619922671382969252613399888466709557647948642336238841217295122913830803695016809704209418578023017076133109437843422646980792487107934634929053697419046728791934655666364442441987137102643239892387101635211921249530496197022023106846226760898212531758694177241548960167525681484132304747568944318778358824109890806727213047564953750308328278180736735618844557778822076620026163921353458182446818005707879515262274049492853329650007312470274199494764564169818213493792784",PREC,10);
  22. #elif N==19
  23. mpf_class lgNfac("0.850946212137333975964698739152232290497722598539741313666693692604975901121704368750355100556296757361027674434366966254232695409553455460266582675428869914901788449345745312405367296319968315301755458488231441913562629865478936247705580959544455805900395715382213953051315219874753989127629876010525986493352917760143107634347329232153136968704005032429683472476441737279430492864799897795421161103213111473027049755441995046895226977637239548704389930860493142953939222255825750125024331612127445734972760487230900554774279206065609336164912792207562521921239080623868310552493760140533286470371507126924693208875992718787157791189617033373710137143314727609236238021956995679939251060307241823914231119900032391106105110220784245009795447192405122788054763703283817434933612341252442070872044427165258901230672760014360685856995672268660367727305603634683999334284035610283323291409916227886668235038196216294120445299666276693725670018678583631147313926708913141583311503508563545202440476109853",PREC,10);
  24. mpf_class lgNa1fac("0.850946212137334011666463992635473017402022809256750220737861524054429426162957988448086157530142948743240175214089472758943864734211757559115942781584924627189602119767804517487908707188706866080996391202445104629289279116780140621800543064392583297389498524127554633491432612521879808951797395192773274626108121637922373967633150528453834243397052728406531400140933818626750411877448954615069796431410071089448644430916519085552547279013177482992247692193340541697318928258921090801558411359441072061296119022333336106441771688203560867153984955737185433184231163076835813224488354949874194526196835468058718901108089685414327129155434322515165593723807471477594063056542970829434397835657297407073741221224631357775615305590561173886526541598978688679660006761838421748232726397296559474365752013910730711858376006742372605072590736670392430375878077408090879964545098504783583010168341301522609626203382232438247050908091114547388778128713452523527236708133296030169824838113168989754469585938634",PREC,10);
  25. #elif N==20
  26. mpf_class lgNfac("0.3861246168777145928102087686397162558179621413160826726770968303876057793865949463619627621738158477767872446348676920045180374220888690520263506012475825665546818091887679339023564639435997750313594487470017268085107025184098226711244427954746848914396520340786216586177261434663212462310456602850349248458780852836456283117856600621802606146561694212937583548471922615434027638715119546571682386026121402654936564745413711245500003742637917753883515504146779826373943251460924120970746580567070421356469832933827986691670371425160436948308827757489914200186162228664145772984935393255920669063410337520731025370178383237592979110346471650302763152194633997463021122209140196337474945744236959245234340178868940896962351520092968213776280845754756735561081152982066425933300422644336756226480183863342004389797307005395598025181220368212150969847750646602958606219069496767131007397847895728352638152783853757158138909379122806238066747426972583376548051143725715090473866360419257548956557375744010",PREC,10);
  27. mpf_class lgNa1fac("0.3861246168777145929887175949071324601495252181672115152818198022008222004395554496309888961514840651227181918292491239971646978829541678280147882487908780519279071159475925601717825671007693588388143793740872039767949076953442499405918578293298110342634026529426519809897818518540372898789685500888436524258561912356229588508388910290403332071475486729438676877746021434896113867529820767474291007966504510514854509030221494294355526144413841024937376481900151941868611432731369746195119962537514045556969381524416626882450179542296909561384449333825480232271476902929035217283337791738458055071505728019525696619451689115463219560152158237202423875997258403756554282079477181234954814166089359805628686258773581211437487821467498563112263652129216427051475641881981439358069136383153772058116033738477452281314742904068044649590404266651099230974623905000418365722049906679343546253507263642664440750568060705944203780888807454343315673233624312214812497009352371806258760301258405700042849746787873",PREC,10);
  28. #elif N==21
  29. mpf_class lgNfac("0.7083439116116338608174529304874677585016634018307439440105974244110883720723352354710426456689207226286762474784740834843226127349977153284024918911075540534857278346960033858534919736816032091793173639668921642499371024129713089418241166945509966731920360611112857395408199647397863462054936410754157134707462506214219181775593226522157584953132423474614011990453527772041106584128779909504715378125711095182843033184057619197025025206546830337151755395843535849283219009324341208759337780975150232909948694363363967359462353257034824462208241017340123464884088871757918507158292547722369349608462269403568496047003974795523276990428779577937557310933803703396554978189510040080328063295413897417845672139242983437251326902053889349259633769531650585826417278902702522323669801299789441323786485885290851827784269485137072448966779430177264844394328955994303930618867310922488237581705315771239385502821052336539379907746809947454172477031302224594564413246938374575554812437427976064806997301229909",PREC,10);
  30. mpf_class lgNa1fac("0.7083439116116338608259533507859161491857825686978583113858082562894771842694836106144295023694352613953568719982557373715769482855184894123049656278059451700755419465432019135764935761452745598698692628163697004837606217507752347004727136544150519480403290235627216350774073737729120657117034295982600646069141082824021863364384164515037066967361140539961674558920516945007955304175027616909413795100480148289437536915139702760456965964329389619689088068726496274873647591615808193916982786378979663662216085230280087689963415921005592928489219546292801346003866631820753400792104234146677593860411077439615396728215099026211737381790297532012435252625929536710804843607711149349146326998354436540304405945616868365818824500903879022893545380764256090451925500696644134737957733688230063512352696930397653095991304905401991590327051566157673057069663901221161000331836050871724441037602485330827469241641241413647597303090226097426607660522460713609997993069056846075785954268441652706331506449457059",PREC,10);
  31. #elif N==22
  32. mpf_class lgNfac("0.507665924338400967813917964549850269765554737593055799175639042179495935800090925292646841728380736728355376621212604575899591358462462701707735680340424076699709295187262253528985370057082031953548524476389094254484499263117036288438467239951510111717696096130680780689683707209679146489302795442511606235665306316697654628285694800713483382689265486107718460897728244117898012842903535136153814736243891549920506582720053596317978001531737968475837102251259335025949858233179566725292890908179781684441286179607913979092354373268856899555402152480433332938499727966311535634872088166928447865962537690074314169931621809977201268830951540066876170046342968613261371584684424718145276566711345675171430013566252595321528254841081817163187633757021162008038837121622899393153385880371094940035803939348271983419278292080114874727426406187530360188611886267262917295420385134912817184976696568544175168922688236087203876738000296127882074349236394608928948045945321680971235389925457681829548892671801",PREC,10);
  33. mpf_class lgNa1fac("0.507665924338400967817781791958235901894735907891494429875408695014005781253015777179339903334583566457910780500105476340342410046865492605963814304185530501759933054400328474759717017051100416306028174609002347179513061103906624525417262859914969207191661521051690789571731927136198741417089817894074987107433590362415957359421919234231606505362080973212162133106050397788093588784238701078505170746074242079348252415942488065376003407016218519402362625222687829035020963435602589454319249317257713335726029321489623861278773197543680983203580527922117945983225111269479875949207406567531028859495369518488677711709380119749482493357870009476270512136692838602767113141883334252320267117159918133595781490791964739441067736840331075881288406924426485475309899007220711654817719878666015284037384197242327205891903124654352881641425055238355741498348857324982011029246925864764726016429946438942147510417443254214072137212402498131766929072868896846711142540231782470725502647047923276143034634464661",PREC,10);
  34. #elif N==23
  35. mpf_class lgNfac("0.4124944284514329756491689087061745766740665080954017681931693908326382299768155652028488750000528792674493491487376427652086925627989968590428446272504291013606166847476553342673170658155470290257835339035602696746170669074228839086674247552085763293961699294187977670301544015623976211152241762364623583602876810671583095731046641493364920578558378534727342334125463142657099843917919296316158989154785204781112584862608412434319402855498124503008017387348688202582954407809467654015007317990880652748150508476880007613782167346183269395854697255222027975286740577753359281454207323581643318083979851596888052795254347937575177373688382291887092926241593426884147327919315541475361690091379203866016257710697697535484850403257479159786558505847912263856263850257995189455368479522139367775721966300282858285715545510614066451919056253888410033832896816054605273693889760167892942935196429417828141871898376019843078016117444886267196632351955662528673184527995416828356923646255027630496466189711366",PREC,10);
  36. mpf_class lgNa1fac("0.4124944284514329756491857079557762533354977377230193788226556774174234726954069799941999603021127430377707619774282352698134559308882770618523427888451851006055215366145341734384306598628055724466545007131151442383432465464216843057555595184734331622658403913798367595093333513986329783827998482119032946522903665646474511483605387510677869317172067727263776117557989754814193199168300276190561119334458525367732870582018060431135898247407805298728545872686054648434464084172267831958629695094536861509754262921665601258105543859050054808102684785781674389909905033979193749759792159167468276359490975525971289518257526895912190325689623285309928465824003857969588297055259005047594934604337790713569407147643430933212732533910498371480759075873547066983096860387962045539608579691572087147621049012704481364140282978935513598691699972662878016308649259507240141176675039355559320093700335980135403584860924895030151249691774192086667760961928853021635978047314777572071685783192330726274307246848767",PREC,10);
  37. #elif N==24
  38. mpf_class lgNfac("0.7927056701630389985854134961347689661787650166724232569543174145191919505837502167065989870022276968999674734912826747960559653484247947283322653344360452364873685293634068997425614730458326587079877996669248389953127223636983093062293394417636878795577296381658194101563509159617466799709652327581647525747993368980825580407360750193847923555319502877730372527838055588547741835720938163242358185464639599695128971143797486815291648959939806367812999434929017634016018836201952046005648595324164123388493687987787681243521995394105887445754940133248247021664623790610337539246602644525181839646621683566269463569752575140539668077092692659370556260945604070768044769925446087915279723672411786533062924390011187347849142612134256336095501190980303324339120481802440851201829930416311083147782271567911627342770304410061844917122858753300052932921240997357625805297306361621820005466692290593365223045739865106608431558618634093469632129910793399714846675743012582213251677662540752061790064675973892",PREC,10);
  39. mpf_class lgNa1fac("0.7927056701630389985854141961035023693729913179198812319853919857741659525656538103466421884778803171955659891418020464350044179071771563933350532269386528540990162162476491293491011026248996367790258219577182163054884572572800634839527435089641676993064474465628702366629118339112380171664931944983396691527816453687750249305025615938976194014154569501208452374891867751119126575955251021314384625925857895916079375182154831072300292302338067102413291938161422625728385615970433044012808082250769446269221540477369911430869948600926636421513749728118042740092673928597173751167336669821635711558289873758805003826878809305693075188998559501756264912335233477496050781464715270987387883044415543894242613820340835894441873848249904117961225759444081080607894622657599230743160928607657162714828649814073237715617644580066033684531187177303465984488276456238993888124996470513936887022412294758854933191923614934071064133210494542109405225493410557871649685153706688028414612970657473251909807419751801",PREC,10);
  40. #elif N==25
  41. mpf_class lgNfac("0.1906456788350766081579357066857829126423852537482061743334624922649755720349011977327444827658553528185985191084206840378664295861577477163328806670266540863583626008550200944189220044226267717056199938705676737610038432999744512135215700427232092578585047130850188835311616730040549853343998873902001001817137498828199229440342207414548985199604124516714572375847093832238547544020298865689835735618823017331309941163854054418082095429938449957454748288356444267176010777791755304314204097414278172145499543094593508969729810955916232221467110202683543665134777494529782614381719378294407074460544022778656799247147794102928025438778983425512450230842965531060475001551079686600208334304552351690422706272073370336136649792389888398562530393855600098793067703244875634205096309810142454836565992695558136365637035919299370238474409361413071728680350911421076591527735439308124637253816331592433283210248550024880394630459721034380949975094205400224045201308978978315470567948719364054281830876705577",PREC,10);
  42. mpf_class lgNa1fac("0.1906456788350766081579357346845322487701543057981261541201048701359421563799936706258242360769176343922873563690126396821757048655784298176212838543278503553384355945133862658318608082925051144833934695962993247907920952034962242885917190366331163206594630191482193331232191352133804753818522127066594416157236310803435814872087892687685445990528298002125134931474062684297650648235558201355262528778009926529223835320633013999706845439750910807290406686991861280862277774800507521136533601745058070343537269078708938401870147183144977328084883810000916680735439389032227945185758620957227339716122581839796044400768014321501319685083853565584827686863710577506268567582798689153110858164925948595244049067203343993618934408388359134393375008311961720009784946418993624172055247559408976248620775250690930036178925287384177981977689851798500025522803074174250000060737348891989952387439757080007018711900066940769132134620935703857444233940694460890562723444867342449708829506824720011806099330795229",PREC,10);
  43. #elif N==26
  44. mpf_class lgNfac("0.6056190268058945725781797593526043702403043235973829874450809283896368529395199480273191891429411845322609517865808229481981575115028493782486719780631244960911887277326457506547257997976702049640017503817799948939374953830269654635944954054540570260266783277809043576876195434850683191953114536011390875496029145790526130084620248835477742920619708190413476299442626525473007645491146792687923808381421165093255185131719883881693843827390604530218009188253050398280592674243678932776453870627664911597222642708012958900734888954637478535678687939385890944643306827424625043463535151885225511635617238524598034871409944859145696591790973173455201823726660850796858545484718174884251920852348703292346234204955886032468621023789996927752657063460157138914921816028798676669364015495806440549492269846033564173958132147270601180604109992514089613073625030359614458785057162920828340618135800076218470917160650095845319484617554996776376325830506939515011317933439217117990858784110155251388258807115743",PREC,10);
  45. mpf_class lgNa1fac("0.6056190268058945725781797604294793447067569794454568661991797837068830523501809486750049667381416308612062397832069284518581999590893874561107042618789326476190551868116712727519086117898934498803910927595571068545412635561726109378294290413446591463426019381085551423956193856322055309077486808509977028445704330735418106911574275165790645669076711896733592280695616563098657496562177459543958249452873503738961611438256198238494050860561675125233300523940981565825438383874402301306553632095593219864448907679378390973686425724110106070744263571988352403296904702014604571816991820861053165305773893252511372191136003675100955244938187929271615054843923286655428825328748470468534717546482628056851392387476284068703037260353068106041339229495530526325418722066592419680474842087185731766479009693997286103705658882397551973703308936607142957092826387166000180460145448154094144476384939344744641609579151587103701494773971991828592543332404274796145346960708188282761695544363887932186405118003832",PREC,10);
  46. #elif N==27
  47. mpf_class lgNfac("0.369827909648818844632634691179502978406909161694705819346778493053243112905033171562242560857900890636550300913369606288870649381785314321197030962777127258909175832921598231240814126841006024999594215882664593016345059650958802390952771694372378788648452911583669172528574933735028698979905690104053344242527405031435004761979132430081529250183878594851466196621075959753559308247326554480170373004809729198118709065541661237167591255709546278020825171417458541799749696575246771136877462033028542924778529221362204557889543126551884176084181885923382979392648133633062654948746416677354632986193260888485277743183111798786776829592212723544114692376819353764486913818461008286184769440078994501605252712176178916855279961570672737781116531796165805406254207153092185423447075450562759265199930723329161948022700637182800522297495227257602891254652187554204591632237816870980371254681568098635985178974199488365213935406013852484613121832374153364925156625339495732660127120771082511482002169195902",PREC,10);
  48. mpf_class lgNa1fac("0.369827909648818844632634691578345561542632367564362811954096273683980352522332653716764149115865322484661283514180728713564589965717876229398202718366938737293665438095595609608509918781922423580876068576798688330689286284982075375195792874057642478897380481594126461728026219768290453382268724837917701123110818129898842423913609211208316986377280389367278993012924758647905366068707507801523927317065199501532875930675281051897564853046089144880485906989358749336684923315805897695568936430221799130614659999171647145219345708463765154225586367958110285907218999940487057128146692690509491006171780459820630015768760804419996849333144913019173298805447079057854667668392718622277484193465318637401503145319354994312155283705216386322841446878696411898681959642665438697939594007321497587523692722735848687626019300380288752536079898889383091193713050545403512109169472508529185920678842806008228225991542273846374329744777428688679935875454224500434425950393805737623679116104962767135144918745887",PREC,10);
  49. #elif N==28
  50. mpf_class lgNfac("0.4841408223071011056029575171595725448606430754176530710591677252777941297414315021961902881696510064864482945495532969733182068711291973692048853474755459530177265854546107229090711867498331244716325855356852405081374480626915863440024598795596346647037319221460106716127279076418481476081108566989953598173232175481521296160948006578200850969281929416064364824867493061242981991109566294368310293409086482721556366172855075800744250440180168022379710881965318550375981172071671321163438629882539761945399573372975942280422158005977791700802927094159131589037307304071244443491916415777201704933797512770972998452022428747781690602054483368340710591668574765083047748728616308876798326829871937113770463288081274174061195184015736540787953763252426659053865647923195662559189054727265299768486624658532967771788075444382683177500375269528526877315111530483902418852546335356648959439233024468315500943898196685301193711669794434219988740919336056123771713598618999577299093466047546614809428106590062",PREC,10);
  51. mpf_class lgNa1fac("0.4841408223071011056029575171609969826575563725814732746042569321843043051433003920070292383271026572267168308739063819122485559028078769532398027311019667157872401865581440394690318427457087676479026712032962413793913772266620353925252468433493715588675278318633095707067605327668904753638168384297807662535324091423874660660977624023684706169134012040881161025186369461892041270179143925469585432798849459853515214097124295473795933487085351495732036564604560193754969209694540495487113274218719184889276134147189412812232184536559762076211613458662758399836811195791673038225462151016791945503711676373855090745814478193599815600760567779904023154744295779828264748880616248580390509491536702078891786710436729906028649260833295738640380185747654135268144843127200328037241915316203277225416283766476433512440734681170876328804502154885090451375361577180300966114202422405861138631103588191041221039530997230382922421019087575660958990448385802992033491943972456676938156809946900625893069485033086",PREC,10);
  52. #elif N==29
  53. mpf_class lgNfac("0.9465388202060571929358042801288275361149374933048069121245646329303800617009681750206329726120811354999988382281411826129147561793116155732499125301548691071329097656467050813413582891123401963742807816957810993291570486303584578581073647622406771260327287826210478936496658926770510862435779635717066586073357940491680936002471190265229749315563175416748606483198438197011901765779832218943086665302656434645373152365745773678548698504268224473363383923250678797518553128721210385213163380213411426327664772799974959444538498600588550101332124588018545736714499539121581872708128818774747201852208191822912848686117406357376338966389962576441070527688662121851814243202906803862760755210798397615413325784260362403221222816909442065387744255351667318087369139211979697695360599598075190100263661651666724191009474601247163466318121610078187234726543710140292007499071073318313968225002843193600815282340085307275664294250986125793444439919696692904922183796287872538499749465587650183478105914686309",PREC,10);
  54. mpf_class lgNa1fac("0.9465388202060571929358042801288766546596586414828696777640505454253427022315073893460997909680634441790657520708865811702540947044771666735732608947597439706995252893691550579598932253535910005316408534231151677810850573214923664637171999339225959102737089611516335274192061793214616602233300069880298784806798585263546353449354666761390469363348719909421576782559806502186752681171752666315221815123528308372627686615386695332661139136200729451848176065960172563569016866081863099805594048592057326062188869777029389128812047189078819000285685287860810875104082234131310311079347880853167727687941235140803235341446254890377034394151628160043494131067877033867770319898625603697914884417425621901665209683210114240000048418809176892667965623900215818681254010523483242328629369568969431731926088378617610464085935988119518528800241223464103501108170735372422041947511940982829624528703201740552930013490298905429486116933686837647033191613439634107509021479619500666863752618355513301191350902692464",PREC,10);
  55. #elif N==30
  56. mpf_class lgNfac("0.4236600749257196302308321833839428453150663574955027769544302732356092144846292980636013282596974370104635309963932285064777253215368429245402562362263985170661527174998764388311434934078169955529333387646099207983893854910480961166076253502350707436454511037468687468380785426398626031444710020414620742322190693571983894228257484796764344758751232218227936448924588008438752320031892139540502186843785956013660993677019699463706614313707871722630989250972148178691605469498399664666637910681865970103516734971091374663590049991226685314800622570197709748297613307857727743203199240372123575635733532610875262976708462003923365712323709093137374817238714956174357032647487748330071704740041828018499665286667126698016775796169667335397230744797003540251146602920077533946721619582993963005499548610765256782364330764551229913549250021659358327453552762538488718448131291301697978437184765867739986702944601771448962444513805744362856705253652430854893463360254632076722838911141292603509353702046362",PREC,10);
  57. mpf_class lgNa1fac("0.4236600749257196302308321833839444825998903957681048691424131370749431837086625007762608451726385977657879260742362062627506401243542058318605684852290201585078302950298347257452193161120352184712831417978879178667771909273476572361102190632116777638247494994003546319408833719770279101888636213387063897269149810532399125615048821288064112123286536087281548857362853056056347188066159588484177367816943016431593332840398926487838082982434673092531957568172322091332820510169253615442587101306128848608222947234179084878407799474289013169638872949752393756101414525521849089812694382642740010506026416829941518805937460373257422623947912720097443907374260197456518949555221794921771834245462883136719655139620953239608732668089523318410390132477906791551331358227234303399097808157226399397288355308223078914045771435370728644699306989107269876865852818115073943647119013021486131166640894838353930739060512718134909055904281567777361459630180219700675026300663454384718757443909278959023463639454132",PREC,10);
  58. #endif
复制代码


代码里面宏CHECK_COUNT会预先检查部分边缘上的候选数是否合法来减少搜索范围。现在的经验是对于越大的目标,需要使用越大的CHECK_COUNT.
比如对于16!已经需要使用5000了。
如果能够调整代码,根据不同位置的取值特征不同动态调整CHECK_COUNT应该会有更好的结果。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 昨天 16:13 | 显示全部楼层
N=17
51727020012214
(使用CHECK_COUNT=50000,共5分钟左右)
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 3 小时前 | 显示全部楼层
上面分析过程关键不等式为\(\frac{LOW}s \le \frac{(p_hq_{h-1}-q_np_{h-1})(u-x^*v)}{q_{h-1}x^*+q_h}\lt \frac{UPP}s\)
我们也可以写成\(LOW \le \frac{(p_hq_{h-1}-q_np_{h-1})(su-x^*sv)}{q_{h-1}x^*+q_h}\lt UPP\)
于是分两种情况,
1. \(\frac{p_hq_{h-1}-q_np_{h-1}}{q_{h-1}x^*+q_h}\gt 0\)时,变为
\(\frac{LOW(q_{h-1}x^*+q_h)}{p_hq_{h-1}-q_np_{h-1}}\le su-x^*sv\lt \frac{UPP(q_{h-1}x^*+q_h)}{p_hq_{h-1}-q_np_{h-1}}\)
2.\(\frac{p_hq_{h-1}-q_np_{h-1}}{q_{h-1}x^*+q_h}\lt 0\)时,变为
\(\frac{UPP(q_{h-1}x^*+q_h)}{p_hq_{h-1}-q_np_{h-1}}\lt su-x^*sv\le \frac{LOW(q_{h-1}x^*+q_h)}{p_hq_{h-1}-q_np_{h-1}}\)
其中我们记左边为\(L_i\),右边为\(R_i\), 那么相当于找最小整数k,使得区间\([\frac{L_i}k,\frac{R_i}k]\) (左右区间可能分别半开半闭)包含\(x-x^*\),其中\(x=\frac {su}{k}=\frac{su}{sv}\)
而其中我们搜索的目标是使得\(q_{h-1}su+q_hsv\)最小,
所以如果我们在尝试\(sv=k\)为某个指定值时,那么必然有整数su满足不等式\(L_i+x^*k \le su \le R_i+x^*k\), 也就是这个区间里面有整数时才是一个解,也就是通常左边的取整小于右边的取整。但是这个区间里面还可能存在多个整数,那么也就是su有多种取值,这时由于目标是\(q_{h-1}su+q_hsv\)最小,我们应该取su最小的一个(因为sv都相同),但是前面代码里面取了最大一个,这是一个潜在小错误。
另外,如果在相同参数\(L_i,R_i,x^*\)下搜索时有多个\(sv=k\)有整数解,那么我们这是总是值需要取最小的一个即可,因为根据上面分析,如果有整数解,我们总会选择\(su=\lceil L_i+x^* k\rceil\), 所以k越小,对应的su只可能越小, 不可能更大,对应的\(q_{h-1}su+q_hsv\)也会更小。所以代码里面对应的循环找到一个整数解后会继续循环也没有必要,只要保证是从小到大搜索,找到一个解即可停下,退出循环。


毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 3 小时前 | 显示全部楼层
另一个问题是现在代码里面的CHECK_COUNT选择没有任何依据,这会严重影响到性能。
实际上是我们已知k可能在一个区域[KL,KU], 其中KU有可能可以是无穷大。
其中每个k对应\(x=\frac uv\)一个可能的候选区间\([x^*+\frac{L_i}k,x^*+\frac{R_i}k] \cap [0,1)\) ,
要么我们依次对部分k判断区间\([L_i+x^*k,R_i+x^*k]\)是否包含整数,
要么对于余下的k, 将它们对应的x连成一片,然后穷举判断1/x的整数部分的所有可能取值,转化为下一级连分数问题。
通常对于x比较小的部分,1/x的整数部分变化更快,不利于转化为连分数处理,我们需要优先对k进行穷举,而对于x较大部分,比如如果大于1/2,那么1/x的整数部分就唯一了,适合转化为下一级连分数问题。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2025-9-13 15:30 , Processed in 0.035440 second(s), 15 queries .

Powered by Discuz! X3.5

© 2001-2025 Discuz! Team.

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