aimisiyou 发表于 2025-6-19 18:17:55

递推求数

这个递推数列怎么求？

wayne 发表于 2025-6-22 22:38:51

整数部分是 759903.
这题最快的方法 可能就是高精度浮点数了
Block[{n=2025},IntegerPart@Nest[#+1/FractionalPart[#]&,N,2*n],n]]

aimisiyou 发表于 2025-6-22 22:40:38

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

怎么算的？

王守恩 发表于 2025-6-23 08:14:47

凑个热闹！
Table) &, N, 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,
——长得慢一点的。
Table) &, N, (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 发表于 2025-6-23 13:47:57

的确4000位计算精度就够了，不过我计算结果位759901
下面是对应定点计算代码
#include <gmpxx.h>
#include <iostream>

#define TM 4000
int main()
{
    mpz_class M=1;
    int i;
    for(i=0;i<TM;i++)M*=10;//generte 10^TM
    mpz_class M2=M*M;
    mpz_class sqrt2("14142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358314132226659275055927557999505011527820605714701095599716059702745345968620147285174186408891986095523292304843087143214508397626036279952514079896872533965463318088296406206152583523950547457502877599617298355752203375318570113543746034084988471603868999706990048150305440277903164542478230684929369186215805784631115966687130130156185689872372352885092648612494977154218334204285686060146824720771435854874155657069677653720226485447015858801620758474922657226002085584466521458398893944370926591800311388246468157082630100594858704003186480342194897278290641045072636881313739855256117322040245091227700226941127573627280495738108967504018369868368450725799364729060762996941380475654823728997180326802474420629269124859052181004459842150591120249441341728531478105803603371077309182869314710171111683916581726889419758716582152128229518488472089694633862891562882765952635140542267653239694617511291602408715510135150455381287560052631468017127402653969470240300517495318862925631385188163478001569369176881852378684052287837629389214300655869568685964595155501644724509836896036887323114389415576651040883914292338113206052433629485317049915771756228549741438999188021762430965206564211827316726257539594717255934637238632261482742622208671155839599926521176252698917540988159348640083457085181472231814204070426509056532333398436457865796796519267292399875366617215982578860263363617827495994219403777753681426217738799194551397231274066898329989895386728822856378697749662519966583525776198939322845344735694794962952168891485492538904755828834526096524096542889394538646625744927556381964410316979833061852019379384940057156333720548068540575867999670121372239475821426306585132217408832382947287617393647467837431960001592188807347857617252211867490424977366929207311096369721608933708661156734585334833295254675851644710757848602463600834449114818587655554286455123314219926311332517970608436559704352856410087918500760361009159465670676883605571740076756905096136719401324935605240185999105062108163597726431380605467010293569971042425105781749531057255934984451126922780344913506637568747760283162829605532422426957534529028838768446429173282770888318087025339852338122749990812371892540726475367850304821591801886167108972869229201197599880703818543332536460211082299279293072871780799888099176741774108983060800326311816427988231171543638696617029999341616148786860180455055539869131151860103863753250045581860448040750241195184305674533683613674597374423988553285179308960373898915173195874134428817842125021916951875593444387396189314549999906107587049090260883517636224749757858858368037457931157339802099986622186949922595913276423619410592100328026149874566599688874067956167391859572888642473463585886864496822386006983352642799056283165613913942557649062065186021647263033362975075697870606606856498160092718709292153132368281356988937097416504474590960537472796524477094099241238710614470543986743647338477454819100872886222149589529591187892149179833981083788278153065562315810360648675873036014502273208829351341387227684176678436905294286984908384557445794095986260742499549168028530773989382960362133539875320509199893607513906444495768456993471276364507163279154701597733548638939423257277540038260274785674172580951416307159597849818009443560379390985590168272154034581581521004936662953448827107292396602321638238266612626830502572781169451035379371568823365932297823192986064679789864092085609558142614363631004615594332550474493975933999125419532300932175304476533964706627611661753518754646209676345587386164880198848497479264045065444896910040794211816925796857563784881498986416854994916357614484047021033989215342377037233353115645944389703653166721949049351882905806307401346862641672470110653463493916407146285567980177933814424045269137066609777638784866238003392324370474115331872531906019165996455381157888");
    mpz_class sqrt2u=sqrt2+1;
    mpz_class u0=sqrt2-M;
    mpz_class u1=u0+1;
    mpz_class s1=1;
    for(i=2;i<=2025;i++){
        mpz_class v0=u0+M2/u0;
        mpz_class v1=u1+M2/u1;
        mpz_class r0=v0%M;
        mpz_class r1=v1%M;
        if(v0-r0!=v1-r1){
                std::cout<<"Overflow in round "<<i<<std::endl;
                return -1;
        }
        s1+=(v0-r0)/M;
        u0=r0;
        u1=r1;
    }
    std::cout<<s1<<std::endl;
    return 0;
}

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

wayne 发表于 2025-6-23 15:03:44

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

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

#include <stdio.h>
#include "flint/arb.h"

int main() {
    // 1. 定义计算精度（以比特为单位）。256位对于此任务足够高。
    slong prec = 12000;

    // 2. 声明并初始化 arb 变量。
    // 'a' 将持有序列 a(n) 的当前值。
    // 'inv_a' 是一个临时变量，用于存储 1/a(n-1)。
    arb_t a,temp_floor, temp_frac, inv_a;
    arb_init(a);
    arb_init(temp_floor);
    arb_init(temp_frac);
    arb_init(inv_a);
    // 3. 设置初始条件 a(1) = sqrt(2)
    // arb_sqrt_ui 计算无符号整数的平方根。
    arb_sqrt_ui(a, 2, prec);

    // 4. 从 n=2 循环到 2025 进行迭代计算
    for (int n = 2; n <= 2025; ++n) {
      // 计算 a(n-1) 的倒数，并存入 inv_a
      // arb_inv(result, source, precision)
      arb_floor(temp_floor, a, prec);
      arb_sub(temp_frac, a, temp_floor, prec);
      arb_inv(temp_frac, temp_frac, prec);
      // 更新 a: a = a + 1/(小数部分)
      arb_add(a, a, temp_frac, prec);
      printf("Value of a(%d): ",n); arb_printn(a, 100, 0); printf("\n");

    }

    // 5. 打印最终结果 a(2025)
    // arb_printd 将结果以十进制形式打印到指定的位数。
    printf("The value of a(2025) is:\n");
    arb_printd(a, prec); // 打印约80位小数
    printf("\n");

    // 6. 清理内存，防止泄漏
    arb_clear(a);
    arb_clear(inv_a);
    arb_clear(temp_floor);
    arb_clear(temp_frac);

    return 0;
}

wayne 发表于 2025-6-23 15:10:01

如果我们拿到 $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$是有理数

王守恩 发表于 2025-6-23 17:02:55

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

Table) &, N, 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 发表于 2025-6-25 12:26:16

