Может ли кто-нибудь дать мне отзыв о моем математическом движке C ++?

Я создавал математический движок для C ++ и решил, что сейчас хорошее время, чтобы попросить мнение сообщества. Сейчас основные «хорошие» вещи, которые у него есть, — это комплексные числа, исчисление и класс дробей. Я планирую добавить теорию множеств когда-нибудь в будущем, хотя не знаю, насколько это будет полезно для программистов. Вот мой код:

#include <iostream>
using namespace std;
namespace math{
const double mini = double (numeric_limits<double>::min());
const double e = 2.718281828459045;
const double pi = 3.141592653589793238462643383;
int def_prec = 16;
}

template <class t = double>
t exp(t base, int exponent){
if(exponent<0){
return 1/exp(base,-exponent);
} else if(exponent>0){
double output = 1;
for(int i=1;i<=exponent;i++){output *= base;}
return output;
} else {return 1;}
    
}

double nsqrt(double x, int precision = math::def_prec){
    double sum = 1;
    for (int i=1;i<=precision;i++){
        sum -= ((sum*sum)-x)/(2*sum);
    }
    return sum;
}

double esqrt(double num, double precision = math::def_prec){
    double h = num/precision;
    double y = 1;
    for (double x=1;x!=num;x+=h){
        y += h*(1/(2*y));
    }
    return y;
}

double sqrt(double num, int precision = math::def_prec){
    return nsqrt(num,precision);
}

double root(double num, int root, int precision = math::def_prec){
 double sum = 1;
  for (int i=1;i<=precision;i++){
      sum -= ((exp(sum,root))-num)/(root*exp(sum,root-1));
  }
  return sum;
}


class complex{
private:
    double r,i;
public:
    complex(double r_in,double i_in): r(r_in), i(i_in) {};
    double get_re() {return r;}
    double get_im(){return i;}
    double magnitude(int precision = math::def_prec){ return sqrt((r*r)+(i*i),precision);}
    complex operator + (complex a){
        r += a.get_re();
        i += a.get_im();
        return (*this);
    }
    complex operator + (double a){r += a;return (*this);}
    complex operator - (complex a){
        r -= a.get_re();
        i -= a.get_im();
        return (*this);
}
    complex operator - (double a){r -= a;return (*this);}
    complex operator * (complex a){
        r = (r*a.get_re())-(i*a.get_im());
        i = (r*a.get_im())+(i*a.get_re());
        return (*this);
    }
    complex operator * (double a){r *=a;i *= a;return (*this);}
    complex operator / (complex a){
           r = (r/a.get_re())-(i/a.get_im());
           i = (r/a.get_im())+(i/a.get_re());
           return (*this);
       }
    complex operator / (double a){r /=a;i /= a;return (*this);}
    void operator = (double a[2]){r=a[0];i=a[1];}
    void operator = (complex a){r=a.get_re();i=a.get_im();}
    bool operator == (complex a){return (r==a.get_re() && i==a.get_im());}
    bool operator == (double a[2]){return (r==a[0] && i==a[1]);}
    bool operator != (complex a){return (r!=a.get_re() || i!=a.get_im());}
    bool operator != (double a[2]){return (r!=a[0] || i!=a[1]);}
    void operator += (complex a){
        r += a.get_re();
        i += a.get_im();
    }
    bool operator > (complex a){return (r==a.get_re())? (i>a.get_im()) : (r>a.get_re());}
    bool operator < (complex a){return (r==a.get_re())? (i<a.get_im()) : (r<a.get_re());}
    bool operator <= (complex a){return (r==a.get_re())? (i<=a.get_im()) : (r<=a.get_re());}
    bool operator >= (complex a){return (r==a.get_re())? (i>=a.get_im()) : (r>=a.get_re());}
    void operator += (double a){r += a;}
    void operator -= (complex a){
        r -= a.get_re();
        i -= a.get_im();
    }
    void operator -= (double a){r -= a;}
    void operator *= (complex a){
        r = (r*a.get_re())-(i*a.get_im());
        i = (r*a.get_im())+(i*a.get_re());
    }
    void operator *= (double a){r *= a;i *= a;}
    void operator /= (complex a){
        r = (r/a.get_re())-(i/a.get_im());
        i = (r/a.get_im())+(i/a.get_re());
    }
    void operator /= (double a){r /= a;i /= a;}
    void operator ++ (){r++;}
    void operator -- (){r--;}
    
