Программа на СИ

[F1mka]

Была написана программа, которая много чего ищет. По заданию, все функции нужно считать без библиотеки math.h Все функции должны считаться приближенными методами вычисления. И тут Загвостка, Натуральный логариф(9) нужно вычислить комбинированным методом хорд и секущих, а не так как я посчитал на калькуляторе. И ума никак не приложу, как сделать это.

#include <stdio.h>
double p=3.14159;
double exp(double x);
double cos(double x);
double fabs(double x);
double sqrt(double x);
int main()
{
double ln01=-2.30258;
double A0=0.75; // Объявление энергию сигала
double f0=10000;// Объявление частоты несущего колебания
double delf;// Ввод девиации частоты радиосигнала
double Tc=0.015;// Объявление общей длительности радиосигнала
double bigW;// Ввод частоты модулирующей функции
double t; // Ввод времени
double T=0.007; // Объявление общей длительности радиосигнала
double alf; // Ввод коэффициента длительности радиосигнала
double fd=80000; // Объявление частоты дискретизации
double delt; // Ввод периода дискретизации
double si;
double sj;
double sk;
double s1;
double s2;
double s3;
double x;
double y;
double E; // Ввод Энергии радиосигнала
double ts; // Ввод среднего времени радиосигнала
double Te; // Ввод эффективной длительность радиосигнала
double Ed; // Ввод Энергии производной радиосигнала
int i; // Ввод Счетчика итераций
int N=Tc*fd; // Объявление Числа отсчетов
delf=0.05*f0; // Присваивание значения девиации частоты
bigW=1.0/Tc; // Присваевание значения частоты модулирующей функции
alf=-4*ln01/(T*T); // Присваевание значения коэффициента длительности радиосигнала
delt=1/fd; // Присваевание значения периоду дискретизации
s1=0; // Объявление значения интеграла в заданный момент
s2=0; // Объявление значения интеграла в заданный момент
s3=0; // Объявление значения интеграла в заданный момент
E=0; // Объявление значения энергии радиосигнала
ts=0; // Объявление значения среднего времени радиосигнала
Te=0; // Объявление значения эффективной длительности радиосигнала
Ed=0; // Объявление значения энергии производной
// Открытие файла "out.txt" на запись
FILE *out_file;
out_file = fopen("out.txt", "w");
FILE *out1_file;
out1_file = fopen("out1.txt", "w");
for (i=0; i < N; i=i+1) // Цикл перебора всех i-х точек на оси времени
{
t=delt*i-Tc/2; // Нахождение t(i)
x=A0*exp(-alf*t*t); // Нахождение x(i)
s1=s1+2*delf*delt*p*cos(p/2-bigW*t); // Нахождение интегралф в период времени t(i)
y=cos(2*f0*p*t+s1); // Нахождение y(i)
si=x*y; // Вычисление цифрового отсчета сигнала в данный момент времени i
E=E+delt*si*si; // Вычисление энергии сигнала
ts=ts+t*delt*si*si; // Вычисление средней продолжительности радиосигнала
fprintf(out_file, "%4.8f\n", si); // Вывод i-ого отсчета в файл "out.txt"
fprintf(out1_file, "%4.8f\n", s1);
}
ts=ts/E; // Вычисление средней продолжительность радиосигнала
for (i=0; i < (N-1); i=i+1) // Цикл перебора всех i-х точек на оси времени
{
t=delt*i-Tc/2; // Нахождение t(i)
x=A0*exp(-alf*t*t); // Нахождение x(i)
s1=s1+2*delf*delt*p*cos(p/2-bigW*t); // Нахождение интеграла в период времени t(i)
y=cos(2*f0*p*t+s1); // Нахождение y(i)
si=x*y; // Вычисление цифрового отсчета сигнала в данный момент времени i
t=delt*(i+1)-Tc/2; // Нахождение t(i+1)
x=A0*exp(-alf*t*t); // Нахождение x(i+1)
s2=s2+2*delf*delt*p*cos(p/2-bigW*t); // Нахождение интеграла в период времени t(i+1)
y=cos(2*f0*p*t+s2); // Нахождение y(i+1)
sj=x*y; // Вычисление цифрового отсчета сигнала в данный момент времени k
t=delt*(i+2)-Tc/2; // Нахождение t(i+2)
x=A0*exp(-alf*t*t); // Нахождение x(i+2)
s3=s3+2*delf*delt*p*cos(p/2-bigW*t); // Нахождение интеграла в период времени t(i+2)
y=cos(2*f0*p*t+s3); // Нахождение y(i+2)
sk=x*y; // Вычисление цифрового отсчета сигнала в данный момент времени k
Te=Te+si*si*(i*delt-Tc/2-ts)*(i*delt-Tc/2-ts); // Нахождение эффективной длительносттирадиосигнала
Ed=Ed+(sk-sj)*(sk-sj)/delt; // Нахождение эффективной энергии производной
}
Te=sqrt(Te*delt/E); // Вычисление эффективной длительности радиосигнала
printf("E= %f\n", E); // Вывод на экран значения энергии сигнала
printf("ts= %f\n", ts); // Вывод на экран значения средней продолжительности радиосигнала
printf("Te= %f\n", Te); // Вывод на экран значения эффективной длительности радиосигнала
printf("Ed= %f\n", Ed); // Вывод на экран значения энергии производной
fclose(out_file); // Закрытие файла "out.txt"
fclose(out1_file);
printf("\nPress ENTER key to exit..."); // Вывод на экран строки запроса нажатия клавиши Enter
scanf("%*c"); // Закрытие файла при нажатии на клавишу Enter
return 0; // Возвращение значения функции
}

[F1mka]

double fabs (double x) // Функция возвращает модуль числа с плавующей точкой
{
if (x < 0) // Проверка Условия, является ли входной аргумент отрицательным числом
  {
  x = -x;
  }
return x; // Возвращение значения функции
}
double exp(double x)
{
  double yn, uk, uk_1;// Объявление переменных, куда записываются знчения последовательных приближений
  int k;// Объявление переменной, определяющей номер текущей итерации
  double e = 1.0e-8;// Предельная абсолютная погрешность вычислений
  k = 1;
  uk = 1.0;// Установка начального условия
  yn = uk;// значение искомой функции на текущей итерации
  while (fabs(uk) > e) // Цикл вычисления итераций
    {
    uk_1 = uk; // Переход к новой итерации
    uk = x/k*uk_1; // Вычисление n-й итерации
    k = k + 1; // увеличение номера текущей итерации
    yn = yn + uk; // значение искомой функции на текущей итерации
    }
  x=yn;
  return x; // Возвращение значения функции
}
double cos(double x)
{
    double yn, uk, uk_1,uk1; // Объявление переменных, куда записываются знчения последовательных приближений
    int k; // Объявление переменной, определяющей номер текущей итерации
    double e = 1.0e-8;// Предельная абсолютная погрешность вычислений
    double a=0;// Установка начального условия
    int b=0;// Установка начального условия
    k = 1; 
    uk = 1.0; // Установка начального условия
    yn = uk;// значение искомой функции на текущей итерации
    uk1=1.0;
    a=x/(p*2);
    b=a;
    x=(a-b)*2*p;
    while (uk1 > e) // Цикл вычисления итераций
      {
        uk_1 = uk; // Переход к новой итерации
        uk = -x*x/(2*k)/(2*k-1)*uk_1; // Вычисление n-й итерации
        k = k + 1; // увеличение номера текущей итерации
        yn = yn + uk; // значение искомой функции на текущей итерации
        if(uk<=0)  // Задание условий при uk меньше равным нулю
          {
            uk1=-uk;
          }
      }
    x=yn;
    return x; // Возвращение значения функции
}
double sqrt(double x)
{
  double yn, yn_1; // Объявление переменных, куда записываются знчения последовательных приближений
  double e = 1.0e-8;// Предельная абсолютная погрешность вычислений
  yn_1 = x; // Установка начального условия
  yn = 0.5*(yn_1 + x/yn_1); // Выполнение первой итерации
  while (fabs(yn-yn_1) > e) // Цикл вычисления итераций
  {
    yn_1 = yn; // Переход к новой итерации
    yn = 0.5*(yn_1 + x/yn_1); // Вычисление n-й итерации
  }
  return yn; // Возвращение значения функции
 }