GilgaLab

8Jul/112

Project Euler – Problem 28

Posted by Henrique

Alright!
I have been away for a while as there are too many things happening right now and I have been pretty busy.
Well, this time I will resolve Problem 28 from Projet Euler. For this, I will use mostly a manual method for deducing the formulas and then a small C program to finish the sum for us :)
The problem states:

Starting with the number 1 and moving to the right in a clockwise direction a 5 by 5 spiral is formed as follows:

21 22 23 24 25
20  7  8  9 10
19  6  1  2 11
18  5  4  3 12
17 16 15 14 13

It can be verified that the sum of the numbers on the diagonals is 101.

What is the sum of the numbers on the diagonals in a 1001 by 1001 spiral formed in the same way?

According to the problem for this spiral the sum is 101 because: 1 + 3 + 5 + 7 + 9 + 13 + 17 + 21 + 25 = 101
Good! So in order to resolve this problem all we have to do is figure out what are the rules that guide the diagonal numbers.

I have divided this spiral in 4 diagonals:

  • Center to right upper corner
  • Center to right lower corner
  • Center to left lower corner
  • Center to left upper corner

If we expand this spiral a a little bit more, we will have:

73 74 75 76 77 78 79 80 81
72 43 44 45 46 47 48 49 50
71 42 21 22 23 24 25 26 51
70 41 20 07 08 09 10 27 52
69 40 19 06 01 02 11 28 53
68 39 18 05 04 03 12 29 54
67 38 17 16 15 14 13 30 55
66 37 36 35 34 33 32 31 56
65 64 63 62 61 60 59 58 57

Based on this we can deduce the formulas that generates each diagonal.

1. Center to right upper corner (1, 9, 25, 49, 81, ...)

For this one, we can see that for each level that expands, we have a squared number. So if we put a little bit of thinking on it, we can come up with the following formula:

      f(n) = (2*n+1)^2

Where n is the number of the line counting from the center up (The first line being 0)
You can try it!

      f(0) = (2*0+1)^2 = (1)^2 = 1\\      f(1) = (2*1+1)^2 = (3)^2 = 9\\      f(2) = (2*2+1)^2 = (5)^2 = 25\\      f(3) = (2*3+1)^2 = (7)^2 = 49

2. Center to right lower corner (1, 3, 13, 31, 57, ...)

Comparing the results from the center to right upper corner results, we notice that the value here is the same as that diagonal minus 6*n. So the formula goes:

      f(n) = ((2*n+1)^2) - 6*n


Where n is the number of the line counting from the center down (The first line being 0)
You can try it!

      f(0) = ((2*0+1)^2) - 6*0 = (1)^2 - 0  = 1\\      f(1) = ((2*1+1)^2) - 6*1 = (3)^2 - 6  = 3\\      f(2) = ((2*2+1)^2) - 6*2 = (5)^2 - 12 = 13\\      f(3) = ((2*3+1)^2) - 6*3 = (7)^2 - 18 = 31

3. Center to left lower corner (1, 5, 17, 37, 65, ...)

Comparing the results from the center to right upper corner results, we notice that the value here is the same as that diagonal minus 4*n. So the formula goes:

      f(n) = (2*n+1)^2-4n

Where n is the number of the line counting from the center down (The first line being 0)
You can also try this one!

      f(0) = (2*n+1)^2-4n = (1)^2 - 0 = 1\\      f(1) = (2*n+1)^2-4n = (3)^2 - 4 = 5\\      f(2) = (2*n+1)^2-4n = (5)^2 - 8 = 17\\      f(3) = (2*n+1)^2-4n = (7)^2 - 12 = 37


4. Center to left upper corner (1, 7, 21, 43, 73, ...)

Comparing the results from the center to right upper corner results, we notice that the value here is the same as that diagonal minus 2*n. So the formula goes:

      f(n) = (2*n+1)^2-2n

Where n is the number of the line counting from the center up (The first line being 0)
And once again you can try this one!

      f(0) = (2*n+1)^2-2n = (1)^2 - 0 = 1\\      f(1) = (2*n+1)^2-2n = (3)^2 - 2 = 7\\      f(2) = (2*n+1)^2-2n = (5)^2 - 4 = 21\\      f(3) = (2*n+1)^2-2n = (7)^2 - 6 = 43

Now we need to figure out from what value to what value our 'n' should range. We know the following about our grid:

  • When the grid is 3x3, we have 9 (3^2) at the right top corner
  • When the grid is 5x5, we have 25 (5^2) at the right top corner
  • When the grid is 7x7, we have 49 (7^2) at the right top corner
  • When the grid is 9x9, we have 81 (9^2) at the right top corner
  • When the grid is 1001x1001, we have 1002001 (1001^2) at the right top corner

Based on this we can conclude that we will have 500 lines up, 500 lines down and one line in the center, each line with 1001 rows.
Our sums should look like this then:

      \sum_{n = 0}^{500} (2*n+1)^2 \\      \\      \sum_{n = 0}^{500} (2*n+1)^2 - 6*n \\      \\      \sum_{n = 0}^{500} (2*n+1)^2 - 4*n \\      \\      \sum_{n = 0}^{500} (2*n+1)^2 - 2*n

We calculate them all and sum them all together. After this is done, we subtract 3 from the sum because in them all we have added the element in the middle of the grid (1), so we have summed the number 1 four times, and the exercise asks us to sum all the numbers in the diagonals, meaning we should add each one of them only once.
We can than use WolframAlpha to calculate this sum, or just write a simple C code for it:

#include <stdio.h>
 
int main()
{
    int n;
    unsigned long long sum = 0;
 
    for(n = 0; n <= 500; n++)
    {
        sum += ((2*n+1)*(2*n+1));       // Same as (2*n+1)^2
        sum += ((2*n+1)*(2*n+1)) - 6*n; // Same as (2*n+1)^2 - 6*n
        sum += ((2*n+1)*(2*n+1)) - 4*n; // Same as (2*n+1)^2 - 4*n
        sum += ((2*n+1)*(2*n+1)) - 2*n; // Same as (2*n+1)^2 - 2*n
    }
 
    printf("Result is: %llu", sum-3);
 
    return 0;
}

I will leave it up to the user to execute the above code to get the result, or to calculate the formulas using some tool.
Hope you enjoyed! :)

3May/112

Python C API – Second Step

Posted by Henrique

As promissed, here is the second step on playing with the Python C API. If you haven't read the first post, I strongly suggest you to do so: Python C API - First Step

This will be a much more exciting read with many new concepts. If you find something to be confusing or is having a hard time understanding it, please let me know and we can discuss!

For this post I will show a code that does many things, such as:

  • Look for the 'sys' module in the default loaded modules
  • Add a new path to the sys.path list in the sys module
  • Import a module wirtten by us
  • Access the members of the module that we loaded
  • Call a function from the module that we loaded
  • And more!

I hope that from this introduction you feel excited about what is coming ahead!

So, let's start looking at the module we will import:

''' Save this code in a file called example.py '''
 
def doubleValue(x):
    y = 2 * x
    return y
 
myString = "Hello World!"

And finally the code that will do all those wonderful things mentioned in the list above:

/**
 * Save the code in a file called main.c
 * Compile with one of the following:
 *    gcc -Wall `python3.2-config --libs` `python3.2-config --include` main.c -o pyx
 *    gcc -Wall -lpthread -ldl -lutil -lm -lpython3.2mu -I/usr/include/python3.2mu main.c -o pyx
 * 
 * Run with: ./pyx
 * 
 * Use the code as you wish. If any help is needed, please let me know.
 * If you are going to republish this code somewhere, I would be grateful
 * if you can keep the credits of the code and/or a link to the original.
 *
 * Author: Henrique M. D.
 * Email: typoon at gmail dot com
 */
 
#include "Python.h"
#include <stdio.h>
#include <stdlib.h>
 
