From 2ee06a30e175d4571139b6f90b5571b7aec8f1a0 Mon Sep 17 00:00:00 2001 From: Piyush Sethia Date: Tue, 2 Oct 2018 04:01:29 +0530 Subject: [PATCH 1/2] Adds complex class and fast fourier transform --- mathx/algebra.h | 169 +++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 147 insertions(+), 22 deletions(-) diff --git a/mathx/algebra.h b/mathx/algebra.h index 5281de5..b4cee07 100644 --- a/mathx/algebra.h +++ b/mathx/algebra.h @@ -8,22 +8,107 @@ namespace mathx { + double const PI = acos(-1); + + class Complex + { + double imaginary, real; + public: + //Constructor + Complex(double a = 0, double b = 0) : real(a), imaginary(b) { } + //Operators + Complex operator +(const Complex &o) const { + return Complex(this->real + o.real, this->imaginary + o.imaginary); + } + Complex operator -(const Complex &o) const { + return Complex(this->real - o.real, this->imaginary - o.imaginary); + } + Complex operator *(const Complex &o) const { + return Complex(this->real * o.real - this->imaginary * o.imaginary, this->real * o.imaginary + this->imaginary * o.real); + } + Complex operator *(const double &o) const { + return Complex(this->real * o, this->imaginary * o); + } + Complex operator /(const double &o) const { + return Complex(this->real / o, this->imaginary / o); + } + Complex operator /=(const double &o) { + this->real /= o; this->imaginary /= o; + return *this; + } + Complex &operator *=(const Complex &o) { + double r = this->real; + this->real = this->real * o.real - this->imaginary * o.imaginary; + this->imaginary = r * o.imaginary + this->imaginary * o.real ; + return (*this); + } + friend ostream &operator<<(ostream &output, const Complex &o ) { + output << "(" << o.real << ", " << o.imaginary << ")" ; + return output; + } + //Accessors + double Real() const { + return this->real; + } + double Imag() const { + return this->imaginary; + } + + }; + class Polynomial { private: std::string var; //variable For eg. "x","y" - float coefficients[100]; + int max_coefficient; + std::vector coefficients; + + void fft(std::vector& a, bool invert) //sub part which must be abstracted and must not be public + { + unsigned n = a.size(); + for (unsigned i = 0, j = 0; i < n; i++) { + + if (i > j) + swap(a[i], a[j]); + + for (unsigned l = n >> 1; (j ^= l) < l; l >>= 1); + + } + for (unsigned len = 2; len <= n; len <<= 1) { + + double ang = (double)(invert ? -2 : 2) * PI / (double)(len); + Complex wlen((double)cos(ang), (double)sin(ang)); + for (unsigned i = 0; i < n; i += len) { + + Complex w(1); + + for (unsigned j = 0; j < (len >> 1); j++) { + + Complex u = a[i + j], v = a[i + j + (len >> 1)] * w; + + a[i + j] = u + v; + a[i + j + (len >> 1)] = u - v; + w *= wlen; + } + } + } + + if (invert) + for (int i = 0; i < n; ++i) + a[i] /= (double)(n); + } public: Polynomial () { var = "x"; - for(int i = 0 ; i < 100; i++) - coefficients[i] = 0; + coefficients.resize(1); + coefficients[0] = 0; + max_coefficient = 0; } Polynomial (std::string x) { var = "x"; - for(int i = 0 ; i < 100; i++) - coefficients[i] = 0; + max_coefficient = 0; + coefficients[0] = 0; parse(x); } @@ -32,38 +117,45 @@ namespace mathx { bool flag = false; int temp_power, start = 0; float temp_co; - for(int i = 0; i < exp.length(); ) + for (int i = 0; i < exp.length(); ) { start = i; - if(exp[i] == '-' || exp[i] == '+') + if (exp[i] == '-' || exp[i] == '+') i++; - while(exp[i] >= 48 && exp[i] <= 57 ) + while (exp[i] >= 48 && exp[i] <= 57 ) i++; - if ( exp.substr(start,i-start) == "+" || exp.substr(start,i-start) == "-") { - temp_co = std::stof(exp.substr(start,i-start) + "1"); + if ( exp.substr(start, i - start) == "+" || exp.substr(start, i - start) == "-") { + temp_co = std::stof(exp.substr(start, i - start) + "1"); } else - temp_co = std::stof(exp.substr(start,i-start)); - while(!(exp[i] >= 48 && exp[i] <= 57)) + temp_co = std::stof(exp.substr(start, i - start)); + while (!(exp[i] >= 48 && exp[i] <= 57)) i++; start = i; - while(exp[i] >= 48 && exp[i] <= 57 ) + while (exp[i] >= 48 && exp[i] <= 57 ) i++; - temp_power = std::stof(exp.substr(start,i-start)); + temp_power = std::stof(exp.substr(start, i - start)); + if (temp_power > max_coefficient) + { + max_coefficient = temp_power; + coefficients.resize(temp_power + 1); + } coefficients[temp_power] = temp_co; } } float evaluate (float const x) { float sum = 0; - for(int i = 0; i<100; i++) { + for (int i = 0; i <= max_coefficient; i++) { if (coefficients[i] != 0) - sum += std::pow(x,i)*coefficients[i]; + sum += std::pow(x, i) * coefficients[i]; } return sum; } float const getCoeff (int i) { + if (i > max_coefficient) + return static_cast(0); return coefficients[i]; } @@ -73,17 +165,22 @@ namespace mathx { Polynomial derivative () { Polynomial der; - for(int i = 0; i < 99; i++) { - der.coefficients[i] = (i+1)*coefficients[i+1]; - } + der.coefficients.resize(max_coefficient); + for (int i = 0; i < max_coefficient; i++) + der.coefficients[i] = (i + 1) * coefficients[i + 1]; return der; } std::string const display () { std::string exp = ""; - for (int i = 99; i >= 0; i--) { - if(coefficients[i] != 0) { - exp += "+("; + bool first = true; + for (int i = max_coefficient; i >= 0; i--) { + if (coefficients[i] != 0) { + if (!first) + exp += "+"; + else + first = false; + exp += "("; exp += std::to_string(coefficients[i]); exp += ")x^("; exp += std::to_string(i); @@ -92,6 +189,34 @@ namespace mathx { } return exp; } + + Polynomial operator *(const Polynomial &b) const + { + unsigned n = 1; + + while (n < this->max_coefficient + b.max_coefficient + 1) + n <<= 1; + + std::vector fa(this->coefficients.begin(), this->coefficients.end()), fb(b.coefficients.begin(), b.coefficients.end()); + fa.resize(n); + fb.resize(n); + + fft(fa, false); + fft(fb, false); + + for (unsigned i = 0; i < n; ++i) + fa[i] *= fb[i]; + + fft(fa, true); + + Polynomial res; + res.max_coefficient = this->max_coefficient + b.max_coefficient; + + for (unsigned i = 0; i <= res.max_coefficient; ++i) + res.coefficients[i] = round(fa[i].Real()); + + return res; + } }; } From e282154abf5792de0aa503172836ba7050acf797 Mon Sep 17 00:00:00 2001 From: Piyush Sethia Date: Tue, 2 Oct 2018 04:13:20 +0530 Subject: [PATCH 2/2] Minor bug fixes --- mathx/algebra.h | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/mathx/algebra.h b/mathx/algebra.h index b4cee07..a8bde3c 100644 --- a/mathx/algebra.h +++ b/mathx/algebra.h @@ -60,7 +60,7 @@ namespace mathx { private: std::string var; //variable For eg. "x","y" int max_coefficient; - std::vector coefficients; + std::vector coefficients; void fft(std::vector& a, bool invert) //sub part which must be abstracted and must not be public { @@ -107,6 +107,7 @@ namespace mathx { Polynomial (std::string x) { var = "x"; + coefficients.resize(1); max_coefficient = 0; coefficients[0] = 0; parse(x); @@ -116,7 +117,7 @@ namespace mathx { trim(exp); bool flag = false; int temp_power, start = 0; - float temp_co; + double temp_co; for (int i = 0; i < exp.length(); ) { start = i; @@ -125,16 +126,16 @@ namespace mathx { while (exp[i] >= 48 && exp[i] <= 57 ) i++; if ( exp.substr(start, i - start) == "+" || exp.substr(start, i - start) == "-") { - temp_co = std::stof(exp.substr(start, i - start) + "1"); + temp_co = std::stod(exp.substr(start, i - start) + "1"); } else - temp_co = std::stof(exp.substr(start, i - start)); + temp_co = std::stod(exp.substr(start, i - start)); while (!(exp[i] >= 48 && exp[i] <= 57)) i++; start = i; while (exp[i] >= 48 && exp[i] <= 57 ) i++; - temp_power = std::stof(exp.substr(start, i - start)); + temp_power = std::stod(exp.substr(start, i - start)); if (temp_power > max_coefficient) { max_coefficient = temp_power; @@ -144,8 +145,8 @@ namespace mathx { } } - float evaluate (float const x) { - float sum = 0; + double evaluate (double const x) { + double sum = 0; for (int i = 0; i <= max_coefficient; i++) { if (coefficients[i] != 0) sum += std::pow(x, i) * coefficients[i]; @@ -153,9 +154,9 @@ namespace mathx { return sum; } - float const getCoeff (int i) { + double const getCoeff (int i) { if (i > max_coefficient) - return static_cast(0); + return static_cast(0); return coefficients[i]; }