Lisp + Calculus = Identities

R. Wm. Gosper

Symbolics, Inc.

700 East El Camino Real, Suite 120, Mountain View, California 94040

∑_{k>−∞}\frac{∏_{n≥2} 1 + tan \frac{2k + 1}{(−2)^{n}}π}{(2k + 1)^{2}} e^{iπ \frac{2k+1}{260}} = \frac{2π^{2}}{9} ,  (1)

∑_{k>−∞}\frac{∏_{n≥2} 1 + tan \frac{2k + 1}{(−2)^{n}}π}{(2k + 1)^{2}} e^{iπ \frac{2k+1}{261}} = \frac{2π^{2}}{9} + \frac{4474421120 + 52016651039i}{2^{42} + i} \frac{5π^{2}}{36} ,

(2)

∑_{k>−∞} \frac{1 − cos(π(k + ϕ)t) ∏_{n≥1} 1 + s (sin \frac{k + ϕ}{3^{n}}π)^{2}}{k + ϕ} = 0, −1 < t < 1, (3)

∑_{n≥0} \frac{q^{\binom{n+1}{2}}}{a − cq^{n}} \frac{∏_{k≥0} 1 − q^{k−n} − ap^{k} + cp^{k}q^{k}}{(q − 1)(q^{2} − 1) … (q^{n} − 1)} = \frac{(1 − ∏_{n≥0}1 − ap^{n}) ∏_{n≥1} 1 − q^{n}}{−a} , (4)

∑_{n≥0} \frac{(−)^{n} cos π\sqrt{(n + \frac{1}{2})^{2} + a^{2}}}{(n + \frac{1}{2})(n + \frac{1}{2} − a)(n + \frac{1}{2} + a)} = \frac{π}{2} \frac{sin^{2} πa − 2 sin^{2} \frac{πa}{\sqrt{2}}}{a^{2} cos πa} .  (5)

(Identity (3) was derived assuming ϕ real and |s + 2| < 2, but appears to hold more widely. Identity (4) requires |p|, |q| < 1, and these suffice.)

In fathering Lisp, John McCarthy grandfathered computer algebra, all the while maintaining a clandestine affair with mathematical analysis. Thus, he must have been interested to hear that analysis and computer algebra were shacked up at my place. Sadly, all that issued was a bunch of useless identities—proverbial bicycles for John’s proverbial fish.

Yet can we be sure that fish do not find bicycles entertaining? To keep the dosage small, we have selected five sums of infinite products. (The cos √ in (5) comes from an infinite product.)

This selection of identities is a little misleading, since the first three came from Fourier analysis of fractal curves, while the last two came from a (so far) unrelated series derangement scheme. Perhaps the juxtaposition will help someone find a connection.

1 Rotational Symmetry

Let z(t), −1 ≤ t ≤ 1, be an arc in the complex plane, and z_{m}(t) be the period 2m arrangement of m such arcs around a regular m-gon, so that z_{m}(t + 2) = e^{2iπ/m}z_{m}(t). Suppose further that z_{m} is nice enough to equal its Fourier approximation,

z_{m}(t) = ∑_{j} a_{j}e^{iπjt/m}.

Then the Fourier coefficients are:

a_{j} = \frac{1}{2m} ∫_{−m}^{m} e^{−iπjt/m}z_{m}(t)dt = \frac{1}{2m} ∫_{−1}^{2m−1} e^{−iπjt/m}z_{m}(t)dt

= \frac{1}{2m} ∑_{n=0}^{m−1} ∫_{−1}^{1} e^{−iπj(t+2n)/m}z_{m}(t + 2n)dt

= \frac{1}{2m} ∑_{n=0}^{m−1} e^{−iπ(j−1)2n/m} ∫_{−1}^{1} e^{−iπjt/m}z_{m}(t)dt

