=  David Harvey: would like to experiment further with speeding up object construction =

== gmpy lessons??? ==

Do "sage -i gmpy-1.0.1".  Then you can see the following (*ignore* the wall time):



{{{
import gmpy
}}}

{{{
a = gmpy.mpz(9393r); b = gmpy.mpz(1293r)
}}}

{{{
%time
for i in range(10^6):
    c = a * b
///
CPU time: 0.38 s,  Wall time: 0.53 s
}}}

{{{
a = 9393; b = 1293
}}}

{{{
%time
for i in range(10^6):
    c = a * b
///
CPU time: 0.78 s,  Wall time: 1.16 s
}}}

{{{
a = 9393r; b = 1293r
}}}

{{{
%time
for i in range(10^6):
    c = a * b
///
CPU time: 0.37 s,  Wall time: 0.66 s
}}}

====================================


{{{
import gmpy
}}}

{{{
a = ZZ.random_element(2^256); b = ZZ.random_element(2^256)
}}}

{{{
a = gmpy.mpz(int(a)); b = gmpy.mpz(int(b))
}}}

{{{
%time
for i in range(10^6):
    c = a * b
///
CPU time: 0.59 s,  Wall time: 0.86 s
}}}

{{{
a = ZZ(long(a)); b = ZZ(long(b))
}}}

{{{
%time
for i in range(10^6):
    c = a * b
///
CPU time: 1.12 s,  Wall time: 1.41 s
}}}

{{{
a = long(a); b = long(b)
}}}

{{{
%time
for i in range(10^6):
    c = a * b
///
CPU time: 1.62 s,  Wall time: 2.15 s
}}}




== Email announcing the project from David Harvey ==
 
hi people,
 
Just want to float an idea for discussion and possibly a coding  
sprint at SD3.
 
Some background: Today I did some work on speeding up getting  
coefficients out of NTL objects, specifically polynomials in Z[x].  
It's much improved now; when you request a coefficient of an NTL ZZX  
object, it copies the bytes directly into a new Integer object,  
instead of what it used to do (which went via a C string in decimal,  
and a python string, and a python long, etc etc).
 
But still what is taking a lot of time is constructing the Integer  
object. In fact, it's quite embarrasing: it takes us about half as  
long to construct 100,000 Integer objects as NTL takes to *multiply*  
two polynomials with 100,000 small coefficients:
{{{
sage: time for i in range(100000): pass
CPU times: user 0.05 s, sys: 0.00 s, total: 0.05 s
 
sage: time for i in range(100000): x = None
CPU times: user 0.09 s, sys: 0.00 s, total: 0.09 s
 
sage: time for i in range(100000): x = int()
CPU times: user 0.16 s, sys: 0.00 s, total: 0.17 s
 
sage: time for i in range(100000): x = Integer()
CPU times: user 0.36 s, sys: 0.00 s, total: 0.36 s
 
sage: f = PolynomialRing(ZZ, "x")([ZZ.random_element() for _ in range 
(100000)])
sage: time g = f*f
CPU times: user 0.76 s, sys: 0.02 s, total: 0.79 s
}}}
This is despite all the work we put into this at SD2.
 
It would be good to be able to optimise object construction in  
general. Unfortunately I think the general case is a very difficult  
problem. Anyone who worked on this at SD2 will agree, I'm sure 
 
On the other hand, I would wager that construction of Integer objects  
is by far the most important. So perhaps we should give up some  
beauty and unity of code to just get Integers working damn fast. So  
here's what I propose: at SD3, let's try writing an experimental pure  
C function for constructing Integers that gets inserted into whatever  
tp_xyz slot is appropriate. I don't care if it has to deal directly  
with mangled pyrex names or whatever. From memory, all it needs to do  
is:
(1) reference counting on the integer ring (ha ha we could even skip  
this if we could guarantee that no-one else ever resets the parent,  
and that there is always at least one reference to the integer ring  
lying around somewhere)
(2) malloc some space for the actual python object
(3) fill in some fields, like the TypeObject*
(4) mpz_init
 