int main(int argc, char** argv) {
    int i;
    char module_dir[255];
    char *name;
    char *str;
    //wchar_t *module_path;
    PyObject *poModule;
    PyObject *poDictModule;
    PyObject *poKey;
    PyObject *poValue;
    PyObject *poString;
    PyObject *poSys;
    PyObject *poList;
    PyObject *poListItem;
    PyObject *poPath;
    PyObject *poMyString;
    PyObject *poDoubleValue;
    PyObject *poResult;
    Py_ssize_t size;
 
    /* Start Part 1 */
 
    // This is the path where our .py scripts are. We will need to
    // add this path to the sys.path list later on prior to loading
    // it. If we don't do that, Python won't be able to load our script
    // as it won't find it.
    // Change this path accordingly to your setup
    memset(module_dir, 0, sizeof(module_dir));
    sprintf(module_dir, "/home/gilgamesh/codes/Pyx/scripts");
 
    Py_Initialize();
 
    poDictModule = PyImport_GetModuleDict();
 
    poString = PyUnicode_FromString("sys");
    poSys = PyDict_GetItem(poDictModule, poString);
    Py_DECREF(poString); // refcount = 0, it can be freed now
 
    if(poSys == NULL)
    {
        printf("Cannot find the sys item in the dictionary. Aborting...\n");
        return -1;
    }
 
    if(!PyModule_Check(poSys))
    {
        printf("This is not the sys module. Leaving...\n");
        return -1;
    }
 
    poDictModule = PyModule_GetDict(poSys);
    poList = PyDict_GetItemString(poDictModule, "path");
 
    if(!PyList_Check(poList))
    {
        printf("This is not the path list. Aborting...\n");
        return -1;
    }
 
    printf("Path list before:\n");
    for(i = 0; i < PyList_Size(poList); i++)
    {
        poListItem = PyList_GetItem(poList, i);
        if(PyUnicode_Check(poListItem))
        {
            char *value = PyBytes_AsString(PyUnicode_AsEncodedString(poListItem, "utf-8", "Error"));
            printf("\tValue [%d] = %s\n", i, value);
        }
    }
    printf("\n");
 
    poPath = PyUnicode_FromString(module_dir);
    PyList_Append(poList, poPath);
    Py_DECREF(poPath);
 
    printf("Path list after:\n");
    for(i = 0; i < PyList_Size(poList); i++)
    {
        poListItem = PyList_GetItem(poList, i);
        if(PyUnicode_Check(poListItem))
        {
            char *value = PyBytes_AsString(PyUnicode_AsEncodedString(poListItem, "utf-8", "Error"));
            printf("\tValue [%d] = %s\n", i, value);
        }
    }
 
    /* End Part 1 */
 
    /* Start Part 2 */
 
    // Here we will import our script. It is called example.py
    // We just pass the name of the module without the .py as the
    // Python interpreter knows that it should look for the file
    // example.py
    poModule = PyImport_ImportModule("example");
 
    if(poModule == NULL)
    {
        printf("Error importing module\n");
        PyErr_Print();
        return -1;
    }
 
    printf("Module imported succesfully!\n\n");
    poDictModule = PyModule_GetDict(poModule);
 
    size = 0;
    while(PyDict_Next(poDictModule, &size, &poKey, &poValue))
    {
        name = PyBytes_AsString(PyUnicode_AsEncodedString(poKey, "utf-8", "Error"));
        printf("Type: %s - Symbol name: %s \n", (char *)poValue->ob_type->tp_name, name);
    }
    printf("\n");
 
    // Get the objects of the symbols we defined in our script    
    poMyString = PyDict_GetItemString(poDictModule, "myString");
 
    if(poMyString)
    {
        if(PyUnicode_Check(poMyString))
        {
            str = PyBytes_AsString(PyUnicode_AsEncodedString(poMyString, "utf-8", "Error"));
            printf("myString value is: %s\n", str);
        }
        else
        {
            printf("myString is not a PyUnicode object. Are you messing up with the script?\n");
        }
    }
    else
    {
        printf("Symbol 'myString' not found...\n");
    }
 
    poDoubleValue = PyDict_GetItemString(poDictModule, "doubleValue");
 
    if(poDoubleValue)
    {
        if(PyFunction_Check(poDoubleValue))
        {
            poResult = PyObject_CallFunction(poDoubleValue, "i", 2);
            printf("doubleValue(2) = %ld\n", PyLong_AsLong(poResult));
            Py_DECREF(poResult);
        }
        else
        {
            printf("doubleValue is not a PyFunction object. Are you messing up with the script?\n");
        }
    }
    else
    {
        printf("Symbol 'doubleValue' not found...\n");
    }
 
    Py_Finalize();
 
    /* End Part 2 */
 
    return 0;
 
}

I will not repeat myself about things that have already been explained in the first post, so I will just give an overview on what is happening here.

This code will import a Python module (a script written by us), print the value of the variable 'myString' that is declared in the script and execute the function 'doubleValue' that (as the name suggests) takes one integer argument, multiplies it by 2 and returns the result.

In order to achieve this, we need some preparation.

First, for Python to be able to find our module, we need to first provide it with the path to where our script is. This is what 'Part One' in the code is all about.

When one wants to import a module, the Python interpreter will check all the paths defined in the list "sys.path" for the module to be imported. If the module is found, it is then loaded otherwise the interpreter will throw an error about the module not being found.

As our module will be in a place different than the default places that python looks for, we need to add this path to the list of module paths.

    poDictModule = PyImport_GetModuleDict();
    poString = PyUnicode_FromString("sys");
    poSys = PyDict_GetItem(poDictModule, poString);
    Py_DECREF(poString); // refcount = 0, it can be freed now

First, we grab the sys.modules dictionary so we can navigate later to the sys module and grab the path list.

We then create a PyUnicode representing the string "sys" that will be used to search for the "sys" module in the dictionary.

PyDict_GetItem() will two arguments: the dictionary we want to get an item from, and a PyObject with the string (a PyUnicode) representing the key in the dictionary. It then returns a PyObject* of the item associated with that key.

After we have our item, we decrement the reference count of 'poString' so it can be free()'d by the system. I will talk about refcounts in another post that is on its way, so don't worry too much about it for now.

At this point we have the module 'sys' in the poSys variable, and can finally go ahead to grab its dictionary and finding the 'path' list in it.

    poDictModule = PyModule_GetDict(poSys);
    poList = PyDict_GetItemString(poDictModule, "path");

This is the other way of getting an item from the dictionary. Instead of having all the trouble of creating a PyUnicode object, we can use PyDict_GetItemString() to get the item identified by the key passed as a char*.
Not much new stuff here, we grab the dictionary for the "sys" module and then grab the "path" list from it.
I know that "sys.path" is a list, because it is documented in the Python documentation somewhere (don't recall where thou).
So now have the sys.path list in the 'poList' variable.

    printf("Path list before:\n");
    for(i = 0; i < PyList_Size(poList); i++)
    {
        poListItem = PyList_GetItem(poList, i);
        if(PyUnicode_Check(poListItem))
        {
            char *value = PyBytes_AsString(PyUnicode_AsEncodedString(poListItem, "utf-8", "Error"));
            printf("\tValue [%d] = %s\n", i, value);
        }
    }
    printf("\n");

Here we are just printing all the items in the sys.path list. Prior to this point one can see that I have made all the checks to make sure that poList is really a PyList and that we can access it as such.

The sys.path list is a list of Strings, so it is okay to just go ahead and convert each item (that we retrieved using PyList_GetItem()) to a char* and then print it out.

An equivalent code in Python to do this would be:

import sys
 
print(sys.path)

Now we want to add the diretory where our module will be

    poPath = PyUnicode_FromString(module_dir);
    PyList_Append(poList, poPath);
    Py_DECREF(poPath);

As we are working with a list of PyUnicode (strings), we need to create a PyObject that represents our path as a PyUnicode, so again we are using PyUnicode_FromString() for this.

After that, we simply append the new path to the list, and again decremet the refcount of the recently created PyObject to make sure it will be cleaned up by the interpreter when it is time.

    printf("Path list after:\n");
    for(i = 0; i < PyList_Size(poList); i++)
    {
        poListItem = PyList_GetItem(poList, i);
        if(PyUnicode_Check(poListItem))
        {
            char *value = PyBytes_AsString(PyUnicode_AsEncodedString(poListItem, "utf-8", "Error"));
            printf("\tValue [%d] = %s\n", i, value);
        }
    }

Print the list again, to prove that our path is now in the sys.path list! :D

Great! Up to this point we saw some new functions and learned some cool stuff on how to handle a list and repracticed handling dictionaries.

What we want to do now is load our module, print all its dictionary entries, access its 'myString' member and call its 'doubleValue' function.
All this things are identified in the code as Part 2.

    // Here we will import our script. It is called example.py
    // We just pass the name of the module without the .py as the
    // Python interpreter knows that it should look for the file
    // example.py
    poModule = PyImport_ImportModule("example");

I guess the comment in this part says it all. Just to make sure, the script presented in the begining of this post should be saved in a file 'example.py' for this to work. If you have saved it in a file with a different name, just make the changes here :)

    printf("Module imported succesfully!\n\n");
    poDictModule = PyModule_GetDict(poModule);
 
    size = 0;
    while(PyDict_Next(poDictModule, &size, &poKey, &poValue))
    {
        name = PyBytes_AsString(PyUnicode_AsEncodedString(poKey, "utf-8", "Error"));
        printf("Type: %s - Symbol name: %s \n", (char *)poValue->ob_type->tp_name, name);
    }
    printf("\n");

