Serge's World

The greatest common divisor of two numbers is the largest number that divides evenly into both numbers.

The method of calculating this is is to use repeated division and subtraction.

Given two numbers a and b, the first step is to repeatedly divide both numbers by 2 as long as both are even, keeping count of the number of iterations as d.

Now, we loop again while a is not equal to b. If either of the numbers are even we divide the number by 2, without incrementing d in this case. If both numbers are odd, we then replace the bigger of a or b with the difference of the two numbers divided by 2.

Eventually the two numbers will be equal and the loop exits.

Now at this point the greatest common divisor is given by gcd = a (or b) * 2^d

public static int GreatestCommonDivisor(int a, int b)
{
	if (a <= 0 || b <= 0)
		throw new Exception("Parameters need to be greater than 0");
	int d = 0;
	while ((a % 2 == 0) && (b % 2 == 0))
	{
		a = a / 2;
		b = b / 2;
		d += 1;
	}

	while (a != b)
	{
		if (a % 2 == 0)
			a = a / 2;
		else if (b % 2 == 0)
			b = b / 2;
		else if (a > b)
			a = (a - b) / 2;
		else
			b = (b - a) / 2;
	}
	int gcd = a * ((int)Math.Pow(2, d));

	return gcd;
}

The full sourcecode for the MathLib library is available at https://github.com/sjmeunier/mathlib

Originally posted on my old blog, Smoky Cogs, on 23 Oct 2009

Updated 5 Oct 2016: Updated code snippet after refactoring MathLib library

Factorials are defined as a number that is multiplied by every positive number below it, and is denoted by an exclamation mark, ie n! = n * (n – 1) * (n – 2) … 2 * 1.

This is a very easy function to write.

public static long Factorial(int num)  
{  
    long fact = 1;  
    for (int i = 2; i <= num; i++)  
    {  
        fact = fact * i;  
    }  
    return fact;  
}

Finding the roots of a quadratic equation (which is where it crosses the x axis) can be a rather involved calculation. I remember in school having to do them, and the equations they gave usually worked out to nice numbers. Often in the real world, it is not quite so easy.

Fortunately, computers are not scared by numbers that are difficult, and can solve quadratics quite well.

The function below solves quadratics in the form of Ax² + Bx + C = 0. Providing A, B and C, the function finds the two roots of the equation, and also determines whether the roots are real, real repeating or complex roots.

public static void SolveQuadratic(double a, double b, double c, ref double x1Real, ref double x1Im, ref double x2Real, ref double x2Im, ref Enums.RootType rootType)
{
	double z = 0;

	z = Math.Pow(b, 2) - (4 * a * c);
	
	if (Math.Round(z, 2) < 0)
	{
		x1Real = -1 * b / (2.0 * a);
		x2Real = x1Real;

		x1Im = Math.Sqrt(Math.Abs(z)) / (2.0 * a);
		x2Im = -1 * x1Im;

		rootType = Enums.RootType.Complex;
	}
	else if(Math.Round(z, 2) == 0)
	{
		x1Real = -1 * b / (2.0 * a);
		x2Real = x1Real;

		x1Im = 0.0;
		x2Im = 0.0;

		rootType = Enums.RootType.Repeating;
	}
	else
	{
		x1Real = ((-1 * b) + Math.Sqrt(z)) / (2.0 * a);
		x2Real = ((-1 * b) - Math.Sqrt(z)) / (2.0 * a);

		x1Im = 0.0;
		x2Im = 0.0;

		rootType = Enums.RootType.Real;
	}
}

The full sourcecode for the MathLib library is available at https://github.com/sjmeunier/mathlib

Originally posted on my old blog, Smoky Cogs, on 23 Oct 2009

Updated 5 Oct 2016: Updated code snippet after refactoring MathLib library

Complex numbers are a vital part of many mathematical calculations, yet Visual Studio does not natively support complex numbers, but fortunately, this is a very easy oversight to fix.

First, a little maths lesson. Imaginary numbers are numbers that do not exist, like for example the square root of -1. This number is denoted by i and is central to complex numbers. Complex numbers are defined as a number that is made up of an imaginary component as well as a real component, and is denoted by the formula x*i + y.

So some examples of complex numbers are 5i + 4, 2i + 0 (which simplifies to 2i), and 0i + 6 (which simplifies to simply 6).

What the above example shows as well, is that the range of complex numbers also includes all the real numbers.

The complex number class contains two variables to contain the real and imaginary components, namely, re and im. Most of the rest of the class are overloaded arithmetic operators.

The ^ operator takes a complex number and raises it to a specified power, returning a new complex number.

The + and - operators add and subtract two complex numbers, which entails adding or subtracting the real and imaginary parts of the one to the real and imaginary parts of the other complex number.

The * and / operators are a bit more complicated in how they are calculated, but the net result is the multiplication and division of the complex numbers.

The Abs() function returns the absolute value of the complex number, which is the square root of the sum of the squares of the real and imaginary components - very much like the absolute value of a normal vector.

