Add 'half turns' functions (cospi et al) to math module following C23 & IEEE754-2019

Ideas
1

C23 standardized multiple new trig functions in <math.h> that operate on units of “half turns”. These functions are recommended by IEEE754-2019. I propose adding these functions to CPython.

These function can be both more accurate and more performant than the traditional functions which take radian arguments, because inexact range reduction (e.g., modulo 2π) can be replaced by exact calculations (e.g., fmod(x, 2)) that never require excess precision.

Modern Linux systems such as Debian Trixie using glibc 2.41 include an implementation of all these C23 functions (regardless of whether building in C23 mode).

For instance, consider calculating the sine of 500π or 500 half-turns. The result is exactly zero with a proper math.sinpi implementation but nonzero with the traditional math.sin function, because 500π is not an exact FP number:

>>> math.sinpi(500), math.sin(math.pi*500)
(0.0, -1.6070832296378168e-13)

For platforms which do not directly support these functions, for instance because they do not support the C23 standard or do not provide the functions when not building in C23 mode, implementations similar to the existing internal function mathmodule function m_sinpi can be used. These perform range reduction with fmod & special value checking before delegating to a platform radians-based routine (e.g., m_sinpi → sin() or cos()) with necessary scaling by a factor of π.

I have prototyped this addition. You can view it on github: Comparing python:main...jepler:c23-math-pifuns · python/cpython · GitHub

This branch includes new tests & documentation for the new functions.

20 Likes
2

I was a little shocked to discover that sin(n*pi) is 0 only for n = 0. sin(pi) = 1.2246467991473532e-16. No need for n = 500 to find a discrepancy. (I know the example shows than sinpi works for 500, but I presume the same should be true for any integral value.)

+1 for the addition

6 Likes
3

I worry that motivation come from “it’s in the standard”, not from real application. (On this ground, we should consider also other added in the C23 functions, e.g. exp10, exp2m1 and so on.) For instance, scipy has no such functions. mpmath has sinpi/cospi, but nothing else (even internally).

We don’t require C23 yet and can’t rely on presence of those functions; that requires to implement them. Some (inverse) are trivial, some — not. Given that we already have private sinpi (which can be translated to cospi) — I think it’s fine to expose just this pair of functions.

4 Likes
4

It might be worth considering the code provided by the CORE-MATH project.

5

FAQ:

Q1: Why not provide a full library?
The long-term goal of CORE-MATH is to ensure that correctly rounded functions are available to the end-user. We believe that collaborating with existing largely used libraries (instead of competing with them by creating yet another library) is the best way to reach that goal.

I think it should be discussed separately if we should have goal to be better than libm and e.g. bundle core-math code. Doesn’t make much sense to do this just for selected set of functions.

1 Like
6

Based on github code search here are a few projects which have implemented sinpi:

  • pyopencl

  • stdlib-js, math-io and other things in js ecosystem

  • R language

  • JuliaPy

  • OpenCL

  • metal shading language

  • pycrlibm

  • mpfr

github code search just for python files reports 4.1k hits for “sinpi”, most of which look like implementations or use of the sinpi function.

1 Like
7

scipy discussed adding sinpi/cospi in 2024 with many +1s but perhaps it was never actually added: RFC: add sinPi, cosPi, tanPi, and cotPi functions to special - #5 by fancidev - scipy - Scientific Python

1 Like
8

actually it exists but is undocumented or internal:

>>> import scipy.special
>>> scipy.special._ufuncs._sinpi
<ufunc '_sinpi'>
1 Like
9

Yeah, I’m with you on this point. Do we have one convincing example of someone who would use this?

You’re making code harder to read (because most readers aren’t familiar with these functions), and you’re adding a development burden to Python implementations.

Who is going to actually use this? Someone who needs extreme precision? Are they really going to use math instead of numpy? I think if this should be pushed anywhere, it’s either numpy, scipy or array-api-extensions.

Also, if scipy couldn’t be convinced (due to insufficient motivation) to add it yet, I don’t see why CPython should jump to implement it.

1 Like
10

We have an implementation of sinpi() in Modules/mathmodule.c. This is one of the reasons why you can no longer build Python versions older than 3.7 on modern platforms.

1 Like
11

Right. Any integer n other than 0 won’t return 0 for math.sin(n * math.pi), because no product n * pi is exactly representable as a machine float, and modern trig libraries do argument reduction “as if” using an infinite precision value for pi. If such a product were exactly representable, that would imply pi is rational, since all representable machine floats are in fact rationals. That is, if the mathematical n * pi is representable, then n * p = r for some rational r, and so pi = r/n would also be rational.

If not, the implementation of sinpi() would be buggy.

The “use pi as if to infinite precision” argument reduction can be expressive, especially so for “very large” magnitude arguments, all to deliver a nonsense result, like:

>>> import math
>>> math.sin(1e200 * math.pi)
0.9177998132538255

It’s nevertheless (on this box) the closest representable approximation to the true sine of the closest representable float to the product of the closest representable approximation to 1e200 and the closest representable approximation to \pi.

>>> import math, mpmath
>>> mpmath.mp.prec = 1000000
>>> float(mpmath.sin(1e200 * math.pi))
0.9177998132538255  # same thing with a million bits of precision

Most apps don’t care. But some do, and that’s the audience for “half turn” functions. Also for convenience. For example, surveying data is usually expressed in degrees, which have to be converted to radians for most trig functions to process. Nobody remembers how to do that :wink:

But if you have d degrees, then the number of “half rotations” is clearly d/180, so sinpi(d/180) does the job in an easily remembered, faster, and more accurate way.

Of course math.sin(math.radians(d)) also works, but not to hide the accuracy errors:

>>> import math
>>> math.sin(math.radians(180))
1.2246467991473532e-16

