Journal: When I coded a Multivariate Polynomial Operation Tool in C++ (discont. Journal)

I did not have much experience in C++ like I do in python. I started this as a helper tool for Dr. Anirban Nag who only has been very warm and welcoming for anytime I got stuck! Though I could not kepe up with the journalling I still kept up with the codebase and it is still up in my GitHub and I am still updating it. Only caveat is that the progress is intentionally slow because of my office work load!

In this project I went from absolute zero in C++ to handling complicated situations like namespaces and object oriented operator overloading.

1. I have no idea how to setup c++ compiler
2. I use linux so I installed g++ with sudo
3. Created the utter necessary directories /src, /tests, /include
4. Rapidly study, Struct, Classes, Namespaces
5. inside src i created a new atom.h
6. inside i made a total simple class

class Coefficient {
    public:
        float num;
};

7. Then I called it using main.cpp and set a new coefficient and build it with g++ to get a *.out file. Fun lol.
8. Now need to implement proper atom.
9. I create the proper atom class

class Atom {
    public:
        float coeff;
        int exp;
        char var;
};

10. And love how it worked for $3x^2$

int main() {
    Atom a;
    a.coeff = 3;
    a.var = 'x';
    a.exp = 2;

    std::cout << a.coeff << a.var << "^" << a.exp << std::endl;

    return 0;
}

11. Now for simplest polynomial I need to be able to find the value of the above when $x=2$, I add a method

float value(float x) {
    return coeff * pow(x, exp);
}

12. Keep in mind, just the constant is just, the variable with exponent as 0.
13. Now I try to implement a quadratic polynomial $3x^2+5x-2$, then evaluate at $x=2$

int main() {
    Atom a(3, 'x', 2);
    Atom b(5, 'x', 1);
    Atom c(-2, 'x', 0);

    Atom poly[3] = {a, b, c};

    float sol = 0;
    for (int i = 0; i<3 ;i++) {
        sol = sol + poly[i].eval(2);
    };
    
    return 0;
}

14. I need to put together $30x^2yz$ where it is like combination of three or more atoms
15. One simple way is to make an array of atoms.
16. After implementing this, adding a eval method was too easy.
17. Now I need to add a complexity like i need it to understand $2x^2\times 3x = 6x^3$

Particle operator*(const Particle& p2) {
            Particle p3;
            float coeff = 1;
            int count = 0;
            std::map<char, int> variables;

            // calculating coeff
            for (int i = 0; i < this->count; i++) {
                Atom atom = this->atoms[i];
                char var = atom.var;

                coeff *= atom.coeff;

                if (variables.find(var) == variables.end()) {
                    variables[var] = atom.exp;
                    count++;
                } else { 
                    variables[var] += atom.exp;                 
                }
            }
            for (int i = 0; i < p2.count; i++) {
                Atom atom = p2.atoms[i];
                char var = atom.var;
                coeff *= p2.atoms[i].coeff;

                if (variables.find(var) == variables.end()) {
                    variables[var] = atom.exp;
                    count++;
                } else { 
                    variables[var] += atom.exp;                 
                }
            }

            p3.count = count;

            int i = 0;
            for (std::pair<char, int> pair : variables) {
                p3.atoms[i] = Atom(i == 0 ? coeff : 1, pair.first, pair.second);
                i++;
            }
            
            return p3;
        }

18. What I am doing in the above method is, first I multiply all the coefficient of all individual atoms together and save in coeff
19. Then I make a map, whenever I see a same variable i just increment the exponent has i get the {var: exp} pair.
20. Then I just iterate each var as a new atom and add it to our result Particle’s list.
21. Now I need to know how I can add up 2 different particle $2xy + 3xy = 5xy$.
22. But before that I need to understand how to even check if those two particle are addable or not.
23. What is coming to tip of my head is, I cant turn every atom I have into a map as an attribute of a particle. if I see two particle have same map then it can be added by coefficient

void mapVars () {
	coefficient = 1;
	variables.clear();
	
	for (int i = 0; i < this->count; i++) {
		Atom atom = this->atoms[i];
		char var = atom.var;
		if (atom.exp == 0) {
			var = '@';
		}

		coeff *= atom.coeff;
		
		if (variables.find(var) == variables.end()) {
			variables[var] = atom.exp;
		} else { 
			variables[var] += atom.exp;                 
		}
	}
}

24. It is working fine, but we notice something, if the two maps are same then cool, but if it doesn’t then we cant just send an empty particle rather we need to send a polynomial
25. Got to make the Polynomial class now.
26. It is just array of particles.
27. Now we need to add up two polynomial $(x+1) + (2x+3)$. If we just stack up the polynomials it will become -> $(x+1+2x+3)$.
28. But I previously implemented to add up particles if it matches save {var: exp} map
29. So we check for individual elements. When they match up we add, otherwise we stack up like -> $(3x + 4)$

