Dogs and fat tails

Posted on by John

I was reading a blog post on boat names because it was on Hacker News this morning. The post contained a link to a data set on dog names in NYC and I poked around the data a little. The top names were not at all what I expected, but then again this is limited to NYC; it’s not a sample across the US. These were the top 10 names:

  1. Bella
  2. Luna
  3. Max
  4. Charlie
  5. Coco
  6. Lola
  7. Rocky
  8. Milo
  9. Teddy
  10. Lucy

I wondered if the name frequencies might fit a power-law distribution. They do not, but they follow a log-normal distribution remarkably well.

Related posts

Manually unbreakable cryptography

Posted on by John

Suppose you were able to go back in time, to an era before computers, and give someone contemporary cryptography. Encryption methods that are essentially unbreakable now would certainly be unbreakable then. But there’s a catch: not only do attackers not have computers, neither do users.

Manual cryptography

If you told someone about RSA encryption, for example, you’d lose them right after you said “First you find a couple 1000-digit primes.” But there’s no need for using 1000-digit primes if no attacker has a computer. You could use 100 digit primes. Could you use 10 digit primes? If you chose primes just big enough to make the method unbreakable by hand, could someone implement it by hand?

Kirchoff’s principle says the strength of an encryption method should depend only on keeping the key secret, not the method. If you could keep the method secret, RSA would be unbreakable because nobody thought of anything like it before computers. But to make our thought experiment more interesting, let’s suppose that an enemy has also traveled back in time. If you tell your side about RSA, he can tell his side about it as well. So we’re back to Kirchoff’s principle.

An encryption method combinining substitution and permutation would have been practical to carry our manually. The ADFGVX cipher from 1918 was a start in this direction. That idea could been extended further, with a larger substitution set and longer permutations, and with more than one round of substitution and permutation, approaching what would be come the approach used in modern symmetric encryption. Such a method might have been manually implementable without being manually breakable.

Mechanized cryptography

World War II was a time of transition from manual cryptography to computerized cryptography. Encryption machines were attacked by cryptanalysis machines, though these machines were general-purpose computers. If you could implement a symmetric encryption method like AES in a mechanical device, no mechanical device could break it.

You could use something like DES, simpler than AES but still unbreakable at the time. DES is considered broken because now you could throw enough compute power at it to break it by brute force, but that would not be possible with only mechanical devices.

My hunch is that the best approach would be stream ciphers. Maybe it would be practical to implement one of these by hand or with the aid of simple calculating machines. Something like PCG, which is not cryptographically secure today [1], would have been then, though I don’t know how practical it would have been to carry out PCG, say, in the 1940s.

More pre-computer cryptography

[1] In 2020, Charles Bouillaguet, Florette Martinez, and Julia Sauvage were able to break PCG using 20,000 CPU-hours. See their paper Practical seed-recovery for the PCG Pseudo-Random Number Generator. IACR Transactions on Symmetric Cryptology. ISSN 2519-173X, Vol. 2020, No. 3, pp. 175–196.

Learning from historical mistakes

Posted on by John

The following extraordinary paragraph comes from Knuth’s TAOCP Volume 4A, right before the last set of exercises.

Many of the exercises below ask a modern reader to find and/or to correct errors in the literature of bygone days. The point is not to gloat over how smart we are in the 21st century; the point is rather to understand that even the pioneers of a subject can stumble. One good way to learn that a set of ideas is not really as simple as it might seem to today’s computer scientists and mathematicians is to observe that some of the world’s leading thinkers had to struggle with the concepts when they were new.

Inverse differential equations

Posted on by John

In science and engineering classes, you might describe a system using Newton’s laws and end up with a differential equation. You then solve the differential equation, analytically or numerically, to see how the solutions behave.

You might also do the opposite, especially in a mathematics class: look at what differential equation a set of functions satisfy in order to understand those functions.

Bessel functions came out of solving differential equations from astronomy. But then they turned out to be useful, not just in other areas of science, but in pure mathematics as well. So there are Bessel function users who came to the functions first and haven’t seen the differential equation they came from.

