In 1973, oil-exporting Arab countries embargoed shipments to countries that had supported Israel in the Yom Kippur War of that same year. This included the United States. Congress, in an attempt to reduce fuel consumption at home, signed the first Corporate Average Fuel Economy (CAFE) standards into law as a result.

These standards stuck around and evolved to focus on global warming and emissions instead. Around 2007, Congress rewrote the legislation to deal with some (legitimate) concerns. They created different efficiency targets for different vehicle "footprints" (base size, length times width) with larger vehicles having lower mileage per gallon requirements.

While the legislation was intended to reduce the amount of fuel consumed, it actually just incentivized manufacturers to build larger vehicles and stop selling small ones. The CAFE standards killed the station wagon, birthed the modern SUV and luxury pickup truck, and failed to get us to the goal of 50 MPG average efficiency by 2025.

There's a similar story about unintended side effects from the British colonial administration of India. The government was concerned about the number of cobras in Delhi, and offered a bounty for each dead cobra brought in by hunters. The locals hunted cobras and were rewarded.

However, some of them started breeding cobras to take in and claim bounties with. As soon as the British got wind of the scam, they canceled the program. The scammers released their captive cobras, and the city was worse off than before. This idea of regulations creating "perverse incentives" is even called the "cobra effect."

These are both famous examples of policy having unintended side effects. But why did CAFE fail to achieve the (ambitious) goals Congress set in 2007? Was it doomed to fail from the start?

CAFE's story isn't unique. A couple of conversations I've had in the past six months got me thinking: regulations that suffer from the cobra effect have something in common. They're not written like engineering requirements, yet engineers have to follow them. In this post I'll touch on a couple other examples, and hopefully leave you wondering whether engineers could have written them better.

How Lawyers Screwed Up

GDPR and the Cost of Compliance

Since 2018 the European Union has required all organizations which handle EU citizens' data to comply with the General Data Protection Regulation (GDPR). This legislation was promoted as a way to give consumers power over how their data was used, and break the stranglehold of data-collecting/advertising monopolies like Google.

Well, it didn't. In fact, it actually made the advertising space less competitive. Big companies like Google easily absorbed the cost of GDPR compliance, while smaller companies and startups had their margins squeezed out by the extra burden. Companies like Microsoft, which is famously a law firm with a software engineering division attached, were already dealing with large-scale compliance tracking. Smaller firms were not, and were more impacted than the big ones.

GDPR and, to a lesser extent, other regulations like it (in the UK, California, Texas, and some other states) are more examples of a well-intentioned piece of legislation written by lawyers having the opposite effect they wanted. It fortified the position of the monopolies instead of breaking it.

To get a feel for how the GDPR reads, let's look at how Article 32 (Security of processing) opens:

Taking into account the state of the art, the costs of implementation and the nature, scope, context and purposes of processing as well as the risk of varying likelihood and severity for the rights and freedoms of natural persons, the controller and the processor shall implement appropriate technical and organisational measures to ensure a level of security appropriate to the risk, including inter alia as appropriate:
(1) the pseudonymisation and encryption of personal data;
(2) the ability to ensure the ongoing confidentiality, integrity, availability and resilience of processing systems and services;
(3) the ability to restore the availability and access to personal data in a timely manner in the event of a physical or technical incident;
(4) a process for regularly testing, assessing and evaluating the effectiveness of technical and organisational measures for ensuring the security of the processing.

What is a "level of security appropriate to the risk?" What does it mean to take into account "costs of implementation?" How available and resilient should processing systems be? What does "timely manner" mean? All of these questions are (intentionally) left to the courts to decide, but that creates friction and uncertainty.

Remember, GDPR was brought to you by the same geniuses that created cookie pop-ups. We'll come back to this later.

America's EAR, FDPR, and Semiconductors

To slow the PRC's progress in designing advanced chips, the US government has applied export controls and other restrictions, with varying levels of effect. The regulations we'll discuss here were intended to hinder China's progress and protect American companies.

The US uses the Export Administration Regulations (EAR) as the primary framework for export controls. An entire agency, the Bureau of Industry and Security, handles the task of preventing sensitive technology from being sold to the wrong people. This is not a bad thing; our foreign trade shouldn't directly hurt our national interests. However, as a side effect, the regulations also hurt America's position in the semiconductor supply chain.

Before 2020, any foreign-made products containing more than 25% American-origin components were subject to US export controls. This was the de minimis rule, and it meant a device (e.g. a graphics card) made in Taiwan could be export-controlled by the US if more than 25% of the parts were American. Since the US doesn't want China to have access to advanced chips or chipmaking equipment, that card couldn't be sold to China.

However, the PRC has a population of over a billion people. It's a huge market. Of course companies want to sell to China. So they started "designing out" American components and swapping them for Japanese or Korean components instead, to keep the American-origin percentage below 25%. Then they could sell their product to China without being subject to changing export controls.

In 2020 the Foreign Direct Product Rule (FDPR) was expanded to fight the de minimis loophole. Now, all products designed with US software (e.g. Cadence EDA tools) or made with US-designed machines (e.g. semiconductor fabs) are subject to US export controls regardless of what percentage of the finished product is American.

This greatly expanded the scope of US export controls, but it also drove companies doing business with China to transition away from US toolchains. Huawei is swapping US-sourced RF chips for Chinese ones, SMIC stockpiled Japanese lithography equipment, and European companies like STMicroelectronics gained market share by positioning themselves as "independent."

Asian and European companies benefitted from a "de-risking" spending spree by Chinese companies while American firms lost revenue and Japan kept supplying "mature node" (read: old) chip tech to China. The policies are slowing China's growth in semiconductors, but they've also harmed the competitiveness of American firms in the process.

I'm not going to pretend I know better on this one. National security and economic policy are full of cost-benefit analyses and weighing of trade-offs, and I have nowhere near a complete picture of what's involved. But surely a more coherent effort could have stopped advanced tech sales without incentivizing "de-risking" and "designing out" strategies in friendly countries.

How Engineers Could've Done Better

Engineers write product requirements that lay out expectations for a completed system. These requirements have some features that CAFE, the GDPR, and America's export control regime all lack.

Some relevant properties of engineering requirements are these:

• Atomicity: Compound statements are disallowed. An individual statement contains a single atomic requirement.
• Unambiguity: Subjective words like "timely" or "safe" aren't helpful; the courts just wind up defining them later.
• Verifiability: There is a feasible way to check that the requirement is met. This means you need a metric you can test.
• Quantifiability: The targets are specified precisely, so that they're unambiguous and verifiable.
• Consistency: Individual requirements don't contradict each other without defining a trade-off. You generally can't design a part with 300% margins that's also weight-efficient.
• Completeness: There are no unhandled edge cases, and negative requirements (what the system should not do) are included.

Would regulations be improved if their authors tried to apply these principles? I think so.

ADA & Requirements Done Right

The Americans with Disabilities Act (ADA) changed the way our buildings and cities look by requiring accessibility be baked into designs from the start. It's had a massive positive impact, even for people without disabilities. It's a case study in writing requirements-based regulation.

As an example, let's look at an excerpt of Title IV of the ADA.

The [Federal Communications] Commission shall, not later than 1 year after July 26, 1990, prescribe regulations to implement this section, including regulations that—
(A) establish functional requirements, guidelines, and operations procedures for telecommunications relay services;
(B) establish minimum standards that shall be met in carrying out subsection (c) of this section;
(C) require that telecommunications relay services operate every day for 24 hours per day;
(D) require that users of telecommunications relay services pay rates no greater than the rates paid for functionally equivalent voice communication services with respect to such factors as the duration of the call, the time of day, and the distance from point of origination to point of termination;
(E) prohibit relay operators from failing to fulfill the obligations of common carriers by refusing calls or limiting the length of calls that use telecommunications relay services;
(F) prohibit relay operators from disclosing the content of any relayed conversation and from keeping records of the content of any such conversation beyond the duration of the call; and
(G) prohibit relay operators from intentionally altering a relayed conversation.

This act empowers the FCC to set up hard engineering requirements for public infrastructure (in this case, telecomms equipment), so it's accessible to everyone. It's a brilliant and thoughtful piece of legislation.

One of the critiques of regulations that look like engineering requirements is that they're "brittle." If a piece of legislation had, say, mandated SHA-1 as a hashing algorithm, it would be outdated today since SHA-1 is insecure. An amendment would be needed. By assigning an agency the power to update the specifics, you get clear engineering requirements that can be updated as technology evolves.

You can read the regulations implemented by the FCC here. No matter where you scroll, it feels a lot like an actual Product Requirements Document, but look at the section I linked. Each line lays out an atomic, verifiable, quantifiable, and unambiguous requirement. An engineer can test a system or building and determine whether it's compliant.

I'm not saying the ADA hasn't created some perverse incentives. It has, especially around hiring decisions. But in general, the technical aspects are well-written and builders can reliably verify compliance.

Back to the GDPR

The GPDR, on the other hand, violates every principle in the list, but perhaps it violates "unambiguity" most egregiously. Many terms in it lack hard definitions, and in the absence of rulings from the Court of Justice of the European Union (CJEU), lower authorities have issued directives and advice.

In September the CJEU delivered a judgement on whether "pseudonymized" data is personal data. The court considers a dataset "personal data" if the intended recipient is able to reconstruct the identities of the people involved. However, the current guidance from the European Data Protection Board (EDPB) considers pseudonymized data "personal data" in all cases. The EDPB directive is more strict than the current legal precedent, so the Board and the Court disagree.

If the term "personal data" was better defined, this conflict would be impossible, but at the moment, it smells like a loophole. I could offload data processing to a third party without having to keep up with GDPR compliance, as long as they don't have access to the identifiers or look-up table. According to the court, this is okay. According to the EDPB, I am about to be fined for non-compliance. Who's right?

I'm convinced the ADA offers a better model. If all the GDPR did was empower the Data Protection Board to write specific requirements, the Board's definition would be the only definition. But that's not what happened, and now the EU is stuck with ambiguous definitions that must be hammered out in court.

Fixing CAFE

Like GDPR, I think CAFE violates the properties of good requirements, but in a different way. The issue with CAFE is more about "consistency" than "unambiguity." You cannot simultaneously ask for a more efficient vehicle fleet and allow vehicles to be less efficient if they are larger. Companies will just build larger vehicles to avoid the extra engineering effort.

CAFE failed because it "punished" the wrong thing. Congress wanted to reduce fuel usage. Yet the standards made it harder to build smaller vehicles. Shouldn't the regulations have made it harder to guzzle gas instead?

I'm no economist, but it seems like making more efficient cars cheaper, and making less efficient cars more expensive, would incentivize people to buy more efficient cars. Maybe a flat tax on less efficient vehicles would have the desired effect? I don't see it leading to ever-larger SUVs as a side effect, so perhaps it's already an improvement.

Good Requirements & Bad Requirements

This comes down to a fundamental idea in systems engineering: requirements should be descriptive, not prescriptive.

I'll say it again. Requirements should be descriptive, not prescriptive.

