/*
    MathArray - Stores and operates on arrays of values and numbers
    
    Copyright (C) 1995, Adam Fedor
    
    $Id: MathArray.m,v 1.18 1996/08/16 23:19:54 adam Exp $
*/
#ifdef NEXT_FOUNDATION
#import <foundation/NSArray.h>
#import <foundation/NSString.h>
#import <foundation/NSCoder.h>
#import <foundation/NSUtilities.h>
#import "MathArray/NSInvocation.h"
#else
#import <Foundation/NSArray.h>
#import <Foundation/NSString.h>
#import <Foundation/NSCoder.h>
#import <Foundation/NSUtilities.h>
#import <Foundation/NSInvocation.h>
#endif
#include "MathArrayPrivate.h"
#include "MathArray/MAValueData.h"
#include "MathArray/MAValue.h"
#include "MathArray/array_encoding.h"
#include <signal.h>
#include <errno.h>

#ifdef GNU_FOUNDATION
#  import <stdio.h>
#endif

#define DEFAULT_OBJC_TYPE @encode(float)

/* Set to non-zero for some types of math errors */
int ma_fpe_errno = 0;

/* Class variables */
// holds the list of subclasses
static NSMutableArray* numberTypes = nil;
// maximize precision
static BOOL max_precision = NO;
// handling signals and traps
static int   sig_installed = 0;
static int   sig_old_mask;
/*
static struct sigvec sig_math;
static struct sigvec sig_old;
*/
static ExceptionMask trap_vec[10];
// Exceptions
NSString *MAFloatingPointException = @"MathArray FloatingPointException";
NSString *MAParameterException = @"MathArray ParameterException";
NSString *MAArrayMismatchException = @"MathArray ArrayMismatchException";
NSString *MARangeException = @"MathArray RangeException";
NSString *MACastException = @"MathArray CastException";
NSString *MAPerformException = @"MathArray PerformException";

/* Our singal handler */
static void array_sig_handler(int sig);

/* This is a dummy function so that we can force the methods in the
   MathArrayExtensions category to be loaded by non-NeXT linkers. */
extern const char * MathArrayExtensions();
extern const char * MathArrayNSObjectExtra();

@interface NSObject (Transmuting)
- (Class)transmuteClassToSibling:(Class)aClassObject;
@end

@interface MathArray (SubclassResponsibility)
- convertFromObjCType:(const char *)old_type;
- _maOperate:(ma_operator_t)operator with:(MathArray *)value;
- _maPerform:(double_func_t)mathFunction;
- _maPerformFunction:(perform_func_t)perform_func userInfo:(void *)info;
@end

/* We have to do this because GNU gcc doesn't load in unused classes by
   default */
@interface MathUCharArray : MathArray
@end
@interface MathCharArray : MathArray
@end
@interface MathUShortArray : MathArray
@end
@interface MathShortArray : MathArray
@end
@interface MathUIntArray : MathArray
@end
@interface MathIntArray : MathArray
@end
@interface MathULongArray : MathArray
@end
@interface MathLongArray : MathArray
@end
@interface MathULongLongArray : MathArray
@end
@interface MathLongLongArray : MathArray
@end
@interface MathFloatArray : MathArray
@end
@interface MathDoubleArray : MathArray
@end
@interface MathComplexFloatArray : MathArray
@end
@interface MathComplexDoubleArray : MathArray
@end

@implementation MathArray

#ifdef NEXT_FOUNDATION
+ initialize
#else
+ (void)initialize
#endif
{
    /* Register our known subclasses */
    if (self == [MathArray class]) {
    	[self registerArraySubclass:[MathUCharArray class]];
    	[self registerArraySubclass:[MathCharArray class]];
    	[self registerArraySubclass:[MathUShortArray class]];
    	[self registerArraySubclass:[MathShortArray class]];
    	[self registerArraySubclass:[MathUIntArray class]];
    	[self registerArraySubclass:[MathIntArray class]];
    	[self registerArraySubclass:[MathULongArray class]];
    	[self registerArraySubclass:[MathLongArray class]];
    	[self registerArraySubclass:[MathULongLongArray class]];
    	[self registerArraySubclass:[MathLongLongArray class]];
    	[self registerArraySubclass:[MathFloatArray class]];
    	[self registerArraySubclass:[MathDoubleArray class]];
    	[self registerArraySubclass:[MathComplexFloatArray class]];
    	[self registerArraySubclass:[MathComplexDoubleArray class]];
		
    }
    
    [self setTrap:FLT_IMP action:RaiseException];
    [self setTrap:!FLT_IMP action:LogException];
#ifdef NEXT_FOUNDATION
    return self;
#endif
}

