Evaluation/Interpolation (II)Backtrack from the GCD ideas a bitSlide 3How does this work in some more detailSlide 5Chinese Remainder TheoremChinese Remainder ExampleConversion to modular CRA formConversion to Normal Form (Garner’s Alg.)Slide 10Slide 11Slide 12Interpolationpolynomial interpolation (Lagrange)polynomial interpolationAnother, incremental, approach (Newton Interpolation, CRA like)An incremental approach (Newton Interpolation)Slide 18Sparse Multivariate InterpolationSparse Multivariate Interpolation: Basic Idea / exampleFirst stage: Find a “skeleton” in variable xWhat if the first stage is wrong?Second stageSecond stage assumption of sparsenessThird stageHow good is this?Polynomial evaluationSingle polynomial evaluationTwo-point polynomial evaluation, at r, -rShort break... Factoring a polynomial h(x) using interpolation and factoring integersTwo possible directions to go from hereRichard Fateman CS 282 Lecture 8 1Evaluation/Interpolation (II)Lecture 8Richard Fateman CS 282 Lecture 8 2Backtrack from the GCD ideas a bit•We can do some operations faster in a simpler domain, even if we need to do them repeatedlyGCD(F(x,y,z),G(x,y,z))GCD(Fp(x,0,0),Gp(x,0,0))GCD(Fp(x,y,z),Gp(x,y,z))H(x,y,z)Hp(x,y,z)Hp(x,0,0)Richard Fateman CS 282 Lecture 8 3Backtrack from the GCD ideas a bit•Some of the detailsGCD(F(x,y,z),G(x,y,z))GCD(Fp(x,0,0),Gp(x,0,0))GCD(Fp(x,y,z),Gp(x,y,z))H(x,y,z)Hp(x,y,z)Hp(x,0,0)reduce mod p1, p2, ...evaluate at z=0,1,2,... evaluate at y=0,1,2,... compute (many) GCDs interpolate OR sparse interpolateCRA p1, p2, ...Richard Fateman CS 282 Lecture 8 4How does this work in some more detail•How many primes p1, p2, ...?–Bound the coeffs by p1*p2 *...pn ? (or +/- half that)•bad idea, the bounds are poor. given ||f|| what is ||h|| where h is factor of f ? What is a bad example? A pyramid X a difference operation? Are there worse?•(x-1)*(1+2x+3x2+....+nxn-1+ (n-1)xn+..+2x2n-1+x2n) is•-1-x-x2-x3...-xn+xn+1 + ... +x2n+1 ... •Some bound.... ||h|| · (d+1)1/22d max(||f||) •Don’t bound but try an answer and test to see if it divides?•See if WHAT divides? –compute cofactors A, B, A*H=F, B*H=G, and when A*H= F or B*H=G, you are done.Richard Fateman CS 282 Lecture 8 5How does this work in some more detail–The process doesn’t recover the leading coefficient since F modulo p etc might as well be monic.–The inputs F and G are assumed primitive; restore the contents.–There may be unlucky primes or evaluation points.Richard Fateman CS 282 Lecture 8 6Chinese Remainder Theorem•(Integer) Chinese Remainder Theorem:•We can represent a number x by its remainders modulo some collection•of relatively prime integers n1 n2 n3...•Let N = n0*n1* ...*nk. Then the Chinese Remainder Thm. tells us that we can represent any number x in the range -(N-1)/2 to +(N-1)/2 by its residues modulo n0, n1, n2, ..., nk. {or some other similar sized range, 0 to N-1 would do}Richard Fateman CS 282 Lecture 8 7Chinese Remainder ExampleExample n0=3, n2=5 , N=3*5=15x x mod 3 x mod 5-7 -1 -2 note: if you hate balanced notation -7+15=8. mod 3 is 2->-1-6 0 -1-5 1 0-4 -1 1-3 0 2-2 1 -2 note: x mod 5 agrees with x, for small x 2 [-2,2], +-(n-1)/2-1 -1 -1 note:0 0 0 note:1 1 1 note:2 -1 23 0 -24 1 -15 -1 0 note: symmetry with -56 0 17 1 2 note: also 22, 37, .... and -8, ...Richard Fateman CS 282 Lecture 8 8Conversion to modular CRA formConverting from normal notation to modular representation is easy in principle;you do remainder computations (one instruction if x is small, otherwise a software routine simulating long division)Richard Fateman CS 282 Lecture 8 9Conversion to Normal Form (Garner’s Alg.)Converting to normal rep. takes k2 steps. Beforehand, computeinverse of n1 mod n0, inverse of n2 mod n0*n1, and also the products n0*n1, etc. Aside: how to compute these inverses:These can be done by using the Extended Euclidean Algorithm.Given r= n0, s= n1, or any 2 relatively prime numbers, EEA computes a, b such that a*r+b*s=1 = gcd(r,s)Look at this equation mod s: b*s is 0 (s is 0 mod s) and so we have a solution a*r=1 and hence a = inverse of r mod s. That’s what we want. It’s not too expensive since in practice we can precompute all we need, and computing these is modest anyway. (Proof of this has occupied quite a few people..)Richard Fateman CS 282 Lecture 8 10Conversion to Normal Form (Garner’s Alg.)Here is Garner's algorithm:Input: x as a list of residues {ui}: u0 = x mod n0, u1 = x mod n1, ...Output: x as an integer in [-(N-1)/2,(N-1)/2]. (Other possibilities include x in another range, also x as a rational fraction!)Consider the mixed radix representation x = v0 + v1*n0 + v2*(n0*n1)+ ... +vk*(n0*...nk-1) [G]if we find the vi, we are done, after k more mults. These products are small X bignum, so the cost is more like k2/2.Richard Fateman CS 282 Lecture 8 11Conversion to Normal Form (Garner’s Alg.)x = v0 + v1*n0 + v2*(n0*n1)+ ... +vk*(n0*...nk-1) [G]It should be clear that v0 = u0, each being x mod n0Next, computing remainder mod n1 of each side of [G]u1 = v0 +v1*n0 mod n1, with everything else dropping outso v1 = (u1-v0)* n0-1 all mod n1in general, vk = ( uk - [v0+v1*n0 +...+vk-1* (n0...nk-2) ]) * (n0*...nk-1)-1 mod nk. mod nkRichard Fateman CS 282 Lecture 8 12Conversion to Normal Form (Garner’s Alg.)Note that all the vk are “small” numbers and the items in green are precomputed.Cost: if we find all these values in sequence, we have k2 multiplication operations, where k is the number of primes needed. In practice we pick the k largest single-precision primes, and so k is about 2* log (SomeBound/231)Richard Fateman CS 282 Lecture 8 13Interpolation•Abstractly THE SAME AS CRA–change primes p1, p2, ... to (x-x1) ...–change residues u1= x mod p1 for some integer x to yk=F(xk) for a polynomial F in one variable–Two ways it is usually presented, sometimes more (solving linear vandermonde system...)Richard Fateman CS 282 Lecture 8 14polynomial interpolation (Lagrange)•The problem: find the (unique) polynomial f(x) of degree k-1 given a set of evaluation points {xi}i=1,k and a set of values {yi=f(xi)}Solution: for each i=1,...,kfind a polynomial pi(x) that takes on the value yi at xi, and is zero for all other abscissae..x1, ...,xi-1,..xi+1,..xkThat’s easy here’s one solution: pi(x)=yi +(x-x1)(x-x2)¢ ...¢(x-xi-1)(x-xi+1)...Richard Fateman CS 282 Lecture 8