递推求数

aimisiyou

这个递推数列怎么求？

wayne

整数部分是 759903.
这题最快的方法 可能就是高精度浮点数了

1. Block[{n=2025},IntegerPart@Nest[#+1/FractionalPart[#]&,N[Sqrt[2],2*n],n]]

复制代码

aimisiyou

wayne 发表于 2025-6-22 22:38
整数部分是 759903.
这题最快的方法 可能就是高精度浮点数了

怎么算的？

王守恩

凑个热闹！

1. Table[Floor@Nest[# + 1/(# - Floor[#]) &, N[Sqrt[2], Floor[(GoldenRatio + 1) n]], n], {n, 2025}]

复制代码

{3, 5, 33, 38, 40, 61, 63, 65, 67, 69, 3311, 3315, 3317, 3803, 3825, 3833, 3836, 3838, 3977, 3979, 3981, 4183, 4204, 4206, 4209, 4211, 4218, 4221, 4231, 4239, 4241, 4244, 4246, 4253, 4259, 4261, 4263, 4505, 4508, 4510, 4512, 4514,
——长得慢一点的。

1. Table[Floor@Nest[# + 1/(# - Floor[#]) &, N[Sqrt[94], (GoldenRatio + 1) n], n], {n, 2025}]

复制代码

{11, 18, 20, 24, 26, 34, 72, 76, 78, 81, 83, 87, 89, 91, 93, 95, 104, 106, 113, 129, 131, 133, 135, 144, 146, 152, 154, 237, 239, 260, 262, 264, 266, 269, 271, 279, 287, 291, 295, 299, 302, 305, 308, 310, 312, 371, 373, 391, 395, 397, 400, 402,  

mathe

的确4000位计算精度就够了，不过我计算结果位759901
下面是对应定点计算代码

1. #include <gmpxx.h>
   
2. #include <iostream>
   
3. 
   
4. #define TM 4000
   
5. int main()
   
6. {
   
7.     mpz_class M=1;
   
8.     int i;
   
9.     for(i=0;i<TM;i++)M*=10;//generte 10^TM
   
10.     mpz_class M2=M*M;
    
11.     mpz_class sqrt2("14142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095523292304843087143214508397626036279952514079896872533965463318088296406206152583523950547457502877599617298355752203375318570113543746034084988471603868999706990048150305440277903164542478230684929369186215805784631115966687130130156185689872372352885092648612494977154218334204285686060146824720771435854874155657069677653720226485447015858801620758474922657226002085584466521458398893944370926591800311388246468157082630100594858704003186480342194897278290641045072636881313739855256117322040245091227700226941127573627280495738108967504018369868368450725799364729060762996941380475654823728997180326802474420629269124859052181004459842150591120249441341728531478105803603371077309182869314710171111683916581726889419758716582152128229518488472089694633862891562882765952635140542267653239694617511291602408715510135150455381287560052631468017127402653969470240300517495318862925631385188163478001569369176881852378684052287837629389214300655869568685964595155501644724509836896036887323114389415576651040883914292338113206052433629485317049915771756228549741438999188021762430965206564211827316726257539594717255934637238632261482742622208671155839599926521176252698917540988159348640083457085181472231814204070426509056532333398436457865796796519267292399875366617215982578860263363617827495994219403777753681426217738799194551397231274066898329989895386728822856378697749662519966583525776198939322845344735694794962952168891485492538904755828834526096524096542889394538646625744927556381964410316979833061852019379384940057156333720548068540575867999670121372239475821426306585132217408832382947287617393647467837431960001592188807347857617252211867490424977366929207311096369721608933708661156734585334833295254675851644710757848602463600834449114818587655554286455123314219926311332517970608436559704352856410087918500760361009159465670676883605571740076756905096136719401324935605240185999105062108163597726431380605467010293569971042425105781749531057255934984451126922780344913506637568747760283162829605532422426957534529028838768446429173282770888318087025339852338122749990812371892540726475367850304821591801886167108972869229201197599880703818543332536460211082299279293072871780799888099176741774108983060800326311816427988231171543638696617029999341616148786860180455055539869131151860103863753250045581860448040750241195184305674533683613674597374423988553285179308960373898915173195874134428817842125021916951875593444387396189314549999906107587049090260883517636224749757858858368037457931157339802099986622186949922595913276423619410592100328026149874566599688874067956167391859572888642473463585886864496822386006983352642799056283165613913942557649062065186021647263033362975075697870606606856498160092718709292153132368281356988937097416504474590960537472796524477094099241238710614470543986743647338477454819100872886222149589529591187892149179833981083788278153065562315810360648675873036014502273208829351341387227684176678436905294286984908384557445794095986260742499549168028530773989382960362133539875320509199893607513906444495768456993471276364507163279154701597733548638939423257277540038260274785674172580951416307159597849818009443560379390985590168272154034581581521004936662953448827107292396602321638238266612626830502572781169451035379371568823365932297823192986064679789864092085609558142614363631004615594332550474493975933999125419532300932175304476533964706627611661753518754646209676345587386164880198848497479264045065444896910040794211816925796857563784881498986416854994916357614484047021033989215342377037233353115645944389703653166721949049351882905806307401346862641672470110653463493916407146285567980177933814424045269137066609777638784866238003392324370474115331872531906019165996455381157888");
    
12.     mpz_class sqrt2u=sqrt2+1;
    
13.     mpz_class u0=sqrt2-M;
    
14.     mpz_class u1=u0+1;
    
15.     mpz_class s1=1;
    
16.     for(i=2;i<=2025;i++){
    
17.         mpz_class v0=u0+M2/u0;
    
18.         mpz_class v1=u1+M2/u1;
    
19.         mpz_class r0=v0%M;
    
20.         mpz_class r1=v1%M;
    
21.         if(v0-r0!=v1-r1){
    
22.                 std::cout<<"Overflow in round "<<i<<std::endl;
    
23.                 return -1;
    
24.         }
    
25.         s1+=(v0-r0)/M;
    
26.         u0=r0;
    
27.         u1=r1;
    
28.     }
    
29.     std::cout<<s1<<std::endl;
    
30.     return 0;
    
31. }
    

复制代码

在计算了一下发现3000位精度也够了，但是2000位就不行了，第1643轮发现精度不够了

wayne

前几天在QQ群里 学到了 flint这个库, https://flintlib.org/, 让大模型 帮我写一个, 然后我稍微改改, 也能用了.
sqrt2直接用 arb_sqrt_ui 产生.

然后,我发现flint 需要我预留12000位的浮点精度才能计算出2025项.

1. #include <stdio.h>
   
2. #include "flint/arb.h"
   
3. 
   
4. int main() {
   
5.     // 1. 定义计算精度（以比特为单位）。256位对于此任务足够高。
   
6.     slong prec = 12000;
   
7. 
   
8.     // 2. 声明并初始化 arb 变量。
   
9.     // 'a' 将持有序列 a(n) 的当前值。
   
10.     // 'inv_a' 是一个临时变量，用于存储 1/a(n-1)。
    
11.     arb_t a,temp_floor, temp_frac, inv_a;
    
12.     arb_init(a);
    
13.     arb_init(temp_floor);
    
14.     arb_init(temp_frac);
    
15.     arb_init(inv_a);
    
16.     // 3. 设置初始条件 a(1) = sqrt(2)
    
17.     // arb_sqrt_ui 计算无符号整数的平方根。
    
18.     arb_sqrt_ui(a, 2, prec);
    
19. 
    
20.     // 4. 从 n=2 循环到 2025 进行迭代计算
    
21.     for (int n = 2; n <= 2025; ++n) {
    
22.         // 计算 a(n-1) 的倒数，并存入 inv_a
    
23.         // arb_inv(result, source, precision)
    
24.         arb_floor(temp_floor, a, prec);
    
25.         arb_sub(temp_frac, a, temp_floor, prec);
    
26.         arb_inv(temp_frac, temp_frac, prec);
    
27.         // 更新 a: a = a + 1/(小数部分)
    
28.         arb_add(a, a, temp_frac, prec);
    
29.         printf("Value of a(%d): ",n); arb_printn(a, 100, 0); printf("\n");
    
30. 
    
31.     }
    
32. 
    
33.     // 5. 打印最终结果 a(2025)
    
34.     // arb_printd 将结果以十进制形式打印到指定的位数。
    
35.     printf("The value of a(2025) is:\n");
    
36.     arb_printd(a, prec); // 打印约80位小数
    
37.     printf("\n");
    
38. 
    
39.     // 6. 清理内存，防止泄漏
    
40.     arb_clear(a);
    
41.     arb_clear(inv_a);
    
42.     arb_clear(temp_floor);
    
43.     arb_clear(temp_frac);
    
44. 
    
45.     return 0;
    
46. }

复制代码

wayne

如果我们拿到 $a+b\sqrt{2}+\frac{1}{a+b\sqrt{2}} = a-\frac{a}{2 b^2-a^2} + (\frac{b}{2 b^2-a^2}+b)\sqrt{2}$的整数部分的表达式的话(a,b是有理数). 这题目是可以存在完全是有理数表示的递推公式,进而高效的计算出a(n)的精确值的.
因为$a(n)=x_n+y_n\sqrt{2},x_n,y_n$是有理数

王守恩

凑个热闹！——这个通项公式可以的。

Table[Floor@Nest[# + 1/(# - Floor[#]) &, N[Sqrt[2], 2 n], n - 1], {n, 2025}]

{1, 3, 5, 33, 38, 40, 61, 63, 65, 67, 69, 3311, 3315, 3317, 3803, 3825, 3833, 3836, 3838, 3977, 3979, 3981, 4183, 4204, 4206, 4209, 4211, 4218, 4221, 4231, 4239, 4241, 4244, 4246, 4253, 4259, 4261, 4263, 4505, 4508, 4510, ..., 759901}

Gongwen0519

用maple程序跑起来还是挺快的：

1. > restart;
   
2. ># 创建一个函数来计算并输出每一步迭代的结果
   
3. calculate_and_display := proc(n)
   
4.     local initial_value, result, i, step_results, output_file, file_path;
   
5.    
   
6.     # 设置足够高的精度
   
7.     Digits := 20*n;
   
8.    
   
9.     # 计算初始值
   
10.     initial_value := evalf(sqrt(2), Digits);
    
11.     result := initial_value;
    
12.    
    
13.     # 存储每一步的结果
    
14.     step_results := [floor(result)];
    
15.    
    
16.     # 执行迭代并记录每一步结果
    
17.     for i from 1 to n-1 do
    
18.         result := result + 1/(result - floor(result));
    
19.         step_results := [op(step_results), floor(result)];
    
20.     end do;
    
21.    
    
22.     # 输出前10步和后10步的结果
    
23.     printf("迭代步骤结果（前10步和后10步）:\n");
    
24.     for i from 1 to min(10, n) do
    
25.         printf("步骤 %3d: %d\n", i-1, step_results[i]);
    
26.     end do;
    
27.     if n > 20 then
    
28.         printf("...\n");
    
29.     end if;
    
30.     for i from max(n-9, 11) to n do
    
31.         printf("步骤 %3d: %d\n", i-1, step_results[i]);
    
32.     end do;
    
33.    
    
34.     # 构建D盘根目录下的文件路径
    
35.     file_path := "D:\\iteration_results_n" || n || ".txt";
    
36.    
    
37.     # 将完整结果保存到D盘根目录
    
38.     output_file := fopen(file_path, WRITE);
    
39.     if output_file = 0 then
    
40.         printf("错误: 无法打开文件 %s 进行写入。请确保D盘存在且有写入权限。\n", file_path);
    
41.         return FAIL;
    
42.     end if;
    
43.    
    
44.     fprintf(output_file, "迭代步骤结果 (n = %d)\n", n);
    
45.     fprintf(output_file, "步骤, 结果\n");
    
46.     for i from 1 to n do
    
47.         fprintf(output_file, "%d, %d\n", i-1, step_results[i]);
    
48.     end do;
    
49.     fclose(output_file);
    
50.    
    
51.     printf("完整结果已保存到文件: %s\n", file_path);
    
52.    
    
53.     # 返回最终结果
    
54.     return step_results[n];
    
55. end proc:
    
56. 
    
57. # 计算n=2025的结果
    
58. final_result := calculate_and_display(2025);
    
59. 
    
60. # 输出最终结果
    
61. printf("当n=2025时的最终结果: %a\n", final_result);

复制代码

运行结果：

1. 迭代步骤结果（前10步和后10步）:
   
2. 步骤   0: 1
   
3. 步骤   1: 3
   
4. 步骤   2: 5
   
5. 步骤   3: 33
   
6. 步骤   4: 38
   
7. 步骤   5: 40
   
8. 步骤   6: 61
   
9. 步骤   7: 63
   
10. 步骤   8: 65
    
11. 步骤   9: 67
    
12. ...
    
13. 步骤 2015: 759780
    
14. 步骤 2016: 759782
    
15. 步骤 2017: 759784
    
16. 步骤 2018: 759786
    
17. 步骤 2019: 759852
    
18. 步骤 2020: 759855
    
19. 步骤 2021: 759858
    
20. 步骤 2022: 759895
    
21. 步骤 2023: 759898
    
22. 步骤 2024: 759901
    
23. 完整结果已保存到文件: D:\iteration_results_n2025.txt
    
24. 
    
25.                         final_result := 759901
    
26. 
    
27. 当n=2025时的最终结果: 759901

复制代码

结果附件：
iteration_results_n2025.txt (26.21 KB, 下载次数: 3)

昨天 12:25 上传

点击文件名下载附件