Transcript Document

Basic Domains of Interest used in
Computer Algebra Systems
Lecture 2b
Richard Fateman CS 282 Lecture 2b
1
Some ideas just for representing integers
Integers are sequences of characters, 0..9.
Integers are sequences of words modulo 109 which is the
largest power of 10 less than 231.
Integers are sequences of hexadecimal digits.
Integers are sequences of 32-bit words.
Integers are sequences of 64-bit double-floats (with 8
bits wasted).
Sequences are linked lists
Sequences are vectors
Sequences are stored in sequential disk locations
Richard Fateman CS 282 Lecture 2
2
Aside: using half the bits in a word is
common strategy
Integers are sequences of 32-bit words, but only
the bottom 16 are used: There is some lack of
uniformity in architectural support for unsigned
32X32 bit multiply. But 16X16  32 bit is
supported, so programs can be portable.
Downside: if you multiply two n-bigit numbers in
time n2 then with half-length bigits you need
twice as many of them and so you take time
(2n)2, or 4 times slower and 2X the space.
Richard Fateman CS 282 Lecture 2
3
Yet more ideas
Integers are stored in redundant form a+b+…
Integer are stored modulo set of different
primes
Integers are stored in p-adic form as a sequence
of x mod p, x mod p2, …
Richard Fateman CS 282 Lecture 2
4
Redundant (big precision) floats
• Sequences (xn + ….+x2 +x1)
• Non-overlapping means each the lsb of xk is more
significiant than the msb of xk-1. May be big gaps.
• Not unique.
• (binary values..) 1100 + -10.1 = 1001 + 0.1 = 1000+1+0.1
• Easy to tell the sign. Look at the leading term.
• Adding must restore non-overlapping property.
• Important use by Jonathan Shewchuk (UCB) in
geometric predicate calculations.
Richard Fateman CS 282 Lecture 2
5
Modular (mod a set of primes q1, q2, …qn)
• images unique only within multiple of product of
the primes Q = q1 ¢ q2 ¢ …¢ qn.
• CRA (Chinese Remainder Algorithm) provides a
way of going from modular to conventional
positional notation, but takes O(n2) in practice.
• This and its generalizations heavily used in
computer algebra systems and this course.
Richard Fateman CS 282 Lecture 2
6
Modular arithmetic is really fast…
all the arithmetic can be done without carry, in
parallel.
Not usually used because
(a) You can’t tell for sure if a number is +, -, 0 or
Q or 3 Q….
(b) Parallelism is almost always irrelevant 
(c) If you must see the answer converted to
decimal, the conversion is O(n2)
(d) Conversion to decimal may be very common if
your application is a bignum calculator.
Richard Fateman CS 282 Lecture 2
7
What’s a p-adic integer?
For the moment assume p is a prime number.
Consider representing an integer a =
a0+a1.p+a2¢ p2 + … where ai are chosen from
integers in the range 0 · ai < p
For any finite positive integer a, either all the ak
are zero in which case a = 0 and ||a||p is 0,
or some initial set of the ai are zero. Let ar+1
be the first non-zero term. ||a||p = p-r. This
replaces the absolute-value valuation |a|
where we consider that x and y are close if
x=y mod pk for many values of k=0,1,2,….
Richard Fateman CS 282 Lecture 2
8
p-adic ordering is odd.
14 3-adically is 2*30+1*31+1*32= 2 + 3 + 9 = 14
Which is very close to 5, 3-adically because
5 is 2*30+1*31, and they are the same modulo 30 and 31.
||15-4||3 = 3-1
What about negative numbers? pk-1 is possible, but
consider
pk-1 = (p-1)¢(pk-1+pk-2…+p+1)=
(p-1)*1+(p-1)*p+(p-1)*p2+… so
-1 3-adically is 2*30+2*31+…. Infinite number of terms
(2,2,2,2,2,….)
Richard Fateman CS 282 Lecture 2
9
How to multiply p-adically
(a,b,c)
(d,e,f) 
Compute a*d , a*e, a*f;
add columns
b*d, b*e; b*f
c*d, c*e, c*f
Add modulo p, p2, p3 … with a carry.
Richard Fateman CS 282 Lecture 2
10
How to multiply p-adically approximately
(a,b,c …)
(d,e,f …) 
Compute a*d , a*e, a*f;
add columns
b*d, b*e; b*f
c*d, c*e, c*f
Ignore
these
Add modulo p, p2, p3 … with a carry. Dropping off
extra terms is like 3.14159 vs 3.14, but with
respect to p-adic distances.
Richard Fateman CS 282 Lecture 2
11
What other p-adic numbers are there??
-1/2 3-adically .. is 1*30+1*31+1*32= (1,1,1,1,1…). Proof:
Multiply by 2, which is (2,0,0,0,….) to get (2,2,2,2…) which
is -1
What about p 7 ?
(a0+3*a1+…)2-7=0 (mod 3i+1)
Mod 30, a02-7=0 so a0 is 1 or –1. Let’s choose 1.
Next solve
(1+3*a1+…)2-7=0 (mod 31)
= 1+2*3*a1 –7 =0 so a1=1
Eventually, p7 = (1,1,1,2,…)
Richard Fateman CS 282 Lecture 2
12
What does this buy us??
Not a great deal for integers, but…
We’ll use p-adic representation where p is not a
number, but an indeterminate (say x), or a
polynomial (x^2-3). Or a polynomial in several
variables (x+y+1). If we compute a gcd of 2
polynomials p-adically approximately to high
enough degree, we will know the exact GCD. If
we can do this computation faster than other
means, we have a winner. (This is the case.)
Richard Fateman CS 282 Lecture 2
13
Multiplication, usual representation
Extremely well studied.
The usual method takes O(n2),
Karatsuba style O(n1.585)
or FFT style O(n log n).
These will be studied in the context of multiplying
polynomials.
Note that 345 can be mapped to p(x)=3x2+4x+5
where p(10) is 345.
Except for the “carry”, the operation is the same.
Richard Fateman CS 282 Lecture 2
14
Integer Division
• This is too tedious to present in a lecture.
• Techniques for guessing the next big digit
(bigit) of a quotient within +1 are available
• For exact division, consider Newton iteration
is an alternative
• FFT / fast multiplication helps
Richard Fateman CS 282 Lecture 2
15
GCD
• Euclid’s algorithm is O(n2 log n) but is hard to
beat in practice, though see analysis of HGCD
(Yap) for an O(n log2 n) algorithm..
• HGCD is portrayed as a winner for
polynomials, but only by complexity analysts
who (especially in this case) assume that
–
certain costs are constant when in fact they grow
exponentially.
– Multiplication cost n log n when, for relevant cases
n2 algorithms are much faster than the n log n ones
Richard Fateman CS 282 Lecture 2
16
Reminder… A Ring R is Euclidean
If there is a function y
Richard Fateman CS 282 Lecture 2
17
Shows the tendency to obfuscate…
What Rings do we use, and what is y?
For integers, absolute value | | = y
For p-adic numbers, || ||p norm
For polynomials in x, degree in x = y
Richard Fateman CS 282 Lecture 2
18
Where next?
• We could spend a semester on integer
arithmetic, but this does not accomplish any
higher goals of CAS
• We can proceed to build models of real
numbers by approximation (e.g. as limit of
rational intervals). We may return to this..
• We proceed to polynomials, typically with
integer coefficients or finite field coeffs.
Richard Fateman CS 282 Lecture 2
19
Interlude, regarding your homework (the last
problem)
• Usual representation in Lisp is trivial.
• 3*x^2+4 is something like (+ (* 3 (^ x 2) 4)
• In Lisp, * is a symbol. Times is a symbol..
Richard Fateman CS 282 Lecture 2
20
Trivial parts..
• Given a, b, a program to produce a product is
(defun prod(a b)(list ‘* a b)) ;common lisp syntax
(define (prod a b)(list ‘* a b)) ;scheme syntax
function make_prod(a:tree,b:tree):tree
var temp:^tree; {hm, something like this..}
begin
temp:= new(tree); temp^.head:= asterisk; temp^.left:=a;
temp^.right:=b, make_prod:=temp
end {Pascal?? Similarly for C, Java, …}
Richard Fateman CS 282 Lecture 2
21
Harder parts..
• If you need to know that (* x 0), (+ x (* -1 x))
and 0 are the same, how much work must you
do?
• (Schwartz-Zippel polynomial identity testing)
Richard Fateman CS 282 Lecture 2
22
Reminder: The Usual Coefficients
•
•
•
•
•
Z: natural numbers (ring, integers, +-* not closed under
division!)
Zp :integers modulo p, p a prime usually.
Q: rational numbers (a field). What is the difference
between 3/4 and 6/8?
R: real includes irrationals p 2 ,transcendentals (e, p)
{mention continued fracs, intervals}
C: complex
– approximation
– special case of algebraic extension
– Analysis: functions of a complex variable
Richard Fateman CS 282 Lecture 2
23
Preview: Extensions
if x is an indeterminate, and
if D is a domain, we can talk about D[x],
polynomials in x with coefficients in D
or D(x) ratios of polys in x.
Or D[[x]] : (truncated) power series in x over D.
Or Matrices over D.
Richard Fateman CS 282 Lecture 2
24
Beware of notation
A particular polynomial expression might be
referred to as p(x), which notation is also used
to denote a function or mapping p:D1D2 or a
function application if x is not an indeterminate
but a element in a field as p(3). Confused? The
notation is unfortunate. Not confused? You are
probably used to it.
Richard Fateman CS 282 Lecture 2
25
Aside 2.1: canonical forms vs. mathematical
equivalence
Mathematically, we don’t distinguish between
two equivalent elements in Z(x), say 1/(x+1) and
(x-1)/(x2-1).
Computationally these can be distinguished, and
generally they must be distinguished. Often we
must compute a canonical form for an
expression by finding a particular “simplest”
form in an equivalence class.
Richard Fateman CS 282 Lecture 2
26
Aside 2.2: Simplification is almost
everything in this business..
Trivial reduction. All computational problems in
computer algebra can be reduced to
simplification:
simplify (ProblemStatement) to CanonicalSolution
Richard Fateman CS 282 Lecture 2
27
Aside 2.3: Computer representation and
canonical forms
A computer might distinguish between two strings
“abc” and “abc” if they are stored in different
locations in memory. Or might not.
Usually it is advantageous to store an object only
once in memory, but not always. (Should we
store 43 just once? How about 3.141592654 ?
How about ax2+bx+c?)
Richard Fateman CS 282 Lecture 2
28
Back to extensions
if r is a root of an irreducible polynomial p, that is,
p(r)=0, we will also talk about a ring or field extended
by r: Q[r]. E.g. p(r)=r2-1=0 means r = p(-1) or i, and we
have just constructed the complex rationals Q[r].
Z[i] is called “Gaussian integers" The set of elements
a+bi, with a, b, integers.
Q[i] would allow rational a, b. (remember rationalizing
denominators?)
Given such a field, you can extend it again.
IF you want to represent Q extended by sqrt(2),
and then THAT extended by sqrt(3), you can do so. Don't
extend it again by sqrt(6). (why?)
Richard Fateman CS 282 Lecture 2
29
More on extensions
In nice cases (primitive element), algebraic
arithmetic can be done by "reducing" modulo r. This is
accomplished by dividing by p(r) and discarding the
remainder: if E = a+b*p(r) then E´ a
You may need reminders of shortcuts. e.g. remainder of
p(x) / (x-a) is the same as substituting a for x in p.
(other terms we will use on occasion: Euclidean domains,
unique factorization domains, ideals, differential fields,
algebraic curves. We’ll motivate them when needed)
Richard Fateman CS 282 Lecture 2
30
Other extensions
Differential fields have what amounts to log() and
exp() extensions.
And an operation of differentiation such that
D(exp(x)) = exp(x), D(log(x)) =1/x.
Exp and log can be nested, and you can make trig
functions:
Richard Fateman CS 282 Lecture 2
31
There’s more… Other kinds of symbolic
computation
Whole careers have been made out of other kinds of
symbolic computation:
theorem proving, string manipulation, group
representations, geometric computation, type
theory/programming language representations, etc.
(J. Symbolic Computation publishes broadly…)
We will not probably not get to any of these areas in this
course, although I could be swayed by student
interest… also projects involving these topics are
generally appropriate.
Richard Fateman CS 282 Lecture 2
32