Giant Steps

Posted on by John

John Coltrane’s song Giant Steps is known for its unusual and difficult chord changes. Although the chord progressions are complicated, there aren’t that many unique chords, only nine. And there is a simple pattern to the chords; the difficulty comes from the giant steps between the chords.

Giant Steps chords

If you wrap the chromatic scale around a circle like a clock, there is a three-fold symmetry. There is only one type of chord for each root, and the three notes not represented are evenly spaced. And the pattern of the chord types going around the circle is

minor 7th, dominant 7th, major 7th, skip
minor 7th, dominant 7th, major 7th, skip
minor 7th, dominant 7th, major 7th, skip

To be clear, this is not the order of the chords in Giant Steps. It’s the order of the sorted list of chords.

For more details see the video The simplest song that nobody can play.

Related posts

Tritone substitution

Posted on by John

Big moves in roots can correspond to small moves in chords.

Imagine the 12 notes of a chromatic scale arranged around the hours of a clock: C at 12:00, C♯ at 1:00, D at 2:00, etc. The furthest apart two notes can be is 6 half steps, just as the furthest apart two times can be is 6 hours.

Musical clock

An interval of 6 half steps is called a tritone. That’s a common term in jazz. In classical music you’d likely say augmented fourth or diminished fifth. Same thing.

The largest possible movement in roots corresponds to almost the smallest possible movement between chords. Specifically, to go from a dominant seventh chord to another dominant seventh chord whose roots are a tritone apart only requires moving two notes of the chord a half step each.

For example, C and F♯ are a tritone apart, but a C7 chord and a F♯7 chord are very close together. To move from the former to the latter you only need to move two notes a half step.

Musical clock

Replacing a dominant seventh chord with one a tritone away is called a tritone substitution, or just tritone sub. It’s called this for two reasons. The root moves a tritone, but also the tritone inside the chord does not move. In the example above, the third and the seventh of the C7 chord become the seventh and third of the F♯7 chord. On the diagram, the dots at 4:00 and 10:00 don’t move.

Tritone substitutions are a common technique for making basic chord progressions more sophisticated. A common tritone sub is to replace the V of a ii-V-I chord progression, giving a nice chromatic progression in the bass line. For example, in the key of C, a D min – G7– C progression becomes D min – D♭7 – C.

Related posts

Bitcoin mining difficulty

Posted on by John

The previous post looked at the Bitcoin network hash rate, currently around one zettahash per second, i.e. 1021 hashes per second. The difficulty of mining a Bitcoin block adjusts over time to keep the rate of block production relatively constant, around one block every 10 minutes. The plot below shows this in action.

Bitcoin hash rate, difficulty, and ratio of the two

Notice the difficulty graph is more quantized than the hash rate graph. This is because the difficulty changes every 2,016 blocks, or about every two weeks. The number 2016 was chosen to be the number of blocks that would be produced in two weeks if every block took exactly 10 minutes to create.

The ratio of the hash rate to difficulty is basically constant with noise. The noticeable dip in mid 2021 was due to China cracking down on Bitcoin mining. This caused the hash rate to drop suddenly, and it took a while for the difficulty level to be adjusted accordingly.

Mining difficulty

At the current difficulty level, how many hashes would it take to mine a Bitcoin block if there were no competition? How does this compare to the number of hashes the network computes during this time?

To answer these questions, we have to back up a bit. The current mining difficulty is around 1014, but what does that mean?

The original Bitcoin mining task was to produce a hash [1] with 32 leading zeros. On average, this would take 232 attempts. Mining difficulty is defined so that the original mining difficult was 1 and current mining difficulty is proportional to the expected number of hashes needed. So a difficulty of around 1014 means that the expected number of hashes is around

1014 × 232 = 4.3 × 1023.

At one zetahash per second, the number of hashes computed by the entire network over a 10 minute interval would be

1021 × 60 × 10 = 6 × 1023.

