- 注册时间
- 2007-12-27
- 最后登录
- 1970-1-1
- 威望
- 星
- 金币
- 枚
- 贡献
- 分
- 经验
- 点
- 鲜花
- 朵
- 魅力
- 点
- 上传
- 次
- 下载
- 次
- 积分
- 40157
- 在线时间
- 小时
|
发表于 2014-1-10 17:24:31
|
显示全部楼层
今天试着写了个Pari/gp的代码
- min_val()={1e-5;}
- polremove_zero(p)={
- local(q,d);
- q=0;d=poldegree(p);
- q=polcoeff(p,0);
- for(u=1,d,
- if(abs(polcoeff(p,u))>min_val(),
- q=q+polcoeff(p,u)*r^u)
- );
- q
- }
- find_best_3(a,b)={
- local(c,best);
- local(ca,sa,cb,sb,cc,sc);
- local(w1,w2,v,w3);
- local(u1,u2,z,u3);
- local(a1,a2,a3);
- local(m1,m2,m3,pr,rr,ts,tt,tr,smu,tmv);
- c=180-a-b;best=0.0;
- print("A="a";B="b";C="c);
- a=a*Pi/180;b=b*Pi/180;c=c*Pi/180;
- ca=cos(a);sa=sin(a);cb=cos(b);sb=sin(b);cc=cos(c);sc=sin(c);
- w1=2*sb*sa*(sc*r-sb*ca);w2=(sc*r-sb*ca)^2-(sb*sa)^2;
- v=(sc*r-sb*ca)^2+(sb*sa)^2;
- w3=w1*cb-w2*sb;
- u1=2*sc*sa*r*(sb-sc*ca*r);u2=(sb-sc*ca*r)^2-(sc*sa*r)^2;
- z=(sc*sa*r)^2+(sb-sc*ca*r)^2;
- u3=u1*cc-u2*sc;
- a1=((-v*sb-w3)*z*sc-v*sb*u3)*r;
- a2=((v*sb+w3)*z*sc*r+(v*sb*z*sc+v*sb*u3));
- a3=-v*sb*z*sc;
- m1=(4*a1*a3-a2^2);m2=(4*a1*w3*u3);
- m3=deriv(m1,r)*m2-m1*deriv(m2,r);
- m3=polremove_zero(m3);
- if(poldegree(m3)>0,
- pr=polroots(m3);
- for(u=1,length(pr),
- if(abs(imag(pr[u]))<min_val(),
- rr=real(pr[u]);
- if(subst(a1,r,rr)!=0,
- tt=subst(a2/(-2*a1),r,rr);
- ts=tt*rr;
- tr=subst(m1/m2,r,rr);
- smu=subst(sb*v/w3*(1-ts),r,rr);
- tmv=subst(sc*z/u3*(1-tt),r,rr);
- if(0<=tmv&&tmv<=tt&&tt<=1&&0<=smu&&smu<=ts&&ts<=1&&sb*sa/(sc*rr-sb*ca)>tan(b/2)&&sc*sa*rr/(sb-sc*ca*rr)>tan(c/2),print("t="tt";s="ts";u="ts-smu";v="tt-tmv";r="tr);if(tr>best,best=tr))
- )
- )
- )
- );best
- }
- find_best_bisector(A,B)=
- {
- local(C,tr,best);
- C=180-A-B;
- A=A*Pi/180;B=B*Pi/180;C=C*Pi/180;
- if(A<B,tr=A;A=B;B=tr);
- if(A<C,tr=A;A=C;C=tr);
- if(B<C,tr=B;B=C;C=tr);
- tr=sin(B)/(sin(A)+sin(B));best=sin(C)/(sin(B)+sin(C));
- if(tr>best,best=tr);
- best
- }
- find_best_perp(A,B)=
- {
- local(C,tr,best);
- C=180-A-B;
- A=A*Pi/180;B=B*Pi/180;C=C*Pi/180;
- if(A<B,tr=A;A=B;B=tr);
- if(A<C,tr=A;A=C;C=tr);
- if(B<C,tr=B;B=C;C=tr);
- tr=(tan(A)+tan(B))/(3*tan(A)+tan(B));
- best=(tan(B)+tan(C))/(3*tan(B)+tan(C));
- if(tr>best,best=tr);
- best
- }
- find_best_outside_triangle(A,B)=
- {
- local(C,best);
- C=180-A-B;best=0;
- A=A*Pi/180;B=B*Pi/180;C=C*Pi/180;
- if(A<B,tr=A;A=B;B=tr);
- if(A<C,tr=A;A=C;C=tr);
- if(B<C,tr=B;B=C;C=tr);
- if(C<Pi/4&&A<Pi/2,
- best = 1/(1+1/tan(A)+1/tan(C))
- );
- best
- }
- find_best_by_angle(A,B)=
- {
- local(C,tr,best);
- C=180-A-B;
- best=find_best_bisector(A,B);
- print("best bisector result "best);
- tr=find_best_perp(A,B);
- print("best perpandicular result "tr);
- if(tr>best,best=tr);
- tr=find_best_outside_triangle(A,B);
- if(tr>0,print("best outside triangle result "tr));
- if(tr>best,best=tr);
- tr=find_best_3(A,B);
- if(tr>best,best=tr);
- tr=find_best_3(B,C);
- if(tr>best,best=tr);
- tr=find_best_3(C,A);
- if(tr>best,best=tr);
- best
- }
- find_best_by_edge(a,b,c)=
- {
- local(A,B,C,tr,best);
- A=acos((b^2+c^2-a^2)/(2*b*c));
- B=acos((c^2+a^2-b^2)/(2*c*a));
- find_best_by_angle(A,B)
- }
复制代码 |
|