IONF Math Functions

This is a brief document describing the general usage of math functionality provided for the IONF software. To use these functions, simply #include "util/IONF_Maths.h" in your source file.

General functions and macros

double simpson(EvalFn f, double a, double b, int N, double parameters[])

Performs integration over a function using Simpson's rule. f is the function to be integrated over from x=a to x=b, N is the number of intervals to use in the Simpson's rule integration. N must be even.
f must be a function of the form:

double f(double x, double[] params)

where x is the point at which the function is being evaluated, and params is an array of parameters that you can pass to the function through the parameters argument to simpson. Pass NULL if the function does not require any extra parameters.

EQ(a,b)

tests a and b for equality. NOTE: This method of comparison is preferable to using the C equality operator (==) because rounding errors and misrepresentations of decimals are inevitable when using C's double (or, more commonly float) datatype. EQ is actually just a macro that tests if the difference between the two numbers is less than a specified epsilon, which is scaled accordind to the size of the numbers. therefore, when comparing something agains zero, always pass zero as the second argument.

MIN(a,b)

Macro that evaluates to the least of two values.

MAX(a,b)

Macro that evaluates to the greates of two values.

IONF_Maths.h also includes the standard header file math.h, so any functions accessable from that header file will accessable using IONF_Maths.

Vectors and Matrices

Standard C provides no support for matrices or matrix operations, and has less than optimal support for higher level datatypes such as those. Therefore, we must use several tricks to get the compiler to do what we want. I've split the information below into two sections, the first part consists of general usage and should be (more than) enough for most people, the second part provides details on the implementation and reasons why this implementation was chosen.

Creating Vectors and Matrices

Matrices are created in the following way:

MATRIX(3, 3, matrix1);

The above creates a 3 x 3 matrix named "matrix1". To access elements of the matrix, you would then write:

matrix1.d[1][2] = 3.0;

The first subscript is row number, and the second is column number.NOTE:Indices into the matrix start with zero. This means, in a 3 x 3 matrix, you will be able to access rows 0, 1, and 2 and columns 0, 1, and 2. The matrix elements would look like this:

matrix1.d[0][0]	matrix1.d[0][1]	matrix1.d[0][2]
matrix1.d[1][0]	matrix1.d[1][1]	matrix1.d[1][2]
matrix1.d[2][0]	matrix1.d[2][1]	matrix1.d[2][2]

The result of accessing an element outside of these bounds is undefined (thus, bad). Matrix elements are doubles.
You can also obtain the width and height of the matrix from the following members:

height = matrix1.h;
width = matrix1.w;

Please treat these members as read-only. Changing them will have unexpected results.

Vectors are implemented the same way as matrices but have slightly different interface. In fact, the only difference is that elemnts can be accessed using only 1 subscript. For now, only column vectors can provide this interface, row vectors must be created as a 1 x n matrix. The following example shows how to create a column vector.

CVECTOR(3, vector);
vector.d[0] = 2.0;
vector.d[2] = 7.0;
height = vector.h;

CVECT and MATRIX are both macros that declare a variable, and therefore (as mandated by the C language) must be place at the top of the function before any computation steps.

Also, because they generate a type, the dimensions of the matrix must be known at compile time. Therefore, the width and height parameters must be constants or expressions that evaluate to a constant during compile. If you don't know what size matrix you want at compile time, then you must use a VMATRIXPTR (see below)

Dynamically Allocating Matrices

These are all examples of statically allocated vectors and matrices. That means that storage is allocated for them when you use the CVECTOR or MATRIX macro. When they are created, all elements are initialized to zero. IONF_Maths also provides a facility for creating dynamic matrices. First we must create a pointer to a matrix that we can use to point to a matrix we are going to create.

MATRIXPTR(3,3, ptr);

This creates a pointer named "ptr" to a 3 x 3 matrix. But since no memory is allocated, we cannot use it yet. To allocate memory for a matrix, do the following:

NEW_MATRIX(3, 3, ptr);

This allocates memory for the new matrix and sets ptr to point to it. Once this command is carried out, memory is allocated for a 3 x 3 matrix and the h and w members are initialized. Thus, we can access elements of the matrix as follows:

ptr->d[2][1] = 42.0;
width = ptr->w;

The only problem with this method is that the pointer we create can only point to 3 x 3 matrices. If we want to create a pointer that can point to any matrix, then we have to do things a bit differently.

First, you must declare the pointer. Do this with the VMATRIXPTR macro:

VMATRIXPTR(vptr);

Once you have created the pointer, you can set it to point to a matrix the same way you would set a fixed-size matrix pointer.

NEW_MATRIX(4,4,vptr);

Here, however, is where the C language has its limitations. We can no longer access element using the subscripts. But, this isn't all that much of a problem. Instead, use the ELEM macro. It is used as follows:

double val = ELEM(vptr, 1,1);

ELEM(vptr, 2,3) = 8.0;

ELEM(vptr, 0,0) = ELEM(vptr,1,0) * ELEM(vptr, 2,0);

One danger of using dynamically allocated arrays is that you must be sure to deallocate them when finished. To do this, use

