Eli Bendersky's website

Monte-Carlo simulations

2026-10-07T19:34:00-07:00

Monte Carlo simulations (or methods) is the technique of applying randomness and the Law of large numbers to the solution of various scientific and engineering problems. One of its first documented uses was by Stanislaw Ulam and John von Neumann for nuclear weapon simulations after WWII [1].

In this post I want to provide examples of some simple uses of Monte Carlo simulations. We'll start with the classical example of calculating the value of \pi by throwing darts.

Estimating pi

Suppose we take a square board and inscribe a quarter of a circle into it. We then proceed to throw darts at the board and record whether each dart hits inside or outside the quarter circle. Having thrown many such darts, we calculate the ratio of the darts inside the circle to the total number thrown.

Assuming our darts are distributed uniformly over the square, by the Law of large numbers this ratio should approach the ratio of areas of the quarter circle A_{circ} to the full square A_{square}.

With a square side length of 1, we have:

\[\frac{A_{circ}}{A_{square}}=\frac{\pi/4}{1}=\frac{\pi}{4}\]

Therefore:

\[\pi\approx 4\frac{\text{hits inside}}{\text{total throws}}\]

Here's a visualization:

Samples 200

101000

π ≈ —

Change the number of samples (dart throws) and click "Run" to regenerate. The code is very simple - here's a slightly sanitized version:

const total = Number(samples.value);
let inside = 0;
for (let i = 0; i < total; i++) {
  const x = Math.random();
  const y = Math.random();
  const inCircle = x * x + y * y <= 1;
  if (inCircle) inside++;
}
estimateValue = 4 * inside / total;

You'll notice that the estimate is relatively poor - even with 1000 samples - if you click "Run" several times, some numbers will be way off mark. While this method does estimate \pi, it's not a particularly good estimate. I find that running ~10 billion samples is necessary to estimate it to 4 digits after the decimal with reasonable reliability.

In general, for independent trials like these, the typical sampling error decreases in proportion to 1/\sqrt{N}, where N is the number of trials. This means that halving the error requires four times as many samples.

While the \pi estimation may seem whimsical, it's an example of an important class of problems to which Monte Carlo simulation is applied: numerical integration. Our simulation estimates the area under the quarter-circle curve, which is a definite integral.

Many integrals are very difficult to solve analytically, and much research has been done in the area of numerical analysis to develop methods to calculate integrals. Monte Carlo methods are particularly useful for high-dimensional integrals, where other numerical methods can become prohibitively expensive.

Combinatorial simulation - the game of SET

A common use of Monte Carlo methods is estimating complex combinatorial calculations. These often don't have analytical solutions, and enumerating all options is intractable due to the scale of the numbers involved. As an example, let's consider the game of SET. Each SET card has four attributes:

Number of shapes (1, 2 or 3)
Color (Red, Green or Purple)
Shape type (Oval, Diamond or Squiggle)
Shading (Empty, Striped or Solid)

And the goal is to find a "set" - three cards that are either all different or all the same for each attribute separately. As an example, here's a hand with a single set; see if you can find it [2]:

SET hand with a single set

And the next hand doesn't have any sets:

SET hand with no hands

Here's a question: given a freshly shuffled SET deck, what are the odds that the first 12 cards drawn will have no sets among them? This question is difficult to answer without using a computer.

It's easy to calculate the number of ways to deal a 12-card hand from a deck of 81:

\[\binom{81}{12}=70,724,320,184,700\]

But how many of these hands have no sets? Enumerating 70 trillion SET hands and checking each one can take quite a while, and there is no straightforward counting formula to answer this question. Some clever methods can be employed to leverage symmetries and other mathematical properties of SET to cut down this search space considerably. Donald Knuth himself worked on this problem and came up with a neat program (setset-all on his programs page) that found 2,284,535,476,080 such hands. Therefore, the answer to our question is:

\[\frac{2,284,535,476,080}{70,724,320,184,700}\approx 0.0323\]

There's a 3.23% chance that a randomly drawn hand of 12 cards from a full deck of SET will have no set in it.

Let's see how we can use a Monte Carlo simulation to answer this question with relatively small effort, without deep knowledge of the mathematical properties of SET that enable cutting down the search space Knuth-style. We can use the following pseudo-code:

C = 0
run N times:
  draw a random 12-card hand from a fresh deck
  count sets in the hand
  if no sets:
    C += 1

Estimated probability = C / N

After running 10 million simulated draws, I got an answer of 0.0323, which matches the real answer very closely.

The Monte Carlo approach lets us solve rather complicated problems in a very simple way. Suppose we want to answer the same question for a hand of 15 cards; this would blow up the search space considerably - there are about 100x more ways to select 15-card hands than there are to select 12-card hands. But for a Monte Carlo simulation, we adjust one small parameter and get a very reliable [3] answer (about 0.00037, in case you were wondering).

Retirement projection

One domain where Monte Carlo simulations are ubiquitous is projections for retirement portfolios. Suppose someone prepares to retire with a total sum of 1 million dollars in their portfolio; they'd like to be able to draw $30,000 a year from the portfolio for their living expenses. Would that work?

There's a large number of factors to take into account when analyzing this question, but for simplicity let's focus on just two: portfolio return and inflation. We can run a naive estimate, assuming average values: suppose an average yearly portfolio return of 4%, and average yearly inflation of 2% [4]

Let's denote our portfolio return as r=0.04, and inflation as q=0.02. Then the real return each year is:

\[r^{\mathrm{real}}=\frac{1+r}{1+q}-1\approx0.0196\]

Starting with $1,000,000, at the end of the year we'll have $1,019,600 and then draw $30,000 for living expenses [5], ending with $989,600. If we continue this way, the money runs out after ~55 years, which means that a person retiring at the age of 65 should be reasonably safe, right?