You can learn a lot about Bessel functions, and other special functions, by looking at their defining differential equation even if you’re not directly interested in the differential equation or the physical problem that motivated it.

Bessel functions satisfy

x² y″ + x y′ + (x² − n²) y = 0.

You can tell a lot about Bessel functions just by inspecting this equation without solving it. If we divide by x² and write the equation in the form

y″ + (p(x)/x) y′ + (q(x)/x²) y = 0

then p(x) = 1 and q(x) = x² − n². The indicial equation

r(r − 1) + p(0) r + q(0) = 0

reduces to

r² = n²

and so r = ±n. That alone tells us there are two solutions, one analytic at zero and one singular at zero. These are Jn and Yn respectively. It also tells us the behavior of these functions as x goes to zero and as x goes to infinity. To find out more, look up “method of Frobinius.”

DNA and Bessel functions

Posted on by John

I was reading a book on the history of the discovery of the structure of DNA [1] and was surprised by a few passing references to Bessel functions.

According to Claude,

When X-rays are diffracted by a helical structure, the resulting diffraction pattern breaks into a series of horizontal “layer lines.” Cochran, Crick, and Vand showed mathematically that the diffracted amplitude on the n-th layer line is proportional to a Bessel function of the first kind, order n:

Jn(2π r R)

where r is the radius of the helix and R is the distance out from the center (the meridian) in the diffraction pattern.

The citation for this paragraph is a paper from 1952 [2] that amazingly is behind a paywall.

Update: Evan pointed out something fishy in the expression above: “Isn’t there a dimensional problem with J(2π r R)? Maybe 2πr/R instead?”

This was Claude’s response:

  • r is the radius of the helix, a real-space length (say, in Å).
  • R is the radial coordinate in reciprocal space — it has units of inverse length (Å⁻¹), because that’s how X-ray diffraction/Fourier space works (spatial frequency, not distance).

So Evan and Claude were both right. There is an implicit division happening, but crystallographers do write r R.

Related posts

[1] Watson and Crick didn’t “discover DNA” as is commonly said. DNA was discovered in 1878. Watson and Crick discovered the structure of DNA in 1953.

[2] Cochran, W., Crick, F. H. C., & Vand, V. (1952). “The Structure of Synthetic Polypeptides. I. The Transform of Atoms on a Helix.” Acta Crystallographica, 5(5), 581–586.

A simple range reduction method

Posted on by John

At the end of my post on how not to calculate cosine I said that the first step in calculating cosine, particularly cosine of a large number, would be to do range reduction. This post will present a simple range reduction method by Cody and Waite that is adequate for moderately large arguments.

If you want to compute the sine or cosine of an angle x you could start by reducing x mod 2π since that would not change the result. However, accurately reducing a number mod 2π is not trivial; that’s why range reduction is an area of algorithm development.

Range reduction mod π/2

Even better would be to reduce x mod π/2. Reducing to a smaller range means that power series method, and other methods such as rational approximation, will be more efficient.

So suppose you can find an integer k such that

x − k π/2 = y

where 0 ≤ y ≤ π/2. Then sin(x) is ±sin(y) or ±cos(y), depending on k mod 4 equals 0, 1, 2, or 3.

from math import *

def reduced_sin(x, k):
     match k % 4:
        case 0: return sin(x)
        case 1: return cos(x)
        case 2: return -sin(x)
        case 3: return -cos(x)

Naive range reduction

Now let’s set x = 500. Then k = 318 because that’s the multiple of π/2 we need to subtract to bring x into range, and the sine of x should be the negative of the sine of the reduced value y because 318 = 2 mod 4.

The following code computes sin(x) with naive range reduction

def naive_sin(x):
    k = floor(x / (pi/2))
    y = x % (pi/2)
    return reduced_sin(y, k)

and when x = 500 the error is on the order of 1.7 × 10−14.

Better range reduction

The value of k above is fine, but we’d like to calculate y more accurately. The following code is much better.

def Cody_Waite_sin(x):
    C1 = 1686629713 / 2**30
    C2 = 4701928774853425 / 2**86

    k = floor(x / (pi/2))
    y = (x - k*C1) - k*C2
    return reduced_sin(y, k)

