[compalg] Les 3 en 4 (snelle vermenigvuldiging)

1866 days ago by J.Demeyer

Oefening 2 (DFT)

Algemene functies
def DFT(omega, f, n): R = parent(f) L = [f(omega^i) for i in range(n)] return R(L) 
       
def punt_product(f, g): R = parent(f) n = min(f.degree(), g.degree())+1 L = [f[i]*g[i] for i in range(n)] return R(L) 
       

Over een eindig veld
F = FiniteField(113) omega = F(40) n = 16 
       
P.<x> = F[] f = 256*x^8 - 448*x^6 + 224*x^4 - 32*x^2 + 1 g = x^7 + 12*x^5 - 35*x^3 + 1 
       
ff = DFT(omega, f, n) ff 
        
43*x^15 + 58*x^14 + 75*x^13 + 57*x^12 + 65*x^11 + 8*x^10 + 40*x^9 + x^8
+ 43*x^7 + 58*x^6 + 75*x^5 + 57*x^4 + 65*x^3 + 8*x^2 + 40*x + 1
43*x^15 + 58*x^14 + 75*x^13 + 57*x^12 + 65*x^11 + 8*x^10 + 40*x^9 + x^8 + 43*x^7 + 58*x^6 + 75*x^5 + 57*x^4 + 65*x^3 + 8*x^2 + 40*x + 1
gg = DFT(omega, g, n) gg 
        
40*x^15 + 8*x^14 + 2*x^13 + 13*x^12 + 3*x^11 + 102*x^10 + 35*x^9 +
23*x^8 + 75*x^7 + 107*x^6 + 102*x^4 + 112*x^3 + 13*x^2 + 80*x + 92
40*x^15 + 8*x^14 + 2*x^13 + 13*x^12 + 3*x^11 + 102*x^10 + 35*x^9 + 23*x^8 + 75*x^7 + 107*x^6 + 102*x^4 + 112*x^3 + 13*x^2 + 80*x + 92
ffgg = punt_product(ff, gg) ffgg 
        
25*x^15 + 12*x^14 + 37*x^13 + 63*x^12 + 82*x^11 + 25*x^10 + 44*x^9 +
23*x^8 + 61*x^7 + 104*x^6 + 51*x^4 + 48*x^3 + 104*x^2 + 36*x + 92
25*x^15 + 12*x^14 + 37*x^13 + 63*x^12 + 82*x^11 + 25*x^10 + 44*x^9 + 23*x^8 + 61*x^7 + 104*x^6 + 51*x^4 + 48*x^3 + 104*x^2 + 36*x + 92
DFT(omega^-1, ffgg, n)/n 
        
30*x^15 + 25*x^13 + 13*x^11 + 30*x^9 + 30*x^8 + 26*x^7 + 4*x^6 + 2*x^5 +
111*x^4 + 78*x^3 + 81*x^2 + 1
30*x^15 + 25*x^13 + 13*x^11 + 30*x^9 + 30*x^8 + 26*x^7 + 4*x^6 + 2*x^5 + 111*x^4 + 78*x^3 + 81*x^2 + 1

Controle:
f*g 
        
30*x^15 + 25*x^13 + 13*x^11 + 30*x^9 + 30*x^8 + 26*x^7 + 4*x^6 + 2*x^5 +
111*x^4 + 78*x^3 + 81*x^2 + 1
30*x^15 + 25*x^13 + 13*x^11 + 30*x^9 + 30*x^8 + 26*x^7 + 4*x^6 + 2*x^5 + 111*x^4 + 78*x^3 + 81*x^2 + 1

Over de gehele getallen
n = 16 
       
P.<x> = ZZ[] f = 256*x^8 - 448*x^6 + 224*x^4 - 32*x^2 + 1 g = x^7 + 12*x^5 - 35*x^3 + 1 
       
maxf = max(abs(c) for c in f.coefficients()) maxg = max(abs(c) for c in g.coefficients()) B = (maxf * maxg * (g.degree()+1)) * 2 + 1 B 
        
250881
250881
B += (1-B) % n B 
        
250881
250881
p = B while not p.is_prime(): p += n p 
        
