This is MathArrayExtensions.m in view mode; [Download] [Up]
/*
MathArrayExtensions - More methods for MathArray
Copyright (C) 1995, Adam Fedor
$Id: MathArrayExtensions.m,v 1.17 1996/05/01 02:50:39 adam Exp $
*/
#include <math.h>
#include <string.h>
#ifdef NEXT_FOUNDATION
#import <foundation/NSString.h>
#else
#import <Foundation/NSString.h>
#endif
#include "MathArrayPrivate.h"
#include "MathArray/MAValueData.h"
int ma_random_seed = 0;
/* This is a dummy function so that we can force the methods in this
file to be loaded by non-NeXT linkers. */
const char *
MathArrayExtensions()
{
return "MathArrayExtensions";
}
@implementation MathArray (ArrayOperators)
/* MathArray is a general class for performing mathematical operations on array (1-D vectors, 2-D matrices, etc) of values. See MathArray for a better description of the class. This document describes some convienent extensions to the 'spartan' class definition. Most of these methods are implemented in terms of other methods.
*/
/* Replaces the receiver by (receiver)^exp. */
- maExponent:value
{
return [self maOperate:ma_exponent with:value];
}
/* Multiplies the receiver by factor. */
- maMultiply:factor
{
return [self maOperate:ma_multiply with:factor];
}
/* Divides the receiver by divisor */
- maDivide:divisor
{
return [self maOperate:ma_divide with:divisor];
}
/* Adds value to the receiver. */
- maAdd:value
{
return [self maOperate:ma_add with:value];
}
/* Subtracts value from the receiver. */
- maSubtract:value
{
return [self maOperate:ma_subtract with:value];
}
/* Replaces each element of the receiver with the mathematical
MOD of the receiver and value.*/
- maMod:value
{
return [self maOperate:ma_mod with:value];
}
/* Replaces each element of the receiver by the minimum of either the
receiver or value. */
- maMinimum:value
{
return [self maOperate:ma_minimum with:value];
}
/* Replaces each element of the receiver by the maximum of either the
receiver or value. */
- maMaximum:value
{
return [self maOperate:ma_maximum with:value];
}
/* Performs a logical NOT on the receiver. */
- maNot
{
return [self maOperate:ma_not with:nil];
}
/* Logically tests for the receiver and value being equal */
- maEqual:value
{
return [self maOperate:ma_add with:value];
}
/* Logically test for the receiver and value being not equal. */
- maNotEqual:value
{
return [self maOperate:ma_not_equal with:value];
}
/* Logically test for the receiver being less than or equal to value. */
- maLessOrEqual:value
{
return [self maOperate:ma_less_or_equal with:value];
}
/* Logically test for the receiver being less than value. */
- maLess:value
{
return [self maOperate:ma_less with:value];
}
/* Logically test for the receiver being greater than or equal to value */
- maGreaterOrEqual:value
{
return [self maOperate:ma_greater_or_equal with:value];
}
/* Logically test for the receiver being greater than value. */
- maGreater:value
{
return [self maOperate:ma_greater with:value];
}
/* Takes the logical AND of value and the receiver. */
- maAnd:value
{
return [self maOperate:ma_and with:value];
}
/* Takes the logical OR of value and the receiver. */
- maOr:value
{
return [self maOperate:ma_or with:value];
}
/* Takes the logical XOR of value and the receiver. */
- maXor:value
{
return [self maOperate:ma_xor with:value];
}
/* Performs matrix or vector multiplication between the receiver and
otherArray. Raises an MAArrayMismatchException if the dimension or
sizes of the arrays are not compatible. */
- maMatrixMultiply:(MathArray *)otherArray
{
const unsigned *a_size;
const unsigned *b_size;
if (dimension > 2 || [otherArray dimension] > 2) {
[NSException raise:MAArrayMismatchException
format:@"Unable to multiply arrays of dimension > 2"];
}
if (dimension == 0 || [otherArray dimension] == 0) {
[NSException raise:MAArrayMismatchException
format:@"Request to matrix multiply scalar value."];
}
a_size = [self sizes];
b_size = [otherArray sizes];
if ((dimension == 2 && a_size[1] != b_size[0])
|| (dimension ==1 && b_size[0] != 1)) {
[NSException raise:MAArrayMismatchException
format:@"Matrix multiply with incompatible arrays."];
}
[self _updateLazyArray];
/* Subclasses should implement the actual multiplication */
return self;
}
@implementation MathArray (ArrayFunctions)
/* Takes the absolute value of the receiver. */
- maAbs
{
return [self maPerform: (double_func_t)fabs];
}
/* Takes the inverse or 'arc' cosine of the receiver. */
- maArcCos
{
return [self maPerform: (double_func_t)acos];
}
/* Takes the inverse or 'arc' sine of the receiver. */
- maArcSin
{
return [self maPerform: (double_func_t)asin];
}
/* Takes the inverse or 'arc' tangent of the receiver. */
- maArcTan
{
return [self maPerform: (double_func_t)atan];
}
/* Takes the cosine of the receiver. */
- maCos
{
return [self maPerform: (double_func_t)cos];
}
/* Replaces the receiver by exp(receiver). */
- maExp
{
return [self maPerform: (double_func_t)exp];
}
/* Take the logarithm (base 10) of the receiver. */
- maLogBase10
{
return [self maPerform: (double_func_t)log10];
}
/* Take the natural logarithm (base e) of the receiver. */
- maNaturalLog
{
return [self maPerform: (double_func_t)log];
}
/* Take the tangent of the receiver. */
- maTan
{
return [self maPerform: (double_func_t)tan];
}
/* Take the sin of the receiver. */
- maSin
{
return [self maPerform: (double_func_t)sin];
}
/* Take the square root of the receiver. */
- maSqrt
{
return [self maPerform: (double_func_t)sqrt];
}
/* Invert the array */
- maInvert
{
[self maSubclassResponsibility:_cmd];
return self;
}
/* Transpose the array. */
- (void)_doTranspose
{
MAMutableValueData *tMatrix;
MAMutableValueData *tSize;
NSRange *range;
const void *data;
void *t_data;
const unsigned *sizes;
unsigned *t_sizes;
unsigned *t_index;
unsigned max_elements;
int bpe;
int i;
/* FIXME: There must be a way to do this "in place" */
tMatrix = [arrayData mutableCopy];
tSize = [size mutableCopy];
t_sizes = (unsigned *)[tSize mutableBytes];
t_index = [[MAMutableValueData dataWithCount:dimension
objCType:@encode(unsigned)] mutableBytes];
range = [[MAMutableValueData dataWithCount:dimension
objCType:@encode(NSRange)] mutableBytes];
sizes = [self sizes];
for (i=0; i < dimension; i++) {
t_sizes[i] = sizes[dimension - i - 1];
range[i] = NSMakeRange(0, sizes[i]);
}
data = [arrayData bytes];
t_data = [tMatrix mutableBytes];
bpe = array_aligned_sizeof_elements([arrayData objCType]);
max_elements = array_num_elements(dimension, [size bytes]);
start_index_from_range(dimension, range, t_index);
for (i=0; i < max_elements; i++) {
memcpy((t_data+inverted_ordered_index(dimension, tSize, t_index)*bpe),
(data+i*bpe), bpe);
increment_index_in_range(dimension, range, t_index, 1);
}
[arrayData release];
[size release];
arrayData = tMatrix;
size = tSize;
}
- maTranspose
{
[self _updateLazyArray];
/* No real need for this - Should we enforce it anyway?
if (dimension > 2) {
[NSException raise:MAArrayMismatchException
format:@"Unable to transpose arrays of dimension > 2"];
}
*/
if (dimension == 0)
return self;
if (dimension == 1) {
unsigned new_sizes[2];
new_sizes[1] = 1;
new_sizes[2] = [(NSNumber *)[size valueAtIndex:0] unsignedIntValue];
return [self reformArrayToDimension:2 size: new_sizes];
}
if (dimension == 2
&& [(NSNumber *)[size valueAtIndex:0] unsignedIntValue] == 1) {
return [self reformArrayToDimension:0 size: NULL];
}
[self _doTranspose];
return self;
}
/* Performs a Fast-Fourier Transform on the receiver in the specified
direction: -1 for a forward transform, +1 for an inverse transform. */
- maFFt:(int)direction
{
if (![self isComplex] && math_flags.promote == 1) {
if ([[self class] willPromoteToDouble])
[self castToObjCType:@encode(complex_double)];
else
[self castToObjCType:@encode(complex_float)];
return [self maFFt:direction];
}
[NSException raise:MAPerformException
format:@"Unable to perform FFt on non-complex arrays."];
return self;
}
/* Shift each value of the array by an amount specified in the array
positions. Values wrap around in the array. */
- maShift:(int *)positions
{
int i;
int array_size;
void *data;
unsigned *index = [[[size mutableCopy] autorelease] mutableBytes];
const unsigned *sizes = [size bytes];
NSRange *range;
[self _updateLazyArray];
range = [[MAMutableValueData dataWithCount:dimension
objCType:@encode(NSRange)] mutableBytes];
for (i=0; i < dimension; i++) {
range[i].location = 0;
range[i].length = sizes[i];
}
array_size = array_num_elements(dimension, sizes);
data = [arrayData mutableBytes];
for (i = 0; i < dimension; i++) {
array_size /= sizes[i];
if (positions[i] != 0) {
int shift;
MAMutableValueData *temp;
shift = (positions[i] % sizes[i]);
if (shift < 0) shift = sizes[i] - shift;
temp = [MAMutableValueData dataWithCount:array_size*shift
objCType:[self objCType]];
index = start_index_from_range(dimension, range, index);
do {
NSRange trange;
/* Copy end of array into temp array */
index[i] = sizes[i] - shift;
trange = NSMakeRange(ordered_index(dimension, size, index),
array_size*shift);
[arrayData getValues:[temp mutableBytes] range:trange];
/* Shift rest of array over */
/* Here we assume the replace operation correctly handles
the case where the replacing buffer overlaps the replaced
buffer */
index[i] = shift;
trange = NSMakeRange(ordered_index(dimension, size, index),
array_size*(sizes[i]-shift));
index[i] = 0;
[arrayData replaceValuesInRange:trange
withValues: ((char *) data
+ordered_index(dimension,size,index)
*array_aligned_sizeof_elements([self objCType]))];
/* Now copy back temp array to beggining or the array */
trange = NSMakeRange(ordered_index(dimension, size, index),
array_size*shift);
[arrayData replaceValuesInRange:trange
withValues:[temp bytes]];
} while (!(increment_index_in_range(dimension,
range, index, array_size*sizes[i])));
} /* if */
} /* for dimension */
return self;
}
+ (void) setRandomSeed: (int)seed
{
ma_random_seed = seed;
}
- maRandomUniformMin: (double)min toMax: (double)max
{
[self maSubclassResponsibility:_cmd];
return self;
}
- maRandomNormalMean: (double)mean deviation: (double)dev
{
[self maSubclassResponsibility:_cmd];
return self;
}
/* Returns a new number which contains the minimum value of the array. */
- (id <NSNumber,ComplexNumber>)maMinimumValue
{
/* If we're not a subclass, we must be un-initialized */
return [NSNumber numberWithFloat:0.0];
}
/* Returns a new number which contains the maximum value of the array. */
- (id <NSNumber,ComplexNumber>)maMaximumValue
{
/* If we're not a subclass, we must be un-initialized */
return [NSNumber numberWithFloat:0.0];
}
/* Returns a new number which contains the sum of all the elements in
the array. */
- (id <NSNumber,ComplexNumber>)maTotal
{
/* If we're not a subclass, we must be un-initialized */
return [NSNumber numberWithFloat:0.0];
}
/* Replaces the receiver by an array which contains the orderedIndex values
of where the receiver was nonzero. You can use this in more complex
expressions to, for instance, return the location of where the receiver
is greater than 2:
number2 = [NSNumber numberWithDouble:2.0];
c = [[[MA_TEMPORARY(b) maGreater:number2] maWhere]
*/
- maWhere
{
[self maSubclassResponsibility:_cmd];
return self;
}
@end
@implementation MathArray (ComplexExtensions)
/* Sneaky method to take a real array and cast it to a complex array that
is unmeshed (imaginary numbers comes after all real numbers). After
that, the values are meshed together. Used by
maComplexArrayWithReal:imaginary: */
- _reformToComplexMeshed
{
unsigned i, elements;
unsigned bpe;
unsigned *new_sizes;
char *data;
const char *saved;
MAMutableValueData* tempData;
/* Set arrayData to nil to avoid having the cast function recast
the values in the array - since they are already the correct type. */
tempData = arrayData;
arrayData = nil;
if (strcmp([tempData objCType], @encode(float)) == 0) {
[self castToObjCType:@encode(complex_float)];
} else {
[self castToObjCType:@encode(complex_double)];
}
arrayData = tempData;
if (strcmp([tempData objCType], @encode(float)) == 0) {
[arrayData setObjCType:@encode(complex_float)];
} else {
[arrayData setObjCType:@encode(complex_double)];
}
new_sizes = [[[size mutableCopy] autorelease] mutableBytes];
new_sizes[0] = new_sizes[0]/2;
[size release];
size = [[MAValueData dataWithValues:new_sizes
count:dimension
objCType:@encode(unsigned)] retain];
/* Now mesh it -
This method uses a lot of space. There is another method which
only invloves swaping values, but this takes time because it's
basically a weird permutation problem. Which is better? */
bpe = array_aligned_sizeof_elements([arrayData objCType]) / 2;
elements = array_num_elements(dimension, [size bytes]);
saved = [[[arrayData copy] autorelease] bytes];
data = [arrayData mutableBytes];
for (i=0; i < elements; i++) {
memcpy(data+2*i*bpe, saved+i*bpe, bpe);
memcpy(data+(2*i+1)*bpe, saved+(elements+i)*bpe, bpe);
}
return self;
}
/* Returns a new complex array made from combining two non-complex
arrays specifing the real and imaginary parts of the array */
+ (MathArray *)maComplexArrayWithReal:(MathArray *)realArray
imaginary:(MathArray *)imagArray
{
MathArray *array;
unsigned long length;
const char* theType;
precision_t r_prec, i_prec;
if (!realArray && !imagArray) {
MA_RAISE_PARAMETER;
/* NOT REACHED */
}
theType = (realArray) ? [realArray objCType] : [imagArray objCType];
if ([realArray isComplex] || [imagArray isComplex]) {
MA_RAISE_ILLEGAL;
/* NOT REACHED */
}
if (realArray && imagArray) {
// FIXME: should arrays have same dimensions?
length = array_num_elements([realArray dimension], [realArray sizes]);
if (length
!= array_num_elements([imagArray dimension],
[imagArray sizes])) {
MA_RAISE_ARRAY_MISMATCH;
/* NOT REACHED */
}
r_prec = [realArray precision];
i_prec = [imagArray precision];
if (r_prec > i_prec) {
imagArray = MA_TEMPORARY(imagArray);
[imagArray castToObjCType:[realArray objCType]];
} else if (i_prec > r_prec) {
realArray = MA_TEMPORARY(realArray);
[realArray castToObjCType:[imagArray objCType]];
}
}
if (realArray)
array = [realArray copy];
else
array =[[[self class] alloc] initArrayOfDimension:[imagArray dimension]
size:[imagArray sizes]
objCType:theType
zero:YES];
if (imagArray)
[array concatArray:imagArray];
i_prec = [[[self class] classForObjCType:@encode(float)] precision];
if ([array precision] < i_prec)
[array castToObjCType:@encode(float)];
[array _reformToComplexMeshed];
return [array autorelease];
}
/* Returns YES if the receiver is a complex array */
- (BOOL)isComplex
{
if ([self objCType] == NULL)
return NO;
if (strcmp([self objCType], @encode(complex_float)) == 0
|| strcmp([self objCType], @encode(complex_double)) == 0)
return YES;
return NO;
}
/* Replaces the receiver by the real portion of itself. If the receiver is
not a complex array, it remains unchanged */
- maReal
{
[self _updateLazyArray];
if (strcmp([self objCType], @encode(complex_float)) == 0)
[self castToObjCType:@encode(float)];
else if (strcmp([self objCType], @encode(complex_double)) == 0)
[self castToObjCType:@encode(double)];
return self;
}
/* Replaces the receiver by the imaginary portion of itself. If the
receiver is not a complex array, a zero-valued array is returned */
- maImaginary
{
[self _updateLazyArray];
NSAssert(!([self isComplex]), @"Invalid call to super from complex array");
return [self maMultiply:[NSNumber numberWithInt:0]];
}
/* Replaces the receiver by the magnitude portion of itself. If the
receiver is not a complex array, it remains unchanged */
- maMagnitude
{
[self _updateLazyArray];
NSAssert(!([self isComplex]), @"Invalid call to super from complex array");
return self;
}
/* Replaces the receiver by the phase portion of itself. If the
receiver is not a complex array, a zero-valued array is returned */
- maPhase
{
[self _updateLazyArray];
NSAssert(!([self isComplex]), @"Invalid call to super from complex array");
return [self maMultiply:[NSNumber numberWithInt:0]];
}
/* Replaces the receiver by its complex conjugate. If the
receiver is not a complex array, it is converted into a complex array
with a zero-values imaginary part. */
- maConjugate
{
[self _updateLazyArray];
NSAssert(!([self isComplex]), @"Invalid call to super from complex array");
if (!math_flags.promote)
return self;
if (strcmp([self objCType], @encode(complex_double)) == 0)
[self castToObjCType:@encode(complex_double)];
else
[self castToObjCType:@encode(complex_float)];
return self;
}
@end
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.