So the number of hashes computed by the entire network is only about 40% greater than what would be necessary to mine a block without competition.

Related posts

[1] The hash function used in Bitcoin’s proof of work is double SHA256, i.e. the Bitcoin hash of x is SHA256( SHA256( x ) ). So a single Bitcoin hash consists of two applications of the SHA256 hash function.

Exahash, Zettahash, Yottahash

Posted on by John

When I first heard of cryptographic hash functions, they were called “one-way functions” and seemed like a mild curiosity. I had no idea that one day the world would compute a mind-boggling number of hashes every second.

Because Bitcoin mining requires computing hash functions to solve proof-of-work problems, the world currently computes around 1,000,000,000,000,000,000,000 hashes, one zettahash, per second. Other cryptocurrencies uses hash functions for proof-of-work as well, but they contribute a negligible amount of hashes per second compared to Bitcoin.

The hashrate varies over time because the difficulty of Bitcoin mining is continually adjusted to keep new blocks being produced at about one every ten minutes. As hardware has gotten faster, the difficulty level of mining has gotten higher. When the price of Bitcoin drops and mining becomes less profitable, the difficulty adjusts downward. There are other factors too, and hashrate is variable.

Bitcoin network hashes per second over time

The prefix giga- (109) wasn’t widely known until computer memory and storage entered the gigabyte range. Then the prefix tera- (1012) became familiar when disk drives entered terabyte territory. The prefixes for larger units such as peta- (1015) and exa- (1018) are still not widely known. The prefixes zetta- (1021) and  yotta- (1024) were adopted in 1991 and ronna- (1027) and quetta (1030) were adopted in 2022.

When the Bitcoin hashrate is relatively low, it’s in the range of hundreds of exahashes per second. At the time of writing the hashrate is 1.136 ZH/s, 1.136 × 1021 hashes per second. This puts the hashrate per day in the tens of yottahashes and the number of hashes per year in tens of ronnahashes.

Related posts

10,000,000th Fibonacci number

Posted on by John

I’ve written a couple times about Fibonacci numbers and certificates. Here the certificate is auxiliary data that makes it faster to confirm that the original calculation was correct.

This post puts some timing numbers to this.

I calculated the 10 millionth Fibonacci number using code from this post.

n = 10_000_000    
F = fib_mpmath(n)

This took 150.3 seconds.

As I’ve discussed before, a number f is a Fibonacci number if and only if 5f² ± 4 is a square r². And for the nth Fibonacci number, we can take ± to be positive if n is even and negative if n is odd.

I computed the certificate r with

r = isqrt(5*F**2 + 4 - 8*(n%2))

and this took 65.2 seconds.

Verifying that F is correct with

n = abs(5*F**2 - r**2)
assert(n == 4)

took 3.3 seconds.

About certificates

So in total it took 68.5 seconds to verify that F was correct. But if someone had pre-computed r and handed me F and r I could use r to verify the correctness of F in 3.3 seconds, about 2% of the time it took to compute F.

Sometimes you can get a certificate of correctness for free because it is a byproduct of the problem you’re solving. Other times computing the certificate takes a substantial amount of work. Here computing F and r took about 40% more work than just computing F.

What’s not typical about this example is that the solution provider carries out exactly the same process to create the certificate that someone receiving the solution without a certificate would do. Typically, even if the certificate isn’t free, it does come at a “discount,” i.e. the problem solver could compute the certificate more efficiently than someone given only the solution.

Proving you have the right Fibonacci number

Now suppose I have you the number F above and claim that it is the 10,000,000th Fibonacci number. You carry out the calculations above and say “OK, you’ve convinced me that F is a Fibonacci number, but how do I know it’s the 10,000,000th Fibonacci number?

The 10,000,000th Fibonacci number has 2,089,877 digits. That’s almost enough information to verify that a Fibonacci number is indeed the 10,000,000th Fibonacci number. The log base 10 of the nth Fibonacci number is very nearly