Lastly, the Arg() function returns the arg of the complex number in degrees, which allows the complex number to be expressed in polar notation.

namespace MathLib
{

    public class Complex
    {
        public double Real;
        public double Imaginary;

        public Complex(double real, double imaginary)
        {
            Real = real;
            Imaginary = imaginary;
        }

        public static Complex operator ^(Complex arg1, int arg2)
        {
            int i = 0;
            Complex x = new Complex(0.0, 0.0);

            if (arg2 == 0)
            {
                return x;
            }
            else
            {
                x = arg1;
                for (i = 1; i < arg2; i++)
                {
                    x = x * arg1;
                }
                return x;
            }
        }


        public static Complex operator +(Complex arg1, Complex arg2)
        {
            return (new Complex(arg1.Real + arg2.Real, arg1.Imaginary + arg2.Imaginary));
        }

        public static Complex operator -(Complex arg1)
        {
            return (new Complex(-arg1.Real, -arg1.Imaginary));
        }

        public static Complex operator -(Complex arg1, Complex arg2)
        {
            return (new Complex(arg1.Real - arg2.Real, arg1.Imaginary - arg2.Imaginary));
        }

        public static Complex operator *(Complex arg1, Complex arg2)
        {
            return (new Complex(arg1.Real * arg2.Real - arg1.Imaginary * arg2.Imaginary, arg1.Real * arg2.Imaginary + arg2.Real * arg1.Imaginary));
        }

        public static Complex operator /(Complex arg1, Complex arg2)
        {
            double c1, c2, d;
            d = arg2.Real * arg2.Real + arg2.Imaginary * arg2.Imaginary;
            if (d == 0)
            {
                return (new Complex(0, 0));
            }
            c1 = arg1.Real * arg2.Real + arg1.Imaginary * arg2.Imaginary;
            c2 = arg1.Imaginary * arg2.Real - arg1.Real * arg2.Imaginary;
            return (new Complex(c1 / d, c2 / d));
        }

        public double Abs()
        {
            return (Math.Sqrt(Real * Real + Imaginary * Imaginary));
        }

        //Arg of complex number in degrees

        public double Arg()
        {
            double ret = 0;
            if (Real != 0)
                ret = (180 / Math.PI) * Math.Atan(Imaginary / Real);
            return (ret);
           
        }

        

        public override string ToString()
        {
            return (String.Format("Complex: ({0}, {1})", Real, Imaginary));
        }

        public string ToComplexString(int iRounding)
        {
            string ComplexNumber = Convert.ToString(Math.Round(Real, iRounding)) + "+" + Convert.ToString(Math.Round(Imaginary, iRounding)) + "i";
            return ComplexNumber;
        }
    }
}

The full sourcecode for the MathLib library is available at https://github.com/sjmeunier/mathlib

Originally posted on my old blog, Smoky Cogs, on 23 Oct 2009

Updated 5 Oct 2016: Updated code snippet after refactoring MathLib library

Matrices are a very useful thing in computer programming, especially in 3D graphics programming. Out of all the mathematics I learned at University, matrices are the one thing that keep rearing their head over and over.

I am not going to go into major detail over what a matrix is - wikipedia can tell you that much better than I can - so will focus mainly on how to implement one in C#.

Essentially, a matrix is a 2-dimensional array of numbers, often representing a set of equations. The number of columns and rows is flexible.

Using the Matrix class below, the constructor can either take an x and y value or an array to determine how big to make the internal array for the matrix.

The Zero() function sets all entries in the matrix to zero by looping through the entire array.

The SetIdentity() function makes the matrix an identity matrix, which means that all the values are 0, except for the values along the diagonal from top left to bottom right being set to 1.

The Transpose function transforms the columns of the matrix into the rows, and vice versa.

The Scale function multiplies each value in the matrix by a scaling factor.

The determinant of an array is found by adding all the products of the entries of each diagonal sloping in a right direction, and subtracting the products of the entries of each diagonal sloping in a opposite direction. This is done with the Determinant function

The Trace function returns the sum of the main diagonal of the matrix.

The matrix class has also got a few overloaded operators, to make handling matrices easier.

The + operator is overloaded so that when passed two matrices, each value in the one matrix is added to the equivalent value in the second matrix. Both matrices need to be the same size for this to work.

A second overload for the + operator takes a scalar value and adds it to the matrix, which just adds that value to each value in the matrix.

The - operator performs identically to the + operator, except that it subtracts instead of adds the matrices.

The** * **operator multiplies the two matrices together, which is a little more involved than addition.
Read More »

A vector is simply a 1-dimensional array of numbers that often refer to coordinates, or directions. The theory of vectors is best left to wikipedia.

Implementing a vector class in C# is pretty easy. The vector is stored as an array within the class. The constructor can either set this array directly, or initialises it to a particular size.

The Zero() function sets the entire vector to 0. The Abs() function calculates the absolute value of the vector, which is the square root of the sum of the squares of the components of the vector. The Dot() function returns the dot product of the vector, with another vector. The dot product is the sum of the products of each component in the vectors. The Scale() function scales the vector by the scaling factor supplied. The + and - operators are overloaded to make it easy to add and subtract vectors.

