state from 2014-09-07 13:30 MESZ

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