But this is very simplistic; assuming just average returns is risky, because they do a poor job of representing reality, and many factors have uncertainty. For example, the sequence of returns matters a lot; a bad year (-10%) followed by a great year (+18%) would still count as "4% on average" but produces significantly less money than two consecutive +4% years. Inflation is also unpredictable, and sometimes correlated with portfolio returns; there could be bad years of high inflation and low / volatile returns that can wreak havoc on a portfolio.

As we add factors (variance in yearly draws, mixed portfolios of stocks, bonds, real estate, life expectancy, unexpected events, changing tax laws etc.), relying on a single average estimate becomes increasingly more fraught. This is why Monte Carlo simulations are very popular in this domain: by drawing from reasonable distributions based on historical data, a Monte Carlo simulation can easily run a million different scenarios and provide estimates: for example, what are the odds of money running out before death.

Here's a useful chart from a simulation I ran:

In the top chart:

The dashed line shows the constant assumptions mentioned before: what happens when yearly return is always 4% and inflation is always 2%.
The shaded blue areas demonstrate the outcomes of 1,000,000 simulations where inflation and return numbers are drawn from reasonable normal distributions based on historical data. We see that in 25% of the cases, all money ran out by roughly 22 years.

In the bottom chart:

It's even easier to see how long the funds last; if we're interested in knowing, say, what are the odds that this plan will have enough money for 30 years - the chart shows it's about 60% (since in 40% of the simulations the funds were depleted at this point).

Looking a this simulation, under the current assumptions the plan sounds much riskier than the average assumption makes it appear. Assuming that a 65-y.o. person would plan for 25 years of retirement until death, the ~30% odds of not having sufficient funds for this duration of time are sobering. Perhaps a change in plans is needed (such as a more frugal lifestyle or securing additional funds in some way).

Retirement projection is only one of may ways in which Monte Carlo simulations are used for financial and economical applications; given the high uncertainty of these domains, it's very difficult to plan using analytical calculations. Company sales projections, growth projections, stock offering prices and much more uses Monte Carlo simulations to arrive at estimates with reasonable error bars.

Code

All the code for the explorations in this post is available on GitHub.

[1]

While these techniques were conceptually understood much earlier, it's not surprising that their first real applications coincided with the development of the first digital computers. As we'll see later in the post, Monte Carlo simulations benefit from running large numbers of trials to get reasonable accuracy.

[2]	The answer is: three empty red diamonds, three solid purple ovals, three striped green squiggles. Note that for each of the 4 SET attributes, these three cards are either all the same or all different.

[3]	For the same number of trials, this estimate has greater relative uncertainty than the 12-card estimate, because we encounter far fewer hands without a set.

[4]	For the examples in this post, we'll be using real dollars, or today's value of the money. We'll assume that the $30,000 yearly draw doesn't change and will instead apply inflation to the portfolio itself.

[5]	More sophisticated simulations would let us control these parameters; for example, are the expense funds drawn at the beginning or end of each year, or distributed on a monthly basis?

Summary of reading: July - September 2026

2026-09-30T17:05:00-07:00