+ (precision_t)precision
{
#ifdef GNU_FOUNDATION
    fprintf(stderr, "Warning (MathArray): Requesting precision on untyped array\n");
#else
    NSLog(@"Warning (MathArray): Requesting precision on untyped array");
#endif
    return 0;
}

+ (const char *)objCType
{
#ifdef GNU_FOUNDATION
    fprintf(stderr, "Warning (MathArray): Requesting objCType on untyped array\n");
#else
    NSLog(@"Warning (MathArray): Requesting objCType on untyped array");
#endif
    return NULL;
}

+ (void)registerArraySubclass:mathArraySubclass
{
    BOOL added;
    int i, count;
    
    if (mathArraySubclass == nil)
    	return;
#if 0
    /* You cannot use this to check for where a class object is in the
       heirarchy. Returns YES only for isKindOf:NSObject . FIXME */
    if (![mathArraySubclass isKindOf:[MathArray class]])
    	return;
#endif
    if (numberTypes == nil)
    	numberTypes = [[NSMutableArray alloc] init];
    added = NO;
    count = [numberTypes count];
    for (i=0; i < count; i++) {
    	Class array;
        array = [numberTypes objectAtIndex:i];
    	// FIXME: Check for precision being equal?
    	if ([array precision] > [mathArraySubclass precision]) {
	    [numberTypes insertObject:mathArraySubclass atIndex:i];
	    added = YES;
	    break;
	}
    }
    if (!added)
	[numberTypes addObject:mathArraySubclass];
}

/* If flag is yes, then when MathArray is forced to promote an array to float (complex) type, it will promote it instead to double (complex double) type. This increases the precision of the calculations at the expense of doubling the size of the array. */
+ (void)promoteToDouble:(BOOL)flag
{
    max_precision = flag;
}

+ (BOOL)willPromoteToDouble
{
    return max_precision;
}

+ (void)setTrap:(unsigned)mathTrap action:(ExceptionMask)action
{
    if (mathTrap & FLT_INT_DIV0) trap_vec[0] = action;
    if (mathTrap & FLT_INT_OVER) trap_vec[1] = action;
    if (mathTrap & FLT_INEXACT) trap_vec[2] = action;
    if (mathTrap & FLT_DIV0) trap_vec[3] = action;
    if (mathTrap & FLT_UNDERFLOW) trap_vec[4] = action;
    if (mathTrap & FLT_OPERAND) trap_vec[5] = action;
    if (mathTrap & FLT_OVERFLOW) trap_vec[6] = action;
    if (mathTrap & FLT_OTHER) trap_vec[7] = action;
}


/* Returns the subclass of MathArray (as registered in registerArraySubclass:)
   that stores the Objective-C type 'theType' */
+ (Class)classForObjCType:(const char *)theType
{
    NSEnumerator* enumerator;
    Class array;
    NSString *wantType = [NSString stringWithCString:theType];
    
    if (numberTypes == nil)
    	return Nil;
    enumerator = [numberTypes objectEnumerator];
    while ((array = [enumerator nextObject])) {
	NSString *subType;
        subType = [NSString stringWithCString:[array objCType]];
    	if ([wantType isEqual:subType]) {
	    return array;
	}
    }
    
    return Nil;
}

+ (Class)classForPrecision:(precision_t)precision
{
    NSEnumerator* enumerator;
    Class array;
    
    if (numberTypes == nil)
    	return Nil;
    enumerator = [numberTypes objectEnumerator];
    while ((array = [enumerator nextObject])) {
    	if ([array precision] == precision)
	    return array;
    }
    
    return Nil;
}