This will compute sin(500) to full machine precision. What kind of magic is this?

The trick is that the exact value of C1 + C2 equals π/2 to more precision than is possible in a single float [1]. You can confirm, with bc or some other extended precision software, that the difference between C1 + C2 and π/2 is roughly 2−88, while the limit of float precision is 2−52.

If we compute

y = x - k*(C1 + C2)

then we’re doing the same calculation as naive_sin and will get the same error. But if we compute

y = (x - k*C1) - k*C2

we will get a more accurate result, provided x isn’t too large.

You can use the following code to play around and see how large x can be before errors start to creep in. For small enough x, like 500, the Cody and Waite sine returns full precision. For larger x it’s better than naive sine but does not return full precision. And for large enough x it completely breaks down.

def compare(x):
    y0 = naive_sin(x) 
    y1 = Cody_Waite_sin(x)
    y2 = sin(x)
    print("Naive error:     ", y2 - y0)
    print("Cody Waite error:", y2 - y1)

Now this may seem circular since we’re using math.sin as our gold standard. However, this function is calling the sine function on your CPU, which is using sophisticated range reduction to compute its result accurately down to the last bit, assuming you run the code on a computer that’s less than 40 years old.

The Cody and Waite algorithm is inadequate for large x, but it’s a good place to begin studying range reduction. It shows there are clever ways of squeezing out more precision than seems possible.

 

[1] The numerator n1 of C1 is ⌊230 π/2⌋. The numerator n2 of C2 is the solution to

286−30 n1 + n2 = ⌊286 π/2⌋.

Corrupted apostrophes

Posted on by John

I have a program that shares files between my laptop and my phone. It works well, except for apostrophes.

When I type an apostrophe ' on my laptop, it becomes ’ on my phone. And when I type 's on my phone, it becomes 痴 on my laptop.

Apparently the phone turns the apostrophe (U+0027) into a right single quote (U+2019), then bungles bytes in the UTF-8 encoding of U+2019 as three Windows-1252 characters. The bytes E28099hex are interpreted as â (E2hex), € (80hex), and ™ (99hex).

When I type 's on my phone, it is encoded as two Windows-1252 characters 92hex and 73hex. Then by the time the text appears on my laptop, the bytes 9273hex are interpreted as a Shift-JIS encoding of the CJK character 痴 (U+75F4).

Here’s Python code to reproduce the problem.

def mojibake(s: str, n: int, bad_encoding: str = 'cp1252') -> str:
    for _ in range(n):
        s = s.encode('utf-8').decode(bad_encoding, errors='replace')
    return s

print(mojibake("’s", 1))
print(mojibake("’s", 1, 'shift_jis'))

The code won’t corrupt ASCII text. The problem started with an ASCII character being replaced by a similar non-ASCII character. Text containing a non-ASCII character gets more corrupted with each round.

How not to calculate cosine

Posted on by John

Calculus professors with no experience in numerical computing will tell students that computers calculate trig functions with power series. They don’t. I worked on the implementation of trig functions in hardware, and I can assure you we didn’t just use power series.

Power series are an excellent way to calculate functions near the center of the series, such as computing sine for small angles. But the further you get from the center, the less useful power series are.

Let’s suppose you want to calculate cos(200) using the power series for cosine. The nth term of that series is

(−1)n x2n / (2n)!

This is an alternating series, and so the error in truncating the series after n terms is bounded by the size of the n+1 term, if you’ve gone far enough out in the series that the terms are monotonically decreasing in absolute value.

To calculate cos(200) to machine precision, i.e. with an error of less than 2−52, we’d need to sum the series up to n where

| 2002n+2 / (2n + 2)! | < 2−52

Actually, that will ensure that the absolute error is small enough, but not that the relative error is small enough; if the value of cos(200) is small, we’d need more terms. Let’s ignore that and assume we’re only concerned with absolute error.

Turns out we’d need 287 terms. That’s a lot of terms. But you might say “That’s fine. I’m not in a hurry, and it’s just more work for the computer, not for me.” OK, so let’s try.

