如何消除这段程序的计算误差（200分不犹豫！）？(200分)

下面是我的一段计算程序。
计算的正确结果应该是h=58.4505,sitac=160,t=55.89997。
可是这段程序的计算结果却是：h=53.5254,sitac=158.1714,t=74.4416。
如果把h=58.4505,sitac=160代入公式中，可以得出t=55.89997，可见公式没有错误，那么造成这种结果的原因只能是h和sitac的计算误差过大。
这该如何处理呢？请各位高手朋友赐教。
多谢！

longdo
uble A1,B1,C1,D1,r1,h1,s1,r2,h2,s2,r3,h3,s3,r4,h4,s4,r5,h5,s5,t,t1,t2,t3;
longdo
uble pi=3.1415926/180.0;
longdo
uble A2,B2,C2,D2,A3,B3,C3,D3;
longdo
uble X1,X2,X3,X4,Y1,Y2,Y3,Y4,E,F,G,h,ss,sitac;

r1=60.82608;h1=55.4160;s1=158.8167*pi;
r2=52.3740;h2=56.5904;s2=159.3864*pi;
r3=57.6858;h3=61.0784;s3=160.8904*pi;
r4=52.9642;h4=61.0030;s4=161.1023*pi;
r5=56.3365;h5=51.0672;s5=157.2643*pi;

A1=r1*r2*h3*sin(s2-s1)+r2*r3*h1*sin(s3-s2)-r1*r3*h2*sin(s3-s1);
B1=r1*r2*sin(s2-s1)+r2*r3*sin(s3-s2)-r1*r3*sin(s3-s1);
C1=r1*(h3-h2)*cos(s1)-r2*(h3-h1)*cos(s2)+r3*(h2-h1)*cos(s3);
D1=r1*(h3-h2)*sin(s1)-r2*(h3-h1)*sin(s2)+r3*(h2-h1)*sin(s3);

A2=r1*r2*h4*sin(s2-s1)+r2*r4*h1*sin(s4-s2)-r1*r4*h2*sin(s4-s1);
B2=r1*r2*sin(s2-s1)+r2*r4*sin(s4-s2)-r1*r4*sin(s4-s1);
C2=r1*(h4-h2)*cos(s1)-r2*(h4-h1)*cos(s2)+r4*(h2-h1)*cos(s4);
D2=r1*(h4-h2)*sin(s1)-r2*(h4-h1)*sin(s2)+r4*(h2-h1)*sin(s4);

A3=r1*r2*h5*sin(s2-s1)+r2*r5*h1*sin(s5-s2)-r1*r5*h2*sin(s5-s1);
B3=r1*r2*sin(s2-s1)+r2*r5*sin(s5-s2)-r1*r5*sin(s5-s1);
C3=r1*(h5-h2)*cos(s1)-r2*(h5-h1)*cos(s2)+r5*(h2-h1)*cos(s5);
D3=r1*(h5-h2)*sin(s1)-r2*(h5-h1)*sin(s2)+r5*(h2-h1)*sin(s5);

X1=A1*C2-A2*C1;
X2=A1*D2-A2*D1;
X3=B1*C2-B2*C1;
X4=B1*D2-B2*D1;
Y1=A1*C3-A3*C1;
Y2=A1*D3-A3*D1;
Y3=B1*C3-B3*C1;
Y4=B1*D3-B3*D1;
E=X1*Y3-X3*Y1;
F=X4*Y1+X3*Y2-X1*Y4-X2*Y3;
G=X2*Y4-X4*Y2;
ss=-F+sqrt(F*F-4*E*G);
sitac=atan(ss/(2*E))+3.1415926;
h=(Y1*sin(sitac)-Y2*cos(sitac))/(Y3*sin(sitac)-Y4*cos(sitac));

t1=A2-B2*h;
t2=C2*sin(sitac);
t3=D2*cos(sitac);
t=t1/(t2-t3);
sitac=sitac/pi;

1, 中间最好计算少点, 简化公式,
2, 看看

如此复杂的计算，编程是次要的，主要是数学。应该简化公式。将三角函数用
级数展开计算，只用加减乘除。这样，可以自己写高精度的计算公式。我曾看过
一个计算Pi的程序，就是这样。
要不然，用构件，好象有个150位的计算构件，忘记在哪了。要收费的。

呵呵， 我只做过整数的，
这样找数学系的人帮忙吧。

wo tou yun

将pi取到1000为试一试?
3.14159265358979323846264338327950288419716939937510......
呵呵,我只能背到50位

用一行写完你的算式,减少截断误差?

