12.5 Иррациональные и трансцендентные функции

Common Lisp не содержит тип данных, который точно отображает иррациональные числовые значения. В данном разделе функции описаны так, как если бы результаты были математически точными, но фактически все они возвращают число с плавающей точкой приблизительно равное настоящему математическому значению. В некоторых местах изложены математические тождества, связанные со значениями функций, однако, два математически идентичных выражения могут быть различными по причине ошибок процесса аппроксимации при вычислениях чисел с плавающей точкой.

Когда все аргументы для функции данного раздела являются рациональными и Ъ математический результат также является (математически) рациональным, тогда, если не указано иное, реализация может вернуть или точный результат типа rational или приближенное значение с плавающей точкой одинарной точности. Если все аргументы рациональны, но результат не может быть рациональным, тогда возвращается число с плавающей точкой одинарной точности.

Если все переданные в функцию аргументы принадлежат типу (or rational (complex rational)) и Ъ математический результат является (математически) комплексным числом с рациональными действительной и мнимой частями, тогда, если не указано иное, реализация может вернуть или точный результат типа (or rational (complex rational)) или приближенный типа с плавающей точкой с одинарной точностью single-float (только если мнимая часть равна нулю) или (complex single-float). Если все аргументы типа (or rational (complex rational)), но результат не может быть выражен рациональными или комплексным рациональным числом, тогда он будет принадлежать типу single-float (только если мнимая часть равна нулю) или (complex single-float).

Все функции за исключением expt подчиняются правилам неявного приведения плавающей точки или комплексного числа. Когда, возможно после приведения типов, все аргументы становятся с плавающей точкой или комплексным числом с плавающей точкой, тогда результат будет такого же типа, что и аргумента, если не указано иное. ____________________________

Заметка для реализации: Для понимания работы функций из данного раздела может быть полезной «поваренная книга чисел с плавающей точкой» от Cody и Waite [14].

___________________________________________________________________________________________________________

12.5.1 Экспоненциальные и логарифмические функции

Наряду с обычными одно-аргументными и двух-аргументными экспоненциальными и логарифмическими функциями, sqrt рассматривается как экспоненциальная функция, потому что она возводит число в степень 1/2.

[Функция] exp number

Возвращает е, возведённое в степень number, где е является основанием натурального логарифма.

[Функция] expt base-number power-number

Возвращает base-number, возведённый в степень power-number. Если base-number принадлежит типу rational и power-number — integer, тогда результат вычислений будет принадлежать типу rational, в противном случае результат будет приближенным числом с плавающей точкой.

X3J13 voted in March 1989 to clarify that provisions similar to those of the previous paragraph apply to complex numbers. If the base-number is of type (complex rational) and the power-number is an integer, the calculation will also be exact and the result will be of type (or rational (complex rational)); otherwise a floating-point or complex floating-point approximation may result.

Если power-number равен 0 (ноль целочисленного типа), тогда результат всегда будет значение 1, такого же типа что и base-number, даже если base-number равен нулю (любого типа). То есть:

(expt x 0) ≡ (coerce 1 (type-of x))

Если power-number является нулём любого другого типа данных, тогда результат также равен 1 такого же типа, что и аргументы после применения правил приведения числовых типов, с одним условием: если base-number равен нулю, когда power-number является нецелочисленным нулём, то это является ошибкой.

Реализация expt может использовать различные алгоритмы для случаев с рациональным и с плавающей точкой аргументом power-number. Суть в том, что в большинстве случаев боле аккуратный результат может быть достигнут для рационального power-number. Например, (expt pi 16) и (expt pi 16.0) могут вернуть слегка разные результаты, если в первом случае алгоритм «повторных квадратов», а во втором использование логарифмов.

Результатом expt может быть комплексное число, даже если ни один аргумент не был комплексным. Такой результат получается если base-number отрицательное число и power-number не является целочисленным. Следует отметить, что (expt -8 1/3) не может вернуть -2, хотя и -2 несомненно является одним из кубических корней для -8, но основным корнем является аппроксимированное комплексное число #C(1.0 1.73205).

[Функция] log number &optional base

Функция возвращает логарифм числа number c основанием base, которое по умолчанию равно e (эпсилон, основание для натурального логарифма). Например:

(log 8.0 2) ⇒ 3.0
(log 100.0 10) ⇒ 2.0

В зависимости от реализации, результатом (log 8 2) может быть как 3, так и 3.0. Например:

(log -1.0) ≡ (complex 0.0 (float pi 0.0))

X3J13 voted in January 1989 to specify certain floating-point behavior when minus zero is supported. As a part of that vote it approved a mathematical definition of complex logarithm in terms of real logarithm, absolute value, arc tangent of two real arguments, and the phase function as

Logarithm log |z| + i phase z

This specifies the branch cuts precisely whether minus zero is supported or not; see phase and atan.

[Функция] sqrt number

Функция возвращает квадратный корень числа number. Если number не является комплексным числом и отрицательно, тогда результат будет комплексным числом. Например:

(sqrt 9.0) ⇒ 3.0
(sqrt -9.0) ⇒ #c(0.0 3.0)

В зависимости от реализации результат (sqrt 9) может быть как 3, так и 3.0. Результат (sqrt -9) может быть как #(0 3) или #c(0.0 3.0).

X3J13 voted in January 1989 to specify certain floating-point behavior when minus zero is supported. As a part of that vote it approved a mathematical definition of complex square root in terms of complex logarithm and exponential functions as

Square root e(log z)∕2

This specifies the branch cuts precisely whether minus zero is supported or not; see phase and atan.

[Функция] isqrt integer

Целочисленный квадратный корень: аргумент должен быть неотрицательным целым, и результат является наибольшим целым числом, которое меньше или равно точному положительному квадратному корню аргумента.

(isqrt 9) ⇒ 3
(isqrt 12) ⇒ 3
(isqrt 300) ⇒ 17
(isqrt 325) ⇒ 18

12.5.2 Тригонометрические и связанные с ними функции

Некоторые из функций в данном разделе, такие как abs и signum, несомненно не относятся к тригонометрическим функциям, когда рассматриваются как функции только для действительных чисел. Однако, путь, с помощью которого они расширены для операций на комплексных числах, делает их связь с тригонометрическими функциями явной.

[Функция] abs number

Возвращает абсолютное значение аргумента. Для некомплексных чисел x,

(abs x) ≡ (if (minusp x) (- x) x)

и результат всегда имеет тот же тип, что и аргумент.

Для комплексных чисел z, абсолютное значение может быть вычислено как

(sqrt (+ (expt (realpart z) 2) (expt (imagpart z) 2)))

__________________________________________________________________________

Заметка для реализации: Аккуратные разработчики не будут напрямую использовать эту формулу для всех комплексных чисел. Очень большие и очень маленькие части комплексных чисел будут обрабатываться специализированно для избежания выходов за верхние и нижние границы значений.

___________________________________________________________________________________________________________

Например:

(abs #c(3.0 -4.0)) ⇒ 5.0

Результатом (abs #(3 4)) может быть или 5 или 5.0, в зависимости от реализации.

[Функция] phase number

Аргументом (так называется функция phase) числа z (arg(z) ≡(phase z)) является угол φ (в радианах) радиус-вектора точки, соответствующей комплексному числу z.

X3J13 voted in January 1989 to specify certain floating-point behavior when minus zero is supported; phase is still defined in terms of atan as above, but thanks to a change in atan the range of phase becomes − π inclusive to π inclusive. The value − π results from an argument whose real part is negative and whose imaginary part is minus zero. The phase function therefore has a branch cut along the negative real axis. The phase of + 0 + 0i is + 0, of + 0 − 0i is − 0, of − 0 + 0i is + π, and of − 0 − 0i is − π.

Если аргумент является комплексным числом с частями из чисел с плавающей точкой, то результатом является число с плавающей точкой того же типа. Если аргумент является числом с плавающей точкой, то результатом является число с плавающей точкой того же типа. Если аргумент является комплексным числом с частями из рациональных чисел, то результатом является число с плавающей точкой одинарного типа.

[Функция] signum number

По определению,

(signum x) ≡ (if (zerop x) x (/ x (abs x)))

Для рационального числа, signum будет возвращать один из вариантов: -1, 0 или 1, в соответствии с тем, является ли число отрицательным, нулём или положительным. Для числа с плавающей точкой, результатом будет число с плавающей точкой того же типа и значением: − 1, 0 или 1. Для комплексного числа z, (signum z) является комплексным числом с таким же аргументом (phase), но с единичным модулем. Но если z равен комплексному нулю, результатом является само число z. Например:

(signum 0) ⇒ 0
(signum -3.7L5) ⇒ -1.0L0
(signum 4/5) ⇒ 1
(signum #C(7.5 10.0)) ⇒ #C(0.6 0.8)
(signum #C(0.0 -14.7)) ⇒ #C(0.0 -1.0)

Для некомплексных рациональных чисел signum является рациональной функцией, но для комплексных чисел может быть иррациональной.

[Функция] sin radians
[Функция] cos radians
[Функция] tan radians

sin возвращает синус аргумента, cos — косинус, tan — тангенс. Аргумент рассматривается как угол в радианах. Аргумент может быть комплексным числом.

[Функция] cis radians

Функция вычисляет ei⋅radians. Имя cis обозначает «cos + i sin», потому что ei𝜃 = cos 𝜃 + i sin 𝜃. Аргумент рассматривается как угол в радианах и может быть любым некомплексным числом. Результатом является комплексное число, у которого действительная часть это косинус аргумента, и мнимая — синус аргумента. Другими словами, результат это комплексное число, у которого аргумент (фаза) равна (mod 2∕pi) и модуль равен единице. _________________________________________________________________

Заметка для реализации: Чаще всего дешевле вычислить синус и косинус угла вместе, чем выполнять два раздельных вычисления.

___________________________________________________________________________________________________________

[Функция] asin number
[Функция] acos number

asin возвращаем арксинус аргумента, и acos — арккосинус аргумента. Результат будет в радианах. Аргумент может быть комплексным числом.

Функции арксинус и арккосинус могут быть математически определены для аргумента z следующим образом:

Арксинус −i log ( ) iz + √1--−-z2
Арккосинус −i log ( √ ------) z + i 1 − z2

Следует отметить, что результат asin или acos может быть комплексным числом, даже если аргумент не являлся комплексным. Такое происходит, когда абсолютное значение аргумента превышает 1.

Kahan [25] suggests for acos the defining formula

Arc cosine 2 log ( ∘ 1+z- ∘ -1−-z) 2 + i 2 i

or even the much simpler (π∕2) − arcsin z. Both equations are mathematically equivalent to the formula shown above.
__________________________________________________________________________

Заметка для реализации: These formulae are mathematically correct, assuming completely accurate computation. They may be terrible methods for floating-point computation. Implementors should consult a good text on numerical analysis. The formulae given above are not necessarily the simplest ones for real-valued computations, either; they are chosen to define the branch cuts in desirable ways for the complex case.

___________________________________________________________________________________________________________

[Функция] atan y &optional x

Функция вычисляет арктангенс и возвращает результат в радианах.

Ни один из двух аргументов y и x не может быть комплексным. Знаки y и x используются для определения квадранта. x может быть нулём при условии, что y не ноль. Значение atan лежит между − pi (невключительно) и pi (включительно). Следующая таблица описывает различные специальные случаи.

Условие
Декартово местоположение Промежуток результата
y = +0 x > 0 Точно выше положительной оси x + 0
y > 0  x > 0 Квадрант I + 0 < result < π∕2
y > 0  x = ±0Положительная ось y π∕2
y > 0  x < 0 Квадрант II π∕2 < result < π
y = +0 x < 0 Точно ниже негативной оси x π
y = −0 x < 0 Точно выше отрицательной оси x π
y < 0  x < 0 Квадрант III − π < result < −π∕2
y < 0  x = ±0Отрицательная ось y − π∕2
y < 0  x > 0 Квадрант IV − π∕2 < result < −0
y = −0 x > 0 Точно ниже положительной оси x − 0
y = +0 x = +0Рядом с центром + 0
y = −0 x = +0Рядом с центром − 0
y = +0 x = −0Рядом с центром π
y = −0 x = −0Рядом с центром − π
 

Следует отметить, что случай y = 0,x = 0 при отсутствии минус ноля является ошибкой, но четыре случая y = ±0,x = ±0 определены при условии существования минус ноля.

Если не указывать x, y может быть комплексным. Результатом функции является арктангенс y, который может быть определён следующей формулой:

Арктангенс log(1+iy)−log(1−iy) 2i

Для некомплексного аргумента y, результат будет некомплексным и будет лежать между − π∕2 и π∕2 (оба невключительно). _______________

Заметка для реализации: This formula is mathematically correct, assuming completely accurate computation. It may be a terrible method for floating-point computation. Implementors should consult a good text on numerical analysis. The formula given above is not necessarily the simplest one for real-valued computations, either; it is chosen to define the branch cuts in desirable ways for the complex case.

__________________________________________________________________________

[Константа] pi

Данная глобальная переменная содержит значение наиболее приближенное к pi в длинном формате числа с плавающей точкой. Например:

(defun sind (x)     ;Аргумент в градусах
  (sin (* x (/ (float pi x) 180))))

Приближение к pi с другими точностями может быть выполнено с помощью (float pi x), где x является числом с плавающей точкой необходимой точности, или с помощью (coerce pi type), где type является именем необходимого типа, как например short-float.

[Функция] sinh number
[Функция] cosh number
[Функция] tanh number
[Функция] asinh number
[Функция] acosh number
[Функция] atanh number

Данные функции вычисляют гиперболический синус, косинус, тангенс, арксинус, арккосинус и арктангенс, которые для аргумента z математически определены следующим образом:

Гиперболический синус  (e z −e−z)∕2
Гиперболический косинус  (e z + e−z)∕2
Гиперболический тангенс  (e z − e−z)∕(ez + e−z)
Гиперболический арксинус  log( √ ------) z + 1 + z2
Гиперболический арккосинус  log( ∘ ---------------) z + (z + 1 ) (z − 1)∕(z + 1)
Гиперболический арктангенс  log( ∘ ---------) (1 + z) 1∕(1 − z2)

Следует отметить, что результат acosh может быть комплексным, даже если аргумент комплексным числом не был. Это происходит если аргумент меньше 1. Также результат atanh может быть комплексным, даже если аргумент не был комплексным. Это происходит если абсолютное значение аргумента было больше 1. _______________________________________________

Заметка для реализации: Эти формулы математически корректны, предполагая совершенно точные вычисления. Но они могут быть неудобны для вычислений чисел с плавающей точкой. Разработчики реализация должны ознакомится с хорошей книгой по численному анализу. Вышеприведённые формы не обязательно являются самыми простыми для вещественных вычислений. Они выбраны для определения точек ветвлений для случаев использования комплексных чисел.

__________________________________________________________________________

12.5.3 Точки ветвления, главные значения и краевые условия на комплексной плоскости

Многие иррациональные и трансцендентные функции на комплексной плоскости многозначны. Например, в общем случае, логарифмическая функция имеет бесконечное количество комплексных значений. В каждом таком случае для возврата должно быть выбрано одно главное значение.

Common Lisp определяет точки ветвления, главные значения и краевые условия для комплексных функции в соотвествие с предложениями для комплексных функций в APL [36]. Содержимое раздела по большей части скопировано из этих предложений.

Indeed, X3J13 voted in January 1989 to alter the direction of continuity for the branch cuts of atan, and also to address the treatment of branch cuts in implementations that have a distinct floating-point minus zero.

The treatment of minus zero centers in two-argument atan. If there is no minus zero, then the branch cut runs just below the negative real axis as before, and the range of two-argument atan is (−π,π]. If there is a minus zero, however, then the branch cut runs precisely on the negative real axis, skittering between pairs of numbers of the form −x ± 0i, and the range of two-argument atan is [−π,π].

The treatment of minus zero by all other irrational and transcendental functions is then specified by defining those functions in terms of two-argument atan. First, phase is defined in terms of two-argument atan, and complex abs in terms of real sqrt; then complex log is defined in terms of phase, abs, and real log; then complex sqrt in terms of complex log; and finally all others are defined in terms of these.

Kahan [25] treats these matters in some detail and also suggests specific algorithms for implementing irrational and transcendental functions in IEEE standard floating-point arithmetic [23].

Remarks in the first edition about the direction of the continuity of branch cuts continue to hold in the absence of minus zero and may be ignored if minus zero is supported; since all branch cuts happen to run along the principal axes, they run between plus zero and minus zero, and so each sort of zero is associated with the obvious quadrant.

With these definitions, the following useful identities are obeyed throughout the applicable portion of the complex domain, even on the branch cuts:

sin iz = i sinh z   sinh iz = i sin z arctan iz = i arctanh z
cos i z = cosh z   cosh i z = cos z arcsinh i z = i arcsin z
tan iz = i tanh z  arcsin iz = i arcsinh z arctanh iz = i arctan z

I thought it would be useful to provide some graphs illustrating the behavior of the irrational and transcendental functions in the complex plane. It also provides an opportunity to show off the Common Lisp code that was used to generate them.

Imagine the complex plane to be decorated as follows. The real and imaginary axes are painted with thick lines. Parallels from the axes on both sides at distances of 1, 2, and 3 are painted with thin lines; these parallels are doubly infinite lines, as are the axes. Four annuli (rings) are painted in gradated shades of gray. Ring 1, the inner ring, consists of points whose radial distances from the origin lie in the range [1∕4, 1∕2]; ring 2 is in the radial range [3∕4, 1]; ring 3, in the range [π∕2, 2]; and ring 4, in the range [3,π]. Ring j is divided into 2j+1 equal sectors, with each sector painted a different shade of gray, darkening as one proceeds counterclockwise from the positive real axis.

We can illustrate the behavior of a numerical function f by considering how it maps the complex plane to itself. More specifically, consider each point z of the decorated plane. We decorate a new plane by coloring the point f(z) with the same color that point z had in the original decorated plane. In other words, the newly decorated plane illustrates how the f maps the axes, other horizontal and vertical lines, and annuli.

In each figure we will show only a fragment of the complex plane, with the real axis horizontal in the usual manner ( −∞ to the left, + ∞ to the right) and the imaginary axis vertical ( −∞i below, + ∞i above). Each fragment shows a region containing points whose real and imaginary parts are in the range [−4.1, 4.1]. The axes of the new plane are shown as very thin lines, with large tick marks at integer coordinates and somewhat smaller tick marks at multiples of π∕2.

Figure 12.1 shows the result of plotting the identity function (quite literally); the graph exhibits the decoration of the original plane.

Figures 12.2 through 12.20 show the graphs for the functions sqrt, exp, log, sin, asin, cos, acos, tan, atan, sinh, asinh, cosh, acosh, tanh, and atanh, and as a bonus, the graphs for the functions √ ------ 1 − z2, √ ------ 1 + z2, (z − 1)∕(z + 1), and (1 + z)∕(1 −z). All of these are related to the trigonometric functions in various ways. For example, if f(z) = (z − 1)∕(z + 1), then tanh z = f(e2z), and if g(z) = √ -----2 1 − z, then cos z = g(sin z). It is instructive to examine the graph for √ ------ 1 − z2 and try to visualize how it transforms the graph for sin into the graph for cos.

Each figure is accompanied by a commentary on what maps to what and other interesting features. None of this material is terribly new; much of it may be found in any good textbook on complex analysis. I believe that the particular form in which the graphs are presented is novel, as well as the fact that the graphs have been generated as PostScript [1] code by Common Lisp code. This PostScript code was then fed directly to the typesetting equipment that set the pages for this book. Samples of this PostScript code follow the figures themselves, after which the code for the entire program is presented.

In the commentaries that accompany the figures I sometimes speak of mapping the points ±∞ or ±∞i. When I say that function f maps + ∞ to a certain point z, I mean that

z = lim x→+∞f(x + 0 i)

Similarly, when I say that f maps −∞i to z, I mean that
z = lim y→−∞f(0 + yi)

In other words, I am considering a limit as one travels out along one of the main axes. I also speak in a similar manner of mapping to one of these infinities.

Изображение 12.1: Initial Decoration of the Complex Plane (Identity Function)

PIC
This figure was produced in exactly the same manner as succeeding figures, simply by plotting the function identity instead of a numerical function. Thus the first of these figures was produced by the last function of the first edition. I knew it would come in handy someday!


Изображение 12.2: Illustration of the Range of the Square Root Function

PIC
The sqrt function maps the complex plane into the right half of the plane by slitting it along the negative real axis and then sweeping it around as if half-closing a folding fan. The fan also shrinks, as if it were made of cotton and had gotten wetter at the periphery than at the center. The positive real axis is mapped onto itself. The negative real axis is mapped onto the positive imaginary axis (but if minus zero is supported, then −x + 0i is mapped onto the positive imaginary axis and −x − 0i onto the negative imaginary axis, assuming x > 0). The positive imaginary axis is mapped onto the northeast diagonal, and the negative imaginary axis onto the southeast diagonal. More generally, lines are mapped to rectangular hyperbolas (or fragments thereof) centered on the origin; lines through the origin are mapped to degenerate hyperbolas (perpendicular lines through the origin).


Изображение 12.3: Illustration of the Range of the Exponential Function

PIC
The exp function maps horizontal lines to radii and maps vertical lines to circles centered at the origin. The origin is mapped to 1. (It is instructive to compare this graph with those of other functions that map the origin to 1, for example (1 + z)∕(1 −z), cosz, and √1--−-z2-.) The entire real axis is mapped to the positive real axis, with −∞ mapping to the origin and + ∞ to itself. The imaginary axis is mapped to the unit circle with infinite multiplicity (period 2π); therefore the mapping of the imaginary infinities ±∞i is indeterminate. It follows that the entire left half-plane is mapped to the interior of the unit circle, and the right half-plane is mapped to the exterior of the unit circle. A line at any angle other than horizontal or vertical is mapped to a logarithmic spiral (but this is not illustrated here).


Изображение 12.4: Illustration of the Range of the Natural Logarithm Function

PIC
The log function, which is the inverse of exp, naturally maps radial lines to horizontal lines and circles centered at the origin to vertical lines. The interior of the unit circle is thus mapped to the entire left half-plane, and the exterior of the unit circle is mapped to the right half-plane. The positive real axis is mapped to the entire real axis, and the negative real axis to a horizontal line of height π. The positive and negative imaginary axes are mapped to horizontal lines of height ± π∕2. The origin is mapped to −∞.


Изображение 12.5: Illustration of the Range of the Function (z−1)∕(z+1)

PIC
A line is a degenerate circle with infinite radius; when I say “circles” here I also mean lines. Then (z − 1)∕(z + 1) maps circles into circles. All circles through − 1 become lines; all lines become circles through 1. The real axis is mapped onto itself: 1 to the origin, the origin to − 1, − 1 to infinity, and infinity to 1. The imaginary axis becomes the unit circle; i is mapped to itself, as is −i. Thus the entire right half-plane is mapped to the interior of the unit circle, the unit circle interior to the left half-plane, the left half-plane to the unit circle exterior, and the unit circle exterior to the right half-plane. Imagine the complex plane to be a vast sea. The Colossus of Rhodes straddles the origin, its left foot on i and its right foot on −i. It bends down and briefly paddles water between its legs so furiously that the water directly beneath is pushed out into the entire area behind it; much that was behind swirls forward to either side; and all that was before is sucked in to lie between its feet.


Изображение 12.6: Illustration of the Range of the Function (1+z)∕(1−z)

PIC
The function h(z) = (1 + z)∕(1 −z) is the inverse of f(z) = (z − 1)∕(z + 1); that is, h(f(z)) = f(h(z)) = z. At first glance, the graph of h appears to be that of f flipped left-to-right, or perhaps reflected in the origin, but careful consideration of the shaded annuli reveals that this is not so; something more subtle is going on. Note that f(f(z)) = h(h(z)) = g(z) = −1∕z. The functions f, g, h, and the identity function thus form a group under composition, isomorphic to the group of the cyclic permutations of the points − 1, 0, 1, and ∞, as indeed these functions accomplish the four possible cyclic permutations on those points. This function group is a subset of the group of bilinear transformations (az + b)∕(cz + d), all of which are conformal (angle-preserving) and map circles onto circles. Now, doesn’t that tangle of circles through − 1 look like something the cat got into?


Изображение 12.7: Illustration of the Range of the Sine Function

PIC
We are used to seeing sin looking like a wiggly ocean wave, graphed vertically as a function of the real axis only. Here is a different view. The entire real axis is mapped to the segment [−1,1] of the real axis with infinite multiplicity (period 2π). The imaginary axis is mapped to itself as if by sinh considered as a real function. The origin is mapped to itself. Horizontal lines are mapped to ellipses with foci at ± 1 (note that two horizontal lines equidistant from the real axis will map onto the same ellipse). Vertical lines are mapped to hyperbolas with the same foci. There is a curious accident: the ellipse for horizontal lines at distance ± 1 from the real axis appears to intercept the real axis at ± π∕2 ≈±1.57… but this is not so; the intercepts are actually at ± (e + 1∕e)∕2 ≈±1.54….


Изображение 12.8: Illustration of the Range of the Arc Sine Function

PIC
Just as sin grabs horizontal lines and bends them into elliptical loops around the origin, so its inverse asin takes annuli and yanks them more or less horizontally straight. Because sine is not injective, its inverse as a function cannot be surjective. This is just a highfalutin way of saying that the range of the asin function doesn’t cover the entire plane but only a strip π wide; arc sine as a one-to-many relation would cover the plane with an infinite number of copies of this strip side by side, looking for all the world like the tail of a peacock with an infinite number of feathers. The imaginary axis is mapped to itself as if by asinh considered as a real function. The real axis is mapped to a bent path, turning corners at ± π∕2 (the points to which ± 1 are mapped); + ∞ is mapped to π∕2 −∞i, and −∞ to − π∕2 + ∞i.


Изображение 12.9: Illustration of the Range of the Cosine Function

PIC
We are used to seeing cos looking exactly like sin, a wiggly ocean wave, only displaced. Indeed the complex mapping of cos is also similar to that of sin, with horizontal and vertical lines mapping to the same ellipses and hyperbolas with foci at ± 1, although mapping to them in a different manner, to be sure. The entire real axis is again mapped to the segment [−1,1] of the real axis, but each half of the imaginary axis is mapped to the real axis to the right of 1 (as if by cosh considered as a real function). Therefore ±∞i both map to + ∞. The origin is mapped to 1. Whereas sin is an odd function, cos is an even function; as a result two points in each annulus, one the negative of the other, are mapped to the same shaded point in this graph; the shading shown here is taken from points in the original upper half-plane.


Изображение 12.10: Illustration of the Range of the Arc Cosine Function

PIC
The graph of acos is very much like that of asin. One might think that our nervous peacock has shuffled half a step to the right, but the shading on the annuli shows that we have instead caught the bird exactly in mid-flight while doing a cartwheel. This is easily understood if we recall that arccosz = (π∕2) − arcsinz; negating arcsinz rotates it upside down, and adding the result to π∕2 translates it π∕2 to the right. The imaginary axis is mapped upside down to the vertical line at π∕2. The point + 1 is mapped to the origin, and − 1 to π. The image of the real axis is again cranky; + ∞ is mapped to + ∞i, and −∞ to π −∞i.


Изображение 12.11: Illustration of the Range of the Tangent Function

PIC
The usual graph of tan as a real function looks like an infinite chorus line of disco dancers, left hands pointed skyward and right hands to the floor. The tan function is the quotient of sin and cos but it doesn’t much look like either except for having period 2π. This goes for the complex plane as well, although the swoopy loops produced from the annulus between π∕2 and 2 look vaguely like those from the graph of sin inside out. The real axis is mapped onto itself with infinite multiplicity (period 2π). The imaginary axis is mapped backwards onto [−i,i]: + ∞i is mapped to −i and −∞i to + i. Horizontal lines below or above the real axis become circles surrounding + i or −i, respectively. Vertical lines become circular arcs from + i to −i; two vertical lines separated by (2k + 1)π for integer k together become a complete circle. It seems that two arcs shown hit the real axis at ± π∕2 = ±1.57… but that is a coincidence; they really hit the axis at ± tan1 = 1.55….


Изображение 12.12: Illustration of the Range of the Arc Tangent Function

PIC
All I can say is that this peacock is a horse of another color. At first glance, the axes seem to map in the same way as for asin and acos, but look again: this time it’s the imaginary axis doing weird things. All infinities map multiply to the points (2k + 1)π∕2; within the strip of principal values we may say that the real axis is mapped to the interval [−π∕2,+π∕2] and therefore −∞ is mapped to − π∕2 and + ∞ to + π∕2. The point + i is mapped to + ∞i, and −i to −∞i, and so the imaginary axis is mapped into three pieces: the segment [−∞i,−i] is mapped to [π∕2,π∕2 −∞i]; the segment [−i,i] is mapped to the imaginary axis [−∞i,+∞i]; and the segment [+i,+∞i] is mapped to [−π∕2 + ∞i,−π∕2].


Изображение 12.13: Illustration of the Range of the Hyperbolic Sine Function

PIC
It would seem that the graph of sinh is merely that of sin rotated 90 degrees. If that were so, then we would have sinhz = isinz. Careful inspection of the shading, however, reveals that this is not quite the case; in both graphs the lightest and darkest shades, which initially are adjacent to the positive real axis, remain adjacent to the positive real axis in both cases. To derive the graph of sinh from sin we must therefore first rotate the complex plane by − 90 degrees, then apply sin, then rotate the result by 90 degrees. In other words, sinhz = isin(−i)z; consistently replacing z with iz in this formula yields the familiar identity sinhiz = isinz.


Изображение 12.14: Illustration of the Range of the Hyperbolic Arc Sine Function

PIC
The peacock sleeps. Because arcsinh iz = iarcsinz, the graph of asinh is related to that of asin by pre- and post-rotations of the complex plane in the same way as for sinh and sin.


Изображение 12.15: Illustration of the Range of the Hyperbolic Cosine Function

PIC
The graph of cosh does not look like that of cos rotated 90 degrees; instead it looks like that of cos unrotated. That is because coshiz is not equal to icosz; rather, coshiz = cosz. Interpreted, that means that the shading is pre-rotated but there is no post-rotation.


Изображение 12.16: Illustration of the Range of the Hyperbolic Arc Cosine Function

PIC
Hmm—I’d rather not say what happened to this peacock. This feather looks a bit mangled. Actually it is all right—the principal value for acosh is so chosen that its graph does not look simply like a rotated version of the graph of acos, but if all values were shown, the two graphs would fill the plane in repeating patterns related by a rotation.


Изображение 12.17: Illustration of the Range of the Hyperbolic Tangent Function

PIC
The diagram for tanh is simply that of tan turned on its ear: itanz = tanhiz. The imaginary axis is mapped onto itself with infinite multiplicity (period 2π), and the real axis is mapped onto the segment [−1,+1]: + ∞ is mapped to + 1, and −∞ to − 1. Vertical lines to the left or right of the real axis are mapped to circles surrounding − 1 or 1, respectively. Horizontal lines are mapped to circular arcs anchored at − 1 and + 1; two horizontal lines separated by a distance (2k + 1)π for integer k are together mapped into a complete circle. How do we know these really are circles? Well, tanhz = ((exp2z) − 1)∕((exp2z) + 1), which is the composition of the bilinear transform (z − 1)∕(z + 1), the exponential expz, and the magnification 2z. Magnification maps lines to lines of the same slope; the exponential maps horizontal lines to circles and vertical lines to radial lines; and a bilinear transform maps generalized circles (including lines) to generalized circles. Q.E.D.


Изображение 12.18: Illustration of the Range of the Hyperbolic Arc Tangent Function

PIC
A sleeping peacock of another color: arctanh iz = iarctanz.


Изображение 12.19: Illustration of the Range of the Function ------ √ 1 − z2

PIC
Here is a curious graph indeed for so simple a function! The origin is mapped to 1. The real axis segment [0,1] is mapped backwards (and non-linearly) into itself; the segment [1,+∞] is mapped non-linearly onto the positive imaginary axis. The negative real axis is mapped to the same points as the positive real axis. Both halves of the imaginary axis are mapped into [1,+∞] on the real axis. Horizontal lines become vaguely vertical, and vertical lines become vaguely horizontal. Circles centered at the origin are transformed into Cassinian (half-)ovals; the unit circle is mapped to a (half-)lemniscate of Bernoulli. The outermost annulus appears to have its inner edge at π on the real axis and its outer edge at 3 on the imaginary axis, but this is another accident; the intercept on the real axis, for example, is not really at π ≈ 3.14… but at ∘ --------- 1 − (3i)2 = √ --- 10 ≈ 3.16….


Изображение 12.20: Illustration of the Range of the Function ------ √ 1 + z2

PIC
The graph of q(z) = √ ----2- 1+ z looks like that of p(z) = √ ----2- 1− z except for the shading. You might not expect p and q to be related in the same way that cos and cosh are, but after a little reflection (or perhaps I should say, after turning it around in one’s mind) one can see that q(iz) = p(z). This formula is indeed of exactly the same form as coshiz = cosz. The function √ -----2 1 + z maps both halves of the real axis into [1,+∞] on the real axis. The segments [0,i] and [0,−i] of the imaginary axis are each mapped backwards onto segment [0,1] of the real axis; [i,+∞i] and [−, −∞i] are each mapped onto the positive imaginary axis (but if minus zero is supported then opposite sides of the imaginary axis map to opposite halves of the imaginary axis—for example, q(+0 + 2i) = √ -- 5i but q(−0 + 2i) = −√ -- 5i).

Here is a sample of the PostScript code that generated figure 12.1, showing the initial scaling, translation, and clipping parameters; the code for one sector of the innermost annulus; and the code for the negative imaginary axis. Comment lines indicate how path or boundary segments were generated separately and then spliced (in order to allow for the places that a singularity might lurk, in which case the generating code can “inch up” to the problematical argument value).

The size of the entire PostScript file for the identity function was about 68 kilobytes (2757 lines, including comments). The smallest files were the plots for atan and atanh, about 65 kilobytes apiece; the largest were the plots for sin, cos, sinh, and cosh, about 138 kilobytes apiece.


% PostScript file for plot of function IDENTITY
% Plot is to fit in a region 4.666666666666667 inches square
%  showing axes extending 4.1 units from the origin.

40.97560975609756 40.97560975609756 scale
4.1 4.1 translate
newpath
  -4.1 -4.1 moveto
  4.1 -4.1 lineto
  4.1 4.1 lineto
  -4.1 4.1 lineto
  closepath
clip
% Moby grid for function IDENTITY
% Annulus 0.25 0.5 4 0.97 0.45
% Sector from 4.7124 to 6.2832 (quadrant 3)
newpath
  0.0 -0.25 moveto
  0.0 -0.375 lineto
  %middle radial
  0.0 -0.375 lineto
  0.0 -0.5 lineto
  %end radial
  0.0 -0.5 lineto
  0.092 -0.4915 lineto
  0.1843 -0.4648 lineto
  0.273 -0.4189 lineto
  0.3536 -0.3536 lineto
  %middle circumferential
  0.3536 -0.3536 lineto
  0.413 -0.2818 lineto
  0.4594 -0.1974 lineto
  0.4894 -0.1024 lineto
  0.5 0.0 lineto
  %end circumferential
  0.5 0.0 lineto
  0.375 0.0 lineto
  %middle radial
  0.375 0.0 lineto
  0.25 0.0 lineto
  %end radial
  0.25 0.0 lineto
  0.2297 -0.0987 lineto
  0.1768 -0.1768 lineto
  %middle circumferential
  0.1768 -0.1768 lineto
  0.0922 -0.2324 lineto
  0.0 -0.25 lineto
  %end circumferential
  closepath
currentgray   0.45 setgray   fill   setgray

[2598 lines omitted]

% Vertical line from (0.0, -0.5) to (0.0, 0.0)
newpath
  0.0 -0.5 moveto
  0.0 0.0 lineto
0.05 setlinewidth   1 setlinecap  stroke
% Vertical line from (0.0, -0.5) to (0.0, -1.0)
newpath
  0.0 -0.5 moveto
  0.0 -1.0 lineto
0.05 setlinewidth   1 setlinecap  stroke
% Vertical line from (0.0, -2.0) to (0.0, -1.0)
newpath
  0.0 -2.0 moveto
  0.0 -1.0 lineto
0.05 setlinewidth   1 setlinecap  stroke
% Vertical line from (0.0, -2.0) to (0.0, -1.1579208923731617E77)
newpath
  0.0 -2.0 moveto
  0.0 -6.3553 lineto
  0.0 -6.378103166302659 lineto
  0.0 -6.378103166302659 lineto
  0.0 -6.378103166302659 lineto
0.05 setlinewidth 1 setlinecap stroke

[84 lines omitted]

% End of PostScript file for plot of function IDENTITY

Here is the program that generated the PostScript code for the graphs shown in figures 12.1 through 12.20. It contains a mixture of fairly general mechanisms and ad hoc kludges for plotting functions of a single complex argument while gracefully handling extremely large and small values, branch cuts, singularities, and periodic behavior. The aim was to provide a simple user interface that would not require the caller to provide special advice for each function to be plotted. The file for figure 12.1, for example, was generated by the call (picture ’identity), which resulted in the writing of a file named identity-plot.ps.

The program assumes that any periodic behavior will have a period that is a multiple of 2π; that branch cuts will fall along the real or imaginary axis; and that singularities or very large or small values will occur only at the origin, at ± 1 or ±i, or on the boundaries of the annuli (particularly those with radius π∕2 or π). The central function is parametric-path, which accepts four arguments: two real numbers that are the endpoints of an interval of real numbers, a function that maps this interval into a path in the complex plane, and the function to be plotted; the task of parametric-path is to generate PostScript code (a series of lineto operations) that will plot an approximation to the image of the parametric path as transformed by the function to be plotted. Each of the functions hline, vline, -hline, -vline, radial, and circumferential takes appropriate parameters and returns a function suitable for use as the third argument to parametric-path. There is some code that defends against errors (by using ignore-errors) and against certain peculiarities of IEEE floating-point arithmetic (the code that checks for not-a-number (NaN) results).

The program is offered here without further comment or apology.


(defparameter units-to-show 4.1)
(defparameter text-width-in-picas 28.0)
(defparameter device-pixels-per-inch 300)
(defparameter pixels-per-unit
  (* (/ (/ text-width-in-picas 6)
        (* units-to-show 2))
     device-pixels-per-inch))

(defparameter big (sqrt (sqrt most-positive-single-float)))
(defparameter tiny (sqrt (sqrt least-positive-single-float)))

(defparameter path-really-losing 1000.0)
(defparameter path-outer-limit (* units-to-show (sqrt 2) 1.1))
(defparameter path-minimal-delta (/ 10 pixels-per-unit))
(defparameter path-outer-delta (* path-outer-limit 0.3))
(defparameter path-relative-closeness 0.00001)
(defparameter back-off-delta 0.0005)

(defun comment-line (stream &rest stuff)
  (format stream "~%% ")
  (apply #’format stream stuff)
  (format t "~%% ")
  (apply #’format t stuff))

(defun parametric-path (from to paramfn plotfn)
  (assert (and (plusp from) (plusp to)))
  (flet ((domainval (x) (funcall paramfn x))
         (rangeval (x) (funcall plotfn (funcall paramfn x)))
         (losing (x) (or (null x)
                         (/= (realpart x) (realpart x))  ;NaN?
                         (/= (imagpart x) (imagpart x))  ;NaN?
                         (> (abs (realpart x)) path-really-losing)
                         (> (abs (imagpart x)) path-really-losing))))
    (when (> to 1000.0)
      (let ((f0 (rangeval from))
            (f1 (rangeval (+ from 1)))
            (f2 (rangeval (+ from (* 2 pi))))
            (f3 (rangeval (+ from 1 (* 2 pi))))
            (f4 (rangeval (+ from (* 4 pi)))))
        (flet ((close (x y)
                 (or (< (careful-abs (- x y)) path-minimal-delta)
                     (< (careful-abs (- x y))
                        (* (+ (careful-abs x) (careful-abs y))
                           path-relative-closeness)))))
          (when (and (close f0 f2)
                     (close f2 f4)
                     (close f1 f3)
                     (or (and (close f0 f1)
                              (close f2 f3))
                         (and (not (close f0 f1))
                              (not (close f2 f3)))))
            (format t "~&Periodicity detected.")
            (setq to (+ from (* (signum (- to from)) 2 pi)))))))
     (let ((fromrange (ignore-errors (rangeval from)))
          (torange (ignore-errors (rangeval to))))
      (if (losing fromrange)
          (if (losing torange)
              ’()
              (parametric-path (back-off from to) to paramfn plotfn))
          (if (losing torange)
              (parametric-path from (back-off to from) paramfn plotfn)
              (expand-path (refine-path (list from to) #’rangeval)
                           #’rangeval))))))

(defun back-off (point other)
  (if (or (> point 10.0) (< point 0.1))
      (let ((sp (sqrt point)))
        (if (or (> point sp other) (< point sp other))
            sp
            (* sp (sqrt other))))
      (+ point (* (signum (- other point)) back-off-delta))))

(defun careful-abs (z)
  (cond ((or (> (realpart z) big)
             (< (realpart z) (- big))
             (> (imagpart z) big)
             (< (imagpart z) (- big)))
         big)
        ((complexp z) (abs z))
        ((minusp z) (- z))
        (t z)))

(defparameter max-refinements 5000)

(defun refine-path (original-path rangevalfn)
  (flet ((rangeval (x) (funcall rangevalfn x)))
    (let ((path original-path))
      (do ((j 0 (+ j 1)))
          ((null (rest path)))
        (when (zerop (mod (+ j 1) max-refinements))
              (break "Runaway path"))
        (let* ((from (first path))
               (to (second path))
               (fromrange (rangeval from))
               (torange (rangeval to))
               (dist (careful-abs (- torange fromrange)))
               (mid (* (sqrt from) (sqrt to)))
               (midrange (rangeval mid)))
          (cond ((or (and (far-out fromrange) (far-out torange))
                     (and (< dist path-minimal-delta)
                          (< (abs (- midrange fromrange))
                             path-minimal-delta)
                          ;; Next test is intentionally asymmetric to
                          ;;  avoid problems with periodic functions.
                          (< (abs (- (rangeval (/ (+ to (* from 1.5))
                                                  2.5))
                                     fromrange))
                             path-minimal-delta)))
                 (pop path))
                ((= mid from) (pop path))
                ((= mid to) (pop path))
                (t (setf (rest path) (cons mid (rest path)))))))))
  original-path)

(defun expand-path (path rangevalfn)
  (flet ((rangeval (x) (funcall rangevalfn x)))
    (let ((final-path (list (rangeval (first path)))))
      (do ((p (rest path) (cdr p)))
          ((null p)
           (unless (rest final-path)
             (break "Singleton path"))
           (reverse final-path))
        (let ((v (rangeval (car p))))
          (cond ((and (rest final-path)
                      (not (far-out v))
                      (not (far-out (first final-path)))
                      (between v (first final-path)
                                 (second final-path)))
                 (setf (first final-path) v))
                ((null (rest p))   ;Mustn’t omit last point
                 (push v final-path))
                ((< (abs (- v (first final-path))) path-minimal-delta))
                ((far-out v)
                 (unless (and (far-out (first final-path))
                              (< (abs (- v (first final-path)))
                                 path-outer-delta))
                   (push (* 1.01 path-outer-limit (signum v))
                         final-path)))
                (t (push v final-path))))))))

(defun far-out (x)
  (> (careful-abs x) path-outer-limit))

(defparameter between-tolerance 0.000001)

(defun between (p q r)
  (let ((px (realpart p)) (py (imagpart p))
        (qx (realpart q)) (qy (imagpart q))
        (rx (realpart r)) (ry (imagpart r)))
    (and (or (<= px qx rx) (>= px qx rx))
         (or (<= py qy ry) (>= py qy ry))
         (< (abs (- (* (- qx px) (- ry qy))
                    (* (- rx qx) (- qy py))))
            between-tolerance))))

(defun circle (radius)
  #’(lambda (angle) (* radius (cis angle))))

(defun hline (imag)
  #’(lambda (real) (complex real imag)))

(defun vline (real)
  #’(lambda (imag) (complex real imag)))

(defun -hline (imag)
  #’(lambda (real) (complex (- real) imag)))

(defun -vline (real)
  #’(lambda (imag) (complex real (- imag))))

(defun radial (phi quadrant)
  #’(lambda (rho) (repair-quadrant (* rho (cis phi)) quadrant)))

(defun circumferential (rho quadrant)
  #’(lambda (phi) (repair-quadrant (* rho (cis phi)) quadrant)))

;;; Quadrant is 0, 1, 2, or 3, meaning I, II, III, or IV.

(defun repair-quadrant (z quadrant)
  (complex (* (+ (abs (realpart z)) tiny)
              (case quadrant (0 1.0) (1 -1.0) (2 -1.0) (3 1.0)))
           (* (+ (abs (imagpart z)) tiny)
              (case quadrant (0 1.0) (1 1.0) (2 -1.0) (3 -1.0)))))

(defun clamp-real (x)
  (if (far-out x)
      (* (signum x) path-outer-limit)
      (round-real x)))

(defun round-real (x)
  (/ (round (* x 10000.0)) 10000.0))

(defun round-point (z)
  (complex (round-real (realpart z)) (round-real (imagpart z))))

(defparameter hiringshade 0.97)
(defparameter loringshade 0.45)

(defparameter ticklength 0.12)
(defparameter smallticklength 0.09)

;;; This determines the pattern of lines and annuli to be drawn.
(defun moby-grid (&optional (fn ’sqrt) (stream t))
  (comment-line stream "Moby grid for function ~S" fn)
  (shaded-annulus 0.25 0.5 4 hiringshade loringshade fn stream)
  (shaded-annulus 0.75 1.0 8 hiringshade loringshade fn stream)
  (shaded-annulus (/ pi 2) 2.0 16 hiringshade loringshade fn stream)
  (shaded-annulus 3 pi 32 hiringshade loringshade fn stream)
  (moby-lines :horizontal 1.0 fn stream)
  (moby-lines :horizontal -1.0 fn stream)
  (moby-lines :vertical 1.0 fn stream)
  (moby-lines :vertical -1.0 fn stream)
  (let ((tickline 0.015)
        (axisline 0.008))
    (flet ((tick (n) (straight-line (complex n ticklength)
                                    (complex n (- ticklength))
                                    tickline
                                    stream))
           (smalltick (n) (straight-line (complex n smallticklength)
                                         (complex n (- smallticklength))
                                         tickline
                                         stream)))
      (comment-line stream "Real axis")
      (straight-line #c(-5 0) #c(5 0) axisline stream)
      (dotimes (j (floor units-to-show))
        (let ((q (+ j 1))) (tick q) (tick (- q))))
      (dotimes (j (floor units-to-show (/ pi 2)))
        (let ((q (* (/ pi 2) (+ j 1))))
          (smalltick q)
          (smalltick (- q)))))
    (flet ((tick (n) (straight-line (complex ticklength n)
                                    (complex (- ticklength) n)
                                    tickline
                                    stream))
           (smalltick (n) (straight-line (complex smallticklength n)
                                         (complex (- smallticklength) n)
                                         tickline
                                         stream)))
      (comment-line stream "Imaginary axis")
      (straight-line #c(0 -5) #c(0 5) axisline stream)
      (dotimes (j (floor units-to-show))
        (let ((q (+ j 1))) (tick q) (tick (- q))))
      (dotimes (j (floor units-to-show (/ pi 2)))
        (let ((q (* (/ pi 2) (+ j 1))))
          (smalltick q)
          (smalltick (- q)))))))

(defun straight-line (from to wid stream)
  (format stream
          "~%newpath  ~S ~S moveto  ~S ~S lineto  ~S ~
           setlinewidth  1  setlinecap  stroke"
          (realpart from)
          (imagpart from)
          (realpart to)
          (imagpart to)
          wid))

;;; This function draws the lines for the pattern.
(defun moby-lines (orientation signum plotfn stream)
  (let ((paramfn (ecase orientation
                   (:horizontal (if (< signum 0) #’-hline #’hline))
                   (:vertical (if (< signum 0) #’-vline #’vline)))))
    (flet ((foo (from to other wid)
             (ecase orientation
               (:horizontal
                (comment-line stream
                              "Horizontal line from (~S, ~S) to (~S, ~S)"
                              (round-real (* signum from))
                              (round-real other)
                              (round-real (* signum to))
                              (round-real other)))
               (:vertical
                (comment-line stream
                              "Vertical line from (~S, ~S) to (~S, ~S)"
                              (round-real other)
                              (round-real (* signum from))
                              (round-real other)
                              (round-real (* signum to)))))
             (postscript-path
               stream
               (parametric-path from
                                to
                                (funcall paramfn other)
                                plotfn))
             (postscript-penstroke stream wid)))
      (let* ((thick 0.05)
             (thin 0.02))
        ;; Main axis
        (foo 0.5 tiny 0.0 thick)
        (foo 0.5 1.0 0.0 thick)
        (foo 2.0 1.0 0.0 thick)
        (foo 2.0 big 0.0 thick)
        ;; Parallels at 1 and -1
        (foo 2.0 tiny 1.0 thin)
        (foo 2.0 big 1.0 thin)
        (foo 2.0 tiny -1.0 thin)
        (foo 2.0 big -1.0 thin)
        ;; Parallels at 2, 3, -2, -3
        (foo tiny big 2.0 thin)
        (foo tiny big -2.0 thin)
        (foo tiny big 3.0 thin)
        (foo tiny big -3.0 thin)))))

(defun splice (p q)
  (let ((v (car (last p)))
        (w (first q)))
    (and (far-out v)
         (far-out w)
         (>= (abs (- v w)) path-outer-delta)
         ;; Two far-apart far-out points.  Try to walk around
         ;;  outside the perimeter, in the shorter direction.
         (let* ((pdiff (phase (/ v w)))
                (npoints (floor (abs pdiff) (asin .2)))
                (delta (/ pdiff (+ npoints 1)))
                (incr (cis delta)))
           (do ((j 0 (+ j 1))
                (p (list w "end splice") (cons (* (car p) incr) p)))
               ((= j npoints) (cons "start splice" p)))))))

;;; This function draws the annuli for the pattern.
(defun shaded-annulus (inner outer sectors firstshade lastshade fn stream)
  (assert (zerop (mod sectors 4)))
  (comment-line stream "Annulus ~S ~S ~S ~S ~S"
                (round-real inner) (round-real outer)
                sectors firstshade lastshade)
  (dotimes (jj sectors)
    (let ((j (- sectors jj 1)))
      (let* ((lophase (+ tiny (* 2 pi (/ j sectors))))
             (hiphase (* 2 pi (/ (+ j 1) sectors)))
             (midphase (/ (+ lophase hiphase) 2.0))
             (midradius (/ (+ inner outer) 2.0))
             (quadrant (floor (* j 4) sectors)))
        (comment-line stream "Sector from ~S to ~S (quadrant ~S)"
                      (round-real lophase)
                      (round-real hiphase)
                      quadrant)
        (let ((p0 (reverse (parametric-path midradius
                                            inner
                                            (radial lophase quadrant)
                                            fn)))
              (p1 (parametric-path midradius
                                   outer
                                   (radial lophase quadrant)
                                   fn))
              (p2 (reverse (parametric-path midphase
                                            lophase
                                            (circumferential outer
                                                             quadrant)
                                            fn)))
              (p3 (parametric-path midphase
                                   hiphase
                                   (circumferential outer quadrant)
                                   fn))
              (p4 (reverse (parametric-path midradius
                                            outer
                                            (radial hiphase quadrant)
                                            fn)))
              (p5 (parametric-path midradius
                                   inner
                                   (radial hiphase quadrant)
                                   fn))
              (p6 (reverse (parametric-path midphase
                                            hiphase
                                            (circumferential inner
                                                             quadrant)
                                            fn)))
              (p7 (parametric-path midphase
                                   lophase
                                   (circumferential inner quadrant)
                                   fn)))
          (postscript-closed-path stream
            (append
              p0 (splice p0 p1) ’("middle radial")
              p1 (splice p1 p2) ’("end radial")
              p2 (splice p2 p3) ’("middle circumferential")
              p3 (splice p3 p4) ’("end circumferential")
              p4 (splice p4 p5) ’("middle radial")
              p5 (splice p5 p6) ’("end radial")
              p6 (splice p6 p7) ’("middle circumferential")
              p7 (splice p7 p0) ’("end circumferential")
              )))

        (postscript-shade stream
                          (/ (+ (* firstshade (- (- sectors 1) j))
                                (* lastshade j))
                             (- sectors 1)))))))

(defun postscript-penstroke (stream wid)
  (format stream "~%~S setlinewidth   1 setlinecap  stroke"
          wid))

(defun postscript-shade (stream shade)
  (format stream "~%currentgray   ~S setgray   fill   setgray"
          shade))

(defun postscript-closed-path (stream path)
  (unless (every #’far-out (remove-if-not #’numberp path))
    (postscript-raw-path stream path)
    (format stream "~%  closepath")))

(defun postscript-path (stream path)
  (unless (every #’far-out (remove-if-not #’numberp path))
    (postscript-raw-path stream path)))

;;; Print a path as a series of PostScript "lineto" commands.
(defun postscript-raw-path (stream path)
  (format stream "~%newpath")
  (let ((fmt "~%  ~S ~S moveto"))
    (dolist (pt path)
      (cond ((stringp pt)
             (format stream "~%  %~A" pt))
            (t (format stream
                       fmt
                       (clamp-real (realpart pt))
                       (clamp-real (imagpart pt)))
               (setq fmt "~%  ~S ~S lineto"))))))

;;; Definitions of functions to be plotted that are not
;;; standard Common Lisp functions.

(defun one-plus-over-one-minus (x) (/ (+ 1 x) (- 1 x)))

(defun one-minus-over-one-plus (x) (/ (- 1 x) (+ 1 x)))

(defun sqrt-square-minus-one (x) (sqrt (- 1 (* x x))))

(defun sqrt-one-plus-square (x) (sqrt (+ 1 (* x x))))

;;; Because X3J13 voted for a new definition of the atan function,
;;; the following definition was used in place of the atan function
;;; provided by the Common Lisp implementation I was using.

(defun good-atan (x)
  (/ (- (log (+ 1 (* x #c(0 1))))
        (log (- 1 (* x #c(0 1)))))
     #c(0 2)))

;;; Because the first edition had an erroneous definition of atanh,
;;; the following definition was used in place of the atanh function
;;; provided by the Common Lisp implementation I was using.

(defun really-good-atanh (x)
  (/ (- (log (+ 1 x))
        (log (- 1 x)))
     2))

;;; This is the main procedure that is intended to be called by a user.
(defun picture (&optional (fn #’sqrt))
  (with-open-file (stream (concatenate ’string
                                       (string-downcase (string fn))
                                       "-plot.ps")
                          :direction :output)
    (format stream "% PostScript file for plot of function ~S~%" fn)
    (format stream "% Plot is to fit in a region ~S inches square~%"
            (/ text-width-in-picas 6.0))
    (format stream
            "%  showing axes extending ~S units from the origin.~%"
            units-to-show)
    (let ((scaling (/ (* text-width-in-picas 12) (* units-to-show 2))))
      (format stream "~%~S ~:*~S scale" scaling))
    (format stream "~%~S ~:*~S translate" units-to-show)
    (format stream "~%newpath")
    (format stream "~%  ~S ~S moveto" (- units-to-show) (- units-to-show))
    (format stream "~%  ~S ~S lineto" units-to-show (- units-to-show))
    (format stream "~%  ~S ~S lineto" units-to-show units-to-show)
    (format stream "~%  ~S ~S lineto" (- units-to-show) units-to-show)
    (format stream "~%  closepath")
    (format stream "~%clip")
    (moby-grid fn stream)
    (format stream
            "~%% End of PostScript file for plot of function ~S"
            fn)
    (terpri stream)))