The module was succesfully imported, we get its dictionary and start printing everything on it!

    poMyString = PyDict_GetItemString(poDictModule, "myString");
 
    if(poMyString)
    {
        if(PyUnicode_Check(poMyString))
        {
            str = PyBytes_AsString(PyUnicode_AsEncodedString(poMyString, "utf-8", "Error"));
            printf("myString value is: %s\n", str);
        }
        else
        {
            printf("myString is not a PyUnicode object. Are you messing up with the script?\n");
        }
    }
    else
    {
        printf("Symbol 'myString' not found...\n");
    }

We retrieve the object that refers to the 'myString' variable in our Python code. If the symbol 'myString' is not present in the dictionary, the function PyDict_GetItemString() will return NULL, so that is why we check to see if there is anything in the poMyString variable. After that, we just want to make sure it is a PyUnicode object, so we can print its value. The printing part is pretty much the same thing we have been doing for every other PyUnicode.

    poDoubleValue = PyDict_GetItemString(poDictModule, "doubleValue");
 
    if(poDoubleValue)
    {
        if(PyFunction_Check(poDoubleValue))
        {
            poResult = PyObject_CallFunction(poDoubleValue, "i", 2);
            printf("doubleValue(2) = %ld\n", PyLong_AsLong(poResult));
            Py_DECREF(poResult);
        }
        else
        {
            printf("doubleValue is not a PyFunction object. Are you messing up with the script?\n");
        }
    }
    else
    {
        printf("Symbol 'doubleValue' not found...\n");
    }

Now we retrieve the object that refers to the function 'myDouble' that we declared in the script. So far, this is pretty much the same thing we have been doing all along. After that we check to make sure that what we got is really a function, and then the interesting part!
PyObject_CallFunction() is one of the many ways to execute a function. Let's have a look at the prototype for PyObject_CallFunction() then:

PyObject* PyObject_CallFunction(PyObject *callable, char *format, ...)

Let's talk about the arguments first.
callable is the PyObject that represents the function we want to call.

format just like printf() is a function that receives a 'format' parameter, this parameter works the same way, the difference is that the characters used to represent thigs are different. This string will have the format of what kind of parameters will come right after. For our code, we have PyObject_CallFunction(poDoubleValue, "i", 2) - Here what happens is, the parameter right after the format is taken as an int, because the format string is "i". If the format string was "iii", we would need to pass 3 parameters, like this: PyObject_CallFunction(poDoubleValue, "iii", 2,3,4);

... all other parameters that are needed depending on the 'format' string.

The "format" parameter shall reflect as many parameters as the function is expecting.

After calling the function, simply print the result (which is an int, treated as a PyLong in Python).
We decrement the reference count for cleanup.

After all this, we just execute a PyFinalize() to finish the interpreter, and we are done!

References:

Python/C API Reference Manual

Filed under: C, Programming, Python 2 Comments
3May/111

Python C API – First Step

Posted by Henrique

I am once again in that time of the year where I start to look back at projects I started and did not finish. I am trying to make a personal commitment to select one or two of those and have them rolling. One of these projects involve embedding Python into it to make it extensible. As it has been a long time since I started this project, I barely remember how that stuff works, so I am once again tackling the problem to see where I can get.

I plan to make this post a little bit simple (but can't guarantee I will be able to), showing some code and explaining parts of it. Where I see the need (and have the knowledge to do so) I will try to discuss a little bit about how certain thing work.

For everything discussed here, I will be using Python 3.2. It is possible some of this stuff won't work with Python 2.x but will probably work fine with Python 3.1.

So, let's start with a simple program that will just initialize Python and list all modules that are loaded as soon as the initialization is finished.

#include "Python.h"
#include <stdio.h>
#include <stdlib.h>
 
int main(int argc, char** argv) {
    char *name;
    PyObject *dictModule;
    Py_ssize_t size;
    PyObject *poKey;
    PyObject *poValue;
 
    Py_Initialize();
    dictModule = PyImport_GetModuleDict();
 
    size = 0;
    while(PyDict_Next(dictModule, &size, &poKey, &poValue))
    {
        name = PyBytes_AsString(PyUnicode_AsEncodedString(poKey, "utf-8", "Error"));
        printf("Type: %s - Name: %s\n", (char *)poValue->ob_type->tp_name, name);
    }
 
    Py_Finalize();
 
    return 0;
}

To compile this code, you can use the following command (Considering you saved the source in a file called main.c):
gcc -Wall `python3.2-config --include` `python3.2-config --libs` -o pyx main.c

If python3.2-config is not available at you platform, this is what it expands to in my computer:
gcc -Wall -lpthread -ldl -lutil -lm -lpython3.2mu -I/usr/include/python3.2mu -I/usr/include/python3.2mu -o pyx main.c

It is up to you to figure out the include paths :) If you have trouble, let me know in the comments.

Since this is a new topic, I will talk a little bit about the new types we see: PyObject and Py_ssize_t.

Py_ssize_t is a typedef to ssize_t or to a long int. It is recommended to use Py_ssize_t as it might change in the future and this guarantees compatibility.

PyObject is an opaque type that seems to be the base to pretty much everything in Python. A module, a dictionary, an object, a function... they are all represented by a PyObject. There are functions and/or macros that can be used to identify what kind of information is stored in a PyObject. For example, if one wants to know if the PyObject in hand is a dictionary, the macro PyDict_Check(PyObject *p) can be used. This macro returns true if the PyObject 'p' is a dictionary.

I won't be discussing much of the inner workings of PyObject, because that would take a long time, and I am still not too familiar with that yet (Hopefully I will have a pot in the future just to discuss it). I guess it is enough to know that most functions that handle a python object will deal with this structure.

PyInitialize() simply initializes the interpreter and loads everything that python needs. I couldn't express this better than the documentation, so let me quote it:

Initialize the Python interpreter. In an application embedding Python, this should be called before using any other Python/C API functions; with the exception of Py_SetProgramName() and Py_SetPath(). This initializes the table of loaded modules (sys.modules), and creates the fundamental modules builtins, __main__ and sys. It also initializes the module search path (sys.path). It does not set sys.argv; use PySys_SetArgvEx() for that. This is a no-op when called for a second time (without calling Py_Finalize() first). There is no return value; it is a fatal error if the initialization fails.

It is safe to ignore the information regarding other functions for now.

So, as soon as Python has been initialized, I was curious about which modules were loaded to see if they would match what the documentation mentioned. Digging into it, I found the function PyImport_GetModuleDict() that will return the dictionary of the default loaded module sys.modules. With this in hand I am able to iterate through the dictionary and print everything that is available in it.

In order to navigate through the dictionary, we can use the function PyDict_Next() which has the following prototype:

int PyDict_Next(PyObject *p, Py_ssize_t *ppos, PyObject **pkey, PyObject **pvalue);

PyObject *p - The dictionary we want to iterate through.
Py_ssize_t *ppos - This will tell the offset inside the dictionary where the current data is found. This value HAS to be initialized to 0 prior to calling this function and should not be changed inside the loop calling PyDict_Next (as this is the variable that keeps track of the offset being used by the function).
PyObject **pkey - This will be a string object containing the name of the dictionary entry (the key of the dictionary, which represents the symbol being analyzed).
PyObject **pvalue - The value to which the current key refers to. This can be anything: a module, another dictionary, a function, a long, a string and so on...

Strings in Python are all Unicode and are represented by a PyObject as well. For that reason, in order to print the string value to the terminal, we need to first convert the PyUnicode object to a PyBytes object which can then be converted to a regular string (char *). This is what is happening in the following line:

name = PyBytes_AsString(PyUnicode_AsEncodedString(poKey, "utf-8", "Error"));

When trying to determine the type of the PyObject I could not find a function that would return a string telling the type, so I started analyzing the PyObject structure and found a field tp_name (which I presume means "type name") inside the ob_type structure contained in it. I guess accessing PyObject's members directly isn't recommended, but as this is being used for learning purposes I might be forgiven.

After everything is completed, we just finalize the Python interpreter by calling PyFinalize(). This will make sure Python does its cleanup and frees all memory it was using.

As stated in the post title, this is just the first step playing with the API. I plan to show some more code in the posts to follow, where we will be able to execute python code without interacting with the interpreter, and a code to import a Python module (a .py file with Python code) and then interact with the interpreter, by reading the data in the file and executing its functions.

For now, this is it!

References:

Python/C API Reference Manual

Filed under: C, Programming, Python 1 Comment
26Jan/110

Conjunto de Mandelbrot

Posted by Henrique

Aprendi sobre fractais quando era um pouco mais novo e conheci o conjunto de Mandelbrot naquela época.
Algum tempo atrás quis plotar o conjunto, porém não conseguia encontrar informações completas e/ou sequenciais que explicassem bem o assunto de uma maneira simples de entender. Sempre havia alguma coisa faltando.
Escrevi esse texto com o intuito de fornecer um raciocínio lógico de se seguir para entender como o conjunto de Mandelbrot é formado e como plota-lo.
Para aqueles que desejam conhecer mais sobre o conjunto, recomendo uma leitura no artigo da wikipedia sobre o assunto: Conjunto de Mandelbrot
O código que será apresentado mais afrente é escrito em linguagem C e utiliza a lib SDL.

