IEEE-754 floats have the concept of a “signed zero”; 0.0 has a different bit representation to -0.0. In most cases, -0.0 behaves the same way as 0.0, and it compares equal in arithmetic operations. It becomes more obviously distinct in floating-point operations that involve some form of limiting behaviour. For example, x / 0.0 and x / -0.0 are opposite-signed infinities.

Along with the other IEEE-754 special values, like quiet/signalling NaN and infinities, these sorts of behaviours prevent compilers from making certain mathematical rewrites that would appear to be completely sound in regular arithmetic. $-(a-b)$ in regular arithmetic over the reals can be written as $b-a$, and despite floats having a symmetrical range of positive and negative values (unlike two’s complement integers), this is an invalid transformation in IEEE-754 arithmetic regardless of whether a symmetric rounding mode is in effect. In all rounding modes, the problem occurs at a = b; all finite IEEE-754 floats satisfy x - x = 0.0, therefore a - b and b - a are both 0.0, and negating one to make -0.0 makes it distinct.

In most uses, this is largely a curiosity and has little impact. It can be useful in general when the only information needed from the result of a long calculation is its sign. If the result were to underflow, beyond even the subnormal floats, the resulting sign of the zero would still be able to indicate the correct direction. The signed zero, then, can be thought of as representing the behaviour in the limit.

Where it becomes far more than a curiosity, with major impactful results, is when the signed zero becomes involved in a calculation with a discontinuity at zero. For real numbers, the most obvious of these is atan2(y, x), which is the floating-point version of $\arctan (y/x)$ including the quadrant of the rotation. Since the float -0.0 is then interpreted as ${\lim }_{x\to 0^{-}}x$—the limit as $x$ approaches zero by becoming less negative—there is a natural distinction between atan2(0.0, 0.0) and atan2(0.0, -0.0). Treating these as being calculated within ${\lim }_{y\to 0^{+}}$, it becomes sensible that atan2(0.0, 0.0) would be a zero rotation, while atan2(0.0, -0.0) approximates $\pi$.

I originally licensed the text of all of my blog posts on this website under the Creative Commons Attribution Share-Alike version 4.0 (CC BY-SA 4.0). Concurrent with this blog post, I am changing the licence I distribute my posts on this website to Creative Commons Attribution version 4.0 (CC BY 4.0), dropping the share-alike requirement, including all previous posts. You can still use blog posts dated before the 1st of April, 2024 under the terms of the CC BY-SA 4.0 that they were originally distributed with if you so choose, though to my understanding, CC BY 4.0 is strictly more permissive for the licensee. Code samples remain licensed under the MIT licence.

The longer I have been a professional software developer, the less interest I have had in reciprocity conditions in licensing. I originally used CC BY-SA 4.0 simply because that was what Wikipedia did at the time I needed to choose a licence, and that seemed fine to me. Thinking again now, I see no benefit in enforcing the reciprocity condition.

Python scoping can be deceptive, especially since it can appear fairly simple at first. Generally, name resolution starts from an inner-most scope, and proceeds outwards until some outer scope either contains the name, or we run out of scopes to search. This procedure is familiar to users of many programming languages. Newcomers from other languages are often surprised that various control-flow constructs like for and with don’t open new scopes, so new variables defined with them are still accessible after the block has completed, but this is relatively simple to accept; Python often seems very permissive. This post presents a few corners of Python’s scoping rules that are more unusual than this, some of which continue to trip me up.

Tuples are as immutable as things get in Python. They have special treatment being a core builtin, so unlike user-defined classes, there’s no private attribute that you can mutate if you really want to mess around. However, Python is an interpreted language with no baked-in access control, and (almost) everything really is mutable if you try hard enough.

All Python objects at the C level can be addressed by a pointer of type PyObject *, which (for most regular CPython builds) has the elements

1
2
3
4struct _object {
    Py_ssize_t ob_refcnt;
    PyTypeObject *ob_type;
}

Tuples are a type of “variable-length” object, in the sense that the same C struct definition is used for all sizes. With the head components inlined, this looks like

1
2
3
4
5
6struct _tuple {
    Py_ssize_t ob_refcnt;
    PyTypeObject *ob_type;
    Py_ssize_t ob_size;
    PyObject *ob_item[1];
}

From Python space, tuple.__setitem__ raises a TypeError, but in principle at the C level, there is no such immutability. As of Python 3.11, the tupleobject.h header file even has a comment explaining how tuple can be used as a temporary array, and there is a PyTuple_SetItem function implementation, although this has slightly different reference-counting semantics to normal.

We can use the Python C FFI ctypes, and the CPython implementation detail that id returns the memory address of the Python object to cheekily define our own Python-space tuple_setitem. The PyObject type in ctypes doesn’t expose the internal fields of the object, but instead we can change the items we want by casting the integer output of id into an size_t (assuming, as in many systems that pointers have the same size as size_t), and manually assigning into the ob_item array. To be better citizens of Python—rampant pointer abuse aside—we should also increment the reference count of the object we are inserting into the tuple, and decrement the count for the object we are removing. This leaves us with the Python-space function

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18from ctypes import cast, c_size_t, c_ssize_t, POINTER
from typing import Any

def tuple_setitem(tup: tuple, index: int, new: Any):
    if -len(tup) <= index < 0:
        index = len(tup) + index
    if not 0 <= index < len(tup):
        raise IndexError("tuple index out of range")

    # Decrement the old item's refcount.
    old = tup[index]
    cast(id(old), POINTER(c_ssize_t))[0] -= 1

    # Increment the new item's refcount.
    cast(id(new), POINTER(c_ssize_t))[0] += 1

    # Replace the item.
    cast(id(tup), POINTER(c_size_t))[3 + index] = id(new)