Each component of first vector is added (or subtracted) to the equivalent component of the second vector, and then the resulting vector is returned.
Read More »

We have already looked at several other ways of doing a least squares fit to find an quation representing a set of data. We now look at the least squares fit using full logs, which tries to match the data with the equation y = B*x^M, using the least squares method.

As with the previous least squares functions, the function below returns the calculated values for M and B, and if no solution exists returns 0 for both of them.

public static void LeastSquaresFitLogFull(List<Point> points, ref double m, ref double b)
{
	double x1, y1, xy, x2, j, lx, ly;

	x1 = y1 = xy = x2 = 0.0;

	foreach (Point point in points)
	{
		lx = Math.Log10(point.X);
		ly = Math.Log10(point.Y);
		x1 += lx;
		y1 += ly;
		xy += ly * lx;
		x2 += lx * lx;
	}

	j = (points.Count * x2) - (x1 * x1);
	if (j != 0.0)
	{
		m = ((points.Count * xy) - (x1 * y1)) / j;
		b = ((y1 * x2) - (x1 * xy)) / j;
	}
	else
	{
		m = 0;
		b = 0;
	}
}

The full sourcecode for the MathLib library is available at https://github.com/sjmeunier/mathlib

Originally posted on my old blog, Smoky Cogs, on 23 Oct 2009

Updated 5 Oct 2016: Updated code snippet after refactoring MathLib library

This algorithm for the least squares fit uses a logarithmic abscissa, and tries to find an equation matching the data set with the form y = M * log(X)/log(10) + B, using the least squares method.

The function below finds M and B, and if no solution exists, it returns 0 for both these values.

public static void LeastSquaresFitLogAbscissa(List<Point> points, ref double m, ref double b)
{
	double x1, y1, xy, x2, j, lx;

	x1 = y1 = xy = x2 = 0.0;

	foreach (Point point in points)
	{
		lx = Math.Log10(point.X);
		x1 += lx;
		y1 += point.Y;
		xy += point.Y * lx;
		x2 += point.X * point.X;
	}

	j = (points.Count * x2) - (x1 * x1);
	if (j != 0.0)
	{
		m = ((points.Count * xy) - (x1 * y1)) / j;
		b = ((y1 * x2) - (x1 * xy)) / j;
	}
	else
	{
		m = 0;
		b = 0;
	}
}

The full sourcecode for the MathLib library is available at https://github.com/sjmeunier/mathlib

Originally posted on my old blog, Smoky Cogs, on 23 Oct 2009

Updated 5 Oct 2016: Updated code snippet after refactoring MathLib library

We have already looked at the linear least squares fit, but that does not always produces a nice fit to the data. The least squares fit with a log ordinate tries to match up the data to the equation y = B * 10^Mx using the least squares method.

In the code below, if there is no solution, the function returns M and B as 0.

public static void LeastSquaresFitLogOrdinate(List<Point> points, ref double m, ref double b)
{
	double x1, y1, xy, x2, j, ly;

	x1 = y1 = xy = x2 = 0.0;

	foreach (Point point in points)
	{
		ly = Math.Log10(point.Y);
		x1 += point.X;
		y1 += ly;
		xy += point.X * ly;
		x2 += point.X * point.X;
	}

	j = (points.Count * x2) - (x1 * x1);
	if (j != 0.0)
	{
		m = ((points.Count * xy) - (x1 * y1)) / j;
		b = ((y1 * x2) - (x1 * xy)) / j;
	}
	else
	{
		m = 0;
		b = 0;
	}
}

The full sourcecode for the MathLib library is available at https://github.com/sjmeunier/mathlib

Originally posted on my old blog, Smoky Cogs, on 23 Oct 2009

Updated 5 Oct 2016: Updated code snippet after refactoring MathLib library

When analysing data it is often useful to find an equation to show the relationship in a certain set of data. The linear least squares fit tries to do exactly that.

Given a set of data points, it tries to find a straight line with the equation y = Mx + B which best fits the data, using the least squares technique, which works out the squares of the deviations of each point and then tries to minimize those.

public static void LeastSquaresFitLinear(List<Point> points, ref double m, ref double b)
{
	double x1, y1, xy, x2, j;

	x1 = y1 = xy = x2 = 0.0;

	foreach(Point point in points)
	{
		x1 += point.X;
		y1 += point.Y;
		xy += point.X * point.Y;
		x2 += point.X * point.X;
	}

	j = (points.Count * x2) - (x1 * x1);
	if (j != 0.0)
	{
		m = ((points.Count * xy) - (x1 * y1)) / j;
		b = ((y1 * x2) - (x1 * xy)) / j;
	}
	else
	{
		m = 0;
		b = 0;
	}
}

The full sourcecode for the MathLib library is available at https://github.com/sjmeunier/mathlib

Originally posted on my old blog, Smoky Cogs, on 23 Oct 2009

Updated 5 Oct 2016: Updated code snippet after refactoring MathLib library