A descriptive requirement describes what things the system shall and shall not do. A prescriptive requirement prescribes how the system shall do "the things."

I picked up an example at Trane. One of the mechanical engineers was drafting some requirements for firmware in our lab automation set-up. He prescribed a seach algorithm: "iterate through the list until you've found the key." One of the software guys explained that we have better search algorithms, but his prescriptive requirement would block us from using them. He refactored it into a descriptive requirement, "the system shall locate the key in the list," and everyone was happy.

CAFE used prescriptive requirements to tell companies how efficient cars of certain sizes must be, instead of describing the objective: "less efficient cars shall be less attractive to operate." Congress tried to be too clever, and wound up creating a perverse incentive.

Consistent Export Controls

As with CAFE, Congress tried to be too clever with the 25% de minimis threshold. If the objective was to block sales of advanced technologies to China, that should have been the law. Instead, Congress said, "any advanced technology with more than 25% American components cannot be exported to China."

That created a perverse incentive to swap American-made components with foreign ones, so companies could keep the American-origin share below 25%. For years, that impacted America's position in the semiconductor industry. Congress and the federal government eventually course-corrected by implementing FDPR and getting other nations like the Netherlands and Japan to agree on sanctions. The question is, why didn't we go for that from the start?

The failure here was an example of Goodhart's Law:

When a measure becomes a target, it ceases to be a good measure.

The de minimis threshold (measure) became a target, and ceased to fulfil the purpose of blocking advanced tech from entering China. The FDPR is a better method of determining what is export-controlled and what is not, because there's no arbitrary measure companies can game. So the first Trump administration and the Biden administration actually fixed this regulatory issue by moving to more descriptive, more consistent requirements.

Again, national security and trade policy are especially complicated, and I'm not claiming to know better. I just think a more coordinated strategy could have minimized side effects.

Conclusion

The cobra policy, CAFE, the GDPR, and the EAR's de minimis rule all created perverse incentives. CAFE and the EAR were too prescriptive. The GDPR was too ambiguous. In contrast, ADA is generally unambiguous and consistent.

To be clear, I'm not saying a "product requirements mindset" would fix every piece of legislation, or that every piece of legislation should look like the ADA. I'm saying that the product requirements mindset can improve regulations which engineers have to follow.