So convenient for some, an actual help for some others, but most people won’t use them. So I’m just +0 for the core. Against it is the burden of writing, testing, documenting, and maintaining plenty of new code, for something that’s certainly not a pressing need in the core distro. You’re going to want at least new versions of sin, cos, tan, and their inverses.

4 Likes
12

Yes, but it might be nice not to be worse than glibc on platforms that don’t have glibc.
But glibc is a moving target since they’re adopting more and more functions from CORE-MATH.
So, if CPython decides to ship its own implementation of sinpi(), etc. then I would suggest at least considering the CORE-MATH code (which is MIT licensed).

1 Like
13

Well, kinda, but actually GCC documents that for the math library functions, “The accuracy is unknown.”
Unless a math library says it’s correctly rounded, one must assume it is not. For example, glibc 2.41 (and also glibc 2.43) have input values for which sinpi() and cospi() have an error of more than 1.85 ulps. (The worst error found for Intel’s math library was more than 0.531 for sinpi and more than 0.502 for cospi, so better, but still not correctly rounded, which would be 0.500.) See the double precision table in this paper [Note: Paul Zimmermann is part of CORE-MATH team].
If you’re interested in complex math functions, see this paper [Note: I’m a co-author here, but I am not part of CORE-MATH team].

2 Likes
14

In context, @tjreedy was talking about passing sinpi an argument that’s an exactly representable integer. And that’s what I was addressing. sinpi(int) is exactly 0, and any sinpi() implementation that didn’t deliver 0 would be missing the point of sinpi() rather spectacularly.

The point was about special cases, not about general error bounds. It doesn’t require anything to deliver 0 when passed an int beyond noting that the argument doesn’t have a non-zero fractional part,

1 Like
15

Yes, thanks Tim, of course, you’re right.
Getting back to the original point of the thread, I was surprised that Python did not already have sinpi() and cospi(). I’m not a Python expert, so I cannot speak to how useful these functions would be to the community. I will say that a nice thing I’ve found about using Python is that it seems to me that Python has often anticipated my needs. For example, I felt joy in finding str.capitalize() when I needed it. In that spirit, I think the argument, “If C or C++ has it, then Python should have it” has some merit…
Just my thoughts as somewhat of a newbie here.
I also acknowledge that it’s easy for me to give a +1 to a new feature that I don’t have to implement or maintain. :wink:

4 Likes
16

That’s entirely what drove Python’s math module at the start: it aimed to be a thin wrapper around all and only the “standard parts” of the platform C’s libm.

Life has gotten more complicated ever since :wink: For example, many platform libms got fine points of IEEE-754 wrong for years, so we grew ever more wrapper code to hide platform differences in those. And sometimes platforms had unbearably sloppy implementations, :Last of those I recall was horrid numerical behavior on some platform’s implementations of the gamma functions. So @mdickinson wrote new onrs for Python, and ignore the platform.

Which isn’t a pure win either. Some libms have better behavior for those functions than Python’s have now, but now we’re overriding them. We don’t have enough people with enough time, interest, and expertise to dig into which platforms should be left alone, And when we do hide platform differences, programmers can be baffled by that they get different results from C and Python on the same platform.

And the C/C++/POSIX standards have - as such things always do - bloated and bloated. A sticking point for us is that Windows is very important platform, but Microsoft’s C typically lags the shiny new standards by years. If sinpi (etc) were supplied by all major platform libms, great, it’s easy to write little Python wrappers around them. But it’s too much of a trick to “just wrap” something Microsoft doesn’t supply :wink:.

More generally, Python stays away from trying to be “a pioneer” in numerical work - we aim to follow there, not lead. The people with the right stuff for that tend to contribute to other projects (like numpy, mpmath, sympy …). We attract someone with serious numeric chops maybe once or twice per decade.

The decimal module is an exception. but that was driven by a single person who was extraordinarily dedicated to it.

5 Likes
17

As Serhiy noted above, Python does already have sinpi, hidden as a private function within mathmodule.c. I have not looked, but is it is used, for anything tested, it is already maintained, so the question is whether to expose it. (And perhaps the fairly trivial cospi, …)

1 Like
18

Not so simple. It’s an unadvertised “internal function” for reasons: it ignores 754 endcases, It’s instead used to reduce (not eliminate) numeric error in Python’s implementations of some gamme-related functions. It only handles cases that arise in that context.

It will return 0.0 for any integer argument, though, because fmod() is always exact, and returns exactly 0.0 (or 1.0) if and only if the argument is an exact even (or odd) integer. In both cases it ends up returning sin(pi * +0.0), and I expect all platform libm sin functions return 0.0 for sin(0.0)

The standards require 7:

  • sinpi(x)
  • cospi(x)
  • tanpi(x)
  • asinpi(x)
  • acospi(x)
  • atanpi(x)
  • atan2pi(y, x)

Let’s do the whole job or not bother with any :smile:.

6 Likes
19

Thank you Tim for the reporting the issues with mathmodule’s private sinpy. You may have been the best qualifies to do so. I mentioned it because two core dev mentioned it here and on the gh issue without such details. I will not personally pursue it further.

I realize from the discussion that wrapping sin, etc, and using fmod to detect int values and return exact 0.0 or 1.0 when appropriate will likely serve any personal uses I can currently imagine.

3 Likes
20

I deal with geospatial and survey data all day and have been burned by the trig conversions so many times, it’s especially bad when library functions (ESRI) are inconsistent with usage of radians and degrees meaning you end up with compounding conversion errors lol.

There are definitely some practical uses for sinpi in that sphere though, especially when I’m recursively applying some rotation as I traverse through a point array and don’t want to have to worry about the compounding error eventually breaking free of my if val <= tol check.

3 Likes