Use libj2 for double-max precision integersHEADmaster

Signed-off-by: Mattias Andrée <m@maandree.se>

8 files changed, 74 insertions, 190 deletions

diff --git a/Makefile b/Makefile
index 26127d3..9bf77c2 100644
--- a/Makefile
+++ b/Makefile
@@ -16,7 +16,7 @@ LIB_VERSION = $(LIB_MAJOR).$(LIB_MINOR)
 LIB_NAME = quanta
 
 
-OBJ_PUBLIC =\
+OBJ =\
 	libquanta_palette_size.o\
 	libquanta_malloc_palette.o\
 	libquanta_vquantise.o\
@@ -24,13 +24,6 @@ OBJ_PUBLIC =\
 	libquanta_vquantise_wu.o\
 	libquanta_quantise_wu.o
 
-OBJ_PRIVATE =\
-	libquanta_bigint_divmod_small__.o
-
-OBJ =\
-	$(OBJ_PUBLIC)\
-	$(OBJ_PRIVATE)
-
 HDR =\
 	libquanta.h\
 	common.h
@@ -38,7 +31,7 @@ HDR =\
 LOBJ = $(OBJ:.o=.lo)
 
 MAN3 =\
-	$(OBJ_PUBLIC:.o=.3)\
+	$(OBJ:.o=.3)\
 	struct_libquanta_image.3 libquanta_image.3\
 	struct_libquanta_channel.3 libquanta_channel.3\
 	struct_libquanta_palette.3 libquanta_palette.3\
diff --git a/common.h b/common.h
index 3c4726e..3054683 100644
--- a/common.h
+++ b/common.h
@@ -11,82 +11,3 @@
 #define PALETTE_VALUE_TYPE      uint64_t
 #define PALETTE_VALUE_SIZE      sizeof(PALETTE_VALUE_TYPE)
 #define PALETTE_VALUE_MAX_BITS  64U
-
-
-struct bigint {
-	uintmax_t high, low;
-};
-
-
-#if defined(__GNUC__)
-__attribute__((__visibility__("hidden")))
-#endif
-uintmax_t libquanta_bigint_divmod_small__(struct bigint *big, uintmax_t small);
-#define bigint_divmod_small libquanta_bigint_divmod_small__
-
-
-static inline void
-bigint_add_small(struct bigint *big, uintmax_t small)
-{
-#if defined(__GNUC__)
-	if (__builtin_add_overflow(big->low, small, &big->low))
-		big->high += 1U;
-#else
-	if (big->low > UINTMAX_MAX - small)
-		big->high += 1U;
-	big->low += small;
-#endif
-}
-
-
-static inline void
-bigint_sub_small(struct bigint *big, uintmax_t small)
-{
-#if defined(__GNUC__)
-	if (__builtin_sub_overflow(big->low, small, &big->low))
-		big->high -= 1U;
-#else
-	if (big->low < small)
-		big->high -= 1U;
-	big->low -= small;
-#endif
-}
-
-
-static inline void
-bigint_add_big(struct bigint *res, const struct bigint *restrict other)
-{
-	bigint_add_small(res, other->low);
-	res->high += other->high;
-}
-
-
-static inline void
-bigint_sub_big(struct bigint *res, const struct bigint *restrict other)
-{
-	bigint_sub_small(res, other->low);
-	res->high -= other->high;
-}
-
-
-static inline void
-bigint_rsub_big(struct bigint *res, const struct bigint *restrict other)
-{
-	res->high = other->high - res->high;
-#if defined(__GNUC__)
-	if (__builtin_sub_overflow(other->low, res->low, &res->low))
-		res->high -= 1U;
-#else
-	if (res->low > other->low)
-		res->high -= 1U;
-	res->low = other->low - res->low;
-#endif
-}
-
-
-static inline uintmax_t
-bigint_div_small(const struct bigint *big, uintmax_t small)
-{
-	struct bigint big_copy = *big;
-	return bigint_divmod_small(&big_copy, small);
-}
diff --git a/config.mk b/config.mk
index 037fe35..71198ec 100644
--- a/config.mk
+++ b/config.mk
@@ -5,4 +5,4 @@ CC = c99
 
 CPPFLAGS = -D_DEFAULT_SOURCE -D_BSD_SOURCE -D_XOPEN_SOURCE=700 -D_GNU_SOURCE
 CFLAGS   =
-LDFLAGS  = -lm
+LDFLAGS  = -lm -lj2
diff --git a/libquanta.7 b/libquanta.7
index d68b73e..4999b8b 100644
--- a/libquanta.7
+++ b/libquanta.7
@@ -9,6 +9,7 @@ libquanta \- Colour quantisation library
 .PP
 Link with
 .I -lquanta