    ostream& operator<<(ostream& os)
    {
        os << r << "+" << i << "i";
        return os;
    }
    
    void operator = (int a){r=a;}
};

complex exp(complex base, int exponent){
if(exponent<0){
complex one(1,0);
return one/exp(base,-exponent);
} else if(exponent>0){
complex output(1,0);
for(int i=1;i<=exponent;i++){output *= base;}
return output;
} else {return base/base;}
    
}


template <class t = int>
t gcd(t X, t Y) {
    t pre = 0;
    if (X%Y == 0) {return (X>Y)? Y : X;}
    pre = X%Y;
    X = Y;
    Y = pre;
    while (X%Y != 0) {
        pre = X%Y;
        X = Y;
        Y = pre;
    }
    return pre;
}

int round(double x){return (x-int(x)>0.5)? int(x)+1 : int(x);}

class fraction{
private:
    int n,d;
public:
    fraction(int n_in,int d_in): n(n_in),d(d_in){};
    void simplify(){n/=gcd(n,d);d/=gcd(n,d);}
    double get_double(){return n/d;}
    int get_num(){return n;}
    int get_den(){return d;}
    fraction operator + (fraction a){
        n = (n*a.get_den())+(a.get_num()*d);
        d *= a.get_den();
        simplify();
        return (*this);
    }
    fraction operator + (int a){n += a*d;return (*this);}
    fraction operator - (fraction a){
        n = (n*a.get_den())-(a.get_num()*d);
        d *= a.get_den();
        simplify();
        return (*this);
    }
    fraction operator - (int a){n -= a*d;return (*this);}
    fraction operator * (fraction a){
        n *= a.get_num();
        d *= a.get_den();
        simplify();
        return (*this);
    }
    fraction operator * (int a){n *= a;return (*this);}
    fraction operator / (fraction a){
        n /= a.get_num();
        d /= a.get_den();
        simplify();
        return (*this);
    }
    fraction operator / (int a){n /= a;return (*this);}
    void operator += (fraction a){
        n = (n*a.get_den())+(a.get_num()*d);
        d *= a.get_den();
        simplify();
    }
    void operator += (int a){n += a*d;}
    void operator -= (fraction a){
        n = (n*a.get_den())-(a.get_num()*d);
        d *= a.get_den();
        simplify();
    }
    void operator -= (int a){n -= a*d;}
    void operator *= (fraction a){
        n *= a.get_num();
        d *= a.get_den();
        simplify();
    }
    void operator *= (int a){n*=a;}
    void operator /= (fraction a){
        n /= a.get_num();
        d /= a.get_den();
        simplify();
    }
    void operator /= (int a){n /= a;}
    void operator = (fraction a){n = a.get_num();d=a.get_den();}
    void operator = (int a){n=a;d=1;}
    
    bool operator < (fraction a){return (n*a.get_den()<a.get_num()*d);}
    bool operator > (fraction a){return (n*a.get_den()>a.get_num()*d);}
    bool operator >= (fraction a){return (n*a.get_den()>=a.get_num()*d);}
    bool operator <= (fraction a){return (n*a.get_den()<=a.get_num()*d);}
    bool operator == (fraction a){return (n*a.get_den()==a.get_num()*d);}
    bool operator != (fraction a){return (n*a.get_den()!=a.get_num()*d);}
    
    bool operator < (int a){return (n<a*d);}
    bool operator > (int a){return (n>a*d);}
    bool operator >= (int a){return (n>=a*d);}
    bool operator <= (int a){return (n<=a*d);}
    bool operator == (int a){return (n==a*d);}
    bool operator != (int a){return (n!=a*d);}
    operator double () {return n/d;}
};

fraction tofrac(double x, int precision = 1000000000){
    long gcd_ = gcd(round(x-int(x) * precision), precision);
    fraction output(round(x-int(x)*precision)/gcd_,precision/gcd_);
    return output;
}

template <class t = double>
t limit(double (*f)
    return (f(to+math::mini)+f(to-math::mini))/2;
}

template <class t = double>
t derivative(double (*f)
    return ((f(x+math::mini)-f(x))/math::mini+(f(x-math::mini)-f(x))/(-math::mini))/2;
}

template <class t = double>
t sigma(double (*f)
    t out = 0;
    for (int n=from;n<=to;n++){out += f(n);}
    return out;
}

template <class t = double>
t pi(double (*f)
    t out = 1;
    for (int n=from;n<=to;n++){out *= f(n);}
    return out;
}

