回到十六世纪，让我用初等数学方法编算三角正弦函数表

 

 

                      一   前     言

 

    古人是怎样编算常用对数表的？古人是怎样编算三角函数表的？古人是怎样计算圆周率π的？这几个问题一直绕在我脑子里，总想以古人的思维和方法，自己动手来编算一下。

常用对数表总算编出来了，要问，是否摸到古人的边？还是不知道。但总算了了一桩心事。

    用割圆法计算圆周率π，也已经完工。是否就是古人之方法？也是不知道。但也总算了了一桩心事。

    剩下一愿，编三角函数表，能实现吗？说到编三角函数表，都说用高等数学，把三角函数展成级数就可以算了。这当然不错。但三角学之初是这样编算的吗？把三角函数展为级数，是十八世纪欧拉时代的事。而最初的三角函数表，在十六世纪早就有了。

    电脑《百科》上有说法：

    “欧洲的「文艺复兴时期」，﹝14世纪-16世纪﹞伟大的天文学家哥白尼﹝1473-1543﹞提倡地动学说，他的学生利提克斯（1514--1574）见到当时天文观测日益精密，认为推算更精确的三角函数值表刻不容缓。于是他定圆的半径为1015，以制作每隔10"的正弦、正切及正割值表。当时还没有对数，更没有计算器。全靠笔算，任务十分繁重。利提克斯和他的助手们以坚毅不拔的意志，勤奋工作达12年之久，遗憾的是，他生前没能完成这项工作，直到1596年，才由他的学生鄂图﹝1550-1605﹞完成并公布于世，1613年海得堡的彼提克斯﹝1561-1613﹞又修订了利提克斯的三角函数表，重新再版。后来英国数学家纳皮尔发现了对数，这就大大地简化了三角计算，为进一步造出更精确的三角函数表创造了条件。”

    “文艺复兴后期，法国数学家韦达成为三角公式的集大成者．他的《应用于三角形的数学定律》（1579年）是较早系统论述平面和球面三角学的专著之一．其中第一部分列出６种三角函数表，有些以分和度为间隔。给出精确到５位和１０位小数的三角函数值。第二部分给出造表的方法，解释了三角形中诸三角线量值关系的运算公式．除汇总前人的成果外，还补充了自己发现的新公式．如正切定律、和差化积公式等等．他将这些公式列在一个总表中，使得任意给出某些已知量后，可以从表中得出未知量的值。”

 

    我想，当时只能是用初等数学方法来编算的，当然也可能不太完备，而且编算细节也不太清楚了。我想带一只仅能作加减乘除开方的计算器，回到十六世纪，步利提克斯、韦达们的后尘，试编一下三角函数表。

 

二   编 算 的 基 础

 

    我以为：由直角三角形或单位圆定义的三角函数，最初只能通过勾股定理，解直角三角形来得到三角函数值的。这样得到的函数值是原始的、精确的，但为数极少只有六个，它们的正弦值为： 

Sin 0°=0                                 Sin 30°=1/2=0.5 

Sin 45°= (√2 )/2 =0.707106781           Sin 60°= (√3)/2=0.866025403  

Sin 72°= √(10+2*√5 ) /4=0.951056516    Sin 90° =  1  

    从这几个精确值出发，再利用以下公式

1  二倍角公式    sin 2α=2 sinα cosα

2  和差公式    sin (α±β) = sinα cosβ ± cosα sinβ

3  半角公式sinα/2=√(( 1- cos α) / 2 )

就可以推算其他角度的Sin，但它们是属于派生的，列一些于下：

Sin 3°=0.052335956       Sin 9°=0.156434465     Sin 15°=0.258819045      

Sin 18°=0.309016994      Sin 27°=0.453990499     Sin 36°=0.587785252        

Sin 45°=0.707106781      Sin 54°=0.809016994     Sin 63°=0.891006524 

Sin 75°=0.965925826      Sin 81°=0.987688340 

    当然还有很多，并请注意，它们都是3的倍数。这样看来，要编3度间隔的正弦表，立马可得。但3度间隔的正弦表，又有多少用途呢？之所以列这些角的正弦，是因为在造1度间隔的正弦表时，要用它们作控制、检验的。

 

 