n log10 φ − 0.5 log10 5

based on Binet’s formula. From this we can determine that the nth Fibonacci number has 2,089,877 digits if n = 10,000,000 + k where k equals 0, 1, 2, or 3.

Because

log10 F10,000,000 = 2089876.053014785

and

100.053014785 = 1.129834377783962

we know that the first few digits of F10,000,000 are 11298…. If I give you a 2,089,877 digits that you can prove to be a Fibonacci number, and its first digit is 1, then you know it’s the 10,000,000th Fibonacci number.

You could also verify the number by looking at the last digit. As I wrote about here, the last digits of Fibonacci number have a period of 60. That means the last digit of the 10,000,000th Fibonacci number is the same as the last digit of the 40th Fibonacci number, which is 5. That is enough to rule out the other three possible Fibonacci numbers with 2,089,877 digits.

Computing big, certified Fibonacci numbers

Posted on by John

I’ve written before about computing big Fibonacci numbers, and about creating a certificate to verify a Fibonacci number has been calculated correctly. This post will revisit both, giving a different approach to computing big Fibonacci numbers that produces a certificate along the way.

As I’ve said before, I’m not aware of any practical reason to compute large Fibonacci numbers. However, the process illustrates techniques that are practical for other problems.

The fastest way to compute the nth Fibonacci number for sufficiently large n is Binet’s formula:

Fn = round( φn / √5 )

where φ is the golden ratio. The point where n is “sufficiently large” depends on your hardware and software, but in my experience I found the crossover to be somewhere 1,000 and 10,000.

The problem with Binet’s formula is that it requires extended precision floating point math. You need extra guard digits to make sure the integer part of your result is entirely correct. How many guard digits you’ll need isn’t obvious a priori. This post gave a way of detecting errors, and it could be turned into a method for correcting errors.

But how do we know an error didn’t slip by undetected? This question brings us back to the idea of certifying a Fibonacci number. A number f Fibonacci number if and only if one of 5f2 ± 4 is a perfect square.

Binet’s formula, implemented in finite precision, takes a positive integer n and gives us a number f that is approximately the nth Fibonacci number. Even in low-precision arithmetic, the relative error in the computation is small. And because the ratio of consecutive Fibonacci numbers is approximately φ, an approximation to Fn is far from the Fibonacci numbers Fn − 1 and Fn + 1. So the closest Fibonacci number to an approximation of Fn is exactly Fn.

Now if f is approximately Fn, then 5f2 is approximately a square. Find the integer r minimizing  |5f2 − r²|. In Python you could do this with the isqrt function. Then either r² + 4 or r² − 4 is divisible by 5. You can know which one by looking at r² mod 5. Whichever one is, divide it by 5 and you have the square of Fn. You can take the square root exactly, in integer arithmetic, and you have Fn. Furthermore, the number r that you computed along the way is the certificate of the calculation of Fn.

Visualizing orbital velocity

Posted on by John

The shape of a planet’s orbit around a star is an ellipse. To put it another way, a plot of the position of a planet’s orbit over time forms an ellipse. What about the velocity? Is its plot also an ellipse? Surprisingly, a plot of the velocity forms a circle even if a plot of the position is an ellipse.

If an object is in a circular orbit, it’s velocity vector traces out a circle too, with the same center. If the object is in an elliptical orbit, it’s velocity vector still traces out a circle, but one with a different center. When the orbit is eccentric, the hodograph, the figure traced out by the velocity vector, is also eccentric, though the two uses of the word “eccentric” are slightly different.

The eccentricity e of an ellipse is the ratio c/a where c is the distance between the two foci and a is the semi-major axis. For a circle, c = 0 and so e = 0. The more elongated an ellipse is, the further apart the foci are relative to the axes and so the greater the eccentricity.

The plot of the orbit is eccentric in the sense that the two foci are distinct because the shape is an ellipse. The hodograph is eccentric in the sense that the circle is not centered at the origin.