Let's skip all the function calls, all the crap that pyrex puts in,  
etc etc. Basically the only stuff left that will really suck up time  
is the two malloc calls. We could even try writing a buffering system  
for mallocing space for a whole bunch of Integers at once, if that  
proves to be taking up time.
 
Similarly we would need a destructor.
 
I'm sure it's not quite as simple as what I've put in the above list,  
but let's just see what we can do. I would like it to be as close as  
possible to the int() construction time. We still are a factor of > 2  
away.
 
David


== Email back from William ==
{{{
> (1) reference counting on the integer ring (ha ha we could even skip
> this if we could guarantee that no-one else ever resets the parent,
> and that there is always at least one reference to the integer ring
> lying around somewhere)

There is my integer with the integer ring -- it should be immutable and
create at module load time and there should only ever be exactly one
copy of it.  I think we should definitely be allowed to forgot about
reference counting for it. 

> (2) malloc some space for the actual python object
> (3) fill in some fields, like the TypeObject*
> (4) mpz_init

You should put 
  (0) or (5) object pool
as an important step -- this "object pool" idea is one of the tricks
that Python uses for its ints.  For example, in your benchmark:

  time for i in range(10^5): x = int()

Python is looking up and returning exactly the same int (the 0) every time:
  sage: int() is int()
  True
  sage: a = 999038r; b=999038r
  sage: a is b
  True

In contrast, when you do Integer(), so is creating a new integer object
every time, and probably (?) also freeing one:
  sage: Integer() is Integer()
  False
  sage: a = 999038; b=999038
  sage: a is b
  False


> Let's skip all the function calls, all the crap that pyrex puts in,
> etc etc. Basically the only stuff left that will really suck up time
> is the two malloc calls. We could even try writing a buffering system
> for mallocing space for a whole bunch of Integers at once, if that
> proves to be taking up time.
>
> Similarly we would need a destructor.

Or return objects to the pool -- this can also speed up desctructing,
since you just don't do it!

> I'm sure it's not quite as simple as what I've put in the above list,
> but let's just see what we can do. I would like it to be as close as
> possible to the int() construction time. We still are a factor of > 2
> away.

This project is very very well worth pursing. 

  -- William
}}}

== Some Ideas from Robert Bradshaw ==
{{{
Some thoughts: I think a pool is a very important idea to consider. I
can think of two instances where 1000's of integer objects would be
created: first, in some large object such as a matrix or polynomial
(in which case there should be a specialized type) and second in some
huge loop (in which case a pool would help immensely).

Also, it'd be interesting to look at the distribution, but I wouldn't
be surprised if the majority of integers (ephemerally) created were
relatively small--say < 100s. Zero and one especially are used all
over. Similar to the pool idea, it might be worth allocating the
first 100 integers and whenever you want to create a "small" integer,
it would simply return one of these. (I think small one-limb mpz_t's
can be detected very easily with mpz_size and a bit mask.) Of course,
using python ints might be in order for many of these cases too.

A related idea came up in the discussion we had here on linear
algebra. Right now if one wants to optimize linear algebra over a new
ring one must re-implement matrix multiplication, addition, etc. The
generic algorithms request an entries (as a Python object), perform
the arithmetic, then store the resulting python object. This can be
hugely inefficient. Rather, what if the matrix had void* methods
_get_unsafe_raw(i,j) and _set_unsafe_raw(i,j), and the corresponding
ring had _add_raw(), _mul_raw(), etc. Also, the ring could have
_get_raw() and _create_from_raw(). For the integer ring, these would
return mpz_t* and, for instance, _mul_raw() could even be a macro to
mpz_mul. The generic base case would just pass around python objects.
The "reference counting" for these raw results would have to be done
manually. I would suggest giving them the same semantics as gmp. This
way one could implement generic polynomial/matrix/etc algorithms that
would be able to operate efficiently on any ring with the above
methods. For some cases (such as the integers) one would want actual
specialized matrices, etc. but it would greatly reduce the work to
get significant speedup for objects containing many elements of a
given generic ring. Also, it would make implementations for specific
ring elements easier to swap in and out (without having to change all
the types that access the element internals).
}}}