Browsing Tag

# math

A few months back I posted some code on GitHub to convert Latex markup code into PNG files using Python named it PyTex2png. I know realize that I didn’t explain one bit of the code. I did however write a very clear (to me) README.md file explaining how to use the code if you wanted to download and use it. I even included a couple of cool examples. In this post I want to dive into the code that make it work.

Lets start by explaining exactly what Latex is and what this module does. Latex is a markup language used for scientific notation and papers. I came across it in my second year in college and science then I have never looked back. It is by far my favorite word processor. I even have my resume and letters typed up in Latex. It just make pretty things. Among other things, Latex can create beautiful Mathematical typography like the one below. Latex does have a bit of learning curve, but it shouldn’t be too difficult to master. We are only looking at the math portion of it anyway. Here is some Latex code and the image that got generated from it. Note that this post will not teach you Latex. There are plenty of books and information about Latex online already.

Code:

\cos (2\theta) = \cos^2 \theta - \sin^2 \theta \\
\lim_{x \to \infty} \exp(-x) = 0 \\
a \bmod b \\
x \equiv a \pmod b

[\code]

Pretty Output:

So how does that happen? Well there are a few things to look at. The first is a C++ code that was written by Bruno Bachelet. It is available here (and here as backup). There is also and online version of the code running here. This code is the main power behind the software. Everything else is just a Python wrapper to add some fancy functions like multiple files, remove background and invert colors.

Within pytex2png.py you will see a function name covert(). This function calls the C++ module to convert images, creates the PNG, removes any background and finally converts the math to white. The last portion is not really required. I was coding this up for a larger project and in that one the background is not white, so I need white fonts.


def convert(source,output,display=False):
if(display): print &amp;amp;quot;Coverting LaTeX from &amp;amp;quot; + source + &amp;amp;quot; to &amp;amp;quot; + output
data = get_file(source).replace(&amp;amp;quot;\\&amp;amp;quot;,&amp;amp;quot;\\\\&amp;amp;quot;)
make_png(data,output,display)
make_transparent_bg(output,display)
make_white_text(output,display)



This doesn’t tell us anything about how this program works. It is just a plain function that makes calls to other functions. We will have to look at the other functions to understand how this one ties everything together. Lets take a look at the code:


# read file content as a single variable
def get_file(filename):
file = open(filename, "r")
file.close()
return lines

# create png from latex
def make_png(data,image,display):
if(display): print "Making image..."
command_line = "./tex2png -r* -lq \""+data+"\" " + image
exe_command(command_line,display)

# remove white background
def make_transparent_bg(image,display):
if(display): print "Removing white background..."
command_line = "convert "+image+" -channel rgba -fill \"rgba(255,255,255,0)\" -opaque \"rgb(255,255,255)\" " + image
exe_command(command_line)

# convert output text to white
def make_white_text(image,display):
if(display): print "Converting output image text to white..."
command_line = "convert "+image+" -channel rgb -fill \"rgba(255,255,255)\" -opaque \"rgb(0,0,0)\" " + image
exe_command(command_line)



After a careful inspection of these functions, you should see a pattern. The first one is a get_file function. All it does is open a file and reads it’s content. Then returns the data from the file as a single object. You can read more about reading a writing from files in the Python Documentation. The other 3 functions, make_png, make_tansparent_bg, and make_white_text follow a templet:

display output
build command line command
execute command

That’s all folks! it is pretty straight forward. The commands each one execute are based off of what we need done. The first one, make_png calls the executable made from the C++ module. The remanning 2 functions, make_transparent_bg and make_white_text use a command line tool called convert. This is a Unix utility and you can read all about it in the convert manual file or on this site.

You will notice that there is one function we haven’t mentioned yet, exe_command. Let’s look at that function:


# runs a command line command and waits for it to complete
# by default it will not display any output
def exe_command(cmd, display=False):
args = shlex.split(cmd)
if(display):
p = subprocess.Popen(args)
else:
p = subprocess.Popen(args,stdout=subprocess.PIPE)
p.wait()



This one should be taken as is. It takes in a command line command and excites it using a Python module called subprocess. It can be used to execute and command in the command line and can retune the results if needed. It is a good little function to keep handy for other programs in the future...

That is all the magic that ties this thing together. There are some more files and notable mentions.

• Examples folder holds text files with Latex markup to be converted.
• Output folder will contain PNG’s that have been converted already.
• Source has the C++ file need to make the PNG files in the first place.
• Makefile is used to compile the C++ file and create and executable for your computer architecture.
• examples.py is a simple Python module to load the example text files from the folder and run them through PyTex2png converter.
• gpl-3.0.txt is a GPL License under which the code is distributed.

Here is a PyTex2png zip file with all the files together. It is also available via GitHub (pytex2png).

For the most part, all of us take for granted that the cpu can do multiplication in a timely matter. However, there are some systems that do not have multiplication. Furthermore, how is multiplication actually done? Multiplication is nothing more than repeated applications of additions. For example, if we want to compute 5 x 4 we can add 5 4 times or add 4 5 times and keep track of the sum. Let’s take a look how this might look like in Python 2.7.


import sys

def main(argv):

if len(argv) != 2:
sys.exit('Usage: simple_multi.py <a> <b>')

a = int(sys.argv[1])
b = int(sys.argv[2])
sum = 0

for i in range(a):
sum += b

print "The result of " + str(a) + " times " + str(b) + " is: " + str(sum)

if __name__ == "__main__":
main(sys.argv[1:])



And the output:

Sweet. Now for most of you this will be it. Trust me, for me the multiplication symbol is it. I never thought I will have to code multiplication. That is until it came up. I have a friend who happens to be a mathematician with a particular fascination to numbers and counting. So he came up to me claiming that he has this fast way to do multiplication that has to do with the ten’s places and basic multiplication. Basic multiplication being the 10 x 10 grid (or 12 x 12) we all had to memorize at some point. So I decided to put him up for the test, I wrote his algorithm in Python as well as the one above and Russian Multiplication. I won’t spend time trying to explain the algorithms. If you want to, you will find the link helpful enough to explain Russian Multiplication. If you want to find out more about my friend’s Ten’s Multipication, send him a message on his blog DeadEndMath. All in python, running on the same computer with the same Linux time command measuring their performance. Before we get to the data, let’s take a look at the code:

Ten’s Multiplication:


import sys

def main(argv):

if len(argv) != 2:
sys.exit('Usage: tens_multi.py <a> <b>')

a = sys.argv[1]
b = sys.argv[2]
sum = 0;
c = len(a) - 1;

for i in a:
d = len(b) - 1;
for j in b:
sum += int(i)*int(j)*(10**c)*(10**d);
d -= 1;
c -= 1;

print "The result of " + a + " times " + b + " is: " + str(sum)

if __name__ == "__main__":
main(sys.argv[1:])



And the output for verification:

Well, at least we got the same results. Let us look at the Russian Multiplication in Python:


import sys

def main(argv):

if len(argv) != 2:
sys.exit('Usage: russian_multi.py <a> <b>')

a = int(sys.argv[1])
b = int(sys.argv[2])

sum = 0;

if a == 0 or b == 0:
print 0;
exit();

if b%2 == 1:
sum += a;

while b != 1:
a = a*2;
b = b/2;
if b%2 == 1:
sum += a;

print "The result of " + str(a) + " times " + str(b) + " is: " + str(sum)

if __name__ == "__main__":
main(sys.argv[1:])



The output:

Again, good thing we got the same results. Now we know we have all three versions of multiplication and they are all computing products of 2 numbers correctly. Now let us look at their time compression. I have organized the data in a table for convince:

 Simple Ten's Russian 4 x 5 0.013 0.013 0.014 22 x 54 0.014 0.014 0.014 5142 x 9810 0.014 0.014 0.014 716234 x 847263 0.040 0.014 0.014 817263548 x  827461930 NaN 0.014 0.014

So what can we see? Well first of all the data gathered is from the linux time command. When execute any program in a Linux shell, you can add the command time before the program. This will return you with an output similar to the one below after the command is complete. the times shown are subjective to your computer specifications and current status (i.e. running programs). Here is an example of a Linux time output:

So what does all this mean? The real time is total execution time. The user time is how long it took the user to key in the input. The sys time is the time your command was running on the cpu. In the table above, the time was taken from the sys value of execution of each command. You will notice that simple multiplication gets more cpu consuming as the length of the input increases, while the time for Ten’s and Russian Multiplication remained constant. This means that either algorithm could work as a substitute to multiplication if you happen to work in an environment that does not support multiplication natively. As a final note let me add that even a basic multiplication in Python operates within the same time as Ten’s and Russian multiplication, if not a second faster.

How math is done

Over the past few days I was debating what I should write about next. I really enjoy sharing my thoughts and have been getting great responses from my readers. I thanks you for your comments. At your request I have decided to cover another math topic, solving equations using python. Equation solving is something fundamental in computing, mathematics and other scientific applications. Using computers enables users to solve a wider range of equations than we could by non-computer methods. However, there is a cost at hand, accuracy. Although computers can solve equations using a few methods, we will witness that accuracy will change between methods. In some cases you might not get to a solution at all.

Before we dive into code behind solving equations, let us step back. A good programmer should always consult other resources to see what else has been done already. The Internet (and Google) have given us great tools to research any topic our heart desires. There are many tools available to solve equations and do almost any other mathematical calculations. Matlab and Matamatica are both great mathematical softwares that can do amazing things. However, they are paid distributions. I do not have anything against paid software, but for the sake of educational useage I would like to keep everything in the realm of free (open source preferred). I like to post code that anyone can download, run, modify, distribute or do anything else they would like to do with it. Using paid languages will be counter productive to such goals and limit the readers of this blog. That being said, if there is a request I will put up Matlab code for the algorithms discussed in this post.

On another note, we will be using Python2.6 for all the examples that are to follow. We will not be using NumPy or SciPy. Both are great scientific math tools but an over kill for what we are trying to do. Our mission is simple, solving an equations. More precisely, we want to figure out where a given function intersect the x-axis without using “the big guns”, including Wolfram-Alpha. So let us get to it.

The Bisection Method:

At out first stop, the most fundamental approach to solving equations is the Bisection Method. Here is the theory in a nutshell:

“If f is a continues function on [a,b], where f(a)f(b) < 0, there must be a point r between a and b where f(r) = 0”

So what does this means? It means that if we have a function that cross the axises between point a and b, then the value of the function at f(a), multiplied by the value of the function at b, f(b), must be smaller that 0. Further more between a and b there is a point where the function f is 0. Here is a picture describing such event:

Here is a general solution in code that does exactly that:


def bisection(a,b,tol):
c = (a+b)/2.0
while (b-a)/2.0 > tol:
if f(c) == 0:
return c
elif f(a)*f(c) &lt; 0:
b = c
else :
a = c
c = (a+b)/2.0

return c



Notice that we have a while loop that runs while we are within tolerance and that we assume function f(x). The latter part can be solved by writing a function f that outputs the desired function. Here is the function we are going to work with:

And here is the python equivalent:


def f(x):
return x**3 + x -1



Once we add the code together we can run it and get a solution within tolerance. What is the tolerance? The Bisection method literally cuts in half the solution at each iteration, assuming the mid point between a and b is the solution. So in theory we can keep asking for a more accurate solution by setting a lower tolerance. The “real” solution to the function according to Wolfram Alpha is 0.6823278... We will use this solution as a guide to compare the accuracy of solutions we compute. Here are a few outputs for our code, for a solution between 0 and 1 using various tolerances:

(The full code for the Bisection Methon in Python)

As you can see, the first solution is only accurate to 2 decimal places. If we tried a lower tolerance we would have got maybe one or none of the digits correct. When we increased our tolerance we arrive to a more accurate solution, within tolerance range. Bisection is a good method to use when it applies. The Bisection method will not work on all functions, but is a quick solution to many equations.

Fixed Point Iteration (FPI)

Fixed point iteration is another method of solving equations. The definition of a fixed point r is when g(r)=r. Using this definition, we start with some initial guess, x and plug it into the function. However, unlike bisection, this method does require us to do some math work before code. We have to transform our function to be in the form x=g(x). This is where the heart of FPI is. Based on your setup and initial guess, you may reach a solution really fast, or not at all. Here is how the code for FPI looks like:


def fpi(x, k):
for i in range(k):
x = g(x)
return x



You will note that now we are no longer using a while loop. FPI could be setup in a smilier way to Bisection, but for simplicity is setup in a fixed i-loop. Using the same function as before, let us say we transform the function to be like this:

And in python (using python’s math library):


def g(x):
return math.pow(1-x,1/3.0)



You will notice that python only has a square root function. For any other roots we use python’s power function. Please remember that 1/3 in integer division is 0. Therefore, we have to use 1/3.0 to get the cubic root. Running the program for 10 iterations starting at 0.5 yields the following output:

(The full code for Fixed Point Iteration in Python)

This is not very close, we only have 1 digit and we started close to the solution. If we let the program run for 100 iterations we get really close to the solution:

How many iteration do you need? It depends on your function. A variation on this function might solve it in 7 iterations (can you find the function?). Again, in FPI it depends on the setup and the initial guess. As a rule of thumb, it will either conversation very fast, or not at all. But then again, sometimes it just needs one or more steps.

Newton’s Method

Our next stop, and the basis for our next method, is Newton’s method. In this case, like FPI we start with an initial guess and keep redefining x, however, we use this function to find the next x :

Newton’s method adds another math task we will have to do outside of our code, taking the derivative. We will see that the last method will take care of this, but at another cost. Here is how Newton’s method will look like in code:


def newt(x,n):
for i in range(n):
x = x - f(x)/f_prime(x)
return x



This means we will have to define our original function as well as a derivative (prime) function. Not a problem, we can use the same function we have been solving all along. This is what it will look like:

And in Python:


def f(x):
return x**3 + x - 1

def f_prime(x):
return 3*x**2 + 1



Before we run our code, we need to take care of one special case. What if the derivative is 0 or end up computing to 0? Than we will have a division by zero error. This can be simply avoided by a couple of if statements. Like so:


def newt(x,n):
for i in range(n):
if f_prime(x) == 0:
return x
x = x - f(x)/f_prime(x)
return x



This should take care of any division by zero errors. However, do not take my word for it, try it out for yourself!

And now running the code, using -0.7 as our initial x for 10 iterations will yield:

(The full code for Newton’s Method in Python)

If you add some debug statements and examine the output you will note that the solutions was found after only 7 iterations. That is really fast and a lot better than 100 iterations of FPI. You might have output like this:

Secant Method

Although Newton’s method works really well, taking derivatives is not always an option. For that, there is another method, mathematically based on Newton’s method and the definition of derivatives. This time we start off with 2 initial guesses and move forward by using one of the guess and a new point. The new point is calculated from the formula:

In Python the Secant Method might look like this:


def secant(x0,x1,n):
for i in range(n):
x_temp = x1 - (f(x1)*(x1-x0)*1.0)/(f(x1)-f(x0))
x0 = x1
x1 = x_temp
return x1



Note that there is a chance that we will be dividing by 0, again we must add a check to ensure that this does not happen. Also note that we use a temp variable to swap values around. This is so we do not require much memory. A new, modified version of the Secant method might look like this:


def secant_2(x0,x1,n):
for i in range(n):
if f(x1)-f(x0) == 0:
return x1
x_temp = x1 - (f(x1)*(x1-x0)*1.0)/(f(x1)-f(x0))
x0 = x1
x1 = x_temp
return x1



As before, we are using for loops, which might be swapped to while loops. I recommend you try it your self, but remember, while loops must terminate at some point. Running the Secant code starting with 0 and 1 for 10 iterations, on out “famous” function yields the following output:

(The full code for Secant Method in Python)

We have gone over 4 basic methods of solving equations. There are many more approaches to solve equations, this is by no means a complete survey. I trust that if you are in need for more information you can find it using a friendly search engine. Note that there is still much to say about the accuracy of the solutions computed by these methods. Also note that there are some programming “tricks” that can be done to shorten calculation. For more accurate mathematical definitions and information please seek out a mathematician.

This post was made possible by Numerical Analysis by Timothy Sauer

To mathematicians prime numbers represent a certain Nobel challenge. They have been studying prime numbers for a while, trying to learn their properties. Some mathematicians actually made progress. Others have gone insane trying to prove a pattern or how to generate the nth prime. Some are actively looking for the largest prime they can find. For computer scientists and programmers prime numbers have properties that can be utilized in many areas. Cryptography and pseudo-random number generators are 2 examples applications of prime numbers. Let us take a dive into the world of prime numbers. We will look at how we can generate and validate prime numbers in python.

Any valid math discussion must begin with a definition. Here is a definition of prime numbers I like to use:

“A natural number p > 2 is called prime if and only if the only natural numbers that divide p [without reminder] are 1 and p. A natural n > 1 that is not prime is called composite. Thus,, n > 1 is composite if n = ab, where a and b are natural numbers with a > 1 and b < n. “ Discrete Mathematics with Graph Theory

What does this means? It means that a number, for example 17, is prime if and only if it can be divided only by 17 and 1. On the other hand, 18 is not prime. 18 is a composite number that can be written as 2 x 9 or 3 x 6. This will prove to be a very useful property.

So how do we find prime number? Well, we do not. Or not exactly. Since the natural number are infinite and a prime number can be any number greater then 2, by definition there are infinite prime numbers. Many mathematicians have tried to prove a formula or correlation among prime numbers. None to date have succeeded. The best we can do is to generate a number and check if it is prime. However, there are some tricks we can do to speed things up, but we will get to that.

Here is our basic problem. We want to generate all the prime numbers between x and y. For starts, let us think about the numbers between 2 and 100. The straight forward, brute force, approach will be to check every number in that range and see if any number up to it can divide it without a reminder. If we find such number then it is not prime. Let us take a look at a simple program that can do this for us:


import sys

def is_prime(num):
for j in range(2,num):
if (num % j) == 0:
return False
return True

def main(argv):

if (len(sys.argv) != 3):
sys.exit('Usage: prime_numbers.py <lowest_bound> <upper_bound>')

low = int(sys.argv[1])
high = int(sys.argv[2])

for i in range(low,high):
if is_prime(i):
print i,

if __name__ == "__main__":
main(sys.argv[1:])



And the output for prime numbers between 2 and 100:

Now for some tricks. First, we don’t have to check every second number. All the even numbers will be divisible by 2 by definition of even numbers. So our first loop can be rewritten to jump by 2 and start at an odd number. That way we are just checking the odd numbers. By this hypothesis, all odd number are potentially a prime number. This takes away half the numbers we need to check. If we were to ask for all the prime numbers between 100 and 10,000 it might take a while. Lets look at how long it will take my computer (Intel core duo 2 2.53GHz with 8GB ram) to run this program using the unix time command.

The original script from before:

Modified script to take out the evens numbers:


import sys
import math

def is_prime(num):
for j in range(2,num):
if (num % j) == 0:
return False
return True

def main(argv):

if (len(sys.argv) != 3):
sys.exit('Usage: prime_numbers2.py <lowest_bound> <upper_bound>')

low = int(sys.argv[1])
high = int(sys.argv[2])

if (low == 2):
print 2,

if (low % 2 == 0):
low += 1

for i in range(low,high,2):
if is_prime(i):
print i,

if __name__ == "__main__":
main(sys.argv[1:])



This modification cut the program time by nearly half, but we can do better. According to the Sieve of Eratosthenes we do not need to check if all the numbers will divide our potential prime, we only need to check up to the floor of the square root of the number. So our next version of the code will look like this:


import sys
import math

def is_prime(num):
for j in range(2,math.sqrt(num)):
if (num % j) == 0:
return False
return True

def main(argv):

if (len(sys.argv) != 3):
sys.exit('Usage: prime_numbers3.py <lowest_bound> <upper_bound>')

low = int(sys.argv[1])
high = int(sys.argv[2])

if (low % 2 == 0):
low += 1

for i in range(low,high,2):
if is_prime(i):
print i,

if __name__ == "__main__":
main(sys.argv[1:])


And the time to find all the primes between 100 to 10,000:

This is really great. We have improved our program by another half and nearly a quarter than the original run. Let us examine 2 (2 things) variations of this program. First we will calculate the nth prime number. Second we will try to get a random prime number for a given length.

In order to calculate the nth prime we will have to modify them main() loop into a while loop. We will also modify the program to accept only 1 number, the nth prime we want to calculate. Once we reach the nth prime we will print it out.


import sys
import math

def is_prime(num):
for j in range(2,int(math.sqrt(num)+1)):
if (num % j) == 0:
return False
return True

def main(argv):

if (len(sys.argv) != 2):
sys.exit('Usage: prime_numbers4.py <nth_prime>')

i = 0
num = 2
nth = int(sys.argv[1])

while i < nth:
if is_prime(num):
i += 1
if i == nth:
print 'The ' + str(nth) + ' prime number is: ' + str(num)
num += 1

if __name__ == "__main__":
main(sys.argv[1:])



Sample output:

For many application computing a prime number is not enough. For some we want a specific length of prime number. Lets say we need a 5 digit prime number. How are we to calculate this? The simple approach is to generate all 5 digits prime numbers and choose 1 of them randomly. Let us look at how to do this:


import sys
import math
import random

def is_prime(num):
for j in range(2,int(math.sqrt(num)+1)):
if (num % j) == 0:
return False
return True

def main(argv):

if (len(sys.argv) != 2):
sys.exit('Usage: prime_numbers5.py <num_digits>')

digits = int(sys.argv[1])
low = int('1' + '0' * (digits-1))
high = int('9' * digits)
prime_number_list = []

if (low == 1):
prime_number_list.append(2)
low = 3

if (low % 2 == 0):
low += 1

for i in range(low,high,2):
if is_prime(i):
prime_number_list.append(i)

print 'A random ' + str(digits) + ' digits prime number is: ' + str(random.choice(prime_number_list))

if __name__ == "__main__":
main(sys.argv[1:])



Output:

If you were to run the program, you will notice that the it will take longer and loner as n gets larger. I found that for any n grater than 6 digits, it is very demanding on the hardware. I do not recommend letting this program try to calculate a high digit prime. You will also notice if you try to ask for 10 digits, python will raise an error in the range function. So how do we compute n digit prime for a large n? One way I found is by using a random number and checking if it is prime. This seems to work reasonably well with numbers less than 15, but then again it is random.


import sys
import math
import random

def is_prime(num):
for j in range(2,int(math.sqrt(num)+1)):
if (num % j) == 0:
return False
return True

def main(argv):

if (len(sys.argv) != 2):
sys.exit('Usage: prime_numbers6.py <num_digits>')

digits = int(sys.argv[1])
low = int('1' + '0' * (digits-1))
high = int('9' * digits)
done = False

if (low == 1):
low = 2

while not done:
num = random.randint(low,high)
if is_prime(num):
print 'A random ' + str(digits) + ' digits prime number is: ' + str(num)
done = True

if __name__ == "__main__":
main(sys.argv[1:])



So how about really large primes? That is where time becomes an enemy and computational complexity goes up. I found a few modules that can help compute large primes in python but I have yet to try to run anything. Once I get a chance I will try to run them and put a post up. If you really need prime numbers I would advise you to either download a list of them or generate a list of primes and save them to a text file. Prime numbers will show up every once in a while. They come up in cryptography and other encryption methods. Most usage of prime numbers is due to the nature of the numbers. For more information about prime numbers I would suggest seeking out a mathematician, which I am clearly not.