This is MathDecimalArray.m in view mode; [Download] [Up]
/*
MathDecimalArray - MathArray subclass of NSDecimal numbers
Copyright (C) 1997, Adam Fedor
*/
#import <Foundation/NSString.h>
#import <Foundation/NSDecimalNumber.h>
#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 */
static precision_t precision = 0x320108;
/* A minor non-thread safe variable */
static int rounding;
inline static NSDecimal find_min(NSDecimal a, NSDecimal b);
inline static NSDecimal find_max(NSDecimal a, NSDecimal b);
inline static NSDecimal find_sum(NSDecimal sum, NSDecimal num);
typedef NSDecimal (*caster_f_t)(void *, unsigned long);
typedef NSDecimal (*operate_f_t)(NSDecimal, NSDecimal);
static caster_f_t cast_function(const char *);
static operate_f_t operate_function(ma_operator_t);
static NSDecimal
NSSetDecimal(double num)
{
NSDecimal result;
//FIXME: find a better way of setting
result = [[NSNumber numberWithDouble: num] decimalValue];
return result;
}
static ma_trap_t
trap_error_from_decimal(NSCalculationError err)
{
ma_trap_t trap;
trap = 0;
switch(err)
{
case NSCalculationLossOfPrecision: trap = FLT_INEXACT; break;
case NSCalculationUnderflow: trap = FLT_UNDERFLOW; break;
case NSCalculationOverflow: trap = FLT_OVERFLOW; break;
case NSCalculationDivideByZero: trap = FLT_DIV0; break;
default:
trap = 0;
}
return trap;
}
@interface MAMutableValueData(MutableData)
- (void) setLength: (unsigned)length;
@end
@interface MathDecimalArray : MathArray
{
}
@end
@implementation MathDecimalArray : MathArray
+ (precision_t)precision
{
return precision;
}
+ (const char *)objCType
{
return @encode(NSDecimal);
}
- convertFromObjCType:(const char *)oldType
{
ordered_index_t i, max_elements;
const char* newType;
NSDecimal *data;
NSDecimal (*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;
NSDecimal *a_data;
NSDecimal *c_data;
NSDecimal (*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] = NSSetDecimal(0);
for (k=0; k < a_size[1]; k++)
{
NSCalculationError calc_error;
NSDecimal result, temp;
temp = caster(b_data, k*b_size[1]+j);
calc_error = NSDecimalMultiply(&result, &a_data[i*a_size[1]+k],
&temp, math_flags.round);
ma_fpe_errno |= trap_error_from_decimal(calc_error);
calc_error = NSDecimalAdd(&c_data[loc], &c_data[loc], &result,
math_flags.round);
ma_fpe_errno |= trap_error_from_decimal(calc_error);
}
}
}
[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;
NSDecimal *data;
NSDecimal (*operate)(NSDecimal, NSDecimal)
= operate_function(operator);
NSDecimal (*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++) {
rounding = math_flags.round;
data[i] = operate(data[i], caster(otherData, j));
j += increment;
}
return self;
}
- _maPerform: (double_func_t)mathFunction
{
unsigned long i, max_elements;
NSDecimal *data;
/* Can't do functions on NSDecimal numbers */
MA_RAISE_PERFORM
max_elements = array_num_elements(dimension, [size bytes]);
data = [arrayData mutableBytes];
for (i=0; i < max_elements; i++)
{
NSDecimalNumber* number;
number = [NSDecimalNumber decimalNumberWithDecimal: data[i]];
data[i] = NSSetDecimal(mathFunction([number doubleValue]));
}
/* Raise an error since we are not really operating on NSDecimal numbers */
ma_fpe_errno |= FLT_OTHER;
return self;
}
- _maPerformFunction:(perform_func_t)perform_func userInfo:(void *)info
{
unsigned long i, max_elements;
unsigned *index = [[[size mutableCopy] autorelease] mutableBytes];
NSDecimal *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++) {
NSDecimalNumber *number;
NSDecimal temp;
[[arrayData valueAtIndex:i] getValue: &temp];
number = [NSDecimalNumber decimalNumberWithDecimal: temp];
// FIXME: Catch the exceptions here
number = perform_func(number, index, info);
data[i] = [number decimalValue];
increment_index_in_range(dimension, range, index, 1);
}
return self;
}
- maInvert
{
[self maNotImplemented:_cmd];
return self;
}
// FIXME: use NSDecimal to calculate random numbers...
/* Generate a uniform random number n, such that min <= n < max */
- maRandomUniformMin: (double)min toMax: (double)max
{
unsigned long i, max_elements;
NSDecimal *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] = NSSetDecimal(GEN_RANDOM * (max - min) + min);
}
return self;
}
/* Generate a normal 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;
NSDecimal *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] = NSSetDecimal(r * cos(theta) + mean);
data[i*2+1] = NSSetDecimal(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] = NSSetDecimal(r * cos(theta) + mean);
}
return self;
}
- (id <NSNumber,ComplexNumber>)_loopArrayWith:
(NSDecimal (*)(NSDecimal, NSDecimal))loopFunc
{
unsigned long i, max_elements;
const NSDecimal *data;
NSDecimal 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 [NSDecimalNumber decimalNumberWithDecimal: number];
}
- (id <NSNumber,ComplexNumber>)maMinimumValue
{
return [self _loopArrayWith:find_min];
}
- (id <NSNumber,ComplexNumber>)maMaximumValue
{
return [self _loopArrayWith:find_max];
}
- (id <NSNumber,ComplexNumber>)maTotal
{
rounding = math_flags.round;
return [self _loopArrayWith:find_sum];
}
- maWhere
{
ordered_index_t i, max_elements;
unsigned j;
NSDecimal *data;
NSDecimal zero;
j = 0;
max_elements = array_num_elements(dimension, [size bytes]);
data = [arrayData mutableBytes];
zero = NSSetDecimal(0);
for (i=0; i < max_elements; i++)
if (NSDecimalCompare(&data[i], &zero) != NSOrderedSame)
data[j++] = NSSetDecimal(i);
[arrayData setCount:j];
[size release];
dimension = 1;
size = [[MAValueData dataWithValues:&j
count:dimension
objCType:@encode(unsigned)] retain];
return self;
}
@end
inline static
NSDecimal find_min(NSDecimal a, NSDecimal b)
{
return (NSDecimalCompare(&a, &b) == NSOrderedAscending) ? a : b;
}
inline static
NSDecimal find_max(NSDecimal a, NSDecimal b)
{
return (NSDecimalCompare(&a, &b) == NSOrderedAscending) ? b : a;
}
inline static
NSDecimal find_sum(NSDecimal sum, NSDecimal num)
{
NSDecimal result;
NSCalculationError calc_error;
calc_error = NSDecimalAdd(&result, &sum, &num, rounding);
ma_fpe_errno |= trap_error_from_decimal(calc_error);
return result;
}
/* All the functions which know how to cast from any known type to our type */
static inline NSDecimal
cast_from_uchar(void *data, unsigned long loc)
{return NSSetDecimal(*((unsigned char*)data + loc));}
static inline NSDecimal
cast_from_char(void *data, unsigned long loc)
{return NSSetDecimal(*((char*)data + loc));}
static inline NSDecimal
cast_from_ushort(void *data, unsigned long loc)
{return NSSetDecimal(*((unsigned short*)data + loc));}
static inline NSDecimal
cast_from_short(void *data, unsigned long loc)
{return NSSetDecimal(*((short*)data + loc));}
static inline NSDecimal
cast_from_uint(void *data, unsigned long loc)
{return NSSetDecimal(*((unsigned int*)data + loc));}
static inline NSDecimal
cast_from_int(void *data, unsigned long loc)
{return NSSetDecimal(*((int*)data + loc));}
static inline NSDecimal
cast_from_ulong(void *data, unsigned long loc)
{return NSSetDecimal(*((unsigned long*)data + loc));}
static inline NSDecimal
cast_from_long(void *data, unsigned long loc)
{return NSSetDecimal(*((long*)data + loc));}
static inline NSDecimal
cast_from_ulonglong(void *data, unsigned long loc)
{return NSSetDecimal(*((unsigned long long*)data + loc));}
static inline NSDecimal
cast_from_longlong(void *data, unsigned long loc)
{return NSSetDecimal(*((long long*)data + loc));}
static inline NSDecimal
cast_from_float(void *data, unsigned long loc)
{return NSSetDecimal(*((float*)data + loc));}
static inline NSDecimal
cast_from_double(void *data, unsigned long loc)
{return NSSetDecimal(*((double*)data + loc));}
static inline NSDecimal
cast_from_complexfloat(void *data, unsigned long loc)
{return NSSetDecimal(*((float*)data + 2*loc));}
static inline NSDecimal
cast_from_complexdouble(void *data, unsigned long loc)
{return NSSetDecimal(*((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 */
/*
static inline NSDecimal
op_exponent(NSDecimal a, NSDecimal b)
{
NSDecimal result;
NSCalculationError calc_error;
calc_error = NSDecimalPower(&result, &a, &b, rounding);
ma_fpe_errno |= trap_error_from_decimal(calc_error);
return result;
}
*/
static inline NSDecimal
op_multiply(NSDecimal a, NSDecimal b)
{
NSDecimal result;
NSCalculationError calc_error;
calc_error = NSDecimalMultiply(&result, &a, &b, rounding);
ma_fpe_errno |= trap_error_from_decimal(calc_error);
return result;
}
static inline NSDecimal
op_divide(NSDecimal a, NSDecimal b)
{
NSDecimal result;
NSCalculationError calc_error;
calc_error = NSDecimalDivide(&result, &a, &b, rounding);
ma_fpe_errno |= trap_error_from_decimal(calc_error);
return result;
}
static inline NSDecimal
op_add(NSDecimal a, NSDecimal b)
{
NSDecimal result;
NSCalculationError calc_error;
calc_error = NSDecimalAdd(&result, &a, &b, rounding);
ma_fpe_errno |= trap_error_from_decimal(calc_error);
return result;
}
static inline NSDecimal
op_subtract(NSDecimal a, NSDecimal b)
{
NSDecimal result;
NSCalculationError calc_error;
calc_error = NSDecimalSubtract(&result, &a, &b, rounding);
ma_fpe_errno |= trap_error_from_decimal(calc_error);
return result;
}
/*
static inline NSDecimal
op_mod(NSDecimal a, NSDecimal b)
{
ma_fpe_errno = FLT_OTHER;
return NSSetDecimal(0);
}
*/
static inline NSDecimal
op_minimum(NSDecimal a, NSDecimal b)
{
return (NSDecimalCompare(&a, &b) == NSOrderedAscending) ? a : b;
}
static inline NSDecimal
op_maximum(NSDecimal a, NSDecimal b)
{
return (NSDecimalCompare(&a, &b) == NSOrderedAscending) ? b : a;
}
static inline NSDecimal
op_not(NSDecimal a, NSDecimal b)
{
NSDecimal zero = NSSetDecimal(0);
return (NSDecimalCompare(&a, &zero) == NSOrderedSame)
? NSSetDecimal(1) : zero;
}
static inline NSDecimal
op_equal(NSDecimal a, NSDecimal b)
{
return (NSDecimalCompare(&a, &b) == NSOrderedSame)
? NSSetDecimal(1) : NSSetDecimal(0);
}
static inline NSDecimal
op_not_equal(NSDecimal a, NSDecimal b)
{
return (NSDecimalCompare(&a, &b) != NSOrderedSame)
? NSSetDecimal(1) : NSSetDecimal(0);
}
static inline NSDecimal
op_less_or_equal(NSDecimal a, NSDecimal b)
{
return (NSDecimalCompare(&a, &b) <= NSOrderedSame)
? NSSetDecimal(1) : NSSetDecimal(0);
}
static inline NSDecimal
op_less(NSDecimal a, NSDecimal b)
{
return (NSDecimalCompare(&a, &b) == NSOrderedAscending)
? NSSetDecimal(1) : NSSetDecimal(0);
}
static inline NSDecimal
op_greater_or_equal(NSDecimal a, NSDecimal b)
{
return (NSDecimalCompare(&a, &b) >= NSOrderedSame)
? NSSetDecimal(1) : NSSetDecimal(0);
}
static inline NSDecimal
op_greater(NSDecimal a, NSDecimal b)
{
return (NSDecimalCompare(&a, &b) == NSOrderedDescending)
? NSSetDecimal(1) : NSSetDecimal(0);
}
// FIXME: How do you do logical ops on floats?
static inline NSDecimal
op_and(NSDecimal a, NSDecimal b)
{
NSDecimal zero = NSSetDecimal(0);
return (NSDecimalCompare(&a, &zero) != NSOrderedSame
&& NSDecimalCompare(&b, &zero) != NSOrderedSame) ? b : zero;
}
static inline NSDecimal
op_or(NSDecimal a, NSDecimal b)
{
NSDecimal zero = NSSetDecimal(0);
return (NSDecimalCompare(&a, &zero) != NSOrderedSame
|| NSDecimalCompare(&b, &zero) != NSOrderedSame) ? b : zero;
}
static inline NSDecimal
op_xor(NSDecimal a, NSDecimal b)
{
NSDecimal zero = NSSetDecimal(0);
return ((NSDecimalCompare(&a, &zero) != NSOrderedSame
|| NSDecimalCompare(&b, &zero) != NSOrderedSame)
&& !(NSDecimalCompare(&a, &zero) != NSOrderedSame
&& NSDecimalCompare(&b, &zero) != NSOrderedSame)) ? b : zero;
}
static operate_f_t
operate_function(ma_operator_t operator)
{
switch(operator) {
case ma_exponent:
MA_RAISE_PERFORM
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:
MA_RAISE_PERFORM
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.