しょき - プログラム倉庫

#include <stdio.h>
#include <math.h>
double mlx(double w1,double w2,double tr){
double l1,l2,mls;
l1=1-(w1+w2+tr)/(3*w1);
l2=1-(w1+w2+tr)/(3*w2);
if (l1<0)  l2=(w2-tr)/(2*w2);
if (l1<0)  l1=0;
if (l2<0)  l1=(w1-tr)/(2*w1);
if (l2<0)  l2=0;
mls=l1;
return mls;
}

double flx(double w1,double w2,double tr){
double l1,l2,fls;
l1=1-(w1+w2+tr)/(3*w1);
l2=1-(w1+w2+tr)/(3*w2);
if (l1<0)  l2=(w2-tr)/(2*w2);
if (l1<0)  l1=0;
if (l2<0)  l1=(w1-tr)/(2*w1);
if (l2<0)  l2=0;
fls=l2;
return fls;
}



int main(void){
double th1[11],th2[11];
double w1,w2,tk,tl,tr,c1,l1,l2;
double tr1,tr2,tr3,b1,b2,wel;
int s1,s2,t,n1,n2;
double maxtk,maxtr,maxwel,maxtl;
double y1[11][11],y2[11][11],c[11][11];
for (s1=1;s1<=10;s1++){
th1[s1]=0.2*s1;
}
for (s2=1;s2<=10;s2++){
th2[s2]=0.1*s2;
}
maxwel=-999;
for (n1=10;n1<=40;n1++){
for (n2=10;n2<=40;n2++){
tk=0.01*n1;
tl=0.01*n2;
tr1=0.1;
tr2=0.05;
b1=0;
for (s1=1;s1<=10;s1++){
for (s2=1;s2<=10;s2++){
w1=(1-tk)*th1[s1];
w2=(1-tl)*th2[s2];
l1=mlx(w1,w2,tr1);
l2=flx(w1,w2,tr1);
c1=w1*l1+w2*l2+tr1;
b1=b1+th1[s1]*l1+th2[s2]*l2-c1;
}
}
t=0;
while (t<100){
b2=0;
for (s1=1;s1<=10;s1++){
for (s2=1;s2<=10;s2++){
w1=(1-tk)*th1[s1];
w2=(1-tl)*th2[s2];
l1=mlx(w1,w2,tr2);
l2=flx(w1,w2,tr2);
c1=w1*l1+w2*l2+tr2;
b2=b2+th1[s1]*l1+th2[s2]*l2-c1;
}
}
tr3=tr2-b2*(tr2-tr1)/(b2-b1);
tr1=tr2;
tr2=tr3;
b1=b2;
if (b2*b2<0.001)t=1000;
t=t+1;
}
tr=tr2;
wel=0;
for (s1=1;s1<=10;s1++){
for (s2=1;s2<=10;s2++){
w1=(1-tk)*th1[s1];
w2=(1-tl)*th2[s2];
l1=mlx(w1,w2,tr);
l2=flx(w1,w2,tr);
c1=w1*l1+w2*l2+tr;
wel=wel+log(c1)+log(1-l1)+log(1-l2);
}
}
if (wel>maxwel)maxtk=tk;
if (wel>maxwel)maxtl=tl;
if (wel>maxwel)maxtr=tr;
if (wel>maxwel)maxwel=wel;
}
}
printf("%f",maxtk);
printf("\n");
printf("%f",maxtl);
printf("\n");
printf("%f",maxtr);
tk=maxtk;
tl=maxtl;
tr=maxtr;
for (s1=1;s1<=10;s1++){
for (s2=1;s2<=10;s2++){
w1=(1-tk)*th1[s1];
w2=(1-tl)*th2[s2];
y1[s1][s2]=th1[s1]*mlx(w1,w2,tr);
y2[s1][s2]=th2[s2]*flx(w1,w2,tr);
c[s1][s2]=(1-tk)*y1[s1][s2]+(1-tl)*y2[s1][s2]+tr;
}
}
return 0;
}