The two kinds of eccentricity are related: the displacement of the center of the hodograph from the origin is proportional to the eccentricity of the ellipse.

Imagine the the orbit of a planet with its major axis along the x-axis and the coordinate of its star positive.  The hodograph is a circle shifted up by an amount proportional to the eccentricity of the orbit e. The top of the circle corresponds to perihelion, the point closest to the star, and the bottom corresponds to aphelion, the point furthest from the star. For more details, see the post Max and min orbital speed.

Race between primes of the forms 4k + 1 and 4k + 3

Posted on by John

The last few posts have looked at expressing an odd prime p as a sum of two squares. This is possible if and only if p is of the form 4k + 1. I illustrated an algorithm for finding the squares with p = 2255 − 19, a prime that is used in cryptography. It is being used in bringing this page to you if the TLS connection between my server and your browser is uses Curve25519 or Ed25519.

World records

I thought about illustrating the algorithm with a larger prime too, such as a world record. But then I realized all the latest record primes have been of the form 4k + 3 and so cannot be written as a sum of squares. Why is p mod 4 equal to 3 for all the records? Are more primes congruent to 3 than to 1 mod 4? The answer to that question is subtle; more on that shortly.

More record primes are congruent to 3 mod 4 because Mersenne primes are easier to find, and that’s because there’s an algorithm, the Lucas-Lehmer test, that can test whether a Mersenne number is prime more efficiently than testing general numbers. Lucas developed his test in 1878 and Lehmer refined it in 1930.

Since the time Lucas first developed his test, the largest known prime has always been a Mersenne prime, with exceptions in 1951 and in 1989.

Chebyshev bias

So, are more primes congruent to 3 mod 4 than are congruent to 1 mod 4?

Define the function f(n) to be the ratio of the number of primes in each residue class.

