How to return array from function

We want function which dynamically creates array of integers, so no static declaration is possible because we do not know how big array needs to be.

#include <stdio.h> #include <stdlib.h> int * allocator(int size) { int * array = malloc(size * sizeof(int)); printf ("%p\n", (void*)array); return array; } int main() { int * i = allocator(2); i[0] = 1; i[1] = 2; printf ("%d, %d\n", i[0], i[1]); printf ("%p\n", (void*)i); free(i); return 0; }

We print address of dynamically allocated array in function and in main to be sure it points to the same place and use it to show it is usable. We must use pointer to integer in main because expression

int i[] = allocator(2);

is illegal initializer. It doesn't prevent us from using pointer to integer later as array. Finally we free allocated memory. That what we did is called manual memory allocation and such memory is placed on heap. If we want to use pointer arithmetic we will do something like this:

*i = 1;
*(i + 1) = 2;
printf ("%d, %d\n", *i, *(i+1));

Please note that *(i + 1) is not the same as *i + 1 even if it produces the same output in our example. Dereferencing pointer got higher order of precedence than addition. Because of precedence we are using braces. What if we need two-dimensional array? Two-dimensional array is array of arrays and we will need pointer to pointer.

#include <stdio.h> #include <stdlib.h> int ** allocator(int columns, int rows) { int i; int ** array = malloc(rows * sizeof(int)); printf ("%p\n", (void*)array); for(i = 0; i < columns; i++) array[i] = malloc(columns * sizeof(int)); return array; } int main() { int ** i = allocator(2, 3), j; i[0][0] = 1; i[1][1] = 2; printf ("%d, %d\n", i[0][0], i[1][1]); printf ("%p\n",(void*)i); for(j = 0; j < 3; j++) free(i[j]); free(i); return 0; }

As we see allocating is more complicated and feeing also. It goes in two steps, first we allocate array of rows and after that we allocate arrays of columns. As in previous case we can use array syntax. Now we can go mixing pointer arithmetic and array syntax or go completely pointer arithmetic

*(array+i) = malloc(columns * sizeof(int)); … **i = 1; *(*(i+1)+1) = 2; … *i[0] = 1; *(i[1]+1) = 2; printf ("%d, %d\n", i[0][0], i[1][1]);

Maybe to add that *(*i+1) corresponds to i[0][1]. That was how to do it, now how not to do it. Going heroic and forcing through pointer to array as return type.

typedef int (*array_t)[]; array_t justdont() { static int array[] = { 0, 1}; return &array; } int (*fun())[] { static int array[] = { 0, 1}; return &array; } int main() { int* i = justdont(); printf ("%d, %d\n", i[0], i[1]); printf ("%p\n",(void*)i); i = fun(); printf ("%d, %d\n", i[0], i[1]); printf ("%p\n",(void*)i); return 0; }

This compiles and executes with few warnings, if we change fun return type by adding another pointer int* (*fun())[], that also works. If something looks difficult, like returning address of array, than design needs additional work. For example

int * justdont() { static int array[] = { 0, 1}; return array; }

this is much nicer and easy to understand. Please note that returning static array is perfectly fine and that we do not use free on it. Statically declared variables are allocated automatically and available through execution of program. If we forget to declare array as static than we will return dangling pointer, local variable goes out of scope and returned pointer points to who knows what.

bgfx

Not sure into which category falls this blog entry, maybe review? I never did any 3D programming, so tutorial is out of question. Anyway I was fiddling with bgfx library and here are my experiences.
According to author, Branimir Karadzic, bgfx is cross-platform, graphics API agnostic, "Bring Your Own Engine/Framework" style rendering library. Since it works on Linux I was interested to take a look at it. Best way to get info on project is to visit https://github.com/bkaradzic/bgfx and read README.md. If you are curious to try it you will need also https://github.com/bkaradzic/bx to build it. Library uses custom tool GENie so that is reason why bx is needed, we will need it for projects as well. I have built it and attempted running examples. Some were working, some not. I was attempting to run them from directory where compiled assemblies are but they need runtime. So to run examples go to /examples/runtime and execute ../../.build/linux64_gcc/bin/example-00-helloworldDebug.
How to generate project make files or codeblocks or codelite projects. Since this is cross-platform tool Windows and Mac OS X targets are supported but I will focus on Linux only. Go to bgfx project root and run

../bx/tools/bin/linux/genie --help

To generate make files, do

../bx/tools/bin/linux/genie --gcc=linux-gcc --with-tools gmake

That will produce examples as targets plus tools, like shaderc. Now cd to .build/projects/gmake-linux and execute make help to see options. For 64 bit Linux we may build using this command:

make config=debug64 example-XY-name/nothing for all

If we are using IDE, we go to .build/projects/codeblocks or codelite, open generated project set target and build it. Now we are ready to start fiddling with provided examples.
For somebody who never did any 3D coding, like me, examples provided with bgfx are too complicated so I dumbed them down.

#include "common.h" #include <bgfx.h> int _main_(int /*_argc*/, char** /*_argv*/) { uint32_t width = 1280; uint32_t height = 720; uint32_t debug = BGFX_DEBUG_NONE; uint32_t reset = BGFX_RESET_VSYNC; bgfx::init(); bgfx::setDebug(debug); bgfx::setViewClear(0 , BGFX_CLEAR_COLOR|BGFX_CLEAR_DEPTH , 0x0000ffff , 1.0f , 0 ); while (!entry::processEvents(width, height, debug, reset) ) { bgfx::setViewRect(0, 0, 0, width, height); bgfx::submit(0); bgfx::frame(); } bgfx::shutdown(); return 0; }

This is example-00-helloworld reduced to creating and showing blue window. There is huge amount of plumbing code which makes those examples work in very few lines of code. For example entry::processEvents. In original shape example will print some text and place image on the screen.
Next to be slaughtered is cubes example. It draws now single triangle:

#include "common.h" #include "bgfx_utils.h" #include <stdlib.h> struct PosColorVertex { float m_x; float m_y; float m_z; uint32_t m_abgr; static void init() { ms_decl .begin() .add(bgfx::Attrib::Position, 3, bgfx::AttribType::Float) .add(bgfx::Attrib::Color0, 4, bgfx::AttribType::Uint8, true) .end(); }; static bgfx::VertexDecl ms_decl; }; bgfx::VertexDecl PosColorVertex::ms_decl; static PosColorVertex s_triangleVertices[3] = { { 0.0f, 1.0f, 0.0f, 0xff0000ff }, {-1.0f, -1.0f, 0.0f, 0xffff0000 }, { 1.0f, -1.0f, 0.0f, 0xffff00ff }, }; static const uint16_t s_triangleIndices[3] = { 0, 1, 2 }; float * view; float * proj; int _main_(int /*_argc*/, char** /*_argv*/) { uint32_t width = 1280; uint32_t height = 720; uint32_t debug = BGFX_DEBUG_TEXT; uint32_t reset = BGFX_RESET_VSYNC; bgfx::init(); bgfx::reset(width, height, reset); // Enable debug text. bgfx::setDebug(debug); // Set view 0 clear state. bgfx::setViewClear(0 , BGFX_CLEAR_COLOR|BGFX_CLEAR_DEPTH , 0x303030ff , 1.0f , 0 ); // Create vertex stream declaration. PosColorVertex::init(); // Create static vertex buffer. bgfx::VertexBufferHandle vbh = bgfx::createVertexBuffer( // Static data can be passed with bgfx::makeRef bgfx::makeRef(s_triangleVertices, sizeof(s_triangleVertices) ) , PosColorVertex::ms_decl ); // Create static index buffer. bgfx::IndexBufferHandle ibh = bgfx::createIndexBuffer( // Static data can be passed with bgfx::makeRef bgfx::makeRef(s_triangleIndices, sizeof(s_triangleIndices) ) ); // Create program from shaders. bgfx::ProgramHandle program = loadProgram("vs_cubes", "fs_cubes"); view = (float*)calloc(16,sizeof(float)); view[0] = view[5] = view[10] = view[15] = 1; proj = (float*)calloc(16,sizeof(float)); proj[0] = proj[5] = 0.9; proj[15] = 1; while (!entry::processEvents(width, height, debug, reset) ) { // Use debug font to print information about this example. bgfx::dbgTextClear(); bgfx::dbgTextPrintf(0, 1, 0x4f, "bgfx/examples/01-cube"); bgfx::dbgTextPrintf(0, 2, 0x6f, "Description: Rendering simple static mesh."); bgfx::setViewTransform(0, view, proj); const bgfx::HMD* hmd = bgfx::getHMD(); if (NULL != hmd) { bgfx::setViewRect(0, 0, 0, hmd->width, hmd->height); } else { bgfx::setViewRect(0, 0, 0, width, height); } // Set vertex and fragment shaders. bgfx::setProgram(program); // Set vertex and index buffer. bgfx::setVertexBuffer(vbh); bgfx::setIndexBuffer(ibh); // Set render states. bgfx::setState(BGFX_STATE_DEFAULT); // Submit primitive for rendering to view 0. bgfx::submit(0); // Advance to next frame. Rendering thread will be kicked to // process submitted rendering primitives. bgfx::frame(); } // Cleanup. bgfx::destroyIndexBuffer(ibh); bgfx::destroyVertexBuffer(vbh); bgfx::destroyProgram(program); free(proj); free(view); // Shutdown bgfx. bgfx::shutdown(); return 0; }

