`math.exp` and `math.log` delegate to `__exp__` and `__log__` methods?

Ideas
1

SciPy 1.15 is slated to support basic mathematical operations on random variables. For instance, if X is an object representing a normal random variable, it will support operations like 3*abs(X)**2 + 1; the resulting random variable follows a shifted scaled chi-squared distribution with one degree of freedom. However, if I understand correctly, there is currently no way to override the behavior of math.exp or math.log. To that end, I’d propose that math.exp and math.log attempt to delegate to __exp__ and __log__ methods of the argument.

Surprisingly, I haven’t found much discussion on this topic.
math.log and math.log10 domain error on very large Fractions · Issue #87052 · python/cpython · GitHub suggests:

I guess another approach might be to change the math module so that all its functions support franctions (and decimal, and other esoteric number types). But that would require having log and many others as methods on the numeric types.

and later

There’s the rub. We could create a dunder method for every math function and require that every numeric type to support every dunder method.

I don’t mean to propose going that far. math.floor, math.ceil, and math.trunc already delegate to __floor__, __ceil__, and __trunc__. Allowing math.exp and math.log to delegate, too, might open the floodgates a little further for similar requests (e.g. math.pow to behave like pow), but the precedent is already there.

The exponential and logarithm functions are defined for many mathematical objects, including matrices and quaternions, so there are several other potential use cases besides random variables.

If I understand correctly, there may be hesitation to add the overhead of delegation. However, built-in pow (which does delegate to __pow__) seems just as fast as math.pow for floats, so perhaps this can be done efficiently.

Thanks for considering the suggestion!
Matt

P.S. Mailman 3 Double "at" operator for matmul exponentiation - Python-ideas - python.org mentions that there could be a double “at” operator @@ with associated method __exp__, but this discussion does not appear to be active any more.

6 Likes
2

It might make sense to add __exp__ and __log__ but I think that you will still find that quite limiting: you would immediately find yourself wanting more functions like sin, cos etc.

Ideally there would be a version of the math module that uses @singledispatch for all functions making them extensible to third-party types. That would also hit up against atan2 which requires multiple dispatch though. Dispatchable mathematical functions are much more extensible and could handle everything from int and float to numpy arrays, scipy random variables, mpmath, gmpy2, python-flint’s arb, sympy expressions, jax and automatic differentiation and so on.

The multiple dispatch approach works very well in Julia and makes it easy for many small third-party libraries to build up a big collection of useful stuff that has a coherent user interface rather than every library defining their own functions and not quite compatible interfaces like math.sin, cmath.sin, np.sin, mpmath.sin, sympy.sin, etc.

5 Likes
3

You could imagine going the Array API path and adding a __math_namespace__ property that returns a namespace which has the math functions for the given object.

Then you can have fully generic math functions as well:

def exp(x):
    math = x.__math_namespace__
    return math.exp(x)
2 Likes
4

Currently, most math functions implicitly convert the argument to float, and cmath functions convert the argument to complex. With your proposition, math function will start supporting complex arguments, if complex.__log__ was implemented. What should math.log(-1) return?

2 Likes
5

If you dispatch on type then the caller needs to pass the correct type. This is what Julia does:

julia> log(-1)
ERROR: DomainError with -1.0:
log was called with a negative real argument but will only return a complex result if called with a complex argument.

julia> log(-1+0im)
0.0 + 3.141592653589793im

Instead of having many log functions that accept many types and coerce to their own type you have one log function but the caller needs to pass the intended type.

2 Likes
6

Surely it makes more sense for scipy to just have its own exp()/log() functions? I’ve been using numpy.{exp,log,sin,cos,...}() for years and I’ve never for a second wished that the math module became some weird abstract redirect system so that I could write math.exp() instead (and I don’t think I would write math.exp() even if I could).

6 Likes
7

Surely it makes more sense for scipy

Maybe, but not surely (I think). Why do math.floor, math.ceil, and math.trunc delegate? Is there a difference between those and other math functions that makes it important for only those to delegate, or would it have been better if those did not delegate?

I’m open to the idea that math functions should not be extendable to other input types, but I would like to get a clearer mental model of why some math functions are extendable.

1 Like
8

The numpy functions and the data types that they can handle are sufficient for many people’s use cases. The problem though is that there is no “system” for going anywhere beyond that outside of the existing dunder methods. Consider for example that you can put SymPy expressions into a numpy array:

In [2]: import sympy as sym

In [3]: import numpy as np

In [4]: a = np.array([sym.Symbol('x'), sym.sqrt(2)])

In [5]: a
Out[5]: array([x, sqrt(2)], dtype=object)

In [6]: a + a[::-1]
Out[6]: array([x + sqrt(2), x + sqrt(2)], dtype=object)

So numpy can take advantage of __add__ to add an array of symbolic expressions. There is no dunder for the sin function though:

In [7]: np.sin(a)
...
TypeError: loop of ufunc does not support argument 0 of type Symbol which has no callable sin method

Apparently numpy tries to find a .sin() method but there is no actual convention that a sin method is used to compute sin (such a convention would only be useful in general if int and float have the method but they don’t). Instead if you want the sin of a sympy expression you have to use sympy’s sin function and if you want the sin of a numpy array you have to use numpy’s sin function, and if you use mpmath or jax or … So now if you want the sin of an array of mpmath numbers then mpmath has to have its own array type and so on. Not having a general way to delegate a function like sin limits composability: you can’t just have a numpy array of mpmath numbers etc and you can’t reuse the scipy numerical integration algorithms with mpmath types and so on.