Números Complexos

Uma pequena introdução (mais uma revisão na verdade) do que é um número complexo.
O número complexo é aquele que não pode ser representado no plano cartesiano real e é dividido em duas partes:

z = x + yi

x -> eixo X (Real)
y -> eixo Y (Imaginário)
i -> Número imaginário

i =  \sqrt{-1}


Sua representação é feita no plano complexo (também conhecido por plano de Argand):

Plano Complexo

Logo, quando se deseja plotar um número complexo utiliza-se a parte real como valor a ser marcado no eixo X, e a parte imaginária para ser marcada no eixo Y, representando assim um ponto no plano.

PS: Há muito mais sobre os números complexos do que apenas isso. Para conhecer mais afundo e entender melhor sobre o conjunto de mandelbrot, é bastante interessante estudar os números complexos... pode ser uma boa hora para desenterrar aqueles seus livros de matemática do ensino fundamental/médio.

O conjunto de Mandelbrot

O conjunto de Mandelbrot é definido pela seguinte equação:

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

Onde 'z' e 'c' são números imaginários do formato:

z = u + wi

u -> Parte Real (Coordenada do eixo X)
w -> Parte Imaginária (Coordenada do eixo Y)

c = a + bi

a -> Parte Real (Coordenada do eixo X)
b -> Parte Imaginária (Coordenada do eixo Y)

Como pode-se perceber pela equação, esse conjunto é definido através de uma fórmula recursiva, onde o resultado atual depende do resultado anterior (esse comportamento é observado também na série de Fibonacci).
A equação por si só não é suficiente para definir o conjunto. Pontos cuja distância da origem tenham valor maior que 2, não fazem parte do conjunto (Pode-se perceber (e provar) que quando um ponto atinge uma distância maior que 2, nas iterações que se seguem farão com que ele cresça rumo ao infinito).
Para os que não se lembram, a seguinte fórmula calcula a distância entre dois pontos:

d = \sqrt{(x_1 - x_0)^2 + (y_1 - y_0)^2}

Considerando a origem do plano no ponto O(0,0), a fórmula fica da seguinte maneira:

  d = \sqrt{(x_1 - 0)^2 + (y_1 - 0)^2}\\  d = \sqrt{x_1^2 + y_1^2}

De volta a equação que define a presença dos elementos do conjunto a expanão da recursividade da fórmula fica da seguinte maneira:

  Z_0 = 0 + 0i    => Z_0 = 0\\  Z_1 = Z_0^2 + c => Z_1 = c\\  Z_2 = Z_1^2 + c => Z_2 = c^2 + c\\  Z_3 = Z_2^2 + c => Z_3 = (c^2 + c)^2 + c\\  ...\\  Z_n+1 = (Z_n)^2 + c\\

Vamos voltar um pouco a teoria de números complexos então e fazer uma análise mais a fundo da equação que define o conjunto de Mandelbrot, para que possa ser aplicada em nosso algoritmo:

  z^2 = (u + wi)^2\\  z^2 = u^2 + 2*u*wi + wi^2\\  z^2 = u^2 + 2*u*wi + (w^2 * i^2)\\  z^2 = u^2 + 2*u*wi + (w^2 * -1)\\  z^2 = u^2 + 2*u*wi - w^2\\  z^2 = u^2 - w^2 + 2*u*wi\\

A parte real do número que representa o eixo X, é dado pela seguinte parte da equação:

u^2 - w^2

A parte imaginária do número que representa o eixo Y, é dado pela seguinte parte da equação:

2*u*wi

Portanto:

  	z_{n+1} = (z_n)^2 + c\\  	z_{n+1} = (u_n^2 - w_n^2 + 2*u_n*w_n*i) + (a + bi)

Em termos de código seria algo parecido com:

void mandelbrot(float u, float w, float a, float b)
{
	eixoX = u^2 - w^2 + a;
	eixoY = 2 * u * w + b;
}

Vamos supor alguns valores para as variáveis envolvidas, apenas para fim de demonstrar como ficaria a conta:
a = 2
b = 3
u = 5
w = 1
Substituindo na equação:

  	z_{n+1} = (5^2 - 1^2 + 2*5*1*i) + (2 + 3i)\\  	z_{n+1} = (25 - 1 + 10i) + 2 + 3i\\  	z_{n+1} = 24 + 10i + 2 + 3i\\  	z_{n+1} = 26 + 13i\\

Logo temos o ponto (26,13) no plano imaginário.

É claro que os valores dados acima são apenas um exemplo para demonstrar como é feito o calculo da equação e não se encontram no conjunto de Mandelbrot.

Eis então uma breve explicação de como funcionará o algoritmo para determinar se um número faz ou não parte do conjunto:

  • A primeira iteração tem valor 0 para 'z' e 'c' receberá o valor do ponto que desejamos testar.
  • Executa-se recursivamente o cálculo até um limite definido pelo programador
  • A recursão ocorre para cada ponto que desejamos testar para verificar se está ou não contido no conjunto de Mandelbrot.
  • O conjunto de pontos que podemos testar é infinito e os valores de x devem estar no intervalo [-2; 1] e os valores de y estejam no intervalo [-1; 1]. Não se esqueça, entre dois pontos existem infinitos números quando consideramos o conjunto dos números complexos.

Vamos trabalhar com um exemplo!
Suponha que desejamos testar se o ponto ((0,3), (0,4)) faz parte do conjunto. Teremos então:

  c = 0,3 + 0,4i\\  z_{n0} = 0 + 0,3 + 0,4i\\

Tem-se o ponto ((0,3), (0,4)). É necessário calcular a distância desse ponto à origem.

  d = \sqrt{0,3^2 + 0,4^2}\\  d = 0.5

A distância é menor que 2, logo deve-se realizar o cálculo da próxima iteração

  z_{n1} = z_{n0}^2 + c\\  z_{n1} = (0,3 + 0,4i)^2 + (0,3 + 0,4i)\\  z_{n1} = (0,3^2 + 2*0,3*0,4i + 0,4i^2) + (0,3 + 0,4i)\\  z_{n1} = (0,09 + 0,24i - 0,16i) + (0,3 + 0,4i)\\  z_{n1} = 0,09 + 0,24i - 0,16i + 0,3 + 0,4i\\  z_{n1} = 0,39 + 0,48i\\

Tem-se o ponto ((0,39), (0,48)). É necessário calcular a distância desse ponto à origem.

  d = \sqrt{0,39^2 + 0,48^2}\\  d = 0.6184667

Novamente a distância é menor que 2, logo deve-se realizar o cálculo da próxima iteração.

  z_{n2} = z_{n1}^2 + c\\  z_{n2} = (0,39 + 0,48i)^2 + (0,3 + 0,4i)\\  ...

Esse processo se repetirá por n iterações. No caso do código abaixo, isso se repetirá 1024 vezes ou até que a distância do ponto à origem seja maior que 2.
Caso as 1024 iterações ocorram sem que o resultado de nenhuma delas extrapole a distância de 2 da origem, então marcamos o ponto 'c' como fazendo parte do conjunto.
Por exemplo, quando iteramos o ponto p(-0.75, -0.01562) ao chegarmos na iteração 201, o valor estoura para fora do limite permitido pelo conjunto, atingindo o ponto q(1.11917, 3.12583) que está a uma distância maior do que 2 da origem. Logo, o ponto p(-0.75, -0.01562) não está no conjunto de Mandelbrot.

Implementando

A seguir então uma implementação do algoritmo para plotar o gráfico referente ao conjunto de Mandelbrot.
O código está escrito em C e utiliza a lib SDL para os gráficos.
Boa parte do código está comentada e acredito que sejam suficientes para explicar seu funcionamento.
O código não está otimizado, então ele demora um pouco a executar (no meu computador leva coisa de 5 segundos, sendo que o processador é um Core 2 Duo E7400 de 2.8ghz).
Sinta-se a vontade para otimiza-lo e adicionar cor ao gráfico.

 
/**
 * Salve o código em um arquivo chamado mandelbrot.c
 * Compile com: gcc -Wall -lm -lSDL -o mandelbrot mandelbrot.c
 * Execute com: ./mandelbrot
 * Use o código e altere a vontade. Se precisar de alguma ajuda
 * ou tiver alguma duvida, fique a vontade para me contactar.
 * Se for republicar o código em algum lugar, agradeço se mantiver
 * os créditos ou um link para o código original.
 *
 * Autor: Henrique M. D.
 * Email: typoon at gmail dot com
 */
 
#include <stdio.h>
#include <math.h>
#include <SDL/SDL.h>
 