= \{ \begin{matrix}\frac{1}{2} ∫_{−1}^{1} e^{−iπjt/m}z_{m}(t)dt, & if m|j − 1; \\ 0, & otherwise. \end{matrix}

Let us further stipulate that the arcs remain congruent as we vary m, and each of the z_{m}(t) is rotated so as to differ from z(t) by only a constant c_{m} := c cot(π/m) = c cot π(k+1/m), for −1 ≤ t ≤ 1. (For fixed side-length, m-gons are proportional to cot(π/m).) Then we can write:

z_{m}(t) = ∑_{k} a_{km+1}e^{iπ(k+1/m)t} = ∑_{k} A(k + 1/m)e^{iπ(k+1/m)t},

A(x) := \frac{1}{2} ∫_{−1}^{1} e^{−iπxt}(z(t) + c cot(πx))dt

=: b(x) + c\frac{cos πx}{πx}.

Thus, if we can compute the function b(x) := ∫_{−1}^{1} e^{−iπxt}z(t)dt/2, we know the Fourier series that arranges the arcs around an arbitrary m-gon. This is actually rather magic, and is not equivalent to taking the mth root of some z_{1}(t), which would distort the arcs, or to taking every mth term of some series free of m. In fact, m can be fractional, which arranges the arcs around a “star”.

2 A Fractal Arc

The arc z(t) I wish to repeat around the m-gon is continous, and joins −i to i, so that c = 1 and z_{m}([−1, 1]) = cot π/m + z([−1, 1]), and the m-gon has side length 2. This will make z_{m} continuous, but we will give it unbounded variation, just to annoy the lawyers.

The arc is to be dyadically self-similar: the segment joining z(−1) = −i to the “midpoint” s := z(0) will be z(t) shifted and scaled so that −i maps to s and i maps to −i. Symmetrically, the segment from s to i will mimic z(t) from i to −i.

Koch’s snowflake curve is the case m = 3, s = 1/\sqrt{3}, or alternatively, m = 6, s = −1/\sqrt{3}, of z_{m}(t) [Man, pp. 42–44]. Cesàro’s squarefilling function is m = 2, s = 1, or m = 4, s = −1 [Man, p. 57, 64-65].

For simplicity, temporarily assume s real, and remove the quotes from “midpoint”. Then (using homogeneous coordinates),

\binom{z(t)}{1} = \{ \begin{matrix}M\binom{z(−1 − 2t)}{1}, & −1 ≤ t ≤ 0; \\ \overline{M}\binom{z(1 − 2t)}{1}, & 0 ≤ t ≤ 1; \end{matrix}

where

M := ( \begin{matrix}\frac{is − 1}{2} & \frac{s − i}{2} \\ 0 & 1 \end{matrix} ) .

Thus, the “balanced negabinary” expansion of t defines z(t) as an infinite product of M and \overline{M}, and this product converges to

( \begin{matrix}0 & z(t) \\ 0 & 1 \end{matrix} )

for max{|s − i|, |s + i|} < 2, sufficing to recursively define z(t), without need of termination conditions, such as z(0) = s or z(1) = i. Furthermore, for rational t, where the negabinary expansion repeats, we can get exact valuations of z:

\binom{z(1/3)}{1} = \overline{M}\binom{z(1/3)}{1} = \binom{\frac{−is−1}{2}z(1/3) + \frac{s+i}{2}}{1},

that is,

z(1/3) = \frac{s + i}{is + 3}.

Similarly,

\begin{aligned}\binom{z(0)}{1} & = M\binom{z(−1)}{1} = MM\binom{z(1)}{1} = MM\overline{M}\binom{z(−1)}{1} \\ & = M\binom{\frac{s^{2}+1}{4}z(−1) + \frac{s^{2}−3}{4}i}{1}, \end{aligned}

so that z(−1) = −i and z(0) = s, as expected. Thus, we have the balanced negabinary expansions

1/3 = (.\overline{1}\overline{1}\overline{1}\overline{1}\overline{1}\overline{1}\overline{1} …)_{−2} = (1.1\overline{1}\overline{1}\overline{1}\overline{1}\overline{1} …)_{−2},
0 = (.11\overline{1}1\overline{1}1\overline{1} …)_{−2} = (.\overline{1}\overline{1}1\overline{1}1\overline{1}1 …)_{−2}

= (1.1\overline{1}1\overline{1}1\overline{1} …)_{−2} = (\overline{1}.\overline{1}1\overline{1}1\overline{1}1 …)_{−2} ⋮

where \overline{1} means the digit −1. So this system has infinitely many representations of every number, none of which terminates, for lack of a 0 digit. (But then, what do they do on the left?). Whoever named this “balanced” must have subcontracted on the Tower of Pisa.

3 The Arc’s Fourier Series

Define ω := e^{−iπx} and consider the sequence

\binom{z(0)}{1}  =  \binom{z(0)}{1}

\frac{ω^{−1/2}M + ω^{1/2}\overline{M}}{2} \binom{z(0)}{1}  = ( \frac{ω^{−1/2}z(−\frac{1}{2}) + ω^{1/2}z(\frac{1}{2})}{2} ; cos(−πx/2) )

\frac{ω^{−1/2}M + ω^{1/2}\overline{M}}{2} \frac{ω^{1/4}M + ω^{−1/4}\overline{M}}{2} \binom{z(0)}{1} =

( \frac{ω^{−3/4}z(−\frac{3}{4}) + ω^{−1/4}z(−\frac{1}{4}) + ω^{1/4}z(\frac{1}{4}) + ω^{3/4}z(\frac{3}{4})}{4} ; cos(−πx/2) cos(πx/4) )

⋮

In the limit, the sum on the right becomes an integral:

( \begin{matrix}0 & \frac{1}{2} ∫_{−1}^{1} ω^{t}z(t)dt \\ 0 & d(x) \end{matrix} )

= \frac{ω^{−1/2}M + ω^{1/2}\overline{M}}{2} \frac{ω^{1/4}M + ω^{−1/4}\overline{M}}{2} \frac{ω^{−1/8}M + ω^{1/8}\overline{M}}{2} ⋯

= \frac{ω^{−1/2}M + \overline{ω^{−1/2}M}}{2} \frac{ω^{1/4}M + \overline{ω^{1/4}M}}{2} \frac{ω^{−1/8}M + \overline{ω^{−1/8}M}}{2} ⋯

= ℜ(ω^{−1/2}M) ℜ(ω^{1/4}M) ℜ(ω^{−1/8}M) ℜ(ω^{1/16}M) …,

where

d(x) := ∏_{j≥1} cos (πx(−2)^{−j}) = ∏_{j≥1} cos \frac{πx}{2^{j}}.

So if we define

W_{n} : = ℜ(ω^{(−2)^{−n}}M)

= ( \begin{matrix}\frac{s}{2} sin \frac{πx}{(−2)^{n}} − \frac{1}{2} cos \frac{πx}{2^{n}} & \frac{b}{2} cos \frac{πx}{2^{n}} − \frac{a}{2} sin \frac{πx}{(−2)^{n}} \\ 0 & cos(πx2^{−n}) \end{matrix} ) ,

a : = 1,
b : = s,

then

( \begin{matrix}0 & b(x) \\ 0 & d(x) \end{matrix} ) = W_{1}W_{2}W_{3} ….

That is, our coveted Fourier coefficient b(x) is just the upper right-hand element of the infinite matrix product of the W_{n}. This product converges for all s and x, and, miraculously, telescopes to a product of scalars. Telescopy in matrices means that you can find a nonzero matrix T_{n} so that

W_{n}T_{n+1} = T_{n},  ∀ relevant n,

and thus,

W_{1}W_{2} … W_{n}T_{n+1} = T_{1},  (Telesc)

so as to provide useful relationships among the elements of the product.

With some luck and much computer algebra, we find a successful column vector

T_{n} := ( \begin{matrix}(−2)^{n−1} \frac{(3a − bs) cos (πx(−2)^{1−n}) + (as − b) sin (πx(−2)^{1−n})}{s^{2} − 3} \\ (−2)^{n−1} sin (πx(−2)^{1−n}) \end{matrix} )

= ( \begin{matrix}−(−2)^{n−1} cos (πx(−2)^{1−n}) \\ (−2)^{n−1} sin (πx(−2)^{1−n}) \end{matrix} ) ,

when a and b are specialized for z(t). If we define b_{n} and d_{n} via

W_{1}W_{2} … W_{n} =: ( \begin{matrix}(−2)^{−n} ∏_{j=1}^{n} (cos \frac{πx}{(−2)^{j}} − s sin \frac{πx}{(−2)^{j}} ) & b_{n}(x) \\ 0 & d_{n}(x) \end{matrix} ) ,

where

\binom{b(x)}{d(x)} = lim_{n→∞} \binom{b_{n}(x)}{d_{n}(x)},

then equation (Telesc) becomes

(−2)^{n}b_{n}(x) sin \frac{πx}{(−2)^{n}} − cos \frac{πx}{2^{n}} ∏_{j=1}^{n} (cos \frac{πx}{(−2)^{j}} − s sin \frac{πx}{(−2)^{j}} )

(−2)^{n} sin(πx(−2)^{−n})d_{n}(x)

= ( \begin{matrix}− cos πx \\ sin πx \end{matrix} ) .

Letting n → ∞, and solving for b(x) and d(x),

b(x) = \frac{∏_{j≥1} (cos \frac{πx}{(−2)^{j}} − s sin \frac{πx}{(−2)^{j}} )}{πx} − \frac{cos πx}{πx}

= \frac{sin πx}{π^{2}x^{2}} ∏_{j≥1} (1 − s tan \frac{πx}{(−2)^{j}} ) − \frac{cos πx}{πx},

d(x) = \frac{sin πx}{πx},

where we divided the infinite product by

\frac{πx}{sin πx} ∏_{j≥1} cos (πx(−2)^{−j}) = \frac{d(x)}{d(x)} = 1.

4 Proof of (1) and (2)

When c = 1, the cos(πx)/πx magically cancels from A(x), and the Fourier identity is, at last,

z_{m}(t) = \frac{sin π/m}{π^{2}} ∑_{k}(−)^{k} \frac{e^{iπ(k+1/m)t}}{(k + 1/m)^{2}} ∏_{n≥1} 1 − s tan \frac{π(k + 1/m)}{(−2)^{n}}.

(If c ≠ 1, we still get (the numerator of) m congruent arcs, but displaced radially, and separated by step discontinuities, by virtue of adding the Fourier series of the “m-gon” which is “swept” by lingering at each vertex for 1/m^{th} of the period, and then jumping to the next.)

So now we can write down mysterious-looking things like

∑_{k}(−)^{k} \frac{e^{iπ(k+1/m)/3}}{(k + 1/m)^{2}} ∏_{n≥1} 1 − s tan \frac{π(k + 1/m)}{(−2)^{n}} = \frac{π^{2} (\frac{s + i}{is + 3} + cot \frac{π}{m})}{sin π/m},

from the value of z(1/3).

The graph of z_{−m}(t) is that of z_{m}(t) traced backward, and inside out, which can be seen by negating k in the bilateral series (figure 1, small):

z_{m}(s, t) = z_{−m}(−s, −t).

When s = tan(π/2m), bulging the sides of an m-gon is the same as indenting the sides of a (properly scaled and rotated) 2m-gon (figure 1, large):

z_{m}(s, t) = \frac{e^{−iπ/2m}}{2} sec \frac{π}{2m}z_{2m}(−s, 2t + 1),  s = tan \frac{π}{2m}.

For this s, the infinite product vanishes in the odd terms of the Fourier series for z_{m}(s, t), and simple trigonometry shows that the even terms coincide with the right-hand side, above.

For rational t, we repeat the procedure we used for z(0), recursing until the inevitable recurrence of some previous value (figure 2). Modern Lisps are ideally suited when s is a (perhaps complex) rational, since arbitrary precision complex rationals are a built-in data type. Likewise, some Lisps now have elaborate hash mechanisms to facilitate loop detection. Thus, a Lisp function, looking at s = 1, quickly ran off z(t) for a few hundred rational t, revealing the interesting contrast between z(1/130) and z(2/261), hence identities (1) and (2). (The (−)^{k} is omitted, since the odd terms vanish, since s = 1 = tan(π/4) = tan(π/2m).)

The right-hand side of (2) can also be written

\frac{13 39841 33217702433766728023 + 2^{7} 3^{2} 5 11 71 28252454442382517i}{2^{86} + 4} π^{2}.

For general s, we will need computer algebra. A Macsyma program (next section) finds that identity (1) is the s = 1 case of

∑_{k}(−)^{k} \frac{e^{iπ(k+1/m)/130}}{(k + 1/m)^{2}} ∏_{n≥1} 1 − s tan \frac{π(k + 1/m)}{(−2)^{n}}

= \frac{π^{2}}{sin π/m} (cot \frac{π}{m} + \frac{128s − i (s^{2} − 1) (s^{2} + 1)^{3}}{2(s^{2} + 1)^{3} + 128} ) .

Lisp+Calculus=Identities 239

Figure 1, small: z_{5}(−1,t)=z_{−5}(1,−t), a pentagon indented by Cesàro’s spacefilled halfsquares. Shown is the sum of Fourier terms (rotating vectors) −32≤k<31, for −5≤t≤5, a full period. This is a low pass filtration of motion, rather than space, so the appearance of cusps merely indicates a momentarily vanishing time derivative of the partial sum.

Large: z_{4}(tan π/8,t)∼z_{8}(− tan π/8,2t),i.e., a bulged square = an indented octagon, if |s|=\sqrt{2}−1. Superposed is a Fourier approximation, partially swept (−1≤t≤108/23) by the sum of 18 rotating vectors, the first few of which are discernable at about “8 o’clock” (t= 108/23). From the center, their indices run 0,−1,1,−2,...,8,−9.

240 Gosper

Figure 2: Another way to approximate Cesàro’s halfsquare filler. The values −z_{3}(1,(7k−2557)/2560) for k∈{0,1,...,730}, along with the origin, connected in a self-contacting, but noncrossing polygon. The values were initially computed as exact complex rationals by a Lisp program analogous to the FOU function. But for graphics, floating point approximation proves much faster, providing n/2 bits of both x and y for n levels of recursion.

The constant term in the denominator is 130 because, for s = 0, the arc degenerates to the line segment z(t) = it, −1 ≤ t ≤ 1.

As you might expect, the generalization of (2) is somewhat worse, even with m = 2:

\frac{4}{π^{2}} ∑_{k>−∞}(−)^{k} \frac{e^{iπ(2k+1)/261}}{(2k + 1)^{2}} ∏_{n≥2} 1 + s tan \frac{2k + 1}{(−2)^{n}}π

= \frac{\begin{aligned}(s − i)^{5}(s + i)^{3} ( & s^{77} + 2is^{76} + 34s^{75} + 76is^{74} + ⋯ \\ & −1123939790088261808331s − 2i\frac{2^{84} − 1}{261}) + 2^{84}s \end{aligned}}{(s − i)^{41}(s + i)^{43} + 2^{84}} .

We can get an alternative Fourier expansion for z(t) by reverting to the original product form of b(k + 1/m) in

z(t) = z_{m}(t) − cot \frac{π}{m}

= ∑_{k} \frac{e^{iπ(k+1/m)t}}{π(k + 1/m)} ∏_{n≥1} (cos \frac{π(k + 1/m)}{(−2)^{n}} − s sin \frac{π(k + 1/m)}{(−2)^{n}}) − cot \frac{π}{m},

and then letting m → ∞. The k = 0 term offsets the cot π/m, leaving

z(t) = it + \frac{s}{3} + ∑_{k≠0} \frac{e^{iπkt}}{πk} ∏_{n≥1} (cos \frac{πk}{(−2)^{n}} − s sin \frac{πk}{(−2)^{n}} ) ,  −1 ≤ t ≤ 1,

= \frac{s}{3} + ∑_{k≠0} \frac{e^{iπkt}}{πk} ((−1)^{k−1} + ∏_{n≥1} cos \frac{πk}{(−2)^{n}} − s sin \frac{πk}{(−2)^{n}} ) , |t| < 1,

using the “sawtooth” Fourier series

t = i ∑_{k≠0}(−)^{k} \frac{e^{iπkt}}{πk} = −2 ∑_{k≥1} \frac{sin(πk(t + 1))}{πk},  −1 < t < 1.

This z(t) coincides with z_{2}(t) for −1 < t < 1, but then discontinuously jumps from i back to −i instead of repeating the arc rotated end for end. Thus the centroid of the arc is the constant term, s/3.

Because Fourier approximations interpolate the midpoint of jump discontinuities, we can say that

∑_{k≠0} \frac{(−)^{k}}{k} ∏_{n≥1} (cos \frac{πk}{(−2)^{n}} − s sin \frac{πk}{(−2)^{n}} ) = −\frac{πs}{3}.

We can also say weird things like

∑_{k≠0} \frac{e^{iπk/3}}{k} ∏_{n≥1} (cos \frac{πk}{(−2)^{n}} − s sin \frac{πk}{(−2)^{n}} ) = π(s + i) (\frac{1}{is + 3} − \frac{1}{3}) ,

again using our value of z(1/3).

5 The Macsyma Function

It would have been clearer to name it Z(T), but, as we well know, all programs are named FOO, in honor of Foorier.

(C50) dispfun(fou,bar)$
(E50) FOU(T, [V]) := (V : IF V = [] THEN [1, 0] ELSE V ,
                                                      1
                              GFACTOR(SUBST(V, 'V, BAR )))
                                                      T
                      V  - '(V)
                       2       2
(E51) BAR  := (BAR  : ----------,
         T        T   '(V)  - V
                          1    1
               BAR  :
                  T
                                            [ i S - 1  S - i ]
                                            [ -------  ----- ]
              (i S - 1) FOU(- 2 T - 1, (V . [    2       2   ]) )
                                            [                ]
                                            [    0       1   ]
                                                               1
IF T < 0 THEN ---------------------------------------------------
                                       2
                                                        S - i
                                                      + -----

                                    [ - i S - 1  S + i ]
                                    [ ---------  ----- ]
      (- i S - 1) FOU(1 - 2 T, (V . [     2        2   ]) )
                                    [                  ]
                                    [     0        1   ]
                                                         1
 ELSE -----------------------------------------------------
                                2
                                                      S + i
                                                    + -----)
                                                        2

(To convince (normally case-insensitive) Macsyma to say i instead of %I, alias(\i,%i). You still have to input \i or %i, but the output is prettier.)

The optional argument v maintains the top row of the cumulative product of the matrices M. BAR is the memoizing hash array. Normally, such array elements redefine themselves automatically on the first call, but in this case, they need to redefine themselves twice: first to set the trap for looping back to the same t, and again to record the answer when either the trap springs or some subanswer was found prerecorded. Springing the trap is just evaluating the difference quotient in the temporary definition of BAR. This effectively simulates infinite recursion by computing the fixed point of the pending matrix transformations. In older Macsymas, the simplifier (but not the evaluator) neglects to perform the part selection when given a subscripted list, for example, '[a,b][1]. So in-

stead of subst(v,'v,bar[t]), one might write eval(bar[t]) in the def-

inition of FOU. However, as Lisp teaches us, indiscriminate evaluation is a poor idea (what if s happened to be bound?), so a cleaner program might

subst(map("=",'[v[1],v[2]],v),bar[t]) instead.

Testing it:

(C52) fou(1);
Time= 122 msec.
(D52)                         i
(C53) fou(0);
Time= 61 msec.
(D53)                         S
(C54) fou(1/130);
Time= 10108 msec.
                  8      6      2
              i (S  + 2 S  - 2 S  + 128 i S - 1)
(D54)       - ----------------------------------
                     2        4      2
                 2 (S  + 5) (S  - 2 S  + 13)

Of course, subsequent evaluations of fou(1/130) would be quicker than the exhibited (first) evaluation of fou(0).

Programming quiz: What would happen if you aborted an evaluation with several BAR traps unsprung? Answer: Something bad, the next time you tried to evaluate one of those booby-trapped entries. If you retried the same top-level call to FOU, you would divide 0 by 0. Much worse, if you reach a half-computed BAR entry with a different v accumulated in FOU, you will simply get a wrong answer, and poison the BAR cache.

Like modern Lisps, Macsyma provides the indispensable UNWIND PRO-TECT form with which a more carefully written BAR can detect an abnormal exit and unset the unsprung trap.

6 Adumbration of (3)

Instead of a dyadic fractal arc from −i via s to i, let Z(t) be the triadic from −i via w − i and i − w to i, for some complex w. For convergence, all three segments must have length < 2, that is, |w−i| < 1. For w = i±1/\sqrt{3}, Z_{3}(t) is the triadic dragon that sweeps the “fudgeflake” [Man, p. 73]. [Man, p. 67] incorrectly claims (modulo typos) that if w = (i ± 1)/2, Z_{2}(t) runs back and forth along the “twindragon river”, which is the skin of one side of the traditional (dyadic) dragon curve. In fact, this skin requires a tetradic “isosceles” construction, rather than a “scalene” triadic one.

If w is pure imaginary, the arc is simply linear, but swept in different fashion depending on |w|. For w = 2i/3, the sweep is linear and Z(t) = it. But if w = i, the central segment vanishes, and Z(t)/i (plotted versus t) becomes the devil’s staircase! [Man, p. 82]. Thus, Z_{m}(t) would continuously trace out an ordinary regular m-gon, but with 0 instantaneous velocity “almost always”.

The recursive definition is

Z(t) := \{ \begin{matrix}−\frac{iw}{2} Z(3t + 2) + \frac{w}{2} − i, & −1 ≤ t ≤ −1/3; \\ (iw + 1)Z(3t), & −1/3 ≤ t ≤ 1/3; \\ −\frac{iw}{2} Z(3t − 2) − \frac{w}{2} + i, & 1/3 ≤ t ≤ 1. \end{matrix}

The coefficients of Z() for each case determine the three corresponding matrices M_{−}, M_{0}, M_{+}, which map [−i, i] onto the three segments. Then, analogous to the dyadic case,

2 ∏_{j≥1} \frac{ω^{−3^{−j}}M_{−} + M_{0} + ω^{3^{−j}}M_{+}}{3}

provides the Fourier coefficient A(x) of e^{iπxt}, when ω := e^{−iπx}. (Since the arc is symmetric under 180^{∘} rotation, we need not worry about traversing subarcs backward, hence the lack of (−3)^{−j}.)

Another telescopy miracle reduces the matrix product to a scalar one, and

Z_{m}(t) = \frac{1}{π} ∑_{k>−∞} \frac{∏_{n≥1} 1 + 2iw(sin \frac{k + 1/m}{3^{n}}π)^{2}}{k + 1/m} e^{iπ(k+1/m)t}

repeats the fractal around an m-gon. Since, for −1 ≤ t ≤ 1, Z_{m}(t) = cot(π/m) + Z(t),

∑_{k>−∞} \frac{1 − e^{iπ(k+ϕ)t} ∏_{n≥1} 1 + s (sin \frac{k + ϕ}{3^{n}}π)^{2}}{k + ϕ} ,  −1 ≤ t ≤ 1,

is independent of ϕ, and an odd function of t, given the arc’s aforementioned symmetry.

Identity (3) simply states that the even part is 0. (For t ≠ 0, the series is conditionally convergent, and needs to be summed sufficiently symmetrically that ∑_{k} 1/(k + ϕ) = π cot πϕ.)

More generally, Z(rational t) will be a rational function of w whose degree depends on the balanced ternary expansion of t. For example,

Z(−1) = −i, Z(0) = 0, Z(1) = i,
Z(\frac{1}{3}) = i − w, Z(\frac{1}{2}) = −i − \frac{4}{w + 2i}, Z(\frac{2}{3}) = i − \frac{w}{2} ,

Z(\frac{1}{5}) = −i − \frac{2}{w + i}, Z(\frac{2}{5}) = i + \frac{4w}{w^{2} − 4}, Z(\frac{3}{5}) = −\frac{1}{w + i},

Z(\frac{1}{7}) = i + \frac{4w}{w^{2} + iw − 2}, Z(\frac{3}{7}) = \frac{w − 2i}{w^{2} + iw − 2},

Z(\frac{1}{9}) = i − 2w − iw^{2}, Z(\frac{5}{9}) = i − w − \frac{i}{2}w^{2},

Z(\frac{1}{11}) = −i\frac{(w − i)(w^{2} − 2iw + 4)}{w^{4} − w^{2} − 2iw + 4}, Z(\frac{3}{13}) = \frac{w − i}{w^{2} − 1},

Z(\frac{1}{15}) = \frac{(w − i)^{2}}{w + 1} , Z(\frac{2}{15}) = −\frac{(w − i)(w − 2i)^{2}}{w^{2} − 4} ,

Z(\frac{1}{18}) = i\frac{(w − 2i)(w − i)^{2}}{w + 2i} , Z(\frac{1}{27}) = (w − i)^{3},

Z(\frac{27}{40}) = −\frac{w − 2i}{w(w − i)^{3} + 2}, Z(\frac{27}{41}) = \frac{1}{w^{3} − iw^{2} − w − i},

etc.

It would be interesting to know what subspace of the rational functions this represents, but since Z(t) is continuous, we know that for arbitrarily close t, these functions must arbitrarily well approximate each other in the disk |w − i| < 1. For example, 27/41 < 2/3 < 27/40, and

Z(\frac{27}{40}) − Z(\frac{2}{3}) = \frac{w}{2} \frac{(w − 2i)(w − i)^{3}}{w(w − i)^{3} + 2}

= \frac{(w − i)^{3}}{4} + \frac{(w − i)^{5}}{4} − i\frac{(w − i)^{6}}{8} − \frac{(w − i)^{7}}{8} + ⋯,

Z(\frac{2}{3}) − Z(\frac{27}{41}) = −\frac{w}{2} \frac{(w − i)^{3}}{w^{3} − iw^{2} − w − i}

= \frac{(w − i)^{3}}{4} + \frac{(w − i)^{5}}{4} + i\frac{(w − i)^{6}}{8} + \frac{(w − i)^{7}}{8} + · · · .

We know that, for |t| ≤ 1/3, there will be a factor of w − i, because of the devil’s staircase property; and, of course, Z(t) − it will contain a factor of 3w − 2i. Also, the degree in w (expansion about ∞) appears to be −α_{3} in the rational factorization |t| = 2^{α_{2}}3^{α_{3}}5^{α_{5}} ….

As w describes a continuous arc from i2/3 to i+1/\sqrt{3}, the arc Z(t) continuously deforms from a one-dimensional line segment to a two-dimensional “triadic dragon”. The functions of w above are the (smooth) trajectories of specific fixed t during the dimension boost. Only a Cantor set of the t wind up on the boundary of the dragon.

This entire discussion of rational functions from rational t can be repeated for the dyadic z(t) underlying identities (1) and (2).

We can also repeat the trick of taking m → ∞ in Z(t) = Z_{m}(t) − cot(π/m) to get a nicer but discontinuous

Z(t) = 2i ∑_{k≥1} \frac{sin πkt}{πk} ((−1)^{k−1} + ∏_{n≥1} 1 + 2iw sin^{2} (\frac{πk}{3^{n}} )) , −1 < t < 1.

Lisp+Calculus=Identities 247

Figure 3, frizzy: The trajectories of Z(k/243),k∈{−243,−242,...,243}, as w sweeps from 2i/3 (line segment =Z(t)) to 1/\sqrt{3}+i (spacefilling triadic dragon =Z(t)), along the arc |w|= 2|w−i|, that is, |Z(−1)−Z(−1/3)|=|Z(−1/3)−Z(1/3)|, the “isosceles” Koch constructions. The gaps show Z(t) at intermediate dimensions. Despite the requirement that the points flow from the segment to completely cover an area, there are many confluences—spacefilling functions visit infinitely many points three times.

Knobby: The filled (triadic dragon) area, this time approximated by Fourier terms 0,−1,1,−2,...,63,−64 of Z_{2}(t), −1≤t≤1. During this half-period, the fundamental swept from “5 o’clock” to “11 o’clock”, and the tip of the −127^{th} harmonic swept the meandering curve from “6 o’clock” to “12 o’clock”. Z_{2}(0)=0 via term k cancelling term −k−1, but t=±1 introduces an alternating sign, effectively doubling each term. These are visible as the rotationally symmetric (t=−1 and t=1) broken lines from the center of the diagram. Each segment is hinged in the middle, and folds into two counter-rotating vectors as t departs either endpoint.

248 Gosper

Figure 4: The regions swept by harmonics 1, 5, 13, 29, and 61 during a full period of Z_{2}(t). Observe that they have one lobe (circle), 5 lobes, 13 lobes, ..., respectively.

7 Nonlocal Derangement

While investigating a so-called q-analogue of sin, I conjectured the generalized addition formula,

sin ∑_{k=1}^{n} a_{k} = ∑_{k=1}^{n} sin a_{k} ∏_{{j=1 \atop j≠k}}^{n} (cos a_{j} + sin a_{j} cot(b_{j} − b_{k})),  (n angles)

for arbitrary n, a_{k} and b_{k}, provided the b_{k} are distinct mod π. (Note the similarity of the product to the first formula for b(x).) The special case a_{k} := π/2 appears in Whittaker and Watson [5]. More generally, with a_{k} := arctan y/x,

∑_{k=1}^{n} ∏_{{j=1 \atop j≠k}}^{n} (x + y cot(b_{j} − b_{k})) = \frac{(x^{2} + y^{2})^{n/2}}{y} sin (n arctan \frac{y}{x});

or, if a_{k} := z + πk/n,

∑_{k=1}^{n} ∏_{{j=1 \atop j≠k}}^{n} (cot(z + \frac{πj}{n} ) + cot(b_{j} − b_{k})) = −2^{n−1} \frac{cos(n (z + \frac{π}{2} ))}{sin(nz)}, n ≥ 1.

In (n angles), replacing a_{k} and b_{k} by ϵa_{k} and ϵb_{k}, expanding at ϵ = 0, and equating first order terms:

∑_{k=1}^{n} a_{k} = ∑_{k=1}^{n} a_{k} ∏_{{j=1 \atop j≠k}}^{n} 1 + \frac{a_{j}}{b_{j} − b_{k}} .  (n terms)

As we shall see, this is actually equivalent to (n angles), so in proving the latter, Vilmos Totik proved both [Ism].

Totik observes that the coefficient of (b_{1} − b_{2})^{−1} is

a_{1}a_{2} (∏_{j=3}^{n} (1 + \frac{a_{j}}{b_{j} − b_{2}} ) − ∏_{j=3}^{n} (1 + \frac{a_{j}}{b_{j} − b_{1}} )) ,

which is O(b_{1} − b_{2}). Therefore, the right side of (n terms), as a function of b_{1}, has no poles, and is therefore constant; likewise for the other b_{k}. To determine that constant, choose the b_{k} arbitrarily large, reducing (n terms) to a tautology.

To recover (n angles), range the sums and products in (n terms) over [0, n] instead of [1, n], then subtract the original (n terms), leaving

a_{0} = a_{0} ∏_{j=1}^{n} (1 + \frac{a_{j}}{b_{j} − b_{0}} ) + ∑_{k=1}^{n} a_{k} \frac{a_{0}}{b_{0} − b_{k}} ∏_{{j=1 \atop j≠k}}^{n} 1 + \frac{a_{j}}{b_{j} − b_{k}} .

Then change variables: a_{0} := 1, a_{k} := α_{k}β_{k}, b_{0} := 0, b_{k} := −β_{k}:

1 − ∏_{j=1}^{n} 1 − α_{j} = ∑_{k=1}^{n} α_{k} ∏_{{j=1 \atop j≠k}}^{n} 1 − \frac{α_{j}β_{j}}{β_{j} − β_{k}} .  (n factors)

Then change variables again: α_{k} := 1−e^{−2ia_{k}}, β_{k} := e^{ib_{k}}, multiply through by

\frac{∏_{j=1}^{n} e^{ia_{j}}}{2i} ,

combine the products in the summand, using

\frac{β_{j}}{β_{j} − β_{k}} = \frac{1}{2} + \frac{1}{2} \frac{β_{j} + β_{k}}{β_{j} − β_{k}} ,

and you are done.

But there ain’t nothin’ new under the sun. [Whi] gives (as an exercise!):

cos ∑_{k=1}^{n} (a_{k} − b_{k}) = (W&W)

∏_{j=1}^{n} \frac{sin(x − b_{j})}{sin(x − a_{j})} − ∑_{k=1}^{n} cot (x − a_{k}) sin (a_{k} − b_{k}) ∏_{{j=1 \atop j≠k}}^{n} \frac{sin(a_{k} − b_{j})}{sin(a_{k} − a_{j})}.

For small a_{k}, b_{k}, the expansion to first order is equivalent to (n factors), under a rational function change of variables. This transformation is invertible, making it likely that (W&W) is merely equivalent to all the rest.

8 Proof of (4) and (5)

Identity (4) is just (n factors), with n → ∞, α_{k} := ap^{k}, and β_{k} := c + q^{k}.

Identity (5) also comes from (n factors), and is actually a bilateral series folded in half by coalescing term n with term −n−1. It is the special case ϕ = 1/2, b = i of

∑_{n>−∞} \frac{(−)^{n} cos π\sqrt{(n + ϕ + \frac{a}{b}) (n + ϕ + ab)}}{(n + ϕ)(n + ϕ − a)}

= \frac{π}{a} \frac{\frac{sin(πa) cos(π(ϕ − a))}{sin(πϕ)} − 2 sin^{2} \frac{πa(b + 1)}{2\sqrt{b}}}{sin(π(ϕ − a))} .

This in turn caps a longish sequence of specializations and renamings of coefficients, with both α_{k} and β_{k} initially general reciprocal quadratics, that is, of the form c_{1}/(k + c_{2})(k + c_{3}). These create rational functions in the bilateral infinite product, which simplify to the cosine via the factorial (Γ) reflection formula.

Identity (5) was sufficiently peculiar to merit a Taylor expansion for small a. Indeed, both sides begin

−\frac{π^{5}}{12} a^{2} − \frac{π^{7}}{40} a^{4} − \frac{π^{9}}{96} a^{6} − \frac{17π^{11}}{4032} a^{8} − \frac{409267π^{13}}{239500800}a^{10} − · · · .

For some analogous formulas from (n terms), and generalizations of the cosines to Bessel functions, see Gosper [1]. George Gasper was first to recognize that identities like (5) are discrete analogs of Sonine’s Bessel integerals [Wat].

References

[Gos] R. Wm. Gosper, Mourad Ismail, and Ruiming Zhang. On some strange summation formulas. (To appear.)

[Ism] Mourad Ismail. Electronic mail, January 1991.

[Man] Benoit Mandelbrot. The Fractal Geometry of Nature, W. H. Freeman and Company, San Francisco, 1983.

[Wat] G. N. Watson. A Treatise on the Theory of Bessel Functions. Cambridge University Press, 1944 (second edition), page 415.

[Whi] E. T. Whittaker and G. N. Watson. Modern Analysis. Cambridge University Press, 1963 (fourth edition), page 140.