Math Notes 2026

I enjoy tinkering with math. Unfortunately, I often forget the little things I discover. On this page you will find all the tidbits I found interesting enough to record.

I don't strive for academic perfection, however I also do not want to spread blatant errors. If you notice any errors in these notes, please email me.

See all math notes

Newton’s method allows finding the ‘root’ of a function (where it is zero), given conditions expressed in the wikipedia article no doubt.

More than this however, it can be used to find the inverse of a function. The common example is $\sqrt{2}$. If you know $x$ where $x^2-2=0$, then you know $\sqrt{2}$.

The rough gist of the algorithm is as follows.

1. Take a function like $f(x) = x^2-2$
2. Make an arbitrary guess on the value of $x$ where $f(x)=0$, $g_0$
3. Evaluate $f(x)$ at $g_0$ and find the line tangent to $f(x)$ at $g_0$
4. Find the x-intercept of the tangent line. This x-intercept value is $g_1$
5. Repeat steps 2…5 (incrementing $n$ throughout) until a satisfactory precision is obtained.

Algebraically, we start with the function and its derivative:

And make a guess $g_0$ and create the tangent line at that guess location on $f(x)$

find $b$ by evaluating at $g_0$:

Which finishes the line equation as:

Next, find the x-coordinate by evaluating at $y=0$

As stated above, this $x$ value is actually $g_1$, or generally

Going back the case of $\sqrt{2}$, we substitute in the function

With $g_n=5$ we see the first four approximations as:

Generally…

This process is easily generalized to finding inverses instead of roots — just as in the $\sqrt{2}$ example.

$f(x)$ was previously defined as $x^2-2$, but more generically is described as $f(x)=x^2$. Then, inside the recursive equation, the $-2$ is placed outside the function

The above $g_{n+1}$ approximates $f^{-1}(2)$. the $2$ is replaced with some constant $c$ to find $f^{-1}(c)$:

The generic case for $\sqrt{c}$ then is:

Or

One way to loosely verify this is to: consider that when $n$ approaches $\infty$, the differences between $g_n$ and $g_{n+1}$ will grow increasingly smaller (assuming the limit does indeed converge) and indeed will be come so close that we can consider them to be equal. So, we set $g_{n+1} = g_n = x$ and solve:

Other inverses

This method can be extended to other functions. Some examples:

Approximate Functions

By expanding the recursion a few steps, we can create an obnoxious approximation for $\ln(x)$:

And one for $\sqrt{x}$:

In a way, this approximation for $\sqrt{x}$ is quite reasonably constructed. It has plenty of detail around the origin, but as it goes further the higher-order polynomial components divide out and all that is left is the slope asymptote.

These examples were generated with the Python library SymPy. (I got exhausted of writing it by hand around third recursion)

from sympy import symbols, simplify, print_latex

x = symbols("x")

expr = 1
for i in range(0, 3):
    expr = 0.5 * expr + x / (2 * expr)

print_latex(simplify(expr))

Suppose $f(x)=\sin(x)$, then $f(\pi) = 0$. Using Newton’s method, we arrive at the recursion:

Which gives π accurate to 100 digits in only 5 iterations.

from sympy import N, tan, pi

expr = 3
for i in range(0, 5):
    expr = expr - tan(expr)

print(N(expr, 100))
print(N(pi - expr, 100))

The behavior of $g_{n+1} = g_n - \tan(g_n)$ depends highly the initial value $g_0$. This makes some intuitive sense because $\sin(x)$ has many different roots so the sequence would likely converge to the nearest root to the starting point.

This behavior is clear by plotting $f(x) = x - \tan(x)$:

Around each multiple of pi, the graph flattens out. to the right of each flat point, the value of the function is decreases slightly. To the left, the value increases slightly. When the function is interpreted as the mapping from the previous value $x$ to the next value $f(x)$, the slight decreases and increases near each multiple of π makes the recursion stably approach the nearest multiple of π. Further away from the multiples of π, the slope becomes steep enough to ‘push’ the next step of the recursion away from the nearest multiple of pi. The steep region is not well-behaved. for these notes, ‘well-behaved’ means that the recursion $g_{n+1} = g_n - \tan(g_n)$ has a value which approaches the multiple of π nearest to $g_0$.

