论坛 › 算法交流 › 对反正弦级数系数的一点儿研究

liangbch 发表于 2016-4-14 17:30:06

对反正弦级数系数的一点儿研究

在 http://bbs.emath.ac.cn/forum.php?mod=viewthread&tid=8882，我提到计算正弦函数需要用到反正弦函数的泰勒级数。和其他泰勒级数相

比，反正弦级数的表达式看起来有点复杂。为了实现我的三角函数的任意精度计算，今天特意研究了一下反正弦级数，并在这里分享出来，希望

可以帮助到别人，另外，我自己查找起来一比较方便。


(一）预备知识：
1.阶乘
定义：
   $n! =n(n-1)(n-2) ...2.1$
翻译自：http://mathworld.wolfram.com/Factorial.html

例：$10! = 10xx9xx8xx7xx6xx5xx4xx3xx2xx1$

2.双阶乘
1.定义
$n!! = n.(n-2)...5.3.1$ , 当n是正奇数时
$n!! = n.(n-2)...6.4.2$ , 当n是正奇数时
$n!! = 1$,当n=-1或者0
翻译自：http://mathworld.wolfram.com/DoubleFactorial.html
例:
$10!! =10xx8xx6xx4xx2$
$9!! =9xx7xx5xx3xx1$

3.反双曲正切函数
$arctanh(x)= x + 1/3 x^3 + 1/5 x^5 + 1/7 x^7 + 1/9 x^9 +... $
第n项为 $1/{2n-1} x^(2n-1)$

$ln({1+x}/{1-x}) =2 (x + 1/3 x^3 + 1/5 x^5 + 1/7 x^7 + 1/9 x^9 + ...) $
故 $ln({1+x}/{1-x})= 2arctanh(x)$


（二）反正弦级数

反正弦函数的泰勒级数展开式：
$ arcsin(x)= x + 1/2 {x^3}/3 +{1.3}/{2.4}{x^5}/5 + {1.3.5}/{2.4.6}{x^7}/7 + ...$

写成双阶乘形式
$ arcsin(x)= x + {1!!}/{2!!}{x^3}/3 + {3!!}/{4!!}{x^5}/5 + {5!!}/{6!!}{x^7}/7 + ... $

第n项
${(2n-3)!!}/{(2n-2)!!}{x^(2n-1)}/{2n-1}$

前5项的常系数表示
$ arcsin(x)=x + 1/6 x^3 + 3/40 x^5 + 5/112 x^7 + 35/1152 x^9 + ... $

liangbch 发表于 2016-4-14 18:14:04

当x是偶数时，令x=2y,则arcsin(x)的级数展开式变为
   $ arcsin(x)= 2y + 2^3 {1!!}/{2!!}{y^3}/3 + 2^5 {3!!}/{4!!}{y^5}/5 + 2^7 {5!!}/{6!!}{y^7}/7 + ... $

注：偶数一般指整数，我们这里对“偶数”的概念拓展一下，若x的某种进制表示的最低单元是偶数，我们可将x称为偶数。


第n项变为
$ 2^(2n-1){(2n-3)!!}/{(2n-2)!!}{y^(2n-1)}/{2n-1}$
$ =2^n {(2n-3)!!}/{(n-1)!}{y^(2n-1)}/{2n-1}$

事实上，第n项的分母部分$(n-1)!$可完全约掉，第n项最终可表示为如下形式,其系数的分子部分$p_i$一个正整数，且是4的倍数。分母部分是奇数。
$p_i/{2n-1}y^(2n-1)$

可以证明，分母$(n-1)!$可被完全约掉，我们这里仅仅证明素因子2可被完全约掉。对此问题感兴趣者可自行证明$(n-1)!$的其他素因子也可被完全约掉。
$n!$ 包含 $\lfloor n/2 \rfloor + \lfloor n/4 \rfloor + \lfloor n/8 \rfloor + ...$个因子2，其素因子2的个数一定小于n,当n是2的方幂时，其素因子2的个数达到最大值$n-1$。因分子包含项$2^n$,共n个因子2，而分母$(n-1)!$至多包含$n-2$个因子，故上面的$p_i$一定是4的倍数。

注： $\lfloor n/2 \rfloor $ 表示将n/2向下取整。

我们再来审视一下反正弦级数的第n项的形式$p_i/{2n-1}y^(2n-1)$，它和 $arctanh$的第n项 $1/{2n-1} x^(2n-1)$ 是多么相似

liangbch 发表于 2016-4-14 20:36:45

arcsin(x) 级数的系数求出来了，下面给出前31项的值
$a_1=2/1$
$a_2=4/3$
$a_3=12/5$
$a_4=40/7$
$a_5=140/9$
$a_6=504/11$
$a_7=1848/13$
$a_8=6864/15$
$a_9=25740/17$
$a_10=97240/19$
$a_11=369512/21$
$a_12=1410864/23$
$a_13=5408312/25$
$a_14=20801200/27$
$a_15=80233200/29$
$a_16=310235040/31$
$a_17=1202160780/33$
$a_18=4667212440/35$
$a_19=18150270600/37$
$a_20=70690527600/39$
$a_21=275693057640/41$
$a_22=1076515748880/43$
$a_23=4208197927440/45$
$a_24=16466861455200/47$
$a_25=64495207366200/49$
$a_26=252821212875504/51$
$a_27=991837065896208/53$
$a_28=3893878851296224/55$
$a_29=15297381201520880/57$
$a_30=60134532999082080/59$
$a_31=236529163129722848/61$

以上各项的分子来自于一个已知的数列 http://oeis.org/A028329

liangbch 发表于 2016-4-14 20:39:42

这里给出C语言的源代码，可在VC下编译并运行。
#include <stdlib.h>
#include <stdio.h>
#include <limits.h>

typedef unsigned __int64 UINT64;
typedef unsigned long DWORD;

/*
Get more details fromhttp://bbs.emath.ac.cn/thread-8885-1-1.html
The numerator of a_i = http://oeis.org/A028329

let x=2y, then
arcsin(x)
= 2y + (2^3)(1!!/2!!)(y^3)/3 + 2^5(3!!/4!!)(y^5)/5 + (2^7)(5!!/6!!)(y^7/7) +...
= 2y + (2^2)(1!!/1!)(y^3)/3 + 2^3(3!!/2!)(y^5)/5 + (2^4)(5!!/3!)(y^7/7) +...

a_n
   = (2^n) ((2n-3)!!/(n-1)!)(y^(2n-1))/(2n-1)
*/

void print_arcsin_cofficient( int len)
{
        UINT64 p0,p1;
        UINT64 d,n;
        UINT64 c1,c2,c4,c6;
        UINT64 max_uint64;
        int i;

        max_uint64=ULLONG_MAX;
        c1=1; c2=2;        c4=4; c6=6;
       
        i=1;p0=2;d=1;
        printf("$a_%d=%I64u/%I64u$\n",i,p0,d);
       
        for (i=2;i<=(UINT64)len;i++)
        {
                UINT64 r;
               
                n=i;
                if ( p0 > max_uint64/(c4*n-c6) )//the next step will overflow
                        break;

                p1= p0 * ( c4*n-c6); // p1=p0 * 2 * (2*n-3)=p0*(4*n-6)
                r= p1 % (n-c1);
                if ( r != 0 )
                {
                        printf("Error when i=%u64\n", i); return ;
                }
               
                p1/=(n-c1);//p1=p0*2*(2*n-3)/(n-1);
                d=c2*n-c1;       //d=2n-1
               
                p0=p1;
                printf("$a_%d=%I64u/%I64u$\n",i,p0,d);
        }
}

int main(int argc, char* argv[])
{
        print_arcsin_cofficient(40);
        return 0;
}

落叶 发表于 2016-5-28 19:02:20

本帖最后由 落叶 于 2016-5-28 22:24 编辑

arcsin(x)=x+1/2*x^3/3+1/2*3/4*x^5/5+ 1/2*3/4*5/6*x^7/7+.............   
               
设 Gb 为 246....的最小公倍数   ,Ga 为 357....的最小公倍数
化为:arcsin(x)=x+x*1/2*x^2(1/3+3/4*x^2(1/5+5/6*x^2(1/7......)))
                        
                  =x+(x*x^2*1*Gb/2(Ga/3+x^2*3*Gb/4(Ga/5+x^2*5*Gb/6(Ga/7......)))   )/(Ga*Gb)

最小公倍数的位长小于计算精度时,除法运算开始加速,需要除法有除尽判断功能.
这个公式类似的测试在我的贴子:基于泰勒展开式的高精三角函数实现 - 第3页 - 算法交流 - 数学研发论坛 - Powered by Discuz! http://bbs.emath.ac.cn/thread-8882-3-1.html

arcsin(x)没查到趋近零的细分公式? 难道只能对泰勒公式加速?
                  

落叶 发表于 2016-5-31 20:11:26

本帖最后由 落叶 于 2016-5-31 20:13 编辑

arcsin(x)=x+(x*x^2*1*Gb/2(Ga/3+x^2*3*Gb/4(Ga/5+x^2*5*Gb/6(Ga/7......)))   )/(Ga*Gb)
错了，Gb化简错了，应为：arcsin(x)=x+(x*x^2*1/2(Ga/3+x^2*3/4(Ga/5+x^2*5/6(Ga/7......)))   )/Ga
单纯用这个公式，arcsin(0.5)，需要的级数运算是所需精度的1.5倍多，arcsin(1)那没法算，如果找不到使x趋零的方法，这个方法不能单独实用。
在这里加公倍数Ga没有用，只有副作用。

查看完整版本: 对反正弦级数系数的一点儿研究