This is MathTemplateArray.m in view mode; [Download] [Up]
/*
MathTemplateArray - Template for several number array classes
Copyright (C) 1995, Adam Fedor
Should run this through a preprocessor each time defining the
macro MATH_ORDER to the appropriate values.
$Id: MathTemplateArray.m,v 1.12 1996/09/18 15:45:09 adam Exp $
*/
# line 12 "MathTemplateArray.m"
#ifdef NEXT_FOUNDATION
#import <foundation/NSString.h>
#include <sys/param.h>
#else
#import <Foundation/NSString.h>
#endif
#include <math.h>
#include <time.h>
#include "MathArrayPrivate.h"
#include "MathArray/MAValue.h"
#include "MathArray/MAValueData.h"
#include "MathArray/array_encoding.h"
#ifdef HAVE_DRAND48_DEF
#include <stdlib.h>
#define GEN_RANDOM drand48()
#else
#if (sun && __svr4__) || defined(__hpux) || defined(_SEQUENT_) || defined(mips)
extern double drand48();
extern void srand48();
#define GEN_RANDOM drand48()
#else
extern long random();
extern void srandom();
#define srand48 srandom
#define RANDOM_MAX (2147483647) /* (2**31)-1 */
#define GEN_RANDOM (double)(random() & 0xffffffffl)/RANDOM_MAX
#endif
#endif /* not HAVE_DRAND48_DEF */
/* Note: precsion is defined in terms of three words: 0xaabbcc, where cc
is the storage size of the number. bb is 1 if signed and aa is an
arbitrary order that increases for more "precise" numbers */
#if MATH_ORDER == 1
# define MathTemplateArray MathUCharArray
# define MATH_TYPE unsigned char
# define MATH_CAST_METHOD unsignedCharValue
static precision_t precision = 0x040001;
#elif MATH_ORDER == 2
# define MathTemplateArray MathCharArray
# define MATH_TYPE char
# define MATH_CAST_METHOD charValue
static precision_t precision = 0x080102;
#elif MATH_ORDER == 3
# define MathTemplateArray MathUShortArray
# define MATH_TYPE unsigned short
# define MATH_CAST_METHOD unsignedShortValue
static precision_t precision = 0x0c0002;
#elif MATH_ORDER == 4
# define MathTemplateArray MathShortArray
# define MATH_TYPE short
# define MATH_CAST_METHOD shortValue
static precision_t precision = 0x100102;
#elif MATH_ORDER == 5
# define MathTemplateArray MathUIntArray
# define MATH_TYPE unsigned int
# define MATH_CAST_METHOD unsignedIntValue
static precision_t precision = 0x140004;
#elif MATH_ORDER == 6
# define MathTemplateArray MathIntArray
# define MATH_TYPE int
# define MATH_CAST_METHOD intValue
static precision_t precision = 0x180104;
#elif MATH_ORDER == 7
# define MathTemplateArray MathULongArray
# define MATH_TYPE unsigned long
# define MATH_CAST_METHOD unsignedLongValue
static precision_t precision = 0x1c0004;
#elif MATH_ORDER == 8
# define MathTemplateArray MathLongArray
# define MATH_TYPE long
# define MATH_CAST_METHOD longValue
static precision_t precision = 0x200104;
#elif MATH_ORDER == 9
# define MathTemplateArray MathULongLongArray
# define MATH_TYPE unsigned long long
# define MATH_CAST_METHOD unsignedLongLongValue
static precision_t precision = 0x240008;
#elif MATH_ORDER == 10
# define MathTemplateArray MathLongLongArray
# define MATH_TYPE long long
# define MATH_CAST_METHOD longLongValue
static precision_t precision = 0x280108;
#elif MATH_ORDER == 11
# define MathTemplateArray MathFloatArray
# define MATH_TYPE float
# define MATH_CAST_METHOD floatValue
static precision_t precision = 0x2c0108;
#elif MATH_ORDER == 12
# define MathTemplateArray MathDoubleArray
# define MATH_TYPE double
# define MATH_CAST_METHOD doubleValue
static precision_t precision = 0x300110;
#endif
inline static MATH_TYPE find_min(MATH_TYPE a, MATH_TYPE b);
inline static MATH_TYPE find_max(MATH_TYPE a, MATH_TYPE b);
inline static MATH_TYPE find_sum(MATH_TYPE sum, MATH_TYPE num);
typedef MATH_TYPE (*caster_f_t)(void *, unsigned long);
typedef MATH_TYPE (*operate_f_t)(MATH_TYPE, MATH_TYPE);
static caster_f_t cast_function(const char *);
static operate_f_t operate_function(ma_operator_t);
@interface MAMutableValueData(MutableData)
- (void) setLength: (unsigned)length;
@end
@interface MathTemplateArray : MathArray
{
}
@end
@implementation MathTemplateArray : MathArray
+ (precision_t)precision
{
return precision;
}
+ (const char *)objCType
{
return @encode(MATH_TYPE);
}
- convertFromObjCType:(const char *)oldType
{
ordered_index_t i, max_elements;
const char* newType;
MATH_TYPE *data;
MATH_TYPE (*caster)(void *, unsigned long) = cast_function(oldType);
newType = [self objCType];
max_elements = array_num_elements(dimension, [size bytes]);
if (array_aligned_sizeof_elements(oldType)
< array_aligned_sizeof_elements(newType)) {
[arrayData setLength: max_elements * math_aligned_size(newType)];
data = [arrayData mutableBytes];
for (i=0; i < max_elements; i++)
data[max_elements-i-1] = caster(data, max_elements-i-1);
} else {
data = [arrayData mutableBytes];
for (i=0; i < max_elements; i++)
data[i] = caster(data, i);
[arrayData setLength: max_elements * math_aligned_size(newType)];
}
[arrayData setObjCType: newType];
return self;
}
- maMatrixMultiply:(MathArray *)otherArray
{
unsigned i, j, k;
unsigned a_size[2];
unsigned b_size[2];
void *b_data;
MAMutableValueData *cMatrix;
MATH_TYPE *a_data;
MATH_TYPE *c_data;
MATH_TYPE (*caster)(void *, unsigned long)
= cast_function([otherArray objCType]);
[super maMatrixMultiply:otherArray];
a_size[0] = [self sizes][0];
if (dimension == 2)
a_size[1] = [self sizes][1];
else
a_size[1] = 1;
b_size[0] = [otherArray sizes][0];
if ([otherArray dimension] == 2)
b_size[1] = [otherArray sizes][1];
else
b_size[1] = 1;
cMatrix = [MAMutableValueData dataWithCount:a_size[0]*b_size[1]
objCType:[self objCType]];
a_data = [arrayData mutableBytes];
b_data = [[otherArray mathData] mutableBytes];
c_data = [cMatrix mutableBytes];
[MathArray _startMath];
for (i=0; i < a_size[0]; i++) {
for (j=0; j < b_size[1]; j++) {
unsigned loc = i*b_size[1]+j;
c_data[loc] = 0;
for (k=0; k < a_size[1]; k++) {
c_data[loc] += a_data[i*a_size[1]+k]
* caster(b_data, k*b_size[1]+j);
#ifdef MATH_CHECK_FPE
if (!finite(c_data[loc]))
ma_fpe_errno = FLT_OVERFLOW;
#endif
}
}
}
[arrayData release];
arrayData = [cMatrix retain];
a_size[1] = b_size[1];
if (a_size[1] == 1)
dimension = 1;
[size release];
size = [[MAValueData dataWithValues:a_size
count:dimension
objCType:@encode(unsigned)] retain];
[MathArray _finishMath];
return self;
}
- _maOperate:(ma_operator_t)operator with:(MathArray *)value
{
unsigned long i, j, max_elements;
unsigned long increment;
void *otherData;
MATH_TYPE *data;
MATH_TYPE (*operate)(MATH_TYPE, MATH_TYPE)
= operate_function(operator);
MATH_TYPE (*caster)(void *, unsigned long)
= cast_function([value objCType]);
if (operate == 0) {
MA_RAISE_ILLEGAL;
/* NOT REACHED */
}
if (operator == ma_matrix_multiply) {
return [self maMatrixMultiply:value];
}
data = [arrayData mutableBytes];
increment = 1;
otherData = [[value mathData] mutableBytes];
if ([value dimension] == 0)
increment = 0;
max_elements = array_num_elements(dimension, [size bytes]);
j = 0;
for (i=0; i < max_elements; i++) {
data[i] = operate(data[i], caster(otherData, j));
#ifdef MATH_CHECK_FPE
if (!finite(data[i]))
ma_fpe_errno = FLT_OTHER;
#endif
j += increment;
}
return self;
}
- _maPerform: (double_func_t)mathFunction
{
unsigned long i, max_elements;
MATH_TYPE *data;
max_elements = array_num_elements(dimension, [size bytes]);
data = [arrayData mutableBytes];
for (i=0; i < max_elements; i++) {
data[i] = (MATH_TYPE)mathFunction((double)data[i]);
#ifdef MATH_CHECK_FPE
if (!finite(data[i]))
ma_fpe_errno = FLT_OTHER;
#endif
}
return self;
}
- _maPerformFunction:(perform_func_t)perform_func userInfo:(void *)info
{
unsigned long i, max_elements;
unsigned *index = [[[size mutableCopy] autorelease] mutableBytes];
MATH_TYPE *data;
NSRange *range;
range = [[MAMutableValueData dataWithCount:dimension
objCType:@encode(NSRange)] mutableBytes];
for (i=0; i < dimension; i++) {
range[i].location = 0;
range[i].length = ((unsigned *)[size bytes])[i];
}
index = start_index_from_range(dimension, range, index);
max_elements = array_num_elements(dimension, [size bytes]);
data = [arrayData mutableBytes];
// FIXME: need to set up autorelease pool to get rid of extra objects
for (i=0; i < max_elements; i++) {
id <NSNumber,ComplexNumber> number;
number = (id <NSNumber,ComplexNumber>)[arrayData valueAtIndex:i];
number = perform_func(number, index, info);
data[i] = [number MATH_CAST_METHOD];
#ifdef MATH_CHECK_FPE
if (!finite(data[i]))
ma_fpe_errno = FLT_OTHER;
#endif
increment_index_in_range(dimension, range, index, 1);
}
return self;
}
- maAbs
{
unsigned long i, max_elements;
MATH_TYPE *data;
/* If we're floating, use fabs instead of abs */
if (precision >= [[[self class] classForObjCType:@encode(float)] precision])
return [self maPerform: (double_func_t)fabs];
/* if we're unsigned, what's the point? */
if (MATH_SIGNED(precision) == 0)
return self;
[self _updateLazyArray];
max_elements = array_num_elements(dimension, [size bytes]);
data = [arrayData mutableBytes];
for (i=0; i < max_elements; i++)
data[i] = abs(data[i]);
return self;
}
- maInvert
{
[self maNotImplemented:_cmd];
return self;
}
/* Generate a uniform random number n, such that min <= n < max */
- maRandomUniformMin: (double)min toMax: (double)max
{
unsigned long i, max_elements;
MATH_TYPE *data;
[self _updateLazyArray];
max_elements = array_num_elements(dimension, [size bytes]);
data = [arrayData mutableBytes];
if (ma_random_seed == 0)
ma_random_seed = time(NULL);
srand48(ma_random_seed);
for (i=0; i < max_elements; i++)
{
data[i] = GEN_RANDOM * (max - min) + min;
}
return self;
}
/* Generate a normail random number with the given mean and standard
deviation. Uses the Box-Miller method which produces two numbers
each iteration */
- maRandomNormalMean: (double)mean deviation: (double)dev
{
unsigned long i, max_elements;
MATH_TYPE *data;
[self _updateLazyArray];
max_elements = array_num_elements(dimension, [size bytes]);
data = [arrayData mutableBytes];
if (ma_random_seed == 0)
ma_random_seed = time(NULL);
srand48(ma_random_seed);
for (i=0; i < max_elements/2; i++)
{
double r, theta;
r = dev * sqrt( -2.0 * log(GEN_RANDOM) );
theta = 2.0 * M_PI * GEN_RANDOM;
data[i*2] = r * cos(theta) + mean;
data[i*2+1] = r * sin(theta) + mean;
}
if (max_elements & 0x01)
{
double r, theta;
r = dev * sqrt( -2.0 * log(GEN_RANDOM) );
theta = 2.0 * M_PI * GEN_RANDOM;
data[max_elements - 1] = r * cos(theta) + mean;
}
return self;
}
- (id <NSNumber,ComplexNumber>)_loopArrayWith:
(MATH_TYPE (*)(MATH_TYPE, MATH_TYPE))loopFunc
{
unsigned long i, max_elements;
const MATH_TYPE *data;
MATH_TYPE number;
[self _updateLazyArray];
[MathArray _startMath];
max_elements = array_num_elements(dimension, [size bytes]);
data = [arrayData bytes];
number = data[0];
for (i=1; i < max_elements; i++)
number = loopFunc(number, data[i]);
[MathArray _finishMath];
return (MANumber *)[MANumber value:&number withObjCType:@encode(MATH_TYPE)];
}
- (id <NSNumber,ComplexNumber>)maMinimumValue
{
return [self _loopArrayWith:find_min];
}
- (id <NSNumber,ComplexNumber>)maMaximumValue
{
return [self _loopArrayWith:find_max];
}
- (id <NSNumber,ComplexNumber>)maTotal
{
return [self _loopArrayWith:find_sum];
}
- maWhere
{
ordered_index_t i, max_elements;
unsigned j;
MATH_TYPE *data;
/* Make sure we can handle the ordered index values */
if ([self precision]
< [[MathArray classForObjCType:@encode(unsigned)] precision]
&& math_flags.promote == 1) {
[self castToObjCType:@encode(unsigned)];
return [self maWhere];
}
j = 0;
max_elements = array_num_elements(dimension, [size bytes]);
data = [arrayData mutableBytes];
for (i=0; i < max_elements; i++)
if (data[i])
data[j++] = i;
[arrayData setCount:j];
[size release];
dimension = 1;
size = [[MAValueData dataWithValues:&j
count:dimension
objCType:@encode(unsigned)] retain];
return self;
}
@end
inline static
MATH_TYPE find_min(MATH_TYPE a, MATH_TYPE b)
{
return MIN(a,b);
}
inline static
MATH_TYPE find_max(MATH_TYPE a, MATH_TYPE b)
{
return MAX(a,b);
}
inline static
MATH_TYPE find_sum(MATH_TYPE sum, MATH_TYPE num)
{
return sum+num;
}
/* All the functions which know how to cast from any known type to our type */
static inline MATH_TYPE
cast_from_uchar(void *data, unsigned long loc)
{return *((unsigned char*)data + loc);}
static inline MATH_TYPE
cast_from_char(void *data, unsigned long loc)
{return *((char*)data + loc);}
static inline MATH_TYPE
cast_from_ushort(void *data, unsigned long loc)
{return *((unsigned short*)data + loc);}
static inline MATH_TYPE
cast_from_short(void *data, unsigned long loc)
{return *((short*)data + loc);}
static inline MATH_TYPE
cast_from_uint(void *data, unsigned long loc)
{return *((unsigned int*)data + loc);}
static inline MATH_TYPE
cast_from_int(void *data, unsigned long loc)
{return *((int*)data + loc);}
static inline MATH_TYPE
cast_from_ulong(void *data, unsigned long loc)
{return *((unsigned long*)data + loc);}
static inline MATH_TYPE
cast_from_long(void *data, unsigned long loc)
{return *((long*)data + loc);}
static inline MATH_TYPE
cast_from_ulonglong(void *data, unsigned long loc)
{return *((unsigned long long*)data + loc);}
static inline MATH_TYPE
cast_from_longlong(void *data, unsigned long loc)
{return *((long long*)data + loc);}
static inline MATH_TYPE
cast_from_float(void *data, unsigned long loc)
{return *((float*)data + loc);}
static inline MATH_TYPE
cast_from_double(void *data, unsigned long loc)
{return *((double*)data + loc);}
static inline MATH_TYPE
cast_from_complexfloat(void *data, unsigned long loc)
{return *((float*)data + 2*loc);}
static inline MATH_TYPE
cast_from_complexdouble(void *data, unsigned long loc)
{return *((double*)data + 2*loc);}
static caster_f_t
cast_function(const char *type)
{
switch(*type) {
case _C_CHR:
return cast_from_char;
break;
case _C_UCHR:
return cast_from_uchar;
break;
case _C_SHT:
return cast_from_short;
break;
case _C_USHT:
return cast_from_ushort;
break;
case _C_INT:
return cast_from_int;
break;
case _C_UINT:
return cast_from_uint;
break;
case _C_LNG:
return cast_from_long;
break;
case _C_ULNG:
return cast_from_ulong;
break;
case 'q':
return cast_from_longlong;
break;
case 'Q':
return cast_from_ulonglong;
break;
case _C_FLT:
return cast_from_float;
break;
case _C_DBL:
return cast_from_double;
break;
case _C_STRUCT_B:
// FIXME: need better check for complex
if (math_aligned_size (type)
== math_aligned_size (@encode(complex_float)))
return cast_from_complexfloat;
else if (math_aligned_size (type)
== math_aligned_size (@encode(complex_double)))
return cast_from_complexdouble;
break;
}
MA_RAISE_CAST;
/* NOT REACHED */
return 0;
}
/* All the functions which know how to operate on this type */
// FIXME: need pow for natural numbers
static inline MATH_TYPE
op_exponent(MATH_TYPE a, MATH_TYPE b) {return pow(a,b);}
static inline MATH_TYPE
op_multiply(MATH_TYPE a, MATH_TYPE b) {return a*b;}
static inline MATH_TYPE
op_divide(MATH_TYPE a, MATH_TYPE b) {return a/b;}
static inline MATH_TYPE
op_add(MATH_TYPE a, MATH_TYPE b) {return a+b;}
static inline MATH_TYPE
op_subtract(MATH_TYPE a, MATH_TYPE b) {return a-b;}
static inline MATH_TYPE
op_mod(MATH_TYPE a, MATH_TYPE b)
{
#if MATH_ORDER == 11
return fmod(a,b);
#elif MATH_ORDER == 12
return fmod(a,b);
#else
return a % b;
#endif
}
static inline MATH_TYPE
op_minimum(MATH_TYPE a, MATH_TYPE b) {return MIN(a,b);}
static inline MATH_TYPE
op_maximum(MATH_TYPE a, MATH_TYPE b) {return MAX(a,b);}
static inline MATH_TYPE
op_not(MATH_TYPE a, MATH_TYPE b)
{
#if MATH_ORDER == 11
return (a == 0) ? 1.0 : 0.0;
#elif MATH_ORDER == 12
return (a == 0) ? 1.0 : 0.0;
#else
return ~a;
#endif
}
static inline MATH_TYPE
op_equal(MATH_TYPE a, MATH_TYPE b)
{
return (a == b);
}
static inline MATH_TYPE
op_not_equal(MATH_TYPE a, MATH_TYPE b)
{
return (a != b);
}
static inline MATH_TYPE
op_less_or_equal(MATH_TYPE a, MATH_TYPE b)
{
return (a <= b);
}
static inline MATH_TYPE
op_less(MATH_TYPE a, MATH_TYPE b)
{
return (a < b);
}
static inline MATH_TYPE
op_greater_or_equal(MATH_TYPE a, MATH_TYPE b)
{
return (a >= b);
}
static inline MATH_TYPE
op_greater(MATH_TYPE a, MATH_TYPE b)
{
return (a > b);
}
// FIXME: How do you do logical ops on floats?
static inline MATH_TYPE
op_and(MATH_TYPE a, MATH_TYPE b)
{
#if MATH_ORDER == 11
return (a != 0 && b != 0) ? b : 0.0;
#elif MATH_ORDER == 12
return (a != 0 && b != 0) ? b : 0.0;
#else
return (a & b);
#endif
}
static inline MATH_TYPE
op_or(MATH_TYPE a, MATH_TYPE b)
{
#if MATH_ORDER == 11
return (a != 0 || b != 0) ? 1.0 : 0.0;
#elif MATH_ORDER == 12
return (a != 0 || b != 0) ? 1.0 : 0.0;
#else
return (a | b);
#endif
}
static inline MATH_TYPE
op_xor(MATH_TYPE a, MATH_TYPE b)
{
#if MATH_ORDER == 11
return ((a != 0 || b != 0) && !(a != 0 && b != 0)) ? 1.0 : 0.0;
#elif MATH_ORDER == 12
return ((a != 0 || b != 0) && !(a != 0 && b != 0)) ? 1.0 : 0.0;
#else
return (a ^ b);
#endif
}
static operate_f_t
operate_function(ma_operator_t operator)
{
switch(operator) {
case ma_exponent:
return op_exponent;
break;
case ma_multiply:
return op_multiply;
break;
case ma_matrix_multiply:
NSCAssert(0, @"Invalid operator reference");
break;
case ma_divide:
return op_divide;
break;
case ma_add:
return op_add;
break;
case ma_subtract:
return op_subtract;
break;
case ma_mod:
return op_mod;
break;
case ma_minimum:
return op_minimum;
break;
case ma_maximum:
return op_maximum;
break;
case ma_not:
return op_not;
break;
case ma_equal:
return op_equal;
break;
case ma_not_equal:
return op_not_equal;
break;
case ma_less_or_equal:
return op_less_or_equal;
break;
case ma_less:
return op_less;
break;
case ma_greater_or_equal:
return op_greater_or_equal;
break;
case ma_greater:
return op_greater;
break;
case ma_and:
return op_and;
break;
case ma_or:
return op_or;
break;
case ma_xor:
return op_xor;
break;
}
return 0;
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.