Fun With Binary Palindromes

My recent vicissitudes have led me into the realm of palindromes. Binary palindromes.

Essentially, palindromes are words or numbers that appear the same when read from left or right. The word “Anna” is a palindrome for example, as well as the number 13931.

But my focus has been primarily on binary palindromic numbers, such as 11011 and 01011010. There is nothing inherently useful about them, though they possess a certain beauty which is only evident to someone who likes to fumble with ones and zeroes.

With the help of a friend I conducted some mindless research concerning the Mathematical properties of these numbers and stumbled upon some things.

The first pattern

For a certain amount of bits b, there exists a set of palindromes P, and P_n \in \mathbb{Z}. The following table shows all the palindromes for b=5 in sorted order. Keep in mind that I will continue to use the 5 bit-system for examples hereon.

n P_n Decimal
0 00000 0
1 00100 4
2 01010 10
3 01110 14
4 10001 17
5 10101 21
6 11011 27
7 11111 31

The next thing we asked ourselves was if there lies a sequence behind the pattern of palindromes, or a way to express the nth palindrome?

It’s quite easy to notice a pattern: if you glance at the middle column of each number, you find that it repeats the pattern 0 – 1, thus switching state repeatedly. Similarly, the columns directly nearby repeat the pattern 0 – 0 – 1 – 1, and the columns next to those 0 – 0 – 0 – 0 – 1 – 1 – 1 – 1.

So, the farther one is from the middle column, the longer it takes for the bit to toggle. More specifically, the bit-toggle takes place after 2^i palindromes, whereby i is the offset from the middle. The same holds for an even amount of bits, where the middle is two bits.

With this knowledge it is easy to devise a program to calculate all palindromes for any amount of bits. The following example is Java code.

We need more math!
We now know the algorithm for calculating the set of palindromes for any amount of bits. This leaves us to the next step: a way of expressing the nth palindrome.

Through endless fumbling we realised that each palindrome is the sum of a set of distinctive palindromes in a particular combination. We will call these palindromes elementary palindromes.

Those palindromes have the property of having only two bits set, which are mirrored in the middle, e.g. 01010, 10001, etc. — they are the simplest palindromes in the entire set, and thus it is easy to conclude that every palindrome incorporates them in some way.

For odd bit-systems, they take the form of e_k = 2^{\left \lfloor{\frac{n}{2}}\right \rfloor- k } + 2^{\left \lfloor{\frac{n}{2}}\right \rfloor + k } where k is the offset from the middle. The catch is that the first elementary palindrome is only a single bit, such as 00100. Therefore e_0 = 2^{\left \lfloor{\frac{n}{2}}\right \rfloor}.

The formula is only slightly different for even bit-systems: e_k = 2^{\frac{n}{2} - k - 1} + 2^{\frac{n}{2} + k }.

Now, any palindrome can be written as P_n = a_0 e_0 + a_1 e_1 + a_2 e_2 + ... where a_k are coefficients that are either zero or one. Each palindrome has its own unique combination of these coefficients. But how can this help us compute the nth palindrome save by chance guessing?

Extending the prior table by the combination of coefficients:

n P_n Decimal a_2 a_1 a_0
0 00000 0 0 0 0
1 00100 4 0 0 1
2 01010 10 0 1 0
3 01110 14 0 1 1
4 10001 17 1 0 0
5 10101 21 1 0 1
6 11011 27 1 1 0
7 11111 31 1 1 1

What follows is an interesting pattern: if the coefficients are read like a binary number then that number is equal to n! This allows for feasible computation of the nth palindrome. The only thing required is to calculate the elementary palindromes and to assign to the coefficients the binary representation of n.

Let’s compute for a 5-bit-system. Now:

e_0 = 2^{\left \lfloor{\frac{5}{2}}\right \rfloor} = 2^{2} = 4.

e_1 = 2^{\left \lfloor{\frac{5}{2}}\right \rfloor- 1 } + 2^{\left \lfloor{\frac{5}{2}}\right \rfloor+ 1 } = 10.

e_2 = 2^{\left \lfloor{\frac{5}{2}}\right \rfloor- 2 } + 2^{\left \lfloor{\frac{5}{2}}\right \rfloor+ 2 } = 17.