Polynomial operator+(const Polynomial& P2) {
	Polynomial P;
	
	// copying the array
	for (int i = 0; i < this->count; i++) {
		P.polynom[i] = this->polynom[i];
	}
	P.count = this->count;

	for (int i = 0; i < P2.count ; i++) {
		bool match = false;
		Particle p2 = P2.polynom[i];

		for (int j = 0; j < P.count; j++) {
			Particle p1 = P.polynom[j];

			if (p2.variables == p1.variables) {
				Particle p3 = p1 + p2;
				P.polynom[j] = p3;
				match = true;
				break;
			}
		}

		if (!match) {
			P.polynom[P.count] = p2;
			P.count++;
		}
	}

	return P;
}

30. I have been using hard sized arrays to implement the idea, but I have introduced and replaced ass array to dynamically sized arrays.
31. For multiplication for initial we can just do brute force multiplication

Polynomial operator*(const Polynomial& P2) {
	Polynomial P;

	for (int i = 0; i < this->polynom.size(); i++) {
		for (int j = 0; j < P2.polynom.size(); j++) {
			Polynomial temp;
			temp.addParticle(this->polynom[i] * P2.polynom[j]);
			P = P + temp;
		}
	}

	return P;
}

32. I have to define the temp Polynomial because of two reason:
    • I can’t perform Particle + Polynomial addition
    • IF I do not do polynomial addition x, 2x will be kept separately instead of as 3x
33. Now I need to work on the division of polynomials
34. for simple cases its trivial so I am targeting little bit complex with long division method
35. before that I need to quickly implement the lexicography ordering for polynomial class.
36. Fortunately the variable names follow lexicographic order for particles because of std::map
37. now I need to implement the ordering in Polynomial class.
38. First we check if two particle when compared which one is greater in Lex Ordering:

bool isLexGreater(const Particle& p2) {
	std::vector<float> V = this->getExpVec();
	std::vector<float> W = p2.getExpVec();

	std::vector<float> resultant;

	int maxLen = std::max(V.size(), W.size());

	for (int i = 0; i < maxLen; i++) {
		float v = (i < V.size()) ? V[i] : 0;
		float w = (i < W.size()) ? W[i] : 0;

		resultant.push_back(v - w);
	}

	for (int i = 0; i < resultant.size(); i++) {
		if (resultant[i] != 0) {
			return resultant[i] > 0;
		}
	}
	return false;
}

39. Now every time we want to use division algorithm in the polynomial we just compare and order them with the following:

void sortLexicographic() {
	std::sort(polynom.begin(), polynom.end(), 
		[](const Particle& a, const Particle& b) {
			return a.isLexGreater(b);
		});
}

39. Notice how I am not implementing it during adding new particle every time -> too slow + not needed for other operations.
40. Now onto implementing division
41. I do ground up, so I implement atom division first:

// * Overload / operator, A = a / b
Atom operator/(const Atom& b) {
	float coeff = this->coeff / b.coeff;
	float exp = this->exp - b.exp;

	if (this->var == b.var) {
		return Atom(coeff, this->var, exp);
	}
	return Atom();
}

42. Now I implement particle division:

// * Overload / operator -> p = p1 / p2
Particle operator/(const Particle& p2) {
	std::map<char, float> vars;
	
	float coeff = this->coefficient / p2.coefficient;

	for (const auto& pair: this->variables) {
		char var = pair.first;
		float exp1 = pair.second;
		float exp2 = (p2.variables.count(var) > 0) ? p2.variables.at(var) : 0;
		vars[var] = exp1 - exp2;
	}

	for (const auto& pair : p2.variables) {
		char var = pair.first;
		if (this->variables.count(var) == 0) {
			vars[var] = 0 - pair.second;
		}
	}

	Particle p;
	int i = 0;
	for (const auto pair: vars) {
		p.addAtom(Atom(
			i == 0 ? coeff : 1,
			pair.first,
			pair.second
		));
		i++;
	}

	return p;
}

42. Now onto implementing the sophisticated division algorithm where if we want to divide P with Q then we get, $P = Q*D + R$ , The famous Remainder theorem.
43. After they are sorted Lexicographical, the need is that every variable of leading particle of D must exist in leading particle of P and with exponents >= in lead of P.

// * Overload / Operator -> P = Q*D + R
Polynomial operator/(const Polynomial& divisor) {
	std::cout << "Division started!" << std::endl;
	Polynomial Q;           // Q = 0 -> Quotient
	Polynomial R;           // R = 0 -> Remainder
	Polynomial P = *this;   // P -> copy of this polynom
	Polynomial D = divisor; // D -> copy of divisor

	P.sortLexicographic();
	D.sortLexicographic();

	while (P.polynom.size() > 0) {

		Particle leadP = P.polynom[0];
		Particle leadD = D.polynom[0];
		bool divisible = true;

		for (const auto& pair : leadD.variables) {
			char var = pair.first;
			if (leadP.variables.count(var) > 0) {
				divisible = leadP.variables[var] >= pair.second;
			} else { divisible = false; }

			if (!divisible) {break;}
		}

		if (divisible) {
			Polynomial q;
			q.addParticle(leadP / leadD);
			Q = Q + q;
			P = P - (q * D);
		} else {
			Polynomial leadTerm;
			leadTerm.addParticle(leadP);
			R = R + leadTerm;
			P = P - leadTerm;
		}
		
	}

	return Q;
}

44. Testing multiple polynomial is getting trickier to define them manually
45. The first thing I need to do is add a string2polynom implementation