本帖最后由 Gongwen0519 于 2025-6-25 17:49 编辑

用maple程序跑起来还是挺快的：
> restart;
># 创建一个函数来计算并输出每一步迭代的结果
calculate_and_display := proc(n)
    local initial_value, result, i, step_results, output_file, file_path;
   
    # 设置足够高的精度
    Digits := 20*n;
   
    # 计算初始值
    initial_value := evalf(sqrt(2), Digits);
    result := initial_value;
   
    # 存储每一步的结果
    step_results := ;
   
    # 执行迭代并记录每一步结果
    for i from 1 to n-1 do
      result := result + 1/(result - floor(result));
      step_results := ;
    end do;
   
    # 输出前10步和后10步的结果
    printf("迭代步骤结果（前10步和后10步）:\n");
    for i from 1 to min(10, n) do
      printf("步骤 %3d: %d\n", i-1, step_results);
    end do;
    if n > 20 then
      printf("...\n");
    end if;
    for i from max(n-9, 11) to n do
      printf("步骤 %3d: %d\n", i-1, step_results);
    end do;
   
    # 构建D盘根目录下的文件路径
    file_path := "D:\\iteration_results_n" || n || ".txt";
   
    # 将完整结果保存到D盘根目录
    output_file := fopen(file_path, WRITE);
    if output_file = 0 then
      printf("错误: 无法打开文件 %s 进行写入。请确保D盘存在且有写入权限。\n", file_path);
      return FAIL;
    end if;
   
    fprintf(output_file, "迭代步骤结果 (n = %d)\n", n);
    fprintf(output_file, "步骤, 结果\n");
    for i from 1 to n do
      fprintf(output_file, "%d, %d\n", i-1, step_results);
    end do;
    fclose(output_file);
   
    printf("完整结果已保存到文件: %s\n", file_path);
   
    # 返回最终结果
    return step_results;
end proc:

# 计算n=2025的结果
final_result := calculate_and_display(2025);

# 输出最终结果
printf("当n=2025时的最终结果: %a\n", final_result);

运行结果：

迭代步骤结果（前10步和后10步）:
步骤   0: 1
步骤   1: 3
步骤   2: 5
步骤   3: 33
步骤   4: 38
步骤   5: 40
步骤   6: 61
步骤   7: 63
步骤   8: 65
步骤   9: 67
...
步骤 2015: 759780
步骤 2016: 759782
步骤 2017: 759784
步骤 2018: 759786
步骤 2019: 759852
步骤 2020: 759855
步骤 2021: 759858
步骤 2022: 759895
步骤 2023: 759898
步骤 2024: 759901
完整结果已保存到文件: D:\iteration_results_n2025.txt

                        final_result := 759901

当n=2025时的最终结果: 759901

结果附件：


页: [1]
查看完整版本: 递推求数