+.I -lj2
 .IR -lm .
 
 .SH DESCRIPTION
diff --git a/libquanta_bigint_divmod_small__.c b/libquanta_bigint_divmod_small__.c
deleted file mode 100644
index a2c345f..0000000
--- a/libquanta_bigint_divmod_small__.c
+++ /dev/null
@@ -1,34 +0,0 @@
-/* See LICENSE file for copyright and license details. */
-#include "common.h"
-
-
-uintmax_t
-libquanta_bigint_divmod_small__(struct bigint *big, uintmax_t small)
-{
-	uintmax_t q = 0, hi, lo;
-	int e = 8 * (int)sizeof(small);
-
-#if 0 /* this would overflow (undefined behaviour) and not fit in the result */
-	q = big->high / small;
-	q <<= 8 * (int)sizeof(small);
-#endif
-	big->high %= small;
-
-	while (big->high && --e) {
-		hi = small >> (8 * (int)sizeof(small) - e);
-		lo = small << e;
-
-		if (hi > big->high)
-			continue;
-		if (hi == big->high && lo > big->low)
-			continue;
-
-		q |= (uintmax_t)1 << e;
-		bigint_sub_small(big, lo);
-		big->high -= hi;
-	}
-
-	q += big->low / small;
-	big->low %= small;
-	return q;
-}
diff --git a/libquanta_quantise.3 b/libquanta_quantise.3
index 3a5db49..fbb6dcb 100644
--- a/libquanta_quantise.3
+++ b/libquanta_quantise.3
@@ -33,6 +33,7 @@ struct libquanta_image *\fBlibquanta_vquantise\fP(struct libquanta_palette *\fIp
 .PP
 Link with
 .I -lquanta
+.I -lj2
 .IR -lm .
 
 .SH DESCRIPTION
diff --git a/libquanta_quantise_wu.3 b/libquanta_quantise_wu.3
index d7673e8..264109e 100644
--- a/libquanta_quantise_wu.3
+++ b/libquanta_quantise_wu.3
@@ -12,6 +12,7 @@ struct libquanta_image *\fBlibquanta_vquantise_wu\fP(struct libquanta_palette *\
 .PP
 Link with
 .I -lquanta
+.I -lj2
 .IR -lm .
 
 .SH DESCRIPTION
diff --git a/libquanta_vquantise_wu.c b/libquanta_vquantise_wu.c
index d148c96..25a4ad2 100644
--- a/libquanta_vquantise_wu.c
+++ b/libquanta_vquantise_wu.c
@@ -1,5 +1,6 @@
 /* See LICENSE file for copyright and license details. */
 #include "common.h"
+#include <libj2.h>
 
 /* TODO limitation: 32 values per channel, limits palette size */
 #define WORKING_BITS 5U
@@ -10,10 +11,10 @@ struct channel_data {
 	size_t subindex_multiplier;
 	size_t colour_shift_down;
 	size_t colour_shift_up;
-	struct bigint *histogram;
-	struct bigint whole;
-	struct bigint half;
-	struct bigint base;
+	struct libj2_j2u *histogram;
+	struct libj2_j2u whole;
+	struct libj2_j2u half;
+	struct libj2_j2u base;
 	long double max;
 	size_t cut;
 	int has_cut;
@@ -68,7 +69,7 @@ compute_histogram(struct libquanta_image *image, struct histogram *histogram, st
 				v <<= channels[ch].colour_shift_up;
 				v += 1U;
 				index += v * channels[ch].subindex_multiplier;
-				bigint_add_small(&channels[ch].histogram[index], v);
+				libj2_j2u_add_ju(&channels[ch].histogram[index], v);
 			}
 			histogram[index].occurrences += 1U;
 			histogram[index].square_sum += square_subsum;
@@ -79,37 +80,37 @@ compute_histogram(struct libquanta_image *image, struct histogram *histogram, st
 
 
 static void
-moments_bigint_(size_t channel, const struct channel_data *channels, size_t nchannels,
-                struct bigint *out, const struct bigint *in, size_t index)
+moments_j2u_(size_t channel, const struct channel_data *channels, size_t nchannels,
+             struct libj2_j2u *out, const struct libj2_j2u *in, size_t index)
 {
 	size_t i, mul;
 
 	if (!nchannels) {
 		if (in != out)
 			memcpy(&out[index], &in[index], sizeof(*in));
-		bigint_add_big(&out[index], &out[index - channels[channel].subindex_multiplier]);
+		libj2_j2u_add_j2u(&out[index], &out[index - channels[channel].subindex_multiplier]);
 		return;
 	}
 
 	mul = channels[--nchannels].subindex_multiplier;
 	for (i = 0; i++ < (1U << WORKING_BITS);)
-		moments_bigint_(channel, channels, nchannels, out, in, index + i * mul);
+		moments_j2u_(channel, channels, nchannels, out, in, index + i * mul);
 }
 
 static void
-moments_bigint(size_t channel, struct channel_data *channels, size_t nchannels, struct bigint *buf, size_t bufsize)
+moments_j2u(size_t channel, struct channel_data *channels, size_t nchannels, struct libj2_j2u *buf, size_t bufsize)
 {
-	struct bigint *in = channels[channel].histogram;
+	struct libj2_j2u *in = channels[channel].histogram;
 	size_t i;
 
 	memset(buf, 0, bufsize);
 
 	for (i = nchannels; i-- > 1U;) {
-		moments_bigint_(i, channels, nchannels, buf, in, 0);
+		moments_j2u_(i, channels, nchannels, buf, in, 0);
 		in = buf;
 	}
 
-	moments_bigint_(0, channels, nchannels, channels[channel].histogram, in, 0);
+	moments_j2u_(0, channels, nchannels, channels[channel].histogram, in, 0);
 }
 
 static void
@@ -153,12 +154,12 @@ compute_cumulative_moments(struct histogram *histogram, struct channel_data *cha
 	void *buf;
 	size_t i, bufsize;
 
-	bufsize = bufelems * sizeof(struct bigint);
+	bufsize = bufelems * sizeof(struct libj2_j2u);
 	buf = malloc(bufsize);
 	if (!buf)
 		return -1;
 	for (i = 0; i < nchannels; i++)
-		moments_bigint(i, channels, nchannels, buf, bufsize);
+		moments_j2u(i, channels, nchannels, buf, bufsize);
 	free(buf);
 
 	bufsize = bufelems * sizeof(struct histogram);
@@ -244,38 +245,38 @@ volume_over_occurrences(const struct cube *cube, const struct histogram *data, c
 
 
 static void
-volume_over_bigints_(const struct cube *cube, const struct bigint *data, const struct channel_data *channels, size_t nchannels,
-                     size_t index, int use_addition, struct bigint *result)
+volume_over_j2us_(const struct cube *cube, const struct libj2_j2u *data, const struct channel_data *channels, size_t nchannels,
+                  size_t index, int use_addition, struct libj2_j2u *result)
 {
 	size_t ch;
 
 	if (!nchannels) {
 		if (use_addition)
-			bigint_add_big(result, &data[index]);
+			libj2_j2u_add_j2u(result, &data[index]);
 		else
-			bigint_sub_big(result, &data[index]);
+			libj2_j2u_sub_j2u(result, &data[index]);
 		return;
 	}
 
 	ch = nchannels - 1U;
 
-	volume_over_bigints_(cube, data, channels, ch,
-	                     index + cube->bounds[ch].max * channels[ch].subindex_multiplier,
-	                     use_addition ^ 1, result);
+	volume_over_j2us_(cube, data, channels, ch,
+	                  index + cube->bounds[ch].max * channels[ch].subindex_multiplier,
+	                  use_addition ^ 1, result);
 
-	volume_over_bigints_(cube, data, channels, ch,
-	                     index + cube->bounds[ch].min * channels[ch].subindex_multiplier,
-	                     use_addition, result);
+	volume_over_j2us_(cube, data, channels, ch,
+	                  index + cube->bounds[ch].min * channels[ch].subindex_multiplier,
+	                  use_addition, result);
 }
 
 static void
-volume_over_bigints(const struct cube *cube, const struct bigint *data, const struct channel_data *channels, size_t nchannels,
-                    struct bigint *result)
+volume_over_j2us(const struct cube *cube, const struct libj2_j2u *data, const struct channel_data *channels, size_t nchannels,
+                    struct libj2_j2u *result)
 {
 	int use_addition = (int)(~nchannels & 1U);
 	result->high = 0;
 	result->low = 0;
-	volume_over_bigints_(cube, data, channels, nchannels, 0, use_addition, result);
+	volume_over_j2us_(cube, data, channels, nchannels, 0, use_addition, result);
 }
 
 
@@ -320,41 +321,41 @@ bottom_over_occurrences(const struct cube *cube, size_t channel, const struct hi
 
 
 static void
-bottom_over_bigints_(const struct cube *cube, size_t channel, const struct bigint *data, const struct channel_data *channels,
-                     size_t nchannels, size_t index, int use_addition, struct bigint *result)
+bottom_over_j2us_(const struct cube *cube, size_t channel, const struct libj2_j2u *data, const struct channel_data *channels,
+                  size_t nchannels, size_t index, int use_addition, struct libj2_j2u *result)
 {
 	size_t ch;
 
 	if (!nchannels) {
 		if (use_addition)
-			bigint_add_big(result, &data[index]);
+			libj2_j2u_add_j2u(result, &data[index]);
 		else
-			bigint_sub_big(result, &data[index]);
+			libj2_j2u_sub_j2u(result, &data[index]);
 		return;
 	}
 
 	ch = nchannels - 1U;
 
-	bottom_over_bigints_(cube, channel, data, channels, ch,
-	                     index + cube->bounds[ch].min * channels[ch].subindex_multiplier,
-	                     use_addition, result);
+	bottom_over_j2us_(cube, channel, data, channels, ch,
+	                  index + cube->bounds[ch].min * channels[ch].subindex_multiplier,
+	                  use_addition, result);
 
 	if (ch == channel)
 		return;
 
-	bottom_over_bigints_(cube, channel, data, channels, ch,
-	                     index + cube->bounds[ch].max * channels[ch].subindex_multiplier,
-	                     use_addition ^ 1, result);
+	bottom_over_j2us_(cube, channel, data, channels, ch,
+	                  index + cube->bounds[ch].max * channels[ch].subindex_multiplier,
+	                  use_addition ^ 1, result);
 }
 
 static void
-bottom_over_bigints(const struct cube *cube, size_t channel, const struct bigint *data, const struct channel_data *channels,
-                    size_t nchannels, struct bigint *result)
+bottom_over_j2us(const struct cube *cube, size_t channel, const struct libj2_j2u *data, const struct channel_data *channels,
+                    size_t nchannels, struct libj2_j2u *result)
 {
 	int use_addition = (int)(~nchannels & 1U);
 	result->high = 0;
 	result->low = 0;
-	bottom_over_bigints_(cube, channel, data, channels, nchannels, 0, use_addition, result);
+	bottom_over_j2us_(cube, channel, data, channels, nchannels, 0, use_addition, result);
 }
 
 
@@ -403,45 +404,45 @@ top_over_occurrences(const struct cube *cube, size_t channel, size_t position, c
 
 
 static void
-top_over_bigints_(const struct cube *cube, size_t channel, size_t position, const struct bigint *data,
-                  const struct channel_data *channels, size_t nchannels, size_t index, int use_addition, struct bigint *result)
+top_over_j2us_(const struct cube *cube, size_t channel, size_t position, const struct libj2_j2u *data,
+               const struct channel_data *channels, size_t nchannels, size_t index, int use_addition, struct libj2_j2u *result)
 {
 	size_t ch;
 
 	if (!nchannels) {
 		if (use_addition)
-			bigint_add_big(result, &data[index]);
+			libj2_j2u_add_j2u(result, &data[index]);
 		else
-			bigint_sub_big(result, &data[index]);
+			libj2_j2u_sub_j2u(result, &data[index]);
 		return;
 	}
 
 	ch = nchannels - 1U;
 	
 	if (ch == channel) {
-		top_over_bigints_(cube, channel, position, data, channels, ch,
-		                  index + position * channels[ch].subindex_multiplier,
-		                  use_addition ^ 1, result);
+		top_over_j2us_(cube, channel, position, data, channels, ch,
+		               index + position * channels[ch].subindex_multiplier,
+		               use_addition ^ 1, result);
 
 	} else {
-		top_over_bigints_(cube, channel, position, data, channels, ch,
-		                  index + cube->bounds[ch].min * channels[ch].subindex_multiplier,
-		                  use_addition, result);
+		top_over_j2us_(cube, channel, position, data, channels, ch,
+		               index + cube->bounds[ch].min * channels[ch].subindex_multiplier,
+		               use_addition, result);
 
-		top_over_bigints_(cube, channel, position, data, channels, ch,
-		                  index + cube->bounds[ch].max * channels[ch].subindex_multiplier,
-		                  use_addition ^ 1, result);
+		top_over_j2us_(cube, channel, position, data, channels, ch,
+		               index + cube->bounds[ch].max * channels[ch].subindex_multiplier,
+		               use_addition ^ 1, result);
 	}
 }
 
 static void
-top_over_bigints(const struct cube *cube, size_t channel, size_t position, const struct bigint *data,
-                 const struct channel_data *channels, size_t nchannels, struct bigint *result)
+top_over_j2us(const struct cube *cube, size_t channel, size_t position, const struct libj2_j2u *data,
+                 const struct channel_data *channels, size_t nchannels, struct libj2_j2u *result)
 {
 	int use_addition = (int)(~nchannels & 1U);
 	result->high = 0;
 	result->low = 0;
-	top_over_bigints_(cube, channel, position, data, channels, nchannels, 0, use_addition, result);
+	top_over_j2us_(cube, channel, position, data, channels, nchannels, 0, use_addition, result);
 }
 
 
@@ -460,7 +461,7 @@ maximise(const struct cube *cube, size_t channel, const struct histogram *histog
 	channels[channel].has_cut = 0;
 
 	for (i = 0; i < nchannels; i++)
-		bottom_over_bigints(cube, channel, channels[i].histogram, channels, nchannels, &channels[i].base);
+		bottom_over_j2us(cube, channel, channels[i].histogram, channels, nchannels, &channels[i].base);
 	base_weight = bottom_over_occurrences(cube, channel, histogram, channels, nchannels);
 
 	first = cube->bounds[channel].min + 1U;
@@ -473,8 +474,8 @@ maximise(const struct cube *cube, size_t channel, const struct histogram *histog
 
 		temp = 0;
 		for (i = 0; i < nchannels; i++) {
-			top_over_bigints(cube, channel, pos, channels[i].histogram, channels, nchannels, &channels[i].half);
-			bigint_add_big(&channels[i].half, &channels[i].base);
+			top_over_j2us(cube, channel, pos, channels[i].histogram, channels, nchannels, &channels[i].half);
+			libj2_j2u_add_j2u(&channels[i].half, &channels[i].base);
 
 			v = (long double)channels[i].half.high;
 			v = ldexpl(v, 8 * sizeof(channels[i].half.low));
@@ -489,7 +490,7 @@ maximise(const struct cube *cube, size_t channel, const struct histogram *histog
 
 		temp2 = 0;
 		for (i = 0; i < nchannels; i++) {
-			bigint_rsub_big(&channels[i].half, &channels[i].whole); /* half = whole - half */
+			libj2_j2u_rsub_j2u(&channels[i].half, &channels[i].whole); /* half = whole - half */
 
 			v = (long double)channels[i].half.high;
 			v = ldexpl(v, 8 * sizeof(channels[i].half.low));
@@ -515,7 +516,7 @@ cut(struct cube *cube1, struct cube *cube2, const struct histogram *histogram, s
 	long double max_max;
 
 	for (i = 0; i < nchannels; i++)
-		volume_over_bigints(cube1, channels[i].histogram, channels, nchannels, &channels[i].whole);
+		volume_over_j2us(cube1, channels[i].histogram, channels, nchannels, &channels[i].whole);
 	whole_weight = volume_over_occurrences(cube1, histogram, channels, nchannels);
 
 	for (i = 0; i < nchannels; i++)
@@ -556,13 +557,13 @@ variance(const struct cube *cube, const struct histogram *histogram, const struc
 {
 	uintmax_t weight;
 	long double w, x, y = 0;
-	struct bigint v;
+	struct libj2_j2u v;
 	size_t i;
 
 	x = volume_over_square_sums(cube, histogram, channels, nchannels);
 
 	for (i = 0; i < nchannels; i++) {
-		volume_over_bigints(cube, channels[i].histogram, channels, nchannels, &v);
+		volume_over_j2us(cube, channels[i].histogram, channels, nchannels, &v);
 		w = (long double)v.high;
 		w = ldexpl(w, 8 * sizeof(v.low));
 		w += (long double)v.low;
@@ -658,7 +659,7 @@ libquanta_vquantise_wu(struct libquanta_palette *palette, va_list args)
 	struct libquanta_image *ret = NULL;
 	size_t size, image_size, *tag = NULL;
 	uintmax_t weight, mean;
-	struct bigint sum;
+	struct libj2_j2u sum;
 	struct cube **cubes = NULL;
 	size_t cube_size;
 	size_t ncolours_used;
@@ -767,8 +768,8 @@ libquanta_vquantise_wu(struct libquanta_palette *palette, va_list args)
 		if (!weight)
 			continue;
 		for (j = 0; j < nchannels; j++) {
-			volume_over_bigints(cubes[i], channels[j].histogram, channels, nchannels, &sum);
-			mean = bigint_div_small(&sum, weight);
+			volume_over_j2us(cubes[i], channels[j].histogram, channels, nchannels, &sum);
+			mean = libj2_j2u_div_ju_return(&sum, weight);
 			if (channels[j].ch->bits_out < channels[j].ch->bits_in)
 				mean >>= channels[j].ch->bits_in - channels[j].ch->bits_out;
 			palette->palette[n * nchannels + j] = (uint64_t)mean;