MACM 401/MATH 数值分析

MACM 401/MATH 701/MATH 801
Assignment 3, Spring 2019.
Michael Monagan
Due Monday February 25th at 4pm. Hand in to Dropoff box 1b outside AQ 4100.
Late Penalty: ?20% for up to 48 hours late. Zero after that.
For problems involving Maple calculations and Maple programming, you should submit a printout
of a Maple worksheet of your Maple session.
MATH 801 students should do both questions 4 and 5.
MACM 401 students may do either question 4 or 5. For a 15 mark bonus, you may do both.
Question 1: The Fast Fourier Transform (30 marks)
(a) Let n = 2m and let ω be a primitive n’th root of unity. To apply the FFT recursively, we use
the fact that ω
2
is a primitive m’th root of unity. Prove this.
Also, for p = 97 = 3 × 2
5 + 1, find a primitve 8’th root of unity in Zp. Use the method in
Section 4.8 which first finds a primitive element 1 < α < p ? 1 of Zp.
(b) What is the Fourier Transform for the polynomial a(x) = 1 + x + x
2 + · · · + xn1, i.e., what
is the vector [a(1), a(ω), a(ω2), . . . , a(ωn1)]
(c) Let M(n) be the number of multiplications that the FFT does. A naive implementation of
the algorithm leads to this recurrence:
M(n) = 2M(n/2) + n + 1 for n > 1
with initial value M(1) = 0. In class we said that if we pre-compute the powers ωi for 0 ≤ i ≤ n/2 and store them in an array W, we can save half the multiplications in the FFT
so that
M(n) = 2M(n/2) + n2 for n > 1.
By hand, solve this recurrence and show that M(n) = 12n log2 n.
(d) Using Maple’s rsolve command, solve the following recurrences. Please simplify the output
from rsolve.
(i) T(1) = d, T(n) = 3T(n/2) + cn for n > 1 (Karatsuba),
(ii) M(1) = 0, M(n) = 2M(n/2) + n/2 for n > 1 (optimized FFT) and
(iii) T(1) = 0, T(n) = T(n 1) + (n 1)2
for n > 1 (Gaussian elimination).
(e) Program the FFT in Maple as a recursive procedure. Your Maple procedure should take as
input (n, A, p, ω) where n is a power of 2, A is an array of size n storing the input coefficients
a0, a1, . . . , an1, a prime p and ω a primitive n’th root of unity in Zp. If you want to precompute
an array W = [1, ω, ω2
, . . . , ωn/21
] of the powers of ω to save multiplications you may
do so.
Test your procedure on the following input. Let A = [1, 2, 3, 4, 3, 2, 1, 0], p = 97 and w
be the primitive 8’th root of unity. To check that your output B is correct, verify that
F F T(n, B, p, ω1
) = nA mod p.
1
(f) Let a(x) = x
3 + 3x + 1 and b(x) = 2x
2 + x + 1 be polynomials in Z97[x].
Calculate the product of c(x) = a(x)b(x) using the FFT.
If you could not get your FFT procedure from part (c) to work, use the following one which
computes [a(1), a(ω), . . . , a(ωn1)] using ordinary evaluation.
FourierTransform := proc(n,A,p,omega)
local f,x,i,C,wi;
f := add(A[i]*x^i, i=0..n-1);
C := Array(0..n-1);
wi := 1;
for i from 0 to n-1 do
C[i] := Eval(f,x=wi) mod p;
wi := wi*omega mod p;
od;
return C;
end:
Question 2: The Modular GCD Algorithm (15 marks)
Consider the following pairs of polynomials in Z[x].
a1 = 58 x
4 415 x
3 111 x + 213
b1 = 69 x
3 112 x
2 + 413 x + 113
a2 = x
5 111 x
4 + 112 x
3 + 8 x
2 888 x + 896
b2 = x
5 114 x
4 + 448 x
3 672 x
2 + 669 x 336
a3 = 396 x
5 36 x
4 + 3498 x
3 2532 x
2 + 2844 x 1870
b3 = 156 x
5 + 69 x
4 + 1371 x
3 332 x
2 + 593 x 697
Compute the GCD(ai
, bi) using the modular GCD algorithm. Use primes p = 23, 29, 31, 37, 43, ….
Identify which primes are bad primes and which are unlucky primes.
To compute gcd(φp(a), φp(b)) in Maple, use Gcd(a,b) mod p . Use the Maple commands chrem
for Chinese remaindering, mods to put the coefficients in the symmetric range, and any other Maple
commands that you need.
PLEASE make sure you input the polynomials correctly!
Question 3: Resultants (15 marks)
(a) Calculate the resultant of A = 3x
2 + 3 and B = (x 2)(x + 5) by hand.
Also, calculate the resultant using Maple. See resultant
(b) Let A, B, C be non-constant polynomials in R[x].
Show that res(A, BC) = res(A, B) · res(A, C).
2
(c) Let A, B be two non-zero polynomials in Z[x]. Let A = GAˉ and B = GBˉ where G = gcd(A, B).
Recall that a prime p in the modular gcd algorithm is unlucky iff p|R where R = res(A, ˉ Bˉ) ∈ Z.
Consider the following pair of polynomials from question 2.
Aˉ = 58 x
4 415 x
3 111 x + 213
Bˉ = 69 x
3 112 x
2 + 413 x + 113
Using Maple, compute the resultant R and identify all unlucky primes. For each unlucky
prime p compute Gcd(A, ˉ Bˉ) mod p in Maple to verify that the primes are indeed unlucky.
Question 4: The Chinese remainder theorem in F[y] (15 marks).
Consider the problem of computing GCDs in Zqy, q a prime. If q is large then we can use
evaluation and interpolation for y, i.e., we can evaluate at y = 0, 1, 2, … and interpolate the coeffi-
cients of the GCD in Zq[y]. If q is small, e.g. q = 2, this will not work as there will be insufficient
evaluation points in Zq. Moreover, y = 0 and y = 1 may be bad or unlucky.
But Zq[y] is a Euclidean domain and there are an infinite number of primes (irreducibles) in
Zq[y] which can play the role of integer primes in the modular GCD algorithm for computing GCDs
in Z[x]. For example, here are the irreducibles in Z2[y] up to degree 4.
y, y + 1, y2 + y + 1, y3 + y + 1, y3 + y
2 + 1, y4 + y + 1, y4 + y
3 + 1, y4 + y
3 + y
2 + y + 1.
To do this we need to solve the Chinese remainder problem in Zq[y].
Theorem: Let F be any field (e.g. Zq) and let m1, m2, . . . , mn and u1, u2, . . . , un be polynomials
in F[y] with deg(mi) > 0 for 1 ≤ i ≤ n. If gcd(mi
, mj ) = 1 for 1 ≤ i < j ≤ n then there
exists a unique polynomial u in F[y] s.t.
(i) u ≡ ui (mod mi) for 1 ≤ i ≤ n and
(ii) u = 0 or deg u < Pn
i=1 deg mi
.
Prove the theorem by modifying the proof of the Chinese remainder theorem for Z to work for
F[y]. Now solve the following Chinese remainder problem: find u ∈ Z2[y] such that
u ≡ y
2
(mod y
3 + y + 1) and u ≡ y
2 + y + 1 (mod y
3 + y
2 + 1).
Note, in the statement of the theorem the congruence relation u ≡ ui (mod mi) means mi
|(u ui)
in F[y]. For the extended Euclidean algorithm in Zq[y], use Maple’s Gcdex(…) mod q command
to compute the required inverse.
Question 5: Multivariate Polynomial Division (15 marks)
In assignment 2 question 2 we were given two polynomials A, B ∈ Z[x1, x2, . . . , xn] with B 6= 0, and
I asked you to give pseudo code for the multivariate division algorithm for dividing A by B. I said
the pseudo code should begin like this
Algorithm DIVIDE(A,B)
Inputs A, B ∈ Z[x1, x2, . . . , xn] satisfying B 6= 0 and n ≥ 0.
Output Q ∈ Z[x1, x2, . . . , xn] s.t. A = BQ or FAIL meaning B does not divide A.
3
For this assignment implement your multivariate division algorithm in Maple as the Maple
procedure DIVIDE(A,B) to divide A by B. Test you program on the following inputs

A := (6y^2-5yz+z^2)x^2+(7y^2z-3yz^2)x+2y^2*z^2;
B := (2y-z)x+y*z;
Q := DIVIDE(A,B);
Q := DIVIDE(A+x,B);
Q := DIVIDE(A+2,B);
C := expand(A*B);
Q := DIVIDE(C,B);
The following operations will be helpful.
X := indets(A) union indets(B); # set of all variables
X := {x, y, z}
var := X[1];
var := x
degree(B,var);
2
lcoeff(B,var); # leading coefficient
2y z
I suggest that you get your procedure working with zero variables first, then one variable, then
two variables then three variables. If your procedure is not working you may insert a print(…);
command anywhere in your procedure to print out any value. Also, you may trace the execution
of your procedure by using
trace(DIVIDE);
Maple will display everything that is computed. To turn tracing off use
untrace(DIVIDE);
WX：codehelp

[温馨提示：高防服务器能助您降低 IT 成本，提升运维效率，使您更专注于核心业务创新。]

[图文来源于网络，不代表本站立场，如有侵权，请联系高防服务器网删除]