diff options
Diffstat (limited to 'numlib')
-rw-r--r-- | numlib/Makefile.am | 15 | ||||
-rw-r--r-- | numlib/numsup.c | 333 | ||||
-rw-r--r-- | numlib/numsup.h | 70 |
3 files changed, 400 insertions, 18 deletions
diff --git a/numlib/Makefile.am b/numlib/Makefile.am deleted file mode 100644 index da95596..0000000 --- a/numlib/Makefile.am +++ /dev/null @@ -1,15 +0,0 @@ -include $(top_srcdir)/Makefile.shared - -privatelib_LTLIBRARIES = libargyllnum.la -privatelibdir = $(pkglibdir) - -libargyllnum_la_SOURCES = numlib.h numsup.c numsup.h dnsq.c dnsq.h \ - powell.c powell.h dhsx.c dhsx.h ludecomp.c ludecomp.h svd.c \ - svd.h zbrent.c zbrent.h rand.c rand.h sobol.c sobol.h aatree.c - -LDADD = ./libargyllnum.la - -check_PROGRAMS = dnsqtest tpowell tdhsx LUtest svdtest zbrenttest \ - soboltest - -EXTRA_DIST = License.txt Readme.txt diff --git a/numlib/numsup.c b/numlib/numsup.c index 7648a62..a28dcd2 100644 --- a/numlib/numsup.c +++ b/numlib/numsup.c @@ -293,6 +293,7 @@ a1log *new_a1log( a1loge(g_log, 1, "new_a1log: malloc of a1log failed, calling exit(1)\n"); exit(1); } + log->refc = 1; log->verb = verb; log->debug = debug; @@ -327,7 +328,7 @@ a1log *new_a1log_d(a1log *log) { /* Returns NULL */ a1log *del_a1log(a1log *log) { if (log != NULL) { - if (--log->refc <= 0) { + if (--log->refc == 0) { #ifdef NT DeleteCriticalSection(&log->lock); #endif @@ -518,6 +519,55 @@ error(char *fmt, ...) { /******************************************************************/ +/* Suplimental allcation functions */ +/******************************************************************/ + +#ifndef SIZE_MAX +# define SIZE_MAX ((size_t)(-1)) +#endif + +/* a * b */ +static size_t ssat_mul(size_t a, size_t b) { + size_t c; + + if (a == 0 || b == 0) + return 0; + + if (a > (SIZE_MAX/b)) + return SIZE_MAX; + else + return a * b; +} + +/* reallocate and clear new allocation */ +void *recalloc( /* Return new address */ +void *ptr, /* Current address */ +size_t cnum, /* Current number and unit size */ +size_t csize, +size_t nnum, /* New number and unit size */ +size_t nsize +) { + int ind = 0; + size_t ctot, ntot; + + if (ptr == NULL) + return calloc(nnum, nsize); + + if ((ntot = ssat_mul(nnum, nsize)) == SIZE_MAX) + return NULL; /* Overflow */ + + if ((ctot = ssat_mul(cnum, csize)) == SIZE_MAX) + return NULL; /* Overflow */ + + ptr = realloc(ptr, ntot); + + if (ptr != NULL && ntot > ctot) + memset((char *)ptr + ctot, 0, ntot - ctot); /* Clear the new region */ + + return ptr; +} + +/******************************************************************/ /* Numerical Recipes Vector/Matrix Support functions */ /******************************************************************/ /* Note the z suffix versions return zero'd vectors/matricies */ @@ -1537,3 +1587,284 @@ char *ctime_64(const INR64 *timer) { return rv; } + +/*******************************************/ +/* Native to/from byte buffer functions */ +/*******************************************/ + +/* No overflow detection is done - */ +/* numbers are clipped or truncated. */ + +/* be = Big Endian */ +/* le = Little Endian */ + +/* - - - - - - - - */ +/* Unsigned 8 bit */ + +unsigned int read_ORD8(ORD8 *p) { + unsigned int rv; + rv = ((unsigned int)p[0]); + return rv; +} + +void write_ORD8(unsigned int d, ORD8 *p) { + if (d > 0xff) + d = 0xff; + p[0] = (ORD8)(d); +} + +/* - - - - - - - - */ +/* Signed 8 bit */ + +int read_INR8(ORD8 *p) { + int rv; + rv = (int)(INR8)p[0]; + return rv; +} + +void write_INR8(int d, ORD8 *p) { + if (d > 0x7f) + d = 0x7f; + else if (d < -0x80) + d = -0x80; + p[0] = (ORD8)(d); +} + +/* - - - - - - - - */ +/* Unsigned 16 bit */ + +unsigned int read_ORD16_be(ORD8 *p) { + unsigned int rv; + rv = (((unsigned int)p[0]) << 8) + + (((unsigned int)p[1])); + return rv; +} + +unsigned int read_ORD16_le(ORD8 *p) { + unsigned int rv; + rv = (((unsigned int)p[0])) + + (((unsigned int)p[1]) << 8); + return rv; +} + +void write_ORD16_be(unsigned int d, ORD8 *p) { + if (d > 0xffff) + d = 0xffff; + p[0] = (ORD8)(d >> 8); + p[1] = (ORD8)(d); +} + +void write_ORD16_le(unsigned int d, ORD8 *p) { + if (d > 0xffff) + d = 0xffff; + p[0] = (ORD8)(d); + p[1] = (ORD8)(d >> 8); +} + +/* - - - - - - - - */ +/* Signed 16 bit */ + +int read_INR16_be(ORD8 *p) { + int rv; + rv = (((int)(INR8)p[0]) << 8) + + (((int)p[1])); + return rv; +} + +int read_INR16_le(ORD8 *p) { + int rv; + rv = (((int)p[0])) + + (((int)(INR8)p[1]) << 8); + return rv; +} + +void write_INR16_be(int d, ORD8 *p) { + if (d > 0x7fff) + d = 0x7fff; + else if (d < -0x8000) + d = -0x8000; + p[0] = (ORD8)(d >> 8); + p[1] = (ORD8)(d); +} + +void write_INR16_le(int d, ORD8 *p) { + if (d > 0x7fff) + d = 0x7fff; + else if (d < -0x8000) + d = -0x8000; + p[0] = (ORD8)(d); + p[1] = (ORD8)(d >> 8); +} + +/* - - - - - - - - */ +/* Unsigned 32 bit */ + +unsigned int read_ORD32_be(ORD8 *p) { + unsigned int rv; + rv = (((unsigned int)p[0]) << 24) + + (((unsigned int)p[1]) << 16) + + (((unsigned int)p[2]) << 8) + + (((unsigned int)p[3])); + return rv; +} + +unsigned int read_ORD32_le(ORD8 *p) { + unsigned int rv; + rv = (((unsigned int)p[0])) + + (((unsigned int)p[1]) << 8) + + (((unsigned int)p[2]) << 16) + + (((unsigned int)p[3]) << 24); + return rv; +} + +void write_ORD32_be(unsigned int d, ORD8 *p) { + p[0] = (ORD8)(d >> 24); + p[1] = (ORD8)(d >> 16); + p[2] = (ORD8)(d >> 8); + p[3] = (ORD8)(d); +} + +void write_ORD32_le(unsigned int d, ORD8 *p) { + p[0] = (ORD8)(d); + p[1] = (ORD8)(d >> 8); + p[2] = (ORD8)(d >> 16); + p[3] = (ORD8)(d >> 24); +} + +/* - - - - - - - - */ +/* Signed 32 bit */ + +int read_INR32_be(ORD8 *p) { + int rv; + rv = (((int)(INR8)p[0]) << 24) + + (((int)p[1]) << 16) + + (((int)p[2]) << 8) + + (((int)p[3])); + return rv; +} + +int read_INR32_le(ORD8 *p) { + int rv; + rv = (((int)p[0])) + + (((int)p[1]) << 8) + + (((int)p[2]) << 16) + + (((int)(INR8)p[3]) << 24); + return rv; +} + +void write_INR32_be(int d, ORD8 *p) { + p[0] = (ORD8)(d >> 24); + p[1] = (ORD8)(d >> 16); + p[2] = (ORD8)(d >> 8); + p[3] = (ORD8)(d); +} + +void write_INR32_le(int d, ORD8 *p) { + p[0] = (ORD8)(d); + p[1] = (ORD8)(d >> 8); + p[2] = (ORD8)(d >> 16); + p[3] = (ORD8)(d >> 24); +} + +/* - - - - - - - - */ +/* Unsigned 64 bit */ + +ORD64 read_ORD64_be(ORD8 *p) { + ORD64 rv; + rv = (((ORD64)p[0]) << 56) + + (((ORD64)p[1]) << 48) + + (((ORD64)p[2]) << 40) + + (((ORD64)p[3]) << 32) + + (((ORD64)p[4]) << 24) + + (((ORD64)p[5]) << 16) + + (((ORD64)p[6]) << 8) + + (((ORD64)p[7])); + return rv; +} + +ORD64 read_ORD64_le(ORD8 *p) { + ORD64 rv; + rv = (((ORD64)p[0])) + + (((ORD64)p[1]) << 8) + + (((ORD64)p[2]) << 16) + + (((ORD64)p[3]) << 24) + + (((ORD64)p[4]) << 32) + + (((ORD64)p[5]) << 40) + + (((ORD64)p[6]) << 48) + + (((ORD64)p[7]) << 56); + return rv; +} + +void write_ORD64_be(ORD64 d, ORD8 *p) { + p[0] = (ORD8)(d >> 56); + p[1] = (ORD8)(d >> 48); + p[2] = (ORD8)(d >> 40); + p[3] = (ORD8)(d >> 32); + p[4] = (ORD8)(d >> 24); + p[5] = (ORD8)(d >> 16); + p[6] = (ORD8)(d >> 8); + p[7] = (ORD8)(d); +} + +void write_ORD64_le(ORD64 d, ORD8 *p) { + p[0] = (ORD8)(d); + p[1] = (ORD8)(d >> 8); + p[2] = (ORD8)(d >> 16); + p[3] = (ORD8)(d >> 24); + p[4] = (ORD8)(d >> 32); + p[5] = (ORD8)(d >> 40); + p[6] = (ORD8)(d >> 48); + p[7] = (ORD8)(d >> 56); +} + +/* - - - - - - - - */ +/* Signed 64 bit */ + +INR64 read_INR64_be(ORD8 *p) { + INR64 rv; + rv = (((INR64)(INR8)p[0]) << 56) + + (((INR64)p[1]) << 48) + + (((INR64)p[2]) << 40) + + (((INR64)p[3]) << 32) + + (((INR64)p[4]) << 24) + + (((INR64)p[5]) << 16) + + (((INR64)p[6]) << 8) + + (((INR64)p[7])); + return rv; +} + +INR64 read_INR64_le(ORD8 *p) { + INR64 rv; + rv = (((INR64)p[0])) + + (((INR64)p[1]) << 8) + + (((INR64)p[2]) << 16) + + (((INR64)p[3]) << 24) + + (((INR64)p[4]) << 32) + + (((INR64)p[5]) << 40) + + (((INR64)p[6]) << 48) + + (((INR64)(INR8)p[7]) << 56); + return rv; +} + +void write_INR64_be(INR64 d, ORD8 *p) { + p[0] = (ORD8)(d >> 56); + p[1] = (ORD8)(d >> 48); + p[2] = (ORD8)(d >> 40); + p[3] = (ORD8)(d >> 32); + p[4] = (ORD8)(d >> 24); + p[5] = (ORD8)(d >> 16); + p[6] = (ORD8)(d >> 8); + p[7] = (ORD8)(d); +} + +void write_INR64_le(INR64 d, ORD8 *p) { + p[0] = (ORD8)(d); + p[1] = (ORD8)(d >> 8); + p[2] = (ORD8)(d >> 16); + p[3] = (ORD8)(d >> 24); + p[4] = (ORD8)(d >> 32); + p[5] = (ORD8)(d >> 40); + p[6] = (ORD8)(d >> 48); + p[7] = (ORD8)(d >> 56); +} + diff --git a/numlib/numsup.h b/numlib/numsup.h index 7db7676..b949379 100644 --- a/numlib/numsup.h +++ b/numlib/numsup.h @@ -41,6 +41,8 @@ /* This really needs checking for each different platform. */ /* Using C99 and MSC covers a lot of cases, */ /* and the fallback default is pretty reliable with modern compilers and machines. */ +/* Note that MSWin is LLP64 == 32 bit long, while OS X/Linux is LP64 == 64 bit long. */ +/* so long shouldn't really be used in any code.... */ /* (duplicated in icc.h) */ #if (__STDC_VERSION__ >= 199901L) /* C99 */ @@ -97,7 +99,7 @@ #else /* !_MSC_VER */ /* The following works on a lot of modern systems, including */ -/* LP64 model 64 bit modes */ +/* LLP64 and LP64 models, but won't work with ILP64 which needs int32 */ #define INR8 signed char /* 8 bit signed */ #define INR16 signed short /* 16 bit signed */ @@ -185,7 +187,7 @@ struct _a1log { int refc; /* Reference count */ - char *tag; /* Optional tag name */ + char *tag; /* Optional tag name */ int verb; /* Current verbosity level (public) */ int debug; /* Current debug level (public) */ @@ -277,6 +279,17 @@ extern char cr_char; /* =========================================================== */ +/* reallocate and clear new allocation */ +void *recalloc( /* Return new address */ +void *ptr, /* Current address */ +size_t cnum, /* Current number and unit size */ +size_t csize, +size_t nnum, /* New number and unit size */ +size_t nsize +); + +/* =========================================================== */ + /* Numerical recipes vector/matrix support functions */ /* Note that the index arguments are the inclusive low and high values */ @@ -372,6 +385,59 @@ char *ctime_32(const INR32 *timer); /* A static buffer is used. There is no \n at the end */ char *ctime_64(const INR64 *timer); +/*******************************************/ +/* Native to/from byte buffer functions */ +/*******************************************/ + +/* No overflow detection is done - numbers are clipped */ + +/* be = Big Endian */ +/* le = Little Endian */ + +/* Unsigned 8 bit */ +unsigned int read_ORD8(ORD8 *p); +void write_ORD8(unsigned int d, ORD8 *p); + +/* Signed 8 bit */ +int read_INR8(ORD8 *p); +void write_INR8(int d, ORD8 *p); + +/* Unsigned 16 bit */ +unsigned int read_ORD16_be(ORD8 *p); +unsigned int read_ORD16_le(ORD8 *p); +void write_ORD16_be(unsigned int d, ORD8 *p); +void write_ORD16_le(unsigned int d, ORD8 *p); + +/* Signed 16 bit */ +int read_INR16_be(ORD8 *p); +int read_INR16_le(ORD8 *p); +void write_INR16_be(int d, ORD8 *p); +void write_INR16_le(int d, ORD8 *p); + +/* Unsigned 32 bit */ +unsigned int read_ORD32_be(ORD8 *p); +unsigned int read_ORD32_le(ORD8 *p); +void write_ORD32_be(unsigned int d, ORD8 *p); +void write_ORD32_le(unsigned int d, ORD8 *p); + +/* Signed 32 bit */ +int read_INR32_be(ORD8 *p); +int read_INR32_le(ORD8 *p); +void write_INR32_be(int d, ORD8 *p); +void write_INR32_le(int d, ORD8 *p); + +/* Unsigned 64 bit */ +ORD64 read_ORD64_be(ORD8 *p); +ORD64 read_ORD64_le(ORD8 *p); +void write_ORD64_be(ORD64 d, ORD8 *p); +void write_ORD64_le(ORD64 d, ORD8 *p); + +/* Signed 64 bit */ +INR64 read_INR64_be(ORD8 *p); +INR64 read_INR64_le(ORD8 *p); +void write_INR64_be(INR64 d, ORD8 *p); +void write_INR64_le(INR64 d, ORD8 *p); + #ifdef __cplusplus } #endif |