If policy-makers drew on the principles of good requirements in their regulations, I think verifying compliance would be easier and creating perverse incentives would be harder. But so many of our politicians are lawyers by training (I think there's exactly one engineer in Congress), so they probably aren't used to writing "good requirements." They're used to writing legal documents. I'm not trying to bash the lawyers; they're just not usually exposed to systems engineering (as far as I know!)

This is one of those areas where people with a lot of different backgrounds have to collaborate. Governments aren't going to get tech policy right unless technical people can weigh in. I learned a lot from writing this post, and I'm interested in reading more in this space. If you have an example of regulation done well, or more stories about perverse incentives, I'd love to hear about it. If you think this idea is ridiculous, or that I've overlooked something important, let me know! You can email me at ethan.barry@howdytx.technology or reach out via LinkedIn.

Thanks for reading!

I love coffee; don't you? Even if (somehow) you don't, I think you'll enjoy this article.

I found out about this connection from reading Games & Decisions by Luce and Raiffa. It's a book on game theory, the mathematical study of decision-making. We use phrases like "zero-sum game" and "optimal strategy" in everyday speech, but they come from this field, which has an enormous number of applications.

The history of game theory is fascinating, thanks to the people who had a hand in forming it. Daniel Bernoulli and his brother Nicholas analyzed an early problem in the 18th century, and Émile Borel (of infinite monkey theorem fame) published some (ultimately incorrect) work on it in the early 1920s. By far the most important were John von Neumann, the great mathematician and computer scientist, and Oskar Morgenstern, who worked with him at the Institute for Advanced Study, when the Anschluss prevented him from returning to Austria.

They coauthored a book, Theory of Games and Economic Behavior, which built on von Neumann's 1928 proof of the minimax theorem for two-person zero-sum games (which Borel had discarded). It laid the foundation for work in the '50s by John Nash, who essentially rewrote the theory to give it more "teeth" when dealing with harder problems.

It's John Nash's contribution that I want to share, but to share it, I have to do a bit of setup.

A Simple Game

When I say "game" in this article, I mean this: \(n\) players, each with a set \(S_n\) of possible decisions, and a set of associated outcomes. If player 1 chooses his strategy 3, and player 2 chooses his strategy 5, then there is an outcome for that particular choice, telling what each player gains or loses. Time for an example!

Battle of the Sexes

Say a man and a woman want to see a movie. The man would rather see an action movie, and the woman would rather see a romantic comedy. However, they both prefer going to any movie together over going alone. We can rank the satisfaction or payoff with a numerical ordering.

The woman would most like to see the rom-com with the man, but she'd dislike seeing the rom-com alone. You can assign those the scores 2 and -1. She'd like watching the action movie with him, but would dislike seeing it alone, so assign those 1 and -1.

Now we can do the same for the man. You can swap their places in the ranking above—just reverse the 2 and the 1. We can make a table, called the payoff matrix.

In the table, the man's choice is on the left, and the woman's choice is on top. The ordered pairs are the payoffs \((w, v)\) to the woman and man, respectively. (Each choice, by the way, is an example of a pure strategy in the game.)

A game theorist wants to know what they should choose to maximize their happiness. A little critical thinking rules out either case where they aren't together. That leaves our hypothetical couple with the choices \((1, 2)\) and \((2, 1)\). Which is the "best?"

Of course, you already know the answer intuitively. There is no "best." Someone has to compromise (accept less than a 2), or both will be miserable!

Let's consider what happens when we repeat the experiment multiple times. In real life, people take turns compromising. They go to a rom-com one night, and an action movie the next. That's a great solution! As we repeat the experiment more times, the overall payoff gets closer to \((\frac{3}{2}, \frac{3}{2})\), the average of 1 and 2, and they're both equally happy.

There's also no way either one can get an advantage by changing their strategy, assuming the other player makes no change. If the man suddenly goes to action movies 3 out of 4 times, instead of 1 out of 2, he'll be worse off. This is what we call an equilibrium point, and in a sense, it's the "best solution" to the game.

Cool; so what?

The battle of the sexes was a very simple example, but it's not hard to see how you could model a complicated situation like a barter, labor dispute, or international political squabble in this framework—each side has some choices that lead to some outcomes. We also saw that, even though there was no good solution for either player in terms of a single pure strategy, there was one in terms of mixed strategies.

Question: Is it always the case there is a best solution, or equilibrium point?

The answer is yes—this is John Nash's major contribution. It's called the Nash equilibrium, and the proof that it exists is actually very beautiful.

Coffee & the Existence of Equilibrium Points

Let's put game theory on hold for a moment, and back up to the end of the 19th century. Henri Poincaré was trying to answer an old question—is the Solar System stable? He considered the motion of objects flowing on a level surface, and made a discovery. However, this discovery is usually credited to (and was actually proved by) L. E. J. Brouwer. There's a nice (apocryphal) story associated with it, so I'll tell that instead. Just remember that Poincaré published on it first.

The story goes like this: As Brouwer was stirring sugar into his coffee one morning, he noticed that one point in the coffee always seemed to be stationary. He realized this must always be the case—what would it look like if every point on the surface were moving? It's not physically possible.

The Theorem We Need

This visual intuition is connected to an actual mathematical fact, now known as Brouwer's fixed-point theorem. I'll give you an even more intuitive example, after the mathematical statement.

Brouwer's fixed-point theorem: Every continuous function from a closed disk onto itself has at least one fixed point.

Let's break this down. We have a disk (the surface of the coffee) which is closed (has a boundary), and a function mapping points on the disk, to other points on the disk. (A stirring motion moves the coffee around, but it never leaves the cup.)

Also, this function is continuous, which means all points initially in a neighborhood stay in a neighborhood after applying the function. (The coffee doesn't "jump" to a new location; it moves "continuously.")

This actually applies to any nonempty, convex, compact subset \(U\) of Euclidean space, and continuous transformation \(T : U \rarr U\).

Now for the visual example. Imagine you're holding a map of the state or country you're standing in. You know the transformation from the map to the country is always continuous—you just scale up. Points in a neighborhood remain neighbors. Spread the map on a table, and there will always be a "You Are Here" point that's the same on the map as it is in the country.

Are you convinced yet? It seemed pretty intuitive to me—it's almost one of those things where you say "I wish I'd thought of that!"

Nash's Proof

Let's return to game theory. So far, we have a game, and Brouwer's fixed point theorem. If we can make a transformation from the set of payoffs, to the set of payoffs, in such a way that a fixed point must be an equilibrium point, we have proved the existence of an equilibrium point for all finite games.

Now I can sketch the outline of Nash's proof. I'll follow the notation and structure of Nash's doctoral thesis, found here. You can follow the entire thing if you've taken a bit of calculus, yet it was revolutionary in its field! If you just want the TL;DR you can skip to the summary.

Here are the definitions we need:

Definition: Finite game: A finite game is \(n\) players, ea. with a finite set of pure strategies, and for ea. player \(i\), a payoff function \(p_i\) which maps all \(n\)-tuples of pure strategies into an ordered set, e.g. \(\Bbb{R}\) or \(\Bbb{Z}\).

Definition: Mixed strategy: A mixed strategy of a player \(i\), denoted \(S_i\), is a collection of weights or probabilites (like \(\frac{1}{2}\) in "Battle of the Sexes") on his pure strategies. \[S_i = \sum_\alpha c_{i\alpha} \pi_{i\alpha}\] where \(\sum_\alpha c_{i\alpha} = 1\) & \(c_{i\alpha} \ge 0\) (it's a probability), and \( \pi_{i\alpha}\) is player \(i\)'s \(\alpha\)-th pure strategy.

Definition: Payoff function: We can extend the payoff function in our definition of the finite game to cover all \(n\)-tuples of mixed strategies. We'll use this definition of \(p_i\) from now on. So for any mixed strategy \(𝒮\), we'd write:
\[p_i(𝒮) = p_i(S_1, S_2, \dots , S_n).\] Also note that any \(𝒮\) is a point in \(n\)-dimensional vector space, and the set of all \(𝒮\)s is a convex subset of that space. (Think about our payoff matrix, but for mixed strategies—that's what this represents.)

Definition: Equilibrium point: a point where no player can get an advantage by changing their strategy, assuming the other players make no change. Formally, a mixed strategy \(𝒮\) is an equilibrium point if and only if for every \(i\), \[p_i(𝒮) = \underset{\text{all}\>r_i\text{'s}}\max\,\big[\,p_i(𝒮; r_i)\,\big],\] given that \((𝒮; r_i)\) means "substitute mixed strategy \(r_i\) for the \(i\)-th player's strategy."

Nash actually constructs a sequence of continuous transformations, which we'll look at below. First, let's specify exactly what we're proving.

Theorem: Every finite game has an equilibrium point.

To prove this, we need some continuous functions of a mixed strategy \(𝒮\). Remember that \(p_{i\alpha}(𝒮)\) is player \(i\)'s payoff from pure strategy \(\pi_{i\alpha}\) when the other players use their strategies in \(𝒮\). For each \(\lambda \in \Bbb{Z}\), \[\begin{gather}q_i(𝒮) = \underset{\alpha}\max\>p_{i\alpha}(𝒮) \\ \Phi_{i\alpha}(𝒮, \lambda) = p_{i\alpha} - q_i(𝒮) + \frac{1}{\lambda} \\ \Phi^+_{i\alpha}(𝒮, \lambda) = \max \big(0, \Phi_{i\alpha}(𝒮, \lambda)\big). \end{gather}\]

Function (2) is either \(\frac{1}{\lambda}\), or \(\frac{1}{\lambda}\) plus some negative number. That lets us show the rational function \[ C'_{i\alpha} = \frac{\Phi^+_{i\alpha}(𝒮, \lambda)}{ \sum_{\beta}\Phi^+_{i\beta}(𝒮, \lambda) } \] is also continuous, since the denominator is never zero.

Nash then defines \( S'_i(𝒮, \lambda) = \sum_\alpha \pi_{i\alpha} C'_{i\alpha}(𝒮, \lambda) \), which are different mixed strategies for every player \(i\). We combine all those into a new point \(𝒮'\), given by \( 𝒮'(𝒮, \lambda) = (S'_1, S'_2, \dots , S'_n) \). The transformation \(𝒮\rarr 𝒮'(𝒮, \lambda)\) is continuous because all the operations that compose it are continuous. The space of all mixed strategies is also non-empty, convex, and compact (lacking "holes"). Thus, there will be a fixed point \( 𝒮^* \), by the Brouwer fixed point theorem.

We haven't yet proved that this fixed point must be an equilibrium point. Nash does that by contradiction—assuming it isn't an equilibrium point—which implies that one of its strategies \(\pi_{i\alpha}\) is non-optimal. If that were so then we can show, based on functions (1)–(3), that \(𝒮^*\) is not a fixed point. (I'll leave that to Nash's thesis. See page 6.) This contradicts the fact that \(𝒮^*\) is indeed a fixed point. Therefore, each strategy in \(𝒮^*\) must be optimal against the others. No player can gain by changing to a different strategy, so by definition, \(𝒮^*\) is an equilibrium point.

Proof Summary

Nash created some transformation \(T\), mapping \(n\)-tuples of mixed strategies into \(n\)-tuples of mixed strategies, where \(T\) has two properties:

• \(n\)-tuple \(𝒮\) is an optimal strategy if and only if \(T(𝒮) = 𝒮\), that is, \(𝒮\) is a fixed point under \(T\).
• \(T\) satisfies the conditions of Brouwer's fixed-point theorem, and consequently has a fixed point.

If those are the case, and Nash proved they must be, then any \(n\)-player finite game does indeed have at least one Nash equilibrium. This means there is a solution to any "game"—meaning an abstract model for anything from a couple at the movies to national governments—which no player can improve by switching strategies. To prove that, we used a single theorem in topology that a mathematician allegedly discovered when stirring sugar into his coffee.

If that's not a fascinating connection, I don't know what is.

References

John Nash: Non-Cooperative Games.
R. Duncan Luce & Howard Raiffa: Games & Decisions.
Wikipedia: "Brouwer fixed-point theorem"
Wikipedia: "John Forbes Nash, Jr."

As always, thanks for reading! If you have a comment, I'd love to hear from you! My email is ethan.barry@howdytx.technology. You can also reach me at my LinkedIn or GitHub, linked in the header.

At the end of the year, I always put some thoughts down in one of my notebooks. Writing things out helps me focus on them, and a bulleted list is easier to understand than an abstract thought cloud. I figured it would be nice to do it on my blog this year, especially since no one reads this anyway.

The spring was great—it feels like it was ages ago I finished Calculus II, which was the highlight of that semester, and I'm sure I've forgotten many things already. What I haven't forgotten is the group of friends I made there. I ran into most of them on campus at one point or another in the fall, and it's a blast watching them dive into their majors.

That informed a decision I made at the beginning of the fall semester. I'd never made social activities a priority before, and the spring was over so quickly I didn't get to know more than a few people. But for the fall, I wanted social activity to be on the same footing as my studies. Universities are about more than grades, after all, and we can learn as much from our peers as from our professors.

I'm happy to report that it paid off—I got to do a team project with an amazing group of fellow CS majors, and I got involved in student organizations, mainly as secretary of the Assoc. for Computing Machinery chapter. I'm grateful for the opportunities, and I learned a lot about teamwork from both.

The fall was more fun in some ways, but it was worse in others. Despite being 16 credit-hours instead of 19 like the spring, it was a harder workload. But, because of my choice to make social activities a priority, I still spent plenty of time with friends. We went out to lunches, discussed everything from music to life goals, and generally had a grand old time.

I'll finish with a couple specific things that happened to me, and then the bulleted list of observations I mentioned. First, I found out I have three semesters of compsci courses left, not two. (One would have a single class.) So, I decided to go all-in and double major in mathematics. That also gave me a chance to participate in the Honors College; I'll get to do a two-semester research project in my major.

Second, a couple days after that, I received an offer from Trane Technologies, and will be joining them as a co-op/intern on the controls software team in January. I'm incredibly excited to get a taste of the industry and solve real-world problems. It means I'm not taking classes in the spring, but I'll be back in the fall with a new perspective and, I'm sure, more friends!

The List of Lessons Learned

In no particular order, the lessons I learned (or should have learned!) from this year:

• Classes are only exciting if teachers are excited to teach them, and a good teacher makes hard material enjoyable.
• It always pays to listen to people. Time is precious, and they notice when you invest it.
• It always pays to ask questions. You're guaranteed to be no worse off than before, and you'll have rephrased the problem to ask the question. Oftentimes, that reveals the solution.
• When they tell you a physics lab is one credit-hour, it's a trap. They're lying.
• The Probability & Statistics class (MATH 3351) is evil. It looks like math, but it isn't, because math is fun.
• Social life matters more than grades. I walked out of this semester with a 3.9 GPA, not a 4.0, and it was worth it.
• It's always better to handle problems early. (I keep forgetting and relearning this one.)
• Helping people learn, whether as a TA for Data Structures & Algorithms, or as a friend giving calculus tips, is extremely rewarding. The "Aha!" moment, when someone finally gets it, is worth every minute spent.

I'd love to hear from anyone reading this, if there are such people. Stranger things have happened on the internet. My email is ethan.barry@howdytx.technology.

Because everyone needs more LaTeX in their life.

At least, the easy-to-use kind.

Most folks in STEM fields know about \(\LaTeX\). It's the gold standard for scientific writing, desktop publishing, and math typesetting. Its algorithms still beat MS Word's, Google Docs', or even Adobe's InDesign's. (If you haven't heard of \(\LaTeX\), it's a way to make PDFs with code.)

The website Overleaf.com is an easy-to-use experience for anyone who would rather sit down and write than fiddle with \(\LaTeX\) installation on, say, a Windows computer. Because of that, many thousands of students, academics, and professionals worldwide use it daily to typeset documents. The platform has a huge collection of templates, which can provide your document's formatting and looks.

The Problem

Half the problem is that most of the templates are boring, and only some are exceptional. Of course, any one you choose will produce a better document than MS Word would give you. The issue is that most users never tap into the superpowers of \(\LaTeX\). This is the fault of something called the "document class." A document class defines macros (shortcuts) for the user to set text with. The document classes everyone uses (article, book, &c.) are old, try to be as backwards-compatible as possible, and break easily if you try to be too cute.

Ask me how I know.

Anyway, for a new document, people should ideally be using the memoir class. It is a modern class (from 2001!) which is entirely dedicated to ergonomics. The original developer wanted users to create professional documents without knowing anything about typography. It provides a whole host of deeply integrated functionality that article and book can only dream of.

Side Note: The class has 400 pages of documentation, 500 pages in the user manual, and even comes with a 130 pages of typographical notes by the developer. You will not lack for help resources if you adopt memoir.

Having a general overabundance of okay-ish templates wouldn't be a problem on its own. The other half of the problem is that many universities provide \(\LaTeX\) templates, while mine does not! So, can we create a decent template using memoir that fixes UT Tyler's lack of an Overleaf template?

The Solution

We have a default article document that looks like this:

But, we want to turn it into something pretty. If we simply change from article to memoir as our document class, very little happens. We need to do more. I'll walk through the changes I made, and a beginner-level \(\LaTeX\) user should have no problem following along in the code, available here.

Here is a page from the final result:

So, what changed?

Visible Changes

First, the font. There's nothing wrong with the default latin-modern font, but it is a bit spindly at larger sizes, and we want something classic. You'll find there is no better option than EB Garamond. If you haven't overridden my CSS, this blog is set in EB Garamond as well. You can insert \usepackage{ebgaramond} to set it as the main font. I've included a package that sets it as the math font, too. Notice how nice the Schrödinger equation looks!

Second, the margins are different—I used a wider right margin so I can insert margin notes—which gives the document a more book-like feel. I want to use margin notes to introduce variables. I saw this done in an older journal, and I think it's a very handy idea. It's a shame we don't see it more often. Notice how \(x\) and \(z\) are shown below. If you were skimming the paper for some variable's definition, you could easily find it!

Boldface is used sparingly, with italics and smallcaps (\scshape) providing most of the emphasis. Bold faces are not a major element of classic typography, so it's nice to minimize their use. It goes a long way to making the document more readable.

Finally, my favorite part is the title page and frontmatter. I'll let the images speak for themselves.

The title page is inspired by this template by Mario Vilar on Overleaf, though I didn't use any code from it in this document. I especially like his faint background. It's similar to my research paper's title page here.

Changes "Under the Hood"

The code itself is fanatically organized and heavily commented. I want it to be as easy as possible to do whatever needs to be done. It has sane defaults for document divisions and page styles, with most of the options visible in the preamble.

You can explore it for yourself on GitHub or Overleaf.

Conclusion

I am excited to see whether this gets any notice on Overleaf from my fellow students, but regardless of whether it does, I want to use it for my own papers. Hopefully none of the professors will be too draconian about formatting rules in the coming semesters.

As for extensions to the template, I'm adding some macros for margin notes and variables in the margins. I realized after I took the screenshots that the margin notes should be set \raggedright instead of justified. I think a macro that takes the note as an argument and automatically formats it will be useful.

As always, please send me any feedback you have at ethan.barry@howdytx.technology, or drop me a message on LinkedIn!

Originally all I wanted to do with numerical analysis was string together some functions and plot an ephemeris, given the right data. Planetary motions are fascinating, and I think it will be a really cool project. However, after reading the first chapter of the textbook Celestial Mechanics, where the author lists important numerical methods, I realized that numerical analysis by itself is something interesting.

In the last section of the first chapter he gives a list of programs or routines that are useful to him. One of those was using Gaussian quadrature to compute a definite integral. I already discussed the problem of quadrature in the post on Simpson's Rule, and I mentioned Gaussian quadrature specifically in the post on Lagrange interpolation. But, it's still worth revisiting the question of what quadrature means, before I share what I learned about Gaussian quadrature.

Computing Definite Integrals

In calculus we learn that a definite integral is actually the area beneath a curve. We approximate this with a sum of the areas of many small slices. In fact, the exact value is the number the sum approaches as the number of slices goes to infinity. This is called the Riemann sum approach, and essentially, it chops the integral up in nice rectangular bits. If we wanted to approximate an integral, one way would be to chop it up into lots of rectangles, and add their areas. If we want a closer approximation, then we just make finer slices.

This is certainly a fine method, but it would be a nightmare to compute by hand. Every slice needs its own function evaluation. Even in a computer, this can be slow. And, since we would be using rectangles to approximate a smooth function, this would still be inaccurate.

If, instead of using a flat-topped rectangle, we used a slope-topped trapezoid as our slice shape, we could make the slope line up with the function. That would mean our calculation would be more accurate, and would require fewer function evaluations. This is the trapezoid rule.

And, if we use a quadratic curve as the top of the slice, it can fit the function even more tightly still. That was Simpson's Rule, from my earlier post. In fact, there's a whole class of quadrature rules called Newton-Cotes formulas, which extend this idea, using higher-order polynomials to interpolate the points.

With these methods, we have one degree of freedom. We get to choose the weighting coefficients in equations like Simpson's Rule:

\[\int_a^b f(x)\,dx \approx \frac{1}{3}(f(a) + 4f(c) + f(b))\]

The big idea behind Gaussian quadrature is that it allows us to choose both the weights and the locations at which we evaluate \(f(x)\), giving us two degrees of freedom. With Newton-Cotes formulas, our accuracy increased at \(\mathcal{O}(n)\), as we added more slices. Now, we can converge to the result at \(\mathcal{O}(n^2)\), with the extra degree of freedom.

Let's examine the main idea.

Gaussian Quadrature

Simpson's Rule and the rest of the Newton-Cotes formulas essentially represent the integrand with polynomials. Gaussian quadrature lets us represent our integrand as a polynomial times a weight function. This is good, because some functions aren't nicely interpolated by polynomials. We can also evaluate at different points by applying weights to the evaluations.

This gives us the general Gaussian quadrature rule

\[\int_{-1}^1 W(x)f(x)\,dx \approx \sum_{1=0}^n w_if(x_i)\]

The polynomials we need to use are called Legendre polynomials. So, in our equation, \(n\) is the number of sample points, \(w_i\) are the weights, and \(x_i\) are the roots of the \(n\)th Legendre polynomial. The integral is given on the interval \([-1,1]\) because that is where Legendre polynomials have most of their useful properties, but we can scale up or down to any interval \([a,b]\).

We can compute the roots on the fly for any order \(n\), within reason. There is a technique called the Newton-Raphson method which finds roots of nonlinear functions. It takes a starting guess and iteratively improves it until the difference between the new value and the old one is within tolerance. It's a surprisingly simple technique.

\[x_{new} = x_{old} - \frac{f(x_{old})}{f'(x_{old})}\]

Now we just need the polynomials! There's a recurrence relation for them, so starting with the first two, we can build all the Legendre polynomials after that.

\[\begin{align}p_{-1}(x) &= 0\\ p_0(x) &= 1\\ p_{j+1}(x) &= (x-a_j)p_j(x) - b_jp_{j-1}(x),\quad j = 1,2,...,n\end{align}\]

We need a good guess to use Newton-Raphson. Usually the roots are closely approximated by \[r=\cos\left(\pi\frac{i + \frac{3}{4}}{n + \frac{1}{2}}\right)\quad i = 0,...,\frac{(n+1)}{2}\] where \(n\) is the order of the quadrature.

Once we have the weights \(w_i\), given by \[w_i = \frac{2}{(1-x_i^2)\left[p_n'(x_i)\right]^2}\] where \(p_n'\) is the derivative of the \(n\)th Legendre polynomial. Finally, we can mix all these numbers together in some code. But first, I want to show a result this algorithm gives.

Here's a test problem with the correct approximation to 9 decimal places: \[\int_0^1 \frac{ x^4 }{ \sqrt{ 2(1+x^2) }}\,dx \approx 0.108\,709\,465\,052\,586\,44\] The Rust code I wrote returns this value to nine decimal places with a 7-point Gauss-Legendre quadrature. That means the function only had to be evaluated seven times. I also used Simpson's Rule with seven evaluations, and it gave me an answer of 0.108 710, which is not bad.

Now, if I raise the Gaussian quadrature's order to eight, I get the answer 0.108 709 465 037, which is correct to 10 decimal places. If I raise the Simpson's Rule order to eight as well, I get 0.108 709 862, which is still only correct to six decimal places.

The code is below.

pub fn gaussian_quad<F>(f: F, a: f64, b: f64, n: u32) -> f64
where
    F: Fn(f64) -> f64,
{
    const NEWTON_PRECISION: f64 = 1.0e-10;
    let mut res = 0.0;
    // Scale the bounds of integration a and b to [-1, 1]
    // using the change of variables
    // x = st + c for t ∈ [-1, 1]
    // where
    // s = (b - a) / 2
    // and
    // c = (b + a) / 2.
    let s = (b - a) / 2.0;
    let c = (b + a) / 2.0;
    let m = (n + 1) / 2;
    let mut roots = vec![0.0; n as usize];
    let mut weights = vec![0.0; n as usize];
    for i in 0..m {
        // Approximate the ith root of this Legendre polynomial.
        let mut z = f64::cos(consts::PI * (i as f64 + 0.75) / (n as f64 + 0.5));
        let mut z1 = z - 1.0; // Ensure an initial difference.
        let mut p_prime = 0.0;
        while (z - z1).abs() > NEWTON_PRECISION {
            let mut p1 = 1.0;
            let mut p2 = 0.0;
            // Build the Legendre polynomial using the recurrence relation here...
            for j in 0..n {
                let p3 = p2;
                p2 = p1;
                p1 = ((2.0 * j as f64 + 1.0) * z * p2 - j as f64 * p3) / (j + 1) as f64;
            }
            // p1 is now the desired Legendre polynomial.
            // let p_prime be its derivative, computed by a known relation.
            p_prime = n as f64 * (z * p1 - p2) / (z * z - 1.0);
            z1 = z;
            z = z1 - p1 / p_prime; // Newton-Raphson method.
        }
        // z is now the ith root of the polynomial.
        roots[i as usize] = z;
        roots[(n - 1 - i) as usize] = -z;
        weights[i as usize] = 2.0 / ((1.0 - z * z) * p_prime * p_prime);
        weights[(n - 1 - i) as usize] = weights[i as usize];
    }

    for i in 0..n as usize {
        res += weights[i] * f(s * roots[i] + c);
    }
    res * s // Scale back to normal size.
}

Now, if we have a function that is well-behaved, this will work nicely. But if there's an irregularity in just one or two places, we have to raise the number of points just to cover those two subintervals. Wouldn't it be nice if we could somehow adapt to a problem like that, and only do more evaluations when we need them?

Adaptive Quadrature

Remember that the heart of numerical integration is chopping area into smaller pieces. If there's an interval where the function goes crazy, and one where it's smooth, then let's not worry about the smooth interval too much, and focus on the crazy one. How do we measure "craziness?" The general method is to evaluate the integral on the interval twice—once at low precision, then again at high precision—which gives us a difference. If that difference is within a tolerance, then we leave it alone. If not, we split the interval and repeat the process.

You might think that this will cost us a lot of function evaluations. It will take more, but maybe we can minimize the number needed. And an adaptive routine is already saving us evaluations in comparison to just globally increasing evaluations to boost accuracy, like with Simpson's Rule.

With Gaussian quadrature, a low-precision rule will not have the same \(x\)-coordinates as a high-precision one. That means we can't reuse any of those function values. What we're looking for is a nested quadrature, which "nests" a higher-order quadrature inside a lower-order one. There are two options, Clenshaw-Curtis, and Gauss-Kronrod. Both were invented around 1960; the first by two Western (British, I think) mathematicians, and the second by a Soviet mathematician, Aleksandr Semyonovich Kronrod.

This particular algorithm, from the paper "Adaptive Quadrature—Revisited" by W. Gander and W. Gautschi, is an example of Gauss-Lobatto-Kronrod quadrature. Kronrod's extension of Gaussian quadrature allows us to reuse function evaluations for higher precision, and the "Lobatto" part means we use the endpoints as well.

A Kronrod Rule uses \(2n+1\) points, where \(n\) is the number of points for the Gaussian quadrature. This makes the Kronrod Rule exact for a polynomial up to degree \(3n+1\). The extra \(x\)-values (abscissas) are nowadays calculated with some linear algebra tricks, as described in the paper. They can also be found by solving non-linear equations with Newton's Method, but because either one would be expensive, we'll use tabulated values.

As an example of what Gaussian quadrature and the Kronrod extension are really doing, I produced a little graphic with gnuplot. Here is a 4-point Gauss-Lobatto quadrature of our example \(\int_0^1 \frac{ x^4 }{ \sqrt{ 2(1+x^2) }}\,dx\). The lines pass through the \(y\)-values of the Gaussian abscissas, so the abscissas lie below the intersection of the line and the curve. The Kronrod abscissas would occur somewhere between the Gaussian ones.

The function I created has a pretty simple interface. I didn't include controls on how many iterations to allow, or what precision to calculate the integral to. It simply computes until the answer is stable to within machine precision.

pub fn gaussian_kronrod_quad<F>(func: F, a: f64, b: f64) -> f64
where
    F: Fn(f64) -> f64,
{
/* SNIP; CODE BELOW */
}

This has given very good results on the integrals I've tested. For example, on our problem from earlier, \[\int_0^1 \frac{ x^4 }{ \sqrt{ 2(1+x^2) }}\,dx \approx 0.108\,709\,465\,052\,586\,44\] I get the answer to the ~17 digits allowed by a double precision float f64.

The basic steps are:

1. Evaluate a 4-point Gauss-Lobatto quadrature on the scaled interval.
2. Evaluate a 7-point Kronrod extension using the 4-point values.
3. Evaluate a general estimate for the entire integral with a 13-point Kronrod extension using stored values from the previous two.
4. If the 4- and 7-point values differ by more than machine precision, split the interval into six parts, and perform the same 4- & 7-point procedure on those. Add their results.

This is a very good quadrature method, comparable to Julia's QuadGK.jl package, or some tools from other dedicated systems, like GNU OCTAVE and MATLAB. In fact, in the paper linked above, the authors show that it was actually an improvement over the quadatures available at the time in MATLAB, and the NAG and IMSL libraries.

Conclusion

The speed and type safety of Rust make it a very attractive option for scientific computing, but so far, it hasn't made much headway against more established languages like FORTRAN, Python, and Julia. The numerical analysis and linear algebra packages in Rust are improving, but they're often less-complete than counterparts elsewhere.

I think the killer application of Rust in this space might be concurrency. Rust makes parallelization much more achievable. Rust's SIMD wrappers are nicer than most languages' equivalents as well, so parallelizing and vectorizing sci-comp workloads is a unique benefit Rust can offer. I'm excited to see what happens in this space.

/// An adaptive 4-point Gauss-Lobatto quadrature routine.
/// Kronrod's Rule is used to reduce the computations needed,
/// while using 7- & 13-point Kronrod quadratures for error checking.
///
/// ## Inputs
/// - f: The integrand (must be integrable on \[a,b\])
/// - a: The lower bound of integration.
/// - b: The upper bound of integration.
///
/// NOTE: The algorithm assumes a < b.
/// ## Outputs
/// An `f64` value containing the definite integral calculated to machine precision.
/// ## Examples
/// ```rust
/// use symmetria::quadrature::gaussian_kronrod_quad;
/// use std::f64::consts;
/// let f = |x: f64| x.cos();
/// let val = gaussian_kronrod_quad(f, 0.0, (consts::PI / 2.0));
/// assert!((1.0 - val).abs() < 1.0e-17);
/// ```
pub fn gaussian_kronrod_quad<F>(func: F, a: f64, b: f64) -> f64
where
    F: Fn(f64) -> f64,
{
    // This algorithm was adopted from the paper
    // "Adaptive Quadrature --- Revisited" by Walter Gander and Walter Gautschi.
    // published in BIT Numerical Mathematics, vol. 40, No. 1, pp. 84-101.

    // Scale the bounds of integration a and b to [-1, 1].
    // This converts
    // \int_a^b f(x)dx to \int_{-1}^1 f(sx + c)sdx
    // using the change of variables
    // x = scale * t + coscale for t ∈ [-1, 1]
    // where
    // scale = (b - a) / 2
    // and
    // coscale = (b + a) / 2.
    let scale = (b - a) / 2.0;
    let coscale = (b + a) / 2.0;
    // Now, f must always be evaluated at scale * x + coscale, not simply x, and the answer multiplied by scale.

    let gl4_weights = [1.0 / 6.0, 5.0 / 6.0];
    let gl4_fnevals = [
        func(scale * -1.0 + coscale),
        func(scale * 1.0 + coscale),
        func(scale * -1.0 / 5.0_f64.sqrt() + coscale),
        func(scale * 1.0 / 5.0_f64.sqrt() + coscale),
    ];
    // Next, use the 4-point Gauss-Lobatto quadrature given by
    // \int_{-1}^1 f(x)dx ≈ \frac{1}{6}[f(-1) + f(1)]
    //                      + \frac{5}{6}[f(-\frac{1}{\sqrt{5}}) + f(\frac{1}{\sqrt{5}})]
    let gl4 = gl4_weights[0] * ((gl4_fnevals[0] + gl4_fnevals[1]) * scale)
        + gl4_weights[1] * ((gl4_fnevals[2] + gl4_fnevals[3]) * scale);

    let kron7_fnevals = [
        func(scale * -((2. / 3.0_f64).sqrt()) + coscale),
        func(scale * (2. / 3.0_f64).sqrt() + coscale),
        func(coscale), // scale * 0.0 is 0.
    ];
    // Then we compare this to its 7-point Kronrod extension, given by
    // \int_{-1}^1 f(x)dx ≈ \frac{11}{210}[f(-1) + f(1)]
    //                      + \frac{72}{245}[f(-\sqrt{\frac{2}{3}}) + f(\sqrt{\frac{2}{3}})]
    //                      + \frac{125}{294}[f(-\frac{1}{\sqrt{5}}) + f(\frac{1}{\sqrt{5}})]
    //                      + \frac{16}{35}[f(0)]
    let kron7 = (11. / 210.) * ((gl4_fnevals[0] + gl4_fnevals[1]) * scale)
        + (72. / 245.) * ((kron7_fnevals[0] + kron7_fnevals[1]) * scale)
        + (125. / 294.) * ((gl4_fnevals[2] + gl4_fnevals[3]) * scale)
        + (16. / 35.) * (kron7_fnevals[2] * scale);

    // Finally, we compare all of them to the 13-point Kronrod extension, given by
    // \int_{-1}^1 f(x)dx ≈ A * [f(-1) + f(1)]
    //                      + B * [f(-x_1) + f(x_1)]
    //                      + C * [f(-\sqrt{\frac{2}{3}}) + f(\sqrt{\frac{2}{3}})]
    //                      + D * [f(-x_2) + f(x_2)]
    //                      + E * [f(-\frac{1}{\sqrt{5}}) + f(\frac{1}{\sqrt{5}})]
    //                      + F * [f(-x_3) + f(x_3)]
    //                      + G * [f(0)]
    // where
    // A = .015827191973480183087169986733305510591,
    // B = .094273840218850045531282505077108171960,
    // C = .15507198733658539625363597980210298680,
    // D = .18882157396018245442000533937297167125,
    // E = .19977340522685852679206802206648840246,
    // F = .22492646533333952701601768799639508076,
    // G = .24261107190140773379964095790325635233,
    // and
    // x_1 = .94288241569547971905635175843185720232,
    // x_2 = .64185334234578130578123554132903188354,
    // x_3 = .23638319966214988028222377349205292599.
    let a_weight = 0.015_827_191_973_480_183;
    let b_weight = 0.094_273_840_218_850_05;
    let c = 0.155_071_987_336_585_4;
    let d = 0.188_821_573_960_182_45;
    let e = 0.199_773_405_226_858_52;
    let f = 0.224_926_465_333_339_51;
    let g = 0.242_611_071_901_407_74;

    let x_1 = 0.942_882_415_695_479_7;
    let x_2 = 0.641_853_342_345_781_3;
    let x_3 = 0.236_383_199_662_149_88;

    let kron13 = a_weight * (gl4_fnevals[0] + gl4_fnevals[1]) * scale
        + b_weight * (func(scale * -x_1 + coscale) + func(scale * x_1 + coscale)) * scale
        + c * (kron7_fnevals[0] + kron7_fnevals[1]) * scale
        + d * (func(scale * -x_2 + coscale) + func(scale * x_2 + coscale)) * scale
        + e * (gl4_fnevals[2] + gl4_fnevals[3]) * scale
        + f * (func(scale * -x_3 + coscale) + func(scale * x_3 + coscale)) * scale
        + g * kron7_fnevals[2] * scale;

    let x = vec![a, b];
    let y = batch_eval(&func, &x);

    gkquad_recursive(&func, a, b, y[0], y[1], kron13)
}

/// A helper function for gaussian_kronrod_quad().
fn gkquad_recursive<F>(f: &F, a: f64, b: f64, fa: f64, fb: f64, iest: f64) -> f64
where
    F: Fn(f64) -> f64,
{
    let h = (b - a) / 2.0;
    let m = (a + b) / 2.0;
    let alpha = f64::sqrt(2.0 / 3.0);
    let beta = 1.0 / f64::sqrt(5.0);

    // Our new evaluation points...
    let mll = m - alpha * h;
    let ml = m - beta * h;
    let mr = m + beta * h;
    let mrr = m + alpha * h;
    let x = vec![mll, ml, m, mr, mrr];
    let y = batch_eval(&f, &x);

    let fmll = y[0];
    let fml = y[1];
    let fm = y[2];
    let fmr = y[3];
    let fmrr = y[4];

    let i2 = (h / 6.0) * (fa + fb + 5.0 * (fml + fmr));
    let i1 =
        (h / 1470.) * (77. * (fa + fb) + 432. * (fmll + fmrr) + 625. * (fml + fmr) + 672. * fm);

    if iest + (i1 - i2) == iest || mll <= a || b <= mrr {
        if m <= a || b <= m {
            // At this point, we have exhausted the machine's precision.
            // We should handle this somehow...
        }
        i1
    } else {
        // Sum of the integrals on the ranges [a, mll], [mll, ml], [ml, m], [m, mr], [mr, mrr], [mrr, b];
        gkquad_recursive(f, a, mll, fa, fmll, iest)
            + gkquad_recursive(f, mll, ml, fmll, fml, iest)
            + gkquad_recursive(f, ml, m, fml, fm, iest)
            + gkquad_recursive(f, m, mr, fm, fmr, iest)
            + gkquad_recursive(f, mr, mrr, fmr, fmrr, iest)
            + gkquad_recursive(f, mrr, b, fmrr, fb, iest)
    }
}

NOTE: The function doesn't initially check whether it ought to evaluate recursively or not, but a simple if-else statement should catch that.

At the end of the day, our computers really only understand ones and zeros, the state of logical gates in a piece of finely-engraved silicon, and the layout of charges in chips of memory. Programs, data, noise – they all look the same. Yet, we're able to write instructions in a text-based programming language, that are somehow turned into the bits the computer speaks natively. How does that happen? With a compiled language, like Pascal, C, or Rust, the compiler does the job of turning the language into processor instructions, which are then assembled into bits and stored in special executable formats. Today, I want to explore the deep and beautiful ways in which a programming language can be translated into machine instructions by building a compiler.

What are we compiling?

This is a good question to start with. There are many languages out there – the major, general-purpose programming languages everyone thinks of, and then less common ones, like gnuplot instructions, or sudo configuration, or other Domain-Specific Languages (DSLs). We could write a program to interpret or compile any Turing-complete language ever invented, but most of them are complicated. Instead, we'll choose a simple, classic language that was invented specifically for learning compiler techniques with, PL/0.

PL/0 is an instructional language created by Niklaus Wirth, the same guy who invented Pascal, and contributed so much in the earlier days of computer science. He provided a basic compiler for it in his book Algorithms + Data Structures = Programs, which is a great read if you can find it. He wrote his in Pascal, which is difficult to understand if (like me) you aren't used to Pascal and its lack of abstractions around string processing. I like Dr. Brian Callahan's guide much better – it's written in C, and makes much more sense. I'll be following his guide for most of this series, and I highly recommend checking it out as well, especially when I glaze over details.

PL/0 is a simple language. Unless we extend it, it has one data type, VAR, which is an integer. Here is a simple program, stolen from Wikipedia.

VAR x, squ;

{Square x and load to squ}
PROCEDURE square;
BEGIN
   squ := x * x
END;

BEGIN
   x := 1;
   {Square numbers from 1 to 10}
   WHILE x < 11 DO
   BEGIN
      CALL square;
      x := x + 1
   END
END.

There are not very many nice features here. I haven't even included I/O operations, but we'll deal with that later.

Our first job is to get this string of text into a string of important tokens. We can leave out whitespace and comments (the bits inside the {}). Now, how do we go about turning a string of text from a file into a string of tokens? We need a lexing program.

Building a Lexer

The main.rs file:

The first subtask is to make a program that can read a file. Since I want to make a nice, professional-looking program, I'll use the Clap CLI builder. This will let me run my program like so:

pl0rs -f path/to/source_code.pl0 -o path/to/target

Right now, we only need the -f part, for the file we're reading.

use std::path::PathBuf;
use std::process::exit;

use clap::{Parser, Subcommand};

use pl0rs::{self, lexer::lex, read_file, COMPILER_VERSION};

// Borrowed from Clap's boilerplate tutorial...
#[derive(Parser)]
#[command(version, about, long_about = None)]
struct Cli {
    /// Path to the input file
    #[arg(short, long, value_name = "FILE")]
    file: Option<PathBuf>,

    /// Turn debugging information on
    #[arg(short, long, action = clap::ArgAction::Count)]
    debug: u8,
}

fn main() {
    let cli = Cli::parse();
    let mut file = String::new();
    let mut state = pl0rs::State::default();

    // Print compiler version, &c.
    println!("pl0rs -- PL/0 Compiler version {}", COMPILER_VERSION);
    println!("(c) Ethan Barry, 2024. GPLv3 licensed.");

    // Open the file and pass it into the file string.
    // Program can terminate here with an error code.
    if let Some(file_path) = &cli.file {
        if let Ok(file_string) = read_file(file_path) {
            file = file_string;
        } else {
            eprintln!("Error: No such file found.");
            exit(1);
        }
    }

    match cli.debug {
        0 => {
            println!("Debug mode is off.");
        }
        1 => {
            println!("Debug mode is on.");
            state.debug = true;
        }
        _ => {
            println!("Don't be crazy. Defaulting to debug mode on.");
            state.debug = true;
        }
    }

    // Call the other modules.
    match lex(&mut state, &file) {
        Ok(res) => {
            println!("Lexer succeeded.");
        }
        Err(e) => {
            {
                eprintln!("{e}");
                exit(1);
            } // A returned lexer error.
        }
    }
}

That's a minimum working example. I also want to talk about project organization for a second. We know that this will be a complex program. It will need to be tested. Unit tests, declared in mod test and annotated with #[test] are not enough to validate an entire compiler. We also need integration tests. An integration test acts like an external application calling your program. It should be able to test all the functionality, by essentially acting as another fn main() in the codebase.

That's why you should organize the project like this:

┌ src
├─ main.rs
├─ lib.rs
├─ lexer.rs
╰ tests
Cargo.toml
...

Your main.rs contains the glue logic that calls the real code in lib.rs. This way, your integration tests, in tests, can also call your library and test it against whatever test battery you create. It's a more advanced way of organizing the project.

The lib.rs file:

Now for the library where our compiler's functionality will live. We need to handle some housekeeping first. Let's declare a const VERSION so we can print it out when we want. All the real compilers display version numbers! Let's also declare our first module, the lexer.

//! File lib.rs
use std::{fs::File, io::Read, path::Path, process::exit};

pub mod lexer;

pub const COMPILER_VERSION: &str = "0.1.0";

pub struct State {
    pub debug: bool,
    pub line: u32,
}

impl Default for State {
    fn default() -> Self {
        Self {
            debug: false,
            line: 1,
        }
    }
}

I also created a struct called State. It will store whether we're in debug mode, and what line we're on. It's a good idea to impl Default whenever we create a struct we want to use. I think this trait can also be derived, but this way was more fun.

Next, here's a function to read in a file and return it as a string.

//! File lib.rs
/// Can cause program termination with an error code.
pub fn read_file(filename: &Path) -> Result<String, String> {
    let path = Path::new(filename);

    // TODO: Eliminate unwrap.
    if !path.extension().unwrap_or_default().eq("pl0") {
        eprintln!("Error: File must have a .pl0 extension.");
        exit(1);
    }

    let mut file = match File::open(path) {
        Ok(file) => file,
        Err(e) => return Err(format!("Couldn't open file: {e}")),
    };

    let mut contents = String::new();
    match file.read_to_string(&mut contents) {
        Ok(_) => Ok(contents),
        Err(e) => Err(format!("Couldn't read file: {e}")),
    }
}

I'm enforcing that the file has a .pl0 extension, but that isn't really necessary, just fun. Once we've read in the file, we need to represent different tokens somehow. Dr. Callahan used #defines in his C program, but I think an enum will work nicely for Rust. The different types of tokens are these:

//! File lib.rs
#[derive(Debug, Clone)]
pub enum Token {
    Ident(String),
    Number(i64),
    Const,
    Var,
    Procedure,
    Call,
    Begin,
    End,
    If,
    Then,
    While,
    Do,
    Odd,
    Dot,
    Equal,
    Comma,
    Semicolon,
    Assign,
    Hash,
    LessThan,
    GreaterThan,
    Plus,
    Minus,
    Multiply,
    Divide,
    LParen,
    RParen,
}

impl PartialEq for Token {
    fn eq(&self, other: &Self) -> bool {
        std::mem::discriminant(self) == std::mem::discriminant(other)
    }
}

impl std::fmt::Display for Token {
    // Needs a better format method.
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(f, "{:?}", self)
    }
}

Each keyword, operator, and symbol has its own unique enum variant. Identifiers and numbers also get one. Here's where you choose your integer size; I chose i64, which means that I'll always try to parse a number in PL/0 source code into a 64-bit signed integer.

Brief funny story: I've done impl Display here by using a debug formatter. This is so I don't have to type {:?} every time I want to print a Token. Originally, I had used {}, which I assumed would just work. However, I kept getting stack overflow errors. I didn't realize it right away, but by calling write!(f, "{}", self) I was recursively calling fmt() and overflowing the stack. Just thought that was entertaining...

The lexer.rs file:

We have all the boilerplate. Now, we need to lex everything into the right sort of token. This will mean going through the file, character by character, and deciding what to do on the fly. An iterator would be nice!

Our lexer's primary function will take a string containing our PL/0 code, and a State that we can modify. It should return a Vec<Token>, but that might fail, so we'll return a Result<Vec<Token>, String> instead.

//! File lexer.rs

pub fn lex(state: &mut State, file: &str) -> Result<Vec<Token>, String> {
    let mut chars = file.chars().peekable();
    let mut comment = String::new();
    let mut tokens: Vec<Token> = vec![];

    // TODO!

    Ok(tokens)
}

First, let's process comments. In PL/0, comments are closed with braces, like this:

{this is a comment}

When we see an opening brace, we'll skip through all the characters until we find the closing brace that matches it. If there is none, it's an error. We'll add this after the TODO! comment.

'lexer: loop {
    if let Some(c) = chars.peek() {
        if (*c).eq(&'{') {
            chars.next(); // Consume the opening brace
            'comment: for c in chars.by_ref() {
                if c == '}' {
                    break 'comment;
                }
                if c.eq(&'\n') {
                    state.line += 1;
                }
                comment.push(c);
            }
        }
    } else {
        break 'lexer; // No more characters.
    }
}

There! We simply push all the characters in a comment to the comment string, in case we need them. Also be sure to increment the line count; we'll need that later.

The next trick is to discard whitespace between tokens. We can add some code to handle that in the same way we handle comments. After we've handled whitespace and comments, everything remaining is syntactically significant. A nested set of if-else statements will handle that. The entire lex() function is below.

//! File lexer.rs
pub fn lex(state: &mut State, file: &str) -> Result<Vec<Token>, String> {
    let mut chars = file.chars().peekable();
    let mut comment = String::new();
    let mut tokens: Vec<Token> = vec![];

    'lexer: loop {
        if let Some(c) = chars.peek() {
            if (*c).eq(&'{') {
                chars.next(); // Consume the opening brace
                'comment: for c in chars.by_ref() {
                    if c == '}' {
                        break 'comment;
                    }
                    comment.push(c);
                }
            } else if (*c).is_whitespace() {
                if c.eq(&'\n') {
                    state.line += 1;
                }
                chars.next(); // Consume the whitespace.
            } else if (*c).is_alphabetic() || (*c).eq(&'_') {
                let token = identifier(&mut chars, state)?;
                tokens.push(token);
            } else if (*c).is_numeric() {
                let token = number(&mut chars, state)?;
                tokens.push(token);
            } else if (*c).eq(&':') {
                let token = assignment(&mut chars, state)?;
                tokens.push(token);
            } else {
                let token = match *c {
                    '.' => Token::Dot,
                    '=' => Token::Equal,
                    ',' => Token::Comma,
                    ';' => Token::Semicolon,
                    '#' => Token::Hash,
                    '<' => Token::LessThan,
                    '>' => Token::GreaterThan,
                    '+' => Token::Plus,
                    '-' => Token::Minus,
                    '*' => Token::Multiply,
                    '/' => Token::Divide,
                    '(' => Token::LParen,
                    ')' => Token::RParen,
                    _ => {
                        return Err(format!("Unknown token on line {}", state.line));
                    }
                };
                tokens.push(token);
                chars.next();
            }
        } else {
            break 'lexer; // No more characters
        }
    }

    Ok(tokens)
}

A diagram might make the code easier to understand. I'll draw it as a state machine!

Any of the last four boxes on the right add a token to the Vec<Tokens> that we're storing the tokenized code in. Also note the assignment() block. It is called when the lexer finds a colon. In the PL/0 grammar, the assignment operator is :=, and the only place a colon can occur is in the assignment operator. So, if we find a colon without an equals sign, it's an error. Let's write that method now.

//! File lexer.rs
pub fn assignment(chars: &mut Peekable<Chars<'_>>, state: &mut State) -> Result<Token, String> {
    chars.next(); // Consume the ':' character.

    if let Some(c) = chars.peek() {
        if (*c).eq(&'=') {
            chars.next();
            Ok(Token::Assign)
        } else {
            Err(format!("Unknown token ':' on line {}", state.line))
        }
    } else {
        Err(format!("Unterminated assignment on line {}", state.line))
    }
}

This function is pretty straightforward. The if let syntax is a clean way to deal with the iterator. This method advances the iterator right away, because it is only ever called on a ':'. (If it were called somewhere besides a colon, this would be a logic error. You can make a note of that and move on.)

Words and Numbers

Now there are just two functions left to implement, number() and identifier(). These both need to take one character at a time, until they encounter whitespace, or something that isn't legally a number or identifier. Dr. Callahan suggested extending PL/0 to allow breaking up large integers with underscores, like we can do in modern languages. We'll add that as well.

Also notice that the function lex() calls identifier() whenever it sees an alphabetic character or underscore. This includes keywords and variable names. It will have to distinguish them, and return the correct Token enum variant.

//! File lexer.rs
pub fn identifier(chars: &mut Peekable<Chars<'_>>, state: &mut State) -> Result<Token, String> {
    let mut idt = String::new();

    loop {
        // Push characters onto idt one by one.
        if let Some(c) = chars.peek() {
            if (*c).is_alphanumeric() || (*c).eq(&'_') {
                idt.push(*c);
                chars.next();
            } else {
                break;
            }
        } else {
            // Was a None.
            return Err(format!("Unterminated identifier on line {}", state.line));
        }
    }

    let token = match idt.to_lowercase().as_str() {
        "const" => Token::Const,
        "var" => Token::Var,
        "procedure" => Token::Procedure,
        "call" => Token::Call,
        "begin" => Token::Begin,
        "end" => Token::End,
        "if" => Token::If,
        "then" => Token::Then,
        "while" => Token::While,
        "do" => Token::Do,
        "odd" => Token::Odd,
        _ => Token::Ident(idt),
    };

    Ok(token)
}

pub fn number(chars: &mut Peekable<Chars<'_>>, state: &mut State) -> Result<Token, String> {
    let mut num = String::new();
    loop {
        if let Some(c) = chars.peek() {
            if (*c).is_numeric() || (*c).eq(&'_') {
                num.push(*c);
                chars.next();
            } else {
                break;
            }
        } else {
            // Was a None.
            return Err(format!("Unterminated number on line {}", state.line));
        }
    }

    if let Ok(res) = num.parse::<i64>() {
        Ok(Token::Number(res))
    } else {
        Err(format!("Invalid number at line {}", state.line))
    }
}

Notice that identifier() allows variables with underscores or numbers after an initial letter or underscore, so x_0 and _1000 are legal variable names. number() parses the complete number string into an i64, as we decided earlier.

At the risk of making the article too long for anyone to read, I put the complete lexer.rs file below, as written up to this point.

//! File lexer.rs (complete)
use std::{iter::Peekable, str::Chars};

use crate::{State, Token};

/// ## Preconditions
///
/// - Must only be called for characters which can begin an identifier.
pub fn identifier(chars: &mut Peekable<Chars<'_>>, state: &mut State) -> Result<Token, String> {
    let mut idt = String::new();

    loop {
        // Push characters onto idt one by one.
        if let Some(c) = chars.peek() {
            if (*c).is_alphanumeric() || (*c).eq(&'_') {
                idt.push(*c);
                chars.next();
            } else {
                break;
            }
        } else {
            // Was a None.
            return Err(format!("Unterminated identifier on line {}", state.line));
        }
    }

    let token = match idt.to_lowercase().as_str() {
        "const" => Token::Const,
        "var" => Token::Var,
        "procedure" => Token::Procedure,
        "call" => Token::Call,
        "begin" => Token::Begin,
        "end" => Token::End,
        "if" => Token::If,
        "then" => Token::Then,
        "while" => Token::While,
        "do" => Token::Do,
        "odd" => Token::Odd,
        _ => Token::Ident(idt),
    };

    Ok(token)
}

pub fn number(chars: &mut Peekable<Chars<'_>>, state: &mut State) -> Result<Token, String> {
    let mut num = String::new();
    loop {
        if let Some(c) = chars.peek() {
            if (*c).is_numeric() || (*c).eq(&'_') {
                num.push(*c);
                chars.next();
            } else {
                break;
            }
        } else {
            // Was a None.
            return Err(format!("Unterminated number on line {}", state.line));
        }
    }

    if let Ok(res) = num.parse::<i64>() {
        Ok(Token::Number(res))
    } else {
        Err(format!("Invalid number at line {}", state.line))
    }
}

pub fn assignment(chars: &mut Peekable<Chars<'_>>, state: &mut State) -> Result<Token, String> {
    chars.next(); // Consume the ':' character.

    if let Some(c) = chars.peek() {
        if (*c).eq(&'=') {
            chars.next();
            Ok(Token::Assign)
        } else {
            Err(format!("Unknown token ':' on line {}", state.line))
        }
    } else {
        Err(format!("Unterminated assignment on line {}", state.line))
    }
}

pub fn lex(state: &mut State, file: &str) -> Result<Vec<Token>, String> {
    let mut chars = file.chars().peekable();
    let mut comment = String::new();
    let mut tokens: Vec<Token> = vec![];

    'lexer: loop {
        if let Some(c) = chars.peek() {
            if (*c).eq(&'{') {
                chars.next(); // Consume the opening brace
                'comment: for c in chars.by_ref() {
                    if c == '}' {
                        break 'comment;
                    }
                    if c.eq(&'\n') {
                        state.line += 1;
                    }
                    comment.push(c);
                }
            } else if (*c).is_whitespace() {
                if c.eq(&'\n') {
                    state.line += 1;
                }
                chars.next(); // Consume the whitespace.
            } else if (*c).is_alphabetic() || (*c).eq(&'_') {
                let token = identifier(&mut chars, state)?;
                tokens.push(token);
            } else if (*c).is_numeric() {
                let token = number(&mut chars, state)?;
                tokens.push(token);
            } else if (*c).eq(&':') {
                let token = assignment(&mut chars, state)?;
                tokens.push(token);
            } else {
                let token = match *c {
                    '.' => Token::Dot,
                    '=' => Token::Equal,
                    ',' => Token::Comma,
                    ';' => Token::Semicolon,
                    '#' => Token::Hash,
                    '<' => Token::LessThan,
                    '>' => Token::GreaterThan,
                    '+' => Token::Plus,
                    '-' => Token::Minus,
                    '*' => Token::Multiply,
                    '/' => Token::Divide,
                    '(' => Token::LParen,
                    ')' => Token::RParen,
                    _ => {
                        return Err(format!("Unknown token on line {}", state.line));
                    }
                };
                tokens.push(token);
                chars.next();
            }
        } else {
            break 'lexer; // No more characters
        }
    }

    Ok(tokens)
}

We have a fully functional lexer for our compiler, that tokenizes and understands PL/0! Congratulations! Our next steps come in the next post in this series, where we'll parse the Vec<Token> that our lexer spits out.

As always, thank you for reading! If you have suggestions or comments, please write me at ethan.barry@howdytx.technology. Have a great day!

Interpolation is a very common problem. You need to interpolate any time you have:

• A table of values, in larger increments than you'd like,
• A set of points on some unknown function,
• Known values some time apart, and you'd like to know what's in between.

Actually, these are all the same problem. Let's explore some solutions!

Bad Interpolation

There are a lot of techniques for interpolation. 90% of them are guaranteed to be bad, but I think it's a good idea to show why.

In real life, when I'm trying to interpolate something (maybe on the occasions I mess around with a cookbook) I try to "eyeball" it. For the amount of salt in a fractionally-scaled recipe, that's okay. (There is no such thing as too much salt. Fight me.) But for practical applications, we need to do better.

Back in The Old Days™, people had to consult log tables and/or slide rules for values of functions. This still happens today in astronomy, where orbital values come in tables called ephemerides, or almanacs, like these published by the United States Naval Observatory. If you want something between the given values, you need to interpolate.

Linear Interpolation

The simplest (formal) method is to fit the two points nearest your x-value to a line. If we know \(sin(x)\) at two points, say \(x = 0\) & \(x = 2\), then we might interpolate it like this:

...which is objectively terrible. It could be better if we had a smaller increment than 0 to 2, but linear interpolation is inaccurate unless the thing it's interpolating is fairly linear itself. The upside is that it's easy to compute.

Smooth(er) Interpolation With Cosine

We can get a smoother curve by using weights on the interpolation function. We need something that produces a curve anywhere we like, so let's choose cosine. Now, given two points, we'll say: \[y_i \approx y_1 (1-\frac{1}{2}(1-cos(\pi \frac{x - x_1}{x_1-x_2})) + y_2 (\frac{1}{2}(1-cos(\pi \frac{x - x_1}{x_1-x_2}))\]

That's a mess. The \(\frac{x - x_1}{x_1-x_2}\) terms are our weight. When \(x\) is equal to \(x_1\), the term is zero, and all the weight is on \(y_2\). Conversely when \(x\) equals \(x_2\), the weight is one, and it's all on \(y_1\). Here's our visual example with the cosine interpolation added.

Notice that all this has done is smooth our linear function. It's still not precise. It would make a nice curve between multiple points, but let's try for something better.

Lagrange Polynomials and Lagrangian Interpolation

Let's be more ambitious. Say we have a set of \(n\) points (you decide how many) and we want to fit a polynomial to them. There will be many polynomials of degree \(> n\) that go through them, but only one of degree \(< n\). Usually, this polynomial's degree is \(n-1\). If the points happen to lie on a lower-degree polynomial, then the degree \(n-1\) polynomial's leading coefficients are zero.

Given our set of points \[(x_1, y_1), (x_2, y_2), ..., (x_i, y_i), ..., (x_n, y_n)\] the polynomial we want is this: \[y = \sum_{i=1}^n y_iL_i(x)\] where \(L_i(x)\) is the \(i\)th Lagrange polynomial, given by: \[L_i(x) = \prod_{j=1,j\ne i}^n\frac{x-x_j}{x_i-x_j}\]

In other words, the first Lagrange polynomial is: \[L_1(x) = \frac{(x-x_2)(x-x_3)...(x-x_n)}{(x_1-x_2)(x_1-x_3)...(x_1-x_n)}\]

We have a polynomial that goes through all the points, and it is also the polynomial of least degree that does so. (That's important for maximum accuracy.) The linked textbook walks through this topic in more detail, but I'll jump straight to my code.

To get a computer to do this for us, we'll first need to pull in our table of values, which is a vector of points in Computer-Science-ese, and evaluate the \(i\)th Lagrange polynomial at the \(x\)-value we want to interpolate, for all \(i\) in \([0,n]\). This just amounts to a couple of for-loops.

I used a janky iterator to exclude \(i\) in the internal loop. (Recall the provision \(j\ne i\) in our Lagrange polynomial definition.) I'm sure there's a better way to do this, but that's all I could come up with at the time...

fn lagrangian_interpolation(table: Table, xval: f64) -> f64 {
    let n = table.len();
    let mut sum = 0.0;

    for i in 0..n {
        if let Some(Point(x_i, y_i)) = table.pt_at(i) {
            // Evaluate the ith Lagrange poly'n. at xval and add it to sum.
            let mut product = 1.0;

            // Iterator grabs all but i.
            for j in (0..n)
                .enumerate()
                .filter(|&(pos, _)| (pos != i))
                .map(|(_, e)| e)
            {
                let x_j = table.pt_at(j).unwrap().0;
                product *= (xval - x_j) / (x_i - x_j);
            } // product is now the value of L_i(xval).

            sum += y_i * product;
        } else {
            break;
        }
    }

    sum
}

I set up my main function to print out a tab-separated list of points, and redirected this to a file. Then, I used gnuplot to plot these points in yellow, so we have a comparison of all our methods.

Notice that the yellow tracks the actual function in blue very well, until about 1.5 radians. I ran this test with a table of four values, and began interpolating points with roughly \(x = \frac{\pi}{6}\), in increments of \(0.01\), for 1000 iterations. It's extrapolating past my last tabulated value of \((1.5707963, 1.0)\), or \((\frac{\pi}{2}, 1)\), and doing an acceptable job for a little while.

The reader is welcome to try fixing the table so the interpolation is more accurate over the interval. All the Rust and gnuplot source code are in my repo here.

This technique is fairly fast, efficient, and is useful in Gaussian quadrature, as well as other modern applications. I hope you find it interesting, if not useful. Either way, let me know what you thought on GitHub, LinkedIn, or by email to ethan.barry@howdytx.technology!

References

Paul Bourke's classic 1999 post. His blog is pure gold; do yourself a favor and bookmark it.

Interpolation methods

Numerical Methods Guy, 2008:

High order interpolation is a bad idea? – Numerical Methods Guy

Numerical Methods Guy Numerical Methods for the STEM Undergraduateautarkaw

This textbook on Celestial Mechanics, where I first learned about the Lagrange interpolation algorithm:

1.11: Fitting a Polynomial to a Set of Points - Lagrange Polynomials and Lagrange Interpolation

LibretextsLibretexts

Many problems in calculus textbooks are contrived. In fact, most of them are. In reality, the sorts of integrals people like engineers and applied mathematicians see are not (practically) solvable analytically. You can't use the tricks from Calculus II to antidifferentiate the integrand, because it's simply too complicated. So what do we do?

The Quest for Quadrature

When integrals can't be antidifferentiated, or when it's too time-consuming to do so, the area under the curve can be approximated by a computer. These numerical, rather than analytical, methods are called "quadrature." (There is an algorithm to either symbolically integrate a function, or else show it cannot be expressed in elementary functions, but this is far beyond this post's scope!)

Perhaps the simplest quadrature algorithm is to draw a line between the function's value at the endpoints and calculate the area under it, using simple geometry. This would be both A) "interpolation to a linear function," and B) very inaccurate. But, as always, we can make smaller slices.

Trapezoidal Quadrature

Suppose we take the interval \([a,b]\) and slice it into many smaller pieces, drawing a line between the function values at each. This would give us a series of trapezoids under the curve. It seems like this would nicely approximate the curve, but it turns out that the trapezoidal method is nearly always worse than the technique we originally learned early in calculus for estimating area with sums of rectangles!

Simpson's Rule

So, if we have a gently curved integrand which is too complicated to find an antiderivative for, like this, \(\int_a^b f(x),dx\), and we want to integrate it, then let's try interpolating it to something we can calculate. We'll evaluate the function three times, once at \(a\), then at \(b\), and finally at the midpoint.

This curve can be seen as three points on a piece-wise quadratic. Instead of interpolating to a linear function, we've interpolated to a quadratic, using one extra function evaluation. And the area under this quadratic is approximately: \[\frac{1}{3}(f(a)+4f(c)+f(b))\] where \(c\) is the midpoint of \(a\) and \(b\). Then, to gain accuracy, we can split the interval into \(n\) small sections, and interpolate for each one of those, so that:

\[\int_a^b f(x),dx \approx \sum_{1\leq i\leq n} \frac{x_{i+1} - x_i }{6} (f(x_i) + 4f(\frac{x_i + x_{i+1}}{2}) + f(x_{i+1})) \]

All we need is a function that takes a number of intervals to evaluate (the precision) and the bounds. Since I want to pass in any mathematical operation, I'll also take a closure, which is a way of passing a function as an argument to another function. (Yes, functions are data types too!) Here's my Rust implementation:

fn int_simpson<F>(a: f64, b: f64, n: i32, integrand: F) -> f64
where
    F: Fn(f64) -> f64,
{
    // We want to sum from a to b, in n subintervals, where the width delta_x is (a - b) / n.
    // Each subinterval is [a, a+delta_x] up to b.
    let delta_x = (b - a) / n as f64;
    let mut result = 0.0_f64;
    for i in 1..=n {
        result += delta_x
            * (integrand(a + (i as f64 - 1_f64) * delta_x)
                + 4_f64 * integrand(a - delta_x / 2_f64 + i as f64 * delta_x)
                + integrand(a + i as f64 * delta_x))
            / 6_f64;
    }
    result
}

The complete program can be found here, along with an example of using it, and even asserting that the floating-point result is accurate to a certain precision!

Other Methods

I think it's interesting that we can see the rectangle rule, trapezoid rule, and Simpson's rule all have something in common. They're interpolating to polynomials (zeroeth, first, and second degrees, respectively) and estimating. There are certainly better methods available, especially adaptive quadratures, for functions with large high-order derivatives, but this is a quick-and-dirty result that achieves enough accuracy for most use cases. I hope you enjoyed!

Footnote: I've added LaTeX formatting to this post using the excellent KaTeX JS package. Simply add it to your site header, and it will typeset any LaTeX markup it sees. Instructions are here.

Java is the de facto instructional language of many universities. It's widespread, it runs on most platforms, it makes a good case study in OOP concepts, and there are a lot of teachers available.

But in many cases, the Java that's taught in the classroom isn't much like production Java. We learn the oldest syntax, and the oldest design patterns, that accomplish the task. It makes typing Java feel like drudgery compared to a more modern language like Rust or Go. Seriously, who wants to type this mess:

MySpecialClassName myClass = new MySpecialClassName();

to instantiate a variable?

The var Keyword

In more modern languages, we have type inference. The compiler can tell what a type a variable holds based on what's assigned to it. But wait, Java actually has this too! Instead we can write this:

var myClass = new MySpecialClassName();

which is much more streamlined. var is a special "reserved type name" that was added in Java 10, to make declarations easier to write. So why don't we learn this in university? Personally, the Java courses I've taken have been Java 7 compatible, or older. Java 7 was released in 2011. A lot has changed since then.

The Optional Type

In Rust, there's a special enum included in the compiler's prelude, called Option. It has two variants, Some() and None. The convention is to use it when a function may or may not return something, such as a query to a database.

Java added a similar concept with the Optional<?> in Java 8 (released in 2014). It's a box that might hold a value, or a null pointer. The biggest advantage is that by making this return type explicit, developers have to handle the possibility of a null pointer when it might be returned from a method. Null-checking is simplified. This:

Apple apple = mightBeNull("granny smith");

if (apple != null) {
    System.out.println(apple.toString());
} else {
    System.out.println("apple was null");
}

becomes this (if we use our fancy new type inference):

Optional<Apple> optApple = mightBeNull("granny smith");

var appleVariety = optApple.orElse("No apple in this box...");

System.out.println(appleVariety); // This is safe now!

Record Classes

A common pattern in Java is using an object to hold related data, usually as it's being moved from a source to a consumer. If you had a database of apples, you might represent each apple as a class, like this:

public class Apple {
    private String variety;
    private float pricePerPound;
    
    public String getVariety() {
        return this.variety;
    }
    
    public float getPrice() {
        return this.pricePerPound;
    }
    
    public Apple(String v, float p) {
        this.variety = v;
        this.pricePerPound = p;
    }
    
    @Override
    public String toString() {
        return "Apple(variety=" + this.getVariety() + ", price=" + this.getPrice();
    }
}

That's a lot of code for one apple! Java 14 added a solution in the form of the record class. With modern Java, our Apple is just:

public record AppleRecord(String variety, float pricePerPound) {}

That's literally it.

It implements getters automatically; just call AppleRecord.variety() to get the variety, for example. By default, the fields are immutable, so a record is instantiated like this:

var apple = new AppleRecord("pink lady", 0.20);

If you need custom logic, you can override these default methods in the record's block. Easy!

Conclusion

So why aren't we being taught these concepts? Is it the result of a one-size-fits-all approach to CS education? Are the instructors even familiar with these later developments in the Java ecosystem? I don't know.

What I do know is that being able to think in terms of Optionals, and reduce boilerplate code, will make Java feel like a whole new language for me!

I recently watched a Numberphile video on the Catalan series, which happens to count lot of different things, from partitions of sets, to divisions of polygons, to unique binary trees of a certain order. The video is fascinating and well worth a watch.

One other thing the Catalan numbers count is the number of possible Dyck words of a given semi-length. A Dyck word is a word in a binary (having two elements) alphabet, where no prefix of the word has more of the second element than the first. That's the complicated definition. The simple definition is a balanced set of parentheses, like so: "()()(())" or "((()))()". This should jump out as interesting to any developer - we use parentheses all day!

So, the Catalan numbers count the Dyck words of a given semi-length. For example, there are 42 Dyck words of semi-length 5, since the fifth Catalan number is 42. But how can we generate them all? And more importantly, how can we do it..... in Rust?

It turns out that defining the Catalan numbers and Dyck words in terms of recursion is easy. But with recursion, there's a limit to how far our program can go. Eventually, the call stack won't be able to hold all our recursive function calls, and our program will terminate with a stack overflow. Even if we use the heap with more advanced features, we will slowly eat into memory and run up against the computer's physical limits. (In my case, this is 24 GiB of DDR3 in a mobile workstation from 2013.)

We need a better algorithm than the recursive solution. We need a perfect algorithm.

...and there is one! It is entirely possible to convert recursive definitions and functions into iterative ones, and someone (Dr. Knuth, it seems) has already done so for this particular problem. After a bit of digging, I was able to find a paper on the subject. [EDIT: I was missing this link for a while after I had to rebuild this site, but it's been fixed. My apologies to the author!]

The algorithm goes as close as possible to the metal by representing opening parentheses as a 1 bit, and closing parentheses as a 0 bit. In the abstract you can see the algorithm in C-like pseudo-code. It's simply bitwise integer math. Fast, efficient, loop-less, and branch-less. The only trick is that it requires the "previous" Dyck word to act upon.

You can generate the "next" Dyck word in the sequence by flipping pairs around. Taking the Dyck word of the form "()()()()()()" or "101010101010" to be the least, and the word of form "(((((())))))" or "111111000000" to be the greatest, we can run this algorithm on its output until we reach the greatest Dyck word of the given semi-length. The integer passed in to the algorithm must represent a valid Dyck word, or it's a precondition violation.

The full source of my program is here, but I'll reproduce the Rust port of the algorithm below.

fn next_dyck_word_bin(w: u64) -> u128 {
    // Calculate a, potentially using two's complement for negation.
    let a = w & (!w + 1);

    // Calculate b directly.
    let b = w + a;

    // Calculate c using bitwise XOR.
    let c = w ^ b;

    // Shift c right by 2, ensuring unsigned division.
    let c = (c as u128 / a as u128) >> 2; // Cast to u128 for division.

    // Add 1 to c, handling potential overflow.
    let c = c.wrapping_add(1);

    // Apply mask and bitwise OR with b.
    let mask = 0xaaaa_aaaa_aaaa_aaaa_aaaa_aaaa_aaaa_aaaa; // I think this is long enough...

    // Return the resulting Dyck word as a u128.
    ((c * c - 1) & mask) | b as u128
}

I've run the program for semi-length 20, whose Catalan number is 6,564,120,420. As you can see, the series grows quickly. It ran for 6.5 hours before I terminated it; the output file had reached ~200 GiB and over 4 billion of the 6.5 billion Dyck words, represented as ones and zeroes. But in that time, my memory usage never wavered, and my CPU stayed locked at 100%. It was very interesting to watch in reality what O(1) means - it's an invincible algorithm.

A couple improvements to the rest of the program have been suggested to me; I'll implement them soon. Next, I want to work on some numerical integration techniques, like Simpson's method. This was an exciting peek into the world of combinatorics and theoretical computer science. I'm excited to see more!