From 655d34cb6c42204016facd6b643225e6edf7fd92 Mon Sep 17 00:00:00 2001
From: Florian Prill <florian.prill@dwd.de>
Date: Thu, 25 Apr 2019 11:06:02 +0200
Subject: [PATCH] [divide_timespans] implemented large integer division
algorithm for (days + msecs). - not yet enabled in mtime library.
---
examples/Makefile.in | 22 +++-
examples/int_div_example.c | 99 +++++++++++++++
include/int_div.h | 46 +++++++
src/Makefile.in | 11 +-
src/int_div.c | 253 +++++++++++++++++++++++++++++++++++++
5 files changed, 425 insertions(+), 6 deletions(-)
create mode 100644 examples/int_div_example.c
create mode 100644 include/int_div.h
create mode 100644 src/int_div.c
diff --git a/examples/Makefile.in b/examples/Makefile.in
index 337b5814..8938a376 100644
--- a/examples/Makefile.in
+++ b/examples/Makefile.in
@@ -93,7 +93,8 @@ noinst_PROGRAMS = time_calculus$(EXEEXT) model_integration$(EXEEXT) \
output_control$(EXEEXT) recurrence$(EXEEXT) repetitor$(EXEEXT) \
tas$(EXEEXT) uniq$(EXEEXT) modulo$(EXEEXT) \
comp_weights$(EXEEXT) callback_test$(EXEEXT) iconatm$(EXEEXT) \
- iconoce$(EXEEXT) iconoce_hl$(EXEEXT) simulate_iau$(EXEEXT)
+ iconoce$(EXEEXT) iconoce_hl$(EXEEXT) simulate_iau$(EXEEXT) \
+ int_div_example$(EXEEXT)
subdir = examples
ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
am__aclocal_m4_deps = $(top_srcdir)/m4/ax_prog_doxygen.m4 \
@@ -219,6 +220,12 @@ uniq_DEPENDENCIES = ../src/libmtime.la
uniq_LINK = $(LIBTOOL) $(AM_V_lt) --tag=FC $(AM_LIBTOOLFLAGS) \
$(LIBTOOLFLAGS) --mode=link $(FCLD) $(AM_FCFLAGS) $(FCFLAGS) \
$(uniq_LDFLAGS) $(LDFLAGS) -o $@
+am_int_div_example_OBJECTS = int_div_example.$(OBJEXT)
+int_div_example_OBJECTS = $(am_int_div_example_OBJECTS)
+int_div_example_DEPENDENCIES = ../src/libmtime.la
+int_div_example_LINK = $(LIBTOOL) $(AM_V_lt) --tag=FC $(AM_LIBTOOLFLAGS) \
+ $(LIBTOOLFLAGS) --mode=link $(FCLD) $(AM_FCFLAGS) $(FCFLAGS) \
+ $(int_div_example_LDFLAGS) $(LDFLAGS) -o $@
AM_V_P = $(am__v_P_@AM_V@)
am__v_P_ = $(am__v_P_@AM_DEFAULT_V@)
am__v_P_0 = false
@@ -277,14 +284,16 @@ SOURCES = $(callback_test_SOURCES) $(comp_weights_SOURCES) \
$(model_integration_SOURCES) $(modulo_SOURCES) \
$(output_control_SOURCES) $(recurrence_SOURCES) \
$(repetitor_SOURCES) $(simulate_iau_SOURCES) $(tas_SOURCES) \
- $(time_calculus_SOURCES) $(uniq_SOURCES)
+ $(time_calculus_SOURCES) $(uniq_SOURCES) \
+ $(int_div_example_SOURCES)
DIST_SOURCES = $(callback_test_SOURCES) $(comp_weights_SOURCES) \
$(duration_SOURCES) $(example_SOURCES) $(example_hl_SOURCES) \
$(iconatm_SOURCES) $(iconoce_SOURCES) $(iconoce_hl_SOURCES) \
$(model_integration_SOURCES) $(modulo_SOURCES) \
$(output_control_SOURCES) $(recurrence_SOURCES) \
$(repetitor_SOURCES) $(simulate_iau_SOURCES) $(tas_SOURCES) \
- $(time_calculus_SOURCES) $(uniq_SOURCES)
+ $(time_calculus_SOURCES) $(uniq_SOURCES) \
+ $(int_div_example_SOURCES)
am__can_run_installinfo = \
case $$AM_UPDATE_INFO_DIR in \
n|no|NO) false;; \
@@ -515,6 +524,9 @@ iconoce_hl_LDFLAGS =
simulate_iau_SOURCES = simulate_iau.f90
simulate_iau_LDADD = ../src/libmtime.la
simulate_iau_LDFLAGS =
+int_div_example_SOURCES = int_div_example.f90
+int_div_example_LDADD = ../src/libmtime.la
+int_div_example_LDFLAGS =
all: all-am
.SUFFIXES:
@@ -626,6 +638,10 @@ uniq$(EXEEXT): $(uniq_OBJECTS) $(uniq_DEPENDENCIES) $(EXTRA_uniq_DEPENDENCIES)
@rm -f uniq$(EXEEXT)
$(AM_V_FCLD)$(uniq_LINK) $(uniq_OBJECTS) $(uniq_LDADD) $(LIBS)
+int_div_example$(EXEEXT): $(int_div_example_OBJECTS) $(int_div_example_DEPENDENCIES) $(EXTRA_int_div_example_DEPENDENCIES)
+ @rm -f int_div_example$(EXEEXT)
+ $(AM_V_FCLD)$(int_div_example_LINK) $(int_div_example_OBJECTS) $(int_div_example_LDADD) $(LIBS) -lgmp
+
mostlyclean-compile:
-rm -f *.$(OBJEXT)
diff --git a/examples/int_div_example.c b/examples/int_div_example.c
new file mode 100644
index 00000000..9bbc075d
--- /dev/null
+++ b/examples/int_div_example.c
@@ -0,0 +1,99 @@
+/* --------------------------------------------------------------------------------
+ * EXAMPLE.C
+ *
+ * 04/2019 : F. Prill, DWD
+ *
+ * limitations:
+ *
+ * - the remainder in the mtime struct "_divisionquotienttimespan" is
+ * limited to int64_t msecs.
+ *
+ * -------------------------------------------------------------------------------- */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <inttypes.h>
+#include "int_div.h"
+#include <gmp.h>
+
+int divide_timespan_mpz(const t_int days1, const t_int day_fraction1,
+ const t_int days2, const t_int day_fraction2,
+ t_int *q_decimal, t_int *remainder_days, t_int *remainder_ms);
+
+int main(void) {
+ const int NREPEAT = 3000000;
+
+ int ret, i, cnt = 0;
+ t_int days1, day_fraction1, days2, day_fraction2;
+ t_int q_decimal0 = 0, remainder_days0 = 0, remainder_ms0 = 0;
+ t_int q_decimal = 0, remainder_days = 0, remainder_ms = 0;
+
+ double f1 = 1728000000.0 / (double) RAND_MAX;
+ double f2 = 86400000.0 / (double) RAND_MAX;
+
+ printf("test program\n");
+ printf(" checking %d random divisions against GNU Multiple Precision Arithmetic Library\n",NREPEAT);
+ srand ( 123 );
+
+ for (i=0; i<NREPEAT; i++) {
+ /* test randomly generated input against GNU Multiple Precision
+ Arithmetic Library: */
+ days1 = (t_int) ( f1 * (double) rand() );
+ day_fraction1 = (t_int) ( f2 * (double) rand() );
+ days2 = (t_int) ( f1 * (double) rand() );
+ day_fraction2 = (t_int) ( f2 * (double) rand() );
+
+ ret = divide_timespan_mpz(days1, day_fraction1, days2, day_fraction2,
+ &q_decimal0, &remainder_days0, &remainder_ms0);
+ ret = divide_timespan(days1, day_fraction1, days2, day_fraction2,
+ &q_decimal, &remainder_days, &remainder_ms);
+
+ if ((ret != 0) ||
+ (q_decimal != q_decimal0) ||
+ (remainder_days != remainder_days0) ||
+ (remainder_ms != remainder_ms0))
+ {
+ printf("dividend: %" PRId64 " : %" PRId64 " \n", days1, day_fraction1);
+ printf("divisor : %" PRId64 " : %" PRId64 " \n", days2, day_fraction2);
+ printf("comparison: quotient_out = %" PRId64 ", remainder = %" PRId64 " : %" PRId64 "\n",
+ q_decimal0, remainder_days0, remainder_ms0);
+ printf("quotient_out = %" PRId64 ", remainder = %" PRId64 " : %" PRId64 "\n",
+ q_decimal, remainder_days, remainder_ms);
+ printf("return code = %d\n",ret);
+ cnt++;
+ }
+ }
+
+ printf(" %d discrepancies reported.\n", cnt);
+ printf("test program done after %d tests.\n", NREPEAT);
+ exit(0);
+}
+
+
+int divide_timespan_mpz(const t_int days1, const t_int day_fraction1,
+ const t_int days2, const t_int day_fraction2,
+ t_int *q_decimal, t_int *remainder_days, t_int *remainder_ms)
+{
+ mpz_t a0, a1, a2, b1, b2, q, r;
+
+ mpz_inits(a0,a1,a2,b1,b2,q,r,NULL);
+
+ mpz_set_ui(a0, 86400000);
+ mpz_set_ui(a1,days1);
+ mpz_set_ui(a2,day_fraction1);
+ mpz_set_ui(b1,days2);
+ mpz_set_ui(b2,day_fraction2);
+
+ mpz_addmul(a2,a1,a0);
+ mpz_addmul(b2,b1,a0);
+ mpz_fdiv_qr(q,r,a2,b2);
+
+ *q_decimal = mpz_get_ui(q);
+
+ mpz_fdiv_qr(a2,b2,r,a0);
+ *remainder_days = mpz_get_ui(a2);
+ *remainder_ms = mpz_get_ui(b2);
+
+ mpz_clears(a0,a1,a2,b1,b2,q,r,NULL);
+ return 0;
+}
diff --git a/include/int_div.h b/include/int_div.h
new file mode 100644
index 00000000..480bc4f9
--- /dev/null
+++ b/include/int_div.h
@@ -0,0 +1,46 @@
+/*
+ Large integer division algorithm for (days + msecs).
+
+ 04/2019 : F. Prill, DWD.
+
+ The method uses a basis b=1200 representation for time spans and the
+ division algorithm explained in
+
+ P. B. Hansen: Multiple-Length Division Revisited: A Tour of the Minefield (192)
+ Electrical Engineering and Computer Science Technical Reports. 166.
+*/
+
+#ifndef INT_DIV_H_
+#define INT_DIV_H_
+
+#include <stdint.h>
+
+/* array representation has (w+1) digits: */
+#define W_DIGITS 6
+
+typedef int64_t t_int;
+
+/* data type for a natural number represented by an array of
+ (W_DIGITS+1) digits in radix b. */
+typedef t_int t_number[W_DIGITS + 1];
+
+
+void number_init(t_number *x, const t_int val, const unsigned int m);
+void print_number(const t_number x);
+unsigned int length(const t_number x);
+void product(t_number *x, const t_number y, const t_int k, const t_int b, int* ret);
+void quotient(t_number *x, const t_number y, t_int k, const t_int b);
+void iremainder(t_number *x, const t_number y, t_int k, const t_int b);
+t_int trial(const t_number r, const t_number d, const int k, const int m, const t_int b, int* ret);
+int smaller(const t_number r, const t_number dq, const int k, const int m);
+void difference(t_number *r, const t_number dq, const int k, const int m, const t_int b, int *ret);
+int division(const t_number x, const t_number y, t_number *q, t_number *r, const t_int b);
+
+void convert_basis(const t_int v_in, t_number *v_out, const int digits, const int istart);
+t_int convert_to_decimal(const t_number v_in, const int istart, const int iend);
+
+int divide_timespan(const t_int days1, const t_int day_fraction1,
+ const t_int days2, const t_int day_fraction2,
+ t_int *q_decimal, t_int *remainder_days, t_int *remainder_ms);
+
+#endif /* INT_DIV_H_ */
diff --git a/src/Makefile.in b/src/Makefile.in
index 2cd42561..a8bad531 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -139,7 +139,7 @@ am__objects_1 = mtime_calendar.lo mtime_calendar360day.lo \
mtime_iso8601.lo mtime_date.lo mtime_datetime.lo \
mtime_julianDay.lo mtime_time.lo mtime_timedelta.lo \
mtime_eventList.lo mtime_eventHandling.lo mtime_utilities.lo \
- vsop87.lo kepler.lo orbit.lo
+ vsop87.lo kepler.lo orbit.lo int_div.lo
am__objects_2 = error_handling.lo mtime_constants.lo \
mtime_error_handling.lo mtime_c_bindings.lo libmtime.lo \
libmtime_hl.lo
@@ -175,7 +175,8 @@ am__depfiles_remade = ./$(DEPDIR)/kepler.Plo \
./$(DEPDIR)/mtime_julianDay.Plo ./$(DEPDIR)/mtime_time.Plo \
./$(DEPDIR)/mtime_timedelta.Plo \
./$(DEPDIR)/mtime_utilities.Plo ./$(DEPDIR)/orbit.Plo \
- ./$(DEPDIR)/vsop87.Plo
+ ./$(DEPDIR)/vsop87.Plo \
+ ./$(DEPDIR)/int_div.Plo
am__mv = mv -f
COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \
$(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS)
@@ -408,7 +409,8 @@ c_SOURCES = mtime_calendar.c \
mtime_utilities.c \
vsop87.c \
kepler.c \
- orbit.c
+ orbit.c \
+ int_div.c
# Fortran sources: this list is ordered with respect to the module
@@ -533,6 +535,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/mtime_utilities.Plo@am__quote@ # am--include-marker
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/orbit.Plo@am__quote@ # am--include-marker
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/vsop87.Plo@am__quote@ # am--include-marker
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/int_div.Plo@am__quote@ # am--include-marker
$(am__depfiles_remade):
@$(MKDIR_P) $(@D)
@@ -742,6 +745,7 @@ distclean: distclean-am
-rm -f ./$(DEPDIR)/mtime_utilities.Plo
-rm -f ./$(DEPDIR)/orbit.Plo
-rm -f ./$(DEPDIR)/vsop87.Plo
+ -rm -f ./$(DEPDIR)/int_div.Plo
-rm -f Makefile
distclean-am: clean-am distclean-compile distclean-generic \
distclean-tags
@@ -803,6 +807,7 @@ maintainer-clean: maintainer-clean-am
-rm -f ./$(DEPDIR)/mtime_utilities.Plo
-rm -f ./$(DEPDIR)/orbit.Plo
-rm -f ./$(DEPDIR)/vsop87.Plo
+ -rm -f ./$(DEPDIR)/int_div.Plo
-rm -f Makefile
maintainer-clean-am: distclean-am maintainer-clean-generic
diff --git a/src/int_div.c b/src/int_div.c
new file mode 100644
index 00000000..b947fabc
--- /dev/null
+++ b/src/int_div.c
@@ -0,0 +1,253 @@
+/*
+ Large integer division algorithm for (days + msecs).
+
+ 04/2019 : F. Prill, DWD.
+
+ The method uses a basis b=1200 representation for time spans and the
+ division algorithm explained in
+
+ P. B. Hansen: Multiple-Length Division Revisited: A Tour of the Minefield (192).
+ Electrical Engineering and Computer Science Technical Reports. 166.
+*/
+#include "int_div.h"
+
+const int Fail = -1;
+const t_int iZero = (t_int) 0;
+const t_int msecs_day = 86400000;
+t_int const bb[4] = { 1, 1200, 1440000, 1728000000 }; /* powers of b = 1200 */
+
+
+/* initialize digits m,m+1,...,W_DIGITS by a constant val. */
+void number_init(t_number *x, const t_int val, const unsigned int m) {
+ int i;
+ for (i=m; i<=W_DIGITS; ++i) (*x)[i]=val;
+}
+
+
+/* linear search to determine the length of a natural number. */
+unsigned int length(const t_number x) {
+ int i = W_DIGITS;
+ while ((x[i] == iZero) && (i >= 0)) i--;
+ return (i+1);
+}
+
+
+/* partial arithmetic wrt. radix b: product x=y*k, return ret!=0 if overflow. */
+void product(t_number *x, const t_number y, const t_int k, const t_int b, int* ret) {
+ int i, m=length(y);
+ t_int temp, carry = 0;
+
+ if (*ret != 0) return;
+ for (i=0; i<m; ++i) {
+ temp = y[i]*k + carry;
+ (*x)[i] = temp % b;
+ carry = temp/b;
+ }
+ if (m<=W_DIGITS)
+ (*x)[m] = carry;
+ else
+ *ret = (carry != iZero);
+ for (i=m+1; i<=W_DIGITS; ++i) (*x)[i] = iZero;
+}
+
+
+/* partial arithmetic wrt. radix b: quotient x=y/k. */
+void quotient(t_number *x, const t_number y, t_int k, const t_int b) {
+ int i, m = length(y);
+ t_int temp, carry = iZero;
+ number_init(x, iZero, m);
+
+ for (i=m-1; i>=0; i--) {
+ temp = carry*b + y[i];
+ (*x)[i] = temp / k;
+ carry = temp % k;
+ }
+}
+
+
+/* partial arithmetic wrt. radix b: remainder x of y/k. */
+void iremainder(t_number *x, const t_number y, t_int k, const t_int b) {
+ int m = length(y), i;
+ t_int carry = iZero;
+
+ number_init(x, iZero, 1);
+ for (i=m-1; i>=0; i--)
+ carry = (carry*b + y[i]) % k;
+ (*x)[0] = carry;
+}
+
+
+/* prefix arithmetic wrt. radix b: estimate of a quotient digit, return <0 if error. */
+t_int trial(const t_number r, const t_number d,
+ const int k, const int m, const t_int b, int* ret) {
+ t_int d2,r3,t = iZero;
+ int km = k+m;
+#if DEBUG
+ if ((m < 2) || (km < m) || (km>W_DIGITS)) *ret = Fail;
+#endif
+ if (*ret != 0) return iZero;
+ r3 = (r[km]*b + r[km-1])*b + r[km-2];
+ d2 = d[m-1]*b + d[m-2];
+
+ t = r3/d2;
+ return (t > b-1) ? (b-1) : t;
+}
+
+
+/* prefix arithmetic: prefix comparison r{m+1} < dq, return <0 if error. */
+int smaller(const t_number r, const t_number dq, const int k, const int m) {
+ int i=m;
+
+#if DEBUG
+ if ((k < 0) || (k+m < k) || (k+m>W_DIGITS)) return Fail;
+#endif
+ while ((r[i+k] == dq[i]) && (i>=0)) i--;
+ return r[i+k] < dq[i];
+}
+
+
+/* prefix arithmetic wrt. radix b: prefix subtraction r{m+1} := r{m+1} - dq, return ret!=0 if error. */
+void difference(t_number *r, const t_number dq, const int k, const int m, const t_int b, int *ret) {
+ int i;
+ t_int borrow, diff;
+
+#ifdef DEBUG
+ if ((k < 0) || (k+m < k) || (k+m>W_DIGITS)) *ret = Fail;
+#endif
+ if (*ret != 0) return;
+ borrow = iZero;
+ for (i=0; i<=m; i++) {
+ diff = (*r)[i+k] - dq[i] - borrow + b;
+ (*r)[i+k] = diff % b;
+ borrow = (t_int) 1 - diff / b;
+ }
+ *ret = (borrow != iZero);
+}
+
+
+/* complete algorithm for multiple length division x/y wrt. radix b, return <0 if error. */
+int division(const t_number x, const t_number y,
+ t_number *q, t_number *r, const t_int b) {
+ int i, k, m, n, ret = 0;
+ t_int y1, f, qt;
+ t_number d, dq, rm;
+
+ m = length(y);
+ if (m == 1) {
+ y1 = y[m-1];
+ if (y1 <= 0) return Fail;
+ quotient(q,x,y1,b);
+ iremainder(r,x,y1,b);
+ }
+ else {
+ n = length(x);
+ number_init(q, iZero, 0);
+ number_init(r, iZero, 0);
+ if (m > n)
+ for (i=0; i<n; i++) (*r)[i] = x[i];
+ else {
+ /* long division algorithm */
+ if ((m < 2) || (n < m) || (n>W_DIGITS)) return Fail;
+ f = b / (y[m-1] + (t_int) 1);
+ product(&rm,x,f,b, &ret);
+ product(&d, y,f,b, &ret);
+ for (k=n-m; (k>=0) && (ret == 0); k--) {
+ qt = trial(rm, d, k, m, b, &ret);
+ product(&dq,d,qt,b, &ret);
+
+ if ((ret == 0) && (smaller(rm,dq,k,m)))
+ product(&dq,d,--qt,b, &ret);
+
+ (*q)[k] = qt;
+ difference(&rm,dq,k,m,b, &ret);
+ }
+ quotient(r,rm,f,b);
+ }
+ }
+ return ret;
+}
+
+
+/* convert to basis "b" by a succession of Euclidean divisions, stores to v_out[istart...]. */
+void convert_basis(const t_int v_in, t_number *v_out,
+ const int digits, const int istart) {
+ int i;
+ t_int v = v_in;
+ for (i=digits; i>=0; i--) {
+ (*v_out)[istart+i] = v / bb[i];
+ v %= bb[i];
+ }
+}
+
+
+/* convert to decimal basis. */
+t_int convert_to_decimal(const t_number v_in,
+ const int istart, const int iend) {
+ int i, l = length(v_in);
+ t_int v_out = iZero;
+
+ if (iend < l) l=iend;
+ if ((l-istart) > 3) return (t_int) Fail;
+ for (i=0; i<(l-istart+1); i++)
+ v_out += v_in[i+istart] * bb[i];
+ return v_out;
+}
+
+
+/*
+ convert days and day fraction to b=1200 basis by a succession of
+ Euclidean divisions.
+ By pre-multiplying with a factor of 20, the day fraction can be
+ expressed by three digits wrt. the basis b=1200.
+
+ restrictions:
+
+ - the no. of digits in the algorithm restricts numbers to 1200^3
+ days, ie. approx. 4.7 mio years.
+ - the quotient is limited to 1200^3 in order to fit into int64_t
+ after conversion to decimal base.
+ - all numbers and time spans must be positive.
+*/
+int divide_timespan(const t_int days1, const t_int day_fraction1,
+ const t_int days2, const t_int day_fraction2,
+ t_int *q_decimal, t_int *remainder_days, t_int *remainder_ms)
+{
+ int ret = 0;
+ t_number p1, p2, q, r; /* dividend, divisor, quotient, remainder */
+ t_int qd = 0, rd = 0;
+
+ if ((day_fraction1 > msecs_day) || (day_fraction2 > msecs_day)) return Fail;
+
+ convert_basis(day_fraction1 * (t_int) 20, &p1, 2, 0);
+ convert_basis(days1, &p1, W_DIGITS-3, 3);
+ convert_basis(day_fraction2 * (t_int) 20, &p2, 2, 0);
+ convert_basis(days2, &p2, W_DIGITS-3, 3);
+
+ ret = division(p1, p2, &q, &r, (t_int) 1200); /* perform division. */
+
+ if (ret == 0) {
+ /* convert result back to b=10 basis. */
+ *q_decimal = qd = convert_to_decimal(q,0,W_DIGITS);
+ *remainder_days = convert_to_decimal(r,3,W_DIGITS);
+ *remainder_ms = rd = convert_to_decimal(r,0,2);
+ ret = (qd < iZero) || (rd < iZero);
+
+ *remainder_ms /= ((t_int) 20);
+ }
+ return ret;
+}
+
+
+#ifdef DEBUG
+
+#include <stdio.h>
+#include <inttypes.h>
+
+/* print number to screen. */
+void print_number(const t_number x) {
+ int i, m = length(x) - 1;
+ for (i=m; i>=0; i--)
+ printf("%" PRId64 " ", x[i]);
+ printf("\n");
+}
+#endif
--
GitLab