summaryrefslogtreecommitdiff
path: root/numlib
diff options
context:
space:
mode:
Diffstat (limited to 'numlib')
-rw-r--r--numlib/Makefile.am15
-rw-r--r--numlib/numsup.c333
-rw-r--r--numlib/numsup.h70
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