double factorial1(int x) {
    double out = 1;
    for(int n=x;n>0;n--){out*=n;}
    return out;
}

template <class t = double>
t sin(t x, int precision = math::def_prec){
  t output = 0;
  for(int n=0;n<=precision;n++){
    output += exp(-1,n)*(exp(x,(2*n)+1)/factorial1((2*n)+1));
  }
    return output;
}

template <class t = double>
t cos(t x, int precision = math::def_prec){
    t output = 0;
    for(int n=0;n<=precision;n++){
      output += exp(-1,n)*(exp(x,2*n)/factorial1(2*n));
    }
      return output;
}

template <class t = double>
t tan(t x, int precision = math::def_prec){return sin(x,precision)/cos(x,precision);}

template <class t = double>
t ln(t x, int precision = math::def_prec){
    t output = 0;
    for(int n=1;n<=precision;n++){
        output += (exp(x-2,n)*exp(-1,n-1))/n;
    }
      return output;
}

complex ln(complex x, int precision = math::def_prec){
    complex output(0,0);
    for(int n=1;n<=precision;n++){
        output += (exp(x-1,n)*exp(-1,n-1))/n;
    }
      return output;
}

template <class t = double>
t log(t base, t x, int precision = math::def_prec){return ln(x,precision)/ln(base,precision);}