#define WIDTH  800 /* Largura da tela */
#define HEIGHT 800 /* Altura da tela */
#define MAX_ITER 1024
 
void setPixel(SDL_Surface *screen, int x, int y, SDL_Color *cor);
void scaleXY(int pX, int pY, float *cX, float *cY);
int mandelbrot(float zX, float zY, float cX, float cY, int iter);
float distOrigem(float x, float y);
void putpixel(SDL_Surface *surface, int x, int y, Uint32 pixel);
 
int main(int argc, char **argv)
{
 
    int i, j;
    int iter;
    int exit;
    float cX, cY;
 
    SDL_Surface *screen;
    SDL_Color black;
    SDL_Color white;
    SDL_Event event;
 
    if(SDL_Init(SDL_INIT_VIDEO))
    {
        printf("Error initializing SDL [%s]\n", SDL_GetError());
        return -1;
    }
 
    atexit(SDL_Quit);
    screen = SDL_SetVideoMode(WIDTH, HEIGHT, 32, SDL_HWSURFACE);
 
    if(screen == NULL)
    {
        printf("Error setting video mode [%s]\n", SDL_GetError());
        return -1;
    }
 
    black.r = 0;
    black.g = 0;
    black.b = 0;
 
    white.r = 255;
    white.g = 255;
    white.b = 255;
 
    for(i = 0; i < WIDTH; i++)
    {
        for(j = 0; j < HEIGHT; j++)
        {
            scaleXY(i, j, &cX, &cY); /* Converter o pixel em coordenadas */
            iter = mandelbrot(0.0, 0.0, cX, cY, 0);
            if(iter == MAX_ITER) /* ponto do conjunto, pintar de preto */
            {
                setPixel(screen, i, j, &black);
            }
            else
            {
                setPixel(screen, i, j, &white);
            }
        }
    }
 
    SDL_Flip(screen);
 
    exit = 0;
    while(!exit)
    {
        SDL_PollEvent(&event);
        switch(event.type)
        {
            case SDL_QUIT:
                exit = 1;
            break;
        }
    }
 
    return 0;
 
}
 
/**
 * Função praticamente copiada de:
 * http://www.cpp-home.com/tutorials/154_1.htm
 * Com leves alterações
 * Serve para setar o pixel na superficie com a cor especificada
 *
 * @param screen Superficie onde pintar o pixel
 * @param x Posicao x do pixel
 * @param y Posicao y do pixel
 * @param cor A cor que se deseja pintar o pixel
 *
 */
void setPixel(SDL_Surface *screen, int x, int y, SDL_Color *cor)
{
    int bpp;
    Uint8 *p;
    Uint32 pixel;
 
    pixel = SDL_MapRGB(screen->format, cor->r, cor->g, cor->b);
 
    bpp = screen->format->BytesPerPixel;
    p = (Uint8 *)screen->pixels + y * screen->pitch + x * bpp;
 
    SDL_LockSurface(screen);
 
    switch(bpp)
    {
        case 1:
            *p = pixel;
        break;
 
        case 2:
            *(Uint16 *)p = pixel;
        break;
 
        case 3:
        if(SDL_BYTEORDER == SDL_BIG_ENDIAN)
        {
            p[0] = (pixel >> 16) & 0xff;
            p[1] = (pixel >> 8) & 0xff;
            p[2] = pixel & 0xff;
        } else {
            p[0] = pixel & 0xff;
            p[1] = (pixel >> 8) & 0xff;
            p[2] = (pixel >> 16) & 0xff;
        }
        break;
 
        case 4:
            *(Uint32 *)p = pixel;
        break;
    }
 
    SDL_UnlockSurface(screen);
 
}
 
 
/**
 * Faz a escala que converte o valor do pixel da tela
 * para o equivalente do plano complexo para o conjunto de Mandelbrot
 * que encontra-se em:
 * x = [-2, 1]
 * y = [-1, 1]
 *
 * @param pX - Pixel no eixo X
 * @param pY - Pixel no eixo Y
 * @param cX - Valor em que o pixel é representado no plano complexo entre -2 e 1
 * @param cY - Valor em que o pixel é representado no plano complexo entre -1 e 1
 *
 */ 
void scaleXY(int pX, int pY, float *cX, float *cY)
{
 
    int coordX0;
    int coordY0;
    float scaleX;
    float scaleY;
 
    /* 2.5 para esquerda do eixo */
    /* 1.0 para direita do eixo */
    coordX0 = (WIDTH / 3.0) * 2.0;
 
    /* 1.0 para cima do eixo  */
    /* 1.0 para baixo do eixo */
    coordY0 = (HEIGHT / 2) * 1;
 
    scaleX = 3.0 / WIDTH;
    scaleY = 2.0 / HEIGHT;
 
 
    *cX = -2.0 + (pX * scaleX);
 
    if(pY > coordY0) /* Coordenada Y negativa */
    {
        *cY = -1 + ((HEIGHT - pY) * scaleY);
    }
    else
    {
        *cY = ((pY - (HEIGHT/2)) * scaleY);
    }
 
}
 
/**
 * Essa função verifica se o ponto que lhe foi passado pertence
 * ao conjunto de Mandelbrot para MAX_ITER iteraçãoes.
 * Caso não pertença, ela retorna o valor da iteração na qual o ponto
 * extrapolou o range permitido do conjunto, caso contrário, o valor
 * retornado será igual a MAX_ITER
 *
 * @param zX - Valor do eixo real X para o numero complexo z
 * @param zY - Valor do eixo imaginario Y para o numero complexo z
 * @param cX - Valor do eixo real X para o numero complexo c
 * @param cY - Valor do eixo imaginário Y para o numero complexo c
 * @param iter - Número da iteração na qual estamos (É o valor de n em Z_n)
 * @return A quantidade de iterações recursivas que foram executadas
 *
 */
int mandelbrot(float zX, float zY, float cX, float cY, int iter)
{
 
    float dX;
    float dY;
 
    dX = ((zX*zX) - (zY*zY)) + cX;
    dY = (2*zX*zY) + cY;
 
    if(iter >= MAX_ITER) /* Ponto passou por todas iterações e está no conjunto */
        return iter;
 
    if(distOrigem(dX, dY) <= 2) /* Ponto está no conjunto */
        return mandelbrot(dX, dY, cX, cY, ++iter);
 
    return iter;
}
 
/**
 * Calcula a distancia do ponto passado a origem O(0,0)
 * do plano
 *
 * @param x Valor de X do ponto
 * @param y Valor de Y do ponto
 * @return A distância do ponto a origem
 */
float distOrigem(float x, float y)
{
    return sqrtf((x*x) + (y*y));
}

O resultado desse programa é uma janela com a seguinte imagem
mandelbrot

Filed under: C, Matemática, SDL No Comments
19Jan/110

Project Euler – Problema 1

Posted by Henrique

Ao perguntar a um amigo no IRC sobre uma recomendação de hobby, ele me respondeu 'Matemática!'. Eu tinha em mente um hobby onde não precisasse pensar muito, mas matemática já havia me passado pela mente diversas vezes.
Acabei por arranjar um hobby onde não é preciso pensar muito e resolvi também adicionar a matemática como um hobby. Esse post é para falar sobre este último :)
Faz algum tempo fui apresentado ao 'Project Euler', um site com desafios matemáticos que podem ser resolvidos com ou sem o auxílio de um computador. Decidi que vou tentar resolver os problemas, inicialmente sem o auxílio de um computador (quando possível), e então implementarei uma solução utilizando programação.

O Problema

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.
Find the sum of all the multiples of 3 or 5 below 1000.

Em português:

Se nós listarmos todos os números naturais abaixo de 10 que sejam múltiplos de 3 ou 5, teremos 3, 5, 6 e 9. A soma desses múltiplos é 23.
Encontre a soma de todos os multiplos de 3 ou 5 abaixo de 1000.

A Solução Matemática

A primeira idéia para resolver esse problema é dar uma olhada em parte dos conjuntos formados pelos números menores que 1000 e múltiplos de 3 e de 5

M3 = \{3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36, ..., 996, 999\}

M5 = \{5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, ..., 990, 995\}

M3 terá um total de 333 elementos (999 / 3)
M5 terá um total de 199 elementos (999 / 5)

Algo para se notar é que em ambos os conjuntos existem números repetidos, como por exemplo {15, 30, 45, 60, ...}, logo se realizarmos a soma dos dois conjuntos teremos duplicidade desses números nos dando a resposta incorreta.
Para eliminar esses números de um dos dois conjuntos, precisamos simplesmente subtrair a soma de todos eles do resultado final e para isso precisamos descobrir quais são todos os números múltiplos de 3 E 5 e aplicar o mesmo procedimento de somatória.
Como descobrir esses números então? Simples! Calcula-se o M.M.C. (mínimo múltiplo comum) de 3 e 5, e encontra-se o conjunto dos múltiplos desse valor encontrado para todos os valores abaixo de 1000.

