嵌套高几何类型合理总和所需的二元拆分算法

Question

有一些很好的来源在线使用二进制拆分技术实施快速总和。例如， ch。 34，JörgArndtBook，（2011）， Cheng等。 （2007）和 Haible and Papanikolaou（1997） 和。从上一篇文章中，以下注释适用于对这种线性收敛序列（类型1）的评估，

    S = SUM(n=0,+oo,a(n)/b(n)*PROD(k=0,n,p(k)/q(k)))

其中 a（n），b（n），p（n），q（n），q（n）是 o（log n）位的整数。最常使用的 情况是， a（n），b（n），p（n），q（n）是具有整数系数的 n 中的多项式。考虑到具有以下部分序列的以下序列，给定两个索引边界 [n1，n2] s（n1，n2） s

    S(n1,n2) = SUM(n=n1,n2-1,a(n)/b(n)*PROD(k=n1,n,p(k)/q(k))

不是直接计算。相反，整数的乘积

    P = p(n1) ... p(n2-1)
    Q = q(n1) ... q(n2−1)
    B = b(n1) ... b(n2-1)
    T = B Q S

是通过二进制拆分递归计算的，直到 n2 -n1＆lt; 5 直接计算它们时。在 n1 和 n2 的中间，选择一个索引 nm 属于间隔 n1 =＆lt; N＆lt; nm ，计算属于间隔 nm =＆lt的组件 pr，qr，br，tr ; N＆lt; N2 并将这些产品和总和设置

    P = Pl Pr
    Q = Ql Qr,  
    B = Bl Br
    T = Br Qr Tl + Bl Pl Tr

为最终，该算法应用于 n1 = 0 和 n2 = nmax = o（n），以及最终的流量 - 点分裂

    S = T/(B Q)

要执行

。 s 与 n 精确的位的位复杂性是 o（（（log n）^2 m（n））其中 m（n）是两个 n bit编号的乘法的位复杂性。

在上一篇参考文献中找到的一个稍微修改但更复杂的序列（类型2）也可以通过二进制分割来求和。它具有额外的内部理由总和，其中 c（n） and d（n）是 o（log n）位的

    U = SUM(n=0,+oo,a(n)/b(n) * PROD(k=0,n,p(k)/q(k)) * SUM(m=0,n,c(m)/d(m)))

整数这些部分总和

    U(n1,n2) = SUM(n=n1,n2-1,a(n)/b(n) * PROD(k=n1,n,p(k)/q(k)) * SUM(m=n1,n,c(m)/d(m)))

该算法是上述算法的变体。

    P = p(n1) ... p(n2-1)
    Q = q(n1) ... q(n2−1)
    B = b(n1) ... b(n2-1)
    T = B Q S
    D = d(n1) ... d(n2-1)
    C = D (c(n1)/d(n1) + ... + c(n2-1)/d(n2-1))
    V = D B Q U

如果 n2 -n1 =＆lt; 4 这些值直接计算。如果 n2 -n1＆gt; 4 它们是通过二元分裂计算的。在 n1 和 n2 的中间选择一个索引 nm ，属于间隔 n1 =＆lt的Vl ; N＆lt; nm ，计算组件 pr，qr，br，tr，dr，dr，cr，vr，vr 属于间隔 nm =＆lt; N＆lt; N2 并将这些产品和总和

    P = Pl Pr
    Q = Ql Qr
    B = Bl Br
    T = Br Qr Tl + Bl Pl Tr
    D = Dl Dr
    C = Cl Dr + Cr Dl
    V = Dr Br Qr Vl + Dr Cl Bl Pl Tr + Dl Bl Pl Vr

最终设置为 n1 = 0 和 n2 = nmax = o（n）和最终浮点divisions are performed

    S = T / (B Q)
    V = U / (D B Q)

I have programmed both algorithms in Pari-GP and applied them to compute some mathematical constants using Chudnovsky's formula for Pi,

我想向前一步，通过混合二进制拆分算法和列文型序列变换来加速一些系列。为此，我需要找到这些系列稍微扩展的二元分裂关系。

    W = SUM(n=0,+oo,a(n)/b(n) * PROD(k=0,n,p(k)/q(k)) * SUM(m=0,n,c(m)/d(m) * PROD(i=0,m,f(i)/g(i))))

它在内和内部具有额外的理性产品，其中 f（n）和 g（n）是 o（log n）位。这些系列不是高几幅，但它们是嵌套的超几何类型总和。我认为该算法可能来自这些部分系列

    W(n1,n2) = SUM(n=n1,n2-1,a(n)/b(n) * PROD(k=n1,n,p(k)/q(k)) * SUM(m=n1,n,c(m)/d(m) * PROD(i=n1,m,f(i)/g(i))))

，如果有人可以得出产品和总和来键入此类型的系列。

我将留下Pari-GP代码以获取pari-gp代码 。如解释所述，计算1和2的快速线性收敛系列。使用？sumbinsplit 寻求帮助。还有一些针对2型系列的测试示例。您可以解开其中之一，并用于

    precision(-log10(abs(sumbinsplit(~F)[1]/s-1)),ceil(log10(Digits())));

检查它。

    \\             ANSI COLOR CODES
    {
    DGreen  = Dg = "\e[0;32m";
    Brown   = Br = "\e[0;33m";
    DCyan   = Dc = "\e[0;36m";
    Purple  = Pr = "\e[0;35m"; 
    Gray    = Gy = "\e[0;37m";
    Red     = Rd = "\e[0;91m";
    Green   = Gr = "\e[0;92m";
    Yellow  = Yw = "\e[0;93m";
    Blue    = Bl = "\e[0;94m";
    Magenta = Mg = "\e[0;95m";
    Cyan    = Cy = "\e[0;96m";
    Reset        = "\e[0m";
    White   = Wh = "\e[0;97m";
    Black   = Bk = "\e[0;30m";
    }
    
    eps()={ my(e=1.); while(e+1. != 1., e>>=1); e; }
    addhelp(eps,Str(Yw,"\n SMALLEST REAL NUMBER\n",Wh,"\n eps() ",Gr,"returns the minimum positive real number for current precision"));
    
    log10(x) = if(x==0,log(eps()),log(x))/log(10);
    
    Digits(n) =  if(type(n) == "t_INT" && n > 0,default(realprecision,n); precision(1.), precision(1.));
    addhelp(Digits,Str(Yw,"\n DIGITS\n",Wh,"\n Digits(n)",Gr," Sets global precision to",Wh," n",Gr," decimal digits.",Wh," Digits()",Gr," returns current global precision."));
    
    addhelp(BinSplit2,Str(Yw,"\n SERIES BINARY SPLITTING (TYPE 2)\n\n",Wh,"BinSplit(~F,n1,n2)",Gr," for ",Wh,"F = [a(n),b(n),p(n),q(n),c(n),d(n)]",Gr," a vector of ",Br,"t_CLOSUREs",Gr," whose\n components are typically polynomials, computes by binary splitting method sums of type\n\n",Wh,"S2 = sum(n=n1,n2-1,a(n)/b(n)*prod(k=n1,n,p(k)/q(k))*sum(m=n1,n,c(m)/d(m)))\n\n",Gr,"Output: ",Wh," [P,Q,B,T,D,C,V]",Gr," integer valued algorithm computing parameters"));
    
    addhelp(BinSplit1,Str(Yw,"\n SERIES BINARY SPLITTING (TYPE 1)\n\n",Wh,"BinSplit(~F,n1,n2)",Gr," for ",Wh,"F = [a(n),b(n),p(n),q(n)]",Gr," a vector of ",Br,"t_CLOSUREs",Gr," whose components\n are typically polynomials, computes by binary splitting method sums of type\n\n",Wh,"S1 = sum(n=n1,n2-1,a(n)/b(n)*prod(k=n1,n,p(k)/q(k)))\n\n",Gr,"Output: ",Wh,"[P,Q,B,T]",Gr," integer valued algorithm computing parameters"));
    
    BinSplit2(~F, n1, n2) =
    {
    my(  P = 1,  Q = 1,  B = 1,  T = 0,  D = 1,  C = 0,  V = 0,
        LP, LQ, LB, LT, LD, LC, LV, RP, RQ, RB, RT, RD, RC, RV,
        nm, tmp1 = 1, tmp2, tmp3 );
    \\
    \\ F = [a(n),b(n),p(n),q(n),c(n),d(n)]
    \\
    if( n2 - n1 < 5,
    \\
    \\ then
    \\
    for ( j = n1, n2-1, 
    
    LP    = F[3](j); 
    LQ    = F[4](j);
    LB    = F[2](j);
    LD    = F[6](j);
    LC    = F[5](j);
    \\
    tmp2  = LB * LQ;
    tmp3  = LP * F[1](j) * tmp1;
    T     = T * tmp2 + tmp3;
    C     = C * LD + D * LC;
    V     = V * tmp2 * LD + C * tmp3;
    P    *= LP;
    Q    *= LQ;
    B    *= LB;
    D    *= LD;
    tmp1 *= LP * LB;
    ),
    \\
    \\ else
    \\
    nm = (n1 + n2) >> 1;
    \\
    [RP,RQ,RB,RT,RD,RC,RV] = BinSplit2(~F, nm, n2);
    [LP,LQ,LB,LT,LD,LC,LV] = BinSplit2(~F, n1, nm);
    \\
    tmp1  = RB * RQ;
    tmp2  = LB * LP;
    tmp3  = LC * RD;
    \\
    P     = LP * RP;
    Q     = RQ * LQ;
    B     = LB * RB;
    T     = LT * tmp1 + RT * tmp2;
    D     = LD * RD;
    C     = RC * LD + tmp3;
    V     = RD * LV * tmp1 + ( RT * tmp3 + LD * RV ) * tmp2;
    \\
    \\ end if
    );
    return([P,Q,B,T,D,C,V]);
    }
    
    BinSplit1(~F, n1, n2) =
    {
    my(  P = 1,  Q = 1,  B = 1,  T = 0,
        LP, LQ, LB, LT, RP, RQ, RB, RT,
        tmp1 = 1, nm );
    \\
    \\ F = [a(n),b(n),p(n),q(n)]
    \\
    if( n2 - n1 < 5,
    \\
    \\ then
    \\
    for ( j = n1, n2-1, 
    
    LP    = F[3](j); 
    LQ    = F[4](j);
    LB    = F[2](j);
    \\
    T     = T * LB * LQ + LP * F[1](j) * tmp1;
    P    *= LP;
    Q    *= LQ;
    B    *= LB;
    \\
    tmp1 *= LP * LB;
    ),
    \\
    \\ else
    \\
    nm = (n1 + n2) >> 1;
    \\
    [RP,RQ,RB,RT] = BinSplit1(~F, nm, n2);
    [LP,LQ,LB,LT] = BinSplit1(~F, n1, nm);
    \\
    P    = LP * RP;
    Q    = RQ * LQ;
    B    = LB * RB;
    T    = LT * RB * RQ + RT * LB * LP;
    \\
    \\ end if
    );
    return([P,Q,B,T]);
    }
    
    sumbinsplit(~F, n1 = 1, dgs = getlocalprec()) =
    {
    my( n = #F, P, Q, B, T, D, C, V, [a,b] = F[3..4] );
    my( n2 = 1 + ceil(dgs*log(10)/log(abs(pollead(Pol(b(x),x))/pollead(Pol(a(x),x))))) );
    \\
    if ( n > 4, [P, Q, B, T, D, C, V] = BinSplit2(~F,n1,n2); return(1.*([V/D,T]/B/Q)),\
                [P, Q, B, T] = BinSplit1(~F,n1,n2); return(1.*(T/B/Q)));
    }
    
    addhelp(sumbinsplit,Str(Yw,"\n LINEARLY CONVERGENT SERIES BINARY SPLITTING SUMMATION\n\n",Wh,"sumbinsplit( ~F, {n1 = 1}, {dgs = getlocalprec()} )\n\n",Gr,"for either ",Wh,"F = [a(n),b(n),p(n),q(n)] ",Gr,"or",Wh," F = [a(n),b(n),p(n),q(n),c(n),d(n)]",Gr," vectors of ",Br,"t_CLOSUREs",Gr," whose\n components are typically polynomials. It computes sums of type 1 or type 2 by binary splitting method\n\n (See BinSplit1, BinSplit2 help)\n\n",Wh,"n1",Gr," starting index (default 1),",Wh," dgs",Gr," result's floating precision\n\n",Yw,"OUTPUT:",Gr," either",Wh," S1",Gr," series value (Type 1) or ",Wh," [S2, S1]",Gr," series values [Type 2, Type1]"));
    
    /*    TESTINGS     */
    
    /*
    Digits(100000);
    a = n->1;
    b = n->n;
    p = n->n*(n<<1-1);
    q = n->3*(3*n-1)*(3*n-2);
    c = n->1;
    d = n->n*(n<<1-1)<<1;
    s = log(2)*(3*log(2)+Pi/2)/10-Pi^2/60;
    F = [a,b,p,q,c,d];
    */
    
    /*
    Digits(100000);
    s = -Pi*Catalan/2+33/32*zeta(3)+log(2)*Pi^2/24;
    F = [n->1,n->n^2,n->n*(n<<1-1),n->3*(3*n-1)*(3*n-2),n->1,n->n*(n<<1-1)<<1];
    \\ precision(-log10(abs(sumbinsplit(~F)[1]/s-1)),ceil(log10(Digits()))); 
    */
    
    /*
    Digits(10000);
    a = n->1;
    b = n->n<<1+1;
    p = n->n<<1-1;
    q = n->n<<3;
    c = n->1;
    d = n->(n<<1-1)^2;
    s = Pi^3/648;
    F = [a,b,p,q,c,d];
    */
    
    /*
    Digits(10000);
    a = n->-1;
    b = n->n^3;
    p = n->-n;
    q = n->(n<<1-1)<<1;
    c = n->20*n-9;
    d = n->n*(n<<1-1)<<1;
    s = 2*Pi^4/75;
    F = [a,b,p,q,c,d];
    */
    
    /*
    Digits(10000);
    a = n->2;
    b = n->n^2;
    p = n->n;
    q = n->(n<<1-1);
    c = n->1;
    d = n->n<<1-1;
    s = 7*zeta(3)-2*Pi*Catalan;
    F = [a,b,p,q,c,d];
    */
    
    /*
    Digits(10000);
    a = n->1;
    b = n->n^4;
    p = n->n;
    q = n->(n<<1-1)<<1;
    c = n->36*n-17;
    d = n->n*(n<<1-1)<<1;
    s = 14*zeta(5)/9+5/18*Pi^2*zeta(3);
    F = [a,b,p,q,c,d];
    */
    
    /*
    Digits(10000);
    a = n->1;
    b = n->n^2;
    p = n->n;
    q = n->(n<<1-1)<<1;
    c = n->12*n-5;
    d = n->n*(n<<1-1)<<1;
    s = 5*zeta(3)/3;
    F = [a,b,p,q,c,d];
    */       

Answer 1

对嵌套超几何序列 w 和部分总和进行密切分析 w（n1，n2）我发现了二进制分裂关系，

在中写入内部超几何总和W（n1，n2），因为

Z(n1,n2) = SUM(m=n1,n2-1,c(m)/d(m) * PROD(i=n1,m,f(i)/g(i)))

我们的

W(n1,n2) = = SUM(n=n1,n2-1,a(n)/b(n) * PROD(k=n1,n,p(k)/q(k)) * Z(n1,n+1))

算法是前一个算法的更深层变化。现在，我们有8个多项式函数 a（n），b（n），p（n），q（n），c（n），d（n），f（n），f（n），g（n），g（n）。 设置这些产品和总和将通过二进制分解计算

P = p(n1) ... p(n2-1)
Q = q(n1) ... q(n2−1)
B = b(n1) ... b(n2-1)
T = B Q S
D = d(n1) ... d(n2-1)
F = f(n1) ... f(n2-1)
G = g(n1) ... g(n2-1)
C = D G Z
V = D B Q G W  

如果 n2 -n1 =＆lt ； 4 这些值直接计算。如果 n2 -n1＆gt; 4 它们是通过二进制拆分递归计算的。在 n1 和 n2 的中间选择一个索引 nm ，vl，fl，gl 属于间隔 n1 =＆lt; N＆lt; nm ，计算组件 pr，qr，br，tr，dr，dr，cr，vr，fr，gr，gr N＆lt; N2 使用辅助变量和分解，将这些产品和总和设置为总和

P = Pl Pr
Q = Ql Qr
B = Bl Br
T = Br Qr Tl + Bl Pl Tr
D = Dl Dr
F = Fl Fr
G = Gl Gr
C = Cl Dr Gr + Cr Dl Fl
V = Dr Br Qr Gr Vl + Dr Cl Bl Pl Gr Tr + Dl Bl Pl Fl Vr

，这27个大整数产品可以降低至19个。最后，该算法应用于 n1 = 0 和 n2 = nmax = o（n），并执行最终的浮点划分。算法提供了所有3总和，

    S = T / (B Q)
    Z = C / (D G)
    W = V / (B Q D G)

我将编码并测试此算法以补充此答案。应用于某些缓慢收敛系列的许多收敛加速方法（CAM）具有W系列的结构（例如，某些经典的凸轮，例如Salzer's，Gustavson's，sumalt（来自Pari GP -Cohen，Rodriguez，Rodriguez，Zagier- - 型凸轮）。我认为，Binsplit和Levin型序列转换的合并应为该主题提供强大的推动力。我们会看到。