ObjFW  Artifact [9473f8f25c]

Artifact 9473f8f25c7616cf792a5a56629a423a63b8bcedcd32331fc35c63b29c7ab8f2:

• File src/OFDoubleMatrix.m — part of check-in [e1e7ffa903] at 2011-09-22 23:25:42 on branch trunk — Exceptions are now autoreleased.

  This is safe as an "exception loop" can't happen, since if allocating
  an exception fails, it throws an OFAllocFailedException which is
  preallocated and can always be thrown.

  So, the worst case would be that an autorelease of an exception fails,
  triggering an OFOutOfMemoryException for which there is no memory,
  resulting in an OFAllocFailedException to be thrown. (user: js, size: 10175) [annotate] [blame] [check-ins using]

/*
 * Copyright (c) 2008, 2009, 2010, 2011
 *   Jonathan Schleifer <js@webkeks.org>
 *
 * All rights reserved.
 *
 * This file is part of ObjFW. It may be distributed under the terms of the
 * Q Public License 1.0, which can be found in the file LICENSE.QPL included in
 * the packaging of this file.
 *
 * Alternatively, it may be distributed under the terms of the GNU General
 * Public License, either version 2 or 3, which can be found in the file
 * LICENSE.GPLv2 or LICENSE.GPLv3 respectively included in the packaging of this
 * file.
 */

#include "config.h"

#include <stdlib.h>
#include <string.h>
#include <math.h>

#import "OFDoubleMatrix.h"
#import "OFDoubleVector.h"
#import "OFString.h"

#import "OFInvalidArgumentException.h"
#import "OFNotImplementedException.h"
#import "OFOutOfMemoryException.h"
#import "OFOutOfRangeException.h"

#import "macros.h"

static Class doubleVector = Nil;

@implementation OFDoubleMatrix
+ (void)initialize
{
	if (self == [OFDoubleMatrix class])
		doubleVector = [OFDoubleVector class];
}

+ matrixWithRows: (size_t)rows
	 columns: (size_t)columns
{
	return [[[self alloc] initWithRows: rows
				   columns: columns] autorelease];
}

+ matrixWithRows: (size_t)rows
	 columns: (size_t)columns
	    data: (double)data, ...
{
	id ret;
	va_list arguments;

	va_start(arguments, data);
	ret = [[[self alloc] initWithRows: rows
				  columns: columns
				     data: data
				arguments: arguments] autorelease];
	va_end(arguments);

	return ret;
}

- init
{
	Class c = isa;
	[self release];
	@throw [OFNotImplementedException exceptionWithClass: c
						    selector: _cmd];
}

- initWithRows: (size_t)rows_
       columns: (size_t)columns_
{
	self = [super init];

	@try {
		rows = rows_;
		columns = columns_;

		if (SIZE_MAX / rows < columns ||
		    SIZE_MAX / rows * columns < sizeof(double))
			@throw [OFOutOfRangeException
			    exceptionWithClass: isa];

		if ((data = malloc(rows * columns * sizeof(double))) == NULL)
			@throw [OFOutOfMemoryException
			     exceptionWithClass: isa
				  requestedSize: rows * columns *
						 sizeof(double)];

		memset(data, 0, rows * columns * sizeof(double));

		if (rows == columns) {
			size_t i;

			for (i = 0; i < rows * columns; i += rows + 1)
				data[i] = 1;
		}
	} @catch (id e) {
		[self release];
		@throw e;
	}

	return self;
}

-   initWithRows: (size_t)rows_
	 columns: (size_t)columns_
	    data: (double)data_, ...
{
	id ret;
	va_list arguments;

	va_start(arguments, data_);
	ret = [self initWithRows: rows_
			 columns: columns_
			    data: data_
		       arguments: arguments];
	va_end(arguments);

	return ret;
}

- initWithRows: (size_t)rows_
       columns: (size_t)columns_
	  data: (double)data_
     arguments: (va_list)arguments
{
	self = [super init];

	@try {
		size_t i;

		rows = rows_;
		columns = columns_;

		if (SIZE_MAX / rows < columns ||
		    SIZE_MAX / rows * columns < sizeof(double))
			@throw [OFOutOfRangeException exceptionWithClass: isa];

		if ((data = malloc(rows * columns * sizeof(double))) == NULL)
			@throw [OFOutOfMemoryException
			     exceptionWithClass: isa
				  requestedSize: rows * columns *
						 sizeof(double)];

		for (i = 0; i < rows; i++) {
			size_t j;

			for (j = i; j < rows * columns; j += rows)
				data[j] = (j == 0
				    ? data_ : va_arg(arguments, double));
		}
	} @catch (id e) {
		[self release];
		@throw e;
	}

	return self;
}