zensst：你说的那个计算Pi的程序的程序能不能给我发一份？（xianyu@cmmail.com).
另外，我发现double类型的数据只能取到小数点后面14或15位，如下面一段程序(我用的是C++builder,估计delphi也一样）：
double pi;
pi=3.1415926535897932384626433832795；
最终的输出结果却是pi=3.141592653589793，后面的数位全部丢失。即使换成longdo
uble也没用。这该如何是好？

是的，它有限制。看来用数组自己写吧。

咸鱼：你要计算的那个式子是什么？写出它的表达式出来，看我能不能帮你算算。
源程序太不好看了：） 你不是学软件的吧？
精度当然可以达到。。。。。位。据说这个mathematics可以完成数学专业硕士研究生的所有计算. 可惜我不是：）也许有了它我是半个。
beta:
Pi的1000位值。用mathmetics 4计算出来的。
3.1415926535897932384626433832795028841971693993751058209749445923078164062862/
089986280348253421170679821480865132823066470938446095505822317253594081284811/
174502841027019385211055596446229489549303819644288109756659334461284756482337/
867831652712019091456485669234603486104543266482133936072602491412737245870066/
063155881748815209209628292540917153643678925903600113305305488204665213841469/
519415116094330572703657595919530921861173819326117931051185480744623799627495/
673518857527248912279381830119491298336733624406566430860213949463952247371907/
021798609437027705392171762931767523846748184676694051320005681271452635608277/
857713427577896091736371787214684409012249534301465495853710507922796892589235/
420199561121290219608640344181598136297747713099605187072113499999983729780499/
510597317328160963185950244594553469083026425223082533446850352619311881710100/
031378387528865875332083814206171776691473035982534904287554687311595628638823/
53787593751957781857780532171226806613001927876611195909216420199

已经EMail了，希望有用。回信请写 zensst@yeah.net 。

我要计算的公式如下：
r1、r2....r5、sita1、sita2...sita5、h1、h2....h5都是已知的。
现在要计算rc、sitac、hc。
x1=r1*cos(sita1-sita_c)-rc x2=r2*cos(sita2-sita_c)-rc
y1=r1*sin(sinta1-sita_c) y2=r2*sin(sinta2-sita_c)
z1=h1-hc z2=h2-hc
x3=r3*cos(sita3-sita_c)-rc x4=r4*cos(sita4-sita_c)-rc
y3=r3*sin(sinta3-sita_c) y4=r4*sin(sinta4-sita_c)
z3=h3-hc z4=h4-hc
x5=r5*cos(sita5-sita_c)-rc xc=0
y5=r5*sin(sinta5-sita_c) yc=0
z5=h5-hc zc=0

| x1 y1 z1 1 |
| x2 y2 z2 1 |
| x3 y3 z3 1 | =0 ....................(1)
| xc yc zc 1 |
| x1 y1 z1 1 |
| x2 y2 z2 1 |
| x4 y4 z4 1 | =0 ....................(2)
| xc yc zc 1 |
| x1 y1 z1 1 |
| x2 y2 z2 1 |
| x5 y5 z5 1 | =0 ....................(3)
| xc yc zc 1 |

你这个利用单纯型法求解柱面坐标下的最优解吧？
既然你说公式没有错，那么引起误差的可能有几种：
1。pi的取值。这个好办，取到足够的精度即可。double型的可以取到15位，一定够了。
2。我感觉这个是很关键的：sin，cos的误差！
当a，b很接近时，a-b->0;sin->0
例如：
sin(0.01) = 9.9998333341666646825424382690997e-3
sin(0.001) = 9.9999983333334166666646825397101e-4
sin(0.0001)= 9.9999999833333333416666666646825e-5
当a很小时：sin(ka)=ksin(a);此时的关系近似为线形。
为了区分，此时必须取小数点后13-15位。若此时与一较大的数相加减，如100，double型的数位不够，引起截断误差。
解决办法：（个人意见）
设置一double型数组，记录所有的sin，cos值，计算过程中截取前N位代入计算，
并（截取后面的所有尾数）*10^N;

N的取值，按题意，取8？（试一试，应该好搞定的）
最后添加一减少误差的过程，用修改了的double型数据。用此可能会减少误差，达到所求的精度。
没试过，可能不正确-P，哎，分好难赚啊。。。。

我的数学不是很好。我始终觉得应该不用Sin Cos，最好把它展开。给你的程序
中有算Tan的，也就可以算Sin Cos了。

多人接受答案了。

zenset: 2者其实是一样的。