/*
**************************************************************
 Вычисление интеграла различными методами при разной
 задаваемой в цикле относительной точности с автоматическим
 подбором числа делений отрезка интегрирования.
 Вариант для PASCAL INTEG01.PAS выдает АБСОЛЮТНО те же значения,
 но ГОРАЗДО ДОЛЬШЕ в PASCALABC. В PASCALABCNET скорость хорошая.
 Вариант для BASIC INTEG01.BAS дает АБСОЛЮТНО те же значения,
 но работает медленно (особенно при относительной точности -8).
 Реализованы следующие методы счета:
 - 1-левых прямоугольников
 - 2-правых прямоугольников
 Разница двух методов (левый и правый) состоит в том,
 что в первом случае суммируются значения функции на
 левых границах подоотрезков, а во-втором - на правых.
 Реально это означает, что в первом - сумма всех, кроме
 последнего, а во-втором - всех, кроме первого.
 - 3 -прямоугольников (или среднего). Это метод Гаусса при n=1.
 Третий метод - это сумма всех значений функции в
 серединах подотрезков.
 - 4-трапеций.
 Четвертый метод - трапеций - это сумма всех значений
 функции в точках со второй по предпоследнюю плюс
 полусумма первой и последней.
 - 5-Симпсона.
 Суммирование значений функции во всех нужных
 точках. Таковыми являются все границы подинтервалов
 (кроме начала первого и конца последнего) с коэффициентом
 2, середины всех интервалов с коэф.=4 и значения в точках,
 равных границе всего интервала с коэф.=1.
 В программе для упрощения все вдвое меньше и в конце
 умножается на 2. Сами коэф равны 2 или 1, а граничные на 12.
 В получаемой таблице столбцы имеют следующий смысл:
 ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН
-eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ
 - -eps-значение степени 10 для принятой в этой строке относительной точности
 - ЛЕВЫЙ..СИМПСОН-метод вычисления интеграла
 - ЛП КЛП..СИМ КСИМ-значения вычисленного интеграла и кол-во раз счета функции.
 Ниже приводится принт-скрин экрана.
*****************************************************************************
Введите значение левой и правой границы 0 1
Точное значение равно 0.785398
 ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН
-eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ


3 0.798441 189 0.772125 189 0.786231 14 0.784241 25 0.785398 12
4 0.789771 1652 0.781075 1710 0.785655 44 0.785054 75 0.785398 12
5 0.786794 16109 0.784008 16289 0.785462 170 0.785319 297 0.785398 12
6 0.785841 159894 0.784956 160460 0.785412 779 0.78538 1222 0.785398 21
7 0.785538 1594004 0.785258 1594004 0.785401 3402 0.785394 5457 0.785398 32
8 0.785442 15924545 0.785354 15924545 0.785399 15399 0.785397 24750 0.785398 45
******************************************************************************
 Точное значение интеграла = ARCTG(1)-ARCTG(0)=0.7853981634=PI/4.
*******************************************************************
*/
#include <conio.h>
#include <stdio.h>
#include <math.h>
long int kk;
// Подинтегральная функция.
double f(double q)
{
 kk++;
 return (1/(1+q*q));
}

// Первообразная функция (если она известна).
// Если не известна, то можно оставить как есть 
// и не обращать на это внимания.
double ff(double q)
{
 return (atan(q));
}

// Функция вычисления интеграла.
// Сюда передаю пределы интеграла, число отрезков деления
// всего интервала и тип вычисляемого интеграла:
double integral(double a, double b, long int k, int it)
{
 double x;
 double sum=0; // сюда буду накапливать интеграл
 long int i=0; // будет счетчиком кол-ва просуммированных подинтервалов
 int k1=1; // будет иметь значения 21
// kk=0;
 double y=(b-a)/k; // длина каждого из подотрезков отрезка
 if (it<=2) // Левый и Правый Прямоугольник
 for (x=a,sum=(it==1)?f(a):f(b); ++i<k;sum+=f(x+=y));
 if (it==3) // Прямоугольник
 for (x=a-y/2; ++i<=k; sum+=f(x+=y));
 if (it==4) // Трапеция
 for (x=a,sum=((f(a)+f(b))/2); ++i<k; sum+=f(x+=y));
 if (it==5) // Симпсон
 for (i=1,x=a,sum=(f(a)+f(b))/2,y/=2;i++<2*k; 
 sum+=((k1=3-k1)*f(x+=y)));
 // окончательный интеграл есть полученная сумма,
 // умноженная на длину каждого подотрезка.
 if (it==5) // Для Симпсона
 return 2*sum*y/3;
 else // Для всех других
 return sum*y;
}

