たｃ - プログラム倉庫

#include <stdio.h>
#include <math.h>
double f(double x,double mu,double sig){
double pi,x1,x2,x3,fx;
pi = 3.1415;
x1=-pow(x - mu,2) / (2*pow(sig,2));
x2=exp(x1);
x3=sig*pow(2*pi,0.5);
fx=x2/x3;
return fx;
}
double g(double y,double mu, double sig){
double gx,h,x;
int n,t;  
gx=0;
h=0.001;
t=floor(y/h);
for (n=-2000;n<t;n++){
x=n*h;
gx=gx+f(x,mu,sig)*h;
}
return gx;
}
double seeky(double p,double mu,double sig){
double g1,g2,y1,y2,y3;
int t;
y1=0.4;
y2=-0.2;
g1=g(y1,mu,sig);
t=0;
while(t<100){
g2=g(y2,mu,sig);
y3=y2+(p-g2)*(y2-y1)/(g2-g1);
y1=y2;
y2=y3;
g1=g2;
if (pow(p-g2,2)<0.0001)t=1000;
t=t+1;
}
return y2;
}
double bud(double tl,double tr,double th[101]){
int s;
double bx,w1,l1,y1,c1;
bx=0;
for (s=1;s<=100;s++){
w1=(1-tl)*th[s];
l1=(w1-tr)/(2*w1);
if (l1<0) l1=0;
y1=th[s]*l1;
c1=w1*l1+tr;
bx=bx+y1-c1;
}
return bx;
}
double trs(double tl,double th[101]){
double tr1,tr2,tr3,b1,b2;
int t;
tr1=0.05;
tr2=0.03;
b1=bud(tl,tr1,th);
t=0;
while(t<100){
b2=bud(tl,tr2,th);
tr3=tr2-b2*(tr2-tr1)/(b2-b1);
tr1=tr2;
tr2=tr3;
b1=b2;
if (b2*b2<0.01) t=1000;
t=t+1;
}
return tr2;
}
double wel(double tl,double tr,double th[101]){
double we,w1,l1,c1;
int s;
we=0;
for (s=1;s<=100;s++){
w1=(1-tl)*th[s];
l1=(w1-tr)/(2*w1);
if (l1<0) l1=0;
c1=w1*l1+tr;
we=we+log(c1)+log(1-l1);
}
return we;
}


int main(void)
{
double p,mu,sig,yy;
double th[101];
double maxwel,wep,maxtl,maxtr;
int n,s;
double tl,tr;
mu=0;
sig=0.39;
for (s=1;s<101;s++){
p=0.01*s-0.005;
yy=seeky(p,mu,sig);
th[s]=exp(yy);
}
maxwel=-999;
for (n=10;n<=40;n++){
tl=0.01*n;
tr=trs(tl,th);
wep=wel(tl,tr,th);
if (wep>maxwel) maxtr=tr;
if (wep>maxwel) maxtl=tl;
if (wep>maxwel) maxwel=wep;
}
printf("%f",maxtl);
return 0;
}