250993
250993
F = FiniteField(p) prim = F.primitive_element() omega = prim^((p-1)//n) omega 
        
33673
33673
ff = DFT(omega, f.change_ring(F), n) ff 
        
19309*x^15 + 198519*x^14 + 222640*x^13 + 961*x^12 + 174628*x^11 +
52540*x^10 + 84389*x^9 + x^8 + 19309*x^7 + 198519*x^6 + 222640*x^5 +
961*x^4 + 174628*x^3 + 52540*x^2 + 84389*x + 1
19309*x^15 + 198519*x^14 + 222640*x^13 + 961*x^12 + 174628*x^11 + 52540*x^10 + 84389*x^9 + x^8 + 19309*x^7 + 198519*x^6 + 222640*x^5 + 961*x^4 + 174628*x^3 + 52540*x^2 + 84389*x + 1
gg = DFT(omega, g.change_ring(F), n) gg 
        
31057*x^15 + 228474*x^14 + 84283*x^13 + 9882*x^12 + 141110*x^11 +
222462*x^10 + 136670*x^9 + 23*x^8 + 219938*x^7 + 22521*x^6 + 166712*x^5
+ 241113*x^4 + 109885*x^3 + 28533*x^2 + 114325*x + 250972
31057*x^15 + 228474*x^14 + 84283*x^13 + 9882*x^12 + 141110*x^11 + 222462*x^10 + 136670*x^9 + 23*x^8 + 219938*x^7 + 22521*x^6 + 166712*x^5 + 241113*x^4 + 109885*x^3 + 28533*x^2 + 114325*x + 250972
ffgg = punt_product(ff, gg) ffgg 
        
57336*x^15 + 237955*x^14 + 28454*x^13 + 209861*x^12 + 17319*x^11 +
162449*x^10 + 65287*x^9 + 23*x^8 + 232275*x^7 + 159083*x^6 + 165833*x^5
+ 43054*x^4 + 80944*x^3 + 193624*x^2 + 103491*x + 250972
57336*x^15 + 237955*x^14 + 28454*x^13 + 209861*x^12 + 17319*x^11 + 162449*x^10 + 65287*x^9 + 23*x^8 + 232275*x^7 + 159083*x^6 + 165833*x^5 + 43054*x^4 + 80944*x^3 + 193624*x^2 + 103491*x + 250972
fg = DFT(omega^-1, ffgg, n)/n fg 
        
256*x^15 + 2624*x^13 + 236881*x^11 + 18336*x^9 + 256*x^8 + 242770*x^7 +
250545*x^6 + 1132*x^5 + 224*x^4 + 250958*x^3 + 250961*x^2 + 1
256*x^15 + 2624*x^13 + 236881*x^11 + 18336*x^9 + 256*x^8 + 242770*x^7 + 250545*x^6 + 1132*x^5 + 224*x^4 + 250958*x^3 + 250961*x^2 + 1
L = [c.centerlift() for c in fg] P(L) 
        
__main__:2: DeprecationWarning: centerlift is deprecated. Please use
lift_centered instead.
See http://trac.sagemath.org/15804 for details.
256*x^15 + 2624*x^13 - 14112*x^11 + 18336*x^9 + 256*x^8 - 8223*x^7 -
448*x^6 + 1132*x^5 + 224*x^4 - 35*x^3 - 32*x^2 + 1
__main__:2: DeprecationWarning: centerlift is deprecated. Please use lift_centered instead.
See http://trac.sagemath.org/15804 for details.
256*x^15 + 2624*x^13 - 14112*x^11 + 18336*x^9 + 256*x^8 - 8223*x^7 - 448*x^6 + 1132*x^5 + 224*x^4 - 35*x^3 - 32*x^2 + 1

Controle:
f*g 
        
256*x^15 + 2624*x^13 - 14112*x^11 + 18336*x^9 + 256*x^8 - 8223*x^7 -
448*x^6 + 1132*x^5 + 224*x^4 - 35*x^3 - 32*x^2 + 1
256*x^15 + 2624*x^13 - 14112*x^11 + 18336*x^9 + 256*x^8 - 8223*x^7 - 448*x^6 + 1132*x^5 + 224*x^4 - 35*x^3 - 32*x^2 + 1

Oefening 3 (Schönhage-Strassen)
def verdeel(pol, m, x, y): ret = pol.parent().zero() for i in range(pol.degree() + 1): xpow = i % m ypow = i // m ret += pol[i] * x^xpow * y^ypow return ret 
       
R = FiniteField(7) P.<x> = R[] f = 2*x^25 + x^2 + 4 g = 3*x^18 + 5*x^14 + 2*x^5 + x^3 + 2*x^2 
       
Q.<y> = P[] 
       
n = 64 m = 8 t = 8 
       
fs = verdeel(f,m,x,y) gs = verdeel(g,m,x,y) 
       
fs(x^m) 
        
2*x^25 + x^2 + 4
2*x^25 + x^2 + 4
gs(x^m) 
        
3*x^18 + 5*x^14 + 2*x^5 + x^3 + 2*x^2
3*x^18 + 5*x^14 + 2*x^5 + x^3 + 2*x^2
D = P.quotient(x^(2*m)+1) D 
        
Univariate Quotient Polynomial Ring in xbar over Finite Field of size 7
with modulus x^16 + 1
Univariate Quotient Polynomial Ring in xbar over Finite Field of size 7 with modulus x^16 + 1
eta = D(x^2) eta 
        
xbar^2
xbar^2
parent(fs) 
        
Univariate Polynomial Ring in y over Univariate Polynomial Ring in x
over Finite Field of size 7
Univariate Polynomial Ring in y over Univariate Polynomial Ring in x over Finite Field of size 7
fss = fs(eta*y) gss = gs(eta*y) fss 
        
2*xbar^7*y^3 + xbar^2 + 4
2*xbar^7*y^3 + xbar^2 + 4
parent(fss) 
        
Univariate Polynomial Ring in y over Univariate Quotient Polynomial Ring
in xbar over Finite Field of size 7 with modulus x^16 + 1
Univariate Polynomial Ring in y over Univariate Quotient Polynomial Ring in xbar over Finite Field of size 7 with modulus x^16 + 1
hss = (fss * gss) % (y^t - 1) hss 
        
6*xbar^13*y^5 + 3*xbar^15*y^4 + (4*xbar^12 + 2*xbar^10 + 4*xbar^9)*y^3 +
(3*xbar^8 + 5*xbar^6)*y^2 + (5*xbar^10 + 6*xbar^8)*y + 2*xbar^7 +
2*xbar^5 + 2*xbar^4 + 4*xbar^3 + xbar^2
6*xbar^13*y^5 + 3*xbar^15*y^4 + (4*xbar^12 + 2*xbar^10 + 4*xbar^9)*y^3 + (3*xbar^8 + 5*xbar^6)*y^2 + (5*xbar^10 + 6*xbar^8)*y + 2*xbar^7 + 2*xbar^5 + 2*xbar^4 + 4*xbar^3 + xbar^2
hs = hss(eta^-1 * y) hs 
        
6*xbar^3*y^5 + 3*xbar^7*y^4 + (4*xbar^6 + 2*xbar^4 + 4*xbar^3)*y^3 +
(3*xbar^4 + 5*xbar^2)*y^2 + (5*xbar^8 + 6*xbar^6)*y + 2*xbar^7 +
2*xbar^5 + 2*xbar^4 + 4*xbar^3 + xbar^2
6*xbar^3*y^5 + 3*xbar^7*y^4 + (4*xbar^6 + 2*xbar^4 + 4*xbar^3)*y^3 + (3*xbar^4 + 5*xbar^2)*y^2 + (5*xbar^8 + 6*xbar^6)*y + 2*xbar^7 + 2*xbar^5 + 2*xbar^4 + 4*xbar^3 + xbar^2
L = [c.lift() for c in list(hs)] haccent = Q(L) haccent 
        
6*x^3*y^5 + 3*x^7*y^4 + (4*x^6 + 2*x^4 + 4*x^3)*y^3 + (3*x^4 +
5*x^2)*y^2 + (5*x^8 + 6*x^6)*y + 2*x^7 + 2*x^5 + 2*x^4 + 4*x^3 + x^2
6*x^3*y^5 + 3*x^7*y^4 + (4*x^6 + 2*x^4 + 4*x^3)*y^3 + (3*x^4 + 5*x^2)*y^2 + (5*x^8 + 6*x^6)*y + 2*x^7 + 2*x^5 + 2*x^4 + 4*x^3 + x^2
h = haccent(x^m) %(x^n + 1) h 
        
6*x^43 + 3*x^39 + 4*x^30 + 2*x^28 + 4*x^27 + 3*x^20 + 5*x^18 + 5*x^16 +
6*x^14 + 2*x^7 + 2*x^5 + 2*x^4 + 4*x^3 + x^2
6*x^43 + 3*x^39 + 4*x^30 + 2*x^28 + 4*x^27 + 3*x^20 + 5*x^18 + 5*x^16 + 6*x^14 + 2*x^7 + 2*x^5 + 2*x^4 + 4*x^3 + x^2

Controle:
f * g 
        
6*x^43 + 3*x^39 + 4*x^30 + 2*x^28 + 4*x^27 + 3*x^20 + 5*x^18 + 5*x^16 +
6*x^14 + 2*x^7 + 2*x^5 + 2*x^4 + 4*x^3 + x^2
6*x^43 + 3*x^39 + 4*x^30 + 2*x^28 + 4*x^27 + 3*x^20 + 5*x^18 + 5*x^16 + 6*x^14 + 2*x^7 + 2*x^5 + 2*x^4 + 4*x^3 + x^2

Nu met DFT. We doen enkel de berekening van hss anders:
omega = eta^2 
       
DFTfss = DFT(omega, fss, t) DFTgss = DFT(omega, gss, t) 
       
DFThss = punt_product(DFTfss, DFTgss) 
       
hss = DFT(omega^-1, DFThss, t)/n hss 
        
6*xbar^13*y^5 + 3*xbar^15*y^4 + (4*xbar^12 + 2*xbar^10 + 4*xbar^9)*y^3 +
(3*xbar^8 + 5*xbar^6)*y^2 + (5*xbar^10 + 6*xbar^8)*y + 2*xbar^7 +
2*xbar^5 + 2*xbar^4 + 4*xbar^3 + xbar^2
6*xbar^13*y^5 + 3*xbar^15*y^4 + (4*xbar^12 + 2*xbar^10 + 4*xbar^9)*y^3 + (3*xbar^8 + 5*xbar^6)*y^2 + (5*xbar^10 + 6*xbar^8)*y + 2*xbar^7 + 2*xbar^5 + 2*xbar^4 + 4*xbar^3 + xbar^2

Dit is dezelfde hss als de vorige berekening.