public class tax13{
public static void main(String[] args) {
double th1[]=new double[11];
double th2[]=new double[11];
double y1[][]=new double[11][11];
double y2[][]=new double[11][11];
double c[][]=new double[11][11];
double up[][]=new double[11][11];
double u[][][][]=new double[101][4][4][4];
double w[][][][]=new double[101][4][4][4];
double v[][][][][]= new double[100][3][3][3][21];
int goto1[][][][][]= new int[100][3][3][3][21];
int goto2[][][][][]= new int[100][3][3][3][21];
int goto3[][][][][]= new int[100][3][3][3][21];
int gotoq[][][][][]= new int[100][3][3][3][21];
int end1[][][]=new int[3][3][3];
int end2[][][]=new int[3][3][3];
int end3[][][]=new int[3][3][3];
int endq[][][]=new int[3][3][3];
int op1[]= new int[101];
int op2[]= new int[101];
int op3[]= new int[101];
int opq[]= new int[101];
double endv[][][]=new double[3][3][3];
int mrank[]=new int[101];
int frank[]=new int[101];
double tk,tl,tr;
double c1,w1,w2,l1,l2;
double we,minu,h;
double maxwel,maxtk,maxtl,maxtr;
double u1,v1,vs;
int s,n,n1,n2,n3,s1,s2,t1,t2,q,qx,pp;
int ns1,ns2,ns3,nx1,nx2,nx3,j,e,t;
for (s1=1;s1<11;s1++){
th1[s1]=0.2*s1;
}
for (s2=1;s2<11;s2++){
th2[s2]=0.1*s2;
}
maxwel=-999;
maxtk=0;
maxtl=0;
maxtr=0;
for (n1=5;n1<40;n1++){
for (n2=5;n2<40;n2++){
tk=0.01*n1;
tl=0.01*n2;
tr=trs(tk,tl,th1,th2);
we=wel(tk,tl,tr,th1,th2);;
if (we>maxwel)maxtk=tk;
if (we>maxwel)maxtl=tl;
if (we>maxwel)maxtr=tr;
if (we>maxwel)maxwel=we;
}
}
tk=maxtk;
tl=maxtl;
tr=maxtr;
for (s1=1;s1<11;s1++){
for (s2=1;s2<11;s2++){
l1=mlx(s1,s2,tk,tl,tr,th1,th2);
l2=flx(s1,s2,tk,tl,tr,th1,th2);
w1=(1-tk)*th1[s1];
w2=(1-tl)*th2[s2];
y1[s1][s2]=th1[s1]*l1;
y2[s1][s2]=th2[s2]*l2;
c[s1][s2]=w1*l1+w2*l2+tr;
up[s1][s2]=ux(c[s1][s2],l1,l2);
}
}
for (n=1;n<101;n++){
minu=999;
mrank[n]=1;
frank[n]=1;
for (s1=1;s1<11;s1++){
for (s2=1;s2<11;s2++){
if (up[s1][s2]<minu)mrank[n]=s1;
if (up[s1][s2]<minu)frank[n]=s2;
if (up[s1][s2]<minu)minu=up[s1][s2];
}
}
s1=mrank[n];
s2=frank[n];
up[s1][s2]=9999;
}
t=0;
h=0.01;
while(t<1000){
for (s=1;s<101;s++){
for (n1=-1;n1<2;n1++){
for (n2=-1;n2<2;n2++){
for (n3=-1;n3<2;n3++){
s1=mrank[s];
s2=frank[s];
c1=c[s1][s2]+n3*h;
l1=(y1[s1][s2]+n1*h)/th1[s1];
l2=(y2[s1][s2]+n2*h)/th2[s2];
u[s][n1+1][n2+1][n3+1]=ux(c1,l1,l2);
}
}
}
}
for (s=1;s<100;s++){
for (n1=-1;n1<2;n1++){
for (n2=-1;n2<2;n2++){
for (n3=-1;n3<2;n3++){
s1=mrank[s];
s2=frank[s];
t1=mrank[s+1];
t2=frank[s+1];
c1=c[s1][s2]+n3*h;
l1=(y1[s1][s2]+n1*h)/th1[t1];
l2=(y2[s1][s2]+n2*h)/th2[t2];
w[s][n1+1][n2+1][n3+1]=ux(c1,l1,l2);
}
}
}
}
for (n1=-1;n1<2;n1++){
for (n2=-1;n2<2;n2++){
for (n3=-1;n3<2;n3++){
for (q=-10;q<11;q++){
v[1][n1+1][n2+1][n3+1][q+10]=-999;
}
}
}
}
for (n1=-1;n1<2;n1++){
for (n2=-1;n2<2;n2++){
for (n3=-1;n3<2;n3++){
q=n2+n1-n3;;
v[1][n1+1][n2+1][n3+1][q+10]=u[1][n1+1][n2+1][n3+1];
}
}
}
for (s=2;s<100;s++){
for (n1=-1;n1<2;n1++){
for (n2=-1;n2<2;n2++){
for (n3=-1;n3<2;n3++){
for (q=-10;q<11;q++){
u1=u[s][n1+1][n2+1][n3+1];
qx=q+n1+n2-n3;
pp=0;
if (qx>10)pp=100;
if (qx<-10)pp=100;
if (pp>50)qx=0;
vs=-999;
ns1=0;
ns2=0;
ns3=0;
for (nx1=-1;nx1<2;nx1++){
for (nx2=-1;nx2<2;nx2++){
for (nx3=-1;nx3<2;nx3++){
v1=u1+v[s-1][nx1+1][nx2+1][nx3+1][qx+10];
if (w[s-1][nx1+1][nx2+1][nx3+1]>u1)v1=-999;
if (v1>vs)ns1=nx1;
if (v1>vs)ns2=nx2;
if (v1>vs)ns3=nx3;
if (v1>vs)vs=v1;
}
}
}
if (pp>50)vs=-999;
goto1[s][n1+1][n2+1][n3+1][q+10]=ns1;
goto2[s][n1+1][n2+1][n3+1][q+10]=ns2;
goto3[s][n1+1][n2+1][n3+1][q+10]=ns3;
gotoq[s][n1+1][n2+1][n3+1][q+10]=qx;
v[s][n1+1][n2+1][n3+1][q+10]=vs;
}
}
}
}
}
for (n1=-1;n1<2;n1++){
for (n2=-1;n2<2;n2++){
for (n3=-1;n3<2;n3++){
u1=u[100][n1+1][n2+1][n3+1];
qx=n1+n2-n3;
vs=-999;
ns1=0;
ns2=0;
ns3=0;
for (nx1=-1;nx1<2;nx1++){
for (nx2=-1;nx2<2;nx2++){
for (nx3=-1;nx3<2;nx3++){
v1=u1+v[99][nx1+1][nx2+1][nx3+1][qx+10];
if (w[99][nx1+1][nx2+1][nx3+1]>u1)v1=-999;
if (v1>vs)ns1=nx1;
if (v1>vs)ns2=nx2;
if (v1>vs)ns3=nx3;
if (v1>vs)vs=v1;
}
}
}
end1[n1+1][n2+1][n3+1]=ns1;
end2[n1+1][n2+1][n3+1]=ns2;
end3[n1+1][n2+1][n3+1]=ns3;
endq[n1+1][n2+1][n3+1]=qx;
endv[n1+1][n2+1][n3+1]=vs;
}
}
}
vs=-999;
ns1=0;
ns2=0;
ns3=0;
for (n1=-1;n1<2;n1++){
for (n2=-1;n2<2;n2++){
for (n3=-1;n3<2;n3++){
if (endv[n1+1][n2+1][n3+1]>vs)ns1=n1;
if (endv[n1+1][n2+1][n3+1]>vs)ns2=n2;
if (endv[n1+1][n2+1][n3+1]>vs)ns3=n3;
if (endv[n1+1][n2+1][n3+1]>vs)vs=endv[n1+1][n2+1][n3+1];
}
}
}
op1[100]=ns1;
op2[100]=ns2;
op3[100]=ns3;
op1[99]=end1[op1[100]+1][op2[100]+1][op1[100]+1];
op2[99]=end2[op1[100]+1][op2[100]+1][op2[100]+1];
op3[99]=end3[op1[100]+1][op2[100]+1][op3[100]+1];
opq[99]=endq[op1[100]+1][op2[100]+1][op3[100]+1];
for (j=1;j<99;j++){
s=99-j;
op1[s]=goto1[s+1][op1[s+1]+1][op2[s+1]+1][op3[s+1]+1][opq[s+1]+10];
op2[s]=goto2[s+1][op1[s+1]+1][op2[s+1]+1][op3[s+1]+1][opq[s+1]+10];
op3[s]=goto3[s+1][op1[s+1]+1][op2[s+1]+1][op3[s+1]+1][opq[s+1]+10];
opq[s]=gotoq[s+1][op1[s+1]+1][op2[s+1]+1][op3[s+1]+1][opq[s+1]+10];
}
e=0;
for (s=1;s<101;s++){
e=e+op1[s]*op1[s]+op3[s]*op3[s]+op2[s]*op2[s];
}
for (s=1;s<101;s++){
s1=mrank[s];
s2=frank[s];
c[s1][s2]=c[s1][s2]+op3[s]*h;
y1[s1][s2]=y1[s1][s2]+op1[s]*h;
y2[s1][s2]=y2[s1][s2]+op2[s]*h;
}
System.out.println(e);
if (e<2)h=h/2;
if (h<0.001)t=10000;
t=t+1;
}
}
static double bud(double tk,double tl,double tr,double th1[],double th2[]){
double b1,w1,w2,l1,l2,c1;
int s1,s2,i;
b1=0;
for (s1=1;s1<11;s1++){
for (s2=1;s2<11;s2++){
w1=(1-tk)*th1[s1];
w2=(1-tl)*th2[s2];
l1=1-(w1+w2+tr)/(3*w1);
l2=1-(w1+w2+tr)/(3*w2);
i=0;
if (l1<0)i=i+1;
if (l2>0)i=i+1;
if (i==2)l1=0;
if (i==2)l2=(w2-tr)/(2*w2);
i=0;
if (l1>0)i=i+1;
if (l2<0)i=i+1;
if (i==2)l1=(w1-tr)/(2*w1);
if (i==2)l2=0;
c1=w1*l1+w2*l2+tr;
b1=b1+th1[s1]*l1+th2[s2]*l2-c1;
}
}
return b1;
}
static double wel(double tk,double tl,double tr,double th1[],double th2[]){
double we,w1,w2,l1,l2,c1;
int s1,s2,i;
we=0;
for (s1=1;s1<11;s1++){
for (s2=1;s2<11;s2++){
w1=(1-tk)*th1[s1];
w2=(1-tl)*th2[s2];
l1=1-(w1+w2+tr)/(3*w1);
l2=1-(w1+w2+tr)/(3*w2);
i=0;
if (l1<0)i=i+1;
if (l2>0)i=i+1;
if (i==2)l1=0;
if (i==2)l2=(w2-tr)/(2*w2);
i=0;
if (l1>0)i=i+1;
if (l2<0)i=i+1;
if (i==2)l1=(w1-tr)/(2*w1);
if (i==2)l2=0;
c1=w1*l1+w2*l2+tr;
we=we+Math.log(c1)+Math.log(1-l1)+Math.log(1-l2);
}
}
return we;
}
static double trs(double tk,double tl,double[] th1,double[] th2){
double tr1,tr2,tr3,b1,b2;
int t1;
tr1=0.1;
tr2=0.05;
b1=bud(tk,tl,tr1,th1,th2);
t1=0;
while(t1<100){
b2=bud(tk,tl,tr2,th1,th2);
tr3=tr2-b2*(tr2-tr1)/(b2-b1);
b1=b2;
tr1=tr2;
tr2=tr3;
if (b2*b2<0.001)t1=1000;
t1=t1+1;
}
return tr2;
}
static double ux(double c1,double l1,double l2){
int pp;
double u1,cp1,lp1,lp2;
pp=0;
cp1=c1;
lp1=l1;
lp2=l2;
if (lp1<0)pp=100;
if (lp1>0.99)pp=100;
if (lp2<0)pp=100;
if (lp2>0.99)pp=100;
if (pp>50)lp1=0.5;
if (pp>50)lp2=0.5;
u1=Math.log(cp1)+Math.log(1-lp1)+Math.log(1-lp2);
if (pp>50)u1=-999;
return u1;
}
static double mlx(int s1,int s2,double tk,double tl,double tr,double[] th1,double th2[]){
int i;
double w1,w2,l1,l2;
w1=(1-tk)*th1[s1];
w2=(1-tl)*th2[s2];
l1=1-(w1+w2+tr)/(3*w1);
l2=1-(w1+w2+tr)/(3*w2);
i=0;
if (l1<0)i=i+1;
if (l2>0)i=i+1;
if (i==2)l1=0;
if (i==2)l2=(w2-tr)/(2*w2);
i=0;
if (l1>0)i=i+1;
if (l2<0)i=i+1;
if (i==2)l1=(w1-tr)/(2*w1);
if (i==2)l2=0;
return l1;
}
static double flx(int s1,int s2,double tk,double tl,double tr,double[] th1,double th2[]){
int i;
double w1,w2,l1,l2;
w1=(1-tk)*th1[s1];
w2=(1-tl)*th2[s2];
l1=1-(w1+w2+tr)/(3*w1);
l2=1-(w1+w2+tr)/(3*w2);
i=0;
if (l1<0)i=i+1;
if (l2>0)i=i+1;
if (i==2)l1=0;
if (i==2)l2=(w2-tr)/(2*w2);
i=0;
if (l1>0)i=i+1;
if (l2<0)i=i+1;
if (i==2)l1=(w1-tr)/(2*w1);
if (i==2)l2=0;
return l2;
}
}
最終更新：2010年05月10日 08:52