"Wuthering Heights" by Emily Brontë - good writing, but the protagonists are quite something. There's barely a likable character in the whole story - they are all either crazy, evil, stupid or a combination of these. It's also remarkable to consider how small the world of some people in those times was (pre-Victorian England); whole lives spent within a tiny area, interacting with a handful of other people. And that's for members of the nobility, not some bumpkin on a remote farm.
"The Pickled City: The Story of New York pickles" by Paul Van Ravenstein and Monique Mulder - yes, a book about pickles, with lots of cool photos! That's what you get while randomly browsing in the local public library while your kids are picking out their books. All you wanted to know about the history of pickles, with particular emphasis on NYC in the early 20th century and the connection to Jewish immigration to that city.
"Rebel Island: the incredible history of Taiwan" by Jonathan Clements - a detailed history of Taiwan, focusing mostly on the early years (pre-modern times). The writing is a bit academic and dry, but it's very comprehensive. I wish the book spent more time on the modern era, though.
"The Old Ways: A Journey on Foot" by Robert Macfarlane - an ode to taking long walks (from hours to multi-day hikes). The writing is very beautiful, but much of it feels pointless; the author describes friends and obscure historical poets in a lot of detail which has little to do with the main thesis of the book. I'd say the first 1/4 of the good is very good, and the rest is unnecessary.
"After Rain" by William Trevor - I'm on a quest to find additional authors that write evocative short stories that stay with you, and Trevor was recommended as someone to try since my current favorites are Steinbeck and Lahiri. "After Rain" is a collection of 12 stories, mostly about love and relationships - between partners but also parents and children; the settings are usually in 1980s Ireland or England. Very good stories overall, I'm glad I tried this author and should try something else by him.
"Spoonhandle" by Ruth Moore - a novel about the lives of Maine natives living through changing times, seeing their picturesque coastal area become a magnet for "summer people" - rich vacationers from the big cities. There's no single protagonist in the book and the writing is very warm and human. I enjoyed it.
"The Almighty Dollar: 500 Years of the World's Most Powerful Money" by Brendan Greeley - a history-oriented exploration of the importance of the dollar in world economy, from its beginnings as the Joachimsthaler in 16th century Bohemia and until modern times. This book has some interesting ideas but I found it to be quite academical and not easy to digest.
"Behave: The Biology of Humans at Our Best and Worst" by Robert M. Sapolsky - the scope of this (very long) book is huge: explaining how human behavior works starting from basic biological principles. Our understanding of the actual biological underpinnings of behavior is so nascent, however, that the vast majority of the book focuses on the usual things - the psychological studies (including all the well known ones like prisoner's dilemma, the Stanford prison experiment, etc) and comparison with other animals. Having read several other books on this topic in the past, I haven't found much new in this one.
"Humorous short stories" by Anton Chekhov (read in Russian) - a collection of short stories with a humorous angle. The stories are typically very short, just 2-3 pages each, but the character building is superb. I should try reading something longer by Chekhov.
"Plainsong" by Kent Haruf - describes a year in the life of several loosely connected people in a remote farm town in eastern Colorado. A wholesome story, without a specific beginning or end - just the lives of ordinary folks.
"Behind the Beautiful Forevers: Life, Death, and Hope in a Mumbai Undercity" by Katherine Boo - investigative journalism of the lives of several connected families in Annawadi - one of Mumbai's most miserable slums - roughly from 2007 to 2010. The only thing more shocking than the squalor in which these people live is the degree of corruption, fraud and exploitation surrounding them. A very poignant read.
"The Dispossessed" by Ursula K. Le Guin - science fiction only in the sense of 1984; it's always amusing to examine how unambitious the futuristic vision of some authors is. The book was written in the early 1970s, and it presents life thousands of years in the future as one that's pretty much identical to what it was in the early 1970s, except for close-to-speed-of-light travel. But science is not the point of this book anyway; it's more of a political or philosophical fiction, examining an anarchist-communist society and comparing it to a capitalist, "propertarian" one. Overall, this book gave me very strong vibes of the famous Soviet kids story "Dunno on the moon" (Незнайка на Луне), which I recall reading as a kid.
"The Lion Women of Tehran" by Marjan Kamali - a multi-generational novel focusing on the life of young women in Iran from the early 1950s to the early 1980s, in the context of the various authoritarian governments of those eras (from Mohammad Reza Pahlavi's rule to after the Islamic revolution of 1979). Nice book showing that a good story can be told without too many unnecessary twists and classical attention-grabbing dramas.

Re-reads:

"The Old Way: A story of the first people" by Elizabeth Marshall Thomas
"Interpreter of Maladies" by Jhumpa Lahiri
"Norwegian Wood" by Haruki Murakami
"Cannery Row" by John Steinbeck

Rusty thoughts on "Parse, don't validate"

2026-09-26T08:25:00-07:00

Like many programmers, I find Alexis King's Parse, don't validate article fascinating, because it gives a name to an idiom that seems familiar and important - one I've observed and used in the past without naming it explicitly. This post is a review of the "Parse, don't validate" pattern applied to the Rust programming language (the original post uses Haskell). I was particularly interested in finding educational examples of this pattern in the Rust standard library and other well-known projects.

Without repeating the original article (please read it first!), here's the gist of it.

Consider the venerable Vec; its first method returns Option<&T>. Why? Because a vector is not guaranteed to have any elements in it, so what to do if first is invoked on an empty one? Returning an Option in this case is idiomatic in Rust [1], with convenient syntax sugar for accepting the result of functions that return Option and deciding what to do next.

So what's the issue?

Imagine we have a function to read some configuration paths from an env var, while enforcing the invariant that the list can't be empty:

use anyhow::{Result, ensure};

fn get_configuration_directories() -> Result<Vec<PathBuf>> {
    let value = env::var("CONFIG_DIRS").context("could not read CONFIG_DIRS")?;

    let directories: Vec<PathBuf> = value
        .split(',')
        .map(str::trim)
        .map(PathBuf::from)
        .collect();

    ensure!(!directories.is_empty(), "empty CONFIG_DIRS");
    Ok(directories)
}

So far, so good. Now let's take a typical usage of this function:

fn main() -> Result<()> {
    let config_dirs = get_configuration_directories()?;

    match config_dirs.first() {
        Some(cache_dir) => initialize_cache(cache_dir),
        None => unreachable!("already checked that CONFIG_DIRS is non-empty"),
    }

    Ok(())
}

Once get_configuration_directories returns a successful result, we are guaranteed that the vector isn't empty. And yet, if we want to get the first element of this vector, we have to use the first method that returns Option<&T>. We are therefore forced - again - to handle a potentially empty case (where the option is None).

As the original article states, this has a number of problems with code clarity, potential performance implications and a ticking time bomb if the invariant is ever changed in get_configuration_directories.

The core issue is that Vec is fundamentally a type that can be empty; we can carry along a "This one can't be empty, pinky promise!" comment on all the relevant code, but it's not formally checked by anything.

A type for "non-empty" vector

The solution is leveraging the type system to enforce a newly established invariant. We can use a separate type for "a vector that cannot be empty"; in fact, such types already exist in several Rust crates - for example nonempty:

pub struct NonEmpty<T> {
    pub head: T,
    pub tail: Vec<T>,
}

This type has no constructor that permits "no elements"; its new takes one element, and its first method returns &T without an Option:

pub const fn new(e: T) -> Self {
    Self::singleton(e)
}

pub const fn singleton(head: T) -> Self {
    NonEmpty {
        head,
        tail: Vec::new(),
    }
}

pub const fn first(&self) -> &T {
    &self.head
}

The rest of the crate deals with making NonEmpty behave as close as possible to a normal Vec, by implementing many useful traits, as well as conversions like:

pub fn from_vec(mut vec: Vec<T>) -> Option<NonEmpty<T>> {
    if vec.is_empty() {
        None
    } else {
        let head = vec.remove(0);
        Some(NonEmpty { head, tail: vec })
    }
}

Let's see how our get_configuration_directories function would look if it returned a NonEmpty instead of a plain Vec:

fn get_configuration_directories() -> Result<NonEmpty<PathBuf>> {
    let value = env::var("CONFIG_DIRS").context("could not read CONFIG_DIRS")?;

    let directories = value
        .split(',')
        .map(str::trim)
        .map(PathBuf::from)
        .collect();

    let Some(directories) = NonEmpty::from_vec(directories) else {
        bail!("CONFIG_DIRS cannot be empty");
    };

    Ok(directories)
}

Note the use of NonEmpty::from_vec here - this is where the invariant is established. Now a successful result is NonEmpty, not just Vec. The client code looks like:

fn main() -> Result<()> {
    let config_dirs = get_configuration_directories()?;

    initialize_cache(config_dirs.first())?;
    Ok(())
}

There's no need to check if the returned value is empty again; this is enforced by the type system!

This is where the parse vs. validate terminology of the original article comes from. When get_configuration_directories returned a Vec, it simply validated it. But when it returns a NonEmpty - the vector is transformed into another entity which carries additional meaning. If we treat the concept of parsing in the most generic sense - "transforming data from one format to another", this fits.

To mention a less artificial example, the Rust rewrite of core POSIX utilities uses NonEmpty in several places [2]. For example, when constructing a shell pipeline:

pub struct Pipeline {
    pub commands: NonEmpty<Command>,
    pub negate_status: bool,
}

The command parser's code:

fn parse_pipeline(&mut self, alias_table: &AliasTable) -> ParseResult<Option<Pipeline>> {
    // pipeline = "!" command ("|" linebreak command)*
    let negate_status = self.match_alternatives(&[CommandToken::Bang])?.is_some();
    let mut commands = if let Some(command) = self.parse_command(alias_table)? {
        NonEmpty::new(command)
    } else {
        return Ok(None);
    };

    // ...

A valid Pipeline is only returned if there are some commands in the parsed AST. Otherwise, it just returns None. Once this is done, the client code can use commands.first() without having to worry about the possibility of it returning None.

Non-zero integers

Rust has a generic type called NonZero, to describe unsigned numeric quantities that are known to be non-zero.

For example, thread::available_parallelism is defined as:

pub fn available_parallelism() -> Result<NonZero<usize>>

If the call is successful, it returns a NonZero<usize>, which is like a normal usize with the restriction that it's not zero. Client code doesn't have to keep checking whether the parallelism is 0 - it's enshrined in the type system.

Rust defines the division operator with NonZero<usize> in the denominator as an operation that "cannot panic".

NonZero has an additional advantage: zero is an invalid value for the type, so Rust can use the zero bit pattern to represent None. Consequently, Option<NonZeroUsize> is guaranteed to have the same size and alignment as NonZeroUsize itself (and as usize). This can avoid the extra storage that an Option<usize> would generally require.

Parsing JSON

A common example of the "parse, don't validate" idiom appears in deserializing data from a JSON string. Rust's serde crate enables us to do the parsing, with validated decisions encoded into the type system, e.g.:

#[derive(Debug, Deserialize)]
struct Config {
    name: String,
    workers: NonZeroUsize,
    mode: Mode,
}

#[derive(Debug, Deserialize)]
#[serde(rename_all = "snake_case")]
enum Mode {
    Fast,
    Safe,
}

And then later:

let input = r#"
     {
         "name": "compiler",
         "workers": 4,
         "mode": "fast"
     }
 "#;

 let config: Config = serde_json::from_str(input)?;

There is a lot happening behind the scenes:

The types of all fields are enforced (e.g. "name" cannot be an array).
The mode is validated to be one of the enum values of Mode.
workers is validated to be a non-zero integer, because of the NonZeroUsize field type.

We take code like this for granted these days, but it's still a great example of the pattern discussed in this post. Once the parser converted mode into the Mode enum, no further validation is required.

In dynamic languages like Python and JS, the process is usually much more manual. Python's json.loads gives us a dictionary, and it's up to the user to validate its contents. Libraries like Pydantic permit an approach closer to Rust's, but they're not universally used.

[1]

Other languages - like Go or Python - don't have a dedicated first method or function, and instead have the programmer rely on a runtime check that raises some sort of exception or panic when lst[0] is accessed on an empty list or slice.

This is not by mistakel these methods / functions were not added to Go and Python precisely because their semantics are tricky to narrow down without a type like Option.

[2]	The project implements its own `NonEmpty`, without relying on the `nonempty` crate, but all the points in this post apply.

Notes on discrete-time Fourier series and transform

2026-09-19T09:39:00-07:00

The following are my notes on discrete-time Fourier series (DTFS), as well as the discrete-time Fourier transform (DTFT). These topics serve as an important theoretical underpinning to the digital processing of signals by computers using the DFT (which will be covered in a future post).

For discrete-time signals, we use the square bracket notation to denote samples: x[n] is the n-th sample of signal x.

We’ll say that a discrete-time (just discrete from this point on, for brevity) signal is periodic with period N if:

\[x[n]=x[n+N]\qquad\forall{n}\]

When talking about continuous-time Fourier series, we started with real trigonometric functions and later moved to complex exponentials. Here, we’ll just go ahead and start directly with (discrete) complex exponentials - converting between the two representations isn’t difficult and can be done as needed (Appendix C in this post demonstrates it).

We’ll be dealing with the following family of signals [1]:

\[\phi_k[n]=e^{ikw_0n}=e^{ik(2\pi/N)n}\qquad \forall k\in\mathbb{Z}\]

These are discrete complex exponentials that are periodic with period N. w_0=\frac{2\pi}{N} is the angular frequency. As discussed in Appendix A, there are only N such distinct signals, since:

\[\phi_k[n]=\phi_{k+N}[n]\qquad\forall{k}\]

Coefficients of discrete-time Fourier series

We’ll want to consider the representation of an arbitrary N-periodic x[n] by a linear combination of \phi_k[n]:

\[x[n]=\sum_{k=\langle N\rangle}a_k\phi_k[n]\]

The notation k=\langle N\rangle means k runs over any sequence of N consecutive integers. Since there are only N distinct signals \phi_k, the order of summation doesn’t matter - as long as all of them are used. So it could go from 0 to N-1, or from 1 to N, or from 2 to N+1, and so on. All of these orders will enumerate all distinct \phi_k.

Going back to our linear combination:

\[x[n]=\sum_{k=\langle N\rangle}a_k\phi_k[n]=\sum_{k=\langle N\rangle}a_ke^{ikw_0n}\]

This is the discrete-time Fourier series (DTFS) representation of x[n] and the coefficients a_k are the Fourier series coefficients. Note that there are no convergence issues as in the continuous case, because we’re dealing with a finite sum.

To find the coefficients a_k, we’ll proceed in a way that’s somewhat similar to the continuous case. Multiply both sides of the equation above by e^{-ir(2\pi/N)n} and sum over N terms:

\[\sum_{n=\langle N\rangle}x[n]e^{-ir(2\pi/N)n}= \sum_{n=\langle N\rangle}\sum_{k=\langle N\rangle}a_ke^{i(k-r)(2\pi/N)n}\]

Interchanging the order of summation on the right-hand side:

\[\sum_{n=\langle N\rangle}x[n]e^{-ir(2\pi/N)n}= \sum_{k=\langle N\rangle}a_k\sum_{n=\langle N\rangle}e^{i(k-r)(2\pi/N)n}\]

By Appendix B, the inner sum on the right-hand side is equal to N when N\mid r-k and to 0 otherwise. Without loss of generality - since we iterate over N consecutive values - let’s just assume this happens in our sum when r=k (and therefore r-k=0, which is an integer multiple of N).

So our equation becomes:

\[\sum_{n=\langle N\rangle}x[n]e^{-ir(2\pi/N)n}=a_r N\]

Or:

\[a_r=\frac{1}{N}\sum_{n=\langle N\rangle}x[n]e^{-ir(2\pi/N)n}\]

To conclude, the Fourier decomposition of a periodic discrete signal x[n] into periodic complex exponentials is:

\[\boxed{ \begin{aligned} x[n]&=\sum_{k=\langle N\rangle}a_k e^{ik(2\pi/N)n}\\ \vspace{2pt}\\ a_k&=\frac{1}{N}\sum_{n=\langle N\rangle}x[n]e^{-ik(2\pi/N)n} \end{aligned} }\]

Example: revisiting the triangular function

Let’s revisit the triangular function t(x) from the post on Fourier series. Here, we’ll be using the following sampled version:

The sampling is done with 12 samples per period of 4. The spacing between two samples is \Delta x=\frac{1}{3}. Another way to express it:

\[t_d[n]=t(n/3)\]

Substituting this into the definition of t(x), we get:

\[t_d[n]= \begin{cases} \frac{n}{3} & 0 \leq n \leq 3 \\ 2-\frac{n}{3} & 3 < n \leq 6 \\ \end{cases}\]

We then do an odd extension [2] to the range n=[6..12), and repeat it with period N=12.

Using considerations similar to the continuous time case, because our function is odd we just need a sine series here:

\[t_d[n]=\sum_{k=1}^{5}b_k sin\frac{k\pi n}{6}\]

The five terms correspond to the five pairs of nonzero-frequency indices; see Appendix C.

And the coefficients are:

\[b_k=\frac{1}{6}\sum_{n=0}^{11}t_d[n]sin\frac{k\pi n}{6}\\\]

We can run this sum over half the period, because the terms at and 12-n are equal (both t_d[n] and the sine change signs). Also, the terms at n=0 and n=6 vanish because t_d[0]=t_d[6]=0:

\[b_k=\frac{1}{3}\sum_{n=1}^{5}t_d[n]sin\frac{k\pi n}{6}\\\]

Looking at the values of t_d:

\[b_1=\frac{1}{3}\left(\frac{1}{3}sin\frac{\pi}{6}+\frac{2}{3}sin\frac{2\pi}{6}+sin\frac{3\pi}{6}+\frac{2}{3}sin\frac{4\pi}{6}+\frac{1}{3}sin\frac{\pi}{6}\right)=\frac{4+2\sqrt{3}}{9}\]

We can similarly calculate all b_k until k=5:

\[\begin{aligned} b_1&=\frac{4+2\sqrt{3}}{9}\\ b_2&=0\\ b_3&=-\frac{1}{9}\\ b_4&=0\\ b_5&=\frac{4-2\sqrt{3}}{9}\\ \end{aligned}\]

The resulting Fourier series is therefore:

\[t_d[n]=\frac{4+2\sqrt{3}}{9}sin\frac{\pi n}{6}-\frac{1}{9}sin\frac{\pi n}{2}+\frac{4-2\sqrt{3}}{9}sin\frac{5\pi n}{6}\]

You can use the following interactive plot to explore this approximation to the discrete triangle function t_d[n]:

Terms in the Fourier series

The dropdown box selects how many of the terms in the series shown above to plot. The orange crosses show where the series’ values at the integer indices are, and the orange line is interpolated to better visualize the sine waves being added. Note that when all the coefficients are used, the DTFS exactly reconstructs the input signal.

Discrete-time Fourier Transform (DTFT)

In the post on the Fourier Transform we’ve seen how non-periodic signals can be represented in the frequency domain by taking the Fourier series and calculating them at the limit L\rightarrow\infty. Here, we’ll do something similar for discrete signals.

Suppose we have a non-periodic discrete signal x[n] of finite duration (it’s zero outside a finite range of indices). Here’s a sample plot:

The bottom half is x_N[n] - a periodic signal, one period (N) of which is equal to x[n]. We choose M large enough that x[n]=0 outside [-M,M] and set N=2M+1.

Since x_N is periodic, we can represent it with DTFS [3]:

\[\begin{aligned} x_N[n]&=\sum_{k=-M}^{M}a_k e^{ik(2\pi/N)n}\\ a_k&=\frac{1}{N}\sum_{n=-M}^{M}x[n]e^{-ik(2\pi/N)n} \end{aligned}\]

We’ll use the notation of angular frequency:

\[\Delta w=\frac{2\pi}{N}\qquad w_k=k\Delta w\]

Then:

\[\begin{aligned} x_N[n]&=\sum_{k=-M}^{M}a_k e^{i w_k n}\\ a_k&=\frac{1}{N}\sum_{n=-M}^{M}x[n]e^{-i w_k n} \end{aligned}\]

Now the critical part; since x[n] is zero outside [-M,M], extending the summation to all integers does not change its value:

\[a_k=\frac{1}{N}\sum_{n=-\infty}^{\infty}x[n]e^{-i w_k n}\]

We’ll define the following continuous function with variable :

\[X(w)=\sum_{n=-\infty}^{\infty}x[n]e^{-i w n}\]

Then at points w_k:

\[a_k=\frac{1}{N}X(w_k)\]

The function X(w) is the discrete-time Fourier transform (DTFT) of x[n]. The inverse DTFT process reconstructs x[n] from X(w), by substituting a_k back into the DTFS:

\[x_N[n]=\sum_{k=-M}^{M}\frac{1}{N}X(w_k) e^{i w_k n}\]

By our definition of \Delta w, we have:

\[\frac{1}{N}=\frac{\Delta w}{2\pi}\]

Therefore:

\[x_N[n]=\frac{1}{2\pi}\sum_{k=-M}^{M}X(w_k) e^{i w_k n}\Delta w\]

Just like in the continuous case, we recognize the sum as a Riemann sum. Let M\rightarrow\infty. The spacing \Delta w tends to zero, and we can rewrite the sum as an integral:

\[x[n]=\lim_{M\rightarrow\infty}x_N[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(w)e^{iwn}dw\]

To conclude, the DTFT and inverse DTFT pair are:

\[\boxed{ \begin{aligned} X(w)&=\sum_{n=-\infty}^{\infty}x[n]e^{-i w n}\\ \vspace{2pt}\\ x[n]&=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(w)e^{iwn}dw \end{aligned} }\]

X(w) is a continuous, periodic function with period 2\pi, because:

\[\begin{aligned} X(w+2\pi)&=\sum_{n=-\infty}^{\infty}x[n]e^{-i (w+2\pi) n}\\ &=\sum_{n=-\infty}^{\infty}x[n]e^{-i w n}e^{-i 2\pi n}\\ &=\sum_{n=-\infty}^{\infty}x[n]e^{-i w n}=X(w)\\ \end{aligned}\]

Therefore, it’s enough to run the integral on any interval of length 2\pi.

The DTFS and DTFT are closely related. Starting with a finite-duration signal x[n], we form a periodic signal by repeating a block of N samples containing it. Its DTFS coefficients are equally spaced, scaled samples of the DTFT:

\[a_k=\frac{1}{N}X\left(\frac{2\pi k}{N}\right)\]

The relation is similar to the one between Fourier series and the Fourier transform in the continuous case.

The DTFT also has useful properties similar to the continuous Fourier transform: linearity, time shifting, the convolution theorem etc, but we won’t spend time on them here.

Appendix A: Discrete-time complex exponentials

Discrete-time complex exponentials have the form e^{iw_0n} where w_0 is the angular frequency. The discrete nature of these functions leads to some interesting outcomes; for example, consider the exponential with angular frequency w_0+2\pi:

\[e^{i(w_0+2\pi)n}=e^{iw_0n}e^{i2\pi n}=e^{iw_0n}\]

Therefore, the exponential at w_0+2\pi is exactly the same as the exponential at w_0. In considering discrete complex exponentials, we need only choose w_0 from a range of length 2\pi.

Another interesting aspect of these functions concerns their periodicity. In order for e^{iw_0n} to be periodic with some positive integer period N, the following must hold:

\[e^{iw_0(n+N)}=e^{iw_0n}\]

In other words:

\[e^{iw_0N}=1\]

So w_0N itself must be an integral multiple of 2\pi. For some integer k:

\[w_0N=k\cdot(2\pi)\Longrightarrow w_0=\frac{k2\pi}{N}\]

For a fixed positive integer N, let’s consider all complex exponentials for which N is a period:

\[\phi_k[n]=e^{ik(2\pi/N)n}\qquad \forall k\in\mathbb{Z}\]

But earlier we’ve said that all complex exponentials with angular frequencies 2\pi apart are identical. This means that the set above consists of only N distinct exponentials. We can choose any starting k and get distinct signals with k, k+1, k+2 and so on until k+N-1. After that, they begin repeating: \phi_N[n]=\phi_0[n] and so on.

All of this is very different from the continuous case, where t is a real number. The signals e^{iw_0t} and e^{i(w_0+2\pi)t} agree only at integer values of t, so they are distinct functions of t\in\mathbb{R}. Consequently, for a fixed positive period , the continuous-time exponentials e^{ik(2\pi/T)t} are distinct for every integer k, giving infinitely many distinct signals.

Appendix B: Sum of consecutive complex exponentials

Let’s consider the following discrete function, defined by a sum:

\[A[k]=\sum_{n=0}^{N-1}e^{i(2\pi/N)kn}\]

We’ll want to consider two cases:

(1) When k is an integer multiple of N: k=rN for some integer r (for example k=0, k=N, k=-2N etc.)

In this case, the sum becomes:

\[A[k]=\sum_{n=0}^{N-1}e^{i(2\pi/N)rNn}=\sum_{n=0}^{N-1}e^{i2\pi rn}\]

Since both r and are integers, the exponent is an integer multiple of 2\pi i, so:

\[A[k]=\sum_{n=0}^{N-1}1=N\]

(2) When k is not an integer multiple of N, let’s do a variable change p=n+1:

\[A[k]=\sum_{p=1}^{N}e^{i(2\pi/N)k(p-1)}\]

We can now use the finite sum formula for a geometric progression with the common ratio r=e^{i(2\pi/N)k}; the sum is:

\[S_N=\frac{1-r^N}{1-r}\]

Note that r\neq 1 because k is not an integer multiple of N. However:

\[r^N=e^{i(2\pi/N)kN}=1\]

Therefore:

\[A[k]=S_N=0\]

While these calculations demonstrate the sums from n=0 to N-1, they work exactly the same for any N consecutive indices, because the summand is periodic in with period N. Therefore:

\[A[k]=\sum_{n=\langle N\rangle}e^{i(2\pi/N)kn}= \begin{cases} N & N\mid k\\ 0 & \text{otherwise} \end{cases}\]

Appendix C: Sine series decomposition of t_d[n]

Using the formula for DTFS developed earlier in the post:

\[\begin{aligned} a_k&=\frac{1}{12}\sum_{n=0}^{11}t_d[n]e^{-ik(\pi/6)n}\\ &=\frac{1}{12}\sum_{n=0}^{11}t_d[n]cos\frac{k\pi n}{6}-i\frac{1}{12}\sum_{n=0}^{11}t_d[n]sin\frac{k\pi n}{6}\\ \end{aligned}\]

Recall that t_d[n] is odd over the full range n=[0..12). Specifically, since the cosine is an even function, we have for any m:

\[\begin{aligned} t_d[12-m]&=-t_d[m]\\ cos\frac{k\pi(12-m)}{6}&=cos(2\pi k -\frac{k\pi m}{6})=cos\frac{k\pi m}{6}\\ \end{aligned}\]

Therefore the cosine terms cancel out throughout the period; term 1 is the same as negative term 11, etc.

We’re left with the sines:

\[a_k=-i\frac{1}{12}\sum_{n=0}^{11}t_d[n]sin\frac{k\pi n}{6}\]

All coefficients are purely imaginary. For convenience, let’s write:

\[\begin{aligned} b_k&=\frac{1}{6}\sum_{n=0}^{11}t_d[n]sin\frac{k\pi n}{6}\\ a_k&=\frac{-i b_k}{2}\\ \end{aligned}\]

Let’s turn back to the Fourier reconstruction formula:

\[t_d[n]=\sum_{k=\langle N\rangle}a_k e^{ik(2\pi/N)n}\]

And pair up exponentials at k with ones at 12-k:

\[e^{i(12-k)\pi n/6}=e^{i2\pi n}e^{-ikn\pi/6}=e^{-ikn\pi/6}\]

Also, since t_d[n] are real, the conjugate of a_k is:

\[a_{12-k}=a_k^*=\frac{ib_k}{2}\]

So adding the paired exponentials with their coefficients and using Euler’s formula:

\[\begin{aligned} a_ke^{i k\pi n/6}+a_{12-k}e^{-i k\pi n/6}&=-\frac{ib_k}{2}e^{i k\pi n/6}+\frac{ib_k}{2}e^{-i k\pi n/6}\\ &=b_k\frac{e^{i k\pi n/6}-e^{-i k\pi n/6}}{2i}\\ &=b_k sin\frac{k\pi n}{6} \end{aligned}\]

Therefore:

\[t_d[n]=\sum_{k=0}^{11}a_k e^{ik(2\pi/N)n}=\sum_{k=1}^{5}b_k sin\frac{k\pi n}{6}\]

We don’t include the terms at k=0 and k=6, because a_0=a_6=0 (sine of an integral multiple of \pi).

[1]	There are a lot of symbols used in this equation. is the "time" index of the signal x[n] - its domain; k enumerates the k-th function \phi_k from a family of functions - each \phi_k[n] is a function of ; N is the period; is the imaginary unit.

[2]	The odd extension may be slightly confusing because we don’t use negative indices. But indices are arbitrary! In our case, within a single period index -n is equivalent to 12-n; if we consider one period around the halfway point of n=6, it’s easy to see that the function is indeed odd.

[3]	Note that we use x[n], not x_N[n] in the summation for a_k. This is because within the chosen range they are equal.

How big are factorials?

2026-08-27T18:52:00-07:00

The other day, I found myself wondering how big 52! (52 factorial) is, and that led me to ponder how these could be estimated without a calculator or a computer.

It turns out there’s some fairly interesting math behind being able to estimate the size (number of digits) of a factorial reasonably accurately. This post will start by stating how to do the estimate, and if you’re curious you can read on for the math background.

Without further ado, the approximation is:

\[\text{number of digits in n!}\approx n\log_{10}\left(\frac{n}{e}\right)+2\]

As an example, let’s use my original question, by estimating this for 52!

Well, 52 divided by e is... 20-ish? And \log_{10}(20) is about 1.3 [1]; therefore our estimate comes out to:

\[\text{number of digits in 52!}\approx 52\cdot 1.3 +2 \approx69\]

The real answer is 68, so this is very close! In estimates like this - when you’re dealing with enormous numbers - being off by a couple of digits usually isn't a big deal.

The Gamma function

The Gamma function for real n>0 is defined [2] as:

\[\Gamma(n)=\int_{0}^{\infty}x^{n-1}e^{-x}dx\]

This integral does not have an analytic expression in the general case, but it does have a very useful property that we can take advantage of. Let’s see what \Gamma(n+1) is:

\[\Gamma(n+1)=\int_{0}^{\infty}x^{n}e^{-x}dx\]

And now use integration by parts with:

\[u=x^n\qquad v=-e^{-x}\]

Then:

\[du=nx^{n-1} dx\qquad dv=e^{-x}dx\]

So:

\[\begin{aligned} \Gamma(n+1)&=\int_{0}^{\infty}x^{n}e^{-x}dx\\ &=\left. -x^n e^{-x}\right|_{0}^{\infty}-\int_{0}^{\infty}-e^{-x}n x^{n-1}dx\\ &=n\int_{0}^{\infty}x^{n-1}e^{-x}dx \end{aligned}\]

But notice that the last integral is just \Gamma(n); therefore, we’ve shown that:

\[\Gamma(n+1)=n\Gamma(n)\]

Let’s also calculate \Gamma(1) - it’s a special case that has an analytical solution:

\[\Gamma(1)=\int_{0}^{\infty}e^{-x}dx=\left. -e^{-x}\right|_{0}^{\infty}=1\]

This helps establish an induction argument:

\[\begin{aligned} \Gamma(2)=1\cdot\Gamma(1)&=1!\\ \Gamma(3)=2\cdot\Gamma(2)&=2!\\ \Gamma(4)=3\cdot\Gamma(3)&=3!\\ \dots\\ \Gamma(n+1)=n\cdot\Gamma(n)&=n! \end{aligned}\]

In other words - the Gamma function is an interpolation of the factorial over all positive reals. Here’s a plot of the Gamma function over a small range; note that the y axis is log-scale because of the function’s fast growth:

Stirling’s approximation

You may have encountered Stirling’s approximation before:

\[n!\approx\sqrt{2\pi n}\cdot\left(\frac{n}{e} \right)^n\]

It’s a great approximation that works reasonably well even for small . This section is a brief overview of how Stirling’s formula is derived from the Gamma function.

Taking:

\[n!=\Gamma(n+1)=\int_{0}^{\infty}x^n e^{-x}dx\]

We’ll start by massaging the integrand a bit:

\[n!=\int_{0}^{\infty}x^n e^{-x}dx=\int_{0}^{\infty}e^{n \ln x} e^{-x}dx\]

And making a change of variables x=ny, which means that dx=ndy:

\[\begin{aligned} n!&=\int_{0}^{\infty}ne^{n\ln(ny) - ny}dy=n\int_{0}^{\infty}e^{n(\ln n+\ln y-y)}dy\\ &=ne^{n\ln n}\int_{0}^{\infty}e^{n(\ln y-y)}dy \end{aligned}\]

These steps make the integral amenable to applying Laplace’s method, which allows us to approximate definite integrals of the form:

\[\int_{a}^{b}e^{nf(x)}dx\]

Where is a twice-differentiable function and some large number. By Laplace’s method, such integrals can be approximated by:

\[\int_{a}^{b}e^{nf(x)}dx\approx\sqrt{\frac{2\pi}{n|f''(x_0)|}}e^{n f(x_0)}\]

Where is the global maximum of .

Let’s see how to apply this method [3] to the latest equation we have for n! (renaming the dummy integration variable back to x):

\[n!=ne^{n\ln n}\int_{0}^{\infty}e^{n(\ln x-x)}dx\]

In our case, f(x)=\ln x - x. It’s easy to show that this function is twice differentiable and has a global maximum at . Moreover:

\[\begin{aligned} f(x_0)&=-1\\ f''(x_0)&=-1\\ \end{aligned}\]

Substituting these into the proper places in Laplace’s approximation, we get:

\[\begin{aligned} n! &\approx ne^{n\ln n}\sqrt{\frac{2\pi}{n}}e^{-n}\\ &\approx \sqrt{2\pi n}\cdot e^{n(\ln n - 1)}\\ &\approx \sqrt{2\pi n}\cdot \left(\frac{n}{e}\right)^n\quad\blacksquare \end{aligned}\]

Number of digits from Stirling’s approximation

We can calculate the number of digits in n! by taking the base-10 logarithm of Stirling’s formula:

\[\begin{aligned} \text{number of digits in n!}&\approx \log_{10} \left(\sqrt{2\pi n}\cdot \left(\frac{n}{e}\right)^n\right)\\ &\approx \log_{10}\left(\sqrt{2\pi n}\right)+ \log_{10}\left(\frac{n}{e}\right)^n\\ &\approx \log_{10}\left(\sqrt{2\pi n}\right)+ n\cdot \log_{10}\left(\frac{n}{e}\right) \end{aligned}\]

Note that the first term is not multiplied by itself; therefore, as grows, it will become less and less noticeable. That said, it still adds a couple of digits - so you should take it into account if you want a more accurate approximation [4]

[1]	Mental tricks for calculating \log_{10} is a different topic, but it really helps to remember that \log_{10}2=0.3, \log_{10}3=0.5, and from here using the various logarithm laws to estimate multiples.

[2]	As a matter of fact, the Gamma function is defined for complex numbers, but for our purpose talking about the reals is sufficient.

[3]	As much as I’d like to dive into why this approximation works, the rabbit hole in this post is already deep enough!

[4]	It adds up to 2 extra digits as long as is less than 1600 or so, and may add more than 2 after that, though no more than 3 until n is 160000. It’s not clear why anyone would like to estimate the number of digits of 1600! (about 4450, in case you were wondering), let alone 160000!

Eli Bendersky's website

Monte-Carlo simulations

Estimating pi

Combinatorial simulation - the game of SET

Retirement projection

Code

Summary of reading: July - September 2026

Rusty thoughts on "Parse, don't validate"

A type for "non-empty" vector

Gradual parsing and type refinement

Non-zero integers

Parsing JSON

Notes on discrete-time Fourier series and transform

Coefficients of discrete-time Fourier series

Example: revisiting the triangular function

Discrete-time Fourier Transform (DTFT)

Appendix A: Discrete-time complex exponentials

Appendix B: Sum of consecutive complex exponentials

Appendix C: Sine series decomposition of t_d[n]

How big are factorials?

The Gamma function

Stirling’s approximation

Number of digits from Stirling’s approximation