MMC(3,5) = 15

Temos então o novo conjunto M15:

M15 = \{15, 30, 45, 60, 75, 90, ..., 975, 990\}

M15 terá um total de 66 elementos (999 / 15)

M15 então é um subconjunto de M3 e também um subconjunto de M5. (Pode-se dizer que M15 está contido em M3 e em M5).

O problema propõe a soma de todos os números múltiplos de 3 OU 5, então precisamos remover o subconjunto M15 de um dos dois conjuntos (M3 ou M5).
Por quê?

Imagine que o problema pedisse a soma de todos os múltiplos de 3 ou 5 abaixo de 20. Teriamos os seguintes conjuntos:

M3 = \{3, 6, 9, 12, {\bf15}, 18\}

M5 = \{5, 10, {\bf15}\}

M15 = \{15\}

Se realizarmos a somatória de todos os elementos de M3 e M5 a resposta estará incorreta, pois o valor 15 foi somado duas vezes a resposta final, uma vez na somatória dos elementos de M3 e uma vez na somatória de M5.
A solução então é somar M3 a M5 e subtrair M15, fazendo com que o subconjunto M15 seja removido de um dos dois conjuntos.
Então a equação para resolução do problema é:

\sum M3 + \sum M5 - \sum M15

Não é muito prático realizarmos essa soma na mão, já que ela é composta de muitos números. Olhando bem para as fórmulas, podemos descrever os conjuntos da seguinte maneira:

  • O conjunto M3 é a somatória dos primeiros 333 números múltiplos de 3.
  • O conjunto M5 é a somatória dos primeiros 199 números múltiplos de 5.
  • O conjunto M15 é a somatória dos primeiros 66 números múltiplos de 15.

Em termos matemáticos:

\sum_{i = 1}^{333}3*i

\sum_{i = 1}^{199}5*i

\sum_{i = 1}^{66}15*i

Se expandirmos a primeira equação, teremos algo como:

(3*1) + (3*2) + (3*3) + (3*4) + ... + (3*332) + (3*333)

Todos os elementos estão sendo multiplicados por 3, logo, pode-se colocar o 3 em evidência e teremos:

3*(1+2+3+4+...+332+333)

Que pode ser representado por:

3*\sum_{i = 1}^{333}i

Esse é um resultado bastante interessante pois foi postulado por Gauss (e várias demonstrações do porque podem ser encontradas aqui) que a somatória dos n primeiros números pode ser escrita da seguinte maneira:

\sum_{i = 1}^n = n*(n+1)/2

Logo, temos a seguinte equação:

(3*\sum_{i = 1}^{333}i + 5*\sum_{i = 1}^{199}) - 15*\sum_{i = 1}^{66}

Que pode ser escrita da seguinte maneira:

(3*(333*(333+1)/2) + 5*(199*(199+1)/2)) - 15*(66*(66+1)/2)

Isso resulta em: 233168

E está ai então a resolução matemática.

A solução com programação

Solução onde não se usam os conceitos matemáticos apresentados acima:

 
#include <stdio.h>
 
int main(int argc, char **argv)
{
    int i;
    int resultado = 0;
    for(i = 1; i < 1000; i++)
    {
        if((i%3 == 0) || (i%5 == 0))
            resultado += i;
    }
 
    printf("Resultado: %d\n", resultado);
 
    return 0;
}

Solução onde os conceitos matemáticos acima são utilizados:

 
#include <stdio.h>
 
int gauss(int n)
{
    return (n*(n+1))/2;
}
 
int main(int argc, char **argv)
{
    int resultado;
 
    resultado = (3*gauss(333) + 5*gauss(199)) - 15*gauss(66);
    printf("Resultado: %d\n", resultado);
 
    return 0;
}
19Sep/092

Radix Sort

Posted by Henrique

Estive lendo sobre alguns algoritmos de organização e me deparei com o Radix Sort.
Pela definição da NIST:


Definition: A multiple pass distribution sort algorithm that distributes each item to a bucket according to part of the item's key beginning with the least significant part of the key. After each pass, items are collected from the buckets, keeping the items in order, then redistributed according to the next most significant part of the key.

Tradução livre:
Definição: Um algorítmo de organização de vários passos que distribui cada item para um bucket de acordo com parte da chave do item, começando pela parte menos significativa da chave. Após cada passo, os items são recolhidos dos buckets, mantendo os itens em ordem e então redistribuindo de acordo com a próxima parte mais significativa da chave.

Pelas leituras que fiz por aí, um método que me pareceu bastante interessante é aplicar esse algorítmo a cada bit da nossa informação. É claro que isso vai depender da aplicação do algorítmo. No caso que irei trabalhar aqui, estaremos ordenando um vetor com números inteiros positivos.

Vamos supor o seguinte array:

int vetor = {12,22,13,16,31}

Descrevendo cada um desses valores em binário:
12 - 01100
22 - 10110
13 - 01101
16 - 10000
31 - 11111

Pela descrição do algorítmo e pelo método adotado, devemos começar organizando nossa informação pelo bit menos significativo de cada um dos números em nosso array, ou seja, pegamos o bit mais a direita do array e distribuímos os valores em 2 buckets. O bucket onde o bit que estamos olhando é 0 ou ou bucket onde o bit que estamos olhando é 1. Ilustrando isso, temos:

01100 - 12
10110 - 22
01101 - 13
10000 - 16
11111 - 31

Começamos olhando para os bits que estão em negrito. Quando encontrarmos um bit 0, jogamos aquele número para o bucket0. Quando encontrarmos um bit 1, jogamos aquele número para o bucket1. Nossos buckets estão dessa forma agora então:

bucket0 = {12, 22, 16}
bucket1 = {13, 31}

Agora nós devolvemos os valores para dentro do nosso vetor na ordem em que se encontram dentro dos bukets. Primeiro colocamos dentro do vetor todos os valores do bucket0 e depois todos os valores do bucket1. Nosso vetor ficará da seguinte maneira:

vetor = {12, 22, 16, 13, 31}

Agora nós repetimos o processo para o nosso novo vetor, utilizando o próximo bit mais significativo.

01100 - 12
10110 - 22
10000 - 16
01101 - 13
11111 - 31

Nossos buckets ficam:

bucket0 = {12, 16, 13}
bucket1 = {22, 31}

Nosso novo vetor:

vetor = {12, 16, 13, 22, 31}

Repetimos o processo para o próximo bit mais significativo em nosso novo vetor:

01100 - 12
10000 - 16
01101 - 13
11010 - 22
11111 - 31

Nossos buckets ficam:

bucket0 = {16, 22}
bucket1 = {12, 13, 31}

Nosso novo vetor:

vetor = {16, 22, 12, 13, 31}

Repetimos o processo para o próximo bit mais significativo em nosso novo vetor:

10000 - 16
11010 - 22
01100 - 12
01101 - 13
11111 - 31

Nossos buckets ficam:

bucket0 = {16}
bucket1 = {22, 12, 13, 31}

Nosso novo vetor:

vetor = {16,22,12,13,31}

Até então não parece nada arrumado não é mesmo? Mas o útimo passo vai nos revelar algo fantástico! :D

Repetimos o processo para o próximo bit mais significativo em nosso novo vetor:

10000 - 16
11010 - 22
01100 - 12
01101 - 13
11111 - 31

Nossos buckets ficam:

bucket0 = {12, 13}
bucket1 = {16, 22, 31}

Nosso novo (e último) vetor:

vetor = {12, 13, 16, 22, 31}

Aha! Fantástico!

O mais interessante disso tudo é a velocidade desse algorítmo. Em alguns testes que fiz aqui, consegui organizar 20 milhoes de números em apenas 9,3 segundos!
Irei demonstrar a seguir uma versão um pouco mais lenta do algorítmo que desenvolvi porém de mais fácil entendimento. Logo após irei demonstrar a versão do mesmo algorítmo que foi capaz de organizar 100 milhões de números em 15,56 segundos.

O algorítmo 'simples'

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
 
typedef struct lista {
    int num;
    struct lista *prox;
} lista;
 
void radix_sort(int *, int);
void insere(lista **, int);
int pop(lista **);
void imprime(int *vetor, int tamanho);
 
int main()
{
    int tamanho = 20;
    int vetor[tamanho];
    int i;
 
    srand(time(NULL));
    for(i = 0; i < tamanho; i++) {
        vetor[i] = rand() % tamanho;
    }
 
    imprime(vetor, tamanho);
    printf("\n");
 
    radix_sort(vetor, tamanho);
 
    imprime(vetor, tamanho);
    printf("\n");
 
    return 0;
}
 
