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

liangbch

在 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

当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

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

这里给出C语言的源代码，可在VC下编译并运行。

1. #include <stdlib.h>
   
2. #include <stdio.h>
   
3. #include <limits.h>
   
4. 
   
5. typedef unsigned __int64 UINT64;
   
6. typedef unsigned long DWORD;
   
7. 
   
8. /*
   
9.   Get more details from  http://bbs.emath.ac.cn/thread-8885-1-1.html
   
10.   The numerator of a_i = http://oeis.org/A028329
    
11. 
    
12.   let x=2y, then
    
13. arcsin(x)
    
14.   = 2y + (2^3)(1!!/2!!)(y^3)/3 + 2^5(3!!/4!!)(y^5)/5 + (2^7)(5!!/6!!)(y^7/7) +...
    
15.   = 2y + (2^2)(1!!/1!)(y^3)/3 + 2^3(3!!/2!)(y^5)/5 + (2^4)(5!!/3!)(y^7/7) +...
    
16. 
    
17. a_n
    
18.    = (2^n) ((2n-3)!!/(n-1)!)(y^(2n-1))/(2n-1)
    
19. */
    
20. 
    
21. void print_arcsin_cofficient( int len)
    
22. {
    
23.         UINT64 p0,p1;
    
24.         UINT64 d,n;
    
25.         UINT64 c1,c2,c4,c6;
    
26.         UINT64 max_uint64;
    
27.         int i;
    
28. 
    
29.         max_uint64=ULLONG_MAX;
    
30.         c1=1; c2=2;        c4=4; c6=6;
    
31.        
    
32.         i=1;p0=2;d=1;
    
33.         printf("$a_%d=%I64u/%I64u$\n",i,p0,d);
    
34.        
    
35.         for (i=2;i<=(UINT64)len;i++)
    
36.         {
    
37.                 UINT64 r;
    
38.                
    
39.                 n=i;
    
40.                 if ( p0 > max_uint64/(c4*n-c6) )  //the next step will overflow
    
41.                         break;
    
42. 
    
43.                 p1= p0 * ( c4*n-c6); // p1=p0 * 2 * (2*n-3)=p0*(4*n-6)
    
44.                 r= p1 % (n-c1);
    
45.                 if ( r != 0 )
    
46.                 {
    
47.                         printf("Error when i=%u64\n", i); return ;
    
48.                 }
    
49.                
    
50.                 p1/=(n-c1);  //p1=p0*2*(2*n-3)/(n-1);
    
51.                 d=c2*n-c1;         //d=2n-1
    
52.                
    
53.                 p0=p1;
    
54.                 printf("$a_%d=%I64u/%I64u$\n",i,p0,d);
    
55.         }
    
56. }
    
57. 
    
58. int main(int argc, char* argv[])
    
59. {
    
60.         print_arcsin_cofficient(40);
    
61.         return 0;
    
62. }
    
63. 
    

复制代码

落叶

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 为 2  4  6  ....的最小公倍数   ,Ga 为 3  5  7  ....的最小公倍数
化为: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)没查到趋近零的细分公式? 难道只能对泰勒公式加速?
                  

落叶

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没有用，只有副作用。