How fast can we multiply?
Christian Tismer
Lately, an old issue about long numbers and bitcounts
showed up again. I argued something about "if at all, then
build an extension with fast multiply".
I want to retract this. Python's C code needs a couple of
helper functions. But fast multiplication itself can be
easily done with plain vanilla Python.
Please try the attached script and enjoy.
The Python multiplication is twice as fast as the internal,
at about 100000 bits.
After this break, back to continuous stacklessness :-) chris
______________________________________________________________________
"""
fastlong.py
Fast long integers.
Christian Tismer
990718
From time to time I have argued on the main list that Python's
long integers should have a faster multiplication. It is a well
known fact that by some factorization, integer multiplication's
complexity can be brought from O(n**2) down to O(n**1.5) .
The method is sometimes called Karatsuba algorithm. But studying
chapter 12.7 of (1), it turns out to be nothing else but a special
case of preconditioning polynomials.
I think to remember that the GMP library uses this algorithm,
but I didn't look it up and tried my own straightforward way.
The reason why I don't aim to modify the C implementation any
longer is: It will not be reasonably faster in C. The real
working horse is still the basic internal long multiplication
(sometimes called "highschool algorithm") which is fast for small
long integers. The only effect of a C port would be to move the
break_even point to something smaller.
My proposal to get something like this into the C code:
For very long numbers to multiply, the C version could
be enabled to call back into the Python implementation.
A function to obtain the bit length of a long would
be most helpful. Some kind of split function as well.
Currently, with a "clean" nbits implementation, the
break even point appears to be 50000 bits.
With a marshal-based inaccurate hack, break even is
around 10000 bits.
For a bit length of 100000, my algorithm outperforms
the built in long numbers by a factor of 2.
(1) Aho/Hopcroft/Ullmann:
"The Design and Analysis of Computer Algorithms"
"""
class fastlong:
"""a class around long numbers which can multiply faster than
long"""
break_even = 50000 # clean nbits
break_even = 10000 # hacked nbits
def __init__(self, val):
if isinstance(val, fastlong):
self.val = val.val
else:
self.val = long(val)
def __repr__(self):
return repr(self.val)
def __add__(self, other):
return self.val + fastlong(other).val
def __sub__(self, other):
return self.val - fastlong(other).val
def __div__(self, other):
# not optimized yet
return self.val / fastlong(other).val
def __pos__(self):
return self
def __neg__(self):
return fastlong(-self.val)
def __lshift__(self, bits):
return self.val << bits
def __rshift__(self, bits):
return self.val >> bits
def __mul__(self, other):
"""Karatsuba alogrithm
(a+bX)(c+dX)=(X^2+X)(bd)+(X)(a-b)(d-c)+(X+1)(ac)
"""
val1 = self.val
val2 = fastlong(other).val
if not val1 and val2:
return 0
bitcount = max(nbits(val1), nbits(val2))
if bitcount < self.break_even:
return val1 * val2
shift = bitcount / 2
(hi1, lo1) = split_num(val1, shift)
(hi2, lo2) = split_num(val2, shift)
sum1 = lo1-hi1
sum2 = hi2-lo2
mul1 = hi1 * hi2
mul2 = lo1 * lo2
mul3 = sum1 * sum2
res = mul2 + ((mul1 + mul2 + mul3) << shift) + (mul1 <<
shift+shift)
return res
def split_num(longval, shift):
scale = 1L << (shift)
mask = scale - 1
lo = longval & mask
hi = longval >> shift
return fastlong(hi), fastlong(lo)
def nbits(x):
"""clean, accurate, not so fast"""
if x <= 2:
if x < 0:
return nbits(-x)+1
return x
sum = 1
while x:
shift = 1
x = x2 = x >> 1
while x2:
x = x2
sum = sum + shift
shift = shift + shift
x2 = x >> shift
return sum
import marshal, struct
def nbits(x):
"""not clean, rather fast, inaccurate"""
s=marshal.dumps(x)
return abs(struct.unpack("i", s[1:5])[0])*15 -7
#---------------------------------------------
def test():
import sys
be = fastlong.break_even
for i in range(be/2, be*16, be/4):
arg = (1l << i+1) -1
fastarg = fastlong(arg)
tim1, res1 = timing(lambda x:x*x, arg)
tim2, res2 = timing(lambda x:x*x, fastarg)
if res1 <> res2:
raise ValueError, "different results fo i=" % i
print i, tim1, tim2
sys.stdout.flush()
def timing(func, args, n=1, **keywords) :
import time
time=time.time
appl=apply
if type(args) != type(()) : args=(args,)
rep=range(n)
before=time()
for i in rep: res=appl(func, args, keywords)
return round(time()-before,4), res
# the end ----------------------------------------------------
______________________________________________________________________
--
Christian Tismer :^) <mailto:tis...@appliedbiometrics.com>
Applied Biometrics GmbH : Have a break! Take a ride on Python's
Kaiserin-Augusta-Allee 101 : *Starship* http://starship.python.net
10553 Berlin : PGP key -> http://wwwkeys.pgp.net
PGP Fingerprint E182 71C7 1A9D 66E9 9D15 D3CC D4D7 93E2 1FAE F6DF
we're tired of banana software - shipped green, ripens at home
Christian Tismer
I detected an oversight, which made me loose a lot of speed.
Version 0.2 is now 4 times faster at 100000 bits.
"""
fastlong.py
Fast long integers.
Christian Tismer
990718
Version 0.2:
Adding some forgotten fastlong calls made it much faster.
the built in long numbers by a factor of 4.
(1) Aho/Hopcroft/Ullmann:
"The Design and Analysis of Computer Algorithms"
"""
class fastlong:
"""a class around long numbers which can multiply faster than
long"""
break_even = 50000 # clean nbits
break_even = 10000 # hacked nbits
def __init__(self, val):
if isinstance(val, fastlong):
self.val = val.val
else:
self.val = long(val)
def __repr__(self):
return repr(self.val)
def __add__(self, other):
return fastlong(self.val + fastlong(other).val)
def __sub__(self, other):
return fastlong(self.val - fastlong(other).val)
def __div__(self, other):
# not optimized yet
return fastlong(self.val / fastlong(other).val)
def __pos__(self):
return self
def __neg__(self):
return fastlong(-self.val)
def __lshift__(self, bits):
return fastlong(self.val << bits)
def __rshift__(self, bits):
return fastlong(self.val >> bits)
import marshal, struct
#---------------------------------------------
--
Tim Peters
> ...
> Version 0.2 is now 4 times faster at 100000 bits.
Factor of 5.7 on my machine.
> ...
> fastlong.py
> ...
> It is a well known fact that by some factorization, integer
> multiplication's complexity can be brought from O(n**2) down to
> O(n**1.5).
The exponent is actually log2(3) ~= 1.58.
> The method is sometimes called Karatsuba algorithm.
> ...
> I think to remember that the GMP library uses this algorithm,
> but I didn't look it up and tried my own straightforward way.
Right, GMP does do this. There's a much faster method (O(n)) based on
Fourier transforms, but it's much harder to code and I've only seen it used
in David Bailey's MPFUN package (a highly optimized Fortran pkg for
arbitrary-but-fixed precision float arithmetic + the std elementary
functions).
> The reason why I don't aim to modify the C implementation any
> longer is: It will not be reasonably faster in C. The real
> working horse is still the basic internal long multiplication
> (sometimes called "highschool algorithm") which is fast for small
> long integers. The only effect of a C port would be to move the
> break_even point to something smaller.
Tee hee. I wonder how much of Python we could speed up by throwing out C
code <wink>? Integer exponentiation can also be done faster in Python, and
much faster with a faster multiply to build on. Seem to recall that
repr(dict) is/was quadratic-time in len(dict), while easy to make
linear-time if coded in Python. C is so clumsy you just settle for the
simplest thing to *spell*.
> ...
> A function to obtain the bit length of a long would
> be most helpful.
The lack of this has gotten unbearable, eh?
> Some kind of split function as well.
That one is much weaker; unlike bit length, the obvious ways to split and
join already have the right O() behavior.
> ...
> (1) Aho/Hopcroft/Ullmann:
> "The Design and Analysis of Computer Algorithms"
Knuth Vol 2 is a better reference for fast integer mult, although it goes
off the deep end on theoretical refinements of no practical value ("what if
the number of bits is bigger than the number of electrons in the
universe?").
> ...
> def __mul__(self, other):
> """Karatsuba alogrithm
> (a+bX)(c+dX)=(X^2+X)(bd)+(X)(a-b)(d-c)+(X+1)(ac)
> """
> val1 = self.val
> val2 = fastlong(other).val
> if not val1 and val2:
> return 0
The test only succeeds when val1 == 0 and val2 != 0 -- not what you
intended. OTOH, it's debatable whether programs multiply by 0 often enough
to pay back the cost of testing for it!
> bitcount = max(nbits(val1), nbits(val2))
min is a better choice here; see below.
> if bitcount < self.break_even:
> return val1 * val2
> shift = bitcount / 2
> (hi1, lo1) = split_num(val1, shift)
> (hi2, lo2) = split_num(val2, shift)
Need a fancier strategy when nbits(val1) and nbits(val2) are far apart.
This strategy takes strong advantage of the special-case-for-0 test above
<wink>, but can easily run 100s of times slower than the builtin mult; e.g.,
change your test case to lambda x:x*5.
This is all the more reason to write in Python, though! Much clumsier to
handle imbalance in C. In the meantime, using "min" instead of "max" to
compute bitcount greatly speeds the worst cases, and doesn't hurt the best
ones.
> ...
> def nbits(x):
> """not clean, rather fast, inaccurate"""
> s=marshal.dumps(x)
> return abs(struct.unpack("i", s[1:5])[0])*15 -7
I've been using the following since 1.5.2 came out. Cleaner and accurate.
The one above is substantially faster, but it's now a small-factor thing
rather than an O() difference:
def nbits(n, hex=hex, long=long, len=len,
d2bits={'0':0,
'1':1,
'2':2, '3':2,
'4':3, '5':3, '6':3, '7':3,
'8':4, '9':4, 'A':4, 'B':4,
'C':4, 'D':4, 'E':4, 'F':4}):
"""nbits(n): return number of bits in n.
n must be >= 0; 0->0, 1->1, 2->2, 3->2, 4->3, 5->3, ...
In 1.5.1 and before this takes time quadratic in the number
of bits. In 1.5.2, linear. "Should be" constant time, but
that requires getting at the internal representation.
"""
if n < 0:
raise ValueError("nbits arg must be >= 0")
# hex representation is ashex[2:-1] (leading 0x and trailing L).
ashex = hex(long(n))
num_hex_digits = len(ashex) - 3
return (num_hex_digits - 1)*4 + d2bits[ashex[2]]
> ...
> def timing(func, args, n=1, **keywords) :
> import time
> time=time.time
Not time.clock?
> ...
> # the end ----------------------------------------------------
looks-more-like-the-beginning-to-me<wink>-ly y'rs - tim
Christian Tismer
Tim Peters wrote:
>
> [Christian Tismer, hiding from continuations <wink>]
Oh, just a break.
> > Version 0.2 is now 4 times faster at 100000 bits.
>
> Factor of 5.7 on my machine.
And the cutoff should better be set to 5000, again better.
> > It is a well known fact that by some factorization, integer
> > multiplication's complexity can be brought from O(n**2) down to
> > O(n**1.5).
>
> The exponent is actually log2(3) ~= 1.58.
I dumbhead thought of 4 mults -> 3 and just said 1.5 - but yes,
sure, sorry.
[Karatsuba, GMP,...]
> Right, GMP does do this. There's a much faster method (O(n)) based on
> Fourier transforms, but it's much harder to code and I've only seen it used
> in David Bailey's MPFUN package (a highly optimized Fortran pkg for
> arbitrary-but-fixed precision float arithmetic + the std elementary
> functions).
Can you give me a pointer? What I know is O(n log n loglog n),
based on convolutions, but I may be 20 years behind state of the art.
[we agree that coding in Python leads to etter solutions]
> > ...
> > A function to obtain the bit length of a long would
> > be most helpful.
>
> The lack of this has gotten unbearable, eh?
This is my old fault. If there is a good algorithm with a better Oh,
that's one thing. But then it hurts so much to see an unnecessary
constant factor left over. The counter (number of words) is even
part of the long number and cries to be used.
> > Some kind of split function as well.
>
> That one is much weaker; unlike bit length, the obvious ways to split and
> join already have the right O() behavior.
Sure. I should stop thinking about it. Just so tempting
to have a functon to break a long into two peaces, at
a word boundary, no shifts... sigh. No, I resist [now].
>
> > ...
> > (1) Aho/Hopcroft/Ullmann:
> > "The Design and Analysis of Computer Algorithms"
>
> Knuth Vol 2 is a better reference for fast integer mult, although it goes
> off the deep end on theoretical refinements of no practical value ("what if
> the number of bits is bigger than the number of electrons in the
> universe?").
I know, also the subject "how..." is from there. But they re-arranged
my office. I found the Knuth corner again five minutes ago.
Yes, there is also the O(n) thing, and I'm not so behind.
> > if not val1 and val2:
> The test only succeeds when val1 == 0 and val2 != 0 -- not what you
> intended. OTOH, it's debatable whether programs multiply by 0 often enough
> to pay back the cost of testing for it!
How could I forget the braces. Well, the "why" was in fact to save
time when my "max" is a bad guess.
...
> Need a fancier strategy when nbits(val1) and nbits(val2) are far apart.
> This strategy takes strong advantage of the special-case-for-0 test above
> <wink>, but can easily run 100s of times slower than the builtin mult; e.g.,
> change your test case to lambda x:x*5.
I was undecided what's best. Obviously you are right. The asymmetric
case has a small and a large side, it already doesn't suffer from
the quadratic "flesh". I see it.
Also I was tempted to use "&" on the numbers, in order to save
one nbits call. It was just unclear what happens with negative
numbers and dropped it. Was happy that negatives went correct without
special treatment.
> This is all the more reason to write in Python, though! Much clumsier to
> handle imbalance in C. In the meantime, using "min" instead of "max" to
> compute bitcount greatly speeds the worst cases, and doesn't hurt the best
> ones.
> > def nbits(x):
> > """not clean, rather fast, inaccurate"""
> > s=marshal.dumps(x)
> > return abs(struct.unpack("i", s[1:5])[0])*15 -7
Yes it is ugly. I just needed to lower the break-even point,
since I was impatient for the result. :-)
[nbits with hex...]
brrr. I can't help it. This should really be a quick log2()
call, this seems quite natural to me.
> > time=time.time
>
> Not time.clock?
One of my very first scripts. Don't ask me, perhaps an old
Python/Win32 problem.
> > # the end ----------------------------------------------------
>
> looks-more-like-the-beginning-to-me<wink>-ly y'rs - tim
It was impressive for me to see that with an hour of
coding, I can outperform C in Python. (Actually is was some
more hours, since I didn't remember how, and undusted some books).
My intent was to encourage people to play with stuff like
this, and to see how it can continue. The last word btw. keeps
me busy enough :)
At the same time, the reacton of the list was a little
disappointing to me. There seems to be near to zero interest.
Is this the hot summer? Holidays? Bad title of my post?
Then I also had a look into NumPy and its matrix multiplication.
I stunned. Again no attempt to optimize big O. It could be done
with Python, the algorithms are known, the fruit is waiting
to be picked.
Challenge: I claim Python could be a competitor to Mathematica,
if only some people were interested.
*** Hey Python, where are your ambitioned Math students? ***
(did they wake up, finally?)
waiting-for-the-beginning-to-show-up-ly y'rs - chris
Ivan Van Laningham
Christian Tismer wrote:
>
> Tim Peters wrote:
>
[bobbit]
> At the same time, the reacton of the list was a little
> disappointing to me. There seems to be near to zero interest.
> Is this the hot summer? Holidays? Bad title of my post?
>
Oh, bad title, most certainly! ;-) *I'm* interested! I don't
understand most of what you guys are talking about, but I can see that
if I could afford to spend the time studying it I would be rewarded.
On topic, though, doesn't Knuth include in his "How Fast can We
Multiply?" section a discussion of using modular arithmetic to speed up
multiplies by a pretty large factor? Now, maybe I've misunderstood the
discussion so far so badly that I missed your use of it, but I don't
recall seeing anything about it.
> Then I also had a look into NumPy and its matrix multiplication.
> I stunned. Again no attempt to optimize big O. It could be done
> with Python, the algorithms are known, the fruit is waiting
> to be picked.
>
> Challenge: I claim Python could be a competitor to Mathematica,
> if only some people were interested.
>
I would be interested--just not competent! But a Python version of
Mathematica could indeed be a killer app. (Or is that a killer rabbit?)
<guess-who-flunked-high-school-algebra>-ly y'rs,
Ivan
----------------------------------------------
Ivan Van Laningham
Callware Technologies, Inc.
iva...@callware.com
iva...@home.com
http://www.pauahtun.org
See also:
http://www.foretec.com/python/workshops/1998-11/proceedings.html
Army Signal Corps: Cu Chi, Class of '70
----------------------------------------------
Tim Peters
> Oh, bad title, most certainly! ;-) *I'm* interested! I don't
> understand most of what you guys are talking about, but I can see that
> if I could afford to spend the time studying it I would be rewarded.
Here's the heart of it. Suppose you want to multiply 10*A+B by 10*C+D,
where A, B, C, D are all in range(10). The answer is
100*A*C + 10*(A*D + B*C) + B*D
Multiplying by 100 or 10 is just tacking zeroes on the end, so that's cheap.
Adding is cheap too. Multiplying is expensive (well, imagine you already
understood what follows so that A, B, C, D each had hundreds of digits
instead of one <wink>).
Now time for a trick: consider (A+B)*(C+D). That's A*C+A*D+B*C+B*D. So
the answer above is the same as
100*A*C + 10*((A+B)*(C+D)-A*C-B*D) + B*D
and there are only *three* unique products there instead of 4: A*C,
(A+B)*(C+D), and B*D. That means we can compute the product of two B-bit
numbers via 3 multiplications of B/2-bit numbers, plus some cheap stuff.
The same trick applies to the B/2-bit numbers, down & down. In the end that
implies we can do multiplication in time proportional to B**log2(3) instead
of B**2.
The method Chris actually used computes (A-B)*(D-C) instead for obscure
technical reasons (history, plus that subtracting instead of adding reduces
the size of the multiplicands a little on average).
> On topic, though, doesn't Knuth include in his "How Fast can We
> Multiply?" section a discussion of using modular arithmetic to speed up
> multiplies by a pretty large factor? Now, maybe I've misunderstood the
> discussion so far so badly that I missed your use of it, but I don't
> recall seeing anything about it.
We didn't mention it. Python doesn't store numbers in modular form, so it's
not that interesting -- it's expensive to convert in to and out of modular
form. If it kept them in modular form, then multiplication would be
essentially linear-time. If you want that-- and can live with numbers in
modular form over long stretches of calculations --it's easy enough to write
a class to do that.
> I would be interested--just not competent! But a Python version of
> Mathematica could indeed be a killer app. (Or is that a killer rabbit?)
I thought Mayan calendar analysis was Python's killer app!
or-is-that-an-app-killer<wink>?-ly y'rs - tim
Tim Peters
> ...
> Also I was tempted to use "&" on the numbers, in order to save
> one nbits call. It was just unclear what happens with negative
> numbers and dropped it.
It's clear enough what happens with positive numbers, though <wink>:
m = 1L << 10000
n = m - 1
print m & n # prints 0
Use "|" instead and you're golden. But that takes time proportional to the
number of bits in the larger input, so-- again --we need to fix nbits. Then
nx = nbits(x)
ny = nbits(y)
if nx == 0 or ny == 0:
return 0
would be constant-time.
> Was happy that negatives went correct without special treatment.
Yup! It's a delight.
> ...
> [nbits with hex...]
>
> brrr. I can't help it. This should really be a quick log2()
> call, this seems quite natural to me.
It should be constant-time.
> ...
> At the same time, the reacton of the list was a little
> disappointing to me. There seems to be near to zero interest.
> Is this the hot summer? Holidays? Bad title of my post?
Na, I think "little interest" is it. Outside of crypto, prime-number fans,
and people studying Calendars from Hell <wink>, nobody cares about really
big ints. The first two groups-- who care a lot --care *so* much they're
using hyper-optimized bigint pkgs anyway. Andrew Kuchling's wrapping gmp
for Python is the most useful thing that's ever been done in this area.
> Then I also had a look into NumPy and its matrix multiplication.
> I stunned. Again no attempt to optimize big O. It could be done
> with Python, the algorithms are known, the fruit is waiting
> to be picked.
Large dense matmult isn't very common. Sing along:
Lemon tree, very pretty,
And de lemon flower is sweet.
But de fruit of de poor lemon,
Is impossible to eat.
> Challenge: I claim Python could be a competitor to Mathematica,
> if only some people were interested.
Give me one bored programmer with a few weeks to spare, and you can keep a
whole web full of people who are interested <wink>.
> *** Hey Python, where are your ambitioned Math students? ***
>
> (did they wake up, finally?)
>
> waiting-for-the-beginning-to-show-up-ly y'rs - chris
i'll-join-you-in-the-bar-ly y'rs - tim
Christian Tismer
Tim Peters wrote:
>
> [Christian Tismer]
> > ...
> > Also I was tempted to use "&" on the numbers, in order to save
> > one nbits call. It was just unclear what happens with negative
> > numbers and dropped it.
Waaah, did I say "&"? I used "|" or course.
> It's clear enough what happens with positive numbers, though <wink>:
>
> m = 1L << 10000
> n = m - 1
> print m & n # prints 0
>
> Use "|" instead and you're golden. But that takes time proportional to the
> number of bits in the larger input, so-- again --we need to fix nbits. Then
But your "min" suggestion makes still a lot more sense to me.
> > Was happy that negatives went correct without special treatment.
>
> Yup! It's a delight.
My first attempt had special "-" treatment, until I recognized
that the shift and mask operations already do all we need. :-)
...
> Na, I think "little interest" is it. Outside of crypto, prime-number fans,
> and people studying Calendars from Hell <wink>, nobody cares about really
> big ints. The first two groups-- who care a lot --care *so* much they're
> using hyper-optimized bigint pkgs anyway. Andrew Kuchling's wrapping gmp
> for Python is the most useful thing that's ever been done in this area.
I see. The gmp makes very much sense. Until now.
On the other hand, if the hard coded "C" basic elements are
made fine and tiny, they give us nuts and bolts to build
the real algorithms on top.
We can therefore do what gmp does, not so much slower, but
readable, in Python. No porting, no nothing, just using.
And again I see a win for Python. It is able to show how
algorithms work, in a very clear way using its clean
syntax, and you can use the implementation of math directly
as highschool text. *And* it just works, and is fast.
Fells very very promising to me.
[matmult]
> Large dense matmult isn't very common. Sing along:
[lemon tree - sung]
It is just so easy. Factoring matrix multiplication
into its submatrices and then save (I think it was one)
multiplications - it would be a dlight to write in
Python. Not doing so appears to me as to know
quicksort, but using bubble sort all the time,
not just at break-even (10-20).
But I'm telling nothing new to you.
> > Challenge: I claim Python could be a competitor to Mathematica,
> > if only some people were interested.
>
> Give me one bored programmer with a few weeks to spare, and you can keep a
> whole web full of people who are interested <wink>.
I'm bored enough, just the few weeks are missing. Sigh...
> > *** Hey Python, where are your ambitioned Math students? ***
> >
> > (did they wake up, finally?)
> >
> > waiting-for-the-beginning-to-show-up-ly y'rs - chris
>
> i'll-join-you-in-the-bar-ly y'rs - tim
Reading the last posts, I think we aren't as alone as I thought.
As another little challenge:
Tim already mentioned Jurjen N.E. Bos' high precision math
library. It is in Python's contrib site.
I have the strong impression that it was not used by anybody
since a year now. python.org fails also to advertize it.
It could need a little rework, to make use of 1.5.2's features.
And if someone marries it with a fast multiplication, that
could make a neat tool.
If I were a student, I would propose my advisor to let me do such
things as my semester homework. What a delight :-)
ciao - chris
Les Schaffer
C> Then I also had a look into NumPy and its matrix multiplication. I
C> stunned. Again no attempt to optimize big O. It could be done with
C> Python, the algorithms are known, the fruit is waiting to be picked.
T> Large dense matmult isn't very common. Sing along:
i dont see the connection between floating point matrix multiplies and
the integer arithmetic tim sketched out in the previous post.
Spock, explain!
les 'my-denseness-is-multiplying' schaffer
--
____ Les Schaffer ___| --->> Engineering R&D <<---
Theoretical & Applied Mechanics | Designspring, Inc.
Center for Radiophysics & Space Research | http://www.designspring.com/
Cornell Univ. scha...@tam.cornell.edu | l...@designspring.com
Christian Tismer
Les Schaffer wrote:
>
> Tim and Christian spoke:
>
> C> Then I also had a look into NumPy and its matrix multiplication. I
> C> stunned. Again no attempt to optimize big O. It could be done with
> C> Python, the algorithms are known, the fruit is waiting to be picked.
>
> T> Large dense matmult isn't very common. Sing along:
>
> i dont see the connection between floating point matrix multiplies and
> the integer arithmetic tim sketched out in the previous post.
>
> Spock, explain!
The floating point property is not the point.
It is just a similar strategy : divide and conquer.
Multplying long ints is breaking them into vectors of smaller
size, and then figure out where a mult can be saved by some
more adds.
The way for matrices is factorizing the multiplicants
into smaller ones and then figuring out which mults
can be saved.
And for both also applies: The low level math for the small
pieces which are below a threshold, are still done by the
fast, "dumb" C routines.
Elegant problem reduction is done with Python which controls
that, and saves calls to the low level.
Jeffrey Chang
> [Christian Tismer]
> > Then I also had a look into NumPy and its matrix multiplication.
> > I stunned. Again no attempt to optimize big O. It could be done
> > with Python, the algorithms are known, the fruit is waiting
> > to be picked.
>
> Large dense matmult isn't very common.
Oh, shoot. I'm 2 days into a week long calculation of normalizing a dense
~6000x6000 matrix. Do you mean I'm the only person doing this kind of
thing? ;)
I've been toying around with the idea of working on optimized matrix
functions for NumPy (while waiting for the computations to finish!).
Since I often work with large matrices, there could be big wins for me.
However, it's a little disheartening to hear that few other people would
benefit from it.
I know the O(n**log2(7)) algorithm for matrix multiply and am aware that
there are faster ones. Does anyone know what they are? Can anyone
recommend a good general book for optimized matrix algorithms? Is Knuth 2
the best source?
Thanks,
Jeff
Tim Peters
> i dont see the connection between floating point matrix multiplies and
> the integer arithmetic tim sketched out in the previous post.
Knuth Volume 2 is a good reference, for this and all other vital questions
in life <wink>.
Turns out it's possible to multiply two 2x2 matrices using 7 multiplies.
Since matmult can be "blocked", a larger matrix can be treated as a 2x2
array of smaller submatrices, and the same trick recursively applied. In
the end you get a method with runtime O(n**log2(7)) ~= O(n**2.81) instead of
O(n**3).
The intmult trick is less work to code, and yields a bigger payoff earlier.
much-like-python-itself-ly y'rs - tim
Tim Peters
>> Andrew Kuchling's wrapping gmp for Python is the most useful thing
>> that's ever been done in this area.
[Christian Tismer]
> I see. The gmp makes very much sense. Until now.
>
> On the other hand, if the hard coded "C" basic elements are
> made fine and tiny, they give us nuts and bolts to build
> the real algorithms on top.
>
> We can therefore do what gmp does, not so much slower, but
> readable, in Python. No porting, no nothing, just using.
I play with & investigate bigint algorithms in Python, & I enjoy it. But
serious work in this area requires serious speed, and GMP doesn't compromise
on that -- from slaved-over assembler to low-level mutable packed storage
limbs, it's doing the right stuff for its task. Factors of 5 count here
<wink>.
You do "long1 & long2" in Python, and it's handled by an all-purpose bitop
function that redispatches on the operation in each iteration of the inner
loop. It's "like that" everywhere: the bigint code was written for
compactness (longobject.c is only twice the size of intobject.c!), and
no-hassle porting to the feeblest machine on earth. The compromises have
consequences (some good! that code has been rock solid since the first
release, while it can take a week to track down all the bugfix patches to
the last official GMP release).
> And again I see a win for Python. It is able to show how
> algorithms work, in a very clear way using its clean
> syntax, and you can use the implementation of math directly
> as highschool text. *And* it just works, and is fast.
> Fells very very promising to me.
For educational purposes, sure -- provided you can find enough people
seeking this kind of education <wink>.
[matmult]
>> Large dense matmult isn't very common. Sing along:
> It is just so easy. Factoring matrix multiplication
> into its submatrices and then save (I think it was one)
> multiplications - it would be a dlight to write in
> Python. Not doing so appears to me as to know
> quicksort, but using bubble sort all the time,
> not just at break-even (10-20).
> But I'm telling nothing new to you.
OTOH, Winograd's refinement of the matmult procedure requires 15 floating
adds, where straight 2x2 matmult requires only 4; float + isn't much cheaper
than float * on many modern processors. So the win here isn't nearly as
pure as for intmult, and the reduction is only from O(n**3) to O(n**2.81)
operations anyway. There are reasons few vendors offer this form of matmult
in their libraries <0.5 wink>.
> ...
> As another little challenge:
>
> Tim already mentioned Jurjen N.E. Bos' high precision math
> library. It is in Python's contrib site.
> I have the strong impression that it was not used by anybody
> since a year now. python.org fails also to advertize it.
>
> It could need a little rework, to make use of 1.5.2's features.
> And if someone marries it with a fast multiplication, that
> could make a neat tool.
I use it often; I'm not sure how practical it is <wink>, but it's great fun!
It would benefit a *lot* from a constant-time nbits. The clearest quick win
is to get rid of all the code that's special-casing for "big" left shift
counts; a patch to remove that restriction went into Python shortly after
Jurjen vanished.
> If I were a student, I would propose my advisor to let me do such
> things as my semester homework. What a delight :-)
It would make a great project for any number-head. OTOH, Jurjen's approach
was both highly clever and highly idiosyncratic, and the internals of the
code are difficult to understand. Some of the examples in the docstring no
longer work. And some of the algorithms are clearly flawed; e.g.,
>>> a = r(1)
>>> b = r(1)
>>> a.mant = 8; a.expon = 0
>>> b.mant = 3; b.expon = 0
>>> a
8.+-1
>>> b
3.+-1
>>> a+b
8.+-4
>>>
That is, the result interval is too small -- a+b *can* be as large as 13. A
similar thing happens whenever the before-rounding internal mantissa of a
sum is congruent to 3 mod 4.
A good semester project would be to come up with-- and prove
correct --algorithms for {+ - * / sqrt exp log input output} for numbers in
this format (a very compact flavor of variable-precision floating-point
intervals).
unary-minus-is-easy<wink>-ly y'rs - tim
Christian Tismer
Tim Peters wrote:
>
> [Tim]
> >> Andrew Kuchling's wrapping gmp for Python is the most useful thing
> >> that's ever been done in this area.
>
> [Christian Tismer]
[wants to rebuild the world using nuts & bolts]
> I play with & investigate bigint algorithms in Python, & I enjoy it. But
> serious work in this area requires serious speed, and GMP doesn't compromise
> on that -- from slaved-over assembler to low-level mutable packed storage
> limbs, it's doing the right stuff for its task. Factors of 5 count here
> <wink>.
Agreed. Still, if we resolve to low level functions which deal
with an optimized number of bits at a time (some reasonable cutoff),
then the Python overhead to operate intelligently with them
is (hopefully) small enough to be justified.
> You do "long1 & long2" in Python, and it's handled by an all-purpose bitop
> function that redispatches on the operation in each iteration of the inner
> loop. It's "like that" everywhere: the bigint code was written for
> compactness (longobject.c is only twice the size of intobject.c!), and
> no-hassle porting to the feeblest machine on earth. The compromises have
> consequences (some good! that code has been rock solid since the first
> release, while it can take a week to track down all the bugfix patches to
> the last official GMP release).
Yes. I think this simplicity is a real advantage.
The Python generalism to do the general &-code with
typechecks and everything on every inner loop appears
to be expensive for small objects only.
If I can write algorithms in a way that they operate
only above some "cutoff" level (say 5000 bit ints),
this is again not expensive.
...
[matmult]
> OTOH, Winograd's refinement of the matmult procedure requires 15 floating
> adds, where straight 2x2 matmult requires only 4; float + isn't much cheaper
> than float * on many modern processors. So the win here isn't nearly as
> pure as for intmult, and the reduction is only from O(n**3) to O(n**2.81)
> operations anyway. There are reasons few vendors offer this form of matmult
> in their libraries <0.5 wink>.
15 adds, sure. Would be expensive for a 2x2 matrix of reals.
But since we recusrively apply this to 2x2 matrices of submatrices,
the complexity of multiplication quickly voidens the real-processors
optimization again, and we win quickly.
Thinking of Jeffrey Chang's 6000x6000 matrices, what can we expect?
I didn't try, but let's assume a cutoff at about 94x94 matrices,
which appears fairly high for me.
Jeffrey can than think of the multiplication of 64x64 matrices,
where the scalars are the 94x94 matrices of above.
Cost for addition of two 94x94 matrices vanishes, compared to
multiplication.
So Jeffrey might expect to go from 64 ** 3 == 262144
down to 64 ** ( log(7)/log(2) ) == 117649.0 which gives him
a speedup factor of 2.23 .
If we are more lucky, and the cutoff enables us to 47*47 matrices,
the speedup will be 2.55 . And for 24*24 matrices, we end up at
a speedup of 2.91 .
May we expect a win for 12x12 matrices, already?
It would give us a great speed gain of 3.33, but we have to
guess the cutoff.
Reading Knuth about that, the necessary scaling for stability
implies that this techniques pay off for n >= 40, so if
we are lucky, the 47*47 case above could be implemented,
hopefully giving a factor of between 2and 2.5 .
One concern I have, and which stops me to dig seriously into
this is matrix stability. Since we are dealing with floats in
the low level, the result may be impacted by more rounding
errors, since more additions are involved.
At least this has to be studied carefully before replacing
one algorithm by another. I don't know what kind of
preconditioning is necessary here, but Knuth said it is.
Not so with fastlong, this is a perfect replacement.
I have not further investigated if Fourier transform
can be applied, and if matrix multiplication can be
expressed in terms of convolutions and Fourier transforms.
Multiplication of very long integers *can* be done this way.
[and about Jurjen's exact math, we'll wait for volunteers]
Les Schaffer
> Turns out it's possible to multiply two 2x2 matrices using 7
> multiplies. Since matmult can be "blocked", a larger matrix can be
> treated as a 2x2 array of smaller submatrices, and the same trick
> recursively applied. In the end you get a method with runtime
> O(n**log2(7)) ~= O(n**2.81) instead of O(n**3).
okay, THAT stuff...
Numerical Recipes attributes this to Strassen -- i've seen some stuff
on the web on it too:
http://www.cs.engr.uky.edu/~lewis/cs-alg/materials/Strassen.html
http://www.compapp.dcu.ie/~aileen/research/matmult4.html
http://www.cs.duke.edu/~luoma/strass/strass.html
i'll grab a look at Knuth ....
> The intmult trick is less work to code, and yields a bigger payoff
> earlier.
yeah, and N has to be large enough to justify using strassen ...
les
tim...@email.msn.com
> Oh, shoot. I'm 2 days into a week long calculation of normalizing
> a dense ~6000x6000 matrix. Do you mean I'm the only person doing
> this kind of thing? ;)
Yes <wink>.
> I've been toying around with the idea of working on optimized
> matrix functions for NumPy (while waiting for the computations to
> finish!). Since I often work with large matrices, there could
> be big wins for me. However, it's a little disheartening to hear
> that few other people would benefit from it.
>
> I know the O(n**log2(7)) algorithm for matrix multiply and am
> aware that there are faster ones. Does anyone know what they are?
> Can anyone recommend a good general book for optimized matrix
> algorithms? Is Knuth 2 the best source?
It's probably the best accessible source for a quick overview of the
theory.
In practice, though, you're pretty much hosed: I don't know what your
6000x6000 matrices *contain*, but if it's doubles or even floats then
each consumes a significant fraction of a gigabyte. So unless you're
running on a real-memory machine, cache effects are likely as least as
important as operation count; on many platforms, much more important.
When high performance requires coding to the most obscure details of a
platform's memory hierarchy, portability goes out the window.
There's an interesting paper about how all this relates to Strassen-
Winograd matmult (the simplest of the "advanced" methods) at:
http://www.cs.duke.edu/~alvy/papers/sc98/index.htm
Here's part of the abstract:
Strassen's algorithm for matrix multiplication gains its lower
arithmetic complexity at the expense of reduced locality of
reference, which makes it challenging to implement the algorithm
efficiently on a modern machine with a hierarchical memory
system. We report on an implementation of this algorithm that
uses several unconventional techniques to make the algorithm
memory-friendly.
...
Performance comparisons of our implementation with that of
competing implementations show that our implementation often
outperforms the alternative techniques (up to 25%). However, we
also observe wide variability across platforms and across matrix
sizes, indicating that at this time, no single implementation is
a clear choice for all platforms or matrix sizes.
...
Par for the course.
buy-an-old-cray-while-they're-still-cheap<wink>-ly y'rs - tim
Sent via Deja.com http://www.deja.com/
Share what you know. Learn what you don't.