void insere(lista **l, int valor)
{
    lista *aux = *l;
 
    if(*l == NULL) {
        *l = (lista *)malloc(sizeof(lista));
        (*l)->num = valor;
        (*l)->prox = NULL;
    } else {
        while(aux->prox != NULL) aux = aux->prox;
        aux->prox = (lista *)malloc(sizeof(lista));
        aux = aux->prox;
        aux->num = valor;
        aux->prox = NULL;
    }
}
 
int pop(lista **l)
{
    int ret;
    lista *aux;
 
    ret = (*l)->num;
    aux = *l;
    *l = (*l)->prox;
    free(aux);
 
    return ret;
}
 
void radix_sort(int *vetor, int tamanho)
{
 
    lista *bucket0 = NULL;
    lista *bucket1 = NULL;
    int i, j;
    int numero_bits;
 
    numero_bits = sizeof(int) * 8;
 
    for(i = 0; i < numero_bits; i++) {
        for(j = 0; j < tamanho; j++) {
            if((((vetor[j] >> i) & 0xf) % 2) == 0)
                insere(&bucket0, vetor[j]);
            else
                insere(&bucket1, vetor[j]);
        }
 
        j = 0;
        while(bucket0 != NULL) {
            vetor[j] = pop(&bucket0);
            j++;
        }
 
 
        while(bucket1 != NULL) {
            vetor[j] = pop(&bucket1);
            j++;
        }
 
    }
 
}
 
void imprime(int *vetor, int tamanho)
{
    int i;
 
    for(i = 0; i < tamanho; i++)
        printf("%d ", vetor[i]);
}

[radix.c]

Essa é a versão 'simples' do algorítmo. Ele utiliza listas para ficar mais simples a inserção e remoção dos registros de nossos buckets. O problema com isso é que conforme aumentamos a quantidade de elementos que queremos organizar o programa fica mais lento. A verdadeira função de ordenação aqui é a função 'radix_sort'. As outras funções fazem parte da implementação da lista.
Vamos então dar uma olhada no programa.

...
    int tamanho = 20;
    int vetor[tamanho];
    int i;
 
    srand(time(NULL));
    for(i = 0; i < tamanho; i++) {
        vetor[i] = rand() % tamanho;
    }
...

Simplesmente declaramos um array com 20 posições e colocamos números randomicos nele que vão de 0 a 19. Esse é o vetor que desejamos ordenar.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
void radix_sort(int *vetor, int tamanho)
{
 
    lista *bucket0 = NULL;
    lista *bucket1 = NULL;
    int i, j;
    int numero_bits;
 
    numero_bits = sizeof(int) * 8;
 
    for(i = 0; i < numero_bits; i++) { // Primeiro loop
        for(j = 0; j < tamanho; j++) { // Segundo loop
            if((((vetor[j] >> i) & 0xf) % 2) == 0)
                insere(&bucket0, vetor[j]);
            else
                insere(&bucket1, vetor[j]);
        }
 
        j = 0;
        while(bucket0 != NULL) {
            vetor[j] = pop(&bucket0);
            j++;
        }
 
 
        while(bucket1 != NULL) {
            vetor[j] = pop(&bucket1);
            j++;
        }
    }
}

Como na explicação sobre o algorítmo, estou usando aqui 2 buckets e o vetor com os dados. O bucket0 recebe os números cujo bit atual sendo verificado é 0, e o bucket 1 recebe os números cujo bit atual sendo verificado é 1.
A variável 'numero_bits' se refere a quantos bits nós iremos verificar para realizar a ordenação. No caso como estamos organizando números inteiros, em teoria (mais tarde discutirei uma maneira melhor de determinar esse número) faz sentido percorremos todos os bits do mesmo. sizeof(int) retorna o número de bytes de um int, se multiplicamos esse valor por 8 temos o número de bits.
O primeiro loop se refere a qual bit estamos olhando e o segundo loop passa por cada elemento para determinar se o bit é 0 ou 1 e jogar o número em seu respectivo bucket.
Acredito que a linha que mais chama a atenção aqui é a linha 13. Vamos a explicação então! Imagine que estou trabalhando com um número de 4 bits apenas que se encontra em vetor[0]. O que essa linha faz é pegar o número inteiro e fazer uma lógica AND com o valor 0xf (eu poderia ter usado 0x1). Suponhamos que o numero no nosso vetor é 4. Então a lógica fica assim:

Deslocar o valor dentro de vetor[0] para a direita em 0 bits (primeira interação).
Em hexa: 0x4 & 0xf
Em binario: 0100 & 1111
Resultado: 0100

Deslocar o valor dentro de vetor[0] para a direita em 1 bit (segunda interação).
Em hexa: 0x2 & 0xf
Em binario: 0010 & 1111
Resultado: 0010

Deslocar o valor dentro de vetor[0] para a direita em 2 bits (terceira interação).
Em hexa: 0x1 & 0xf
Em binario: 0001 & 1111
Resultado: 0001

Deslocar o valor dentro de vetor[0] para a direita em 3 bits (quarta interação).
Em hexa: 0x0 & 0xf
Em binario: 0000 & 1111
Resultado: 0000

O que nos interessa aqui é o bit menos significativo, ou seja, o bit mais a direita. Quando realizamos essa lógica AND, o último bit será 0 ou 1. Se o último bit for 0, temos um número par e então o número deve ser colocado no bucket0. Se o último bit for 1, temos um número ímpar e então o número deve ser colocado no bucket1. O deslocamento que fazemos para a direita é apenas para colocar o bit que desejamos verificar na posição menos significativa do número para que possamos saber se ele deve ir ao bucket0 ou ao bucket1.
A lógica que segue é apenas referente a 'remontar' o vetor com todos os valores que estão no bucket0 e depois todos os valores que estão no bucket1.

O algorítmo um pouco menos simples

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
 
void radix_sort(int *, int, int);
void imprime(int *vetor, int tamanho);
 
int main()
{
    int tamanho = 20;
    int *vetor;
    int i;
    int maior = 0;
    clock_t inicio, fim;
    float periodo;
 
    vetor = (int *)malloc(sizeof(int)*tamanho);
 
    srand(time(NULL));
    for(i = 0; i < tamanho; i++) {
        vetor[i] = (rand() % tamanho)+1;
        if(vetor[i] > maior)
            maior = vetor[i];
    }
 
    imprime(vetor, tamanho);
    printf("\n");
 
    inicio = clock();
    radix_sort(vetor, tamanho, maior);
    fim = clock();
    periodo = (float)(((float)fim-(float)inicio)/CLOCKS_PER_SEC);
 
    printf("Tempo: %.15f\n", periodo);
    imprime(vetor, tamanho);
    printf("\n");
 
    return 0;
}
 
void radix_sort(int *vetor, int tamanho, int maior)
{
 
    int i, j, k, l;
    int numero_bits;
 
    int **temp;
 
    temp = (int **)malloc(sizeof(int *)*2);
    temp[0] = (int *)malloc(sizeof(int)*tamanho);
    temp[1] = (int *)malloc(sizeof(int)*tamanho);
 
    numero_bits = (log(maior)/log(2))+1;
 
    for(i = 0; i < numero_bits; i++) {
        k = 0;
        l = 0;
 
        for(j = 0; j < tamanho; j++) {
            if((((vetor[j] >> i) & 0xf) % 2) == 0) {
                temp[0][k++] = vetor[j];
            } else {
                temp[1][l++] = vetor[j];
            }
        }
 
        temp[0][k] = '0';
        temp[1][l] = '0';
 
        j = 0;
        while(temp[0][j] != '0') {
            vetor[j] = temp[0][j];
            j++;
        }
 
        k=0;
        while(temp[1][k] != '0') {
            vetor[j] = temp[1][k];
            j++;
            k++;
        }
 
    }
}
 
void imprime(int *vetor, int tamanho)
{
    int i;
 
    for(i = 0; i < tamanho; i++)
        printf("%d ", vetor[i]);
}

[radix2.c]

Considerações:

  • Essa implementação do algorítmo é apenas para números maiores que 0
  • Para compilar com o gcc nçao esqueça de compilar o código com -lm para a biblioteca math

As chamadas feitas para a função clock() antes e depois da chamada de radix_sort() são para que possamos visualizar quanto tempo a execução do algorítmo levou.
A primeira pequena diferença que vemos nesse código é na linha 12, onde agora declaramos nossa variável vetor como um ponteiro para que possamos alocar o espaço necessário para ela em runtime. Isso se faz necessário pois quando executarmos o programa para organizar 20 milhões de registros, não seria possível alocar o espaço necessário na stack.
A segunda pequena diferença que vemos nesse código está na linha 22, onde garantimos que o número sorteado é no mínimo 1.
A terceira pequena diferença é que agora nós já identificamos qual é o maior número do nosso vetor. Isso será útil mais a frente.
As grandes diferenças começam dentro da função radix_sort().
A primeira coisa para se notar é que não usamos mais listas. Inserir e retirar items de um array é mais rápido, logo é isso que usaremos. Nosso bucket agora é a variável temp; temp[0] é o bucket0 e temp[1] é o bucket1.
A segunda coisa a se notar é a forma de como calculamos o numero de bits. Essa linha é de extrema imoprtância pois ela reflete diretamente na perfomance do algorítmo. Suponha que estamos trabalhando em uma plataforma onde a variável int tem 32 bits. Suponha também que o maior número que se encontra em nosso array é o número 230.
Vamos representar então o número 230 em binário:
00000000000000000000000011100110