Say we want to compute the fifth palindrome. 5 is 101 in binary. Thus P_5 = 1 \times e_0 + 0 \times e_1 + 1 \times e_2 = 1 \times 4 + 0 \times 10 + 1 \times 17 = 21. Glancing at the above table that result is indeed correct.

The Julia set

Quite like I described the Mandelbrot set, the Julia set is also a Mathematical fractal that is computed with almost the exact iteration. The iteration is

z_{n+1}=z_n^2+c,

a quadratic polynomial in this case. Here, z and c are both complex numbers, but play a rather different role than in the Mandelbrot equation.

Going back to our original understanding of this iteration, c was a constant mapped to a point in the set, and z a value used for computation with an initial value of z_0 = 0. Now, by repeating the iteration it was decided that if |z| > 4 the point was not part of the Mandelbrot set. As the points inside the set will naturally not converge (for they are part of the set), the iteration is restricted to a finite number of recursions.

What separates the Julia set(s) is that the meanings of z and c have been swapped. Now z has an initial value mapped to a point in the set and c is now an arbitrary constant. Thus, there exist an infinite amount of Julia sets.

I will present the Julia set only in its form of a quadratic polynomial, as we will discover semblances to the Mandelbrot set when discussing the many variations.

Some Julia sets

Let us start with the most simple form of the set when we set c = 0, thus reducing the iteration to z_{n+1} = z_n^2. What follows is a simple circle, also called the unit circle. It has a radius of one and sits at the origin.

screen-shot-2016-11-20-at-15-40-38
c = 0

 

A fine explanation for this figure is that any number n \in \mathbb{C} will always to converge to zero (and thus be part of the set) if |n| < 1. The entire set of this form are the points with this property.

When we increase the real value of c we see a fractal of spirals emerge which appear to be mirrored both horizontally and vertically.

screen-shot-2016-11-20-at-15-40-55
c = 0.3

As we increase that value the spirals begin to distance themselves from each other and shrink. Notable is that the vertical distance grows faster than the horizontal one.

For larger imaginary values of c we find another interesting spiral pattern. Hereby the spirals repeat inside each other infinitely. At c=0.76i the spirals begin to fade out and the figure looks more like branches connecting with each other. Similar to the shape before, the figure dissolves as the value becomes larger. If we now make the value negative we see the same pattern only mirrored horizontally.

c = 0.7i
c = 0.7i
c = -0.7i
c = -0.7i

When we assign negative real values to c we find some astonishing similarities to the Mandelbrot set. At c=-0.67 the figures becomes very much resembling of the Mandelbrot shape. The figure is again mirrored on both axes. Increasing that value in the negative direction condenses the fractal into lines, or branches. This looks similar to the tail of the Mandelbrot set.

screen-shot-2016-11-20-at-15-42-19
c=-0.67
c=-1.42
c=-1.42

More interesting shapes can be found if mixing real and imaginary values. I encourage you to explore by yourself!

But what about the Mandelbrot set? All of the figures presented can be indeed found inside the Mandelbrot set, as the Mandelbrot set is actually the sum of all possible Julia sets. Astounding, isn’t it? This is because the Mandelbrot set at each point has a different value of c. And that is why all Julia sets (in the form of a quadratic polynomial) combined form the Mandelbrot set.

I have vaguely hinted before, but the Julia set may have different iteration types. For example there is z_{n+1} = z_n^3+c or even z_{n+1} = \sin(z) \cdot c! The Internet will gift you with astonishing videos and images about these at no cost at all.

Computational Complexity of Recursive Functions

A function is deemed recursive if it satisfies the following three rules:

  1. There must be at least one halting condition.
  2. The function must call itself.
  3. The parameters thereof must be altered such that the halting condition will eventually stop the recursion.

The easiest example of a recursive function would be calculating the factorial of a (natural) number.

As there is only a single recursive call, the recursive function will call itself a total of n times. Therefore its complexity is O(n).

A more interesting example is an algorithm for calculating the nth Fibonacci number.

Now there is a total of two recursive calls with two halting conditions (also called the seed values of the iteration). Noticeable is the curious amount of calls it takes for one recursion to end. If we graph the amount of calls vs. n we find a graph looking similar to an exponential function.

screen-shot-2016-11-06-at-20-57-01

As we now assume that the function determining the amount of calls is exponential, and that the amount of calls of n is equal to the amount of calls of n-1 plus the amount of calls of n-2. We can now, by induction, find the function mathematically.

Let  F(n) = F(n-1) + F(n-2).

Since we assume the function to be exponential we can substitute: x^n = x^{n-1} + x^{n-2}.

Dividing by x^{n-2} leaves us with x^2= x + 1.

Solving, we find x = \frac{1+\sqrt{5}}{2} \approx 1.618, also called the Golden ratio. The second solution can be discarded as it is negative (approaches zero as n approaches infinity). Thus, the computational complexity of the Fibonacci function happens to be O(\phi^n). But as there is bias to powers of two in computer science the complexity may be formally declared O(2^n).

Another method of determining the computational complexity would have been sheer analysis. The iterations of a function exponentially increase by the amount of recursive calls. As the Fibonacci function calls itself two times per iteration, the analysed complexity would be O(2^n), which is roughly equal to our non-rounded solution.

The Mandelbrot Set, Part 2: OpenGL Program

In the last part I tried to give a general overview of how the Mandelbrot set is generated and how it can be rendered using a simple algorithm. This part will focus on programming a simple realtime renderer that can also zoom and move around, as I showed in this video. Additionally, it will serve as small introduction to modern OpenGL programming.

In this tutorial we work with C++, but you may use any other language you feel comfortable with that also binds with OpenGL.

You can also view these sources on Github.

1. Installing Prerequisites

Depending on your choice of programming language these instructions differ, I suggest you to search for them on the internet.

You need: a C++ compiler, GLEW, GLFW3, and OpenGL libraries. Install them with a package manager depending on the operating system you run.

A rough explanation of the libraries:

  • GLFW3 is an OpenGL framework that can create graphical contexts across multiple systems. It handles window functions including input.
  • GLEW (The OpenGL Extension Wranger Library) gives us greater access to hidden OpenGL extension features that may not be accessible with the standard API.

Now we’ll start coding. Open your vim/emacs/notepad/whatever.

2. Includes

3. Main function

Before we do anything in OpenGL we have to initialise GLFW. This takes only a few steps and is mostly straightforward.

We’ll declare three global variables, width, height, and the window-object itself. Then, we’ll initialise GLFW and set the appropriate settings and then finally also initialising GLEW.

The layout is not finished yet as we are missing some kind of render loop. The function glfwWindowShouldClose is a considerable halting-condition for our main loop right now. In the render loop, we clear the screen, swap buffers, and then wait for events (more on that later).

We can try compiling it now. Compile flags sadly differ on systems. Generally, Linux users can use this:

OSX users on the other hand have a slightly different method of linking OpenGL:

Explaining this for Windows would be too scary, so you’ll have to figure that out by yourself.

Screen Shot 2016-07-19 at 08.13.43

Starting this should leave you with a black screen and no output in the terminal. If there is an error, see if your graphics card supports OpenGL 4.1.

4. Creating the Shaders

To harness the full hardware-capabilities of your GPU we shall use shaders to render the Mandelbrot set at high-speed.

Shaders are little programs that can be directly executed on the GPU. They are therefore utterly efficient when it comes to processing data, as there is no latency between the CPU and the GPU.

Rendering requires a vertex and a pixel shader. As we are working with a two-dimensional plane only the pixel shader is of importance to us and I will not give much effort explaining how vertex shaders work, nor will you need to know how they work.

The creation of the vertex shaders as follows: (before main loop)

Much like other software, shaders are essentially written in source code too (in a C-like language), then passed to OpenGL, which compiles it to byte-code that the GPU can execute. All we need to care about is the handle of the shader as determined by the variable vs.

What matters here is the function

which assigns the source code to the shader handle to be compiled with glCompileShader. If length is NULL, then OpenGL expects the source string to be null-terminated.

What is left now is the pixel shader (called fragment shader in the OpenGL API). It will be responsible for rendering the Mandelbrot set, depending on which parameters (which are called uniform variables in OpenGL) we give to the pixel shader.

As you may have previously noted, the source code of the vertex shader is hard-coded into the program. But we will not hard-code the pixel shader, as we may want to play around with colours. The vertex shader will be kept static.

Continuing after the vertex shader:

The program will now open a file shader.glsl (containing the pixel shader code) and create the pixel shader, and additionally check for compile errors. Note that by calling glGetShaderiv we can receive an integer determining the success of compilation. If the compilation was not successful, we let OpenGL supply us with an info log that will be helpful for discovering the error.

Inside of shader.glsl:

we find a language similar to C. If you can recall the pseudo-code which I presented in the first part, you may find some resemblance.

We have at the top a macro defining the GLSL version, the output variable (which is a pixel), and the sum of input variables (uniforms). Note that vec4 is a vector with four components, which, in this instance, are holding information about a RGBA-pixel.

In the main function, we use the Mandelbrot algorithm to determine the value of the pixel we must return. Note that gl_FragCoord is a built-in variable of OpenGL which determines the currently processed pixel’s screen position. If you have a hard time following the structure of this program I recommend you to review the first part, where I explained in detail the Mandelbrot algorithm.

At last, to put our shaders to use we have to link them by creating a shader-program. This would technically allow us to use multiple pixel and vertex shaders but for now we only use two.

We can now finally activate the shader by calling the function:

5. Setting up the VBO and VAO

Before we can do anything with the Mandelbrot set, we will have to fill the screen with a single rectangle. The pixel shader needs a flat 2D-plane which the pixel shader will render to.

A vertex buffer object (VBO) stores vertex-data on the GPU. This allows OpenGL to save time by using the hardware of the GPU to its fullest potential. In our case it is neither necessary nor any more efficient as we only have very little data, but in the modern OpenGL specification there is no other way of rendering something. Also, using shaders requires the use of VBOs.

Let’s declare an float-array before our rendering loop.

These values determine the position of two triangles, which together fill up the whole screen.

Take a look at the coordinate system of OpenGL (which we could actually alter, but that is outside the scope of this tutorial).

ckla6

The center is at (0, 0), while one half of the width and one half of the height have a length of one, such that the corners of the screen rest at the upper left corner(-1, 1), the upper right corner (1, 1), the bottom left corner(-1, -1), and the bottom right corner (1, -1).

If you now look back at the points we declared you see that each pair of three values (X, Y, and Z; as the we render a 2D-plane the Z value is the same for every corner) signify a corner on the screen which sum up to two triangles, covering the whole screen. Why not just make a rectangle? you may ask. The modern OpenGL standard requires the use of triangles, therefore we could not have just made a rectangle.

The creation of the vertex buffer object as follows:

Again we need to use a handle to create the VBO. We then bind the buffer and upload the point data to it.

The function

uploads the data from RAM to the GPU. For the size we pass 2 * 9 * sizeof(float) as we have two pairs of three XYZ coordinates. The usage flag determines the memory location on the GPU. As the data will remain static (it would make no sense to alter it) we may pass the flag GL_STATIC_DRAW.

Now we have the data, but OpenGL still is unable to work with it. It does not know the structure of the data. OpenGL therefore introduces another object called the vertex array object which opengl.org defines as:

…an OpenGL Object that stores all of the state needed to supply vertex data… It stores the format of the vertex data as well as the Buffer Objects… providing the vertex data arrays.

Creating the VAO:

Here we have our handle vao which we bind, and set which the appropriate attributes for the data-structure.

The function which is of significant importance here

determines the format of the data. As our data consists of only XYZ-coordinates, we assign an index of zero (meaning the position attribute), a size of three, and a type of GL_FLOAT. Stride is set to zero as there is no margin between the coordinates.

6. Rendering

After we have struggled with setting up all the annoying OpenGL stuff we find we can now finally tackle the render loop and thus incorporate into the program what we have established (I must note that, when I was learning this by myself that my days were long and many nightmares were haunting me in my sleep).

Our final render loop:

Here we pass the global variables we defined to the shader at each frame. Then, we also bind the VAO and render it, specifying a triangle primitive, and a vertex-count of six (as each triangle has three corners/vertices).

At last, the fruits of our labour:

screen-shot-2016-11-05-at-19-33-47

If you happen to encounter problems, please stop and review your code, as it has been a lot of information. Also compare your code with my original source code.

In the next part, I will explain how to further improve the program, as to allow zooming, movement, and such.

Codesign gdb on OSX

Screen Shot 2016-07-06 at 16.58.50

Trying to use gdb on OSX El Capitan gave me following message.

Solution?

You need to codesign gdb with a custom certificate.

  1. Open Keychain Access
  2. In the menu, go to Certificate Assistant -> Create Certificate
  3. Choose an appropriate name like gdb-cert, set Identity Type to Self Signed Root, set Certificate Type to Codesign, and check Let me override defaults
  4. Optional: extend duration
  5. Click continue until specifying location for the certificate, choose System
  6. In System, click on gdb-cert, choose Get Info, go to Trust and set Code Signing to Always Trust
  7. Restart taskgated
  8. Codesign gdb

The Mandelbrot Set, Part 1: Overview

Recently, I have been experimenting with the Mandelbrot set. I made a video where I showcased a simple Mandelbrot explorer, written in C++ and OpenGL. In this little series I want to explain what the Mandelbrot set really is, how to generate it, and how to program it in OpenGL!

What is the Mandelbrot set?

The Mandelbrot set is a set of points, which are determined by a simple iteration. For a certain amount of iterations, this equation is repeated for each points in the set and coloured appropriately.

The iteration is

z_{n+1} = {z_n}^2 + c

where z and c is a complex number, and c is chosen according to the point in the set. The difficulty is that z and c are complex numbers. If they were not, we would not have this much variety in the Mandelbrot set.

A complex number consists of two parts, a real and imaginary part. They form an extension to the real numbers, and solve the awkward case of taking the square root of negative numbers. They are written by summing the real and imaginary part:

z = a + bi

where i is the imaginary number and defined as

i = \sqrt{-1}

You see that the number i does not really exist. That is why it is called the Imaginary Number.

When plotting complex numbers, the x-axis represents the real part, and the y-axis represents the imaginary part. In this way, numbers are more like points on a surface or two-dimensional vector.

2000px-Complex_number_illustration.svg

Similar to that, the Mandelbrot set is also a plot of each point on the graph, where the aforementioned iteration does not diverge.

Mandelset_hires

By looking at the graph, you may notice that c in fact is very small and lies close to the origin. We have to keep this in mind in order to render it correctly later on.

Additionally, all points in the set do not have an imaginary or real part greater than 2. We can use this fact to determine whether a point has escaped the Mandelbrot set. Given that, one has only to compute the absolute value of z, which is

|z| = \sqrt{a^2 + b^2}

Implementing an algorithm

It is very easy to program the Mandelbrot set. In the following you will see some pseudo-code that I mostly wrote to be easy to understand. If the concept has not clicked yet, maybe now it will!

Colouring

Methods of colouring the Mandelbrot set vary a lot. In general, you must use a palette and colour each point depending on iterations taken (iterations / max_iterations). In the pseudo-code that is the map_to_color function.

For example, if you want to colour the Mandelbrot set in gray-scale, you may write something like this.

In this way the RGB components scale linearly with t, resulting in a smooth black-white graph.

Screen Shot 2016-07-06 at 15.39.06

I hope I helped you grasp the general concept of the Mandelbrot. In the next part I will go on about how to implement this algorithm in OpenGL.

An Overview of the Entity-Component-System

Game programming traditionally revolved around the concept of inheritance, that is, game objects and entities inherit functionality, and become more specific in their functioning as the inheritance gets more complex.

But this type of architecture has evolved to be, in most cases, inefficient. As the game matures, functionality becomes increasingly complicated, which often results to deeply rooted chaos. To combat this cause, game programming has shown that the concept of Composition over Inheritance in such matters have turned out to be (sometimes) more efficient than an architecture based on inheritance.

The point of composition is simple. Rather than inheriting functionality, entities should have the functionality. In such a way behavior is more dynamic than inheritance, because components (which define that functionality) can be added, removed, and interchanged at runtime, whereas inheritance is fixed (you cannot change your parents). These two different concepts are also called has-a and is-a.

Basic Structure

In theory the Entity-Component-System is very simple. As the name says, the Entity-Component System, let’s call it ECS for now, is based on three things: entities, components, and systems.

Entities are things that live. They are the owner of components and are defined only by them. Entities are sole components containers, much like organisms are containers of genes – it is up to the genes how that organism behaves. Super power of programming: we can change the genes!

Typically, an entity would be defined by a unique ID. This would remove the problem of passing pointers, for they could get invalidated after an entity ceases to exist. Additionally it makes it easier to save entities, say, in a save-state file – there would be no need to reconstruct pointers.

Components are simple data-structs which hold some type of data relevant to its purpose. As said before, they could be called the genes of entities. Examples include a PhysicsComponent with coordinates and velocity. Or a HealthComponent which tracks health, etc. The basic idea of components is that they act as storage for relevant data, and nothing more. Combined, they build up the essential functionality of entities, an organism.

Say, you have an entity that builds an NPC, that would have a PhysicsComponent, FriendlyAIComponent, HealthComponent, etc. Now, with the magic of ECS, it is easy to turn this NPC to a zombie. How? Just by swapping FriendlyAIComponent with ZombieAIComponent!

Finally, the concept of a system is what actually brings these data structures alive. Systems make things move. They interact with components according to their data. For example, one may have a PhysicsSystem, HealthSystem, etc. Per frame, these would iterate through existing entities and update their specific components. The PhysicsSystem would for example add velocity to position, the HealthSystem might check for entities which have died, and clean up the entity, or trigger some kind of health regeneration. It is up to the designer how specific the systems are involved.

Summary

One of the major drawbacks of ECS is the tremendous amount of iteration one has to do. Then, as more entities and components are added, the game will get slower. Compared with other design architectures, ECS has no substantial benefits, other than its dynamics in handling components. In the end it is clearly a matter of choice. I hope I gave you an idea on ECS and how clever it is!

The Newspeak of the Digital Age

You probably know that I wrote about Emoji before, and while that post was mostly a rant discussing the mere ridicule of Emoji, my critiques of the actual effects of Emoji have not diminished and I would like to continue the discussion.

My prediction said that Emoji narrow our range of thought by restricting expressions to simple symbols that cannot yield much variety when combined. And as we continue the practice of using them in daily conversation there is a loss of deeper language.

In 1984, the classic dystopian novel, that is exactly the case. The Party (the ruling class of society) administers a new language called Newspeak, which is set to replace Oldspeak (Newspeak for English). The goal of Newspeak: to minimise expression of thought by restricting the range of words and sentences one can make use of. Why do that? To keep the Party in power of course!

For example in Newspeak there are no antonyms such as good and bad, but rather good and ungood. The goal of the language is to minimise vocabulary as such to make an uprising against the Party impossible.

In the end we shall make thoughtcrime literally impossible, because there will be no words in which to express it. Every concept that can ever be needed will be expressed by exactly one word, which its meaning rigidly defined and all its subsidiary meaning rubbed out and forgotten. […] Every year fewer and fewer words, and the range of consciousness always a little smaller.

Is that the case for us as well? Just recently Facebook introduced “reactions”, which allow you to respond to a post by choice of six icons. Of course comments are still allowed, but isn’t this the exact path of communication that George Orwell tried to warn us from? Some Emoji are actually banned in some countries, for example Emoji that depict homosexual couples, take Indonesia as example. It is an actual political problem nowadays. Emoji are more easier to control than normal language.

1984-movie-bb_a1

Especially because the ones who control Emoji are the big corporations such as Microsoft and Apple, and you cannot trust these corporations. Symbolically, they resemble the Party, they control everything, Big Brother is watching you.

I don’t want to appear paranoid here but rather simply raise conscious of these things that nowadays appear so integrated into our daily lives.

Newton’s Method: Explanation and Example Usage

Polynomials are compex to solve. Sometimes they can be solved quite easily and sometimes it seems like solving them is impossible. Fortunately there exists something called Newton’s Method which lets you approximate roots of a polynomial by using the Taylor Series.

The Taylor Series is defined as follows

f(x) = f(a) + f'(a)(x-a) + \frac{f''(a)}{2!}(x-a)^2+ \frac{f'''(a)}{3!}(x-a)^3 + \dots

By using the first two terms of the Taylor Series, we can approximate a function at some differentiate point x_0.

f(x) \approx f(x_0) + f'(x_0)(x - x_0)

You need to first find an approximation point to where a root could possibly be, for example by looking at the function’s graph or by iterating through some values of x and watching for close values to zero or changes of sign.

For example let’s take the function f(x) = 2x^3+5x^2-2

graph1

You can see that there is some root between 0.5 and 0.75, we can take that as approximation.

Screenshot 2016-03-12 17.22.19
Let x_0 = 0.5 be our initial guess. Then, let g(x) be our new function that will get assigned the approximated function that will be solved later.

f(x) \approx g(x) = f(0.5) + f'(0.5)(x - 0.5)

Then we find that

g(x) = 2(0.5)^3+5(0.5)^2-2 + (6(0.5)^2+10(0.5))(x - 0.5) = 6.5x - 3.75

If we plot g(x) you can see that it is tangent to f(x). Furthermore it has a root which is very close to the root we want to calculate. Dazzling!

graph3

Now, if we solve g(x) we get x_1 = \frac{15}{26} \approx 0.57. If we zoom in a little more, we find that the approximation turned out to be quite effective with our error of 0.06!

Screenshot 2016-03-12 17.25.09

If you repeat this process again with x_1, you find that you get a new function h(x)=\frac{2625}{338}x + \frac{85575347}{19307236} which has a root at x_2 \approx 0.5707. As seen, the method is quite effective and converges very quickly.

Iteration

Since we want to iterate this mathematical process with a computer, we might find it simpler if we simplify the equation. We can generalize the equation by noticing the pattern, that is

0 = f(x_n) + f'(x_n)(x_{n+1} - x_n)

when solved for x_{n+1} becomes

x_{n+1} = -\frac{f(x_n)}{f'(x_n)} + x_n

Which is in terms of programming much clearer to our eye than the mess above.

Example: Calculate the square root of a number

Getting the square root of some number  k just means solving the function f(x) = x^2 - k. Since we find the derivative to be f'(x) = 2x we find that implementing such program to be pretty easy.

First, notice how quickly Newton’s Method converges for solving k=2.

Step Starting at 2 Starting at 10 Starting at 100
1 1.5 5.1 50.1
2 1.4166666667 2.746078431 25.024996
3 1.4142135624 1.4442380949 12.55245805
4 1.4142135624 1.4145256551 6.355894
5 1.4142135624 1.4142135968 3.335281609

As you can see Newton’s Method converges pretty fast, even if your approximation is absurdly dumb. Such algorithm can be implemented with ease as seen below (C code).

The program though turns out to have big errors for numbers bigger than 10. Notice how bigger the number the worse the approximation? This happens because our iteration count and approximation is only efficient for very small numbers.

Screenshot 2016-03-14 12.50.20

There are a multitude of improvements that could be done such as making the iteration count depend on the number as well as trying to find a good linear approximation to the square root. I won’t discuss these methods here since my only ambition was to introduce you to Newton’s Method and give you an idea of how it works and how awesome it is.

Trailing Zeros of Factorials

Since I’m bored to hell in programming class, my teacher, who is also my physics teacher, gave me the task to find out the trailing zeros of 100! and then for n!. I then began searching for patterns in factorials starting from 1! to 10!. Interestingly, 5! = 120 has one trailing zero and 10! =3628800 has two trailing zeros. By induction you could determine the trailing zeros of factorials to be expressed by \lfloor \frac{n}{5} \rfloor, where \lfloor x \rfloor is the floor function.

But if you go down the pattern a bit further you find that 25! =15511210043330985984000000 which has six zeros! The heck? This contradicts our idea. Let’s investigate a bit further.

In a factorial, the amount of fives getting multiplied adds a zero. For example 10! = 1\times2\times3\times4\times5\times6\times7\times8\times9\times10. If we take a look we see that there are two 5s hidden here, 5 and 10=5\times2. But 25! contains the factor 25 which is 5^2! So in that case, we find that there are six fives, because 25=5\times5.

With this additional information we can determine the number of trailing zeros of 25! to be \lfloor \frac{25}{5}\rfloor +\lfloor \frac{25}{5^2} \rfloor = 6

But we are not finished yet, there is a pattern here and we should be able to define it for any n!.

Given a number n, the trailing zeros of n! is the sum of n divided by all of its prime factors of 5.

\big f(n) = \sum\limits_{i=1}^{k} \lfloor \frac{n}{5^i} \rfloor

where k has to be chosen such that 5^{k+1} > n.

For example let’s calculate the trailing zeros of 100!. We find 5^3 = 125 > 100 leaving us with f(100) = \lfloor \frac{100}{5}\rfloor +\lfloor \frac{100}{25}\rfloor = 24.

Blogging about a variety of topics. Mostly Software and Mathematics.