elliptic function, rf.

[transoid_z]

I just bought the book, and copied, yet, changed as following;

#ifndef ELLIPTCARL_H
#define ELLIPTCARL_H
#include <cmath>
#include <stdlib.h>
#include <iostream.h>
#include "MINE/utility.h"

double rf(const double x, const double y, const double z)
{
const double error = 0.0025, tiny = 1.5e-38, big = 3.0e37, third = 1.0/3.0;
const double C1 = 1.0/24.0, C2 = 0.1, C3 = 3.0/44.0, C4 = 1.0/14.0;
double alamb, ave, delx, dely, delz, e2, e3, sqrtx, sqrty, sqrtz, xt, yt, zt;

if (min(min(x,y),z) < 0.0 || min(min(x+y, x+z), y+z) < tiny || max(max(x,y),z)>big)
{
cout << "invalid arguments in rf" << endl;
exit(1);
}

xt=x; yt=y; zt=z;

do
{
sqrtx=sqrt(xt); sqrty=sqrt(yt); sqrtz=(zt);
alamb = sqrtx*(sqrty+sqrtz)+sqrty*sqrtz;
xt=0.25*(xt+alamb);
yt=0.25*(yt+alamb);
zt=0.25*(zt+alamb);
ave=third*(xt+yt+zt);
delx=(ave-xt)/ave;
dely=(ave-yt)/ave;
delz=(ave-zt)/ave;
}while(max(max(fabs(delx), fabs(dely)), fabs(delz))>error);

e2=delx*dely-delz*delz;
e3=delx*dely*delz;

return (1.0+(C1*e2-C2-C3*e3)*e2+C4*e3)/sqrt(ave);
};

I changed it such that I don't have to use "nr.h", and in "MINE/utility.h" has the functions for max, min. Now I also wrote the function for "complete elliptic integral of the first kind in Legendre form" which is

double comp_ellk(const double ak)
{
return rf(0.0, (1.0-sqr(ak)), 1.0);
};

It seems everything is in order, by the way, the function "sqr(double) is a function that is also defined in "MINE/utility.h", and it complies fine. The problem is that the it gives WRONG ANSWER!!! For example, comp_ellk(1/2) should gives the answer of 1.8541, but it gives 3.54227. I practically copied whole code from the book. The only thing different is that I didn't use "nr.h". Or do I have to? Thanks in advance.

p.s. I am using VC++ 6.0