template <class t = double>
t derivative_fast(t (*f)

double pi(int precision = math::def_prec){
    double output = 0;
    for(int i=0;i<=precision+(precision%2);i++){
        output += exp(-1.0,i)/((2*i)+1);
    }
    return output*4;
}

//this pi algorithm isn't my work
double spigotpi(long digits){
    
    long LEN = (digits / 4 + 1) * 14;  //nec. array length
    
    long array[LEN];                   //array of 4-digit-decimals
    long b;                        //nominator prev. base
    long c = LEN;                  //index
    long d;                        //accumulator and carry
    long e = 0;                    //save previous 4 digits
    long f = 10000;                //new base, 4 decimal digits
    long g;                        //denom previous base
    long h = 0;                    //init switch
    
        for (; (b = c -= 14) > 0;) {    //outer loop: 4 digits per loop
            for (; --b > 0;) {      //inner loop: radix conv
                d *= b;            //acc *= nom prev base
                if (h == 0)
                    d += 2000 * f;   //first outer loop
                else
                    d += array[b] * f;   //non-first outer loop
                g = b + b - 1;           //denom prev base
                array[b] = d % g;
                d /= g;            //save carry
            }
            h =  (("%d"), e + d / f);//print prev 4 digits
            std::cout << h;
            d = e = d % f;         //save currently 4 digits
                                   //assure a small enough d
        }
        if (getchar())
        {
            return 0;
        }
        return 0;

    
}
double pi2(long precision){
    double output = 1;
    for (long k = precision;k>0;k--){
        output *= 1.0+(k/((2.0*k)+1.0));
        cout << output << endl;
    }
    return output;
}

double e(int precision = math::def_prec){
    double output = 0;
    for(int i=0;i<=precision;i++){
        output += 1/factorial1(i);
    }
    return output;
}

double abs(complex x, int precision){return x.magnitude(precision);}

template <class t = double>
double abs(t x){return (x<0)? -x : x;}


complex sqrt(complex x, int precision = math::def_prec){
    if (x.get_im()!=0){
    complex output((1/sqrt(2,precision))*sqrt(x.magnitude()+x.get_re(),precision),((x.get_im()/abs(x.get_im()))/sqrt(2,precision))*sqrt(x.magnitude()-x.get_re(),precision));
    return output;
    } else {
        complex output(sqrt(x.get_re()),0);
        return output;
    }
}

template <class t = double>
t quadratic(t a,t b, t c,t& otherroot, int precision = math::def_prec){
    otherroot = (-b+sqrt((b*b)-4*a*c,precision))/(2*a);
    return (-b-sqrt((b*b)-4*a*c,precision))/(2*a);
}

complex e_to_xi(double x, int precision = math::def_prec){
   complex output(cos(x,precision),sin(x,precision));
    return output;
}

template <class t>
t e_to_x(t x, int precision = math::def_prec){
    t output = 0;
    for(int i=0;i<=precision;i++){
        output += exp(x,i)/factorial1(i);
    }
    return output;
}

complex e_to_x(complex x, int precision = math::def_prec){return e_to_xi(x.get_im(),precision)*e_to_x(x.get_re(),precision);}

complex exp(complex base, complex exponent, int precision = math::def_prec){return e_to_x(ln(base)*exponent);}

template <class t = double>
t exp(t base, double exponent, int precision = math::def_prec){
    fraction asfrac(0,0);
    asfrac = tofrac(exponent);
    return exp(root(base,asfrac.get_den(),precision),asfrac.get_num());
}

template <class t = double>
t zeta(t s, int precision = math::def_prec){
    t output(0,0);
    for(int n=1.0;n<=precision;n++){output += 1/exp(n,s);}
    return output;
}

template <class t = double>
t eta(t s, int precision = math::def_prec){
    t output = 0;
    for(int n=1;n<=precision;n++){output += exp(-1,(n-1)%2)/exp(n,s);}
    return output;
}


template <class c=double>
c rsum(double (*f)(c), int n,c x[n],c t[n-1]){
    c output = 0;
    for(int i=0;i<=n-1;i++){
        output += f(t[i])*(x[i+1]-x[i]);
    }
    return output;
}

template <class t=double>
t inf(t a, t b){return (a<b)? a:b;}

template <class t=double>
t sup(t a, t b){return (a>b)? a:b;}

template <class c=double>
c ldsum(double (*f)(c), int n,c x[n]){
    c output = 0;
    for(int i=0;i<=n-1;i++){
        output += f(inf(x[i],x[i+1]))*(x[i+1]-x[i]);
    }
    return output;
}

template <class c=double>
c udsum(double (*f)(c), int n,c x[n]){
    c output = 0;
    for(int i=0;i<=n-1;i++){
        output += f(sup(x[i],x[i+1]))*(x[i+1]-x[i]);
    }
    return output;
}
template <class c=double>
c rint(c a, c b, c (*f)(c), int precision=math::def_prec){
    c partition[precision];
    for(int i=0;i<=precision;i++){
        partition[i] = a+(((b-a)/precision)*(i+1));
    }
    return rsum(f, precision, partition, partition); //use the lower bounds as tags
}

template <class c=double>
c udint(c a, c b, c (*f)(c), int precision = math::def_prec){
    c P[precision];
}

Что я должен улучшить, изменить, добавить и т. Д.?

4 ответа
4

using namespace std

Учитывая, что вы используете template и объявите свой class, вы работаете в библиотеке только для заголовков. Однако вам следует никогда использовать using namespace в заголовке. Не только это плохая практика, но в этом случае вы можете столкнуться с std::complex.

Удивительно operator семантика

Кроме того, почти все операторы не следуют обычным правилам и поэтому либо вводят в заблуждение, либо даже явно ошибаются. Давай посмотрим на operator-():

    complex operator - (complex a){
        r -= a.get_re();
        i -= a.get_im();
        return (*this);
}

Давай остановимся. Этот оператор изменяет значение текущего элемента. Это совершенно неожиданно. Если бы у меня был

complex a = foo();
complex b = bar();
complex a_cp = a;
complex b_cp = b;
complex c = a - b;

тогда я ожидаю a == a_cp а также b == b_cp. Однако это не так. Это также мешает нам использовать operator-() на consвысокие ценности.

Дублирование кода

Далее мы видим много дублирования. Остановимся на этом операторе и посмотрим на его вариант составного присваивания:

    complex operator - (double a){r -= a;return (*this);}

    void operator -= (complex a){
        r -= a.get_re();
        i -= a.get_im();
    }
    void operator -= (double a){r -= a;}

Мы сразу видим, что между - а также -=. Это почти тот же код. Обычно все операции составного присваивания возвращаются T& вместо void или же T, так что давайте сначала исправим это:

    complex& operator -= (complex a){
        r -= a.get_re();
        i -= a.get_im();
        return *this;
    }

Скотт Мейерс: предпочитаю функции, не являющиеся членами и не являющиеся друзьями

Сейчас operator- можно переписать с помощью operator-=, который также есть в книге Скотта Мейерса:

complex operator-(const complex &a, const complex &b) {
    complex tmp = a;
    tmp -= b;
    return tmp;
}

Обратите внимание, что указанный выше оператор является автономным. Это не функция-член. Это позволяет нам использовать его с типами, которые можно преобразовать в complex.

Предоставьте конструкторы преобразования

В приведенном выше коде отсутствует double или даже double[2] вариант. Однако это не проблема, если у нас есть подходящий конструктор, например

class complex {
...
public:
     complex(double r)    : r(r), i(0) {}
     complex(double c[2]) : r(c[0]), i(c[1]) {}

     complex(const complex&) = default;
     complex(complex&&)      = default;
     complex& operator=(const complex&) = default;
     complex& operator=(complex&&) = default;
}

Пока мы это делаем, мы также должны добавить обычные конструкторы копирования и перемещения, соответственно. задания.

Использовать const& где возможно

Почти все функции в complex скопируйте их аргумент. Им это не нужно.

Использовать const для функций, которые не могут изменить внутреннюю структуру класса

В get_re() а также get_im() функции не меняют complex среди прочего, отметьте их как const:

class complex {
...
    double get_re() const {
        return r;
    }
    double get_im() const {
        return i;
    }
...

Используйте средство форматирования кода

Форматирование кода и стиль слишком сильно меняются. Придерживайтесь удобочитаемый стиль и применять его вручную (что требует времени) или с помощью такого инструмента, как clang-format. Ваша IDE / редактор может иметь встроенную функцию.

  • Спасибо. Это много значит.

    — NY может

  • 4

    @NYcan Пожалуйста. Надеюсь, этот ответ не покажется слишком резким. В вашей реализации много замечательных вещей, но я считаю, что вам нужно прочитать недавнюю книгу по C ++ с лучшими практиками, чтобы ваш код сиял. Но не волнуйтесь, мы все прошли через этот этап в нашей карьере C ++ :).

    — Зета

  • 1

    В то время как const & обычно предпочтительнее, с типами этот маленький const & часто работает медленнее, но в любом случае все это должно быть встроено и, таким образом, в конечном итоге не будет иметь никакого значения.

    — Джек Эйдли

  • Будут ли учитываться принципы и практика программирования Бьярна Страуструпа с использованием C ++ (второе издание)? Он был опубликован в 2014 году.

    — NY может

  • Когда я попытался предоставить конструкторы преобразования для моего класса дробей, он сказал, что это неоднозначно, когда я назначаю одну дробь другой.

    — NY может

В дополнение к другому ответу,

  1. Поместите весь свой код в пространство имен. Почему внутри только константы math пространство имен? Также, math может быть слишком распространенным именем для использования в качестве пространства имен.

Хорошая идея могла бы выглядеть примерно так:

namespace my_math_library
{
    namespace constants
    {
       // put constants here
    }

    // rest of the stuff here
} 

  1. Вы можете объявить все константы, функции и классы как constexpr.

  1. template<class t = double>. Фактически вам не нужно указывать аргумент шаблона по умолчанию, поскольку аргумент шаблона будет просто выведен из аргумента функции.

  1. Почему некоторые из ваших функций являются шаблонными, а другие нет? (nsqrt, esqrt)

  1. for(double x = 1; x != num; x += h). Сравнение двух чисел с плавающей запятой с использованием == или же != это рецепт катастрофы. Сравните это с допуском (эпсилон), например std::numeric_limits<double>::epsilon()

  1. x-int(x) * precision. Будьте осторожны с приоритетом операторов. int(x) * precision оценивается в первую очередь. Для любого значения x> = ~ 4 вы вызываете поведение undefined, поскольку значение не помещается в int.

  1. t derivate(double (*f)
template<typename Func, typename T>
T derivate(Func&& func, T x)

Таким образом, вы также можете передавать лямбды и функторы.


  1. rsum(double (*f)(c), int n,c x[n],c t[n-1] не действует. Кроме того, используйте последовательные обозначения при использовании шаблонов.

template <class c=double>
c udint(c a, c b, c (*f)(c), int precision = math::def_prec){
    c P[precision];
}

Массивы переменной длины недопустимы в C ++. Кроме того, функция ничего не делает.

  • Когда я пытаюсь превратить свою сложную функцию exp в constexpr, она выдает ошибку «Тип возврата функции Constexpr 'math_lib :: complex' не является буквальным типом» (math_lib - это новое пространство имен, которое я создал.)

    - NY может

  • Ваш complex конструктор должен быть отмечен constexpr.

    - Риш

  • Когда я меняю его на constexpr class complex, он говорит, что класс не может быть помечен как constexpr.

    - NY может

  • Отметьте конструктор как constexpr. Как это: godbolt.org/z/hrG9WPGxP

    - Риш

  • ой, ой, я еще не доделал функцию udint. Вот почему он ничего не делает.

    - NY может

Примечание: это не исчерпывающий отчет, я также добавляю его к предыдущим обзорам.

Числовые константы

Взгляни на числовые константы в C ++. Вам не нужно определять их самостоятельно.

Опасная итерация

    double h = num/precision;
    double y = 1;
    for (double x=1;x!=num;x+=h){
        y += h*(1/(2*y));
    }

Не гарантируется, что x когда-либо достигнет стоимости num. Здесь вы можете использовать стиль итерации for (double x=1; x<=num; x+=h). Однако это похоже на n-е приближение к корню sqare, так почему бы не использовать целочисленный счетчик? В конце концов, вы не используете x вообще внутри петли!

Необычная итерация

  for (int i=1;i<=precision;i++){

Идиоматический способ записи - начать с нуля, а затем сравнить на равенство:

  for (int i=0; i!=precision; ++i) {

Единичные тесты

Это не только часть проверки вашего кода. Тем не менее, для таких библиотек использование Test-Driven Development (TDD), вероятно, является стандартом. Я считаю, что некоторые ошибки в вашем коде стали бы очевидными.

Документация

Что nsqrt(), esqrt() а также sqrt()? Зачем добавлять еще sqrt() перегрузка вообще? Почему бы не использовать тот, который предоставляется стандартной библиотекой? Это всего лишь примеры того, что вам нужно поместить в подходящий контекст, чтобы другие могли использовать ваш код.

Кроме того, если вы отложите несколько своих мыслей, это упростит сопровождение кода для вашего будущего «сейчас + 6 месяцев». Кроме того, иногда вы обнаруживаете, что первоначальная мысль на самом деле является глупой идеей, когда вы объясняете ее кому-то другому. То же самое произойдет, если вы запишете это как комментарий (что должно быть и для кого-то другого).

В дополнение к другим ответам вот несколько кратких комментариев по математической части.

Отказ от ответственности: я не числовой аналитик, поэтому, вероятно, есть более существенная эффективность, чем те, которые я предлагаю. Я бы сказал, что эффективность важна, поскольку обычно математический код является основой числового приложения.

Сначала несколько общих принципов:

  • Как правило, вы хотите попытаться свести к минимуму количество выполняемых вами операций, чтобы свести к минимуму ошибки.

  • Ошибки действительно существуют. Скомпилированная программа на C хранит свои числа в двоичном формате (или, если хотите, в шестнадцатеричном формате). Это имеет неожиданные последствия. Типичным примером является то, что число 1/10 (один больше десяти) является периодическим по основанию 2. Следствием этого является то, что оно сохраняется только как приближение, и что любое вычисление, связанное с этим числом, будет переносить ошибки. Это происходит слева и справа с любой обычной реализацией с плавающей запятой. Очень часто это не создает проблем; но в библиотеке вы хотите минимизировать вычисления, как упоминалось выше.

  • Любая функция, реализованная в вашей библиотеке, вероятно, будет многократно использоваться программой. Это делает эффективность очень важной.

Теперь к (очень) неполному списку конкретных проявлений вышеизложенного:

  • Функция Exp просто умножает base, exponent раз. Это можно немного упростить (по времени), сохранив продукты, и этого можно достичь, используя двоичное представление экспоненты. Например, если вы хотите сделать exp(base, 20), вы заметите, что 20 равно 2 ^ 4 + 2 ^ 2. Так ты можешь сделать (base^2)^4 + (base^2)^2. Это заставляет вас вычислять три продукта и одну сумму вместо 20 продуктов.

  • Вы определенно захотите выполнить сложное умножение в полярной форме, где произведения превращаются в суммы.

  • Предел и производная рассчитываются как (f(x+h)+f(x-h))/2 а также ((f(x+h)-f(x))/h+(f(x-h)-f(x))/(-h))/2. Численно это плохо. В численных расчетах вы хотите избегать очень маленьких знаменателей, так как это значительно увеличивает ваши ошибки.

  • В синусе и косинусе вы повторно вычисляете факториал для каждого члена ряда, вместо того, чтобы вычислять его по ходу: в члене n + 1 делитель равен d_{n+1} = d_n * (n+1).

  • Для касательной функции лучше использовать ее собственный ряд Тейлора, чем вычислять как синусы, так и косинусы и делить.

  • Комментирование структуры кода OP может быть полезно, но почти все используемые численные методы бесполезны для любого реального приложения, и решение этой проблемы здесь, IMO, не по теме.

    - алефзеро

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *