Base of np.log1p

[gokceneraslan]

np.log1p is a great function but not having alternatives with different bases (like log2 and log10) is very limiting. I wonder if it'd be possible to add other functions like log2_1p and log10_1p (and possibly exp2_m1 and exp10_m1).

I played around with some numba-based solutions such as this one based on the original implementation in numpy:

@vectorize([float64(int32),
            float64(int64),
            float32(float32),
            float64(float64)])
def log10_1p(x):
    if np.isinf(x) and x > 0: 
        return x
    
    u = 1. + x
    d = u - 1.
    
    if d == 0:
        return x

    return np.log10(u) * x / d

it seems to be working well (see https://nbviewer.ipython.org/gist/gokceneraslan/2744cfeda702fb9b9e48d0216427372c for some results) but having a numpy-native function would be really great.

PS: Can anybody explain why there is a * x / d in the end? Maybe @khinsen can shed some light.

[gokceneraslan]

GLIBC implementation explains it in detail: https://github.com/lattera/glibc/blob/master/sysdeps/ieee754/dbl-64/s_log1p.c#L16

[WarrenWeckesser]

You could create a one-line function for this, e.g.:

In [17]: def log2_1p(x): 
    ...:     return np.log1p(x) / np.log(2) 
    ...:                                                                                    

In [18]: log2_1p(3)                                                                         
Out[18]: 2.0

In [19]: log2_1p(1e-23)                                                                     
Out[19]: 1.4426950408889635e-23

You can make it faster by using the constant value log(2) = 0.6931471805599453:

In [36]: def log2_1p(x): 
    ...:     return np.log1p(x) / 0.6931471805599453 
    ...:                                                                                    

[seberg]

You can suggest adding it. It probably should be run by the mailing list before we actually would do it. But, I think the general sentiment here will be that this should be in scipy.special (or somewhere else out of NumPy). However, the numpy functions should point to such other ones in Documentation.

[gokceneraslan]

@WarrenWeckesser For small x, unlike the log1p(x) = x trick in log1p, that'd return log2_1p(x)=x/log(2). So I wonder if there is a big difference in terms of accuracy between returning x and x/log2 for small x.

[lewlin]

@gokceneraslan a division by a positive number is well-conditioned therefore I think it is quite safe to to divide by log(2).

[WarrenWeckesser]

@gokceneraslan wrote

For small x, unlike the log1p(x) = x trick in log1p, that'd return log2_1p(x)=x/log(2).

Yes, that is the correct result.

For example, in the following I use mpmath to compute log2(1 + 1.5e-24) very accurately:

In [215]: import mpmath

In [216]: mpmath.mp.dps = 50

In [217]: x = 1.5e-24

In [218]: y = mpmath.log(mpmath.mp.one + x, b=2)

In [219]: y
Out[219]: mpf('2.1640425613334450784151339045869271953211073610684447e-24')

In [220]: float(y)
Out[220]: 2.1640425613334452e-24

Now compute the same result using the function that I suggested above:

In [221]: def log2_1p(x): 
     ...:     return np.log1p(x) / 0.6931471805599453  
     ...:                                                                       

In [222]: log2_1p(x)
Out[222]: 2.1640425613334452e-24

[rossbar]

Closing for now as simple alternatives to adding new functions are given in the above discussion. If interested in pursuing the addition of new log*p functions with different bases, consider bringing it up on the mailing list.