1
2
3
4>>> tup = (1, 2, 3)
>>> tuple_setitem(tup, 0, 7)
>>> tup
(7, 2, 3)

Due to recent problems free CI servers have had with abuse from cryptocurrency miners, Travis CI are scaling back what services are available to open source projects. We used to use Travis for QuTiP, but the new pricing model has meant that we won’t be able to afford it once the old free tier travis-ci.org completely shuts down. We trialled GitHub Actions to build the 28 different wheels we now distribute on PyPI, and then moved all our testing there as well.

Here are few nice features that I’ve been using since we started. These are mostly not unique to GitHub Actions, but all my example configuration code will be for it.

In the second part of this series on series expansions (introduction), we found a way to make a series expansion using any basis of orthogonal polynomials. This is how the standard Fourier series is defined, and then we also used it to make a series expansion out of the Legendre polynomials. Now we’re going to compare how these behave with different numbers of terms.

This is finally the fun part of the seminar! My second-year physics students should generally have already known the previous two articles from first-year courses, so they started here, pretty much. To interact with these series expansions, you should load up my Jupyter notebook on either Colab or myBinder, and try plotting your own functions and seeing how the expansions behave.

In reality we don’t use the infinite form of series expansions, we only use the first few terms to get a good approximation of a function. The notebook is used for investigating how different series expansions work with comparably few terms, and which would be most appropriate in a given situation. The series there are the standard Fourier series, the Fourier–Legendre series, and the Taylor series about $x=0$ (which is formed by a different method to the other two). If you need a reminder of how these series look like, the Taylor and Fourier series forms were given in the introduction, and the Legendre series was derived in part two.

I asked the students to consider a couple of the following questions, and discuss them in groups of around five:

• Which series would you call the “most accurate”? Why? Does it depend on the function?
• What sorts of functions are the different expansions best at approximating? Which are they bad at?
• The Legendre series often seems to “give up” in the middle of some shapes at low orders (e.g. sinusoids). Why is this, and are there any things the series is still useful for?
• The Taylor series almost invariably has the largest pointwise error. Why is this? What is the Taylor series useful for?

In the first part of this series on series expansions (introduction), we defined several concepts from first-year linear algebra that we need to work in a general, abstract setting rather than having to repeat the same derivations over and over again. We are now looking at the topic at the centre of this series; generalising the Fourier series method.

Amazingly, the few definitions I gave in the previous article let us define a whole family of different functional series expansions. Like before, we will consider only continuous, finite functions defined on the interval $[-1,1]$. Let’s say we have a basis of functions that are orthogonal under the inner product

One example of this is the trigonometric functions $\sin (k\pi x)$ and $\cos (k\pi x)$ for all non-negative integers $k$. It’s beyond the scope of undergraduate physics courses to prove that the trigonometric functions span this space, but they do, and you also can verify that the inner-product integral is indeed zero for any unequal pair. We’ll refer to elements of this basis as ${\varphi }_n(x)$, where $n$ is just a unique label.

Now, since the functions span the space of functions we can write down any function $f$ in terms of our basis ${\varphi }_n$ and some scalar coefficients $c_n$ as

and this series is unique. For convenience, we’ll call the series representation $F_f$ to distinguish it from $f$ while we’re still determining the coefficients.

In the introduction to this series of articles, we introduced the idea of a series expansion, and saw the explicit forms of the Taylor and Fourier series. We are building towards a generalised series expansions based on “orthogonal” polynomials, but first we have to define some concepts from linear algebra. This article is roughly at the level of early first-year physics undergradutes.

Students at this level have come across the word “orthogonal” before when talking about Euclidean (“normal”) vectors. Here it means the same thing as “perpendicular” or “at right angles”, at least while you have three or fewer dimensions. Once you have more than that, or you’re dealing with some other type of vector, the definition is a little more abstract.

In all my time as a PhD student I have been involved with teaching undergraduate and postgraduate physics students at Imperial College London (they even gave me a prize in 2018!) as a teaching assistant in classroom-style tutorials and seminars. For my final year, though, I’ve also moved up into helping write the teaching materials, particularly for the differential equations part of the second-year course.

The first differential equations seminar we had was showing off some of the properties and uses of the Legendre polynomials, which the students had just met by solving Legendre’s equation. My part of the seminar was illustrating how any orthogonal basis of functions can be used to make a series expansion, and then getting the students to investigate how different types of functional expansion behave at different orders. If you just want to play with this, load up the Jupyter notebook on Colab or myBinder. The source code for this notebook is available on GitHub. You can see an example plot of these below.

The end result we’ll achieve mathematically is an abstract way of making series expansions. In general, a series expansion approximates a function $f(x)$ by using a (possibly infinite) sequence of terms, where each term is a constant multiplied by some basis function. Different series expansions use different bases and different methods of determining the coefficients.

Perhaps the most familiar example is the Taylor series. The Taylor series $T_f(x)$ that approximates a function $f(x)$ around some point $x_0$ is

This post is the final permalink for the work I’ve done for QuTiP in the Google Summer of Code 2020. My mentors have been Alex Pitchford, Eric Giguère and Nathan Shammah, and QuTiP is under the numFOCUS umbrella.

The main aim of the project was to make Qobj, the primary data type in QuTiP, able to use both sparse and dense representations of matrices and have them interoperate seamlessly. This was a huge undertaking that had far-reaching implications all across the library, but we have now succeeded. There is still plenty of work to be done in additional development documentation and on sanding out the edges to improve the UX, but we are moving towards a public beta of a major version update next year.