pole pod wykresem funkcji kwadratowej

Mam problem z tym:

Funkcja
double trapez(double p, double q, double h)
{
return h*(p + q)/2;
}
oblicza pole powierzchni trapezu o bokach a, b i wysokości h. Zastosować tą funkcję do
obliczenia pola powierzchni figury ograniczonej krzywą
y = ax2 + bx + c oraz prostymi y = 0, x = xp , x = xk z zadaną dokładnością ε.
Dane: a, b, c – współczynniki paraboli,
xp, xk - ograniczenia przedziału wartości x,
ε - dokładność obliczenia
Uwaga: Zastosować metodę polegającą na obliczaniu sumy pól trapezów stanowiących
wynik kolejnych podziałów przedziału <xp, xk> na 2, 4, 8, 16, ... części. Procedurę kolejnych
podziałów zakończyć jeżeli sumy pól trapezów częściowych po n-tym kroku podziału Fn i po
kolejnym n+1 kroku Fn+1 spełniają zależność:

[tex]\frac{|F_{n+1} - F_{n}|}{|F_{n}|} <=\epsilon[/tex]

Za wynik przyjąć wartość Fn+1

Bardzo prosze o szybka pomoc

sum = 0;
steps = 1;
do {
lastsum = sum;
sum = 0;
step = (xk - xp) / steps;
for (i = 0; i < steps; i++)
  sum += trapez(f(i * step + xp), f((i + 1) * step + xp), step);
step <<= 1;
} while (abs(sum - lastsum) / lastsum > epsilon)
cout << sum;

gdzie f(x) zwraca współrzędną y punktu na paraboli dla podanego x.