Bruce Sherwood Bruce_Sherwood@ncsu....
Sat Dec 29 16:02:36 CST 2007

```I found by timing measurements that a faster scheme with less penalty
for the case of sqrt(array) looks like this:

def sqrt(x):
if type(x) is float: return mathsqrt(x)
return numpysqrt(x)

Bruce Sherwood wrote:
> Roman Yakovenko wrote:
>> On Dec 29, 2007 7:47 AM, Bruce Sherwood <Bruce_Sherwood@ncsu.edu> wrote:
>>
>>> I realized belatedly that I should upgrade from Boost 1.33 to 1.34.
>>> Alas, that didn't cure my problem.
>>>
>> Can you post small and complete example of what you are trying to
>> achieve?
>>
> I don't have a "small and complete" example available, but I'll
> summarize from earlier posts. VPython (vpython.org) has its own vector
> class to mimic the properties of 3D vectors used in physics, in the
> service of easy creation of 3D animations. There is a beta version
> which imports numpy and uses it internally; the older production
> version uses Numeric. Boost python and thread libraries are used to
> connect the C++ VPython code to Python.
>
> vector*scalar, both producing vector. With Numeric, sqrt produced a
> float, which was a scalar for the operator overloading. With numpy,
> sqrt produces a numpy.float64 which is caught by vector*scalar but not
> by scalar*vector, which means that scalar*vector produces an ndarray
> rather than a vector, which leads to a big performance hit in existing
> VPython programs. The overloading and Boost code is the same in the
> VPython/Numeric and VPython/numpy versions. I don't know whether the
> problem is with numpy or with Boost or with the combination of the two.
>
> Here is the relevant part of the vector class:
>
> inline vector
> operator*( const double s) const throw()
> { return vector( s*x, s*y, s*z); }
>
> and here is the free function for right multiplication:
>
> inline vector
> operator*( const double& s, const vector& v)
> { return vector( s*v.x, s*v.y, s*v.z); }
>
> Maybe the problem is in the Boost definitions:
>
> py::class_<vector>("vector", py::init< py::optional<double, double,
> double> >())
>    .def( self * double())
>    .def( double() * self)
>
> Left multiplication is fine, but right multiplication isn't.
>
> A colleague suggested the following Boost declarations but cautioned
> that he wasn't sure of the syntax for referring to operator, and
> indeed this doesn't compile:
>
> .def( "__mul__", &vector::operator*(double), "Multiply vector times
> scalar")
> .def( "__rmul__", &operator*(const double&, const vector&), "Multiply
> scalar times vector")
>
> I would really appreciate a Boost or numpy expert being able to tell
> me what's wrong (if anything) with these forms. However, I may have a
> useful workaround as I described in a post to the numpy discussion
> list. A colleague suggested that I do something like this for sqrt and
> other such mathematical functions:
>
> def sqrt(x):
>   try: return mathsqrt(x)
>   except TypeError: return numpysqrt(x)
>
> That is, first try the simple case of a scalar argument, handled by
> the math module sqrt, and only use the numpy sqrt routine in the case
> of an array argument. Even with the overhead of the try/except
> machinery, one gets must faster square roots for scalar arguments this
> way than with the numpy sqrt.
>
> Bruce Sherwood
>
```