Re: Bignum multiplication in Python vs GMP (was: constexpr is really very smart!)

On Wed, 25 Dec 2024 11:05:05 +0200
Michael S <already5chosen@yahoo.com> wrote:

On Tue, 24 Dec 2024 05:21:42 -0800
Tim Rentsch <tr.17687@z991.linuxsc.com> wrote:
 

Michael S <already5chosen@yahoo.com> writes:
 

On Sun, 22 Dec 2024 10:25:59 -0800
Tim Rentsch <tr.17687@z991.linuxsc.com> wrote:
  

Michael S <already5chosen@yahoo.com> writes:
  

On Fri, 20 Dec 2024 14:59:07 -0800
Tim Rentsch <tr.17687@z991.linuxsc.com> wrote:
  

I ran a python version to compute fibonacci( 10000000 ).  This
code ran 30% faster than the previous version, which I would be
is almost entirely due to the expression factoring with 'c'.
Great!   

>
This one does fib(10M) in 89 msec on my very old home PC.
>
[code]   

>
You are the speed king!   

>
But approximately the same code in python is MUCH slower.
2.9 sec on 11 y.o. PC that still serves me as a desktop at work.
2.0 sec on 5.5 y.o. mall server.   

 
That matches my experience with doing this in python.
 

I don't really understand why.
I expected that nearly all time is spend in multiplication of few
huge numbers, probably over 75% of time in last couple of steps (6
multiplications).  And I was under impression that internally
Python uses the same GMP library that I used in C. So I expected
the same speed, more or less.  May be, 1.5x slowdown because of
less efficient memory handling in Python. 35x difference is a big
surprise.   

 
Yes, I have no explanation either.  Your reasoning looks sound to
me.
 

 
I did few more time measurements both with GMP and with Python. It
clearly showed that they use not just different implementations of
bignum multiplication, but different algorithms.
GMP multiplication for operands in range of above ~200,000 bits is
approximately O(N*log(N)) with number of bits, or just a little bit
slower than that.
Python multiplication scale much worse, approximately as (N**1.55).
 
And indeed after binging more I had seen the following discussion on
reddit:
https://www.reddit.com/r/Python/comments/7h62ul/python_largeinteger_computation_speed_a_comparison/
 
A relevant quote:
"Python seems to use two algorithms for bignum multiplication: "the
O(N**2) school algorithm" and Karatsuba if both operands contain more
than 70 decimal digits. (Your example reaches 70 decimal digits
already after the 5th(?) squaring, and after the 22nd there are
~17million decimal digits(?). The wast majority of the time is spent
on the last one or two multiplications i.e. Karatsuba.)
 
Meanwhile GMP implements seven multiplication algorithms: The same(?)
"school" and Karatsuba algorithms, plus various Toom and an FFT method
(Strassen), but I'm unclear on what the exact thresholds are that
select these. It seems these thresholds are in "limbs" ("the part of a
multi-precision number that fits in a single machine word") instead of
decimal digits, with e.g. FFT method being used after "somewhere
between 3000 and 10000 limbs", i.e. >~30k - 100k(?) decimal digits, so
it would probably be used for the relevant steps in your example."
 
 
It seem that that threshold for FFT method mentioned in the quote
applies to number of limbs (i.e. 64-bit groups of bits in our specific
case) in the product rather than in the operands. Also it is possible
that in the newer versions of GMP the threshold is a lower than it
was 7 years ago.
 
I wonder, why python uses its own thing instead of GMP.
Incompatible licenses or plain NIH ?
 

Forgot to mention.
As said in the same redit thread, one can use GMP from python. It's
pretty easy, just not the default.

import gmpy2

def fib(n):
  if n < 1:
    return 0

  tmp = n
  msb = 1
  while tmp > 1:
    msb += msb
    tmp //= 2
  f0=gmpy2.mpz(0); f1=gmpy2.mpz(1);
  while msb > 1:
    ff0 = f0*f0 + f1*f1
    f1 = (f0+f0+f1)*f1
    f0 = ff0
    msb //= 2
    if msb & n:
      f0 = f1
      f1 += ff0
  return f1


For big values of n gmpy2-based variant calculates Fibbonaci numbers at
approximately the same speed as C variant. For smaller numbers - not
quite the same as C and for n < ~10,000 it is slower than default
python math.