// Главная программа
void main()
{
 const kn=2; // начальное кол-во подотрезков. При =1 возможно деление на 0.
 double ts=0,temp,a,b,d;
 long int k=kn;
 int l,i;
 clrscr();
 printf("Введите значение левой и правой границы ");
 scanf("%lf%lf",&a,&b);
 printf("Точное значение равно %lfn",ff(b)-ff(a));
 puts(" ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН");
 puts("-eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ");
 // Внешний цикл по относительной точности счета
 for (d=1e-3,l=3;d>1e-8;d/=10,l++)
 {
 // Цикл по возможным методам счета интеграла
 for (printf("n%1d",l),i=1; i<=5; ts=0, k=kn, i++)
 {
 kk=0; // счетчик обращений к вычислению подинтегральной функции для 
 // данного метода при указанной точности. 
 // вычисляю интеграл при К=1. Для начала это предыдущее.
 // Циклить до того момента, пока относительное отношение
 // интеграла при большем и меньшем числе отрезков
 // не достигнет заданной относительной точности.
 // Если еще не достиг, то для следующего раза
 // ныне новое будет старым.
 while (fabs(((temp=integral(a,b,k++,i))-ts)/(ts=temp))>d);
 // если еще не достиг, то для следующего раза
 // ныне новое будет старым.
 // Сюда по достижении точности.
 // Вывожу сами вычисленные интегралы и полученные числа делений отрезков.
 if (i<3)printf("%9.6lg%9ld", temp,kk);
 if (i==3)printf("%9.6lg%6ld", temp,kk);
 if (i==4)printf("%9.6lg%6ld", temp,kk);
 if (i==5)printf("%9.6lg%3ld", temp,kk);
 }
 }
}

//**************************************************************
// Вычисление интеграла различными методами при разной
// задаваемой в цикле относительной точности с автоматическим
// подбором числа делений отрезка интегрирования.
// Вариант для PASCAL INTEG01.PAS выдает АБСОЛЮТНО те же значения,
// но ГОРАЗДО ДОЛЬШЕ в PASCALABC. В PASCALABCNET скорость хорошая. 
// Вариант для BASIC INTEG01.BAS дает АБСОЛЮТНО те же значения,
// но работает медленно (особенно при относительной точности -8).
// Реализованы следующие методы счета:
// - 1-левых прямоугольников 
// - 2-правых прямоугольников
// Разница двух методов (левый и правый) состоит в том,
// что в первом случае суммируются значения функции на
// левых границах подоотрезков, а во-втором - на правых.
// Реально это означает, что в первом - сумма всех, кроме
// последнего, а во-втором - всех, кроме первого.
// - 3 -прямоугольников (или среднего). Это метод Гаусса при n=1.
// Третий метод - это сумма всех значений функции в 
// серединах подотрезков.
// - 4-трапеций.
// Четвертый метод - трапеций - это сумма всех значений
// функции в точках со второй по предпоследнюю плюс
// полусумма первой и последней.
// - 5-Симпсона.
// Суммирование значений функции во всех нужных
// точках. Таковыми являются все границы подинтервалов
// (кроме начала первого и конца последнего) с коэффициентом
// 2, середины всех интервалов с коэф.=4 и значения в точках,
// равных границе всего интервала с коэф.=1.
// В программе для упрощения все вдвое меньше и в конце
// умножается на 2. Сами коэф равны 2 или 1, а граничные на 12.
// В получаемой таблице столбцы имеют следующий смысл:
// ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН
//-eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ
// - -eps-значение степени 10 для принятой в этой строке относительной точности
// - ЛЕВЫЙ..СИМПСОН-метод вычисления интеграла
// - ЛП КЛП..СИМ КСИМ-значения вычисленного интеграла и кол-во раз счета функции.
// Ниже приводится принт-скрин экрана.
//*****************************************************************************
//Введите значение левой и правой границы 0 1
//Точное значение равно 0.785398
// ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН
//-eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ
//
//
//3 0.798441 189 0.772125 189 0.786231 14 0.784241 25 0.785398 12
//4 0.789771 1652 0.781075 1710 0.785655 44 0.785054 75 0.785398 12
//5 0.786794 16109 0.784008 16289 0.785462 170 0.785319 297 0.785398 12
//6 0.785841 159894 0.784956 160460 0.785412 779 0.78538 1222 0.785398 21
//7 0.785538 1594004 0.785258 1594004 0.785401 3402 0.785394 5457 0.785398 32
//8 0.785442 15924545 0.785354 15924545 0.785399 15399 0.785397 24750 0.785398 45 
//******************************************************************************
// Точное значение интеграла = ARCTG(1)-ARCTG(0)=0.7853981634=PI/4.
//*******************************************************************

uses crt;
var
kk:longint;
 kn:integer;
 ts,temp,a,b,d:real;
 k:longint;
 l,i:integer;
// Подинтегральная функция.
function f(q:real):real;
begin
 inc(kk);
 f:= (1/(1+q*q));
end;

// Первообразная функция (если она известа).
// Если не известна, то можно оставить как есть 
// и не обращать на это внимания.
function ff(q:real):real;
begin
 ff:= (arctan(q));
end;