三   编算的思路：先求Sin1°再“滚雪球”求所有整度的正弦

 

    如果通过和差、倍角、半角，能凑出Sin1°的值，那就太好了。因为有了Sin1°，就可用sin(α+β)=sinαcosβ+cosαsinβ公式算出全部整度的正弦。设β始终为1，即β=1，则从α=0开始，使(α+β)成为新的α，即α_i+1=(α_i+1°)，一步一步，α便为1、2、3、4、…90，它的正弦也就有了。这种方法我称它为“滚雪球”法，越滚越大。所以如果有了Sin1°，那么编算1度间隔的正弦函数表，就可以实现。尽管太简单了一些，但必竟可以只用加、减、乘、除、开平方，就能够编算出来的。1°间隔的正弦函数表编出来后，则10’、1’、10” 间隔的正弦函数表也就可以编算了。

    一句话，编算三角正弦函数表的关键，是先要算出sin 1°。

 

 

          四   用‘迭代法’求起算数sin 1°

 

    单靠和差、倍角、半角这三个公式是凑不出Sin1°来的。后来我就在‘三倍角公式’上想办法。靠过去的‘测量平差’经验，我想起了‘迭代法’。迭代法在我们‘测量平差’时用到，在‘开平方’时也能应用，且收敛也较快。这样就可以用sin 3°反求出sin 1°。解决问题了，乃大喜。

    三倍角公式是  sin3α=3sinα－4sin³α ，这样便得到：

sin3α= sinα(3-4 sin²α)，把Sin 3°=0.052335956作为已知数，sinα作为未知数，反过来便有：

Sinα= sin 3α/ (3-4 sin²α)，由于α=1°，sin 3°=0.052335956，所以

Sin 1°= sin 3° / (3-4 sin²1°)= 0.052335956/(3-4 sin² 1°)……(1)

于是就用迭代法反求sin1°。由于sin 1°是未知量，暂时取一个近似值，先代入sin² 1°，按(1)式计算，便得较好的sin 1°，再将此较好的sin 1°第二次代入sin² 1°，得第三次sin 1°，此后反复代入，直到

sin 1°收敛，不再变动，即为结果。计算表明：即使用误差极大的值，如

sin 1°= 0.1作第一次近似值，也只趋近六次。用较好的初始近似值时，迭代次数可以减少。所谓“较好的初始近似值”是用弧长来代替。弧长值为：α/180*3.1416=1/180*3.1416=0.01745333。在算例中第一次sin1°取0.01745，经过三次迭代收敛，最后sin 1°= 0.017452406。

现将迭代运算Sin 1°= 0.052335956/(3-4 sin² 1°) 结果列于下表，以见趋近情况。

 

算例	趋近	sin 1° 近似值	sin 1° 渐近结果
1	1	0.1000000000	0.0176810659
	2	0.0176810659	0.0174525930
	3	0.0174525930	0.0174524062
	4	0.0174524062	0.0174524060
	5	0.0174524060	0.0174524060
	6	0.0174524060	0.0174524064
	7	0.0174524064	0.0174524064
2	1	0.0174500000	0.0174524045
	2	0.0174524045	0.0174524064
	3	0.0174524064	0.0174524064
	4	0.0174524064	0.0174524064

 

    最后采用sin 1°= 0.017452400，有意在第九位上留有误差6的误差，它将会给运算结果带来某些积累误差。

 

 

五     1 度间隔的(0°—90°)正弦函数表的编算

 

    α由0开始，β=1°，(且始终为1)，按循环公式α_i+1=(α_i+1°)及和差公式sin(α+β)=sinαcosβ+cosαsinβ一步年步向后推算，便得sin1°、 sin2° 、sin3° 、sin4°……sin90°

    手算过程如下：

第一步  己知 α=0， sin 0°= 0，cos0°=1，β=1，

sin 1°= 0.017452400，cos 1°=0.999847695

sin 1°=sin(0+1)=sin0 cos1+cos0 sin1

=0.*0.999847695+1.* 0.017452400= 0.017452400、

cos1°=√（1- 0.017452400* 0.017452400）=0.999847695

 

第二步 己知 α=1， sin 1= 0.017452400，cos 1=0.999847695，β=1，

sin 1= 0.017452400，cos 1=0.999847695

sin2°=sin(1+1)=sin1cos1+cos1sin1

=0.017452400*0.999847695+0.999847695* 0.017452400= 0.034899484

cos2°=√（1- 0.034899484*0.034899484）=0.999390827

 

第三步 己知 α=2，sin 2= 0.034899484，cos 2=0.999390827，β=1，sin 1= 0.017452400，