/* Tries to pick the type with the best precicion.  Currently it picks the type with the higher precision, unless there is a conflict between signed and unsigned types.  In this case, it picks the signed version of the higher precision type. */
+ (const char *)resolvedTypecast:(const char *)firstType 
	:(const char *)secondType
{
    Class class1;
    Class class2;
    precision_t prec1, prec2;
    const char *rType;
    
    if ((class1 = [self classForObjCType:firstType]) == Nil) {
    	MA_RAISE_CAST;
	/* NOT REACHED */
    }
    if ((class2 = [self classForObjCType:secondType]) == Nil) {
    	MA_RAISE_CAST;
	/* NOT REACHED */
    }
    prec1 = [class1 precision];
    prec2 = [class2 precision];
    rType = (prec1 > prec2) ? [class1 objCType] : [class2 objCType];
    if (MATH_SIGNED(prec1) != MATH_SIGNED(prec2)) {
        if (prec1 > prec2 && MATH_SIGNED(prec1) == 0) {
	    prec1 |= 0x0100;
	    rType = [[self classForPrecision:prec1] objCType];
	} else if (prec1 < prec2 && MATH_SIGNED(prec2) == 0) {
	    prec2 |= 0x0100;
	    rType = [[self classForPrecision:prec2] objCType];
	}
	if (rType == NULL) {
	    MA_RAISE_CAST;
	    /* NOT REACHED */
	}
    }  

    return rType;
}

/* Creates and returns a MathArray representing a scalar value.  This 
   is a special case of an array with a dimension of zero (but a size 
   of 1). */
+ (MathArray *)maWithScalar:value
{
    MathArray *array;
    unsigned count = 1;

    if (!value) {
    	[NSException raise:MAParameterException
	    format:@"Cannot create an array with a nil value"];
	/* NOT REACHED */
    }
    if (![value respondsToSelector:@selector(objCType)]) {
    	[NSException raise:MAParameterException
		format:@"Creating array with non-value"];
	/* NOT REACHED */
    }
    array = [[self alloc] initArrayOfDimension:0
    		size:&count
		objCType:[value objCType]
		zero:NO];
    [value getValue:[[array mathData] mutableBytes]];
    return [array autorelease];
}

+ (MathArray *)maWithVector:(unsigned)count objCType:(const char *)theType 
{
    MathArray *array;
    array = [[self alloc] initArrayOfDimension:1
    		size:&count
		objCType:theType
		zero:YES];
    return [array autorelease];
}

+ (MathArray *)maMatrixWithCols:(unsigned)cols rows:(unsigned)rows objCType:(const char *)theType 
{
    MathArray *array;
    unsigned count[2];
    count[1] = cols;
    count[0] = rows;
    array = [[self alloc] initArrayOfDimension:2
    		size:count
		objCType:theType
		zero:YES];
    return [array autorelease];
}

+ (MathArray *)maWithValueData:(MAValueData *)valueData;
{
    unsigned count;
    if (!valueData) {
    	[NSException raise:MAParameterException
		format:@"Cannot create array from nil data"];
	/* NOT REACHED */
    }
    
    count = [valueData count];
    return [[[self alloc] initArrayFrom:valueData
    		ofDimension:1 
		size:&count 
		objCType:[valueData objCType] ] autorelease];
}

/* Updates the array to reflect all the currently known information about
   it. For instance, if we were initialized lazily (i.e with no data or type),
   choose a default type and/or actually allocate the data so that someone
   else can access the information. */
- (void)_updateLazyArray
{
    if (arrayData == nil) {
        unsigned total = array_num_elements(dimension, [size bytes]);
	if ([self isMemberOfClass:[MathArray class]])
    	    [self castToObjCType:DEFAULT_OBJC_TYPE];
    	arrayData = [[MAMutableValueData dataWithCount:total
	    			objCType:[self objCType]] retain];
	math_flags.zero = 0;
    }
}

/* Install and de-install signal handlers.  I guess we have to do this
   each time we go into a method, just in case some other class raises
   a fpe, we don't want to catch it. */
+ (void)_startMath
{
    /* Check if handler is already installed, and increment install count */
    if (sig_installed++)
    	return;
    
/*  FIXME: Need more universal signal handling.
    sig_math.sv_handler = array_sig_handler;
    sig_math.sv_mask    = 0;
    sig_math.sv_flags   = 0;
    sig_old_mask = sigblock(0);
    sigsetmask(sig_old_mask & ~(sigmask(SIGFPE)));
    if (sigvec(SIGFPE, &sig_math, &sig_old)) {
  	[NSException raise:NSGenericException
	    format:@"Internal error: cannot install signal handler"];
    }
*/
    errno = 0;
    ma_fpe_errno = 0;
}