- (void)dealloc
{
	free(data);

	[super dealloc];
}

- (void)setValue: (double)value
	  forRow: (size_t)row
	  column: (size_t)column
{
	if (row >= rows || column >= columns)
		@throw [OFOutOfRangeException exceptionWithClass: isa];

	data[row * columns + column] = value;
}

- (double)valueForRow: (size_t)row
	       column: (size_t)column
{
	if (row >= rows || column >= columns)
		@throw [OFOutOfRangeException exceptionWithClass: isa];

	return data[row * columns + column];
}

- (size_t)rows
{
	return rows;
}

- (size_t)columns
{
	return columns;
}

- (BOOL)isEqual: (id)object
{
	OFDoubleMatrix *otherMatrix;

	if (![object isKindOfClass: [OFDoubleMatrix class]])
		return NO;

	otherMatrix = object;

	if (otherMatrix->rows != rows || otherMatrix->columns != columns)
		return NO;

	if (memcmp(otherMatrix->data, data, rows * columns * sizeof(double)))
		return NO;

	return YES;
}

- (uint32_t)hash
{
	size_t i;
	uint32_t hash;

	OF_HASH_INIT(hash);

	for (i = 0; i < rows * columns; i++) {
		union {
			double d;
			uint8_t b[sizeof(double)];
		} u;
		uint8_t j;

		u.d = of_bswap_double_if_be(data[i]);

		for (j = 0; j < sizeof(double); j++)
			OF_HASH_ADD(hash, u.b[j]);
	}

	OF_HASH_FINALIZE(hash);

	return hash;
}

- copy
{
	OFDoubleMatrix *copy = [[isa alloc] initWithRows: rows
						columns: columns];

	memcpy(copy->data, data, rows * columns * sizeof(double));

	return copy;
}

- (OFString*)description
{
	OFMutableString *description;
	size_t i;

	description = [OFMutableString stringWithFormat: @"<%@, (\n",
							 [self className]];

	for (i = 0; i < rows; i++) {
		size_t j;

		[description appendString: @"\t"];

		for (j = 0; j < columns; j++) {

			if (j != columns - 1)
				[description
				    appendFormat: @"%10f ",
						  data[j * rows + i]];
			else
				[description
				    appendFormat: @"%10f\n",
						  data[j * rows + i]];
		}
	}

	[description appendString: @")>"];

	[description makeImmutable];

	return description;
}

- (double*)cArray
{
	return data;
}

- (void)addMatrix: (OFDoubleMatrix*)matrix
{
	size_t i;

	if (matrix->isa != isa || matrix->rows != rows ||
	    matrix->columns != columns)
		@throw [OFInvalidArgumentException exceptionWithClass: isa
							     selector: _cmd];

	for (i = 0; i < rows * columns; i++)
		data[i] += matrix->data[i];
}

- (void)subtractMatrix: (OFDoubleMatrix*)matrix
{
	size_t i;

	if (matrix->isa != isa || matrix->rows != rows ||
	    matrix->columns != columns)
		@throw [OFInvalidArgumentException exceptionWithClass: isa
							     selector: _cmd];

	for (i = 0; i < rows * columns; i++)
		data[i] -= matrix->data[i];
}


- (void)multiplyWithScalar: (double)scalar
{
	size_t i;

	for (i = 0; i < rows * columns; i++)
		data[i] *= scalar;
}

- (void)divideByScalar: (double)scalar
{
	size_t i;

	for (i = 0; i < rows * columns; i++)
		data[i] /= scalar;
}

