Jeff Epler's blog

Calculating unreasonably large Fibonacci numbers quickly

The Fibonacci sequence is an example of a recursively defined sequence of numbers. The first two numbers in the sequence are defined as 0 and 1, and each following number is defined as the sum of the two previous numbers.

Given this definition, it is very easy to write an extremely inefficient recursive function to calculate the n'th number:

def fib(n):
    if n == 0 or n == 1: return n
    return fib(n-1) + fib(n-2)

In most standard languages (the above is intended to be Python) this is inefficient because the program is approached as a sequence of operations. Each time "fib(13) is encountered, the body of the fib function is executed. This in turn evaluates "fib(12)" and "fib(11)". "fib(12)" evaluates "fib(11)" again. In the end, the total number of times a "fib(n)" is evaluated grows very quickly compared to "n". In fact, I think it is one of those cases of increasing literally exponentially. This is not a good method to find the Fibonacci sequence or individual Fibonacci numbers.

Another way that works better is a loop that calculates each number in the sequence from the previous two, until the desired one is reached:

def fib(n):
    a, b = 0, 1
    for i in range(n):
        a, b = a+b, a
    return a

Now, the n'th number is calculated after only n additions. A big improvement! This method is quite good for finding the entire Fibonacci sequence up to some 'n', or for finding particular values for 'n' in the hundreds or thousands.

Past this, if you don't need an exact number, you can use floating-point arithmetic, using a closed-form expression involving (this sounds funny if you don't know it's math words, don't laugh) "powers of the golden ratio".

import decimal

def fib(n):
    sqrt_5 = decimal.Decimal(5).sqrt()
    phi = (1 + sqrt_5) / 2
    return (phi**n - (-phi)**(-n)) / sqrt_5

.. fib(4784970) is Decimal('1.735729075920066439815208698E+999999'), give or take some of the last decimal places. After that, Python's standard decimal package throws in the towel.

You can make approximations beyond this by doing log-math. First, notice how the (-phi) ** (-n) term goes quickly to zero, leaving just phi ** n / sqrt(5). You can find the log10 of a really huge Fibonacci numbe as:

from math import log

def logfib(n, base=10):
    log_sqrt_5 = log(5, base) / 2
    log_phi = log((1 + 5**.5) / 2, base)
    return log_phi * n - log_sqrt_5

def print_log(n):
    whole, frac = divmod(n, 1)
    print(f"{10 ** frac}e{whole:.0f}")

Here, we can get 1.735729075774171e999999 as an estimate of fib(4784970), but can also easily get 1.4140750583014141e417975280 as an estimate of fib(2_000_000_000). As you'll see shortly, the first 8 digits of that second value are accurate.

The method I really want to talk about, though, is the one which can exactly calculate Fibonacci numbers—all the binary or decimal digits, and can actually find the 2 billionth number in short order and the 16 billionth one with effort (and a bigger swapfile). The program usually skips printing the actual number since converting the number from binary to decimal and printing it actually takes a fair amount of time.

$ time fib4.py 2000000000
2000000000: 417975281 digits: 141407506869875260956516…969828793289064439453125
217.43user 14.56system 0:40.17elapsed 577%CPU…
$ time ./fib4.py 16000000000
16000000000: 3343802244 digits: 446863351203393654391107…465863000407211560546875
1931.98user 205.95system 8:22.67elapsed 425%CPU…

How's this work? As usual, by standing on the shoulders of giants. First, the matrix form of the Fibonacci number, in which case the n'th number can be found by taking a particular element of the n'th power of a 2x2 matrix. Second, the "repeated squaring" algorithm, for finding the n'th power in only about log2(n) operations. Third, the gmp library for high quality multiprecision library and its "gmpy" wrapper for Python.

Besides the script below, I'm using a modified version of gmpy which allows the "+" and "*" operations to be multithreaded. This can give a good speed-up as the matrix multiplication step's individual element multiplcations can take place in parallel. As shown above, this improves the elapsed time by more than a factor of 4.

The main barrier to going further seems to be RAM. Doubling n increases the memory needs by about 1.88x, and for the 2 billionth Fibonacci number we're only at about 418 million digits. While the records are strictly incomparable, this is on a par with the largest number of digits of pi computed…in 1999.

I don't know that anyone's worked to set a record largest Fibonacci number calculated, but let this program be my entry in the contest if so.

[permalink]

Heat-staking in 3D prints

This is not a technique I've perfected, but I think it has some interesting applications and I haven't seen anyone else write about it.

I frequently design 3D prints that include PCBs. Common techniques for retaining PCBs in 3D prints include screws and clips. Clips can interfere with access to PCB pins/headers, while screws need, well, screws and sometimes nuts or threaded inserts.

