I've implemented Fraction and Polynomial classes to simulate the mathematical operations that may be applied to them. The Fraction class was mainly implemented so that numbers like 2 / 3 will be represented as 2 / 3 rather than 0.666 only to be rounded off later.
Here's the Fraction class:
import math as m
class Fraction:
def __init__(self, num, denom=1):
if type(num) == float:
num = floatToFraction(num)
if type(denom) == float:
denom = floatToFraction(denom)
#num, denom must be of type int or Fraction
if type(num) == int and type(denom) == int:
#catches zero divison
if denom == 0:
raise ZeroDivisionError()
gcd = m.gcd(num, denom)
#always in simplest form
self.num = int(num / gcd)
self.denom = int(denom / gcd)
elif type(num) in [int, Fraction] and type(denom) in [int, Fraction]:
#catches zero division
if type(denom) == int:
if denom == 0:
raise ZeroDivisionError()
else:
if denom.getNum() == 0:
raise ZeroDivisionError()
new_f = num / denom
self.num = new_f.getNum()
self.denom = new_f.getDenom()
else:
raise TypeError('Constructor must receive inputs of type int, float or Fraction.')
#sign-correction
if self.num * self.denom > 0:
self.num, self.denom = abs(self.num), abs(self.denom)
else:
self.num, self.denom = -1 * abs(self.num), abs(self.denom)
def __repr__(self):
if self.denom == 1:
return f'{self.num}'
return f'( {self.num} / {self.denom} )'
def getNum(self):
return self.num
def getDenom(self):
return self.denom
def inverse(self):
return Fraction(self.denom, self.num)
def __add__(self, o_f):
if type(o_f) == int:
o_f = Fraction(o_f)
return Fraction(self.num * o_f.getDenom() + self.denom * o_f.getNum(), self.denom * o_f.getDenom())
def __radd__(self, o_f):
return self + o_f
def __sub__(self, o_f):
return self + (o_f * -1)
def __rsub__(self, o_f):
return (self - o_f) * -1
def __mul__(self, o_f):
if type(o_f) == int:
o_f = Fraction(o_f)
return Fraction(self.num * o_f.getNum(), self.denom * o_f.getDenom())
def __rmul__(self, o_f):
return self * o_f
def __truediv__(self, o_f):
#define for operations between numbers and fractions
if type(o_f) == int:
o_f = Fraction(o_f)
return self * o_f.inverse()
def __rtruediv__(self, o_f):
return Fraction(o_f) * self.inverse()
def floatToFraction(x):
#only converts non-recurring decimals, CANNOT handle recurring decimals e.g. 0.333333
d = 1
while x % 1 != 0:
x *= 10
d *= 10
return Fraction(int(x), int(d))
And here's the Polynomial class:
import numpy as np
from fraction import Fraction
class Polynomial:
"""
NOTE:
-> Constructor may accept a Polynomial, NumPy Array, List or Number.
-> Empty iterables will be replaced with an object of an equivalent type holding the number 0, for ease of computation.
"""
def __init__(self, ipt=[]):
#arr is the coefficient array describing the polynomial.
if isinstance(ipt, Polynomial):
#Creates a copy of an existing polynomial.
ipt = ipt.getArr()
elif type(ipt) == np.ndarray:
if ipt.size == 0:
ipt += np.array([Fraction(0)])
ipt = np.array([n if isinstance(n, Fraction) else Fraction(float(n)) for n in ipt])
elif type(ipt) == list:
if len(ipt) == 0:
ipt = [0]
ipt = np.array([Fraction(n) for n in ipt])
elif type(ipt) == int:
ipt = np.array([Fraction(ipt)])
else:
raise TypeError('Constructor must be a Polynomial, NumPy Array, List or Number.')
self.coef = np.array(ipt)
self.dim = self.coef.size - 1
def __getitem__(self, idx):
return self.coef[idx]
def __setitem__(self, n, coef):
self.coef[n] = coef
def __iter__(self):
self.index = 0
return self
def __next__(self):
if self.index <= self.dim:
coef = self[self.index]
self.index += 1
return coef
raise StopIteration
def __repr__(self):
p_str = ''
first = True
for i in reversed(list(range(self.coef.size))):
if self[i] != 0 or i == 0:
p_str += ('+ ' if not first else '') + str(self[i]) + (f' * x^{str(i)} ' if i != 0 else '')
first = False
return p_str
def getArr(self):
return self.coef
def getDim(self):
return self.dim
def _expand(self, n):
#NOTE: Expansion will not affect the original array.
k = n - self.dim
exp_cpy = np.append(self.coef, np.zeros(k))
return Polynomial(exp_cpy)
"""
NOTE: Arithmetic operations accept constants and Polynomials only.
"""
def __add__(self, P2):
if not isinstance(P2, Polynomial):
P2 = Polynomial([P2])
if P2.getDim() < self.dim:
return P2 + self
else:
return Polynomial(self._expand(P2.getDim()).getArr() + P2.getArr())
def __sub__(self, P2):
return self + (P2 * -1)
def __mul__(self, P2):
if type(P2) in [int, float, Fraction]:
#Wraps it with a Polynomial object.
P2 = Polynomial([P2])
elif P2.getDim() < 0:
P2 += Polynomial([0])
product = Polynomial([])
s1, s2 = self.dim + 1, P2.getDim() + 1
for i in range(s2):
temp = {}
p_e = P2[i]
for j in range(s1):
temp[i+j] = self[j] * p_e #Single pass per element.
new_p = Polynomial(np.zeros(max(temp.keys())+1))
for k in temp:
new_p[k] = temp[k]
product += new_p
return product
def __rmul__(self, P2):
return self * P2
def compute(self, x):
val = 0
for i in range(self.dim + 1):
val += self[i] * x ** i
return val
I'm trying to implement a concept called perturbation theory in applied mathematics. The details of this aren't important, what I'm trying to compute are the coefficients for the polynomials that can be used to calculate the first x terms of k-powers. For example, for k=2, the sum of the first x terms of the series: 1^2, 2^2, 3^2, 4^2, ... , x^2 can be calculated using (1 / 3) * x^3 + (1 / 2) * x^2 + (1 / 6) * x. Such a polynomial exists for all k. I've written some code that computes the coefficients and memoizes the results as follows:
from fraction import Fraction
pn_bank = {}
def comb(x, y):
#Computes xCy.
v = Fraction(1)
while y > 0:
v *= Fraction(x, y)
x -= 1
y -= 1
return v
def binom(n):
#Computes the polynomial coefficients for (x + 1)^n.
coef = []
for i in range(n + 1):
coef.append(comb(n, i))
return coef
def perturb(n):
c_arr = pn_bank.get(n)
if c_arr:
return c_arr
coefs = binom(n + 1)
coefs[n + 1] = 0
rhs = Polynomial(1)
for i in range(0, n):
#Skips avoiding perturb(n + 1) and perturb(n) and triggering infinite recursion.
rhs += perturb(i) * coefs[i]
lhs = Polynomial(binom(n + 1)) - rhs
lhs *= Fraction(1, coefs[n])
return lhs
When I tried computing the coefficients for k=17 and upwards, the coefficients being computed are correct, but for some reason, the polynomial class's compute method isn't working correctly (at least that's where I think the problem lies). For example, the sum of the first 10 terms for k=17 is 19179318935377305 (according to Wolfram Alpha). According to a JavaScript for loop, it's 19179318935377304 and according to my code, it's 19179318935377305.666 (which is impossible since the sum of integers has to be an integer). I checked for k=18,19,etc. and this only occurs for the 9th term and onwards. Everything before is still accurate. Is there some floating-point round off error here? What might be going on here.
Here's the output as requested for k=0,1,2,...,19:
POWER: 0
POLYNOMIAL: 1 * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 2
SUM of 3rd term: 3
SUM of 10th term: 10
-------------------------------------------------------
POWER: 1
POLYNOMIAL: ( 1 / 2 ) * x^2 + ( 1 / 2 ) * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 3
SUM of 3rd term: 6
SUM of 10th term: 55
-------------------------------------------------------
POWER: 2
POLYNOMIAL: ( 1 / 3 ) * x^3 + ( 1 / 2 ) * x^2 + ( 1 / 6 ) * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 5
SUM of 3rd term: 14
SUM of 10th term: 385
-------------------------------------------------------
POWER: 3
POLYNOMIAL: ( 1 / 4 ) * x^4 + ( 1 / 2 ) * x^3 + ( 1 / 4 ) * x^2 + 0 * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 9
SUM of 3rd term: 36
SUM of 10th term: 3025
-------------------------------------------------------
POWER: 4
POLYNOMIAL: ( 1 / 5 ) * x^5 + ( 1 / 2 ) * x^4 + ( 1 / 3 ) * x^3 + 0 * x^2 + ( -1 / 30 ) * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 17
SUM of 3rd term: 98
SUM of 10th term: 25333
-------------------------------------------------------
POWER: 5
POLYNOMIAL: ( 1 / 6 ) * x^6 + ( 1 / 2 ) * x^5 + ( 5 / 12 ) * x^4 + 0 * x^3 + ( -1 / 12 ) * x^2 + 0 * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 33
SUM of 3rd term: 276
SUM of 10th term: 220825
-------------------------------------------------------
POWER: 6
POLYNOMIAL: ( 1 / 7 ) * x^7 + ( 1 / 2 ) * x^6 + ( 1 / 2 ) * x^5 + 0 * x^4 + ( -1 / 6 ) * x^3 + 0 * x^2 + ( 1 / 42 ) * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 65
SUM of 3rd term: 794
SUM of 10th term: 1978405
-------------------------------------------------------
POWER: 7
POLYNOMIAL: ( 1 / 8 ) * x^8 + ( 1 / 2 ) * x^7 + ( 7 / 12 ) * x^6 + 0 * x^5 + ( -7 / 24 ) * x^4 + 0 * x^3 + ( 1 / 12 ) * x^2 + 0 * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 129
SUM of 3rd term: 2316
SUM of 10th term: 18080425
-------------------------------------------------------
POWER: 8
POLYNOMIAL: ( 1 / 9 ) * x^9 + ( 1 / 2 ) * x^8 + ( 2 / 3 ) * x^7 + 0 * x^6 + ( -7 / 15 ) * x^5 + 0 * x^4 + ( 2 / 9 ) * x^3 + 0 * x^2 + ( -1 / 30 ) * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 257
SUM of 3rd term: 6818
SUM of 10th term: 167731333
-------------------------------------------------------
POWER: 9
POLYNOMIAL: ( 1 / 10 ) * x^10 + ( 1 / 2 ) * x^9 + ( 3 / 4 ) * x^8 + 0 * x^7 + ( -7 / 10 ) * x^6 + 0 * x^5 + ( 1 / 2 ) * x^4 + 0 * x^3 + ( -3 / 20 ) * x^2 + 0 * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 513
SUM of 3rd term: 20196
SUM of 10th term: 1574304985
-------------------------------------------------------
POWER: 10
POLYNOMIAL: ( 1 / 11 ) * x^11 + ( 1 / 2 ) * x^10 + ( 5 / 6 ) * x^9 + 0 * x^8 + -1 * x^7 + 0 * x^6 + 1 * x^5 + 0 * x^4 + ( -1 / 2 ) * x^3 + 0 * x^2 + ( 5 / 66 ) * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 1025
SUM of 3rd term: 60074
SUM of 10th term: 14914341925
-------------------------------------------------------
POWER: 11
POLYNOMIAL: ( 1 / 12 ) * x^12 + ( 1 / 2 ) * x^11 + ( 11 / 12 ) * x^10 + 0 * x^9 + ( -11 / 8 ) * x^8 + 0 * x^7 + ( 11 / 6 ) * x^6 + 0 * x^5 + ( -11 / 8 ) * x^4 + 0 * x^3 + ( 5 / 12 ) * x^2 + 0 * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 2049
SUM of 3rd term: 179196
SUM of 10th term: 142364319625
-------------------------------------------------------
POWER: 12
POLYNOMIAL: ( 1 / 13 ) * x^13 + ( 1 / 2 ) * x^12 + 1 * x^11 + 0 * x^10 + ( -11 / 6 ) * x^9 + 0 * x^8 + ( 22 / 7 ) * x^7 + 0 * x^6 + ( -33 / 10 ) * x^5 + 0 * x^4 + ( 5 / 3 ) * x^3 + 0 * x^2 + ( -691 / 2730 ) * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 4097
SUM of 3rd term: 535538
SUM of 10th term: 1367428536133
-------------------------------------------------------
POWER: 13
POLYNOMIAL: ( 1 / 14 ) * x^14 + ( 1 / 2 ) * x^13 + ( 13 / 12 ) * x^12 + 0 * x^11 + ( -143 / 60 ) * x^10 + 0 * x^9 + ( 143 / 28 ) * x^8 + 0 * x^7 + ( -143 / 20 ) * x^6 + 0 * x^5 + ( 65 / 12 ) * x^4 + 0 * x^3 + ( -691 / 420 ) * x^2 + 0 * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 8193
SUM of 3rd term: 1602516
SUM of 10th term: 13202860761145
-------------------------------------------------------
POWER: 14
POLYNOMIAL: ( 1 / 15 ) * x^15 + ( 1 / 2 ) * x^14 + ( 7 / 6 ) * x^13 + 0 * x^12 + ( -91 / 30 ) * x^11 + 0 * x^10 + ( 143 / 18 ) * x^9 + 0 * x^8 + ( -143 / 10 ) * x^7 + 0 * x^6 + ( 91 / 6 ) * x^5 + 0 * x^4 + ( -691 / 90 ) * x^3 + 0 * x^2 + ( 7 / 6 ) * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 16385
SUM of 3rd term: 4799354
SUM of 10th term: 128037802953445
-------------------------------------------------------
POWER: 15
POLYNOMIAL: ( 1 / 16 ) * x^16 + ( 1 / 2 ) * x^15 + ( 5 / 4 ) * x^14 + 0 * x^13 + ( -91 / 24 ) * x^12 + 0 * x^11 + ( 143 / 12 ) * x^10 + 0 * x^9 + ( -429 / 16 ) * x^8 + 0 * x^7 + ( 455 / 12 ) * x^6 + 0 * x^5 + ( -691 / 24 ) * x^4 + 0 * x^3 + ( 35 / 4 ) * x^2 + 0 * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 32769
SUM of 3rd term: 14381676
SUM of 10th term: 1246324856379625
-------------------------------------------------------
POWER: 16
POLYNOMIAL: ( 1 / 17 ) * x^17 + ( 1 / 2 ) * x^16 + ( 4 / 3 ) * x^15 + 0 * x^14 + ( -14 / 3 ) * x^13 + 0 * x^12 + ( 52 / 3 ) * x^11 + 0 * x^10 + ( -143 / 3 ) * x^9 + 0 * x^8 + ( 260 / 3 ) * x^7 + 0 * x^6 + ( -1382 / 15 ) * x^5 + 0 * x^4 + ( 140 / 3 ) * x^3 + 0 * x^2 + ( -3617 / 510 ) * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 65537
SUM of 3rd term: 43112258
SUM of 10th term: ( 206902004244165440 / 17 )
-------------------------------------------------------
POWER: 17
POLYNOMIAL: ( 1 / 18 ) * x^18 + ( 1 / 2 ) * x^17 + ( 17 / 12 ) * x^16 + 0 * x^15 + ( -17 / 3 ) * x^14 + 0 * x^13 + ( 221 / 9 ) * x^12 + 0 * x^11 + ( -2431 / 30 ) * x^10 + 0 * x^9 + ( 1105 / 6 ) * x^8 + 0 * x^7 + ( -11747 / 45 ) * x^6 + 0 * x^5 + ( 595 / 3 ) * x^4 + 0 * x^3 + ( -3617 / 60 ) * x^2 + 0 * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 131073
SUM of 3rd term: 129271236
SUM of 10th term: ( 1072613870418395776 / 9 )
-------------------------------------------------------
POWER: 18
POLYNOMIAL: ( 1 / 19 ) * x^19 + ( 1 / 2 ) * x^18 + ( 3 / 2 ) * x^17 + 0 * x^16 + ( -34 / 5 ) * x^15 + 0 * x^14 + 34 * x^13 + 0 * x^12 + ( -663 / 5 ) * x^11 + 0 * x^10 + ( 1105 / 3 ) * x^9 + 0 * x^8 + ( -23494 / 35 ) * x^7 + 0 * x^6 + 714 * x^5 + 0 * x^4 + ( -3617 / 10 ) * x^3 + 0 * x^2 + ( 43867 / 798 ) * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 262145
SUM of 3rd term: 387682634
SUM of 10th term: ( 22227014932144214016 / 19 )
-------------------------------------------------------
POWER: 19
POLYNOMIAL: ( 1 / 20 ) * x^20 + ( 1 / 2 ) * x^19 + ( 19 / 12 ) * x^18 + 0 * x^17 + ( -323 / 40 ) * x^16 + 0 * x^15 + ( 323 / 7 ) * x^14 + 0 * x^13 + ( -4199 / 20 ) * x^12 + 0 * x^11 + ( 4199 / 6 ) * x^10 + 0 * x^9 + ( -223193 / 140 ) * x^8 + 0 * x^7 + 2261 * x^6 + 0 * x^5 + ( -68723 / 40 ) * x^4 + 0 * x^3 + ( 43867 / 84 ) * x^2 + 0 * x^1 + 0
SUM of 1st term: 1
SUM of 2nd term: 524289
SUM of 3rd term: 1162785756
SUM of 10th term: 11506994510201251840
-------------------------------------------------------
numpy.polynomial.Polynomialclass useful: numpy.org/doc/stable/reference/routines.polynomials.htmlfractionsmodule of the standard library?fractionsto verify that your own implementation is working correctly.int(num / gcd)use the integer division operator:num // gcd. There may be cases where they get different results due to floating point roundoff errors.