- (void)multiplyWithMatrix: (OFDoubleMatrix*)matrix
{
	double *newData;
	size_t i, base1, base2;

	if (rows != matrix->columns)
		@throw [OFInvalidArgumentException exceptionWithClass: isa
							     selector: _cmd];

	if ((newData = malloc(matrix->rows * columns * sizeof(double))) == NULL)
		@throw [OFOutOfMemoryException
		     exceptionWithClass: isa
			  requestedSize: matrix->rows * columns *
					 sizeof(double)];

	base1 = 0;
	base2 = 0;

	for (i = 0; i < columns; i++) {
		size_t base3 = base2;
		size_t j;

		for (j = 0; j < matrix->rows; j++) {
			size_t base4 = j;
			size_t base5 = base1;
			double tmp = 0.0;
			size_t k;

			for (k = 0; k < matrix->columns; k++) {
				tmp += matrix->data[base4] * data[base5];
				base4 += matrix->rows;
				base5++;
			}

			newData[base3] = tmp;
			base3++;
		}

		base1 += rows;
		base2 += matrix->rows;
	}

	free(data);
	data = newData;

	rows = matrix->rows;
}

- (void)transpose
{
	double *newData;
	size_t i, k;

	if ((newData = malloc(rows * columns * sizeof(double))) == NULL)
		@throw [OFOutOfMemoryException
		    exceptionWithClass: isa
			 requestedSize: rows * columns * sizeof(double)];

	rows ^= columns;
	columns ^= rows;
	rows ^= columns;

	for (i = k = 0; i < rows; i++) {
		size_t j;

		for (j = i; j < rows * columns; j += rows)
			newData[j] = data[k++];
	}

	free(data);
	data = newData;
}

- (void)translateWithVector: (OFDoubleVector*)vector
{
	OFDoubleMatrix *translation;
	double *cArray;

	if (rows != columns || vector->isa != doubleVector ||
	    vector->dimension != rows - 1)
		@throw [OFInvalidArgumentException exceptionWithClass: isa
							     selector: _cmd];

	cArray = [vector cArray];
	translation = [[OFDoubleMatrix alloc] initWithRows: rows
						  columns: columns];

	memcpy(translation->data + (columns - 1) * rows, cArray,
	    (rows - 1) * sizeof(double));

	@try {
		[self multiplyWithMatrix: translation];
	} @finally {
		[translation release];
	}
}

- (void)rotateWithVector: (OFDoubleVector*)vector
		   angle: (double)angle
{
	OFDoubleMatrix *rotation;
	double n[3], m, angleCos, angleSin;

	if (rows != 4 || columns != 4 || vector->isa != doubleVector ||
	    vector->dimension != 3)
		@throw [OFInvalidArgumentException exceptionWithClass: isa
							     selector: _cmd];

	n[0] = vector->data[0];
	n[1] = vector->data[1];
	n[2] = vector->data[2];

	m = sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);

	if (m != 1.0) {
		n[0] /= m;
		n[1] /= m;
		n[2] /= m;
	}

	angle = (double)(angle * M_PI / 180.0);
	angleCos = cos(angle);
	angleSin = sin(angle);

	rotation = [[OFDoubleMatrix alloc] initWithRows: rows
					       columns: columns];

	rotation->data[0] = angleCos + n[0] * n[0] * (1 - angleCos);
	rotation->data[1] = n[1] * n[0] * (1 - angleCos) + n[2] * angleSin;
	rotation->data[2] = n[2] * n[0] * (1 - angleCos) - n[1] * angleSin;

	rotation->data[4] = n[0] * n[1] * (1 - angleCos) - n[2] * angleSin;
	rotation->data[5] = angleCos + n[1] * n[1] * (1 - angleCos);
	rotation->data[6] = n[2] * n[1] * (1 - angleCos) + n[0] * angleSin;

	rotation->data[8] = n[0] * n[2] * (1 - angleCos) + n[1] * angleSin;
	rotation->data[9] = n[1] * n[2] * (1 - angleCos) - n[0] * angleSin;
	rotation->data[10] = angleCos + n[2] * n[2] * (1 - angleCos);

	@try {
		[self multiplyWithMatrix: rotation];
	} @finally {
		[rotation release];
	}
}

- (void)scaleWithVector: (OFDoubleVector*)vector
{
	OFDoubleMatrix *scale;
	double *cArray;
	size_t i, j;

	if (rows != columns || vector->isa != doubleVector ||
	    vector->dimension != rows - 1)
		@throw [OFInvalidArgumentException exceptionWithClass: isa
							     selector: _cmd];

	cArray = [vector cArray];
	scale = [[OFDoubleMatrix alloc] initWithRows: rows
					    columns: columns];

	for (i = j = 0; i < ((rows - 1) * columns) - 1; i += rows + 1)
		scale->data[i] = cArray[j++];

	@try {
		[self multiplyWithMatrix: scale];
	} @finally {
		[scale release];
	}
}
@end