t,i;
double o,w,h,x,y,k,a,b,c;
double g(N,S)double N,S[][2];{
/* Initially, let k hold the geometric mean of given y-values */
for(t=0;t<N;t++)
k+=S[t][1];
k/=N;
/* We approximate 9 times to ensure accuracy */
for(i=0;i<9;i++){
o=w=h=0;
for(t=0;t<N;t++)
/* Here, we are making running totals of partial derivatives */
/* o is the first, w the second, and h the third*/
x=S[t][0],
y=S[t][1],
a=y-k,
c=k*k-2*k*y+x*x+y*y,
o+=-a/sqrt(x*x+a*a),
w+=x*x/pow(c,1.5),
h+=3*x*x*a/pow(c,2.5);
/* We now use these derivatives to find a (hopefully) closer k */
a=h/2;
b=w-h*k;
c=o-w*k+a*k*k;
k=(-b+sqrt(b*b-4*a*c))/h;
}
return k;
}
/* Our testing code */
int main(int argc, char** argv) {
double test[2][2] = {
{5.7, 3.2},
{8.9, 8.1}
};
printf("%.20lf\n", g(2, test));
return 0;
}
But the first partial derivative of each Di is pretty bad...:
