2011-02-25 Linux Assembler SSE add

SSe programming is whery interesting fromthat point that there are parallely 4 numbers that are porcessed.SSE has registers of size 128 bits. They can handle 4 floats.GCC C there is no default type for 128 bits and we define our ownstructure for that.

typedef struct xmm
{
    float a;
    float b;
    float c;
    float d;
} xmm __attribute__ ((aligned (16)));

movaps xmm0, [eax]
movaps xmm1, [ebx]

addps xmm0,xmm1

movaps [eax], xmm0

typedef struct xmm
{
    float a;
    float b;
    float c;
    float d;
} xmm __attribute__ ((aligned (16)));

extern void sse_add( xmm *, xmm * );

int main( int argc, char **argv)
{
    xmm x0,x1;
    x0.a = 1.0;
    x0.b = 2.0;
    x0.c = 3.0;
    x0.d = 4.0;
    x1.a = x1.b = x1.c = x1.d = 5.0;
    
    printf("%10f %10f %10f %10f\n",x0.a,x0.b,x0.c,x0.d);
    printf("%10f %10f %10f %10f\n",x1.a,x1.b,x1.c,x1.d);
    
    sse_add( &x0 , &x1 );
    
    printf("%10f %10f %10f %10f\n",x0.a,x0.b,x0.c,x0.d);
    printf("%10f %10f %10f %10f\n",x1.a,x1.b,x1.c,x1.d);
    
    return 0;
}

format ELF

section '.text'

public sse_add

align 4
sse_add:
    ;arguments that are pointers for 2 xmm data blocks
    x0 equ [ebp+8]
    x1 equ [ebp+12]
    
    push ebp
    mov ebp, esp
    
    mov eax, x0
    mov ebx, x1
    
    ;load in xmm0 and xmm1 values
    ;if values where not aligned than we would used other instruction
    movaps xmm0, [eax]
    movaps xmm1, [ebx]
    
    ;sum up and save inside xmm0
    addps xmm0,xmm1
    
    ;save value in first argument
    movaps [eax], xmm0
    
    pop ebp
    ret

2011-01-16 FPU catch division by zero

There can occure some problems in C onw of them is divison on zero. For this we setup system exception handler or signal handler.  When is division on zero it works.Also for return in main function there is used setjmp and longjmp

void set_exception_handler()
{
	int err;
	fenv.__control_word &= ~FE_ALL_EXCEPT;
	fenv.__cs_selector &= ~FE_ALL_EXCEPT << 7;
	fesetenv( &fenv );	
	
	sa.sa_sigaction = &exception_handler;
	sa.sa_flags = SA_SIGINFO;
	err = sigaction( SIGFPE, &sa, NULL );
	if (err != 0)
		printf("Cannot set FloatingPoint exception handler\n");
	else
		printf("[OK] SIGFPE is set\n");
}

void exception_handler(int i, siginfo_t *s, void *v )
{
	if (s->si_signo == SIGFPE)
	{
		printf("[SIGFPE] SIGFPE Occure\n");
		printf("[SIGFPE] Error number: %d\n", s->si_errno);
		printf("[SIGFPE] Signal code: %d\n", s->si_code);
		switch (s->si_code)
		{
			case FPE_INTDIV:
				printf("[SIGFPE] Divison by 0\n");
				longjmp( jmp , 1 );
				break;
		}
	}
	abort();
}

2010-11-09 Numerical integration. Simpson method

Simpsons method not to compilcated to calculate numerical value of integral.Main point of this Simpson rules is insert values in to parabole
\[ y(x) = Ax^2+Bx+C \]
this parabole going trought points x_0, x_0+h, x_0+2*h we doubling number of points where we will calculate by 2 and this gives us enought points for this formula.
\[ \int_{x_0}^{x_0+2h} y(x)= \frac{h}{3}(y_0+4y_1+y_2) \]

#include <stdio.h>
#include <math.h>

#define TYPE float
#define NUMBER 10

TYPE integ_simpson( TYPE f(TYPE) , TYPE, TYPE, int);
TYPE fun( TYPE );

int main()
{
	printf("Result: %f\n",integ_simpson( &fun , 0.0 , 1.0 , NUMBER ));
	return 0;
}

TYPE integ_simpson( TYPE f(TYPE) , TYPE a, TYPE b, int n)
{
	int i;
	n=2*n;
	TYPE sum=f(a),h=(b-a)/n;
	for (i=1;i<=n-1;i+=2) sum += 4*f(a+h*i);
	for (i=2;i<=n-2;i+=2) sum += 2*f(a+h*i);
	sum += f(b);
	return h * sum / 3;
}

TYPE fun( TYPE x )
{
	return 1/(1+pow(x,2.0));
}

2010-07-20 Python MIT Binary Trees

After reading chapter form MIT Introduction in algorithms about trees I have implemented same algorithm in pythonI haven't tryed to make best perfomance only easy to understund and one to one like is in pseudo-codeTree python class that are used to represent BinaryTree.

class Tree:
	p = None
	left = None
	right = None
	key = 0

Inorder_Tree_Walk( x )
if x != NIL
    then Inorder_Tree_Walk( left[x] )
        print key[x]
        Inorder_Tree_Walk( right[x] )

def inorder_tree_walk( t ):
    if t != None:
        inorder_tree_walk( t.left )
        print t.key,
        inorder_tree_walk( t.right )

Tree_Search( x , k )
if x = NIL or k = key[x]
    then return x
if k < key[x]
    then return Tree_Search( left[x] , k )
    else return Tree_Search( right[x] , k )

def tree_search( t , k ):
    if (t == None) or (k == t.key):
        return t
    if k < t.key:
        return tree_search( t.left, k )
    return tree_search( t.right, k )

Tree_Minimum( x )
while left[x] != NIL
    do x <- left[x]
return x

def tree_minimum( t ):
    while t.left != None:
        t = t.left
    return t

Tree_Maximum( x )
while right[x] != NIL
    do x <- right[x]
return x

def tree_maximum( t ):
    while t.right != None:
        t = t.right
    return t

def tree_root( t ):
    while ( t.p != None):
        t = t.p
    return t

Tree_Successor( x )
if right[x] != NIL
    then return Tree_Minimum( right[x] )
y <- p[x]
while y != NIL and x = right[y]
    do  x <- y
        y <- p[y]
return y

def tree_successor( t ):
    if t.right != None:
        return tree_minimum( t.right )
    y = t.p
    while (y != None) and (t == y.right):
        t = y
        y = y.p
    return y

Tree_Insert( T , z )
y <- NIL
x <- root[T]
while x != NIL
    do y <- x
        if key[z] < key[x]
            then x <- left[x]
            else x <- right[x]
p[x] <- y
if y = NIL
    then root[T] <- z
    else if key[z] < key[y]
        then left[y] <- z
        else right[y] <- z

def tree_insert( t , z ):
    y = None
    x = tree_root( t )
    while x != None:
        y = x
        if z.key < x.key:
            x = x.left
        else:
            x = x.right
    z.p = y
    if y == None:
        r = tree_root( t )
        r = z
    else:
        if z.key < y.key:
            y.left = z
        else:
            y.right = z

def tree_insert_recrusive( t , z ):
    if t.left == None and t.right == None:
        if z.key < t.key:
            t.left = z
        else:
            t.right = z
        return
    if z.key < t.key:
        tree_insert_recrusive( t.left , z )
    else:
        tree_insert_recrusive( t.right , z )

Tree_Delete( T , z )
if left[z] = NIL or right[z] = NIL
    then y <- z
    else y <- Tree_Successor( z )
if left[y] != NIL
    then x <- left[y]
    else x <- right[y]
if x != NIL
    then p[x] <- p[y]
if p[y] = NIL
    then root[T] <- x
    else if y = left[p[y]]
        then left[p[y]] <- x
        else right[p[y]] <- x
if y != z
    then key[z] <- key[y]
return y

def tree_delete( t , z ):
    if (z.left == None) or (z.right == None):
        y = z
    else:
        y = tree_successor( z )
    if y.left != None:
        x = y.left
    else:
        x = y.right
    if x != None:
        x.p = y.p
    if y.p == None:
        r = tree_root( t )
        r = x
        t = r
    else:
        if y == y.p.left:
            y.p.left = x
        else:
            y.p.right = x
    if y != z:
        z.key = y.key
    return y

keys = [10,6,1,0,3,8,7,9,21,15,11,17,25,23,46]
max_deep = log(len(keys),2)

def create_tree( n=0 , p=None):
    if (len(keys) == 0) or (n >= max_deep):
        return None
    t = Tree()
    t.p = p
    t.key = keys.pop(0)
    t.left = create_tree( n+1 , t )
    t.right = create_tree( n+1 , t)
    return t

def print_tree( t ):
    if (t != None) and (t.key != None):
        if t.left == t.right == None:
            print "Key:%d "%(t.key)
            return
        if t.left.key == None:
            print "Key:%d Right:%d"%(t.key,t.right.key)
            print_tree( t.right )
            return
        if t.right.key == None:
            print "Key:%d Left:%d"%(t.key,t.left.key)
            print_tree( t.left )
            return
        print "Key:%d Left:%d Right:%d"%(t.key,t.left.key,t.right.key)
        print_tree( t.left )
        print_tree( t.right )

t = create_tree()
r = tree_search( t, 10 )
n = Tree()
n.key = 150
tree_insert_recrusive( t , n )
inorder_tree_walk( t )
print ""
tree_delete( t , r )
inorder_tree_walk( t )
print ""
r = tree_root( t )
print r.key

2010-01-16 Math

Different stuff that have some math inside
2D Vector Library

Attractors

Attractor
Lorenz attractor

Fractals

Simple Fractal
Blender Fractal
Blender Towers

Numerical

Simposon integration method

Other

Voronoi
Pascal Triangle

2010-01-16 Vector Library

Library with 2D vector functions.I have tested this library with
ChipMunk vector and with
VL vectors
My implementation performs some +% in speed.
Source

2009-11-15 Pascal Triangle

This is Pascal Triangle. And colours is chosen by reminder of division. There are 2 ways how that is made:1. Simple if division reminder is 0 then red else green2. Reminder is like alpha channel of red colour.In gallery shown images are generated by dividing on (17,19,23,41,42,43). There is easy to see that prime numbers (17,19,41,43) has very structured triangle but not prime (42) is different.

./chessboard div_number -

simply show triangle divided by number

./chessboard div image.tga -

save program screen to image.tga

./chessboard div img.tga ant_parm -

changes to second colouring type. Why it I call it chessboard? I were trying make such with OpenGL bet result is Pascal Triangle.
Source

2009-10-24 Making C executables smaller

There are some simple things that can be done to make C executables as small as possible.
Here is some example code we will work with:

#include <SDL/SDL.h>

char quit = 0;

int main()
{
    SDL_Surface *screen,surface;
    SDL_Event e;
    SDL_Init( SDL_INIT_VIDEO );
    screen = SDL_SetVideoMode( 400, 400, 32, SDL_SWSURFACE );
    while(!quit)
        while(SDL_PollEvent(&e)>0)
        {
            if(e.type==SDL_MOUSEBUTTONDOWN) quit=1;
            if(e.type==SDL_KEYDOWN) quit=1;
        }
    SDL_Quit();
}

-O  -O0  -O1  -O2  -O3  -Os
-ffast-math
-fomit-frame-pointer
-fauto-inc-dec
-mpush-args
-mno-red-zone
-mstackrealign 

void _start()

asm ( \
      "movl $1,%eax\n" \
      "xor %ebx,%ebx\n" \
      "int $128\n" \
    );

a=/tmp/I;tail -n+2 $0|zcat>$a;chmod +x $a;$a;rm $a;exit

gcc -Os -ffast-math -fomit-frame-pointer 
-fauto-inc-dec -mpush-args -mno-red-zone -c small.c;
ld -dynamic-linker /lib/ld-linux.so.2 small.o /usr/lib/libSDL.so -o small;
strip -s -R .comment -R .gnu.version small;sstrip small;
7z a -tGZip -mx=9 small.gz small > /dev/null;
cat unpack.header small.gz > small;
chmod a+x small;rm small.gz small.o

2009-05-07 PyGame Attractors

I have found old article about attractors that is resoult what I have done with playing with attractor values. Attractor formula is very simple but result is interesting.