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