Many users trip up on these things. The reason that the numpy docs always reference np.sin rather than from numpy import sin is because experience showed there was (and still is) so much confusion about having these different inompatible functions of the same name in different modules with people mixing up math.sin, numpy.sin etc.

1 Like
9

What should math.log(-1) return?

I wouldn’t mind it returning 3.141592653589793j, but there is an argument against that at the top of the math documentation. If that is still the preferred behavior, perhaps complex should not implement __log__.

10

Good question. It seems it was for returning an int rather than a float. For more informations read the related PEP:

3 Likes
11

Because @mdickinson & I were both apparently asleep when that change happened.

math was intended to be a thin wrapper around the platform C’s libm, in which ceil / floor / trunc are standard functions, taking and returning double (same as Python float).

And that’s how they worked in Python 2. Why they were changed in Python 3 to return int is a mystery to me. math.ceil(1e300) used to return its argument unchanged, but now returns a 301-digit giant int.

While there are certainly use cases for “get rid of the fractional part in some way and change the type”, such a radical change from the libm behavior for floats wasn’t (IMO) a good way to do that. “Serious” numerical Python 2 code using this generally continued to work, because after “wasting” cycles and memory to construct useless potentially giant ints, in context the result is generally mixed with a float again, and so the giant int would burn more cycles to magically convert back to a float.

Trying to force academic “numeric tower” conceits onto a programming language ignores realities. For example, mathematically, Python’s rationals (fractions.Fraction) are not at all a subset of machine floats. To the contrary, machine floats are a tiny, lumpy subset of Python’s rationals.

That could be repaired if Python’s idea of “float” were changed to implement the constructive reals - but then the language would be so slow that only toy programs would use floats :wink:

17 Likes
12

…I admit I have not understood this part… ^__^’

1 Like
13

It would be useful in many applications to have a generic way to use functions like exp, log etc that was extensible to third party types. I don’t think that trying to add it in the stdlib math module is a good approach though. The idea here that only exp and log could be extensible does not really get us anywhere useful. As long as the math module is mostly inextensible then it is better to leave it like that and use a new module/interface for generic functions.

I think a productive way forward if anyone would like to see some improvement on this would be either the math namespace idea or single/multiple dispatch functions. Either way it would need to be developed outside of the stdlib: what would make it successful is if significant third-party packages could get behind it.

1 Like
14

“lumpy”? I think it’s really off-topic here, although fine for a tiny digression. “Lumpy” is in contrast to “uniform”. For example, there are 2**52 representable floats in range(1, 2), but also the same number in range(2**1000, 2**1001). The spacing between representable floats isn’t uniform.

2 Likes
15

First of all, let me just say that SciPy’s new probability distribution library is one of the most exciting (for me) developments in the Python ecosystem, and if you have a chance, the SciPy people have already done a tremendous job on thorough explanation of their design.

I agree with Tim and others that math is for floats and cmath is for complex numbers. That’s the existing Python design. (I don’t think you should even hook into __floor__, etc. because people doing math.trunc(distribution) by accident might be more convenient if it raised.)

I agree with Benjamin that in an ideal world, there would be some library with dispatched functions that all kinds of types could register on. Function dispatch is truly beautiful. However, that doesn’t exist yet.

I think Nicolas’s suggestion of making SciPy probability distributions an Array API client (or whatever they call it) is interesting, but I don’t see how you could implement many of the Array API functions.

Personally, I think Brénainn’s suggestion of SciPy having its own functions is the most logical way forward:

  • you’ll have complete control over the interface,
  • you can iterate very quickly (nice for a first release), and
  • if you change your mind, you can support some other solution and deprecate your functions fairly painlessly.
8 Likes
16

The idea as I understood it would be to write a whole new “Math API” specification that would be separate from (although likely overlapping and interacting with) the Array API. If the Array API approach works well for downstream consumers to be able to abstract over different array implementations then the Math API could do something similar but focusing on the different implementations of scalars rather than arrays.

1 Like
17

Yes, I was thinking of something that is complementary to and overlapping with the Array API. Maybe it’s a subset of the Array API spec? (I haven’t checked what’s actually in math). Beyond __math_namespace__ (or maybe __math__), does it make sense to have __cmath__ where inputs are coerced to the complex domain? Does it make sense to have the automatic-domain version __emath__, where log(-1.0) returns πi?

Anyway, those are details. I was trying to propose looking at the existing approach, where a well-defined set of functions in a namespace (numpy, here math) is opened up to other implementations by making the namespace part of the object that is being worked on.

1 Like
18

Thanks for considering it, all! It is clear now that math functions were not really intended to be extended to third party types, despite floor/ceil/trunc and in contrast with built-ins. For now, it looks like we’ll stick with exp and log functions in SciPy. Thanks!

4 Likes
19

I think it would by better to make the math functions properties, so you could say x.sin. That way a ‘math’ protocol could be used generically. For functions with an argument like atan2 you would need it to be defined as def atan2(self, x): return x.__ratan2__(self) to do the required double dispatch. Other dunder methods could be reserved for operators only.

20

This isn’t scalable for the large number of functions needed and also isn’t extensible for types you don’t control. Every function you might ever want would need to be a method on float, int, ndarray etc. I doubt that adding even the functions from the math module as float methods would be accepted.

On the other hand if you give up on the idea of generic functions and only want to implement this for one type like SciPy’s random variables then adding .exp and .log methods makes sense.

2 Likes