The On-Line Encyclopedia of Integer Sequences gets regularly updated with new integer sequences. One of the recent updates was contributed by me, A285494.

A285494 is the list of all numbers k so that its digit sum equals its number of distinct prime factors.
A number’s digit sum is the sum of all of its decimal digits. The number 62831853, for example, has a digit sum of 6+2+8+3+1+8+5+3 = 36.
A number’s number of distinct prime factors is the number of different prime numbers that multiply together to result in the original number. As an example, 62831853 = 3^2 \cdot 7 \cdot 127 \cdot 7853, so it has five prime factors of which four are distinct.
Thereby one can conclude that 62831853 is not an entry in this sequence, as 36 \neq 4.

The sequence is certainly infinite, as the number k = 2 \cdot 10^n with n \in \mathbb{N}^* has a digit sum of 2 + (0 \cdot n) = 2 and — because k = 2^{n+1} \cdot 5^n — exactly two distinct prime factors.

In the encyclopedia entry, I provided a Mathematica one-liner to compute the first few entries of this sequence. Since then, I have also written a Python two-liner to achieve the same goal.

(* Mathematica *)
Out = {20, 30, 102, 120, 200, 300, 1002, 1200, 2000, 2001, 2002, 3000, 3010}
# Python 2.7
>>> def p(n):exec"i,f=2,set()\nwhile n>1:\n\tif n%i==0:f.add(i);n/=i;i=1\n\ti+=1";return len(f)
>>> print filter(lambda n:p(n)==sum(map(int,str(n))),range(2,10001))
[20, 30, 102, 120, 200, 300, 1002, 1200, 2000, 2001, 2002, 3000, 3010]

Multibrot Set

The Mandelbrot Set is typically defined as the set of all numbers c \in \mathbb{C} for which — with z_0 = 0, z_{n+1} = f_c(z_n) and f_c(z) = z^2 + c — the limit \lim\limits_{n \to \infty} z_n converges. Visualizations of this standard Mandelbrot Set can be seen in three of my posts (Mandelbrot Set, Mandelbrot Set Miscalculations and Mandelbrot Set II).

f_c(z)=z^2+cHowever, one can extend the fractal’s definition beyond only having the exponent 2 in the function to be f_c(z)=z^\text{exp}+c with \text{exp} \in \mathbb{R}. The third post I mentioned actually has some generalization as it allows for \text{exp} \in \{2,3,4,5\}, although the approach used cannot be extended to real or even rational numbers.

f_c(z)=z^3+cThe method I used in the aforementioned post consists of manually expanding (a+b\cdot i)^n for each n. The polynomial (a+b\cdot i)^3, for example, would be expanded to (a^3 - 3 \cdot a \cdot b^2) + (3 \cdot a^2 \cdot b - b^3) \cdot i.
This method is not only tedious, error-prone and has to be done for every exponent (of which there are many), it also only works for whole-number exponents. To visualize real Multibrots, I had to come up with an algorithm for complex number exponentiation.

f_c(z)=z^4+cLuckily enough, there are two main ways to represent a complex number, Cartesian form z = a+b\cdot i and polar form z = k\cdot e^{\alpha\cdot i}. Converting from Cartesian to polar form is simply done by finding the number’s vector’s magnitude k = \sqrt{a^2+b^2} and its angle to the x-axis \alpha = \mbox{atan2}(\frac{a}{b}). (The function \mbox{atan2} is used in favor of \arctan to avoid having to divide by zero. View this Wikipedia article for more on the function and its definition.)
Once having converted the number to polar form, exponentiation becomes easy as z^\text{exp} = (k \cdot e^{\alpha\cdot i})^\text{exp} = k^\text{exp} \cdot e^{\alpha \cdot \text{exp} \cdot i}. With the exponentiated z^\text{exp} in polar form, it can be converted back in Cartesian form with z^\text{exp} = k^\text{exp} \cdot (\cos{(\alpha \cdot \text{exp})} + \sin{(\alpha \cdot \text{exp})} \cdot i \big).

f_c(z)=z^5+cUsing this method, converting the complex number to perform exponentiation, I wrote a Java program which visualizes the Multibrot for a given range of exponents and a number of frames.
Additionally, I added a new strategy for coloring the Multibrot Set, which consists of choosing a few anchor colors and then linearly interpolating the red, green and blue values. The resulting images have a reproducible (in contrast to randomly choosing colors) and more interesting (in contrast to only varying brightness) look.

f_c(z)=z^6+cThe family of Multibrot Sets can also be visualized as an animation, showing the fractal with an increasing exponent. The animated gif shown below was created using ImageMagick’s convert -delay <ms> *.png multibrot.gif command to stitch together the various .png files the Java application creates. To speed up the rendering, a separate thread is created for each frame, often resulting in 100% CPU-usage. (Be aware of this should you render your own Multibrot Sets!)

