1#!/usr/bin/env python 2 3import polypy 4import polypy_test 5 6import random 7 8polypy_test.init() 9 10def check_factorization(p, debug=False, check_product = True): 11 global polypy_test 12 if debug: 13 print("check_factorization:") 14 print("p = {0}".format(p)) 15 # Do the factorization 16 factorization = p.factor() 17 C, factors = factorization[0], factorization[1:] 18 if debug: 19 print("C = {0}".format(C)) 20 print("factors = {0}".format(factors)) 21 # The constant should always be positive 22 if C <= 0: 23 print("Wrong factorization (constant)") 24 print("p = {0}".format(p)) 25 print("C = {0}".format(C)) 26 print("factors = {0}".format(factors)) 27 polypy_test.check(False) 28 return 29 # Check if factorization multiplies to the input 30 if check_product: 31 product = C 32 for (f, d) in factors: 33 if (f.degree() == 0): 34 print("Wrong factorization (constant factor)") 35 print("p = {0}".format(p)) 36 print("factors = {0}".format(factors)) 37 polypy_test.check(False) 38 product = product * f**d 39 if (p != product): 40 print("Wrong factorization (product mismatch)") 41 print("p = {0}".format(p)) 42 print("product = {0}".format(product)) 43 print("factors = {0}".format(factors)) 44 polypy_test.check(False) 45 return 46 # Done, we're OK 47 polypy_test.check(True) 48 49 50""" 51 Returns a list of cyclotomic polynomials F1...Fn 52""" 53def cyclotomic(n): 54 if (n == 1): 55 P1 = polypy.x - 1 56 return [P1] 57 else: 58 # Get up to 1 polynomials 59 L = cyclotomic(n-1) 60 # Compute Pn 61 Pn = polypy.x**n - 1 62 for d, Pd in enumerate(L): 63 if (n % (d+1) == 0): 64 Pn = Pn / Pd 65 L.append(Pn) 66 return L 67 68polypy_test.start("Factorization in Z (Regressions)") 69 70# polypy.trace_enable("factorization") 71# polypy.trace_enable("hensel") 72# polypy.trace_enable("arithmetic") 73 74x = polypy.x 75 76p = (x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5); 77check_factorization(p) 78 79p = 1*x**4 + 94*x**3 + 107*x**2 + 771*x + (-690) 80check_factorization(p) 81 82p= x**4 + 1 83check_factorization(p) 84 85p = x**8 - x**6 + x**4 - x**2 + 1 86check_factorization(p) 87 88p = (x**2 + 1)*(x**2 + 3) 89check_factorization(p) 90 91p = x**2 + 1 92check_factorization(p) 93 94p = (x**2 + 1)*(x**2 + 2) 95check_factorization(p) 96 97p = (x - 1)*(x**2 + 1) 98check_factorization(p) 99 100p = (x - 1)*(x - 2) 101check_factorization(p) 102 103polypy_test.start("Berlekamp in Z_13 (Knuth)") 104 105K = polypy.CoefficientRing(13) 106x = polypy.x.to_ring(K); 107 108p = x**8 + x**6 + 10*x**4 + 10*x**3 + 8*x**2 + 2*x + 8 109 110check_factorization(p) 111 112polypy_test.start("Square-free modular (simple)") 113 114K = polypy.CoefficientRing(13) 115x = polypy.x.to_ring(K); 116 117p = x**13 118check_factorization(p); 119 120p = x**13 - x 121check_factorization(p) 122 123p = x*(x-1)**2*(x-2)**3 124check_factorization(p) 125 126p = 1*x**6 + 3*x**5 + (-5*x**4) + 4*x**3 + 6*x**2 + 6*x + (-1) 127check_factorization(p) 128 129polypy_test.start("Square-free modular (random)") 130 131for factors_count in range(1, 11): 132 for test in range(1, 10): 133 # Constant 134 p = 1 135 # Factors 136 deg = random.randint(1, 2) 137 for k in range(0, factors_count): 138 degree = random.randint(1, 3) 139 f_k = polypy_test.random_upolynomial(K, degree = 3, M = 10, lc = 1) 140 p *= f_k**deg 141 deg += random.randint(1, 2) 142 # Check it 143 check_factorization(p) 144 145polypy_test.start("Cyclotomic (Z)") 146 147#for p in cyclotomic(100): 148# print(p) 149 150# http://mathworld.wolfram.com/Swinnerton-DyerPolynomial.html 151polypy_test.start("Swinnerton-Dyer (Z)") 152