Specifically, let’s look at $x=0\pi$. The well-behaved region will extend an equal distance positive and negative from zero (because the function is symmetric to an extent around the multiples of π) Therefore, we can solve for this distance by finding $g_n = g_{n+1}$, or $x=-(x-tan(x))$.

Or

Where $g(x) = 0$

Plotting this function, we see the zero points are around $\approx\pm1.166$.

I can’t think of any algebra to solve $g(x) = 0$. Instead, we will use Newton’s method to find some ‘radius of stability’ $S$ such that $g(S)=0$

(there’s some humor in here somewhere; I just know it!)

Starting with $g_0 = 1.16$, we see it level out to $S\approx1.16556118$. I do not believe this constant has any simpler symbolic forms.

from sympy import N, tan, sec

expr = 1.16
for i in range(0, 15):
    expr = expr - (tan(expr) - 2 * expr) / (sec(expr) * tan(expr) - 2)
    # evaluate algebra to numeric value each step to prevent
    # algebraic explosion. Use 200 digits to avoid precision loss
    expr = N(expr, 200)
print(N(expr, 15))

So, around $0\pi$ this well-behaved region is described as

And in general is:

Where:

When $g_n$ is outside the well-behaved region, $g_{n+1}$ will skip to be nearby an entirely different multiple of pi. That makes me curious: is there some value of $g_{0}$ such that all subsequent $g_n$ values will keep skipping away and never converge?

Well, there’s the solution $g_{0} = S$ where the iterations bounce back and forth, but that’s more of a bounce than a skip. In theory, there exists some value $g_n$ such that $g_{n+1}=g_{n}+\pi$. A $g_n$ of this description will never converge.

Of course, this same logic applies for the situation where $g_{n+1} = g_{n} \pm q\pi$ where $q\in\mathbb{Z}$ and $q\ne 0$. So, generally:

will diverge. For example:

I was quite interested in describing the nature of this non-well-behaved region. However, the more I looked, the more detail I saw. Slight changes in input create vast differences in output, nudging the inputs even slightly creates entirely different outputs. This detail seems to persist no matter how precise the inputs are.

This all stinks of a fractal.

After discovering that the details were infinite, I stopped trying to directly define the results of the function. However, one last detail I noticed which I have not yet explained is the converging values around $-\pi/2$ (and all other midpoints between adjacent well-behaved areas)

As seen in the graph, the convergent values around $-\pi/2$ generally tend to be the multiple of π less than a certain hyperbola. Fascinating.

I recommend viewing this graph in Desmos yourself. I’ve included the code to recreate it below (make sure to make the range on the y-axis quite large).

f\left(x\right)=x-\tan\left(x\right)

x=-\frac{\pi}{2}

y=f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(f\left(x\right)\right)\right)\right)\right)\right)\right)\right)\right)\right)\right)\right)\right)\right)\right)

\frac{1}{\left(x+\frac{\pi}{2}\right)}

\operatorname{floor}\left(\frac{\frac{1}{\left(x+\frac{\pi}{2}\right)}}{\pi}\right)\pi

Perhaps this is quite similar to 3Blue1Brown: Newton’s Fractal. It has been a while since I have seen that video.

from sympy import N, pi, floor
best = (10000, 1)

for i in range(1, 1_000_000):
    val = N(pi * i, 100)
    lz = val - floor(val)
    print(i, lz, val)
    lz = lz.evalf()
    if lz < best[0]:
        best = (lz, i)

print(best)

If I were to trust google’s calculator… then I’d say π is rational.

Two equivalent expressions of a standard decimal number:

A mixed-radix number is somewhat like having a different base for each digit. The primary objective of this paper is stating mixed-radix representation’s of $e$ and $\pi$, and converting them to more useful bases.

Example of mixed-radix number:

Which expands to

When $a_i=1$ for all $i$, this is the well-known power series for $e$.

Euler’s number spigot

A. H. J. Sale doi 10.1093/comjnl/11.2.229.

Well-known $e$ power series:

Given in mixed-radix form:

Also can be stated as a recursion:

The paper then multiplies the fractional part of the expression by $\frac{10}{10}$, and distributes the numerator into the recursion:

The goal of this manipulation is to extract base-10 digits from this mixed-radix form, one digit at a time. We first take this recursion to some specific depth. $5$ in this example.

From right to left, we consecutively reduce the inner values mod the outer divisor, moving out the quotient in the process. An example is more clear:

This method of distributing the divisor, diving with remainder, and factoring is repeated:

This shows clearly that the 2nd digit is $7$. To get further digits, we again multiply the mixed-radix part by $10/10$.

And repeat the process:

And again:

And again:

Which now produces an incorrect result.

Why? Well, it is reasonable to assume that, because we limited the recursion to a finite number of levels, the remaining levels begin to become significant as the powers of $10$ rise.

Looking at the last reasonable result, we see at the end:

Let’s call $(\ldots)$ to be $r$ , we see

Which expands to

Multiplying in the $10^3$ creates a value less than $1$, and so it does not effect the calculations.

However, when the power of ten is increased to $10^4$, we see:

Which is greater than one and will likely cause a deleterious effect when it is missing from the calculations. So, the question is: how many terms does the recursion need to be expanded to to obtain $d$ digits? Instead of thinking about the algorithm as algebraic manipulations on a recursive equation, I find it easier to think of as an array of numbers. The first step of the above example becomes:

Where:

Inside the array, there are 5 slots (columns). Each slot represents a level of recursion. More slots => more digits extractable.

Let $d$ represent the number of digits being extracted, and $s$ represent the number of slots. A calculation using $s$ slots will yield $d$ digits when the following is true:

Infinite sums are cumbersome so instead a simpler form is desired. So long as the value of the sum overestimated, a simpler form will still work.

The following form is valid according to Desmos and WolframAlpha, so I’ll trust those sites on that.

Substituting this simpler expression in:

Solve for digits:

Plotting this, we see

The linear approximation of

Fits well under the entire domain. This approximation becomes poorer as more slots are added, but never overestimates. When $s=70$, $d \le \approx100$

A polynomial or piecewise approximation could provide improvement. The original paper provides a different approximation entirely. I have not reviewed that approximation.

When in one of the previous examples the result of $2.7180\ldots$ was obtained, only 5 slots were used, and so only two digits after the decimal place should have been expected with confidence to be correct.

Back again to thinking of the problem algorithmically:

Where:

To find the second digit, update the table such that:

• $v_i=10r_i$ for all $i$
• $q_i$ is cleared for all $i$

Then apply the same rules, stopping after $i=1$:

And so digit #2 is $1$: $2.71\ldots$.

An implementation in Rust is provided.

fn main() {
    for i in 1..20 {
        println!("{:?}", e(i));
    }
}

fn e(d: u32) -> Vec<u8> {
    let mut digits = vec![0; d as usize];
    let slots = d as usize + 3;
    let mut tmp = vec![10; slots];
    for stop in 0..d as usize {
        let mut carry = 0;
        for i in (0..slots).rev() {
            let n = i as u32 + 2;
            let s = carry + tmp[i];
            let q = s / n;
            let r = s % n;
            tmp[i] = r * 10;
            carry = q;
        }
        digits[stop] = carry as u8;
    }
    digits
}

Back to π

After a slight detour, the study of π is resumed.

Todo….

On the sums of natural numbers

Think about each integer as an area, and draw each component integer of $L_n$ ascending, in a square of dimension $n$:

In this case, $n=4$ and $L_4 = 10$.

Notice the diagonal set of squares from the bottom left to the top right of the greater square. There exist $n$ units on this diagonal. Above the diagonal, there are zero units as part of the greater square. Below the diagonal, all units are part of the greater square.

We know the diagonal has a total area of $n$, but what about the area under the diagonal? If we consider the greater square with the diagonal removed, then this under-area is clearly half of the greater square’s (sans diagonal) area.

Thus the area of the under-area is $\frac{n^2 -n}{2}$ units

To find $L_4$, the diagonal must be added back

Proof By Induction

Base condition:

Following values

Expand left side

Both sides are equal.

On the sums of squares of natural numbers

Or recursively:

(Note, this section follows from a YouTube video by Perfect Yarn)

Thing about each squared integer as a square. Arrange them ascending:

In this case, $n = 3$ and $S_3 = 14$.

Just as with the un-squared natural numbers, a closed form for $S_n$ can be found by interpreting this geometric construction.

If a equation for the quantity of units in each row $R_{j}$ is found, then the total $S_n$ is the sum of the row $R_{j}$ values. Where $j$ is the row number, from bottom up, starting at $0$:

For the bottom row ( $j=0$), $R_{0,3} = L_3$ (where $L_n$ is the sum of natural numbers from 1…n) or

For the following values of $j > 0$, $R_{j}$ will decrease. Each increase of $j$ will cause a subsequent quantity to be removed:

Again, recall that

Substituting in, $R_{j}$ we see:

Factor

The sum of natural numbers of the right of the equation can be simplified to $L_{n-1}$ and the sum of squares on the left can be simplified to $S_{n-1}$

Recall that:

Substitute $S_{n-1}$

Factor $S_n$ to one side

Simplify

Proof By Induction

Base condition

Following values

Simplify

Both sides are equal.

The sum of natural numbers

Can be described as a recursive function

Or

We seek $f(x)$. We will start with a guess-and-check method starting with $f(x) = x^2$.

This guess for $f(x)$ does not satisfy the constraint and must be updated. A more clear method for doing this will be explained later, for now, halve current guess:

Now, the guess for $f(x)$ is a lot closer to satisfying the equation. Update guess again

A Bit more formally

This method, as previously shown, works by incrementally adding and removing polynomials to $f(x)$ until

Where $h(x)$ is some polynomial function. A table has been made to better understand the effect which adding a power of $x$ to the guess for $f(x)$ has on the overall equation.

The above table indicates that, if $x^2$ is added to the guess for $f(x)$, then $f(x)-f(x-1)$ will have $2x-1$ added to it as a result.

Using this table, the method will again be applied to the sum of squares, where $h(x) = x^2$. It is not necessary to use the table, but the algebra becomes pretty murky without it.

The starting guess for $f(x)$ is informed by the highest power of $x$ in $h(x)$. Because $h(x)=x^2$, the first guess will be $f(x) = \frac{1}{3}x^3$. The specific $\frac{1}{3}x^3$ is chosen because it will make $f(x) - f(x-1)$ contain $x^2$ which is of course desirable.

Now that the expression contains $x^2$, the unnecessary components $-x + \frac{1}{3}$ must be removed. Remove $-x$ by adding $\frac{1}{2}x^2$ to the guess for $f(x)$

Next, $-\frac{1}{6}$ is removed by adding $\frac{1}{6}x$ to the guess for $f(x)$

Thus, the sum of squares formula is

Solver Algorithm Implementation

A sample algorithm implementation in Python using NumPy is provided:

import fractions as F
from numpy.polynomial import polynomial as P

# Possibly unexpected behavior:
assert P.polytrim([0]) == [0]

def polyset(p, n, v):
    while len(p) <= n:
        p.append(0)
    p[n] = v

# Given f(x) - f(x-1) = g(x) Where `g` is given value
given = [F.Fraction(0), F.Fraction(1)] # x^2

# Populate table
table = []
base = [F.Fraction(-1), F.Fraction(1)]
tmp = base.copy()
for i in range(len(given) + 1):
    x = tmp[:-1]
    x = P.polymul(x, [-1])
    leading_coef = x[-1]
    table.append((leading_coef, x))
    tmp = P.polymul(tmp, base)
tmp = None
base = None

# move given to left side of equation: f(x)-f(x-1) - g(x) = 0
working = P.polymul(given, -1)
out = []

# Work until a satisfying function $f$ is found
# such that the equation is zero on both sides
while len(working) > 1 or working[-1] != 0:
    high_pwr = len(working) - 1
    hp_c = working[high_pwr]
    assert hp_c != 0
    lc, tp = table[high_pwr]
    factor = hp_c / lc * -1
    polyset(out, high_pwr + 1, factor)
    tp = P.polymul(tp, [factor])
    working = P.polyadd(working, tp)
    working = P.polytrim(working)

print("The sum of: ", P.Polynomial(given))
print("Has the closed form: ", P.Polynomial(out))

Additionally, I’ve ported it to JavaScript so you can use it interactively. All polynomials must be ine expanded form.