from math import *

s = 0
for n in range(288):
    s += (-1)**n * 200**(2*n) / factorial(2*n)
print(s)

This prints -3.6840358571084123e+67. You may suspect the answer is incorrect since values of cosine are on the order of 1, not on the order of 1067. Something went spectacularly bad. On closer inspection, it’s remarkable the code didn’t crash.

If you changed 200 to 200.0 above, the code would crash. Calculating 200.0**(2*n) overflows when n = 67. But when we calculate 200**(2*n), the result is an integer. And we’re dividing by factorial(2*n), which is also an integer. Both of these integers become too large to fit in a float, but their ratio has a maximum value of around 1080, smaller than the maximum float, which is on the order of 10308.

When we don’t overflow, we have a different problem: catastrophic cancellation. You can’t calculate a number between −1 and 1 as an alternating sum of numbers as large as 1080. You’d need more than 80 + 16 = 96 decimal places of precision to compute the sum accurately, and floating point only gives you between 15 and 16 decimal places of precision.

So how would you calculate cos(200)? The first step would be to use some sort of range reduction on 200. You could reduce 200 mod 2π to get a smaller number to work with.

>>> from math import cos, pi
>>> x = 200 % (2*pi)
>>> x
5.221255477432827

Using a power series to compute the cosine of 5.221255477432827 is feasible, but not optimal. There’s also another problem: the naive range reduction above loses some precision.

>>> cos(x)
0.48718767500701254
>>> cos(x) - cos(200)
6.661338147750939e-15

The error is small, but it’s still an order of magnitude larger than machine precision. You can’t simply reduce n mod 2π with ordinary float division because the integer part of n / 2π pushes some digits of precision off the right end. I intend to write about how range reduction works in future posts.

Update: See this post for a simple range reduction algorithm that is fine for values of x such as 200, but not adequate for much larger values.

cos(200!)

Posted on by John

In a footnote to the previous post, I said that Python’s math library can calculate the logarithm of extremely large numbers but not the cosine. This post will expand on that comment.

In this post I’ll use n = 200! as my example rather than 1000! because this value of N is larger than the largest representable floating point number but small enough to be more convenient to work with.

Suppose someone calculates 200! for you:

78865786736479050355236321393218506229513597768717326329474253324435\
94499634033429203042840119846239041772121389196388302576427902426371\
05061926624952829931113462857270763317237396988943922445621451664240\
25403329186413122742829485327752424240757390324032125740557956866022\
60319041703240623517008587961789222227896237038973747200000000000000\
00000000000000000000000000000000000

You could now calculate log(n) using

n = 7.886578673647905 × 10374

and so

log(n) = log(7.886578673647905 × 10374)
= log(7.886578673647905) + 374 log(10) = 863.2319871924055.

The key thing that makes this possible is that the least significant digits of n only affect the least significant digits of log(n). In the calculation above I kept the first 16 digits of n. Python couldn’t make use of any more digits, and had no need of any more digits, in order to produce the logarithm to machine precision.

Cosine doesn’t work that way. The cosine of n depends on the remainder when n is divided by 2π, and that remainder depends on every single digit of n. I’ll illustrate that below.

Using bc -l and setting the scale to 400, I can calculated n then calculate

cos(n + 10i)

for i running from 0 to 374, tweaking each digit one at a time. (Except when a digit is a 9 and the addition results in a carry.)

    n = 1
    for (i = 1; i <= 200; i++) n *= i
    scale = 400
    for (i = 1; i <= 374; i++) {
        x = c(n+10^i)
        scale = 16
        print x/1, "\n"
        scale = 400
    }

Here’s what a plot of the results look like.

The value of cos(n) is about −0.985, but the values above are all over the map. We can look at the range by projecting all the points over to the left edge then rotating a quarter turn:

The remarkable thing about this image is that there are a few gaps, i.e. a few values the cosine does not take on.

Here’s a more sophisticated way to look at it. The sequence 10i mod 2π is dense in [0, 2π], and so by going far enough out in the sequence, we can find a value that shifts the phase of n by any desired amount within any given tolerance.

