diff --git a/src/main/python/algorithms/math/ChineseRemainderTheorem.py b/src/main/python/algorithms/math/ChineseRemainderTheorem.py new file mode 100644 index 0000000..c42cc11 --- /dev/null +++ b/src/main/python/algorithms/math/ChineseRemainderTheorem.py @@ -0,0 +1,186 @@ +''' + * @author (original JAVA) William Fiset, william.alexandre.fiset@gmail.com + * (conversion to Python) Hiruna Vishwamith, hvishwamith@gmail.com + * @date 20 Nov 2022 + + Use the chinese remainder theorem to solve a set of congruence equations. + +
The first method (eliminateCoefficient) is used to reduce an equation of the form cx≡a(mod + m)cx≡a(mod m) to the form x≡a_new(mod m_new)x≡anew(mod m_new), which gets rids of the + coefficient. A value of null is returned if the coefficient cannot be eliminated. + +
The second method (reduce) is used to reduce a set of equations so that the moduli become + pairwise co-prime (which means that we can apply the Chinese Remainder Theorem). The input and + output are of the form x≡a_0(mod m_0),...,x≡a_n-1(mod m_n-1)x≡a_0(mod m_0),...,x≡a_n-1(mod + m_n-1). Note that the number of equations may change during this process. A value of null is + returned if the set of equations cannot be reduced to co-prime moduli. + +
The third method (crt) is the actual Chinese Remainder Theorem. It assumes that all pairs of
+ moduli are co-prime to one another. This solves a set of equations of the form x≡a_0(mod
+ m_0),...,x≡v_n-1(mod m_n-1)x≡a_0(mod m_0),...,x≡v_n-1(mod m_n-1). It's output is of the form
+ x≡a_new(mod m_new)x≡a_new(mod m_new).
+
+
+'''
+
+import math
+import queue
+import random
+
+# eliminateCoefficient() takes cx?a(mod m) and gives x?a_new(mod m_new).
+
+def eliminateCoefficient( c, a, m) :
+ d = egcd(c, m)[0]
+ if (a % d != 0) :
+ return None
+ c /= d
+ a /= d
+ m /= d
+ inv = egcd(c, m)[1]
+ m = abs(m)
+ a = (((a * inv) % m) + m) % m
+ return [a, m]
+# reduce() takes a set of equations and reduces them to an equivalent
+# set with pairwise co-prime moduli (or null if not solvable).
+
+def reduce( a, m) :
+ aNew = []
+ mNew = []
+ # Split up each equation into prime factors
+ i = 0
+ while (i < len(a)) :
+ factors = primeFactorization(m[i])
+ factors.sort()
+ iterator = factors.listIterator()
+ while (iterator.hasNext()) :
+ val = iterator.next()
+ total = val
+ while (iterator.hasNext()) :
+ nextVal = iterator.next()
+ if (nextVal == val) :
+ total *= val
+ else :
+ iterator.previous()
+ break
+ aNew.add(a[i] % total)
+ mNew.add(total)
+ i += 1
+ # Throw away repeated information and look for conflicts
+ i = 0
+ while (i < aNew.size()) :
+ j = i + 1
+ while (j < aNew.size()) :
+ if (mNew.get(i) % mNew.get(j) == 0 or mNew.get(j) % mNew.get(i) == 0) :
+ if (mNew.get(i) > mNew.get(j)) :
+ if ((aNew.get(i) % mNew.get(j)) == aNew.get(j)) :
+ aNew.remove(j)
+ mNew.remove(j)
+ j -= 1
+ continue
+ else :
+ return None
+ else :
+ if ((aNew.get(j) % mNew.get(i)) == aNew.get(i)) :
+ aNew.remove(i)
+ mNew.remove(i)
+ i -= 1
+ break
+ else :
+ return None
+ j += 1
+ i += 1
+ # Put result into an array
+ res = [[0] * (aNew.size()) for _ in range(2)]
+ i = 0
+ while (i < aNew.size()) :
+ res[0][i] = aNew.get(i)
+ res[1][i] = mNew.get(i)
+ i += 1
+ return res
+
+def crt( a, m) :
+ M = 1
+ i = 0
+ while (i < len(m)) :
+ M *= m[i]
+ i += 1
+ inv = [0] * (len(a))
+ i = 0
+ while (i < len(inv)) :
+ inv[i] = egcd(M / m[i], m[i])[1]
+ i += 1
+ x = 0
+ i = 0
+ while (i < len(m)) :
+ x += (M / m[i]) * a[i] * inv[i]
+ # Overflow could occur here
+ x = ((x % M) + M) % M
+ i += 1
+ return [x, M]
+
+def primeFactorization( n) :
+ factors = []
+ if (n <= 0) :
+ raise Exception("IllegalArgumentException")
+ elif(n == 1) :
+ return factors
+ divisorQueue = queue.PriorityQueue()
+ divisorQueue.add(n)
+ while (not divisorQueue.isEmpty()) :
+ divisor = divisorQueue.remove()
+ if (isPrime(divisor)) :
+ factors.append(divisor)
+ continue
+ next_divisor = pollardRho(divisor)
+ if (next_divisor == divisor) :
+ divisorQueue.add(divisor)
+ else :
+ divisorQueue.add(next_divisor)
+ divisorQueue.add(divisor / next_divisor)
+ return factors
+
+def pollardRho( n) :
+ if (n % 2 == 0) :
+ return 2
+ # Get a number in the range [2, 10^6]
+ x = 2 + int((999999 * random.random()))
+ c = 2 + int((999999 * random.random()))
+ y = x
+ d = 1
+ while (d == 1) :
+ x = (x * x + c) % n
+ y = (y * y + c) % n
+ y = (y * y + c) % n
+ d = gcf(abs(x - y), n)
+ if (d == n) :
+ break
+ return d
+# Extended euclidean algorithm
+
+def egcd( a, b) :
+ if (b == 0) :
+ return [a, 1, 0]
+ else :
+ ret = egcd(b, a % b)
+ tmp = ret[1] - ret[2] * (a / b)
+ ret[1] = ret[2]
+ ret[2] = tmp
+ return ret
+
+def gcf( a, b) :
+ return a if b == 0 else gcf(b, a % b)
+
+def isPrime( n) :
+ if (n < 2) :
+ return False
+ if (n == 2 or n == 3) :
+ return True
+ if (n % 2 == 0 or n % 3 == 0) :
+ return False
+ limit = int(math.sqrt(n))
+ i = 5
+ while (i <= limit) :
+ if (n % i == 0 or n % (i + 2) == 0) :
+ return False
+ i += 6
+ return True
\ No newline at end of file
diff --git a/src/main/python/algorithms/math/CompressedPrimeSieve.py b/src/main/python/algorithms/math/CompressedPrimeSieve.py
new file mode 100644
index 0000000..162ca91
--- /dev/null
+++ b/src/main/python/algorithms/math/CompressedPrimeSieve.py
@@ -0,0 +1,78 @@
+'''
+ * @author (original JAVA) William Fiset, william.alexandre.fiset@gmail.com
+ * (conversion to Python) Hiruna Vishwamith, hvishwamith@gmail.com
+ * @date 20 Nov 2022
+
+ Generate a compressed prime sieve using bit manipulation. The idea is that each bit represents a
+ boolean value indicating whether a number is prime or not. This saves a lot of room when creating
+ the sieve. In this implementation I store all odd numbers in individual longs meaning that for
+ each long I use I can represent a range of 128 numbers (even numbers are omitted because they are
+ not prime, with the exception of 2 which is handled as a special case).
+ Time Complexity: ~O(nloglogn)
+
+'''
+import math
+
+
+NUM_BITS = 128.0
+NUM_BITS_SHIFT = 7
+# 2^7 = 128
+# Sets the bit representing n to 1 indicating this number is not prime
+
+def setBit( arr, n) :
+ if ((n & 1) == 0) :
+ return
+ # n is even
+ arr[n >> NUM_BITS_SHIFT] |= 1 << ((n - 1) >> 1)
+# Returns true if the bit for n is off (meaning n is a prime).
+# Note: do use this method to access numbers outside your prime sieve range!
+
+def isNotSet( arr, n) :
+ if (n < 2) :
+ return False
+ # n is not prime
+ if (n == 2) :
+ return True
+ # two is prime
+ if ((n & 1) == 0) :
+ return False
+ # n is even
+ chunk = arr[n >> NUM_BITS_SHIFT]
+ mask = 1 << ((n - 1) >> 1)
+ return (chunk & mask) != mask
+# Returns true/false depending on whether n is prime.
+
+def isPrime( sieve, n) :
+ return isNotSet(sieve, n)
+# Returns an array of longs with each bit indicating whether a number
+# is prime or not. Use the isNotSet and setBit methods to toggle to bits for each number.
+
+def primeSieve( limit) :
+ numChunks = int(math.ceil(limit / NUM_BITS))
+ sqrtLimit = int(math.sqrt(limit))
+ # if (limit < 2) return 0; // uncomment for primeCount purposes
+ # int primeCount = (int) Math.ceil(limit / 2.0); // Counts number of primes <= limit
+ chunks = [0] * (numChunks)
+ chunks[0] = 1
+ # 1 as not prime
+ i = 3
+ while (i <= sqrtLimit) :
+ if (isNotSet(chunks, i)) :
+ j = i * i
+ while (j <= limit) :
+ if (isNotSet(chunks, j)) :
+ setBit(chunks, j)
+ j += i
+ i += 2
+ return chunks
+# Example usage.
+
+
+limit = 200
+sieve = primeSieve(limit)
+i = 0
+while (i <= limit) :
+ if (isPrime(sieve, i)) :
+ print("%d is prime!\n" % (i),end ="",sep ="")
+ i += 1
+
diff --git a/src/main/python/algorithms/math/EulerTotientFunction.py b/src/main/python/algorithms/math/EulerTotientFunction.py
new file mode 100644
index 0000000..eae7b7c
--- /dev/null
+++ b/src/main/python/algorithms/math/EulerTotientFunction.py
@@ -0,0 +1,75 @@
+import math
+import random
+import queue
+
+def eulersTotient( n) :
+ for p in set() :
+ n -= (n / p)
+ return n
+
+def primeFactorization( n) :
+ factors = []
+ if (n <= 0) :
+ raise Exception("IllegalArgumentException")
+ elif(n == 1) :
+ return factors
+ divisorQueue = queue.PriorityQueue()
+ divisorQueue.add(n)
+ while (not divisorQueue.isEmpty()) :
+ divisor = divisorQueue.remove()
+ if (isPrime(divisor)) :
+ factors.append(divisor)
+ continue
+ next_divisor = pollardRho(divisor)
+ if (next_divisor == divisor) :
+ divisorQueue.add(divisor)
+ else :
+ divisorQueue.add(next_divisor)
+ divisorQueue.add(divisor / next_divisor)
+ return factors
+
+def pollardRho( n) :
+ if (n % 2 == 0) :
+ return 2
+ # Get a number in the range [2, 10^6]
+ x = 2 + int((999999 * random.random()))
+ c = 2 + int((999999 * random.random()))
+ y = x
+ d = 1
+ while (d == 1) :
+ x = (x * x + c) % n
+ y = (y * y + c) % n
+ y = (y * y + c) % n
+ d = gcf(abs(x - y), n)
+ if (d == n) :
+ break
+ return d
+
+def gcf( a, b) :
+ return a if b == 0 else gcf(b, a % b)
+
+def isPrime( n) :
+ if (n < 2) :
+ return False
+ if (n == 2 or n == 3) :
+ return True
+ if (n % 2 == 0 or n % 3 == 0) :
+ return False
+ limit = int(math.sqrt(n))
+ i = 5
+ while (i <= limit) :
+ if (n % i == 0 or n % (i + 2) == 0) :
+ return False
+ i += 6
+ return True
+
+def main( args) :
+ # Prints 8 because 1,2,4,7,8,11,13,14 are all
+ # less than 15 and relatively prime with 15
+ print("phi(15) = %d\n" % (eulersTotient(15)),end ="",sep ="")
+ print()
+ x = 1
+ while (x <= 11) :
+ print("phi(%d) = %d\n" % (x,eulersTotient(x)),end ="",sep ="")
+ x += 1
+
diff --git a/src/main/python/algorithms/math/EulerTotientFunctionWithSieve.py b/src/main/python/algorithms/math/EulerTotientFunctionWithSieve.py
new file mode 100644
index 0000000..a1d9787
--- /dev/null
+++ b/src/main/python/algorithms/math/EulerTotientFunctionWithSieve.py
@@ -0,0 +1,70 @@
+'''
+ * @author (original JAVA) William Fiset, william.alexandre.fiset@gmail.com
+ * (conversion to Python) Hiruna Vishwamith, hvishwamith@gmail.com
+ * @date 20 Nov 2022
+
+ Computes Euler's totient function
+
+'''
+
+import math
+
+def sieve( limit) :
+ if (limit <= 2) :
+ return [0] * (0)
+
+ numPrimes = int((1.25506 * limit / math.log(float(limit))))
+ primes = [0] * (numPrimes)
+ index = 0
+ isComposite = [False] * (limit)
+ sqrtLimit = int(math.sqrt(limit))
+ i = 2
+ while (i <= sqrtLimit) :
+ if (not isComposite[i]) :
+ primes[index] = i
+ index += 1
+ j = i * i
+ while (j < limit) :
+ isComposite[j] = True
+ j += i
+ i += 1
+ i = sqrtLimit + 1
+ while (i < limit) :
+ if (not isComposite[i]) :
+ primes[index] = i
+ index += 1
+ i += 1
+ return primes[:index]
+
+
+MAX = 1000000
+PRIMES = sieve(MAX)
+
+
+def totient(n) :
+ if (n >= MAX - 1) :
+ raise Exception("MAX not large enough!")
+ ans = n
+ i = 1
+ p = PRIMES[0]
+ while (p * p <= n) :
+ if (n % p == 0) :
+ ans -= int(ans / p)
+ while (n % p == 0) :
+ n /= p
+ p = PRIMES[i]
+ i += 1
+ # Last factor
+ if (n != 1) :
+ ans -= int(ans / n)
+ return ans
+
+
+
+# Prints 8 because 1,2,4,7,8,11,13,14 are all
+# less than 15 and relatively prime with 15
+print(totient(15))
+x = 1
+while (x <= 11) :
+ print("phi(%d) = %d\n" % (x,totient(x)),end ="",sep ="")
+ x += 1
diff --git a/src/main/python/algorithms/math/ExtendedEuclideanAlgorithm.py b/src/main/python/algorithms/math/ExtendedEuclideanAlgorithm.py
new file mode 100644
index 0000000..d8bf156
--- /dev/null
+++ b/src/main/python/algorithms/math/ExtendedEuclideanAlgorithm.py
@@ -0,0 +1,23 @@
+'''
+ * @author (original JAVA) William Fiset, william.alexandre.fiset@gmail.com
+ * (conversion to Python) Hiruna Vishwamith, hvishwamith@gmail.com
+ * @date 20 Nov 2022
+
+ This function performs the extended euclidean algorithm on two numbers a and b.
+ The function returns the gcd(a,b) as well as the numbers x and y such
+ that ax + by = gcd(a,b). This calculation is important in number theory
+ and can be used for several things such as finding modular inverses and
+ solutions to linear Diophantine equations.
+ Time Complexity ~O(log(a + b))
+
+'''
+
+def egcd(a,b):
+ if b==0:
+ return {a,1,0}
+ else:
+ ret=egcd(b,a%b)
+ tmp=ret[1]-ret[2]*(a/b)
+ ret[1]=ret[2]
+ ret[2]=tmp
+ return ret
diff --git a/src/main/python/algorithms/math/FastFourierTransformComplexNumber.py b/src/main/python/algorithms/math/FastFourierTransformComplexNumber.py
new file mode 100644
index 0000000..9704a05
--- /dev/null
+++ b/src/main/python/algorithms/math/FastFourierTransformComplexNumber.py
@@ -0,0 +1,62 @@
+'''
+ * @author (original JAVA) William Fiset, william.alexandre.fiset@gmail.com
+ * (conversion to Python) Hiruna Vishwamith, hvishwamith@gmail.com
+ * @date 20 Nov 2022
+
+ This snippet multiplies 2 complex polynomials very efficiently using the Fast Fourier Transform.
+ Time Complexity: O(nlogn)
+
+'''
+import math
+
+def fft(comp_list):
+ n=len(comp_list)
+ if (n==1):return comp_list[0]
+ arr=comp_list[n//2]
+ for k in range(n//2):
+ arr[k]=comp_list[2*k]
+ q=fft(arr)
+ for k in range(n//2):
+ arr[k]=comp_list[2*k+1]
+ r=fft(arr)
+ y=[1]*n
+ for k in range(n//2):
+ kth=-2*k*math.pi/n
+ wk=math.cos(kth)+1j*math.sin(kth)
+ y[k]=q[k]+wk*r[k]
+ y[k+n//2]=q[k]-wk*r[k]
+
+ return y
+
+def polyMult(a,b):
+ exp=math.ceil(math.log2(len(a)+len(b)-1))
+ length=int(math.pow(2,exp))
+ a=pad(a,length)
+ b=pad(b,length)
+ c,d,e=fft(a),fft(d),[1]*length
+ for i in range(length):
+ e[i]=c[i]*d[i]
+ return ifft(e)
+
+def pad(arr,n):
+ padded=[1]*n
+ for i in range(n):
+ if i