Metoda trapezów

Witam. Mamy do wykonania następujące zadanie:

Funkcja
double trapez(double p, double q, double h)
{
Return h * (p + q)/2;
}
Oblicza pole trapezu o bokach p, q i wysokości h. Zastosować tą funkcję do obliczenia
pola powierzchni figury ograniczonej parabolą 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ść obliczeń
  Uwaga: Do obliczenia pola zastosować tzw. metodę trapezów (całkowania
  numerycznego). Polega ona na obliczaniu sum trapezów powstających w wyniku
  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ść:
  |Fn+1 – Fn|/|Fn| ≤ ε
  Za wynik przyjąć wartość Fn+1\

Czy mógłby ktoś zweryfikować moją implementację (na zadanych powyżej warunkach) tej metody? Będę wdzięczny za pomoc.

#include <iostream>
#include <cmath> // abs
using namespace std;
 
float a, b, c;

double trapez(double p, double q, double h)
{
    return h * (p + q)/2;
}

double func (float x)
{
	return a*x*x + b*x + c;
}

int main()
{
	int n = 2;
	float xp, xk, e, dx;
	double S = 0, S1 = 1;
	
	cout << "Wprowadz wspolczynnik a: "; cin >> a;
	cout << "Wprowadz wspolczynnik b: "; cin >> b;
	cout << "Wprowadz wspolczynnik c: "; cin >> c;
	
	cout << "Podaj xp: "; cin >> xp;
	cout << "Podaj xk: "; cin >> xk;
	
	cout << "Podaj dokladnosc e: "; cin >> e;
	

	S1 = trapez (func(xp), func(xk), (xk - xp));
	
	do
	{
		dx = (xk - xp)/n; // wysokość pojedynczego trapezu
		
		S = S1;
		S1 = 0;
		
		for (int i = 0; i < n; i++) S1 += trapez (func(xp + i*dx), func(xp + (i + 1)*dx), dx); // sumowanie n - trapezów
		
		n*=2;
		
	} while ((abs(S1) - abs(S))/abs(S) > e); // warunek "dokladności"
	
	cout << "\n\nPole S wynosi w przyblizeniu: " << abs(S1);
	cin.ignore();
	cin.get();
}