FREEMATRIX(ptr);
FREEMATRIX(vptr);

Matrices as parameters to functions

(This section is meant for people who want to write their own functions that accept matrices. Most people will want to skip ahead (hopefully). I'm assuming a bit of C knowledge here.)
Occasionally, you may want to pass matrices as parameters to functions. Because they are large and rather clunky, the preferable method is to pass a pointer into the function rather than the whole matrix itself.
So, suppose I want to pass a 3x3 matrix to a function called foo. I would declare the function as follows:

void foo(const MATRIXSTRUCT(3,3) *ptr);

The MATRIXSTRUCT macro expands to the type of an m x n matrix. Once inside the function, you can use ptr just as if you declared it using MATRIXPTR.
Additionally, you can specify a function that takes an MHEADER* to use variable-sized arrays or matrix_ptr to accept either.

Matrix and Vector functions

All functions have several things in common.

Functions take as arguments pointers to matrices or vectors. The type matrix_ptr can accept pointers to all types of matrices an vectors. Note that if the function will not alter the matrix, then the type is a const matrix_ptr. Here is an example, using Zero function:

MATRIX(3,3,m1);
CVECTOR(2,v1);
MATRIXPTR(4,4,mptr1);
VMATRIXPTR(mptr2);

NEWMATRIX(4,4,mptr1);
NEWMATRIX(2,2,mptr2);

Zero(&m1);   /* must use address-of operator (&), since   */
Zero(&c1);   /* these are not pointers                    */

Zero(mptr1); /* don't need address-of, since these are    */
Zero(mptr2); /* already pointers                          */

If a function computes a matrix, then you must pass in a matrix to store the result in. You must also make sure that the matrix you pass in is already initialized and is the proper size.
During debugging, several extra checks will be in place to notify the user when the wrong size or type of matrix has been passed to the function. In the final build, these checks will be removed and the results will be undefined. When I say the function will "fail" I mean that it will generate an assertion failure and halt the program during debug, and have undefined resultsfor the final build. So, having a function "fail" is a Very Bad Thing.
For efficiency, a number have macros have and will be defined for the most commonly used operations. Do not concern yourself with these at the moment.

The functions

double DotP(const matrix_ptr V1, const matrix_ptr V2)

Calculates the Euclidian dot product of the two vectors. Fails if they are not vectors.

void Identity(matrix_ptr M);

Makes M an identity matrix. M must be square.

void MCopy(matrix_ptr R, const matrix_ptr A)

Copies all elemnts of A into R. A nd R must have equal dimension.

void MMult(matrix_ptr R, const matrix_ptr A, const matrix_ptr B)

Multiplies two matrices/vectors and stores the result in R. if R is the wrong size or A and B cannot be multiplied, the function fails.

void Normalize(matrix_ptr V)

Normalizes the Vector V. V must be a vector (1 x n or m x 1) or the function will fail. Does not do anything if the lenth is zero.

void Scale(matrix_ptr R, double scal)

Scales the matrix (multiplies all elements) by factor scal

void SetMrowmaj(matrix_ptr M, ...)

Allows you to set elements of a matrix. To use this function, first pass the address of the matrix you wish to change. Then pass a list of doubles like this:
Suppose we have an m x n matrix, A:

A(1,1)   A(1,2)   ...   A(1,n)
A(2,1)   A(2,2)   ...   A(2,n)
...       ...     ...    ...
A(m,1)   A(m,2)   ...   A(m,n)

You would call

SetMrowmaj(&A, A(1,1), A(1,2), ... , A(2,1), A(2,2), ..., A(m,1), A(m,2), ... , A(m,n));

where A(x,y) are doubles you wish to place in the (x,y) position.
NOTE: It is important that you have exactly the correct number of arguments. Also, it is extremely important that each value is a double. That means you cannot pass in a literal "2" because the compiler treats that as an int, and must instead use "(double)2" or "2.0". (Floats would work as well thanks to the "default argument promotions"; however, no one should be using floats.)

void SetMscal(matrix_ptr M, double scal)

Sets every element of M to scal.

void printMatrix(const matrix_ptr M)

Fof debugging purposes only, writes a description of matrix M to stdout.

Some Macros

Here are some macros that are meant to simplify your life a bit. Be careful with them.

CANMULT(a,b)

Checks if two matrices can be multiplied. a and b must be matrix_ptrs.

EQDIM(a,b)

Checks if two matrices have the same dimension. a, b must be matrix_ptrs.

ISSQUARE(a)

Checks if a is a square matrix. a must be a matrix_ptr.

LEN(a)

Computes length of a. a must be a matrix_ptr pointing to a vector. (expands to

sqrt(DotP((a),(a)))

)

ONE(a)

Sets a to all ones. a must be a matrix_ptr. Calls SetMScal.

ZERO(a)

Sets a to all zeros. a must be a matrix_ptr. Calls SetMscal.

Dave Paul

Last modified: Sat Jun 30 12:53:38 EDT 2001
$Id: math.html,v 1.5 2001/06/30 17:05:41 dpaul Exp $