f_c(z)=z^10+cTo use the program on your own, either copy the source code listed below or download the .java file. The sections to change parameters or the color palette are clearly highlighted using block comments (simply search for ‘/*’).
To compile and execute the Java application, run (on Linux or MacOS) the command javac multibrot.java; java -Xmx4096m multibrot in the source code’s directory (-Xmx4096m tag optional, though for many frames at high quality it may be necessary as it allows Java to use more memory).
If you are a sole Windows user, I recommend installing the Windows 10 Bash Shell.

Multibrot animation (probably loading...)

// Java 1.8 Code
// Jonathan Frech, 11th of September 2016
//          edited 17th of April     2017
//          edited 18th of April     2017
//          edited 20th of April     2017
//          edited 21st of April     2017
//          edited 22nd of April     2017

Continue reading


The On-Line Encyclopedia of Integer Sequences (also known by its acronym, OEIS) is a database hosting hundreds of thousands of — as the name implies — integer sequences. Yet, despite the massive number of entries, I contributed a new integer sequence, A278328.

A278328 describes numbers whose absolute difference to their decimal reverse are square. An example would be 12 or 21 (both are the decimal reverse to each other), since \left|12-21\right|=9 and 9=3^2.

Not a whole lot is known about the sequence, partly due to its definition only resulting in the sequence when using the decimal system, though it is known that there are infinitely many numbers with said property. Since there are infinitely many palindromes (numbers whose reverse is the number itself), \left|n-n\right|=0 and 0=0^2.

Due to there — to my knowledge — not being a direct formula for those numbers, I wrote a Python script to generate them. On the sequence’s page, I posted a program which endlessly spews them out, though I later wrote a Python two-liner, which only calculates those members of the sequence in the range from 0 to 98 (shown below entered in a Python shell).

>>> import math
>>> filter(lambda n:math.sqrt(abs(n-int(str(n)[::-1])))%1 == 0, range(99))
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 15, 21, 22, 23, 26, 32, 33, 34, 37, 40, 43, 44, 45, 48, 51, 54, 55, 56, 59, 62, 65, 66, 67, 73, 76, 77, 78, 84, 87, 88, 89, 90, 95, 98]

Collatz Conjecture

The Collatz conjecture states that every positive integer k>0 will — if you iteratively set k to f(k) — result in 1 (function shown beneath).
The graph beneath shows the path length of numbers from 1 to 10\,000. In this range 6170 is the number with the most steps, 261.
f(k)={\begin{cases}\frac{k}{2}&{\text{if }}k \mod 2 = 0\\3 \cdot k+1&{\text{if }}k \mod 2 = 1\end{cases}}

Path length in Collatz sequence [1 to 10**4]

# Python 2.7.7 Code
# Pygame 1.9.1 (for Python 2.7.7)
# Jonathan Frech 2nd of September, 2016

Continue reading


Interpreting the hour hand on a clock as a two-dimensional object on a plane, the hand’s tip can be seen as a complex number.
This clock converts the hour hand’s position into a complex number, sets the number’s length to the current minutes and displays it in the form a + b \cdot i.
The angle \phi is determined by the hours passed (\frac{2 \cdot \pi \cdot \text{hour}}{12} = \frac{\pi \cdot \text{hour}}{6}) but has to be slightly modified because a complex number starts at the horizontal axis and turns anti-clockwise whilst an hour hand starts at the vertical axis and turns — as the name implies — clockwise.
Thus \phi = (2 \cdot \pi - \frac{\pi \cdot \text{hour}}{6}) + \frac{\pi}{2} = (\frac{15 - \text{hour}}{6}) \cdot \pi.
The complex number’s length is simply determined by the minutes passed. Because the length must not be equal to 0, I simply add 1. |z| = k = \text{minute} + 1.
Lastly, to convert a complex number in the form k \cdot e^{\phi \cdot i} into the form a + b \cdot i, I use the formula k \cdot (\cos{\phi} + \sin{\phi} \cdot i) = a + b \cdot i.


# Python 2.7.7 Code
# Pygame 1.9.1 (for Python 2.7.7)
# Jonathan Frech 29th of July, 2016

Continue reading

Triangular Squares

In a recent video Matt Parker showed a triangular number that also is a square number, 6, and asked if there were more.

A triangular number has the form \frac{n^2+n}{2} — shown by Euler — and a square number has the form m^2.
Triangular squares are those numbers for which \frac{n^2+n}{2} = m^2 with n,m \in \mathbb{N}.
Examples are \{0, 1, 6, 35, 204, 1189, 6930, \dots\} (sequence A001109 in OEIS).

To check if triangular numbers are square numbers is easy (code listed below), but a mathematical function would be nicer.
The first thing I tried was to define the triangular number’s square root as a whole number, \sqrt{\frac{n^2+n}{2}} = \lfloor \sqrt{\frac{n^2+n}{2}} \rfloor. This function does not return the square numbers that are triangular but the triangular numbers that are square.
The resulting sequence is \{0, 1, 8, 49, 288, 1681, 9800, \dots\} (sequence A001108 in OEIS).

# Python 2.7.7 Code
# Jonathan Frech 13th of July, 2016
#         edited 15th of July, 2016

Continue reading

Palindrome Function

To get a number’s palindrome in a programming language like python is easy. There are ways to swap between integer and string and strings can be manipulated.

>>> n = 1234
>>> int(str(n)[::-1])

But I wanted to create a mathematical function p(n), which returns an integer’s palindrome. Thus p(1234) = 4321.

Firstly I needed a way of determining the number’s size. In base 10 the length is calculated using the logarithm to said base.
l(n) = \lfloor \log_{10}{n} \rfloor + 1
l(1234) = \lfloor \log_{10}{1234} \rfloor = \lfloor3.09 \rfloor + 1 = 4

Secondly I need a way to isolate a specific digit. Using the floor function, this function returns the i\text{-th} digit (starting on the right with i=0).
d_i(n) = \lfloor \frac{n}{10^i} \rfloor - \lfloor \frac{n}{10^{i+1}} \rfloor \cdot 10
d_2(1234) = \lfloor \frac{1234}{10^2} \rfloor - \lfloor \frac{1234}{10^{2+1}} \rfloor \cdot 10 = \lfloor 12.34 \rfloor - \lfloor 1.23 \rfloor \cdot 10 = 12 - 1 \cdot 10 = 2

Thirdly both of these functions can be used to split up the number into a sum.
n = \sum\limits_{i=0}^{l(n)-1} \Big[ d_i(n) \cdot 10^{i} \Big] = \sum\limits_{i=0}^{\lfloor \log_{10}{n} \rfloor} \Big[ \big( \lfloor \frac{n}{10^i} \rfloor - \lfloor \frac{n}{10^{i+1}} \rfloor \cdot 10 \big) \cdot 10^{i} \Big]

Fourthly I only need to swap the power of ten at the end to get my palindrome function.
p(n) = \sum\limits_{i=0}^{l(n)-1} \Big[ d_i(n) \cdot 10^{l(n) - 1 - i} \Big] = \sum\limits_{i=0}^{\lfloor \log_{10}{n} \rfloor} \Big[ \big( \lfloor \frac{n}{10^i} \rfloor - \lfloor \frac{n}{10^{i+1}} \rfloor \cdot 10 \big) \cdot 10^{\lfloor \log_{10}{n} \rfloor - i} \Big]

Thus the final function p(n) is defined.
p(n) = \sum\limits_{i=0}^{\lfloor \log_{10}{n} \rfloor} \Big[ \big( \lfloor \frac{n}{10^i} \rfloor - \lfloor \frac{n}{10^{i+1}} \rfloor \cdot 10 \big) \cdot 10^{\lfloor \log_{10}{n} \rfloor - i} \Big]

To check if the formula is correct, I use 1234 (as seen above).
p(1234) = \sum\limits_{i=0}^{\lfloor \log_{10}{1234} \rfloor} \Big[ \big( \lfloor \frac{1234}{10^i} \rfloor - \lfloor \frac{1234}{10^{i+1}} \rfloor \cdot 10 \big) \cdot 10^{\lfloor \log_{10}{1234} \rfloor - i} \Big]
p(1234) = \sum\limits_{i=0}^{3} \Big[ \big( \lfloor \frac{1234}{10^i} \rfloor - \lfloor \frac{1234}{10^{i+1}} \rfloor \cdot 10 \big) \cdot 10^{3 - i} \Big]
p(1234) = d_0(1234) \cdot 10^3 + d_1(1234) \cdot 10^2 + d_2(1234) \cdot 10^1 + d_3(1234) \cdot 10^0
p(1234) = 4000 + 300 + 20 + 1 = 4321


A cycloid is the curve generated by tracing a point on a smaller circle rolling around an bigger circle.
The word cycloid comes from the greek word κύκλος meaning circle. There are epicycloids and hypocycloids (smaller circle located above and beneath the bigger circle).
The cycloid is determined by the ratio between the two circle’s radii k = \frac{R}{r}.
More information can be found on this Wikipedia entry.

Hypocycloid with k = 2.1 Hypocycloid with k = 2.1 Hypocycloid with k = 7.2

# Python 2.7.7 Code
# Pygame 1.9.1 (for Python 2.7.7)
# Jonathan Frech 3rd of June, 2016
#         edited 4th of June, 2016

Continue reading