Every digit in n matters, and changing any digit can change the value of cosine to be essentially any value. You cannot calculate the cosine of an enormous number without using some kind of extended precision arithmetic. There are clever range reduction algorithms that minimize the amount of extended arithmetic necessary, but extended arithmetic cannot be completely eliminated.

Calculating log(1000!)

Posted on by John

The previous post pointed out that the following code such as the following unexpectedly works.

>>> from math import log, factorial
>>> log(factorial(1000))
5912.128178488163

If you don’t find this unexpected, note that if you replace math.log with numpy.log the code will fail [1]. Functions like natural logarithm operate on real numbers. Real numbers are represented as floating point numbers in programming languages, and 1000! factorial is too large to represent as a standard floating point number. (More on that here.)

In this post I’d like to look at how you might calculate log(1000!) with less capable software, and even without software.

One approach would be to sum the logarithms of the numbers 1 through 1000. This will give essentially the same result as above, with a little difference in the last couple decimal places due to rounding error.

If you have a way to calculate 1000! but not a way to cast it to a floating point number, you could do this manually.

>>> s = str(factorial(1000))
>>> s[:16]
'4023872600770937'
>>> len(s)
2568

This tells us 1000! = 4.023872600770937 × 102567. Therefore

log(1000!) = log(4.023872600770937) + 2567 log(10)

which only requires working with numbers of modest size.

Calculating by hand

Now suppose it’s 1964. You don’t have a computer, or even a calculator, but you do have a copy of the recently published Handbook of Mathematical Functions by Abramowitz and Stegun (A&S). You turn to Table 6.6 “Factorials for large arguments.” This has values of factorial for 100, 200, 300, …, 1000, so you can simply look up your answer to 20 decimal places.

That was too easy; I didn’t expect that to be there when I started writing this post. If you wanted to compute log(950!), for example, you’d have to work harder. You could find A&S equation 6.1.41 (Stirling’s series) which says

\begin{align*} \ln \Gamma(z) &\sim (z - \tfrac{1}{2}) \ln z - z + \tfrac{1}{2} \ln 2\pi + \frac{1}{12z} - \frac{1}{360z^3} \\ &+ \frac{1}{1260z^5} - \frac{1}{680z^7} + \cdots \end{align*}

So how would you use this formula to calculate log(1000!)? Since n! = Γ(n + 1), you set z = 1001.

You’d need to decide how many terms you need to use. Assuming the error is on the order of the first term you leave out, you’d reason that you could probably stop with the 1/12z term because the next term is between 10−11 and 10−12.

You find Table 4.2 has natural logarithms, but not for 1001. You can look up log(1.001), however, and at the bottom of the same page is log(10) to 16 decimal places, and you can find log(10) to 24 decimal places in Table 1.1. So you calculate

log(1001) = log(1.001 × 10³) = log(1.001) + 3 log(10).

You can find log(2) and log(π) in Table 1.1, and average them to find ½ log(2π).

Here’s Python code to simulate the hand calculations.

log2     = 0.6931_47180_55994_53094_172321 # Table 1.1
log10    = 2.3025_85092_99404_56840_179915 # Table 1.1
logpi    = 1.1447_29885_84940_01741_43427  # Table 1.1
log1_001 = 0.00099_95003_330835            # Table 4.2

z = 1001
logz = log1_001 + 3*log10
s = (z - 0.5)*logz - z + (log2 + logpi)/2 + 1/(12*z)

print(s)

This result differs from the one at the top of the post only in the last decimal place.

Related posts

Doing calculations with tables is not as simple as “just look it up.” It takes a bit of skill.

[1] The code will also fail if you replace math.log with math.cos. Both logarithm and cosine return moderate sized real numbers when given enormous inputs like 1000!, so representing the output as a float is not the problem. But logarithms of huge numbers can be computed with ordinary precision functions, as above. But computing the cosine of a huge number requires extended precision.

Update: The next post expands on why computing the cosine of a large number is more difficult than computing the log.