Naturally I didn't touch shaders. If we remove completely shaders from code, instances of program, triangle will become black, new achievements of hacking. I also did square, here is the relevant code:

static PosColorVertex s_triangleVertices[4] = { {-1.0f, 1.0f, 0.0f, 0xff0000ff }, {-1.0f, -1.0f, 0.0f, 0xffff0000 }, { 1.0f, -1.0f, 0.0f, 0xffff00ff }, { 1.0f, 1.0f, 0.0f, 0xffff00ff }, }; static const uint16_t s_triangleIndices[6] = { 0, 1, 2, 0, 2, 3 };

Everything else is the same as for triangle. I hope my experience helps somebody to get started. After I spent eternity in programming there are still areas completely unknown to me. If I didn't get friendly suggestion to use identity as view and orthogonal matrix as projection I would stay on Hello World/blue window.

Thread-Specific Data API

Last time it was making function thread safe or reentrant via thread local storage, now more complicated thread-specific data. Everything happens on Linux but should be applicable on any UNIX since we are using POSIX.
This one is also called thread-private data API. It creates per thread and per reentrant function storage. Process starts with key creation when thread for the first time calls function, any thread. In order to have key created once pthread_once is used. Then, once we got key, pthread_getspecific and pthread_setspecific are used to deal with thread-specific data. Here is example:

#include <pthread.h> #include <stdio.h> #include <stdlib.h> #define MAXSSZ 128 static pthread_once_t once = PTHREAD_ONCE_INIT; static pthread_key_t key; static void ctor(void){ int err; err = pthread_key_create(&key, free); if (err != 0){ perror("pthread_key_create failed"); exit(1); } } char * reentrant(char * msg){ char * storage; int err; err = pthread_once(&once, ctor); if (err != 0){ perror("pthread_once failed"); exit(1); } storage = pthread_getspecific(key); if (storage == NULL) { storage = malloc(MAXSSZ); if (storage == NULL){ perror("out of memory?"); exit(1); } err = pthread_setspecific(key, storage); if (err != 0){ perror("pthread_setspecific failed"); exit(1); } } if(msg != NULL) snprintf(storage, MAXSSZ, "Message is '%s'", msg); return storage; } static void * threadFunc(void *arg) { char *str; str = reentrant("Hello I am lumberjack!"); printf("From thread: %s\n", str); sleep(1); str = reentrant(NULL); printf("Again thread: %s\n", str); return NULL; } int main(int argc, char *argv[]) { pthread_t t; int err; char *str; str = reentrant("Hello from main!"); printf("Main message: %s\n", str); err = pthread_create(&t, NULL, threadFunc, NULL); if (err != 0){ perror("pthread_create failed"); exit(1); } str = reentrant("I am not lumberjack!"); err = pthread_join(t, NULL); if (err != 0){ perror("pthread_join failed"); exit(1); } printf("Main message: %s\n", str); return 0; }

To compile it the following is used:

gcc -o tsdexample tsdexample.c -pthread

Result of execution is:

Main message: Message is 'Hello from main!'
From thread: Message is 'Hello I am lumberjack!'
Again thread: Message is 'Hello I am lumberjack!'
Main message: Message is 'I am not lumberjack!'

In order to init static variable once, macro  PTHREAD_ONCE_INIT is used. During key creation destructor, ordinary free, is passed as parameter to pthread_key_create. Main thread and another thread then are calling reentrant function and setting and retrieving messages which are thread specific, what is expected. If we want to delete key it looks like this:

if(&key != NULL){ err = pthread_key_delete(key); if (err != 0){ perror("pthread_key_delete failed"); exit(1); } }

Though in example program it is static variable and it will go out of scope together with program.

Thread safety

Before I start with thread safety, evolving design, if something can be improved don't be lazy.
Now I remembered that GCC supports 128 bit integers and that Fibonacci numbers can be calculated up to forty digits. There is problem with printing those 128 bit integers but somebody already offered utility function to convert them to string, that was Perkins . I only added zero case to his solution.

#include <stdio.h> #include <stdlib.h> void uint128_to_str_iter(__uint128_t n, char *out,int firstiter){ static int offset=0; if (firstiter){ offset=0; } if (n == 0) { return; } uint128_to_str_iter(n/10,out,0); out[offset++]=n%10+0x30; } char* uint128_to_str(__uint128_t n){ if(!n) { char * c = (char*)calloc(sizeof(char),2); c[0] = '0'; return c; } char *out=calloc(sizeof(char),40); uint128_to_str_iter(n, out, 1); return out; } __uint128_t * allocate(uint i){ __uint128_t* res = (__uint128_t*)malloc(2 * sizeof(__uint128_t)); res[0] = 1; res[1] = i; return res; } __uint128_t * recurse(uint i){ if(i<2) return allocate(i); __uint128_t * t = recurse(i/2); __uint128_t t0 = t[0]*t[0]+t[1]*t[1]; t[1] = t[1]*((t[0]<<1)-t[1]); t[0] = t0; if(i%2){ t0 = t[0]+t[1]; t[1] = t[0]; t[0] = t0; } return t; } __uint128_t square(uint n){ __uint128_t * result = recurse(n); __uint128_t r = result[1]; free(result); return r; } uint main(){ uint i; for(i=0; i<187; i++){ char * str = uint128_to_str(square(i)); printf("F(%d) = %s\n", i, str); free(str); } return 0; }

Now we can go up to 186th Fibonacci number without overflow.

Threads

There is few buzz words which people are using in job ads, related to threads, but usually they fail to define them. For example multithreading, concurrency, thread safety, reentrant function and similar.
What is reentrant function or thread safe function? That would be such function that does not use statically allocated or global variables and result depends only on input arguments. For example:

int increase(int i){ static counter = 0; counter += i; return counter; }

is not reentrant or threadsafe. If more than one thread calls that function in the same time, that is concurrency, it will produce unpredictable result. Even if it is called sequentially result will be influenced by previous calls, possibly from different thread. We are on Linux and we are using POSIX threads. So, here is example where two threads are calling nonreentrant function:

#include <stdio.h> #include <stdlib.h> #include <pthread.h> int increase(int i){ static int counter = 0; counter += i; return counter; } void * runnable1(void *arg){ int i; for(i = 1; i < 5; i++){ printf("thread 1 call %d = %d, expecting %d\n", i, increase(2), 2*i); } return((void *)1); } void * runnable2(void *arg) { int i; for(i = 1; i < 5; i++){ printf("thread 2 call %d = %d, expecting %d\n", i, increase(3), 3*i); } pthread_exit((void *)2); } int main(void) { int error; pthread_t tid1, tid2; void * exitCode; error = pthread_create(&tid1, NULL, runnable1, NULL); if (error != 0){ perror("create thread 1 failed"); exit(1); } error = pthread_create(&tid2, NULL, runnable2, NULL); if (error != 0){ perror("create thread 2 failed"); exit(1); } error = pthread_join(tid1, &exitCode); if (error != 0){ perror("join thread 1 failed"); exit(1); } printf("thread 1 exit code %ld\n", (long)exitCode); error = pthread_join(tid2, &exitCode); if (error != 0){ perror("join thread 2 failed"); exit(1); } printf("thread 2 exit code %ld\n", (long)exitCode); return 0; }

To compile it I used:

gcc -o notreentrant -g -gdwarf-2 notreentrant.c -pthread

That -gdwarf-2 story is distraction, I am using gcc version 4.9.2 built from source on old Mint based on Ubuntu 12.04 and new gcc outputs by default debug info which gdb doesn’t understand. If you just built new gcc and gdb doesn't work as expected, that is cure. Linker switch -pthread is required.
Example output from execution looks like this:

thread 1 call 1 = 2, expecting 2
thread 1 call 2 = 7, expecting 4
thread 1 call 3 = 9, expecting 6
thread 1 call 4 = 11, expecting 8
thread 2 call 1 = 5, expecting 3
thread 2 call 2 = 14, expecting 6
thread 2 call 3 = 17, expecting 9
thread 2 call 4 = 20, expecting 12
thread 1 exit code 1
thread 2 exit code 2

How do we remedy that? We make increase reentrant or thread safe:

int increase(int i){ static __thread int counter = 0; counter += i; return counter; }

We requested thread local storage, available since kernel 2.6 and supported by gcc 3.3 and higher. When thread terminates and goes out of scope, local thread storage does to. We could make function reentrant via using thread specific data, but that is much more code, maybe next time.

Similarities and differences between Java and C

As almost everybody knows Java was intended to be sort of C without pointers and to be used in smart home appliances. And then Gosling and his co-workers were drinking too much coffee and it come out more like C++ without pointers, resource management and multiple inheritance. So, how really different can be coding the same task in C and Java. Practically we are addressing just resource management part, we will not try to go OO in C. We implement fast Fibonacci algorithm in Java.

import java.math.BigInteger; public class Example { /* |1 1| <- we are using this row |1 0| */ static BigInteger[] Q = {BigInteger.ONE, BigInteger.ONE}; static BigInteger[] I = {BigInteger.ONE, BigInteger.ZERO}; /* F(2n+1) = F(n+1)^2+F(n)^2 F(2n) = F(n)*(2*F(n+1)−F(n)) */ static BigInteger[] squareIt(int n) { if (n == 0) return I; if (n == 1) return Q; BigInteger[] t = squareIt(n / 2); BigInteger[] res = {t[0].multiply(t[0]).add(t[1].multiply(t[1])), t[1].multiply(t[0].shiftLeft(1).subtract(t[1]))}; if(n%2==1) return new BigInteger[]{res[0].add(res[1]), res[0]}; return res; } public static void main(String[] env) { for(int i = 0; i<1000; i++) System.out.println("F("+i+") = "+squareIt(i)[1]); } }

It is really very fast, logarithmic fast. Now we try to do the same in C, we will skip BigInteger and use long. The first thing which wont work in C is array allocation, declaring them like that will create local variable. So we create function to allocate our arrays in order to return them to caller:

long * allocate(){ return (long*)malloc(2 * sizeof(long)); }

Those Q and I arrays can be global variables and we are ready to go, almost. We can not declare t as array since initializer can't be function call. Not big deal, we will use pointer to long but that introduces another problem. That pointer may point to global variables Q and I or allocated variables. Again we can check to what it points and if it doesn't point to Q or I we can call free. Already we have feeling that those arrays returned by squaring function are somehow polymorphic and in Java they were not like that. Allocating and freeing arrays on every step of recursion and checking should we free or not doesn't look right. Attempt to literally translate Java code may look something like this:

long * recurse(int i){ if(i==0) return I; if(i==1) return Q; long * t = recurse(i/2); long * r = allocate(); r[0] = t[0]*t[0]+t[1]*t[1]; r[1] = t[1]*((t[0]<<1)-t[1]); checknfree(t); if(i%2){ long * odd = allocate(); odd[0] = r[0]+r[1]; odd[1] = r[0]; free(r); return odd; } return r; }

If Q and I arrays are created using malloc we could do reassignments of initial array and free it at the end of calculation. So no global Q and I but they will be allocated on the heap and serve all recursive calls.

#include <stdio.h> #include <stdlib.h> long * allocate(int i){ long* res = (long*)malloc(2 * sizeof(long)); res[0] = 1; res[1] = i; return res; } long * recurse(int i){ if(i<2) return allocate(i); long * t = recurse(i/2); long t0 = t[0]*t[0]+t[1]*t[1]; t[1] = t[1]*((t[0]<<1)-t[1]); t[0] = t0; if(i%2){ t0 = t[0]+t[1]; t[1] = t[0]; t[0] = t0; } return t; } long square(int n){ long * result = recurse(n); long r = result[1]; free(result); return r; } int main(){ int i; for(i=0; i<93; i++) printf("F(%d) = %ld\n", i, square(i)); return 0; }

Now if we run valgrind on our program it will report:

==18896== HEAP SUMMARY:
==18896==     in use at exit: 0 bytes in 0 blocks
==18896==   total heap usage: 93 allocs, 93 frees, 1,488 bytes allocated
==18896== 
==18896== All heap blocks were freed -- no leaks are possible

Looking now at initial Java solution, it also could benefit from array reuse instead of allocating new ones on every step. Garbage collector is nice thing but not thinking at all about resource management may result in lower quality of code.

Functions, pointers, nested functions, clojures

While ago GCC introduced nested functions which are allowing easy solution of funargs problem. With funargs problem solved C can accommodate closures and provide everything what different modern languages can, except funny syntax. In order to demonstrate that we need function pointers and they are not part of modern languages derived from C. So, we start with function pointers refresher.

Function pointers

They are more cryptic than other pointers and they scare away younger generations of programmers. We will declare function and pointer to it:

int add( int i, int j) { return i + j; } int main(){ int (*fptr)(int,int) = add; printf("%d\n", fptr(2, 3)); return 0; }

We have return type, variable name in brackets and type list for parameters. We could also use address of operator during assignment. How we pass it as parameter?

int execute(int (*fptr)(int,int), int a, int b) { return fptr(a, b); } printf("%d\n", execute(fptr, 2, 3));

How we return it, write function which is returning function pointer? There is user friendly way, we declare typedef for function pointer:

typedef int (*fptra)(int,int); fptra get() { return add; } printf("%d\n", get()(2, 3));

But there is also user unfriendly way:

int (*fptr2fptr())(int, int) { return add; } int (*fptr)(int,int) = fptr2fptr(); printf("%d\n", fptr(2, 3)); printf("%d\n", (*fptr2fptr())(2, 2));

Return type of return type goes in front, inside first pair of brackets name and own parameters, inside the second pair of brackets list of parameter types of return type. And just for those who value extreme programming above all, pointer to function which returns function which returns function:

int (*(*f2f2f())(void))(int,int){ return fptr2fptr; } int (*(*f2f)())(int,int) = f2f2f(); int (*fadd)(int,int) = f2f(); printf("%d\n",fadd(4, 3));

Nested Functions

Those are GCC specifics and they are just ordinary functions but inside another function. Standard C will not let you define function within function. To demonstrate we can do trampolining of one notoriously recursive function.

int factorial(int i){ int acc = 1; void * f() { if (i < 2) return 0; acc *= i; --i; return f; } while(f()) continue; return acc; } int main() { int i; for(i = 2; i < 14; i++) printf("%d! is %d\n", i, factorial(i)); return 0; }

Here nested function is accessing variables which are host scoped and returns self or NULL. We got whole function back but we are not really using it, except to exit while loop if it is NULL. To remedy that we will dispatch those nested functions out of host function.

struct closure{ int counter; int visited; } data; int (*pgen(int i))(int){ int increase(int k){ data.counter += k; data.visited++; return data.counter; } int decrease(int k){ data.counter -= k; data.visited++; return data.counter; } if(i) return increase; return decrease; } int main(){ data.counter = 4; data.visited = 4; int (*incr)(int) = (*pgen)(1); int (*decr)(int) = (*pgen)(0); printf("counter = %d visited = %d\n", data.counter, data.visited); printf("counter = %d\n", incr(2)); printf("counter = %d\n", decr(1)); printf("counter = %d visited = %d\n", data.counter, data.visited); printf("counter = %d\n", incr(4)); printf("counter = %d\n", decr(1)); printf("counter = %d visited = %d\n", data.counter, data.visited); return 0; }

Here we have struct which is global variable, it is not really closure but just named like that, and one function generator. If we pass zero to generator we get one function back and for non zero the other one. Both nested functions now operating on global variable. In previous example we were working on variables belonging to host function, in this case if we try the same we will have problems since host's variables went out of scope.

Funarg problem

This one is related to functional languages and nested functions. If we pass nested function as argument to another function that is called downward funarg problem. Returning nested function to caller is upwards funarg problem. If any of passed or returned functions references variable from scope of declarer, those must be available to nested function so it can be executed. Passing nested function as argument to another function is not problem in C since all local variables of declarer are alive.

void execute(void(*fptr)(int), int i){ fptr(i); } void print(int(*fptr)(void)){ printf("counter = %d\n", fptr()); } int main(){ int counter = 0; void increase(int k){ counter += k; } void decrease(int k){ counter -= k; } int getter(){ return counter; } execute(increase, 3); print(getter); execute(decrease, 2); print(getter); return 0; }

Returning nested function is more problematic because declarer went out of scope. That could be addressed like in Java, if anonymous function uses some declarer's variable, that variable must be final. In C there is even more options, we can use static variables or allocate them on the heap.

int (*pgen(int i))(int){ static int counter = 0; int increase(int k){ counter += k; return counter; } int decrease(int k){ counter -= k; return counter; } if(i) return increase; return decrease; } int main(){ int (*incr)(int) = (*pgen)(1); int (*decr)(int) = (*pgen)(0); printf("counter = %d\n", incr(2)); printf("counter = %d\n", decr(1)); printf("counter = %d\n", incr(4)); printf("counter = %d\n", decr(1)); return 0; }

It is easier to use static than malloc and free. Removing static from counter declaration will result in exception during execution of program. So, funargs problem is not problem for C in GCC version. This makes C almost functional language, functional with few extras.

Recursive Squaring

This one neatly solves problem of log n multiplications to calculate nth power of of some number. It goes in two stages, calls self recursively with half power until it reaches power 0 or 1 and then squares results. If power was not even correction is required, additional multiplication with base.

static long pow(long base, int n){ if(n==0) return 1; if(n==1) return base; long acc = pow(base, n/2); acc *= acc; if(n%2 == 0) return acc; return base * acc; }

If we want to see how it works and how many times it calls self we can insert debug code:

static long powDbg(long base, int n, int depth){ for(int i = 0; i < depth; i++) System.out.print("-"); System.out.println(n); if(n==0) return 1; if(n==1) return base; long acc = powDbg(base, n/2,depth+1); acc *= acc; for(int i = 0; i < depth; i++) System.out.print("-"); System.out.println(n); if(n%2==0) return acc; return base * acc; }

Should have been included it in previous blog entries about recursion but better late than never.

Recursion and optimizations

As seen in last instalment calculating Fibonacci numbers using naïve implementation creates binary tree of calls and the same subtrees are calculated over and over again what has huge impact on performance. One of ways to improve performance is caching previous results and that is called memoization. In old Java we will do something like this:

private static final Map<Long, Long> memo = new HashMap<>(); static { memo.put(0L, 0L); memo.put(1L, 1L); } public static void main(String[] args) { System.out.println(fibonacci(0)); System.out.println(fibonacci(43)); System.out.println(fibonacci(92)); } static long fibonacci(long n) { if (memo.get(n) == null) { memo.put(n, fibonacci(n - 1) + fibonacci(n - 2)); } return memo.get(n); }

We initialize map container with result for first two cases, to avoid if statements, and later if there is mapping for input, we return it and if there is no mapping for input we calculate it. Without help from memoization naïve implementation will choke and die on input 92. Using Java 8 it becomes even shorter, here is just fibonacci method, everything else is the same:

public static long fibonacci(long n) { return memo.computeIfAbsent(n, m -> fibonacci(m - 1) + fibonacci(m - 2)); }

Signature of computeIfAbsent is the following:

default V computeIfAbsent(K key, Function<? super K, ? extends V> mappingFunction)

where n is key and function is interface Function, we defined it recursively and compiler was not upset. Memoization brings speed but takes space.
How do we turn tail call version into lambda. We can do this for example:

static Function<Long, Long> y = new Function<Long, Long>() { private long tailRecursive(long a, long b, long n) { if (n > 0) { return tailRecursive(b, a + b, n - 1); } return a; } public Long apply(Long x) { return tailRecursive(0, 1, x); } }; public static long fibonacciFunctional(long n) { return y.apply(n); }

what does not look very functional and very Java 8 because it is anonymous class. Another more Java 8 option is declaring functional interface with sufficient number of arguments:

interface FewArgs<T, U, V, R> { public R apply(T t, U u, V v); } static FewArgs<Long, Long, Long, Long> tailRecursive; static { tailRecursive = (a, b, n) -> { if (n > 0) return tailRecursive.apply(b, a + b, n - 1); return a; }; }

We also can not initialize it in one go, declare class field and initialize it, it must be done from some of host class methods.
Remaining optimization ideas are strictly related to Fibonacci numbers and not applicable to most other recursions. There is not so obvious way to generate Fibonacci numbers:

| 1 1 | n      | F(n+1) F(n)   |
| 1 0 |    =   | F(n)   F(n-1) |

If we rise matrix on the left, lets call it A, to power of n-1 its A[0][0] is Fn. We can prove it using induction. For n=1 our claim obviously holds, now

| F(n+1) F(n)   |   | 1 1 |
| F(n)   F(n-1) | x | 1 0 | =

|F(n+1)+F(n) F(n+1)+0|   |F(n+2) F(n+1)|
|F(n)+F(n−1) F(n)+0  | = |F(n+1) F(n)  |

Instead of doing chain multiplication we can do squaring and reduce number of multiplications to log(n) multiplications.

| F(n+1) F(n)   |2
| F(n)   F(n-1) | =

| F(n+1)^2+F(n)^2         F(n+1)*F(n)+F(n)*F(n−1)|
| F(n)*F(n+1)+F(n−1)*F(n) F(n)^2+F(n−1)^2        |

From where one can pull more interesting relations.
Here is illustrative implementation of this idea:

private static long square(long n) { if(n==0) return n; long test = n-1; ArrayList<long> arr = new ArrayList<>(); long[] A = new long[]{ 1, 1, 1, 0 }; long[] result = new long[4]; System.arraycopy(A, 0, result, 0, 4); while(test>1){ if(test%2!=0) arr.add(1L); else arr.add(0L); test/=2; } for(int i=arr.size()-1;i>-1;i--){ result=multiply(result,result); if(arr.get(i)==1) result=multiply(result,A); } return result[0]; } private static long[] multiply(long[] a, long[] b) { long[] result = new long[]{ a[0] * b[0] + a[1] * b[2], a[0] * b[1] + a[1] * b[3], a[2] * b[0] + a[3] * b[2], a[2] * b[1] + a[3] * b[3] }; return result; } public static void main(String[] arg) { for (int i = 0; i < 93; i++) { System.out.println("F(" + i + ") = " + square(i)); } }

That ArrayList is mapper, consumes some additional space but makes things easier to debug.

Recursion

Got recently homework to do as part of interview process, already described it here. After providing them with working solution, thoroughly tested, they decided not to speak with me. My conclusion is that was not so bad outcome.
While searching web for clues I discovered this Find all paths from source to destination.
Guy submits his code, peers do review, everything sounds nice.
I will refrain from reviewing particular solution but after looking at it I decided to write about recursion. My impression is that young generations of programmers are putting great effort into mastering frameworks and new features of programming languages but in the same time somehow missing basics. It also coincides with my exploration of new functional features of Java 8.
Recursion is when method, in Java, or function, in C, calls itself during execution. Execution of caller is suspended until callee finishes, then it proceeds. So, there must be some natural way for recursive calls to stop propagating themselves into infinity. If we are not traversing some kind of tree structure, or marking vertices as visited when traversing graph, we typically passing variable to control depth of recursion. Trivial example of recursion is calculation of factorial:

n! = n*(n-1)*(n-2)*...*2*1

It is defined for positive integers. Here is natural and naive implementation together with test:

public static void main(String[] args) { System.out.println("9! = "+factorial(9)); } static long factorial(long n){ if(n<0) return 0; if(n<2) return 1; return n*factorial(n-1); }

I could throw exception on negative n and assert is 9! equal to 362880. What stops recursion here from proceeding forever is that if statement in combination with decreasing n in recursive calls. Now in order to visualize how execution looks like we will add some debugging code.

public static void main(String[] args) { System.out.println("9! = " + factorial(9, 9)); } static long factorial(long n, long m) { for (int i = 0; i < m - n; i++) { System.out.print("-"); } System.out.println(n); if (n < 0) return 0; if (n < 2) return 1; long result = n * factorial(n - 1, m); for (int i = 0; i < m - n; i++) { System.out.print("-"); } System.out.println(n); return result; }

Code now prints current value of n and number of dashes is how deep we are into recursion. We can see how stack of frames is growing in debugger as well but this is nicer. Output of execution is:

9
-8
--7
---6
----5
-----4
------3
-------2
--------1
-------2
------3
-----4
----5
---6
--7
-8
9
9! = 362880

As expected it behaves in linear fashion. It can be much more interesting than linear, for that we will calculate Fibonacci numbers. In the case of Fibonacci numbers we have the following recursive definition:

Fn = Fn-1 + Fn-2

Trivial and naive implementation looks like this:

public static void main(String[] args) { System.out.println("F42 = " + fibonacci(42)); } static long fibonacci(long n) { if (n < 1) return 0; if (n < 2) return 1; return fibonacci(n - 1) + fibonacci(n - 2); }

Execution of this code will take  about second or two. Using 42 as function argument illustrates point that naïve implementation is not really most efficient. Mechanics of stopping recursive calls is identical to one in the case of factorial. Let us insert debugging code and see what is going on.

static int indent = 0; public static void main(String[] args) { System.out.println("F5 = " + fibonacci(5, 5)); } static long fibonacci(long n, long m) { for (int i = 0; i < m - n + indent; i++) { System.out.print("-"); } System.out.println(n); if (n < 1) return 0; if (n < 2) return 1; long result = fibonacci(n - 1, m); for (int i = 0; i < m - n + indent; i++) { System.out.print("-"); } System.out.println(n); indent--; result = result + fibonacci(n - 2, m); indent++; for (int i = 0; i < m - n + indent; i++) { System.out.print("-"); } System.out.println(n); return result; }

For reasonably sized argument of 2 we have this output:

2
-1
2
-0
2
F2 = 1

Argument 2 is decreased and recursive call is made with argument 1. When it returns to level 0 with result 1 next recursive call is made with argument 0. The second call also returns with result 0 to level 0 and we have 1+0=1 as result.
We try now 3 as argument:

3
-2
--1
-2
--0
-2
3
-1
3
F3 = 2

We start on level 0 with 3-1 call, then we have pattern from last run elevated for one level up, return to level 0 with result 1 and finally 3-2 recursive call and its return with result 1 to level 0. We can make conclusion that recursive calls are building binary tree. We can try argument 4 to recognize patterns for 3 and 2 and so on.
We could achieve the same on the factorial if we were splitting chain multiplication in half, I was writing about that here.
Recursive call does not have to return anything, work can be performed on one of function arguments. For example we can reverse array by swapping elements.

public static void main(String[] args) { String r = "Hello World!"; char[] arr = r.toCharArray(); reverse(arr, r.length()-1, 0); System.out.println(arr); } public static void reverse(char[] array, int hi, int lo) { if(lo>hi) return; char tmp = array[hi]; array[hi] = array[lo]; array[lo] = tmp; reverse(array, hi - 1, lo + 1); }

Here recursion stops when lo is bigger than hi.

Optimization

We can improve on naïve implementation using tail call elimination. We write method in such way that return statement is pure function call and that allows compiler to perform optimization and practically discard frames without placing them onto the stack. For example factorial can be rewritten like this:

static long tailRecursive(long acc, long n) { if (n == 0) return acc; return tailRecursive(n * acc, n - 1); }

We provide 1 as initial value for accumulator, n is the same as before. For Fibonacci we will need two accumulators:

static long tailRecursive(long a, long b, long n) { if (n > 0) return tailRecursive(b, a + b, n - 1); return a; }

Initial values are a = 0, b = 1 and n as before. From this form is quite easy to switch to iterative form.

static long notRecursive(long n) { if(n == 0) return 0; long a = 0; long b = 1; long c; do { c = b; b = a + b; a = c; --n; } while (n > 0); return a; }

About other forms of optimization like memoization, repeated squaring and lambdas, next time.

Scheduling and Johnson's Algorithm

There is actually paper available on the web where Johnson's scheduling algorithm is nicely described. You can find it here.
So that is not all-pairs shortest paths but scheduling one. For those who are lazy to read it, here is the story. We have programming task paired with testing task. We have one programmer and one tester, tester can test only what is coded, must wait for programmer to finish. In which order jobs should be executed so that whole bunch is finished in shortest time?
Now Johnson's algorithm in plain English:

1. For each task T1, T2, ..., Tn, determine the P1 and P2 times.
2. Establish two queues, Q1 at the beginning of the schedule and Q2 at the end of the schedule.
3. Examine all the P1 and P2 times, and determine the smallest.
4. If the smallest time is P1, then insert the corresponding task at the end of Q1. Otherwise, the smallest time is P2, and the corresponding task is inserted at the beginning of Q2. In case of a tie between a P1 and P2 time, use the P1 time.
5. Delete that task from further consideration.
6. Repeat steps 3 – 5 until all tasks have been assigned.

Very simple. Here is the code:

import java.util.ArrayDeque; import java.util.Deque; public class JavaApplication18 { public static void main(String[] args) { Tasks[] t = new Tasks[6]; for (int i = 0; i < 6; i++) { t[i] = new Tasks(); } t[0].programming = 15; t[0].testing = 4; t[1].programming = 22; t[1].testing = 10; t[2].programming = 3; t[2].testing = 6; t[3].programming = 12; t[3].testing = 8; t[4].programming = 9; t[4].testing = 6; t[5].programming = 12; t[5].testing = 9; executeAll(6, t); } static void executeAll(int n, Tasks[] t) { int[] order = getSchedule(n, t); System.out.println("order of execution is:"); for (int i = 0; i < n; i++) { System.out.println((char)(65+order[i]) + " " + t[order[i]].programming + " -> " + t[order[i]].testing); } System.out.println("Complete duration is " + calculateDuration(n, order, t)); } static int calculateDuration(int n, int[] order, Tasks[] t) { int max = 0, programmingTime = 0; int testingTime = 0, carry; for (int i = 0; i < n; i++) { programmingTime += t[order[i]].programming; carry = 0; if (testingTime > programmingTime) carry = testingTime - programmingTime; testingTime = programmingTime + t[order[i]].testing + carry; if (testingTime > max) max = testingTime; } return testingTime; } static int[] getSchedule(int n, Tasks[] t) { int[] order = new int[n]; int counter = 0, taskIndex = -1, min; Deque<integer> programmingDeque = new ArrayDeque<>(); Deque<integer> testingDeque = new ArrayDeque<>(); for (int j = 0; j < n; j++) { min = 2 ^ 16; kind task = null; for (int i = 0; i < n; i++) { if (!t[i].removed) { if (t[i].programming < min) { min = t[i].programming; task = kind.programming; taskIndex = i; } if (t[i].testing < min) { min = t[i].testing; task = kind.testing; taskIndex = i; } } } if (task == kind.programming) { programmingDeque.addLast(taskIndex); t[taskIndex].removed = true; } if (task == kind.testing) { testingDeque.addFirst(taskIndex); t[taskIndex].removed = true; } } while (programmingDeque.size() > 0) { order[counter++] = programmingDeque.removeFirst(); } while (testingDeque.size() > 0) { order[counter++] = testingDeque.removeFirst(); } return order; } enum kind { programming, testing }; static class Tasks { public int programming; public int testing; public boolean removed; } }

I think that code is self-explanatory and I will not bother you with explanation.
I also had a look at C++ code written by riteshkumargupta here it is slight modification of problem mentioned in paper. If you prefer C++ to Java take a look.

Find all paths from vertex to vertex in directed, weighted graph

Why I am doing homework? I have applied for a job at some company. They did some kind of phone interview and decided to give me “task” before they make up their mind will they ever talk to me about job again. OK, that evening I open Firefox and find in mailbox my tasks. That was few NP hard problems and I can pick and choose. Time for execution is limited, to make everything more interesting. Like reality show but you do not eat snails. Going quickly through them I realized that I am dealing with person who thinks that he or she is lecturer at university and that I do not have clue how to solve any of tasks.
BTW my educational background is physicist. Studied technical physics at school of electrical engineering at University of Belgrade. Never finished. As programmer I am self-educated. While I didn’t get university diploma I have learned to use literature.
Natural choice was the first task, it must be more important than others, that is the reason why it is the first. Task was something about finding all paths from vertex A to vertex B on weighted directed multigraph with cycles. Never encountered graphs before, only trees and they are part of the family. OK, time to google hard and yandex hard. After whole day of searching I got two clues. The first one is complete solution to generate all possible paths up to some length on Perl Monks and the other one is notes from some lecturing by Gunnar. Both clues claiming that adjacency matrix is a way to go and tool of choice is modified Floyd-Warshall. Learned quickly what I need and assembled working program based on degenerated matrix multiplication. After three days almost without sleep, submitted solution and went to bed. Next morning I figured out much simpler solution.
Now friendly warning, what follows is theory by somebody who studied graph theory for one whole day.

Adjacency matrix

Graph is one with directed edges, one way edges, multiple edges between pair of vertices may exist and there may be edge starting and finishing at the same vertex. Aadjacency matrix gives metric on graph space, who is neighbour to whom. Where you can go from some place.

     A    B    C    …    Z
A    1    0    0        1

From vertex A we can travel to A or Z in one move, to B and C we can't go in one move. Now comes curious tourist with question how far is B from A. If we start multiplying adjacency matrix with itself we will find out where we can go in two or more moves. Number of moves corresponds to power of adjacency matrix and if after rising it to power of 4 we see something like this:

     A    B    C    …    Z
A    0    2    1        4

meaning of it is we can't reach A in 4 moves, there is only one path to C but 2 paths to B and 4 paths to Z. Now to find out how many different paths exist from A to B where length of path is up to 4 moves one sums adjacencyMatrix[0][1] for powers from 1 to 4. Those are not shortest pats those are all possible paths including loops and cycles. Good reason to keep binary adjacency matrix around even if we got weights for edges.
So how do we go about finding those two or more moves paths, translating them in chain of edges? Since adjacency matrix is metric on that space we can use it to discover paths. For example where we can go from A in two moves.

     A    B    C    D
A    0    0    1    1
B    1    0    1    0
C    0    1    0    1
D    1    1    0    0

We can reach C or D in single move, from C we can get to B and D in single move and from D we can get to A and B in one move. So paths are ACB, ACD, ADA and ADB. That should be enough to write code for pathfinder.

Java implementation

From some vertex we can go where adjacency matrix says we can go and we move to next vertex. For that next vertex, we are in different row but we do the same thing. Looks like case for recursion. When frames start offloading from the stack they will bring us story where repeated recursive calls were taking them. We just need to send story back. If we do not limit depth of recursion it will never stop. BTW replace depth-1 with depth-- and --depth if you are not sure how suffix and prefix are different.

public ArrayList<string> findPath(int start, int end, int depth) { ArrayList<string> result = new ArrayList<>(); ArrayList<string> tmp = new ArrayList<>(); for (int i = 0; i < adjacencyMatrix.length; i++) { if (adjacencyMatrix[start][i] != 0) { if (i == end) { result.add(start + "," + end); } if(depth>1) tmp.addAll(findPath(i, end, depth-1)); } } if(tmp.size()>0) for (String s : tmp) { result.add(start + "," + s); } return result; }

We have built binary adjacencyMatrix and we are traversing graph. If we find target we store it but we proceed looking for target via calling all vertices which are reachable in single move. With limitation that depth of traversal is bigger than 1 and on each new call we are using reduced incoming depth. When function call comes back we append prefix to all paths. If you print content of that string array you will get something like:

2,3,2,3,2,3
2,3,2,4,1,2,3
2,3,4,1,2,3

It is the same format as used by blokhead from PerlMonks in his solution and you can nicely compare output of your code. Do not forget to make  adjacency matrix match :P
What about weights? Instead of having binary adjacency matrix we replace 1 with corresponding edge, actually weight of that edge. I will use integers for weights.

     A    B    C    D
A    0    0    3    5

Slightly modified traversal code:

public HashMap<String, Integer> findWeightedPath(int start, int end, int depth) { HashMap<String, Integer> result = new HashMap<>(); for (int i = 0; i < adjacencyMatrix.length; i++) { if (adjacencyMatrix[start][i] != 0) { if (i == end) { result.put(start + "," + end, adjacencyMatrix[start][i]); } if (depth > 1) { HashMap<String, Integer> tmp = findWeightedPath(i, end, depth - 1); if (tmp.size() > 0) { for (String s : tmp.keySet()) { result.put(start + "," + s, tmp.get(s) + adjacencyMatrix[start][i]); } } } } } return result; }

When you print content of that HashMap it looks like this:

2,3,2,3 = 23
2,3,4,1,2,3 = 33
2,3,2,4,1,2,3 = 35

It still stops on depth. If you need it to stop on weights sum then send data upstream. Don’t go for matrix multiplication if you do not have to, check other options. Nobody is going to ask you to find all paths between all vertices, to find all paths between two nodes recursive traversal looks much better than matrix multiplication.

PostgreSQL and PL/Pythonu

Beside SQL PostgreSQL has support for creating functions in many other so called procedural languages. We just need to install them and they are available within database. Standard distribution includes those four languages PL/pgSQL, PL/Tcl, PL/Perl, and PL/Python. There are other externally maintained languages like Java, Scheme, PHP, Unix shell and so on. If that is not enough it is possible to write procedural language handler, it is open source and consequently fully customizable.
Hosted Python is very interesting option. It is quite powerful, versatile language and code is very concise. If it is executed within PostgreSQL it enjoys very fast communication avoiding overhead of DB calls which external programs suffer. Both Python 2 and Python 3 are supported. If you are using, for example, Ubuntu then postgresql-plpython-9.x is Python 2 and postgresql-plpython3-9.x is Python 3. It needs to be installed, server restarted and language created.

createlang plpythonu [dbname]

If we do not specify db, then it will be created on db to which are we currently connected. That u suffix have meaning that we are dealing with unrestricted/untrusted language. Power is given to us and it is our responsibility to use it to our benefit without causing damage. Once everything is in place we can write Hello World function in Python.

CREATE OR REPLACE FUNCTION lumberjack() RETURNS text AS $BODY$ a = "I'm a lumberjack, and I'm okay. " return a + "I sleep all night and I work all day." $BODY$ LANGUAGE plpythonu

From psql we try it

testpgdb=# select lumberjack(); lumberjack ----------------------------------------------------------------------- I'm a lumberjack, and I'm okay. I sleep all night and I work all day. (1 row)

Or if you prefer GUI, then use pgAdmin III start Query tool and execute it there.
To enjoy full comfort of Python one need to import modules. How do we use modules in PL/Pythonu? It is no different to ordinary Python code:

create or replace function get_response() returns text as $BODY$ import httplib2 h = httplib2.Http(".cache") resp, content = h.request("http://example.org/", "GET") return resp $BODY$ language 'plpythonu'; select get_response();

we should see response. If we see:

ERROR: ImportError: No module named httplib2 CONTEXT: Traceback (most recent call last): PL/Python function "get_response", line 2, in <module> import httplib2 PL/Python function "get_response"

That is sign that httplib2 is not in path, we can place it in path or if we do not want to reload server we can append to function:

from sys import path path.append('/path to httplib2')

There is significant quantity of examples on Internet to help you further. With some knowledge of Python one can do really huge amount of work directly from PostgreSQL without overhead of DB connection. In role of mailer of monthly report PL/Pythonu function beats Java Application Server and even C/C++ application connecting from outside to DB.

GCC Quad-Precision Math Library example

There is small PDF on GCC website with libquadmath documentation. It contains explanation and example how to convert string to __float128 and how to convert __float128 back to string. It took me some googling to find out how to compile example and to figure out what is width parameter. To compile it we need -lquadmath linker switch and width, fourth argument to quadmath_snprintf, is overall width of output, right anchored. Here is my version of mentioned example:

#include <quadmath.h> #include <stdlib.h> #include <stdio.h> int main(){ __float128 pi = strtoflt128 ("3.1415926535897932384626433832795029", NULL); __float128 dpi = 3.1415926535897932384626433832795; __float128 qpi = 3.1415926535897932384626433832795029Q; char buf[128]; int n = quadmath_snprintf (buf, sizeof buf, "%*.34Qf",10, pi); if ((size_t) n < sizeof buf) printf ("%s\n", buf); n = quadmath_snprintf (buf, sizeof buf, "%*.34Qf", 40, dpi); if ((size_t) n < sizeof buf) printf ("%s\n", buf); n = quadmath_snprintf (buf, sizeof buf, "%*.34Qf", 45, qpi); if ((size_t) n < sizeof buf) printf ("%s\n", buf); return 0; }

If we try to assign number without Q suffix then we will have ordinary double inside quad storage. The fourth argument to quadmath_snprintf will be ignored if it is smaller than length of string representation, no spaces will be added in front of number.
If we need to examine layout in memory, we can use modified example from last article. We do not need here -lquadmath linker switch to compile.

#include <stdio.h> #include <quadmath.h> #include <stdlib.h> typedef struct quad_storage { unsigned long mant_low:64; unsigned long mant_high:48; unsigned int exponent:15; unsigned int sign:1; } Quad; typedef union quad_helper { __float128 number; Quad storage; } Helper; int main(){ Helper quad = {3.1415926535897932384626433832795029Q}; printf("sign = %u, exponent = %u, mantissa_h = %lu, mantissa_l = %lu\n", quad.storage.sign, quad.storage.exponent, quad.storage.mant_high, quad.storage.mant_low); printf("sign = %X, exponent = %X, mantissa_h = %lX, mantissa_l = %lX\n", quad.storage.sign, quad.storage.exponent, quad.storage.mant_high, quad.storage.mant_low); return 0; }

Interesting part from GDB session:

(gdb) print &pi $1 = (__float128 *) 0x7fffffffe650 (gdb) x/4 0x7fffffffe650 0x7fffffffe650: 0xc51701b8 0x8469898c 0xb54442d1 0x4000921f

We compare it with output of program to see what we are looking at.
Examining tests from GCC source, we find out that all operators which are supported for other floats are supported on __float128, so no funny function calls to add or compare two numbers. Comfort of normal syntax and 34 significant digits. For example we do square root round trip:

dpi = 2.0; qpi = sqrtq(dpi); n = quadmath_snprintf (buf, sizeof buf, "%*.34Qe", 45, dpi-qpi*qpi); if ((size_t) n < sizeof buf) printf ("%s\n", buf);

That code produced -3.8518598887744717061119558851698546e-34 error on my box.

Double precision floats and GDB

They are not human readable in memory like integers, so if you are examining memory in GDB session it is not very likely that you will recognize what number is it. Their layout in memory is specified by IEEE 754 standard. They occupy 64 bits and there is one bit for sign, 11 bits for exponent and 52 bits significand. There is really plenty of read about floats on the web and I will not bother you with theory, significance of their discrete spectrum, rounding and similar. Just want to show you how to handle them during debug session.
I have seen in some lecture, can't remember exactly where, solution to examine memory layout via struct and union. Here is my interpretation of it:

#include <stdio.h> typedef struct double_storage { unsigned long mantissa:52; unsigned int exponent:11; unsigned int sign:1; } Double; typedef union double_helper { double number; Double storage; } Helper; int main(){ Helper dbl = {3.141592653589793}; printf("number = %16.15lf sign = %u, exponent = %u, mantissa = %lu\n", dbl.number, dbl.storage.sign, dbl.storage.exponent, dbl.storage.mantissa); printf("number = %16.15g sign = %X, exponent = %X, mantissa = %lX\n", dbl.number, dbl.storage.sign, dbl.storage.exponent, dbl.storage.mantissa); return 0; }

In struct we have specified layout and union is helper so that we can see number in human readable representation. We compile it using GCC -g switch and we start GDB session to examine how layout looks like. Before that we execute program once and save output so that we can recognize our number, the second printf is of interest.

number = 3.14159265358979 sign = 0, exponent = 400, mantissa = 921FB54442D18

Here is debug session, rather most interesting parts:

(gdb) print &dbl $1 = (Helper *) 0x7fffffffe700 (gdb) x/2x 0x7fffffffe700 0x7fffffffe700: 0x54442d18 0x400921fb

After variable is initialized we queried address of it and then we asked for two units to be printed in hex format. Just to be on the safe side, let us print them 32bit by 32bit and 4 bytes by 4 bytes as well

(gdb) x/x 0x7fffffffe704 0x7fffffffe704: 0x400921fb (gdb) x/4b 0x7fffffffe704 0x7fffffffe704: 0xfb 0x21 0x09 0x40 (gdb) x/x 0x7fffffffe700 0x7fffffffe700: 0x54442d18 (gdb) x/4b 0x7fffffffe700 0x7fffffffe700: 0x18 0x2d 0x44 0x54

Where we clearly have little-endian layout in memory, least significant byte at lowest address. The same number in register

(gdb) print/x $rax $4 = 0x400921fb54442d18 (gdb) print/x $eax $5 = 0x54442d18 (gdb) print/x $ax $6 = 0x2d18 (gdb) print/x $al $7 = 0x18

Here we do not really have address so no endianness either. Maybe we could just use f format?

(gdb) x/f 0x7fffffffe700 0x7fffffffe700: 3.1415926535897931

OpenMP, file tree walk and map-reduce example

All happens on Linux and we will use GCC and built in OpenMP support and also thread local storage. We will see how easy is to spread workload across threads using OpenMP.
In man pages we find complete example for nftw or file tree walk. That example will print to standard output about everything what it finds while walking directories. Instead of having it printing everything I made it count files:

#define _XOPEN_SOURCE 500 #include <ftw.h> #include <stdio.h> #include <stdlib.h> #include <string.h> #include <stdint.h> int counter = 0; int count_files(const char *fpath, const struct stat *sb, int tflag, struct FTW *ftwbuf) { if(tflag == FTW_F){ counter++; } return 0; /* To tell nftw() to continue */ } int main(int argc, char *argv[]) { if (nftw((argc < 2) ? "." : argv[1], count_files, 20, 0) == -1) { perror("nftw"); exit(EXIT_FAILURE); } printf("There was %d files in %s\n", counter, (argc > 1)?argv[1]:"."); exit(EXIT_SUCCESS); }

To compile it I used:

gcc -O1 -o simpleftw simpleftw.c

In order to count files across multiple directories we will pass few directories in argument list to program and also turn our static variable counter into thread local storage. We will use pragma OpenMP parallel for directive combined with reduction which is in our case simple addition:

#define _XOPEN_SOURCE 500 #include <ftw.h> #include <stdio.h> #include <stdlib.h> #include <string.h> #include <stdint.h> __thread int counter = 0; int count_files(const char *fpath, const struct stat *sb, int tflag, struct FTW *ftwbuf) { if(tflag == FTW_F){ counter++; } return 0; /* To tell nftw() to continue */ } int main(int argc, char *argv[]) { int total_fc = 0; #pragma omp parallel for reduction(+:total_fc) for (int i=1; i< argc; i++){ if (nftw((argc < 2) ? "." : argv[i], count_files, 20, 0) == -1) { perror("nftw"); exit(EXIT_FAILURE); } total_fc += counter; printf("Thread %d %d\n", i, counter); } printf("There was %d files in ",total_fc); for (int i=1; i< argc; i++) printf("%s ", argv[i]); printf("\n"); exit(EXIT_SUCCESS); }

That print inside for loop should demonstrate that thread local storage works as expected. C99 switch is required, or we will have to specify that variable i is private. Also switch for OpenMP is required. To compile it I used:

gcc -O1 -fopenmp -o dwalker dwalker.c -std=c99

Execution produces on my box:

./dwalker ~/Documents/ ~/Downloads/
Thread 1 3894
Thread 2 7784
There was 11678 files in ...

Everything works as expected.

C pointers, assembler and GDB

Idea here is to do some multitasking, find out how C pointer looks like in assembler and use Linux system calls, GCC and GDB. We start with Hello World:

int main() { char *c = "Hello World!\n"; write(1, c, 13); return 0; }

Execute man 2 write in terminal to find out more about write, 1 is standard output and 12 is number of bytes to be written. I saved it as hello.c and this is how it was compiled:

gcc -O0 -m64 -save-temps -o hlo hello.c

GCC will emit assembly code as temp, I like 64 bit version and no optimizations. Content of temp file is:

.file "hello.c" .section .rodata .LC0: .string "Hello World!\n" .text .globl main .type main, @function main: .LFB0: .cfi_startproc pushq %rbp .cfi_def_cfa_offset 16 .cfi_offset 6, -16 movq %rsp, %rbp .cfi_def_cfa_register 6 subq $16, %rsp movq $.LC0, -8(%rbp) movq -8(%rbp), %rax movl $13, %edx movq %rax, %rsi movl $1, %edi movl $0, %eax call write movl $0, %eax leave .cfi_def_cfa 7, 8 ret .cfi_endproc .LFE0: .size main, .-main .ident "GCC: (Debian 4.7.2-5) 4.7.2" .section .note.GNU-stack,"",@progbits

It pushes base pointer onto stack, there is q at the end of the push and bp is called rbp, so it must be 64 bit code. Loads string into memory, later into rax, prepares call to write and so on. Now we take that temp file and compile it with -g switch so that we can use GDB.

gcc -g -O0 -o hlo hello.s
gdb hlo

We type in start at prompt and here is the GDB session:

Temporary breakpoint 1, main () at hello.s:16 16 subq $16, %rsp (gdb) n 17 movq $.LC0, -8(%rbp) (gdb) i r rbp rbp 0x7fffffffe360 0x7fffffffe360 (gdb) x 0x7fffffffe358 0x7fffffffe358: 0x00000000 (gdb) n 18 movq -8(%rbp), %rax (gdb) x 0x7fffffffe358 0x7fffffffe358: 0x004005ec (gdb) x/13c 0x004005ec 0x4005ec: 72 'H' 101 'e' 108 'l' 108 'l' 111 'o' 32 ' ' 87 'W' 111 'o' 0x4005f4: 114 'r' 108 'l' 100 'd' 33 '!' 10 '\n'

Label .LCO is our string, -8(%rbp) is where pointer is. After info registers in short form, I also restricted info to rbp, I find out that rbp points to  0x7fffffffe360 and at  0x7fffffffe360 - 8 = 0x7fffffffe358 is our pointer. No, 0x60 - 8 is not 0x52. After we send next instruction and line 17 is executed we can see our string in memory, slash 13 c is we want 13 characters from that address.
Since write is Linux system call, we can use strace to see what is going on.

strace -tt ./hlo 08:13:04.407694 execve("./hlo", ["./hlo"], [/* 37 vars */]) = 0 ... 08:13:04.412892 write(1, "Hello World!\n", 13Hello World! ) = 13 08:13:04.413108 exit_group(0) = ?

After quite few lines were replaced with three consecutive dots we see that we managed to write all 13 characters to standard output.
Using system calls is so easy that one may even attempt writing to real file:

#include <sys/types.h> #include <sys/stat.h> #include <fcntl.h> int main() { char *c = "Hello World!\n"; int fd; fd = open("deleteme.blog", O_WRONLY | O_CREAT | O_TRUNC | O_APPEND, S_IRUSR | S_IWUSR); if (fd == -1){ write(1, "Open failed!", 12); return 1; } if(write(fd, c, 13)!=13) write(1, "Write failed!", 13); if(close(fd)==-1){ write(1, "Close failed!", 13); return 1; } return 0; }

We added more system calls. We applied common flags and modes in open, checked for all possible errors, after writing string to file we have closed file. If we compile it using

gcc -O0 -m64 -save-temps -o hello2 hello.c

We have binary and also assembly code. With help of strace

strace -tt ./hello2 09:24:49.906151 execve("./hello2", ["./hello2"], [/* 37 vars */]) = 0 ... 09:24:49.911601 open("deleteme.blog", O_WRONLY|O_CREAT|O_TRUNC|O_APPEND, 0600) = 3 09:24:49.911842 write(3, "Hello World!\n", 13) = 13 09:24:49.912018 close(3) = 0 09:24:49.912386 exit_group(0) = ?

Executing cat deleteme.blog we see that write really worked and strace didn't cheat. Examining assembly code will be interesting to those who are curious how if works on low level. Finally every project manager will take a note that from 08:13:04.407694 to 09:24:49.906151 I have written only 21 line of code.

Extending PostgreSQL with own C functions

Power of open source, if you are not happy with what PostgreSQL currently offers, you write own extension in C. Compile code with your functions into shared library, install it and they will become available from PostgreSQL. OK there are some rules and procedures. Once we are inside PostgreSQL we are using its types, interfaces and utilities. Let us do Hello World example. To build example you will need postgresql-server-dev-9.1 or whatever version you are using, installed on your Linux box.

#include "postgres.h" #include <string.h> #include "fmgr.h" #include "utils/builtins.h" #ifdef PG_MODULE_MAGIC PG_MODULE_MAGIC; #endif PG_FUNCTION_INFO_V1(say_hello_n); Datum say_hello_n(PG_FUNCTION_ARGS) { text *arg1 = cstring_to_text("Hello "); text *arg2 = PG_GETARG_TEXT_P(0); int32 new_text_size = VARSIZE(arg1) + VARSIZE(arg2) - VARHDRSZ; text *new_text = (text *) palloc(new_text_size); SET_VARSIZE(new_text, new_text_size); memcpy(VARDATA(new_text), VARDATA(arg1), VARSIZE(arg1) - VARHDRSZ); memcpy(VARDATA(new_text) + (VARSIZE(arg1) - VARHDRSZ), VARDATA(arg2), VARSIZE(arg2) - VARHDRSZ); PG_RETURN_TEXT_P(new_text); }

As we can see we are using new V1 spec. That is modified concat_text example from PostgreSQL documentation, look for 35.9.4. Version 1 Calling Conventions. We have load of different VAR macros. For example:

VARHDRSZ is the same as sizeof(int4), but it's considered good style to use the macro VARHDRSZ to refer to the size of the overhead for a variable-length type.

That is from PostgreSQL documentation. SET_VARSIZE we can find in includes postgres.h

#define SET_VARSIZE_4B(PTR,len) \ (((varattrib_4b *) (PTR))->va_4byte.va_header = (((uint32) (len)) << 2))

Unless you are on big endian. Going through header file one also can read in comments more about Datum and varlena datatypes. Then we got palloc which corresponds to malloc, memcpy you already know and GET and RETURN macros. It is obvious that for writing extensions one needs to familiarize himself with PostgreSQL internals. Power without knowledge and responsibility exists only in fery tails told by “software evangelists” at annual developers developers developers meetings.
Variables passed around by PostgreSQL may be on the disc, do not change them.
To build shared library I used the following Makefile:

MODULES = say_hello PG_CONFIG = pg_config PGXS := $(shell $(PG_CONFIG) --pgxs) include $(PGXS)

Rather long story to get location of pgxs. That pgxs is location of makefiles for building extensions. It is not trivial build and using provided mk files is right way to do it. After that we can copy say_hello.so to some reasonable location or give full path to it in create function declaration.

CREATE FUNCTION say_hello(text) RETURNS text AS '/opt/pg/ext/say_hello', 'say_hello_n' LANGUAGE C STRICT; ------- select say_hello('Filip');

PostgreSQL already allows Python through untrusted language PL/Python. One can utilize power of Python for functions or triggers without learning much about PostgreSQL internals. But again if you need power and speed, you can use what PostgreSQL speaks internally and that is C.

SQLite parameterized query in C

Still very angry at Discovery Health but that is not reason to stop using C. To execute parameterized query we prepare SQL statement with one or more placeholders. Placeholder could be question mark, alone or followed by number, column, dollar sign or at sign followed by alphanumeric. For example:

select title, full_name from jane where id = @id;

We prepare such statement using sqlite3_prepare_v2, later we bind our parameter and finally execute query. To do binding we will use appropriate function, there is few of them:

int sqlite3_bind_blob(sqlite3_stmt*, int, const void*, int n, void(*)(void*)); int sqlite3_bind_double(sqlite3_stmt*, int, double); int sqlite3_bind_int(sqlite3_stmt*, int, int); int sqlite3_bind_int64(sqlite3_stmt*, int, sqlite3_int64); int sqlite3_bind_null(sqlite3_stmt*, int); int sqlite3_bind_text(sqlite3_stmt*, int, const char*, int n, void(*)(void*)); int sqlite3_bind_text16(sqlite3_stmt*, int, const void*, int, void(*)(void*)); int sqlite3_bind_value(sqlite3_stmt*, int, const sqlite3_value*); int sqlite3_bind_zeroblob(sqlite3_stmt*, int, int n);

All bind functions will return SQLITE_OK if binding is successful or error code if it fails. The first argument is handle to prepared statement. The second argument is index of parameter. The third argument is value to be set. For blob, text we have the fourth argument, size in bytes and the fifth – destructor function. Instead of destructor function we can pass constants SQLITE_STATIC - do nothing or SQLITE_TRANSIENT – make local copy and free it when done. To find out what is index of parameter we are using this function:

int sqlite3_bind_parameter_index(sqlite3_stmt*, const char *zName);

We pass prepared statement and parameter name it returns zero if there is no such parameter or parameter index if it exists. Even if we know that our parameter must have index one, we will still look for it to demonstrate how it is done. Here is the code:

#include <sqlite3.h> #include <stdio.h> #include <string.h> int main(int argc, char **argv) { int rc, i, v, col; sqlite3 *db; sqlite3_stmt *stmt; char *sql; rc = sqlite3_open_v2("surogat.db", &db, SQLITE_OPEN_READWRITE, NULL); if(rc != SQLITE_OK) { printf("Can't open database: %s\n", sqlite3_errmsg(db)); sqlite3_close(db); return 1; } sql = "select title, full_name from jane where id = @id"; rc = sqlite3_prepare_v2(db, sql, strlen(sql), &stmt, 0); if(rc != SQLITE_OK) { printf("SQL error: %s\n", sqlite3_errmsg(db)); sqlite3_close(db); return 1; } i = sqlite3_bind_parameter_index(stmt, "@id"); v = 3; rc = sqlite3_bind_int(stmt, i, v); if(rc != SQLITE_OK) { printf("SQL error: %s\n", sqlite3_errmsg(db)); sqlite3_finalize(stmt); sqlite3_close(db); return 1; } rc = sqlite3_step(stmt); col = sqlite3_column_count(stmt); for(i=0; i < col; i++) { printf("%s\t",sqlite3_column_name(stmt, i)); } printf("\n-----------------\n"); while(rc == SQLITE_ROW) { for(i=0; i < col; i++) { printf("%s\t", sqlite3_column_text(stmt, i)); } printf("\n"); rc = sqlite3_step(stmt); } sqlite3_finalize(stmt); sqlite3_close(db); return 0; }

Database should be loaded with required values in previous examples, if not here is sql to create it:

CREATE TABLE jane ( id INTEGER PRIMARY KEY NOT NULL, title TEXT, full_name TEXT NOT NULL )
INSERT INTO jane VALUES(1,'Mr','Smith');
INSERT INTO jane VALUES(2,'Mrs','Doe');
INSERT INTO jane VALUES(3,'Mr','Doe');

After we build it using:

gcc -g test.c -lsqlite3 -o test

We execute test and see the following output:

title full_name ----------------- Mr Doe

We could also misspell parameter name and rebuild to check is error handling working.

PostgreSQL, libpqxx and prepared statement example

As I promised, if people are finding interesting introductory article,  I will write more about PostgreSQL and libpqxx. But before I start with programming – rant.
I am still looking for work and finding none. Today went to Discovery Health, had very pleasant 45 minutes chat with their architect and BA and promise that we will talk again. Later guy from employment agency, who arranged meeting, calls and says that they do not want me since I do not have ANSI C in CV?! Like somebody from management decided to override decision of interviewers. If ANSI C in CV is precondition, why they wasted my time and invited me for interview? I wish them the same from their customers. BTW I am usually employed as senior developer and not as C developer or, Perl developer.
Back to programming, this time C++, had enough of ANSI C for today. Environment is Linux, I am using the same PostgreSQL and libpqxx as last time and to compile example we will use:

g++ hello_prep.cxx -o hello_prep -I/usr/local/include/ -lpqxx -lpq

This time I used test092.cxx, run dpkg with -L switch on libpqxx3-doc to see where is it. It tests passing binary parameter to prepared statement. Test macros are replaced with printing of tested values to standard output and setup of connection/transaction is included. Here is the code:

#include <iostream> #include <pqxx/pqxx> using namespace std; using namespace pqxx; int main(int argc, char** argv) { connection C("dbname=testpgdb user=testuser password=testpass hostaddr=127.0.0.1 port=5432"); if (C.is_open()) { cout << "We are connected to" << C.dbname() << endl; } else { cout << "We are not connected!" << endl; return 0; } work T(C); const char databuf[] = "Test\0data"; const string data(databuf, sizeof(databuf)); cout << "data.size() = " << data.size() << endl; cout << "databuf size = " << sizeof(databuf) << endl; const string Table = "pqxxbin", Field = "binfield", Stat = "nully"; T.exec("CREATE TEMP TABLE " + Table + " (" + Field + " BYTEA)"); if (!T.conn().supports(connection_base::cap_prepared_statements)) { cout << "Backend version does not support prepared statements. Skipping." << endl; return 0; } T.conn().prepare(Stat, "INSERT INTO " + Table + " VALUES ($1)") ("BYTEA", pqxx::prepare::treat_binary); T.prepared(Stat)(data).exec(); const result L( T.exec("SELECT length("+Field+") FROM " + Table) ); cout << "data.size() = " << data.size() << endl; cout << "L[0][0].as<size_t>() = " << L[0][0].as<size_t>() << endl; const result R( T.exec("SELECT " + Field + " FROM " + Table) ); const binarystring roundtrip(R[0][0]); cout << "roundtrip.str() = " << roundtrip.str() << endl; cout << "data = " << data << endl; cout << "roundtrip.size() = " << roundtrip.size() << endl; cout << "data.size() = " << data.size() << endl; T.exec("CREATE TEMP TABLE tuple (one INTEGER, two VARCHAR)"); pqxx::prepare::declaration d( T.conn().prepare("maketuple", "INSERT INTO tuple VALUES ($1, $2)") ); d("INTEGER", pqxx::prepare::treat_direct); d("VARCHAR", pqxx::prepare::treat_string); pqxx::prepare::invocation i( T.prepared("maketuple") ); const string f = "frobnalicious"; i(6); i(f); i.exec(); const result t( T.exec("SELECT * FROM tuple") ); cout << "t.size()= " << t.size() << endl; cout << "sz = " << 1u << endl; cout << "t[0][0].as<string>() = " << t[0][0].as<string>() << endl; cout << "l = " << 6 << endl; cout << "t[0][1].c_str() = " << t[0][1].c_str() << endl; cout << "f = " << f << endl; return 0; }

After connection is successfully obtained, transaction T is constructed. Temp table is created and after confirmation that prepared statements are supported, testing goes on. We have prepare::declaration and prepare::invocation, available in reference generated with doxygen. Adding parameters is in iterative fashion and feels natural, as they say in documentation like varargs from C. Library is well designed and easy to use, documentation, tutorial and tests are supplied. Lengths and contents should match and test succeeds.