OK, I think I figured this out. If I am not mistaken, basically you can use decomposition into a product of reflections, and then improve this to take the signs into account. This gives a $\mathcal{O}(n^3)$ algorithm for multiplication - as good as naïve matrix multiplication.
Let me first describe an algorithm for representing (and multiplying) elements of the group* $\operatorname{O}(n)$. And then I will explain how I think it can be adapted to $\operatorname{Pin}^\pm(n)$.
A presentation of $\operatorname{O}(n)$: We start with the well-known fact that $\operatorname{O}(n)$ is generated by reflections. More precisely, for each unit vector $u \in \mathbb{S}(\mathbb{R}^n)$, let $\sigma_u$ denote the reflection through the hyperplane orthogonal to $u$. (Of course it is more appropriate to parametrize these hyperplanes by the projective space $\mathbb{P}(\mathbb{R}^n)$; but the unit sphere $\mathbb{S}(\mathbb{R}^n)$ is more convenient for concrete algorithmic representation, and will in fact be necessary for the Pin groups --- see below).
Then $\operatorname{O}(n)$ is generated by this set $\{\sigma_u | u \in \mathbb{S}(\mathbb{R}^n)\}$, and presented by the following relations:
- $\forall u \in \mathbb{S}(\mathbb{R}^n),\quad \sigma_u^2 = \operatorname{Id}$;
- $\forall u \in \mathbb{S}(\mathbb{R}^n),\quad \sigma_{-u} = \sigma_u$;
- $\forall u, v \in \mathbb{S}(\mathbb{R}^n)$ s.t. $\langle u, v \rangle = 0$ and $\forall (\alpha, \beta), (\alpha', \beta') \in \mathbb{S}(\mathbb{R}^2)$,
$$\sigma_{\alpha u + \beta v} \sigma_u = \sigma_{(\alpha \alpha' - \beta \beta')u + (\alpha \beta' + \alpha' \beta)v} \sigma_{\alpha' u + \beta' v}.$$
(The last relation is basically about the fact that the product of two distinct reflections is a rotation, but that the decomposition of a rotation into such a product is not unique.)
Canonical form in $\operatorname{O}(n)$: We now define a canonical form for an element $g \in \operatorname{O}(n)$ as a product
$$g = x_1 x_2 \cdots x_n$$
with each $x_i$ lying in the set $R_i$ (or, if you prefer, the smaller set $R^+_i$), which, by definition, comprises $\operatorname{Id}$ and all the reflections $\sigma_u$ with respect to vectors $u$ whose first $i-1$ coordinates all vanish, and whose $i$-th coordinate is nonzero (resp. positive).
It is then easy to see that this canonical form is unique (up to sign if using $R_i$, and completely unique if using $R_i^+$). Indeed, $x_i$ can be computed as the unique (up to sign) element of $R_i$ that sends the vector $f_i = x_{i-1} \cdots x_1 g \cdot e_i$ (guaranteed by induction to be orthogonal to $e_1, \ldots, e_{i-1}$) to $e_i$. Even more explicitly, it is equal to $\operatorname{Id}$ if $e_i = f_i$, and to the reflection about $\frac{e_i - f_i}{\|e_i - f_i\|}$ otherwise.
The multiplication algorithm for $\operatorname{O}(n)$: To explain how to multiply two canonical forms, it suffices to explain how to multiply on the left a canonical form $x_1 x_2 \cdots x_n$ by an arbitrary reflection $\sigma_u =: y_1$. We convert the product $y_1 x_1 \cdots x_n$ to canonical form step by step, in the following way:
$$y_1 x_1 x_2 \cdots x_n \\
= x'_1 y_2 x_2 \cdots x_n \\
= \vdots \\
= x'_1 x'_2 \cdots x'_n y_{n+1},$$
while ensuring that, at each step, $x'_i$ still belongs to $R_i$, and furthermore $y_i$ belongs to $R_i \cup R_{i+1} \cup \ldots \cup R_{n+1}$ (with the obvious convention $R_{n+1} = \{\operatorname{Id}\}$). In the end, we are thus reduced to $y_{n+1} = \operatorname{Id}$ that we can discard, and we are left with $x'_1 x'_2 \cdots x'_n$ which is in canonical form.
More precisely, the $i$-th step (converting $y_i x_i$ to $x'_i y_{i+1}$) works as follows:
- If $y_i$ is already equal to $\operatorname{Id}$, simply set $x_i' = x_i$ and $y_{i+1} = \operatorname{Id}$.
- If $x_i = \operatorname{Id}$, set $x'_i = y_i$ and $y_{i+1} = \operatorname{Id}$.
- If $y_i x_i = \operatorname{Id}$, set $x'_i = \operatorname{Id}$ and $y_{i+1} = \operatorname{Id}$.
- Otherwise, we have $y_i = \sigma_w$ and $x_i = \sigma_u$ for some linearly independent unit vectors $u \in (0, \ldots, 0, \mathbb{R}^*, \mathbb{R}, \ldots, \mathbb{R})$ and $w \in (0, \ldots, 0, \mathbb{R}, \mathbb{R}, \ldots, \mathbb{R})$. Consequently we can express $w$ as $\alpha u + \beta v$, for some unit vector $v \in (0, \ldots, 0, \mathbb{R}, \mathbb{R}, \ldots, \mathbb{R})$ orthogonal to $u$ and for some $(\alpha, \beta) \in \mathbb{S}(\mathbb{R}^2)$. We then apply the relation 3 above, adjusting the pair $(\alpha', \beta')$ so as to make the $i$-th coordinate of $\alpha' u + \beta' v$ vanish. This choice is then unique (up to sign), and one can verify that the $i$-th coordinate of the new $x'_i$ cannot vanish.
A presentation of $\operatorname{Pin}^\pm(n)$: now, in any of these two double covers, each reflection $\sigma_u = \sigma_{-u} \in \operatorname{O}(n)$ has two lifts. One can easily verify that there exists a natural (one-to-one) parametrization of the set of these lifted reflections by $\mathbb{S}(\mathbb{R}^n)$, i.e. we may call the two lifts $\tilde{\sigma}_u$ and $\tilde{\sigma}_{-u}$ (which are now distinct!).
Now if I am not mistaken (this is the part where I am only like 80% confident: please double-check!!), the group $\operatorname{Pin}^\pm(n)$ is presented, over this set $\{\tilde{\sigma}_u | u \in \mathbb{S}(\mathbb{R}^n)\}$, by a very similar set of relations:
- $\forall u \in \mathbb{S}(\mathbb{R}^n),\quad \tilde{\sigma}_u^2 = +\operatorname{Id}$ (for $\operatorname{Pin}^+$) or $-\operatorname{Id}$ (for $\operatorname{Pin}^-$);
- $\forall u \in \mathbb{S}(\mathbb{R}^n),\quad \tilde{\sigma}_{-u} = -\operatorname{Id} \cdot \tilde{\sigma}_u$;
- exactly verbatim as for $\operatorname{O}(n)$ (just add tildas).
Canonical form in $\operatorname{Pin}^\pm(n)$: We can then similarly define a canonical form for an element $g \in \operatorname{Pin}^\pm(n)$ as a product
$$g = x_1 x_2 \cdots x_n,$$
with each $x_i$ lying in the set
$$\tilde{R}_i := \{\pm \operatorname{Id}\} \cup \{\tilde{\sigma}_u | \sigma_u \in R_i\}.$$
The canonical form is then once again unique, up to changing the sign of an even number of factors. Alternatively, to avoid this ambiguity, you may use the smaller set
$$\tilde{R}^+_i := \{\operatorname{Id}\} \cup \{\tilde{\sigma}_u | \sigma_u \in R^+_i\}:$$
you then get a completely unique canonical form, at the expense of adding an extra bit of information to keep track of the global sign: i.e. it then looks like
$$g = (\pm \operatorname{Id}) x_1 x_2 \cdots x_n.$$
The multiplication algorithm for $\operatorname{Pin}^\pm(n)$ is then a completely straightforward generalization of the algorithm for $\operatorname{O}(n)$ (at least assuming that the presentation I have given is correct): just replace $\operatorname{Id}$ by $\pm \operatorname{Id}$ as needed.
* I apologize for the confusion between the big-O notation $\mathcal{O}(n^3)$ and the name of the orthogonal group $\operatorname{O}(n)$ - the Context Club, as Gro-Tsen calls it, must be laughing really hard on this one.
