John D. CookModify

2026-08-08 09:20:25

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 E28099_hex are interpreted as â (E2_hex), € (80_hex), and ™ (99_hex).

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

2026-08-07 23:11:10

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)ⁿx²ⁿ / (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⁻⁵², we’d need to sum the series up to n where

| 200²ⁿ⁺² / (2n + 2)! | < 2⁻⁵²

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 10⁶⁷. 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 10⁸⁰, smaller than the maximum float, which is on the order of 10³⁰⁸.

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 10⁸⁰. 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.

2026-08-07 21:19:07

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 × 10³⁷⁴

and so

log(n) = log(7.886578673647905 × 10³⁷⁴)
= 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 + 10ⁱ)

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 10ⁱ 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.

2026-08-06 21:23:43

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 × 10²⁵⁶⁷. 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

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⁻¹¹ and 10⁻¹².

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.

• Using a table of logarithms
• Using a trig table
• How much precision can you squeeze out of a table?

[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.

2026-08-06 02:26:25

Last week I wrote a post on hiding cryptographic keys in decks of cards. I wrote some code for that post that shouldn’t work, but before fixing I noticed that it in fact did work.

The code computes logarithms for integers larger than the largest representable float. For example, the largest float is on the order of 10³⁰⁸, and yet the following code works.

>>> import math
>>> math.log10(10**400)
400.0

The log, log2, and log10 functions have some code inside that handles large integers specially. It doesn’t simply convert the integers to floats before taking the logarithm. If it did, it would overflow. If you replace math with numpy above, the code will fail. NumPy’s implementation of logarithms is more what I would expect.

While playing around with this I also noticed that you can define floats larger than the largest float without warnings.

>>> math.log(1e308)
709.1962086421661
>>> math.log(1e309)
inf

This isn’t a feature of math.log but of how Python handles scientific notation. The expression 1e308 is the floating point representation of 10³⁰⁸. It is a float, not an int.

>>> type(1e308)
<class 'float'>

The expression 1e309 is also a float. But since it’s larger than is possible for a float, Python interprets it as inf. The code

math.log(1e309)

returns inf based on the reasoning that log(∞) = ∞.

That explains the following behavior:

>>> 1e309 == 1e310
True

The expressions 1e309 and 1e310 are equal because both are alternate ways of writing inf.

2026-08-05 22:49:35

A few days ago I wrote a post on counting rooted trees. That post looked at the sequence c(n) which counts the number of rooted trees with n nodes. Here one node is distinguished as the root, but the nodes below the root are not distinguished from each other; all that matters is how the nodes are connected.

The number of rooted trees with n nodes is the same as the number of ways to configure n − 1 non-overlapping circles. Not only are the counts the same, there is a natural correspondence between the trees and the circles. It’s not obvious that there should be such a correspondence, with the right notation the correspondence is sort of a pun.

The standard way to represent unlabeled trees is as a multiset of their children. We use a multiset, not a set, because some elements will be repeated. We represent a leaf as a pair of parentheses: ().

There is only one rooted tree with one node: ().

There is only one rooted tree with one two nodes: (()). Here the outer parentheses represent the root node and the inner parentheses represent its child.

There are two rooted trees with three nodes, and we can represent them as ((())) and ((),()). The first is the straight line tree: a node that has a single child node that has a single child node. The second is a node that branches to two nodes. (Here’s where we need multisets.)

The four rooted trees with four nodes can be represented as (((()))), ((((),())), ((),(())), and ((),(),(),()).

Here are the nine rooted trees with five nodes:

((((()))))
((((),())))
(((),(())))
(((),(),()))
((()),(()))
((),((())))
((),((),()))
((),(),(()))
((),(),(),())

The correspondence with non-overlapping circles removes the outer parentheses then joins the rest to form circles, with nested parentheses corresponding to concentric circles. A more geometric way to see the correspondence is to start at the bottom of the tree, replace leaves with circles, then work your way up circling connected components.