f(n) = (# primes p < n with p = 3 mod 4) / (# primes p < n with p = 1 mod 4)

As n goes to infinity, the function f(n) converges to 1. So in that sense the number of primes in each category are equal.

If we look at the difference rather than the ratio we get a more subtle story. Define the lead function to be how much the count of primes equal to 3 mod 4 leads the number of primes equal to 1 mod 4.

g(n) = (# primes p < n with p = 3 mod 4) − (# primes p < n with p = 1 mod 4)

For any n, f(n) > 1 if and only if g(n) > 0. However, as n goes to infinity the function g(n) does not converge. It oscillates between positive and negative infinitely often. But g(n) is positive for long stretches. This phenomena is known as Chebyshev bias.

Visualizing the lead function

We can calculate the lead function at primes with the following code.

from numpy import zeros
from sympy import primepi, primerange

N = 1_000_000
leads = zeros(primepi(N) + 1)
for index, prime in enumerate(primerange(2, N), start=1):
    leads[index] = leads[index - 1] + prime % 4 - 2

Here is a list of the primes at which the lead function is zero, i.e. when it changes sign.

[   0,     1,     3,     7,    13,    89,  2943,  2945,  2947,
 2949,  2951,  2953, 50371, 50375, 50377, 50379, 50381, 50393,
50413, 50423, 50425, 50427, 50429, 50431, 50433, 50435, 50437,
50439, 50445, 50449, 50451, 50503, 50507, 50515, 50517, 50821,
50843, 50853, 50855, 50857, 50859, 50861, 50865, 50893, 50899,
50901, 50903, 50905, 50907, 50909, 50911, 50913, 50915, 50917,
50919, 50921, 50927, 50929, 51119, 51121, 51123, 51127, 51151,
51155, 51157, 51159, 51161, 51163, 51177, 51185, 51187, 51189,
51195, 51227, 51261, 51263, 51285, 51287, 51289, 51291, 51293,
51297, 51299, 51319, 51321, 51389, 51391, 51395, 51397, 51505,
51535, 51537, 51543, 51547, 51551, 51553, 51557, 51559, 51567,
51573, 51575, 51577, 51595, 51599, 51607, 51609, 51611, 51615,
51617, 51619, 51621, 51623, 51627]

This is OEIS sequence A038691.

Because the lead function changes more often in some regions than others, it’s best to plot the function over multiple ranges.

The lead function is more often positive than negative. And yet it is zero infinitely often. So while the count of primes with remainder 3 mod 4 is usually ahead, the counts equal out infinitely often.

Wagon’s algorithm in Python

Posted on by John

The last three posts have been about Stan Wagon’s algorithm for finding x and y satisfying

x² + y² = p

where p is an odd prime.

The first post in the series gives Gauss’ formula for a solution, but shows why it is impractical for large p. The bottom of this post introduces Wagon’s algorithm and says that it requires two things: finding a quadratic non-residue mod p and a variation on the Euclidean algorithm.

The next post shows how to find a quadratic non-residue.

The reason Wagon requires a non-residue is because he needs to find a square root of −1 mod p. The previous post showed how that’s done.

In this post we will complete Wagon’s algorithm by writing the modified version of the euclidean algorithm.

Suppose p is an odd prime, and we’ve found x such that x² = −1 mod p as in the previous posts. The last step in Wagon’s algorithm is to apply the Euclidean algorithm to x and p and stop when the numbers are both less than √p.

When we’re working with large integers, how do we find square roots? Maybe p and even √p are too big to represent as a floating point number, so we can’t just apply the sqrt function. Maybe p is less than the largest floating point number (around 10308) but the sqrt function doesn’t have enough precision. Floats only have 53 bits of precision, so an integer larger than 253 cannot necessarily be represented entirely accurately.

The solution is to use the isqrt function, introduced in Python 3.8. It returns the largest integer less than the square root of its argument.

Now we have everything necessary to finish implementing Wagon’s algorithm.

from sympy import legendre_symbol, nextprime
from math import isqrt

def find_nonresidue(p):
    q = 2
    while legendre_symbol(q, p) == 1:
        q = nextprime(q)
    return q

def my_euclidean_algorithm(a, b, stop):
    while a > stop:
        a, b = b, a % b
    return (a, b)

def find_ab(p):
    assert(p % 4 == 1)
    k = p // 4
    c = find_nonresidue(p)
    x = pow(c, k, p)
    return my_euclidean_algorithm(p, x, isqrt(p))

Let’s use this to find a and b such that x² + y² = p where p = 2255 − 19.

>>> a, b = find_ab(p := 2**255 - 19)
>>> a
230614434303103947632580767254119327050
>>> b
68651491678749784955913861047835464643
>>> a**2 + b**2 - p
0

Finis.

Finding a square root of -1 mod p

Posted on by John

If p is an odd prime, there is a theorem that says

x² = −1 mod p

has a solution if and only if p = 1 mod 4. When a solution x exists, how do you find it?

The previous two posts have discussed Stan Wagon’s algorithm for expressing an odd prime p as a sum of two squares. This is possible if and only if p = 1 mod 4, the same condition on p for −1 to have a square root.

In the previous post we discussed how to find c such that c does not have a square root mod p. This is most of the work for finding a square root of −1. Once you have c, set

x = ck

where p = 4k + 1.

For example, let’s find a square root of −1 mod p where p = 2255 − 19. We found in the previous post that c = 2 is a non-residue for this value of p.

>>> p = 2**255 - 19
>>> k = p // 4
>>> x = pow(2, k, p)

Let’s view x and verify that it is a solution.

>>> x
19681161376707505956807079304988542015446066515923890162744021073123829784752
>>> (x**2 + 1) % p
0

Sometimes you’ll see a square root of −1 mod p written as i. This makes sense in context, but it’s a little jarring at first since here i is an integer, not a complex number.

The next post completes this series, giving a full implementation of Wagon’s algorithm.