While contemplating this problem, I was reminded of the construction of common tactile switches. A small metal plate is retained by little plasic blobs at each corner. Presumably, the blobs were originally part of the main molded plastic switch body, and the manufacturing process first placed the metal plate and then melted the plastic. This process seems to be called "heat staking" or "thermal staking".

The main downside of heat staking is that you can't remove the PCB without destroying the print. Ouch.

There are various staking processes. This source shows some possible geometries. Unfortunately, I don't have the special concave tools required. However, knowing suggested geometries is helpful! For instance, given a hole of diameter 'd', a boss that is initially 1.6d above the surface to bond should flatten to a diameter of 2d and a height of .5d, at least using a flared stake tool with the right geometry.

I don't have a special tool, but I do have some nice broad flat tips for my iron: these heat set insert tips. Unfortunately, the inexpensive iron I use with them can only be adjusted down to 280C and I am pretty sure that's just a power setting, not a temperature setting. Melting PETG even at the "280C" setting seems to release some pretty nasty fumes, I don't recommend it.

There's also the problem of the plastic sticking/melting onto the tip. To mitigate this, I use a piece of polyimide (often called by the brand name Kapton) tape between the iron and the plastic. The same piece of tape can be moved & re-used multiple times.

My procedure:

1. Design the bosses into the 3D printed part
2. Print. My test pieces have been in PETG plastic.
3. Place the PCB on the bosses
4. Put polyimide tape over a boss
5. Heat my iron for a bit, turn it off, then press the tip down onto the tape to melt the boss
6. Repeat the last two steps for each stake

I don't know how much further I'll pursue this technique, but it's a good one to have in the toolbox!

[permalink]

wwvbtk: Visualize wwvb timecodes with Tkinter

I added a real-time visualization program to wwvbpy.

[permalink]

Wrap text according to font width in CircuitPython

I wrote a small function called wrap_text which re-wraps a block of text so that, when rendered in a given font, it is no more than a given number of pixels wide. (There are other functions available for wrapping text according to the number of letters but when the letter "l" is much narrower than the letter "w", this can give unsatisfactory results on specific texts)

The program splits the input string on any whitespace, then uses a simple "greedy" algorithm to add each word to the current line if it fits, otherwise break and continue with the next line.

The indentation of the first and subsequent lines can be specified; the default indents the 2nd line because my use case was to wrap text of a poem.

A new string is returned, with newlines inserted in the required spots, suitable for use with a standard circuitpython label.

The sample program is for the Adafruit MagTag. It uses this function to display a poem (placed in CIRCUITPY/poem.txt) 5 lines at a time, updating once per hour. The first line in poem.txt is always displayed at the bottom, and is intended to hold the poem's name and author. By using deep sleep mode, I expect that I'll only need to charge the device every few days.

I wrote this because I was very moved by Amanda Gorman's poem "The Hill We Climb". I could tell this was a poem that would really "stick with me" and wanted the opportunity to experience the poem over a longer period of time. I can look at the MagTag several times a day and reflect on whatever lines happen to be visible at the moment.

Out of respect for Ms. Gorman's copyright I don't include the poem's full text here. However, many news outlets have quoted it in its entirety. And, of course, you can put whatever you like in poem.txt—Why not something by Maya Angelou or Lin-Manuel Miranda, apparently two people who Ms. Gorman has found inspiring.

"Festive" 3D printed stand on thingiverse

Files currently attached to this page:

[permalink]

Seven Years of Prius (2013) Fuel Efficiency

At the tail end of 2019 I changed to a work from home job, and of course COVID ended most activities outside the house for the bulk of 2020, leading to such low mileage. Fuel economy went up a bit.

[permalink]

Where should CircuitPython go in 2021?

Looking back at 2020

In my list for last year, I focused on timekeeping and time measurement. We did do two out of the four goals I mentioned: the built-in time module now works from a real-time clock peripheral, and I implemented minimal time-zone handling for CircuitPython, though it's not sufficiently polished to add to the community bundle.

Together, these enabled one of my 2020 projects, a CircuitPython clock that is accurate to within better than 1 second per month, handles my local time zone rule, and has a battery back up so it doesn't need to be set after power is lost. I never blogged the full project, but there are bits and pieces scattered about: 1 2 3 4

On the other hand, the things I'm proudest of helping with in 2020 didn't even figure on my list—they were unanticipated. The two items I'm thinking of are bringing Zoltán Vörös's ulab into CircuitPython, and making it much easier for users to contribute to CircuitPython translations by setting up Hosted Weblate.

Looking forward to 2021

This year, I don't want to be so closely focused on features that just a few people may be interested in. Also, many of the 2020 submissions I was most struck by were those that focused more on the human side than on the technical. With that in mind, here are some of the broader goals I think we should set:

• Support at least one new microcontroller family
• Continue to help people grow into the roles of reviewer and contributor
• Spotlight or create more "how to" video content around CircuitPython, especially for the more complex libraries

There are also some pain points that we need to deal with sooner or later, but in many cases the right idea hasn't come up yet. We should keep an eye out for ideas on solving these problems:

• Increasing time & complexity to build CircuitPython core
• Risk of being spread too thin by multiple in-progress ports
• No automatic way to gather items from the Bundle for a particular application
• Wireless-enabled devices deserve a pure Wi-Fi or Bluetooth development workflow
• There's an increasing need for a path manipulation library
• We should standardize the location of common files like fonts, icons, and sounds
• Read/write storage that doesn't interfere with access to CIRCUITPY from a host computer

Finally, some personal goals I have related to CircuitPython and projects in general:

• Add an antialiased font format to adafruit_bitmap_font
• Design and fabricate PCBs for my projects
• Improve my skills at 3D printed enclosure design, including learning FreeCAD
• Share the PCB and 3D printable designs

Spilling the beans on one particular project I'd love to do in 2021, why isn't there a MagTag-like podcast player appliance running CircuitPython? There are a few missing pieces before this becomes possible, but I think we could—and should—do it.

I fully anticipate that once again the best things in 2021 will be the ones I didn't even imagine. I can't wait to see what this community builds together in 2021.

[permalink]

ForkAwesome font converted for CircuitPython

The ttf file is also included in case you need to generate other sizes.

The font uses character codes starting at 0xf000, so using regular characters like "012" or "abc" won't show anything. Instead, to find the code of a character, visit the forkawesome website, click the icon you want, and look for the Unicode (f322 for the fa-python icon). Then, in your CircuitPython code, create a label with the text "\uf322" or chr(0xf322) and the FontAwesome font to show it:

from adafruit_bitmap_font.bitmap_font import load_font
from adafruit_display_text.label import Label
from displayio import Group

font = load_font("/forkawesome-36.pcf")

text = chr(0xf322)

group = Group()
label = Label(font=font, text=text, background_color=0x111111)
label.anchor_point = (0, 0)
label.anchored_position = (0,0)

group.append(label)
board.DISPLAY.show(group)
while True:
    pass

You can also use the attachment icons.py as a reference for the character codes. The python identifiers are the same as the font's "id", except that "-" is replaced with "_","500px" is renamed to "fivehundredpx", and several names that are Python built in functions or reserved words have an underscore appended (file_, filter_, list_, map_, print_, try_).

Files currently attached to this page:

[permalink]

Head Scratching C(++) Bug

/* n is a parameter of type size_t */
/* table is a pointer to int such that table[0] .. table[n-1] are all valid indices */
for(size_t i = 0; i < n; i++) {
    if(table[i] == needle) return i;
}

Instead of terminating after at most n iterations, the loop continued indefinitely, eventually returning some nonsense value where another memory location matched the needle value.

Give up? The problem's nothing to do with the loop, but with what came (didn't come) after it:

size_t find_in_table(int *table, size_t n, int needle) {
    for(size_t i = 0; i < n; i++) {
        if(table[i] == needle) return i;
    }
}

There's no return statement after the loop. If the loop terminates, it reaches the end of a non-void function without returning a value. This is undefined behavior. Therefore, the compiler-author logic goes, it must be the case that the loop can never terminate due to the loop condition; it must only terminate by the return statement inside the loop. (Stated in another way, the "undefined behavior" is to continue to run the loop.)

Without the final return (e.g., return (size_t)-1;), the comparison i < n is completely eliminated at typical optimization levels. Even weirder, you can get the program to print that (say) i==5, n==4, and that i<n is true!

I think gcc has an enabled-by-default warning for this, but either it's disabled by my project's CFLAGS or I overlooked it in a torrent of other compiler output.

Without looking beyond the loop or if you're thinking with the mindset that a missing return statement's undefined behavior is "returns an arbitrary value", you don't think about things like this. This is a very easy mistake to make, since it is also a pretty good model of what a modern compiler will do at "for better debugging" optimization levels, as well as what the "C is portable assembler" moto asks us to falsely believe. No, in reality, C(++) is a slippery language. While it's not a bet you could settle in the negative, I bet there's at least one security bug out there due to this undefined behavior and the way that gcc transforms a finite loop into an infinite one.

My recommendation?

Pay attention to compiler diagnostics. Enable more of them. Understand what the diagnostic is saying, then make an appropriate response. Get your editor's help to parse error lists and make sure you don't miss any.

[permalink]

wwvbgen updates

Pi Zero W USB Proxy

GPS and relativity

Some notes on the Si5351a

Si5351 Frequency Planner in Python

Quad CharliePlex FeatherWing hack

Minimal Time-Zone Handling for CircuitPython

Calibrating the DS3231 and PCF8523 RTCs