Я создавал математический движок для 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 ответа
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 / редактор может иметь встроенную функцию.
В дополнение к другому ответу,
- Поместите весь свой код в пространство имен. Почему внутри только константы
math
пространство имен? Также,math
может быть слишком распространенным именем для использования в качестве пространства имен.
Хорошая идея могла бы выглядеть примерно так:
namespace my_math_library
{
namespace constants
{
// put constants here
}
// rest of the stuff here
}
- Вы можете объявить все константы, функции и классы как
constexpr
.
template<class t = double>
. Фактически вам не нужно указывать аргумент шаблона по умолчанию, поскольку аргумент шаблона будет просто выведен из аргумента функции.
- Почему некоторые из ваших функций являются шаблонными, а другие нет? (
nsqrt
,esqrt
)
for(double x = 1; x != num; x += h)
. Сравнение двух чисел с плавающей запятой с использованием==
или же!=
это рецепт катастрофы. Сравните это с допуском (эпсилон), напримерstd::numeric_limits<double>::epsilon()
x-int(x) * precision
. Будьте осторожны с приоритетом операторов.int(x) * precision
оценивается в первую очередь. Для любого значения x> = ~ 4 вы вызываете поведение undefined, поскольку значение не помещается вint
.
t derivate(double (*f)
template<typename Func, typename T>
T derivate(Func&& func, T x)
Таким образом, вы также можете передавать лямбды и функторы.
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)
.Для касательной функции лучше использовать ее собственный ряд Тейлора, чем вычислять как синусы, так и косинусы и делить.
Спасибо. Это много значит.
— NY может
@NYcan Пожалуйста. Надеюсь, этот ответ не покажется слишком резким. В вашей реализации много замечательных вещей, но я считаю, что вам нужно прочитать недавнюю книгу по C ++ с лучшими практиками, чтобы ваш код сиял. Но не волнуйтесь, мы все прошли через этот этап в нашей карьере C ++ :).
— Зета
В то время как const & обычно предпочтительнее, с типами этот маленький const & часто работает медленнее, но в любом случае все это должно быть встроено и, таким образом, в конечном итоге не будет иметь никакого значения.
— Джек Эйдли
Будут ли учитываться принципы и практика программирования Бьярна Страуструпа с использованием C ++ (второе издание)? Он был опубликован в 2014 году.
— NY может
Когда я попытался предоставить конструкторы преобразования для моего класса дробей, он сказал, что это неоднозначно, когда я назначаю одну дробь другой.
— NY может