cos 1=0.999847695

sin3°=sin(2+1)=sin2 cos1+cos2 sin1

= 0.034899484*0.999847695+0.999390827* 0.017452400=0.052335937

cos3°=√（1- 0.052335937*0.052335937）=0.998629536

 

第四步 己知 α=3，sin 3= 0.052335937，cos 3=0.998629536，β=1，sin 1= 0.017452400，

cos 1=0.999847695

sin4°=sin(3+1)=sin3 cos1+cos3 sin1

= 0.052335937*0.999847695+0.998629536*0.017452400=0.069756448

cos4°=√（1- 0.069756448*0.069756448）=0.997564052

 

第五步α=4，β=1，……循环下去。计算结果整理在表A的sin计算值中。

 

表A   1 度间隔的正弦函数  （ sin 1°= 0.017452400     用 9度控制改正）

α°	Sin 计算值	控制点值	改正值	改正后Sin α	检验 Sin α
0	0.000000000	0.	0	0.000000000	0.000000000
1	0.017452400		6	0.017452406	0.017452406
2	0.034899484		13	0.034899497	0.034899497
3	0.052335937		19	0.052335956	0.052335956
4	0.069756448		25	0.069756473	0.069756474
5	0.087155711		32	0.087155743	0.087155743
6	0.104528425		38	0.104528463	0.104528463
7	0.121869299		44	0.121869343	0.121869343
8	0.139173050		51	0.139173101	0.139173101
9	0.156434408	0.156434465	57	0.156434465	0.156434465
10	0.173648114		63	0.173648177	0.173648178
11	0.190808926		69	0.190808995	0.190808995
12	0.207911615		75	0.207911690	0.207911691
13	0.224950972		81	0.224951053	0.224951054
14	0.241921808		87	0.241921895	0.241921896
15	0.258818951		93	0.258819044	0.258819045
16	0.275637256		99	0.275637355	0.275637356
17	0.292371599		105	0.292371704	0.292371705
18	0.309016883	0.309016994	111	0.309016994	0.309016994
19	0.325568038		116	0.325568154	0.325568154
20	0.342020021		121	0.342020142	0.342020143
21	0.358367822		126	0.358367948	0.358367950
22	0.374606461		131	0.374606592	0.374606593
23	0.390730991		137	0.390731128	0.390731128
24	0.406736501		142	0.406736643	0.406736643
25	0.422618115		147	0.422618262	0.422618262
26	0.438370995		152	0.438371147	0.438371147
27	0.453990343	0.453990500	157	0.453990500	0.453990500
28	0.469471402		161	0.469471563	0.469471563
29	0.484809455		164	0.484809619	0.484809620
30	0.499999831		168	0.499999999	0.500000000
31	0.515037902		172	0.515038074	0.515038075
32	0.529919087		175	0.529919262	0.529919264
33	0.544638855		179	0.544639034	0.544639035
34	0.559192720		183	0.559192903	0.559192903
35	0.573576249		186	0.573576435	0.573576436
36	0.587785062	0.587785252	190	0.587785252	0.587785252
37	0.601814830		192	0.601815022	0.601815023
38	0.615661280		194	0.615661474	0.615661475
39	0.629320193		196	0.629320389	0.629320391
40	0.642787409		198	0.642787607	0.642787610
41	0.656058827		200	0.656059027	0.656059029
42	0.669130402		202	0.669130604	0.669130606
43	0.681998154		204	0.681998358	0.681998360
44	0.694658163		206	0.694658369	0.694658370
45	0.707106573	0.707106781	208	0.707106781	0.707106781
46	0.719339591		208	0.719339799	0.719339800
47	0.731353491		208	0.731353699	0.731353702
48	0.743144615		208	0.743144823	0.743144825
49	0.754709369		208	0.754709577	0.754709580
50	0.766044232		209	0.766044441	0.766044443
51	0.777145750		209	0.777145959	0.777145961
52	0.788010543		209	0.788010752	0.788010754
53	0.798635300		209	0.798635509	0.798635510
54	0.809016785	0.809016994	209	0.809016994	0.809016994
55	0.819151836		207	0.819152043	0.819152044
56	0.829037366		205	0.829037571	0.829037573
57	0.838670363		202	0.838670565	0.838670568
58	0.848047893		200	0.848048093	0.848048096
59	0.857167100		198	0.857167298	0.857167301
60	0.866025205		196	0.866025401	0.866025404
61	0.874619511		193	0.874619704	0.874619707
62	0.882947400		191	0.882947591	0.882947593
63	0.891006335	0.891006524	189	0.891006524	0.891006524
64	0.898793860		185	0.898794045	0.898794046
65	0.906307605		180	0.906307785	0.906307787
66	0.913545279		176	0.913545455	0.913545458
67	0.920504679		171	0.920504850	0.920504853
68	0.927183685		167	0.927183852	0.927183855
69	0.933580262		162	0.933580424	0.933580426
70	0.939692461		158	0.939692619	0.939692621
71	0.945518421		153	0.945518574	0.945518576
72	0.951056367	0.951056516	149	0.951056516	0.951056516
73	0.956304613		142	0.956304755	0.956304756
74	0.961261559		135	0.961261694	0.961261696
75	0.965925696		128	0.965925824	0.965925826
76	0.970295603		121	0.970295724	0.970295726
77	0.974369948		115	0.974370063	0.974370065
78	0.978147491		108	0.978147599	0.978147601
79	0.981627081		101	0.981627182	0.981627183
80	0.984807659		94	0.984807753	0.984807753
81	0.987688254	0.987688341	87	0.987688341	0.987688341
82	0.990267991		77	0.990268068	0.990268069
83	0.992546082		68	0.992546150	0.992546152
84	0.994521835		58	0.994521893	0.994521895
85	0.996194647		48	0.996194695	0.996194698
86	0.997564009		39	0.997564048	0.997564050
87	0.998629503		29	0.998629532	0.998629535
88	0.999390805		19	0.999390824	0.999390827
89	0.999847684		10	0.999847694	0.999847695
90	1.000000000	1	0	1.000000000	1.000000000

 

 

            六  1 度间隔的正弦函数表的误差分配及精度分析

 

    表A第一列为整度数。

    表A第二列，是sin计算值。本想这就是结果了，哪知与已知值一比，出乎意料，发现有差值，还真不小，最大209，即第七位上差2，这是为什么？细看表A就会发现，累次“滚雪球”加1°后，由于sin 1°的原始取位误差，引起了运算结果的误差。且每计算一个循环，都会把误差积累下来。在各个控制点0°、9°、18°…90°上，误差从0、57、111、…到209，再由209、189、…0。数值大小由sin 1°的误差引起。误差由小到大、由大到小的变化，则由计算公式特性引起。

    表A第三列，为sin已知值。这些已知值，不是靠电脑查出来的，而在《二编算的基础》中早已列出了的，在此是作为一个控制值用的，检查“滚雪球”后的计算值是否与已知值相同。结果不符值为：