+ (void)trapException: (ma_trap_t)trap reason: (NSString *)reason
{
    ExceptionMask mask;
    
    switch(trap) {
    case FLT_INT_DIV0: mask = trap_vec[0]; break;
    case FLT_INT_OVER: mask = trap_vec[1]; break;
    case FLT_INEXACT: mask = trap_vec[2]; break;
    case FLT_DIV0: mask = trap_vec[3]; break;
    case FLT_UNDERFLOW: mask = trap_vec[4]; break;
    case FLT_OPERAND: mask = trap_vec[5]; break;
    case FLT_OVERFLOW: mask = trap_vec[6]; break;
    case FLT_OTHER: mask = trap_vec[7]; break;
    default:
    	[NSException raise:NSGenericException
    		format:@"Internal Error: bad trap type"];
	return;
    }
    [MaskedException setMask:mask forException:MAFloatingPointException];
    [MaskedException raise:MAFloatingPointException
    	format:reason];
}

+ (void)_finishMath
{
    /* Check for handler and decrement install count*/
    NSAssert(sig_installed, @"Internal error: no handler installed");
    sig_installed--;

/*
    sigsetmask(sig_old_mask);
    if (sigvec(SIGFPE, &sig_old, NULL)) {
    	[NSException raise:NSGenericException
	    format:@"Internal error: cannot deinstall signal handler"];
    }
*/
    
    /* FIXME: need better idea of what caused the fpe */
    if (errno == EDOM) {
     	[self trapException:FLT_OTHER 
		reason:@"A domain error occured during processing"];
    } else if (errno == ERANGE) {
     	[self trapException:FLT_OTHER 
		reason:@"A range error occured during processing"];
    }
    if (ma_fpe_errno & FLT_DIV0) {
     	[self trapException:FLT_DIV0 
		reason:@"A divide by 0 error occured during processing"];
    }
    if (ma_fpe_errno & FLT_OVERFLOW) {
     	[self trapException:FLT_OVERFLOW 
		reason:@"An overflow error occured during processing"];
    }
    if (ma_fpe_errno & FLT_OTHER) {
     	[self trapException:FLT_OTHER 
		reason:@"A fpe error occured during processing"];
    }
    errno = 0;
    ma_fpe_errno = 0;
}

- initArrayFrom:(NSData *)data ofDimension:(unsigned)numDimensions size:(const unsigned *)sizes objCType:(const char *)type 
{
    const char *category_name;
    self = [super init];

    category_name = MathArrayExtensions();
    category_name = MathArrayNSObjectExtra();
    dimension = numDimensions;
    size = [[MAValueData dataWithValues:sizes
    		count:numDimensions
		objCType:@encode(unsigned)] retain];
    /* Notice we cast ourselves before we set data, so the data won't be
       converted (i.e. its already supposed to be this type) */
    if (data != nil && type == NULL) {
	[NSException raise:MAParameterException
	    format:@"Init array with untyped data: need to specify objCType"];
    }
    if (type != NULL) {
        [self castToObjCType:type];
    }
    if (data) {
	if ([data length] != array_num_bytes(dimension, sizes, type)) {
	    MA_RAISE_PARAMETER;
	    /* NOT REACHED */
	}
	arrayData = [[MAMutableValueData alloc] initWithData:data
			objCType:type];
    }
    math_flags.promote = YES;
    return self;
}

- initArrayOfDimension:(unsigned)numDimensions size:(const unsigned *)sizes objCType:(const char *)type zero:(BOOL)doZero 
{
    self =  [self initArrayFrom:nil
    			ofDimension:numDimensions
			size:sizes
			objCType:type];
    math_flags.zero = doZero;
    return self;
}

- (void)dealloc
{
    [arrayData release];
    [size release];
    [super dealloc];
}

// NSCopying
- deepen
{
    if (arrayData)
	arrayData = [arrayData mutableCopyWithZone:[self zone]];
    size = [size copyWithZone:[self zone]];
    return self;
}

- copyWithZone:(NSZone *)zone
{
    return [(MathArray *)NSCopyObject(self, 0, zone) deepen];
}

- copy
{
    return [self copyWithZone:[self zone]];
}

// NSCoding
// FIXME: do I need to implement classForCoder to encode a MathArray cluster?
- (void)encodeWithCoder:(NSCoder *)coder
{
    unsigned b_promote = math_flags.promote;
    
    [self _updateLazyArray];
    [coder encodeObject:arrayData]; 
    [coder encodeObject:size]; 
    [coder encodeValuesOfObjCTypes:"II", &dimension, &b_promote];
}

- (id)initWithCoder:(NSCoder *)coder
{
    unsigned b_promote = math_flags.promote;
    
    arrayData = [[coder decodeObject] retain]; 
    size = [[coder decodeObject] retain]; 
    [coder decodeValuesOfObjCTypes:"II", &dimension, &b_promote];
    math_flags.promote = b_promote;
    return self;
}

- (void)promoteIfNeeded:(BOOL)doPromote
{
    math_flags.promote = doPromote;
}

/* Returns YES if the array is promotable. */
- (BOOL)isPromotable
{
    return math_flags.promote;
}

/* Forces the array to be cast to the specified type. Precision may be 
   lost if the array is cast to a type with less precision (sort of 
   obvious, huh?). This method works even if the array is not promotable. */
- castToObjCType:(const char *)new_type
{
    Class wantClass;
    NSString *oldType = nil;

    if (!new_type)
	return nil;
    wantClass = [[self class] classForObjCType:new_type];
    if (wantClass == Nil) {
	 MA_RAISE_CAST;
	/* NOT REACHED */
    }
    /* Convert ourselves and our instance variables to the new type */
    [self transmuteClassToSibling:wantClass];
    if (arrayData && [arrayData objCType])
	oldType = [NSString stringWithCString:[arrayData objCType]];
    [arrayData setObjCType:new_type];

    /* Now we can tell our "new" self to convert the data from the old type */
    if (arrayData && oldType)
    	[self convertFromObjCType:[oldType cString]];
    return self;
}

/* Returns the dimension of the array. */
- (unsigned)dimension
{
    return dimension;
}

/* Returns an array which contains the sizes of each dimension of the 
   receiver. */
- (const unsigned *)sizes
{
    return (const unsigned *)[size bytes];
}

/* Returns a unique number which specifies the precision of numbers 
   stored by the receiver. */
- (precision_t)precision
{
    return [[self class] precision];
}

/* Returns the current type of number the receiver represents as an
   Objective-C type (obtained from @encode()). */
- (const char *)objCType
{
    return [[self class] objCType];
}

/* Returns the object which contains the numbers stored in the 
   receiver. */
- (MAMutableValueData *)mathData
{
    [self _updateLazyArray];
    return arrayData;
}

- (MathArray *)arraySubrange:(NSRange *)arrayRange
{
    int i;
    unsigned want_elements;
    unsigned *index;
    NSRange  fromRange, toRange;
    MAMutableValueData  *indexData;
    MAMutableValueData  *newData;
    MAMutableValueData  *newSizes;
    MathArray 	 *newArray;
    unsigned *new_sizes;
    const unsigned *sizes;
    [self _updateLazyArray];
    
    newSizes = [[size mutableCopy] autorelease];
    sizes = (const unsigned *)[size bytes];
    new_sizes = (unsigned *)[newSizes mutableBytes];
    /* Check to see the specified ranges fit within our limits. A range
       that isn't specified (i.e. length = 0) is the same as specifying 
       the maximum range for that dimension */
    want_elements = 1;
    for (i=0; i < dimension; i++) {
    	if (arrayRange[i].length == 0) {
	    arrayRange[i].length = sizes[i];
	    arrayRange[i].location = 0;
	}
    	if (NSMaxRange(arrayRange[i]) > sizes[i]) {
	    MA_RAISE_RANGE;
	    /* NOT REACHED */
	}
	want_elements *= arrayRange[i].length;
	new_sizes[i] = arrayRange[i].length;
    }
    
    /* Now copy the numbers into a new array */
    newData = [MAMutableValueData dataWithCount:want_elements 
    			objCType:[arrayData objCType]];
    indexData = [MAMutableValueData dataWithCount:dimension
    			objCType:@encode(unsigned)];
    index = (unsigned *)[indexData mutableBytes];
    index = start_index_from_range(dimension, arrayRange, index);
    fromRange.location = ordered_index(dimension, size, index);
    fromRange.location *= array_aligned_sizeof_elements([arrayData objCType]);
    fromRange.length   = arrayRange[dimension-1].length;
    toRange = fromRange;
    toRange.location = 0;
    i = 0;
    while (i == 0) {
	[newData replaceValuesInRange:toRange 
    		withValues:([arrayData bytes] + fromRange.location)];
	i = increment_index_in_range(dimension, arrayRange, 
			index, arrayRange[dimension-1].length);
	fromRange.location = ordered_index(dimension, size, index);
	toRange.location += toRange.length;
    }
    
    newArray = [[MathArray alloc] initArrayFrom:newData
    			ofDimension:dimension
			size:new_sizes
			objCType:[arrayData objCType]];
    return [newArray autorelease];
}

/* The array 'arrayLocations' is a MathArray where
   the elements of the array specify locations to pick out values from our
   array. These values get inserted into another array that is exactly like
   'arrayLocations' except it has the same type as us (self). This can
   give strange results, especially if 'arrayLocations' is itself an
   n-dimensional array... */
- (MathArray *)arrayValues:(MathArray *)arrayLocations
{
    int i;
    int count;
    unsigned bpe;
    unsigned long max_elements;
    MathArray *locations;
    MAMutableValueData *values;
    ordered_index_t *loc_data;
    [self _updateLazyArray];
    
    bpe = array_aligned_sizeof_elements([arrayData objCType]);
    locations = arrayLocations;
    if (strcmp([arrayLocations objCType], @encode(ordered_index_t)) != 0) {
	locations = [[arrayLocations copy] autorelease];
	[locations castToObjCType:@encode(ordered_index_t)];
    }
    loc_data = (ordered_index_t *)[[locations mathData] bytes];
    count = array_num_elements([locations dimension], [locations sizes]);
    max_elements = array_num_elements(dimension, [size bytes]);
    values = [MAMutableValueData dataWithCount:count
		objCType:[arrayData objCType]];
    for (i=0; i < count; i++) {
    	memcpy(([values mutableBytes]+i*bpe), 
		([arrayData bytes]+loc_data[i]*bpe), bpe);
    }
    
    return [[[MathArray alloc] initArrayFrom:values
		ofDimension:[locations dimension]
		size:[locations sizes]
		objCType:[arrayData objCType]] autorelease];
}

- (id <NSNumber,ComplexNumber>)arrayValueAtIndex:(unsigned *)index
{
    ordered_index_t ordered;
    const unsigned *sizes;
    
    // FIXME: need to test if index is valid (correctly)
    sizes = (const unsigned *)[size bytes];
    ordered = ordered_index(dimension, size, index);
    if (ordered >= array_num_elements(dimension, sizes)) {
    	MA_RAISE_RANGE
	/* NOT REACHED */
    }
    [self _updateLazyArray];
    return (id <NSNumber,ComplexNumber>)[arrayData valueAtIndex:ordered];
}

- (void)setArray:(MathArray *)otherArray atIndex:(unsigned *)startIndex
{
    int i;
    unsigned other_dim;
    const unsigned *sizes;
    const unsigned *other_sizes;
    unsigned *other_index;
    MAValueData *other_size;
    NSRange  *trange, *frange;
    NSRange  toRange, fromRange;
    if (!otherArray)
    	return;
    if (!startIndex)
    	[NSException raise:MAParameterException
	    format:@"No index specified"];
	
    [self _updateLazyArray];
    other_dim = [otherArray dimension];
    if (other_dim > dimension)
    	[NSException raise:MAArrayMismatchException
	    format:@"Trying to insert larger dimension array into smaller"];
    sizes = (const unsigned *)[size bytes];
    other_sizes = [otherArray sizes];
    other_size = [MAValueData dataWithValues:other_sizes 
			count:other_dim objCType:@encode(unsigned)];
    for (i=other_dim-1; i >= 0; i--)
    	if ((startIndex[i+(dimension-other_dim)]+other_sizes[i]) 
		> sizes[i+(dimension-other_dim)])
    	    [NSException raise:MAArrayMismatchException
		format:@"Trying to insert larger size array into smaller"];

    // FIXME: Correct? Should cast other array to our type
    otherArray = [[otherArray copy] autorelease];
    [otherArray castToObjCType:[self objCType]];
    
    other_index = (unsigned *)[[[size mutableCopy] autorelease] mutableBytes];
    trange = [[MAMutableValueData dataWithCount:dimension
		objCType:@encode(NSRange)] mutableBytes];
    frange = [[MAMutableValueData dataWithCount:dimension
		objCType:@encode(NSRange)] mutableBytes];
    for (i=other_dim-1; i >= 0;  i--) {
	other_index[i] = 0;
	trange[i] = NSMakeRange(startIndex[i], other_sizes[i]);
	frange[i] = NSMakeRange(0, other_sizes[i]);
    }

    /* Insert array */
    do {
	toRange.location = ordered_index(dimension, size, startIndex);
	toRange.length   = other_sizes[other_dim-1];
	fromRange.location = ordered_index(other_dim, other_size, other_index);
	fromRange.location *= array_aligned_sizeof_elements([self objCType]);
	[arrayData replaceValuesInRange:toRange
    	    withValues:([[otherArray mathData] bytes] + fromRange.location)];
	(void)increment_index_in_range(dimension, trange, 
			startIndex, trange[dimension-1].length);
	i = increment_index_in_range(other_dim, frange, 
			other_index, frange[other_dim-1].length);
    } while (i == 0);
    
}

- (void)setValues:(MathArray *)otherArray 
	atLocations:(MathArray *)arrayLocations
{
    [self _updateLazyArray];
    [self maNotImplemented:_cmd];
}

- (void)setValue:value atIndex:(unsigned *)index
{
    MathArray *scalar;
    ordered_index_t ordered;
    const unsigned *sizes;
    
    // FIXME: need to test if index is valid (correctly)
    sizes = (const unsigned *)[size bytes];
    ordered = ordered_index(dimension, size, index);
    if (ordered >= array_num_elements(dimension, sizes)) {
    	MA_RAISE_RANGE
	/* NOT REACHED */
    }
    [self _updateLazyArray];
    scalar = [MathArray maWithScalar:value];
    [scalar castToObjCType:[self objCType]];
    [arrayData setValue:[scalar arrayValueAtIndex:0] atIndex:ordered];
}

- concatArray:(MathArray *)otherArray
{
    int i;
    unsigned *new_sizes;
    const unsigned *sizes = (const unsigned *)[size bytes];
    const unsigned *other_sizes = [otherArray sizes];
    [self _updateLazyArray];
    
    /* FIXME: Should this be an exception? */
    if (!otherArray)
    	return self;
    
    /* All but the highest dimension must match, although we are allowed to
       concat a scalar onto a one-dimensional array (and vice versa) */
    if (dimension > 1 && [otherArray dimension] > 1
    		&& dimension != [otherArray dimension]) {
	[NSException raise:MAArrayMismatchException
		format:@"Dimensions must be same to concat array"];
	/* NOT REACHED */
    }
    for (i=1; i < dimension; i++) {
	if (sizes[i] != other_sizes[i]) {
	    [NSException raise:MAArrayMismatchException
		format:@"Array sizes must be same to concat array"];
	    /* NOT REACHED */
	}
    }
    if (dimension == 0)
    	dimension = 1;

    /* We have to be the same type */
    if ([self precision] < [otherArray precision] && math_flags.promote == 1) {
	[self castToObjCType: [[self class] 
		resolvedTypecast:[self objCType] :[otherArray objCType]]];
    } else if ([self precision] > [otherArray precision] 
    	    || ([self precision] > [otherArray precision] 
		&& math_flags.promote == 1)) {
	otherArray = MA_TEMPORARY(otherArray);
	[otherArray castToObjCType: [[self class] 
		resolvedTypecast:[self objCType] :[otherArray objCType]]];
    }
    
    /* Append the data and change our size */
    [arrayData appendValueData:[otherArray mathData]];
    new_sizes = [[[size mutableCopy] autorelease] mutableBytes];
    new_sizes[0] = sizes[0] + other_sizes[0];
    [size release];
    size = [[MAValueData dataWithValues:new_sizes
    		count:dimension
		objCType:@encode(unsigned)] retain];
    return self;
}

- reformArrayToDimension:(unsigned)newDimension size:(unsigned *)newSizes 
{
    if (newDimension != 0 && array_num_elements(dimension, [size bytes]) 
    		!= array_num_elements(newDimension, newSizes)) {
	MA_RAISE_PARAMETER;
	/* NOT REACHED */
    }
    if (newDimension == 0) {
	int i;
	const unsigned *sizes = (const unsigned *)[size bytes];
	if (!newSizes)
	    newSizes = [[[size mutableCopy] autorelease] mutableBytes];
	for (i=0; i < dimension; i++) {
	    if (sizes[i] > 1)
		newSizes[newDimension++] = sizes[i];
	}
    }
    dimension = newDimension;
    [size release];
    size = [[MAValueData dataWithValues:newSizes
    		count:newDimension
		objCType:@encode(unsigned)] retain];
    return self;
}

/* We got here because our current subclass doesn't handle the message,
   so we look for a class (with greater precision) that does, and recast
   ourselves to that class */
- (void)forwardInvocation:(NSInvocation *)anInvocation
{
    NSEnumerator* enumerator;
    Class array;
    
    if (numberTypes == nil) {
    	MA_RAISE_PERFORM;
	/* NOT REACHED */
    }
    enumerator = [numberTypes objectEnumerator];
    while ((array = [enumerator nextObject])) {
    	if ([self precision] < [array precision]
		&& [array instancesRespondToSelector:[anInvocation selector]]) {
	    [self castToObjCType:[array objCType]];
	    [anInvocation setTarget:self];
	    [anInvocation invoke];
	    return;
	}
    }
    
    /* We raise an exception, rather than cause an error, when we don't find
       a class to perform the message */
    MA_RAISE_PERFORM;
}

/* Performs the logical or mathematical operation (specified by operator)on the reciever using value, which may be an NSValue or a MathArray. If value is an NSValue, then value used on each element of the receiver array.  If value is a MathArray, then element 1 of value  is operated with element 1 of the receiver, and so on. Raises an MASizeException or MADimensionException if the sizes or the the dimensions of the two arrays don't match. Raises an MAFloatingPointException if a floating-point exception occurred during the processing. */
- maOperate:(ma_operator_t)operator with:value 
{
    [self _updateLazyArray];
    if ([value isKindOfClass:[MathArray class]]) {
    	// FIXME: check dimensions
    } else {
	/* A scalar is a special case of MathArray where dimension==0  */
	value = [[self class] maWithScalar:value];
    }
    // FIXME: will resolvedTypecast return an error?
    if ([self precision] < [value precision] && math_flags.promote == 1) {
	[self castToObjCType: [[self class] 
		resolvedTypecast:[arrayData objCType] :[value objCType]]];
    }
    
    [MathArray _startMath];
    [self _maOperate:operator with:(MathArray *)value];
    [MathArray _finishMath];
    
    /* Recast to int (at most) for logical operators FIXME: ok? */
    if ((operator >= ma_not) &&
    	    [self precision] > [[[self class] classForObjCType:@encode(int)]
    		precision]
	    && math_flags.promote == 1) {
	[self castToObjCType:@encode(int)];
    }

    return self;
}

/* Performs the function mathFunction on each element of the receiver. The receiver is automatically promoted to at least a float if necessary (and if possible). Raises an MAFloatingPointException if a floating-point exception occurs during the operation. */
- maPerform:(double (*)(double))mathFunction 
{
    [self _updateLazyArray];
    if ([self precision] < [[[self class] classForObjCType:@encode(float)]
    		precision]
	    && math_flags.promote == 1) {
    	if (max_precision)
	    [self castToObjCType:@encode(double)];
	else
	    [self castToObjCType:@encode(float)];
    }
    
    [MathArray _startMath];
    [self _maPerform:mathFunction];
    [MathArray _finishMath];
    return self;
}

/* For each element of the receiver, calls the function perform_func. The
   functions is passed a number which contains the value of the current
   element, an array containing the current index of the element, and a
   pointer to user supplied values given in info.  Raises an 
   MAFloatingPointException if a floating-point exception occurs 
   during the operation. */
- maPerformFunction:(perform_func_t)perform_func userInfo:(void *)info
{
    [self _updateLazyArray];
    
    [MathArray _startMath];
    [self _maPerformFunction:perform_func userInfo:info];
    [MathArray _finishMath];
    return self;
}

@end

@implementation NSObject (Transmuting)

/* This is a scary method. Basically the whole concept of casting in MathArray
   hinges on being able to switch subclasses on the fly (without having to
   allocating a new class and copying variables - that would change the
   pointer). Call it runtime casting.  It is imperative that the new
   class have the same instance variables as the old one.  In essence, they
   should all be subclasses of one superclass that defines all the instance
   variables.  This really needs better checking for this, though. */
- (Class)transmuteClassToSibling:(Class)aClassObject
{
    if (aClassObject != Nil)
        if ([self isKindOfClass:[aClassObject superclass]])
          {
            Class old_isa = isa;
            isa = aClassObject;
            return old_isa;
          }
    /* We raise an exception, rather than cause an error */
    MA_RAISE_CAST;
    /* NOT REACHED */
    return Nil;
}

@end

/* Signal handler - try to construct a decent handler
   FIXME - make a separate class to handle signals? */
void 
array_sig_handler(int sig)
{
    [MaskedException raise:MAFloatingPointException
    		format:@"A floating-point error occured during processing"];
    /* NOT REACHED */
}