São 32 bits, apenas 8 bits significativos e 24 bits não significativos. A partir do momento que organizamos os números usando todos os bits significativos é inútil que façamos a verificação para os dígitos não significativos pois eles irão sempre ser colocados no bucket0 e sempre na mesma ordem.
Para entender o cálculo utilizado, precisamos lembrar um pouco da matemática da escola.
2^x = 16
Quanto vale x?
log(base 2) de 16 = x
x = 4

Ou seja, para representarmos 16 números precisamos de 4 bits. Sabemos que o a contagem para os números binários começa em 0, logo com esses 4 bits podemos contar de 0 a 15. Para que o número 16 fosse considerado precisamos adicionar mais 1 bit.
Logo, para representar o número 16 temos: (log(base 2) de 16)+1
Como em C não temos uma forma de calcular log em base 2. nós usamos uma propriedade matemática que diz:

log em base x de y = (log y) / (log x)
Para os que não se lembram log sem menção de base é log em base 10.

Voltando então ao nosso exemplo do número 230. Se calcularmos:
log(base 2) de 230 = x
x = 7,85

Não há possibilidade de termos 0.85 bit. Em C variáveis int se arredondam pra baixo. Logo em nosso caso x vale 7.
2^7 = 128
Não é o suficiente, e matemáticamente falando 7,85 seria arredondado de fato para 8. É por essa razão que somamos 1.
2^8 = 256
Agora sim temos espaço o bastante para o valor 230.

Caso não tenha ficado muito clara a explicação, tente ler novamente e se precisar relembre um pouco de logarítmos. Tenho certeza que ficará mais claro.

O resto da lógica que segue é basicamente o mesmo. Como precisamos de uma forma de identificar quando os números de nosso bucket acabaram, eu optei por dizer que o último número do bucket é 0. É por essa razão que o esse algorítmo só é capaz de organizar números maiores que 0.
Todo o resto é basicamente o mesmo já apresentado antes.

Observações
Se você testou esses algorítmos deve ter percebido que são bastante rápidos. Eu quis determinar o quão rápida essa segunda implementação apresentada é quando comparada ao quicksort e os resultados são impressionantes.

Eis então os algorítmos usados:

Quicksort

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
 
void quicksort(int left, int right, int *v);
void swap(int i, int j, int **v);
 
int main(int argc, char **argv)
{
 
	int MAX;
 
	if(argc < 2) {
		printf("Use: %s numero\n", argv[0]);
		return -1;
	}
 
	MAX= atoi(argv[1]);
 
	int i;
	int *vetor;
    clock_t inicio, fim;
    float periodo;
 
	vetor = (int *)malloc(sizeof(int)*MAX);
 
	srand(time(NULL));
	for(i = 0; i < MAX; i++) {
		vetor[i] = rand() % MAX;
	}
 
	inicio = clock();	
	quicksort(0, MAX-1, vetor);
	fim = clock();
 
    periodo = (float)(((float)fim-(float)inicio)/CLOCKS_PER_SEC);
    printf("Tempo: %.15f\n", periodo);
 
	return 0;
}
 
void quicksort(int left, int right, int *v)
{
	int pivot = v[(left+right)/2];
	int i = left;
	int j = right;
 
	if ( left < right ) {
 
		do {
 
			while (v[i] < pivot) i++;
			while (v[j] > pivot) j--;
 
			if( i <= j ) {
				swap(i, j, &v);
				i++;
				j--;
			}
 
 
		} while( i <= j );
 
		quicksort(left, j, v);
		quicksort(i, right, v);
	}
}
 
void swap(int i, int j, int **v) {
	int aux;
 
	aux = (*v)[j];
	(*v)[j] = (*v)[i];
	(*v)[i] = aux;
 
}

Radixsort

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
 
void radix_sort(int *, int, int);
void imprime(int *vetor, int tamanho);
 
int main(int argc, char **argv)
{
    int tamanho;
    int *vetor;
    int i;
    int maior = 0;
    clock_t inicio, fim;
    float periodo;
 
	if(argc < 2) {
		printf("Use: %s numero\n", argv[0]);
		return -1;
	}
 
	tamanho = atoi(argv[1]);
 
    vetor = (int *)malloc(sizeof(int)*tamanho);
 
    srand(time(NULL));
    for(i = 0; i < tamanho; i++) {
        vetor[i] = (rand() % tamanho)+1;
        if(vetor[i] > maior)
            maior = vetor[i];
    }
 
    //imprime(vetor, tamanho);
    //printf("\n");
 
    inicio = clock();
    radix_sort(vetor, tamanho, maior);
    fim = clock();
    periodo = (float)(((float)fim-(float)inicio)/CLOCKS_PER_SEC);
 
    printf("Tempo: %.15f\n", periodo);
    //imprime(vetor, tamanho);
    //printf("\n");
 
    return 0;
}
 
void radix_sort(int *vetor, int tamanho, int maior)
{
 
    int i, j, k, l;
    int numero_bits;
 
    int **temp;
 
    temp = (int **)malloc(sizeof(int *)*2);
    temp[0] = (int *)malloc(sizeof(int)*tamanho);
    temp[1] = (int *)malloc(sizeof(int)*tamanho);
 
    numero_bits = (log(maior)/log(2))+1;
 
    for(i = 0; i < numero_bits; i++) {
        k = 0;
        l = 0;
 
        for(j = 0; j < tamanho; j++) {
            if((((vetor[j] >> i) & 0xf) % 2) == 0) {
                temp[0][k++] = vetor[j];
            } else {
                temp[1][l++] = vetor[j];
            }
        }
 
        temp[0][k] = '0';
        temp[1][l] = '0';
 
        j = 0;
        while(temp[0][j] != '0') {
            vetor[j] = temp[0][j];
            j++;
        }
 
        k=0;
        while(temp[1][k] != '0') {
            vetor[j] = temp[1][k];
            j++;
            k++;
        }
 
    }
}
 
void imprime(int *vetor, int tamanho)
{
    int i;
 
    for(i = 0; i < tamanho; i++)
        printf("%d ", vetor[i]);
}

[Radixsort]
[Quicksort]

Dados da máquina onde o teste foi feito:

$ uname -a
Linux myhost 2.6.30-ARCH #1 SMP PREEMPT Mon Aug 17 16:06:45 CEST 2009 x86_64 Intel(R) Core(TM)2 Duo CPU E7400 @ 2.80GHz GenuineIntel GNU/Linux
4GB Ram DDR2 800MHz

Eis então uma tabela comparativa do tempo de execução dos algorítmos para quando o maior valor que pode existir dentro do vetor é o próprio tamanho do vetor.
O tempo é dado em segundos

Numero elementos Radix Sort (s) Quick Sort (s)
10 0.000000000000000 0.000000000000000
1000 0.000000000000000 0.000000000000000
100 mil 0.019999999552965 0.009999999776483
1 milhão 0.310000002384186 0.219999998807907
10 milhões 3.759999990463257 2.549999952316284
100 milhões 42.380001068115234 28.530000686645508

Ao que parece, quicksort vence em todos os testes.
Agora façamos a seguinte alteração nos dois códigos:

No código do quicksort, trocamos

		vetor[i] = rand() % MAX;

por

		vetor[i] = rand() % 1023;

No código do radixsort, trocamos

		vetor[i] = (rand() % tamanho)+1;

por

		vetor[i] = (rand() % 1022)+1;

Eis os novos resultados:

Numero elementos Radix Sort (s) Quick Sort (s)
10 0.000000000000000 0.000000000000000
1000 0.000000000000000 0.000000000000000
100 mil 0.009999999776483 0.009999999776483
1 milhão 0.159999996423721 0.180000007152557
10 milhões 1.559999942779541 2.059999942779541
100 milhões 15.560000419616699 22.809999465942383

O Radixsort vence todas as comparações!
Logo, dependendo do tipo de dados que se deseja ordenar (caracteres por exemplo), o Radixsort pode ter uma performance superior a do Quicksort.

As implementações do Radixsort aqui apresentadas ainda podem ser otimizadas eu acredito, podendo apresentar uma performance ainda maior!

Referências

http://www.itl.nist.gov/div897/sqg/dads/HTML/radixsort.html
http://www.jimloy.com/computer/radix.htm

Tagged as: , , 2 Comments