sin 0  sin 9   sin 18   sin 27  sin 36  sin 45  sin 54   sin63  sin 72  sin 81  sin 90

  0     57      111      157     190     208     209     189     149     87      0

    57      54       46       33      18       1      -20-    40     -62     -87

       -3       -8       -13     -15      -17     -21      -20    -22-    -25

            5       -5        -2     - 2      -4       -1      -2      -3

    上面四行是误差的各次差分。

分析。这是一条开口向下的抛物误差曲线。共十个大区间，每区间的一次差分已成一条斜线，二次差分变成平缓的斜线，三次差分已近乎一条水平线。而每个大区间中，又有九个小区间。如果在每一个大区间中，按线性内插，则每一点上，每度，其平均内插误差在2以下，最大不会超过25 / 9 = 3的误差，即在54°到90°之间，可能会有3的误差。

于是将上述误差，在每个大区间内，分别作线性内插，把误差分配、分摊掉，算得第四列改正值。

    表A第四列，为改正值，只列第九位上的尾数。

    表A第五列，为改正后Sin α  =  Sin α计算值  + 改正值。

    表A第六列，列出正确的Sin α，这才由电脑提供。

改正后的Sin α，与正确的电脑Sin α相比，改正值残剩误差为3的，有14个，占16﹪，且分布在50°之后，正与上述分析相符。残剩误差为2的，有21个，占24﹪。残剩误差为0、1的，有55个，占60﹪。这样的结果，虽然还有点遗憾，但我只能满意了。如果再想提高其表面精度，只能把控制区间定为6°、甚至3°了。

    九位正弦函数表的精度，到底精确到什么程度？先举个例，例如sin 30°=0.5，而sin 30°0’0.001”=0.500000004，说明：九位上即使有误差4，对应的角度误差影响也仅为0.001秒。这样看来，八位正弦函数表精确到0.01”， 七位正弦函数表精确到0.1” 六位正弦函数表精确到1”。一般来说，使用六位正弦函数表就可以了。

现在，表A中每度的Sin，就可作为已知值，成为编算更小间隔(如1’、10“)的正弦表的控制值了。

 

    在整个运算作业中，我的思路慢慢清晰了。用“滚雪球”法编算正弦表，关键两条：

    一，要有一个较精确的起算值sin 1°。编1°的正弦表，取sin 1°= 0.017452400，已足够精确。但如果没有控制、没有检验，“雪球”就越滚越大，引起计算值不准。非但无法判断其精度，甚至算错了也不知道。所以：

    二，一定要有几个控制点、控制值。用它们检验计算的正确性，并根据不符值，对每一计算点作改正。有了控制，即使sin 1°有误差，也不可怕。没有控制，即使sin 1°再精确，也不放心。至于控制点要多少？控制区间有多宽，依我一个星期的来回反复计算，认为：如果控制区间为9度，则正弦值精度，能正确到第八位。如果控制区间为6度，将正确到第九位。如果不作控制，则只能正确到第五位。

 

             七   表A   1度间隔正弦函数的内插及内插精度

 

    有人要问：用这样的方法，所编算的正弦函数表，内插精度又如何？现分别用实例来说明表A的内插精度。

    例：Sin60°30’=？

查表A、得   Sin 60°=0.866025401    Sin61°=0.874619704  

    Sin60°30’=( 0.866025401+0.874619704)/2=0.870322552

再用0.870322552，按反三角正弦公式β= Arc Sinβ，反算得：

β= Arc (0.870322552)=60°29’46”，这就与出发的 60°30’有14”的误差，也就是表A的内插精度。其他算例如下：

 

度  分	Sin 内插值	反算 度 分 秒	误差(秒)
05  30	0.095842103	05  29  59	1
14  30	0.250370469	14  29  58	2
30  30	0.507519037	30  29  55	5
44 30	0.700882575	44  29  52	8
60  30	0.870322552	60  29  46	14
75 30	0.968110774	75  29  30	30
82  30	0.991407109	82  29  04	56
88 30	0.999619259	88  25  08	248

 

    可见1°度一载的正弦函数，自身的精度虽然很高，达0.001”，但内插精度只能达到秒级。小角的内插，误差在10”以下。角度越大，内插误差也越大。在90度附近，达到4分，根本不能内插。因此我认为，三角函数表最合理的安排是1’一载的六位三角函数表。

 

                        八   说   明

 

    1’间隔、10”间隔这两份表的编算方法，与1°间隔的编算方法完全相仿只是取sin 1’=0.000290887、Sin 10”=0.000048481368为起算。但因为表格太长，所以割爱不载了。

正弦函数表有了，其他余弦、正切、余切函数表也就可以造了，所以也不再编算了。

    噢，应该说一下计算过程中的体味。那就是编算过程中，光是加、减、乘、除、还可以承受，但在用“和差公式”时要算Cos=√（1—s i n 2），必须开平方，这工作量就大了。我知道手工方法开方，但我也不是堂·吉诃德。我是用计算器代劳的，造表时又用“电子表格”公式化成批地代算，辛苦不了什么，反而促进了我解题的兴趣与毅力。所以一个星期内就编出1度、1分、10秒间隔的部份正弦函数表。除了在检验结果时，我必须用计算器内的正弦值外，全部运算过程绝不用计算器、计算机内的现存正弦值代替，否则还算个屁。我要全过程，只用加、减、乘、除、开平方来计算sin，效仿古人、效仿古法，如此而已，岂有他哉。

 

                       九   结 束 语

 

    于是我返回现代。

我发现，浅绿色封面、硬精装的《六位三角函数表》已经绝迹，可能在我原先工作的单位“湖南省地质测绘队”还有，但早也不再使用，只是存放在后勤科的仓库里，并蒙着厚厚的灰尘了。

    我了了三个心愿：用古代人可能采用的只靠加、减、乘、除、开平方的方法，编算了常用对数表、计算了圆周率π、编算了三角正弦函数表。自我感觉良好。

    学习、求知是我的天性。阿门！

2013--08--  20～～28   在佛山 年76