void qromb(double a, double b, double& ss)
{
    double dss,eps = 0.000001;
    const int jmax = 20;
    int jmaxp = jmax + 1;
    int j,k = 5;
    int km = k - 1;
    double s[21], h[21];
    h[1] = 1.0;
    for (j = 1; j<=jmax; j++)
	{
        trapzd(a, b, s[j], j);
        if (j >= k)
		{
            polint(h, s, k, 0.0, ss, dss);
            if (fabs(dss) < eps * fabs(ss))
			{
				return;
			}
        }
        s[j + 1] = s[j];
        h[j + 1] = 0.25 * h[j];
    }
    cout<<"too many steps."<<endl;
}