// Функция вычисления интеграла.
// Сюда передаю пределы интеграла, число отрезков деления
// всего интервала и тип вычисляемого интеграла:
function integral(a,b:real;k, it:integer):real;
var
 x: real;
 sum:real ; // сюда буду накапливать интеграл
 i:longint; // будет счетчиком кол-ва просуммированных подинтервалов
 k1:integer ; // будет иметь значения 21
 y:real;
 begin
 k1:=1;
 y:=(b-a)/k; // длина каждого из подотрезков отрезка
 if (it<=2) then// Левый и Правый Прямоугольник
 begin
 if (it=1) then sum:=f(a)
 else sum:=f(b);
 x:=a+y; 
 for i:=1 to k-1 do
 begin
 sum:=sum+f(x);
 x:=x+y;
 end; 
 end; 
 if (it=3) then // Прямоугольник
 begin
 x:=a+y/2;
 sum:=0;
 for i:=1 to k do
 begin
 sum:=sum+f(x);
 x:=x+y;
 end; 
 end;
 if (it=4) then// Трапеция
 begin
 x:=a+y;
 sum:=((f(a)+f(b))/2);
 for i:=1 to k-1 do
 begin
 sum:=sum+f(x);
 x:=x+y;
 end; 
 end;
 if (it=5) then// Симпсон
 begin
 sum:=(f(a)+f(b))/2;
 y:=y/2;
 x:=a+y;
 for i:=1 to 2*k-1 do
 begin 
 k1:=3-k1;
 sum:=sum+f(x)*k1;
 x:=x+y;
 end; 
 end; 
 // окончательный интеграл есть полученная сумма,
 // умноженная на длину каждого подотрезка.
 if (it=5) then// Для Симпсона
 integral:= 2*sum*y/3
 else // Для всех других
 integral:= sum*y;
end ;

// Главная программа

 begin
 kn:=2; // начальное кол-во подотрезков. При =1 возможно деление на 0.
 k:=kn;
 write('Введите значение левой и правой границы ');
 readln(a,b);
 writeln('Точное значение равно ',ff(b)-ff(a));
 writeln(' ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН');
 writeln('-eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ');
 // Внешний цикл по относительной точности счета
 d:=1e-2;
 for l:=3 to 8 do
 begin
 d:=d/10;
 // Цикл по возможным методам счета интеграла
 write(l);
 for i:=1 to 5 do
 begin
 ts:=0;
 k:=kn;
 kk:=0; // счетчик обращений к вычислению подинтегральной функции для 
 // данного метода при указанной точности. 
 // вычисляю интеграл при К=1. Для начала это предыдущее.
 // Циклить до того момента, пока относительное отношение
 // интеграла при большем и меньшем числе отрезков
 // не достигнет заданной относительной точности.
 // Если еще не достиг, то для следующего раза
 // ныне новое будет старым.
 while(d>0) do
 begin
 temp:=integral(a,b,k,i);
 inc(k);
 if (abs((temp-ts)/temp)<=d) then break;
 ts:=temp;
 end;
 // если еще не достиг, то для следующего раза
 // ныне новое будет старым.
 // Сюда по достижении точности.
 // Вывожу сами вычисленные интегралы и полученные числа делений отрезков.
 if (i<3) then write(temp:9:6,kk:9);
 if (i=3) then write(temp:9:6,kk:6);
 if (i=4) then write(temp:9:6,kk:6);
 if (i=5) then write(temp:9:6,kk:4);
 end
 end
end.

'**************************************************************
' Вычисление интеграла различными методами при разной
' задаваемой в цикле относительной точности с автоматическим
' подбором числа делений отрезка интегрирования.
' Вариант для PASCAL INTEG01.PAS выдает АБСОЛЮТНО те же значения,
' но ГОРАЗДО ДОЛЬШЕ в PASCALABC. В PASCALABCNET скорость хорошая.
' Вариант для BASIC INTEG01.BAS дает АБСОЛЮТНО те же значения,
' но работает медленно (особенно при относительной точности -8).
' Реализованы следующие методы счета:
' - 1-левых прямоугольников
' - 2-правых прямоугольников
' Разница двух методов (левый и правый) состоит в том,
' что в первом случае суммируются значения функции на
' левых границах подоотрезков, а во-втором - на правых.
' Реально это означает, что в первом - сумма всех, кроме
' последнего, а во-втором - всех, кроме первого.
' - 3 -прямоугольников (или среднего). Это метод Гаусса при n=1.
' Третий метод - это сумма всех значений функции в
' серединах подотрезков.
' - 4-трапеций.
' Четвертый метод - трапеций - это сумма всех значений
' функции в точках со второй по предпоследнюю плюс
' полусумма первой и последней.
' - 5-Симпсона.
' Суммирование значений функции во всех нужных
' точках. Таковыми являются все границы подинтервалов
' (кроме начала первого и конца последнего) с коэффициентом
' 2, середины всех интервалов с коэф.=4 и значения в точках,
' равных границе всего интервала с коэф.=1.
' В программе для упрощения все вдвое меньше и в конце
' умножается на 2. Сами коэф равны 2 или 1, а граничные на 12.
' В получаемой таблице столбцы имеют следующий смысл AS
' ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН
'-eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ
' - -eps-значение степени 10 для принятой в этой строке относительной точности
' - ЛЕВЫЙ..СИМПСОН-метод вычисления интеграла
' - ЛП КЛП..СИМ КСИМ-значения вычисленного интеграла и кол-во раз счета функции.
' Ниже приводится принт-скрин экрана.
'*****************************************************************************
'Введите значение левой и правой границы 0 1
'Точное значение равно 0.785398
' ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН
'-eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ
'
'
'3 0.798441 189 0.772125 189 0.786231 14 0.784241 25 0.785398 12
'4 0.789771 1652 0.781075 1710 0.785655 44 0.785054 75 0.785398 12
'5 0.786794 16109 0.784008 16289 0.785462 170 0.785319 297 0.785398 12
'6 0.785841 159894 0.784956 160460 0.785412 779 0.78538 1222 0.785398 21
'7 0.785538 1594004 0.785258 1594004 0.785401 3402 0.785394 5457 0.785398 32
'8 0.785442 15924545 0.785354 15924545 0.785399 15399 0.785397 24750 0.785398 45
'******************************************************************************
' Точное значение интеграла = ARCTG(1)-ARCTG(0)=0.7853981634=PI/4.
'*******************************************************************
' Прототипы используемых функций
DECLARE FUNCTION f# (q AS DOUBLE)
DECLARE FUNCTION ff# (q AS DOUBLE)
DECLARE FUNCTION integral# (a AS DOUBLE, b AS DOUBLE, k AS LONG, it AS LONG)
DIM SHARED kk AS LONG
DIM kn AS LONG, k AS LONG, l AS LONG, i AS LONG
DIM ts AS DOUBLE, temp AS DOUBLE, a AS DOUBLE, b AS DOUBLE, d AS DOUBLE

' Главная программа
 kn = 2' начальное кол-во подотрезков. При =1 возможно деление на 0.
 k = kn
 INPUT "Введите значение левой и правой границы ", a, b
 PRINT "Точное значение равно "; ff(b) - ff(a)
 PRINT " ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН"
 PRINT "-eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ"
 ' Внешний цикл по относительной точности счета
 d = .01
 FOR l = 3 TO 8
 d = d / 10
 ' Цикл по возможным методам счета интеграла
 PRINT USING "##"; l;
 FOR i = 1 TO 5
 ts = 0
 k = kn
 kk = 0' счетчик обращений к вычислению подинтегральной функции для
 ' данного метода при указанной точности.
 ' вычисляю интеграл при К=1. Для начала это предыдущее.
 ' Циклить до того момента, пока относительное отношение
 ' интеграла при большем и меньшем числе отрезков
 ' не достигнет заданной относительной точности.
 ' Если еще не достиг, то для следующего раза
 ' ныне новое будет старым.
 DO WHILE (d > 0)
 temp = integral(a, b, k, i)
 k = k + 1
 IF (ABS((temp - ts) / temp) <= d) THEN EXIT DO
 ts = temp
 LOOP
 ' если еще не достиг, то для следующего раза
 ' ныне новое будет старым.
 ' Сюда по достижении точности.
 ' Вывожу сами вычисленные интегралы и полученные числа делений отрезков.
 IF (i < 3) THEN PRINT USING " #.###### ########"; temp; kk;
 IF (i = 3) THEN PRINT USING " #.###### #####"; temp; kk;
 IF (i = 4) THEN PRINT USING " #.###### #####"; temp; kk;
 IF (i = 5) THEN PRINT USING " #.###### ##"; temp; kk;
 NEXT i
 NEXT l
END

' Подинтегральная функция.
FUNCTION f# (q AS DOUBLE)
 kk = kk + 1
 f = 1 / (1 + q * q)
END FUNCTION

' Первообразная функция (если она известа).
' Если не известна, то можно оставить как есть
' и не обращать на это внимания.
FUNCTION ff# (q AS DOUBLE)
 ff = ATN(q)
END FUNCTION

FUNCTION integral# (a AS DOUBLE, b AS DOUBLE, k AS LONG, it AS LONG)
' Функция вычисления интеграла.
' Сюда передаю пределы интеграла, число отрезков деления
' всего интервала и тип вычисляемого интеграла AS
DIM x AS DOUBLE, sum AS DOUBLE ' сюда буду накапливать интеграл
DIM i AS LONG ' будет счетчиком кол-ва просуммированных подинтервалов
DIM k1 AS LONG ' будет иметь значения 21
DIM y AS DOUBLE
 k1 = 1
 y = (b - a) / k' длина каждого из подотрезков отрезка
 IF (it <= 2) THEN ' Левый и Правый Прямоугольник
 IF (it = 1) THEN sum = f(a) ELSE sum = f(b)
 x = a + y
 FOR i = 1 TO k - 1
 sum = sum + f(x)
 x = x + y
 NEXT i
 END IF
 IF (it = 3) THEN ' Прямоугольник
 x = a + y / 2
 sum = 0
 FOR i = 1 TO k
 sum = sum + f(x)
 x = x + y
 NEXT i
 END IF
 IF (it = 4) THEN ' Трапеция
 x = a + y
 sum = ((f(a) + f(b)) / 2)
 FOR i = 1 TO k - 1
 sum = sum + f(x)
 x = x + y
 NEXT i
 END IF
 IF (it = 5) THEN ' Симпсон
 sum = (f(a) + f(b)) / 2
 y = y / 2
 x = a + y
 FOR i = 1 TO 2 * k - 1
 k1 = 3 - k1
 sum = sum + f(x) * k1
 x = x + y
 NEXT i
 END IF
 ' окончательный интеграл есть полученная сумма,
 ' умноженная на длину каждого подотрезка.
 integral = sum * y
 ' Для Симпсона
 IF (it = 5) THEN integral = 2 * sum * y / 3
END FUNCTION

/***************************************************************
 Вычисление интеграла различными методами при разной
 задаваемой в цикле относительной точности с автоматическим
 подбором числа делений отрезка интегрирования.
 В получаемой таблице столбцы имеют следующий смысл:
 ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН
 - eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ
 - eps- степень 10 для принятой в этой строке отн. точности
 - ЛЕВЫЙ..СИМПСОН-метод вычисления интеграла
 - ЛП КЛП..СИМ КСИМ-значения интеграла и кол-во раз счета функ.
****************************************************************
Введите значение левой границы
0
Введите значение правой границы
1
Точное значение равно 0,78540

 ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН
-eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ


3 0,79844 00000189 0,77212 00000189 0,78623 00000014 0,78424 00025 0,78540 12
4 0,78977 00001652 0,78108 00001710 0,78566 00000044 0,78505 00075 0,78540 12
5 0,78679 00016109 0,78401 00016289 0,78546 00000170 0,78532 00297 0,78540 12
6 0,78584 00159894 0,78496 00160460 0,78541 00000779 0,78538 01222 0,78540 21
7 0,78554 01594004 0,78526 01594004 0,78540 00003402 0,78539 05457 0,78540 32
8 0,78544 15924545 0,78535 15924545 0,78540 00015399 0,78540 24750 0,78540 45
//******************************************************************************
// Точное значение интеграла = ARCTG(1)-ARCTG(0)=0.7853981634=PI/4.
//*******************************************************************
*/
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Integral1
{
 class Program
 {
 static int kk;
 static void Main()
 {
 const int kn=2; // начальное кол-во подотрезков. При =1 возможно деление на 0.
 double ts = 0, temp, a, b, d;
 int k = kn;
 int l, i;
 Console.WriteLine("Введите значение левой границы ");
 a = double.Parse(Console.ReadLine());
 Console.WriteLine("Введите значение правой границы ");
 b = double.Parse(Console.ReadLine());
 Console.WriteLine("Точное значение равно {0:f5}n", ff(b) - ff(a));
 Console.WriteLine(" ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН");
 Console.WriteLine("-eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ");
 // Внешний цикл по относительной точности счета
 for (d = 1e-3, l = 3; d > 1e-8; d /= 10, l++)
 {
 // Цикл по возможным методам счета интеграла
 Console.Write("n" + l);
 for ( i = 1; i <= 5; ts = 0, k = kn, i++)
 {
 kk=0; // счетчик обращений к вычислению подинтегральной функции для 
 // данного метода при указанной точности. 
 // вычисляю интеграл при К=1. Для начала это предыдущее.
 // Циклить до того момента, пока относительное отношение
 // интеграла при большем и меньшем числе отрезков
 // не достигнет заданной относительной точности.
 // Если еще не достиг, то для следующего раза
 // ныне новое будет старым.
 while (Math.Abs(((temp=integral(a,b,k++,i))-ts)/(ts=temp))>d);
 // если еще не достиг, то для следующего раза
 // ныне новое будет старым.
 // Сюда по достижении точности.
 // Вывожу сами вычисленные интегралы и полученные числа делений отрезков.
 if (i < 3) Console.Write(" {0:f5} {1:d8}", temp, kk);
 if (i == 3) Console.Write(" {0:f5} {1:d8}", temp, kk);
 if (i == 4) Console.Write(" {0:f5} {1:d5}", temp, kk);
 if (i == 5) Console.Write(" {0:f5} {1:d2}", temp, kk);
 }
 }
 Console.ReadLine();
 }

 // Подинтегральная функция.
 static double f(double q)
 {
 kk++;
 return (1/(1+q*q));
 }

 // Функция вычисления интеграла.
 // Сюда передаю пределы интеграла, число отрезков деления
 // всего интервала и тип вычисляемого интеграла:
 static double integral(double a, double b, int k, int it)
 {
 double x;
 double sum=0; // сюда буду накапливать интеграл
 int i=0; // будет счетчиком кол-ва просуммированных подинтервалов
 int k1=1; // будет иметь значения 21
 // kk=0;
 double y=(b-a)/k; // длина каждого из подотрезков отрезка
 if (it<=2) // Левый и Правый Прямоугольник
 for (x=a,sum=(it==1)?f(a):f(b); ++i<k;sum+=f(x+=y));
 if (it==3) // Прямоугольник
 for (x=a-y/2; ++i<=k; sum+=f(x+=y));
 if (it==4) // Трапеция
 for (x=a,sum=((f(a)+f(b))/2); ++i<k; sum+=f(x+=y));
 if (it==5) // Симпсон
 for (i=1,x=a,sum=(f(a)+f(b))/2,y/=2;i++<2*k; 
 sum+=((k1=3-k1)*f(x+=y)));
 // окончательный интеграл есть полученная сумма,
 // умноженная на длину каждого подотрезка.
 if (it==5) // Для Симпсона
 return 2*sum*y/3;
 else // Для всех других
 return sum*y;
 }

 // Первообразная функция (если она известа).
 // Если не известна, то можно оставить как есть 
 // и не обращать на это внимания.
 static double ff(double q)
 {
 return (Math.Atan(q));
 }
 }
}

#**************************************************************
# Вычисление интеграла различными методами при разной
# задаваемой в цикле относительной точности с автоматическим
# подбором числа делений отрезка интегрирования.
# В получаемой таблице столбцы имеют следующий смысл:
# ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН
# - eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ
# - eps- степень 10 для принятой в этой строке отн. точности
# - ЛЕВЫЙ..СИМПСОН-метод вычисления интеграла
# - ЛП КЛП..СИМ КСИМ-значения интеграла и кол-во раз счета функ.
#****************************************************************
"""
Введите границы зоны интегрирования a и b 0 1
Точное значение равно 0.78540

 ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН
-eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ

3 0.79844 189 0.77212 189 0.78623 14 0.78424 25 0.78540 12
4 0.78977 1652 0.78108 1710 0.78566 44 0.78505 75 0.78540 12
5 0.78679 16109 0.78401 16289 0.78546 170 0.78532 297 0.78540 12
6 0.78584 159894 0.78496 160460 0.78541 779 0.78538 1222 0.78540 21
7 0.78554 1594004 0.78526 1594004 0.78540 3402 0.78539 5457 0.78540 32
8 0.78544 15924545 0.78535 15924545 0.78540 15399 0.78540 24750 0.78540 45
"""
#******************************************************************************
# Точное значение интеграла = ARCTG(1)-ARCTG(0)=0.7853981634=PI/4.
#*******************************************************************

#import pdb; pdb.set_trace()
import math
# Подинтегральная функция.
def f(q):
 global kk
 kk+=1;
 return (1/(1+q*q));

# Первообразная функция (если она известа).
# Если не известна, то можно оставить как есть 
# и не обращать на это внимания.
def ff(q):
 return math.atan(q)
 
# Функция вычисления интеграла.
# Сюда передаю пределы интеграла, число отрезков деления
# всего интервала и тип вычисляемого интеграла:
def integral(a, b, k, it):
 sum=0 # сюда буду накапливать интеграл
 k1=1 # будет иметь значения 21
 y=(b-a)/k # длина каждого из подотрезков отрезка
 if it<=2: # Левый и Правый Прямоугольник
 sum=f(a) if it==1 else f(b)
 x=a
 for i in range(1,k): 
 x+=y
 sum+=f(x)
 if it==3: # Прямоугольник
 x=a-y/2
 for i in range(1,k+1):
 x+=y
 sum+=f(x)
 if it==4: # Трапеция
 sum=((f(a)+f(b))/2)
 x=a
 for i in range(1,k):
 x+=y
 sum+=f(x)
 if it==5: # Симпсон
 sum=(f(a)+f(b))/2
 x=a
 y=y/2
 for i in range(1,2*k):
 k1=3-k1
 x+=y
 sum+=(k1*f(x))
 # окончательный интеграл есть полученная сумма,
 # умноженная на длину каждого подотрезка.
 # Для Симпсона (it=5) и всех других
 if it==5:
 return 2*sum*y/3
 else:
 return sum*y
 

 
(a,b)= list(map(float, input('Введите границы зоны интегрирования a и b ').split()))
kn=2 # начальное кол-во подотрезков. При =1 возможно деление на 0.
ts = 0
k = kn;
print('Точное значение равно %8.5fn'%(ff(b) - ff(a)));
print(' ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН')
print('-eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ')
# Внешний цикл по относительной точности счета
d=1e-2
for l in range(3,9):
 d=d/10
 # Цикл по возможным методам счета интеграла
 print('%d'%(l),end='')
 for i in range(1,6):
 k=kn
 ts=kk=0 # счетчик обращений к вычислению подинтегральной функции для 
 # данного метода при указанной точности. 
 # вычисляю интеграл при К=1. Для начала это предыдущее.
 # Циклить до того момента, пока относительное отношение
 # интеграла при большем и меньшем числе отрезков
 # не достигнет заданной относительной точности.
 # Если еще не достиг, то для следующего раза
 # ныне новое будет старым.
 temp=integral(a,b,k,i)
 while (math.fabs((temp-ts)/temp)>d):
 # если еще не достиг, то для следующего раза ныне новое будет старым.
 ts=temp
 # увеличиваю кол-во отрезков для следующего счета
 k+=1
 temp=integral(a,b,k,i)
 if i<=3: 
 print('%8.5f %8d'%(temp,kk), end='')
 elif i==4:
 print('%8.5f %6d'%(temp,kk), end='')
 else:
 print('%8.5f %3d'%(temp,kk), end='n')
input() # команда ДЕРЖИТ экран

!**************************************************************
! Вычисление интеграла различными методами при разной
! задаваемой в цикле относительной точности с автоматическим
! подбором числа делений отрезка интегрирования.
! В получаемой таблице столбцы имеют следующий смысл:
! ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН
! - eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ
! - eps- степень 10 для принятой в этой строке отн. точности
! - ЛЕВЫЙ..СИМПСОН-метод вычисления интеграла
! - ЛП КЛП..СИМ КСИМ-значения интеграла и кол-во раз счета функ.
!****************************************************************
! Введите значение левой и правой границ
!0 1
! Точное значеие=
! 0.785398163397448
!
! ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН
! -eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ
! 3 0.79844 189 0.77212 189 0.78623 14 0.78424 25 0.78540 12
! 4 0.78977 1652 0.78108 1710 0.78566 44 0.78505 75 0.78540 12
! 5 0.78679 16109 0.78401 16289 0.78546 170 0.78532 297 0.78540 12
! 6 0.78584 159894 0.78496 160460 0.78541 779 0.78538 1222 0.78540 21
! 7 0.78554 1594004 0.78526 1594004 0.78540 3402 0.78539 5457 0.78540 32
! 8 0.78544 15924545 0.78535 15924545 0.78540 15399 0.78540 24750 0.78540 45
!Press any key to continue
!**********************************************************************
! Точное значение интеграла = ARCTG(1)-ARCTG(0)=0.7853981634=PI/4.
!*******************************************************************
!*******************************
! Основная программа
!*******************************
program INTEGRALL ! заголовок основной программы
 common kk ! переменная, общая для всех
 real*8 ts,temp,a,b,d,ffb
 integer k,kn,kk,l,i,it
 ! объявляю массивы с указанием диапазона индексов
 real*8,dimension(1:5)::tem
 integer*4,dimension(1:5)::ki
 ts=0 
 kn=2 ! начальное кол-во подотрезков. При =1 возможно деление на 0.
 k=kn 
 print *, trim(RuDosWin('Введите значение левой и правой границ', .false.))
 read *,a,b
 ffb= ff(b)-ff(a)
 print *, trim(RuDosWin('Точное значеие=', .false.))
 print *,ffb
 print *
 print *, trim(RuDosWin(' ЛЕВЫЙ ПРАВЫЙ ПРЯМОУГ ТРАПЕЦИЯ СИМПСОН', .false.))
 print *, trim(RuDosWin('-eps ЛП KЛП ПП КПП П КП ТР КТР СИМ КСИМ',.false.))
 ! Внешний цикл по относительной точности счета
 d=1e-3
 l=3
 do while (d>1e-8)
 i=1
 ! цикл по типам алгоритмов вычисления интегралов
 do while (i<=5)
 kk=0 ! счетчик обращений к вычислению подинтегральной функции для 
 ! данного метода при указанной точности. 
 ! вычисляю интеграл при К=1. Для начала это предыдущее.
 ! Циклить до того момента, пока относительное отношение
 ! интеграла при большем и меньшем числе отрезков
 ! не достигнет заданной относительной точности.
 ! Если еще не достиг, то для следующего раза
 ! ныне новое будет старым. 
 do while (1>0)
 ! вычисляю новое значение интеграла при заданных границах (a,b),
 !количестве делений - k и номере типа алгоритма счета интеграла - i 
 temp=integral1(a, b, k, i)
 k = k + 1 ! увеличиваю число делений отрезка интегрирования
 ! если достиг относительной точности вычисления интеграла, то выход
 if (abs((temp - ts) / temp) <= d) exit
 ! иначе новый для следующего будет старым
 ts = temp
 end do
 ! Сюда по достижении точности.
 ! Запоминаю в массивах сами вычисленные интегралы и полученные числа делений отрезков.
 tem(i)=temp
 ki(i)=kk
 ts=0
 k=kn
 ! перехожу к очередному алгоритму счета интеграла
 i=i+1
 end do
 ! вывод всего вычисленного для данной точности по всем типам алгоритмов
 1 format(I4,F9.5,I9,F8.5,I9,F8.5,I6,F8.5,I6,F8.5,I3)
 write(*,1)l,tem(1),ki(1),tem(2),ki(2),tem(3),ki(3),tem(4),ki(4),tem(5),ki(5)
 ! поднимаю точность 
 d=d/10
 l=l+1
 end do
contains ! далее идут описания внутренних функций

! Функция вычисления интеграла.
! Сюда передаю пределы интеграла, число отрезков деления
! всего интервала и тип вычисляемого интеграла:
 real*8 function integral1(a, b, k, it)
 ! тут определены типы передаваемых в функцию параметров
 real*8 a,b,x
 integer k,it
 real*8 sum ! сюда буду накапливать интеграл
 real*8 y ! длина каждого из подотрезков отрезка
 integer i ! будет счетчиком кол-ва просуммированных подинтервалов
 integer k1 ! будет иметь значения 21
 sum=0
 i=1
 k1=1
 y=(b-a)/k
 if (it.LE.2) then ! Левый и Правый Прямоугольник
 x=a 
 if (it.EQ.1) then
 sum=f(a) 
 else
 sum=f(b) 
 end if
 do while (i<k)
 x=x+y 
 sum=sum+f(x) 
 i=i+1 
 end do
 end if 
 if (it.EQ.3) then ! Прямоугольник
 x=a-y/2 
 do while (i<=k)
 x=x+y 
 sum=sum+f(x)
 i=i+1 
 end do
 end if 
 if (it.EQ.4) then ! Трапеция
 x=a
 sum=((f(a)+f(b))/2)
 do while (i<k)
 x=x+y 
 sum=sum+f(x)
 i=i+1 
 end do
 end if 
 if (it.EQ.5) then ! Симпсон
 sum=(f(a)+f(b))/2
 y=y/2
 x=a+y
 do while(i<2*k)
 k1=3-k1
 sum=sum+k1*f(x)
 x=x+y
 i=i+1
 end do
 end if 
! окончательный интеграл есть полученная сумма,
! умноженная на длину каждого подотрезка.
 if (it==5)then ! Для Симпсона
 integral1=2*sum*y/3
 return 
 else ! Для всех других
 integral1= sum*y
 return 
 end if
end function


! Подинтегральная функция.
real*8 function f(q)
! определяется тип передаваемого в функцию параметра
real*8 q
 common kk
 integer kk
 kk=kk+1
 ! прсваивается вычисленное значение имени функции
 f=(1/(1+q*q))
 ! возврат в точку вызова функции
 return 
end function
 
! Первообразная функция (если она известа).
! Если не известна, то можно оставить как есть 
! и не обращать на это внимания.
real*8 function ff(q)
! определяется тип передаваемого в функцию параметра
 real*8 q
 ff= atan(q)
 return 
end function ! конец очередной внутренней функции.
!*****************************************************************
! Модуль (функция), позволяющий вводить и выводить кириллицу в кодировке 866.
!*****************************************************************
!module TextTransfer
! Длина строк - результирующих переменных символьных функций
!integer(4), parameter :: ncresults = 250 
!contains
! Функция преобразовывает текст DOS в текст Windows, если dos_win = .TRUE., 
! и преобразовывает текст Windows в текст DOS, если dos_win = .FASLE.
function RuDosWin(string, dos_win) ! Длина строки string не должна превышать
 character(250) :: RuDosWin ! 250 символов
! Параметры string и dos_win имеют вид связи IN
! и, следовательно, не должны изменяться в RuDosWin
 character(*), intent(in) :: string
 logical(4), intent(in) :: dos_win
! dif - величина, на которую при преобразовании изменяется код буквы
 integer(2) :: i, dos_win_code, dif
 RuDosWin = string
 do i = 1, len_trim(RuDosWin)
! dos_win_code – DOS- или Windows-код символа
 dos_win_code = iachar(RuDosWin(i:i))
 dif = 0 ! dif больше нуля, если символ - русская буква
 if(dos_win) then ! Если преобразование DOS - Windows
 select case(dos_win_code) ! Найдем величину dif
 case(128 : 175) ! DOS-русские буквы от А до Я и от а до п
 dif = 64
 case(224 : 239) ! DOS-русские буквы от р до я
 dif = 16
 end select
 else ! Преобразование Windows - DOS
 select case(dos_win_code)
 case(192 : 239) ! Windows-русские буквы от А до Я и от а до п
 dif = -64
 case(240 : 255) ! Windows-русские буквы от р до я
 dif = -16
 end select
 end if
 ! Выполняем преобразование символа, если он является буквой русского алфавита 
 if(dif /= 0) RuDosWin(i:i) = char(dos_win_code + dif)
 end do
 end function RuDosWin
! end module TextTransfer 
end program INTEGRALL ! конец основной программы .