diff --git a/examples/Makefile.in b/examples/Makefile.in
index 337b58149dfa15ddef081c1bee95ccfbf226fe79..8938a376f0d07a0b8bc60e5ee56849e66bcfceb9 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/duration.f90 b/examples/duration.f90
index 77dc3f03e5a6183f0442d891e9eb1670432bf8d2..94f4d3d11b5c120e814098b57be37415c26b0bbb 100644
--- a/examples/duration.f90
+++ b/examples/duration.f90
@@ -65,7 +65,6 @@ contains
     REAL, intent(in) :: expected
 
     integer  :: ierr, ierrs
-    integer(i8)  :: denominator = 0
 
     TYPE(divisionquotienttimespan) :: quotient
     TYPE(datetime),     POINTER   :: dt1_num => NULL()
@@ -73,7 +72,9 @@ contains
     TYPE(datetime),     POINTER   :: dt1_denom => NULL()
     TYPE(datetime),     POINTER   :: dt2_denom => NULL()
     TYPE(timedelta),    POINTER   :: td => NULL()
-
+    TYPE(timedelta),    POINTER   :: td_num   => NULL()
+    TYPE(timedelta),    POINTER   :: td_denom => NULL()
+    TYPE(juliandelta),  POINTER   :: jd_denom => NULL()
 
     ierrs = 0
     dt1_num => newDatetime(dt1_dividend, errno = ierr)
@@ -85,12 +86,29 @@ contains
     dt2_denom => newDatetime(dt2_divisor, errno = ierr)
     ierrs = ierrs + ierr
 
-    CALL divideTwoDatetimeDiffsInSeconds(dt1_num,dt2_num,dt1_denom,dt2_denom,denominator,quotient)
-    print *, '(', TRIM(dt1_dividend), ' - ', TRIM(dt2_dividend), ') / ',  &
-      &      '(', TRIM(dt1_divisor), ' - ', TRIM(dt2_divisor), '):  ', REAL(quotient%remainder_in_ms) / 1000. / REAL(denominator)
+    td_num   => newTimeDelta('PT0H', errno = ierr)
+    ierrs = ierrs + ierr
+    td_denom => newTimeDelta('PT0H', errno = ierr)
+    ierrs = ierrs + ierr
+    td_num   = dt1_num   - dt2_num
+    td_denom = dt1_denom - dt2_denom
+
+    jd_denom => newJuliandelta("+", 0_i8, 0_i8, ierr)
+    ierrs = ierrs + ierr
+    CALL timeDeltaToJulianDelta(td_denom, dt2_denom, jd_denom)
+
+    CALL divideTimeDelta(td_num, td_denom, dt2_num, quotient)
+    print *, '(', TRIM(dt1_dividend), ' - ', TRIM(dt2_dividend), ') mod ',  &
+      &      '(', TRIM(dt1_divisor), ' - ', TRIM(dt2_divisor), '):  ',    &
+      &      (86400.*REAL(quotient%remainder_days) + 0.001*REAL(quotient%remainder_in_ms)) / &
+      &      (86400.*REAL(jd_denom%day) + 0.001*REAL(jd_denom%ms)), &
+      &      " (rel.)"
     print *, 'expected: ', expected
 
     CALL deallocateTimeDelta(td)
+    CALL deallocateTimeDelta(td_num)
+    CALL deallocateTimeDelta(td_denom)
+    CALL deallocateJulianDelta(jd_denom)
 
   end subroutine testDatetimeDiffComp
   
diff --git a/examples/example_hl.f90 b/examples/example_hl.f90
index f1162a6f34fdff1f35634ee32936a0c6216a80a4..af3d1d3d4d693f14db85fccab1941799360c04e4 100644
--- a/examples/example_hl.f90
+++ b/examples/example_hl.f90
@@ -18,6 +18,7 @@ PROGRAM example
 
   WRITE (0,*) "example_hl : test example"
 
+  WRITE (0,*) " "
   WRITE (0,*) "  setting calendar."
   CALL setCalendar(PROLEPTIC_GREGORIAN)
 
@@ -27,6 +28,7 @@ PROGRAM example
   !              CALL deallocateDatetime(dt)
 
   ! TYPE(t_datetime)  :: dt
+  WRITE (0,*) " "
   WRITE (0,*) "  testing assignment of t_datetime."
   dt = t_datetime("1970-01-01T00:00:00")
 
@@ -36,6 +38,7 @@ PROGRAM example
   !              CALL deallocateTimedelta(td)
 
   ! TYPE(t_timedelta) :: td
+  WRITE (0,*) " "
   WRITE (0,*) "  testing assignment of t_timedelta."
   td = t_timedelta("PT1H1M1S")
 
@@ -43,71 +46,87 @@ PROGRAM example
   !              CALL timedeltaToString(td, dstring2)
   !              WRITE (0,*) "compute ", dstring1, " + ", dstring2
 
-  WRITE (0,*) "compute ", dt%toString(), " + ", td%toString() 
-
   ! Adding time intervals is unchanged:
-
+  WRITE (0,*) " "
+  WRITE (0,*) "  adding time intervals."
+  WRITE (0,*) "  compute ", dt%toString(), " + ", td%toString() 
   dt = dt + td
-  WRITE (0,*) "result: ", dt%toString()
+  WRITE (0,*) "    result: ", dt%toString()
 
   ! --- Further examples:
 
   ! subtraction of two dates
+  WRITE (0,*) " "
+  WRITE (0,*) "  subtraction of two dates."
   dt2 = t_datetime("1970-01-01T00:00:00")
   td = dt-dt2
-  WRITE (0,*) "subtraction of dates: time delta: ", td%toString()
+  WRITE (0,*) "  subtraction of dates: time delta: ", td%toString()
 
   ! comparison of dates
+  WRITE (0,*) " "
+  WRITE (0,*) "  comparison of dates."
   dt  = t_datetime("1970-01-01T00:00:00")
   dt2 = t_datetime("1970-01-01T00:00:00")
-  WRITE (0,*) dt%toString(), " == ", dt2%toString(), ": ", (dt == dt2)
+  WRITE (0,*) "  ", dt%toString(), " == ", dt2%toString(), ": ", (dt == dt2)
   dt3 = t_datetime("1970-01-01T00:00:01")
-  WRITE (0,*) dt%toString(), " == ", dt3%toString(), ": ", (dt == dt3)
+  WRITE (0,*) "  ", dt%toString(), " == ", dt3%toString(), ": ", (dt == dt3)
 
   ! min / max
+  WRITE (0,*) " "
+  WRITE (0,*) "  min / max."
   dt4 = min(dt2,dt3)
-  WRITE (0,*) "MIN(", dt2%toString(), ", ", dt3%toString(), "): ", dt4%toString()
+  WRITE (0,*) "  MIN(", dt2%toString(), ", ", dt3%toString(), "): ", dt4%toString()
   dt4 = max(dt2,dt3)
-  WRITE (0,*) "MAX(", dt2%toString(), ", ", dt3%toString(), "): ", dt4%toString()
+  WRITE (0,*) "  MAX(", dt2%toString(), ", ", dt3%toString(), "): ", dt4%toString()
 
   ! interval assignment with milliseconds
   td = t_timedeltaFromMilliseconds(360000)
-  WRITE (0,*) "interval assignment with milliseconds: ", td%toString()
+  WRITE (0,*) " "
+  WRITE (0,*) "  interval assignment with milliseconds"
+  WRITE (0,*) "  interval assignment with milliseconds: ", td%toString()
 
+  WRITE (0,*) " "
+  WRITE (0,*) "  factor multiplication."
   td = t_timedelta("PT1H")
   td = td*0.5_c_double
-  WRITE (0,*) "PT1H * 0.5 = ", td%toString()
+  WRITE (0,*) "  PT1H * 0.5 = ", td%toString()
   td = t_timedelta("PT1H")
   td = 0.5_c_double*td
-  WRITE (0,*) "PT1H * 0.5 = ", td%toString()
+  WRITE (0,*) "  PT1H * 0.5 = ", td%toString()
 
   td = t_timedelta("PT1H")
   td = td*2
-  WRITE (0,*) "PT1H * 2 = ", td%toString()
+  WRITE (0,*) "  PT1H * 2 = ", td%toString()
   td = t_timedelta("PT1H")
   td = 2*td
-  WRITE (0,*) "PT1H * 2 = ", td%toString()
+  WRITE (0,*) "  PT1H * 2 = ", td%toString()
 
   ! division
+  WRITE (0,*) " "
+  WRITE (0,*) "  division."
   td1 = t_timedelta("PT23H42M")
-  dqts = td1%divideInSecondsBy(td)
-  WRITE (0,*) td1%toString(), " / ", td%toString(), " = ", dqts%quotient, " plus stuff"
+  dqts = td1%divideBy(td,dt)
+  WRITE (0,*) "  ", td1%toString(), " / ", td%toString(), " = ", dqts%quotient, " plus stuff"
 
   ! toSeconds, toMilliSeconds
-  WRITE (0,*) td%toString(), " is in seconds ", td%toSeconds(dt)
-  WRITE (0,*) td%toString(), " is in milliseconds ", td%toMilliSeconds(dt)
+  WRITE (0,*) " "
+  WRITE (0,*) "  toSeconds, toMilliSeconds."
+  WRITE (0,*) "  ", td%toString(), " is in seconds ", td%toSeconds(dt)
+  WRITE (0,*) "  ", td%toString(), " is in milliseconds ", td%toMilliSeconds(dt)
 
   dt   = t_datetime("1970-01-01T00:00:00")
   jday = dt%toJulianDay()
   dt2  = jday%toDateTime()
-  print *, dt%toString(), " = ", dt2%toString()
+  WRITE (0,*) "  ", dt%toString(), " = ", dt2%toString()
 
   ! register an error callback without stopping the application for
   ! our tests:
   CALL  register_finish_mtime_procedure(error_callback)
   
   ! produce errors
-  WRITE (0,*) 'error test: dt = t_datetime("1970--01-01T00:00:00")'
+  WRITE (0,*) " "
+  WRITE (0,*) "  produce errors"
+  WRITE (0,*) '  error test: dt = t_datetime("1970--01-01T00:00:00")'
   lerror = .FALSE.
   dt = t_datetime("1970--01-01T00:00:00")
   IF (.NOT. lerror)  WRITE(0,*) ERR_UNCAUGHT
@@ -120,58 +139,58 @@ PROGRAM example
   ! dt = dt5
   ! IF (.NOT. lerror)  WRITE(0,*) ERR_UNCAUGHT
 
-  WRITE (0,*) 'error test: dt5%toString()'
+  WRITE (0,*) '  error test: dt5%toString()'
   lerror = .FALSE.
   WRITE (0,*) dt5%toString()
   IF (.NOT. lerror)  WRITE(0,*) ERR_UNCAUGHT
 
-  WRITE (0,*) 'error test: dt5%to_posix_string()'
+  WRITE (0,*) '  error test: dt5%to_posix_string()'
   lerror = .FALSE.
   WRITE (0,*) dt5%toString("%s%d%LK")
   IF (.NOT. lerror)  WRITE(0,*) ERR_UNCAUGHT
 
-  WRITE (0,*) 'error test: dt = dt + td2'
+  WRITE (0,*) '  error test: dt = dt + td2'
   lerror = .FALSE.
   dt = dt + td2
   IF (.NOT. lerror)  WRITE(0,*) ERR_UNCAUGHT
 
-  WRITE (0,*) 'error test: dt = dt - td2'
+  WRITE (0,*) '  error test: dt = dt - td2'
   lerror = .FALSE.
   dt = dt - td2
   IF (.NOT. lerror)  WRITE(0,*) ERR_UNCAUGHT
 
-  WRITE (0,*) 'error test: dt = dt - td2'
+  WRITE (0,*) '  error test: dt = dt - td2'
   lerror = .FALSE.
   dt = dt - td2
   IF (.NOT. lerror)  WRITE(0,*) ERR_UNCAUGHT
 
-  WRITE (0,*) 'error test: td = t_timedelta(...)'
+  WRITE (0,*) '  error test: td = t_timedelta(...)'
   lerror = .FALSE.
   td = t_timedelta("P1lK")
   IF (.NOT. lerror)  WRITE(0,*) ERR_UNCAUGHT
 
   td = t_timedeltaFromMilliseconds(HUGE(INT(1)))
-  WRITE (0,*) "t_timedeltaFromMilliseconds(HUGE(INT(1))) : HUGE(INT(1)) = ", HUGE(INT(1))
-  WRITE (0,*) "td%toString() = ", td%toString()
+  WRITE (0,*) "  t_timedeltaFromMilliseconds(HUGE(INT(1))) : HUGE(INT(1)) = ", HUGE(INT(1))
+  WRITE (0,*) "  td%toString() = ", td%toString()
 
-  WRITE (0,*) 'error test: td = td * 0.000001D0'
-  WRITE (0,*) 'td%toString() = ', td%toString()
+  WRITE (0,*) '  error test: td = td * 0.000001D0'
+  WRITE (0,*) '  td%toString() = ', td%toString()
   lerror = .FALSE.
   td = td * 0.000001D0
   IF (.NOT. lerror)  WRITE(0,*) ERR_UNCAUGHT
 
   ! FIXME: This does no longer work, probably it initially worked for the wrong
   ! reasons?
-  WRITE (0,*) 'error test: td2 = td2 * 1'
+  WRITE (0,*) '  error test: td2 = td2 * 1'
   lerror = .FALSE.
   td2 = td2 * 1
   IF (.NOT. lerror)  WRITE(0,*) ERR_UNCAUGHT
 
   ! FIXME: This does no longer work, probably it initially worked for the wrong
   ! reasons?
-  WRITE (0,*) 'error test: td2%toString()'
+  WRITE (0,*) '  error test: td2%toString()'
   lerror = .FALSE.
-  WRITE (0,*) 'td2%toString() = ', td2%toString()
+  WRITE (0,*) '  td2%toString() = ', td2%toString()
   IF (.NOT. lerror)  WRITE(0,*) ERR_UNCAUGHT
 
   ! --------------------------------------------------------------------------------
@@ -179,7 +198,7 @@ PROGRAM example
 
   c_sign = '+'
   jd = t_juliandelta(sign=c_sign, day=99_c_int64_t, ms=123_c_int64_t)
-  WRITE (0,*) "error test: jd = t_juliandelta(sign='r', day='99', ms='123')"
+  WRITE (0,*) "  error test: jd = t_juliandelta(sign='r', day='99', ms='123')"
   lerror = .FALSE.
   c_sign = 'r'
   jd = t_juliandelta(sign=c_sign, day=99_c_int64_t, ms=123_c_int64_t)
diff --git a/examples/iconoce.f90 b/examples/iconoce.f90
index 067cf10f2a58583aa81b9b4080309fdffc53c95e..b711c11f0bbd3b6d310c920a4cfc7b90fecfb6c9 100644
--- a/examples/iconoce.f90
+++ b/examples/iconoce.f90
@@ -1,6 +1,7 @@
 module mo_event_manager
 
   use mtime
+  USE mtime_error_handling
 
   implicit none
 
diff --git a/examples/int_div_example.c b/examples/int_div_example.c
new file mode 100644
index 0000000000000000000000000000000000000000..4a9803e60ec14b3dddef5a4cabdeb6273a394f93
--- /dev/null
+++ b/examples/int_div_example.c
@@ -0,0 +1,111 @@
+/* --------------------------------------------------------------------------------
+ * 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_init(a0);
+  mpz_init(a1);
+  mpz_init(a2);
+  mpz_init(b1);
+  mpz_init(b2);
+  mpz_init(q);
+  mpz_init(r);
+
+  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_clear(a0);
+  mpz_clear(a1);
+  mpz_clear(a2);
+  mpz_clear(b1);
+  mpz_clear(b2);
+  mpz_clear(q);
+  mpz_clear(r);
+  return 0;
+}
diff --git a/examples/output_control.f90 b/examples/output_control.f90
index 3cd0349adc74d71b552745090716d033f1776b5a..ef58f9425ca75e6333bc8aad04ade74049efee70 100644
--- a/examples/output_control.f90
+++ b/examples/output_control.f90
@@ -1,6 +1,7 @@
 program output_control
 
   use mtime
+  USE mtime_error_handling
 
   implicit none
 
diff --git a/examples/tas.f90 b/examples/tas.f90
index fa0e1fe637ee9ac39a17c902e24e8ab7b48a14f4..0d42bb844fdb758507d328689d8ec3102998ae61 100644
--- a/examples/tas.f90
+++ b/examples/tas.f90
@@ -72,7 +72,7 @@ program tas
     call timedeltatostring(tdivisor, ctd)
     print *, ctd
 
-    call divideTimeDeltaInSeconds(tdividend, tdivisor, tq)
+    CALL divideTimeDelta(tdividend, tdivisor, dt1, tq)
 
     print *, tq
 
@@ -171,11 +171,11 @@ contains
     tddiff => newtimedelta('PT0S')
     tddiff = mtime_current - mtime_begin
     
-    call dividetimedeltainseconds(tddiff, vlsec, tq)
-    step = tq%quotient
-    
     mtime_step => newDatetime("0000-01-01T00:00:00.000"); 
-    
+
+    CALL divideTimeDelta(tddiff, vlsec, mtime_step, tq)
+    step = tq%quotient
+      
     if (step >= 0) then
     
       mtime_step = mtime_begin + step * vlsec
diff --git a/include/int_div.h b/include/int_div.h
new file mode 100644
index 0000000000000000000000000000000000000000..9092822eae55a20d85e0fa17b2ad77e23df29802
--- /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 (1992)
+                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/include/mtime_timedelta.h b/include/mtime_timedelta.h
index 72a9016582df62d1f88d4d0b0176c50f2938e97a..ec88e29940f4c8e6781ac1f9eb0f3c49d763ad8c 100644
--- a/include/mtime_timedelta.h
+++ b/include/mtime_timedelta.h
@@ -58,6 +58,7 @@ struct _timedelta
 struct _divisionquotienttimespan
 {
   int64_t quotient;          ///< Quotient of two timedeltas.
+  int64_t remainder_days;    ///< Remainder in days of two timedeltas division.
   int64_t remainder_in_ms;   ///< Remainder in milli seconds of two timedeltas division.
 };
 
@@ -86,14 +87,8 @@ struct _timedelta*
 julianDeltaToTimeDelta(struct _juliandelta* jd, struct _datetime *dt, struct _timedelta* td_return);
 
 struct _divisionquotienttimespan*
-divideTimeDeltaInSeconds(struct _timedelta* dividend, struct _timedelta* divisor, struct _datetime* base_dt, 
-struct _divisionquotienttimespan* quo_ret);
-
-struct _divisionquotienttimespan*
-divideTwoDatetimeDiffsInSeconds(struct _datetime* dt1_dividend, struct _datetime* dt2_dividend,struct _datetime* dt1_divisor, struct _datetime* dt2_divisor, int64_t * denominator_ret,  struct _divisionquotienttimespan* quo_ret);
-
-struct _divisionquotienttimespan*
-divideDatetimeDifferenceInSeconds(struct _datetime* dt1, struct _datetime* dt2, struct _timedelta* divisor, struct _divisionquotienttimespan* quo_ret);
+divideTimeDelta(struct _timedelta* dividend, struct _timedelta* divisor,
+		struct _datetime* base_dt, struct _divisionquotienttimespan* quo_ret);
 
 struct _timedelta*
 getTimeDeltaFromDate(struct _date*, struct _date*, struct _timedelta*);
diff --git a/src/Makefile.in b/src/Makefile.in
index 2cd42561af23e8f2eeebc8dc7e7d2e80dcb8a8f2..a8bad5311976e4fcada9ee9eecfe3e06c52b76e9 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 0000000000000000000000000000000000000000..1373ef3448b1bc49321121bb9cb2a5011d3be87d
--- /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 (1992).
+                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
diff --git a/src/libmtime.f90 b/src/libmtime.f90
index 6e9df0ec0f796272c29c002b400fe9975f8ece6a..5855591b8e9a141431224c429f6d235162ea3257 100644
--- a/src/libmtime.f90
+++ b/src/libmtime.f90
@@ -1092,9 +1092,7 @@ module mtime_timedelta
   public :: getPTStringFromMinutes
   public :: getPTStringFromHours
   public :: timeDeltaToJulianDelta
-  public :: divideTimeDeltaInSeconds
-  public :: divideTwoDatetimeDiffsInSeconds
-  public :: divideDatetimeDifferenceInSeconds
+  public :: divideTimeDelta
   public :: operator(+)
   public :: operator(-)
   public :: operator(*)
@@ -1730,7 +1728,7 @@ contains
   END SUBROUTINE timeDeltaToJulianDelta
 
   !>
-  !! @brief division by seconds.
+  !! @brief timedelta division.
   !!
   !! @param[in]  dividend
   !!         A pointer to type timedelta
@@ -1741,7 +1739,7 @@ contains
   !! @param[out]  quotient
   !!         A pointer to type divisionquotienttimespan
   !!
-  SUBROUTINE divideTimeDeltaInSeconds(dividend, divisor, base_dt, quotient, errna)!OK-UNTESTED.
+  SUBROUTINE divideTimeDelta(dividend, divisor, base_dt, quotient, errna)!OK-UNTESTED.
     type(timedelta), target, intent(in) :: dividend
     type(timedelta), target, intent(in) :: divisor
     TYPE(datetime),  TARGET, INTENT(in) :: base_dt
@@ -1749,66 +1747,11 @@ contains
     INTEGER, INTENT(out), optional :: errna
     type(c_ptr) :: dummy_ptr
     IF (PRESENT(errna)) errna = 0 ! FIXME: no_error
-    dummy_ptr = my_dividetimedeltainseconds(c_loc(dividend), c_loc(divisor), c_loc(base_dt), c_loc(quotient))
+    dummy_ptr = my_dividetimedelta(c_loc(dividend), c_loc(divisor), c_loc(base_dt), c_loc(quotient))
     IF (PRESENT(errna) .AND. .NOT. c_associated(dummy_ptr)) THEN
       errna = errna + 2  ! increment error number by 2, see below for an explanation.
     ENDIF
-  end subroutine divideTimeDeltaInSeconds
-  !>
-  !! @brief division of two differences in datetimes.
-  !!
-  !! @param  dt1_dividend, dt2_dividend, dt1_divisor, dt2_divisor
-  !!         Reference date (a pointer to struct _datetime).
-  !!
-  !! @param  intvlsec
-  !!         Interval given in seconds.
-  !!
-  !! @return result of division. NULL indicates error.
-  subroutine divideTwoDatetimeDiffsInSeconds(dt1_dividend,dt2_dividend, &
-      &                                      dt1_divisor, dt2_divisor,  &
-      &                                      denominator, quotient)
-    type(datetime), target, intent(in) :: dt1_dividend
-    type(datetime), target, intent(in) :: dt2_dividend
-    type(datetime), target, intent(in) :: dt1_divisor
-    type(datetime), target, intent(in) :: dt2_divisor
-    integer(c_int64_t), target, intent(out) :: denominator
-    type(divisionquotienttimespan), target, intent(out) :: quotient
-    type(c_ptr) :: dummy_ptr
-    dummy_ptr = my_dividetwodatetimediffsinseconds(c_loc(dt1_dividend), c_loc(dt2_dividend),  &
-        &                                                c_loc(dt1_divisor), c_loc(dt2_divisor),  &
-        &                                                c_loc(denominator), c_loc(quotient))
-  end subroutine divideTwoDatetimeDiffsInSeconds
-
-  !>
-  !! @brief division of an datetime interval by seconds.
-  !!
-  !! the datetime interval is calculated by dt1-dt2.
-  !!
-  !! @param[in]  dt1
-  !!         A pointer to type datetime
-  !!
-  !! @param[in]  dt2
-  !!         A pointer to type datetime
-  !!
-  !! @param[in]  divisor
-  !!         A pointer to type timedelta
-  !!
-  !! @param[out]  quotient
-  !!         A pointer to type divisionquotienttimespan
-  !!
-  subroutine divideDatetimeDifferenceInSeconds(dt1, dt2, divisor, quotient, errna)
-    type(datetime), target, intent(in) :: dt1
-    type(datetime), target, intent(in) :: dt2
-    type(timedelta), target, intent(in) :: divisor
-    type(divisionquotienttimespan), target, intent(out) :: quotient
-    INTEGER, INTENT(out), optional :: errna
-    type(c_ptr) :: dummy_ptr
-    IF (PRESENT(errna)) errna = 0 ! FIXME: no_error
-    dummy_ptr = my_dividedatetimedifferenceinseconds(c_loc(dt1), c_loc(dt2), c_loc(divisor), c_loc(quotient))
-    IF (PRESENT(errna) .AND. .NOT. c_associated(dummy_ptr)) THEN
-      errna = errna + 2  ! increment error number by 2, see below for an explanation.
-    ENDIF
-  end subroutine divideDatetimeDifferenceInSeconds
+  END SUBROUTINE divideTimeDelta
   !
 end module mtime_timedelta
 !>
diff --git a/src/libmtime_hl.f90 b/src/libmtime_hl.f90
index 3d37c4ef9a2f76085b8d927339c15332933e193f..e8514283e6bc98c7d1f57bd42f535a5df388acab 100644
--- a/src/libmtime_hl.f90
+++ b/src/libmtime_hl.f90
@@ -194,7 +194,7 @@ MODULE mtime_hl
     ! thrown.
     PROCEDURE :: toMilliSeconds            => t_timedelta_toMilliSeconds
 
-    PROCEDURE :: divideInSecondsBy         => t_timedelta_divideInSecondsBy
+    PROCEDURE :: divideBy                  => t_timedelta_divideBy
 
 
     ! --- overloaded operators
diff --git a/src/mtime_c_bindings.f90 b/src/mtime_c_bindings.f90
index 6f3da716f9c9d95d83d094732a85ae22acf0dcb2..0c53dbc3d98918986423aa7ffb74d169faf05fce 100644
--- a/src/mtime_c_bindings.f90
+++ b/src/mtime_c_bindings.f90
@@ -1,4 +1,4 @@
-module mtime_c_bindings
+MODULE mtime_c_bindings
   !
   use, intrinsic :: iso_c_binding, only: c_bool, c_int, c_char, c_null_char, c_ptr, c_int64_t, &
        &                                 c_float, c_double, c_associated
@@ -84,6 +84,7 @@ module mtime_c_bindings
   !
   type, bind(c) :: divisionquotienttimespan
     integer(c_int64_t) :: quotient;
+    integer(c_int64_t) :: remainder_days;
     integer(c_int64_t) :: remainder_in_ms;
   end type divisionquotienttimespan
   !
@@ -510,39 +511,15 @@ module mtime_c_bindings
       type(c_ptr), value :: jd
     end function my_timedeltatojuliandelta
     !
-    FUNCTION my_divideTimeDeltaInSeconds(dividend, divisor, base_dt, quotient) &
-RESULT(ret_quotient) BIND(c,name='divideTimeDeltaInSeconds') 
+    FUNCTION my_divideTimeDelta(dividend, divisor, base_dt, quotient)  &
+      RESULT(ret_quotient) BIND(c,name='divideTimeDelta')
       import :: c_ptr
-      type(c_ptr) :: ret_quotient 
+      type(c_ptr) :: ret_quotient
       type(c_ptr), value :: dividend
       type(c_ptr), value :: divisor
       type(c_ptr), value :: base_dt
       type(c_ptr), value :: quotient
-    end function my_divideTimeDeltaInSeconds
-    !
-    FUNCTION my_divideDatetimeDifferenceInSeconds(dt1, dt2, divisor, quotient) RESULT(ret_quotient) &
-         &                                                                     bind(c,name='divideDatetimeDifferenceInSeconds')
-      import :: c_ptr
-      type(c_ptr) :: ret_quotient
-      type(c_ptr), value :: dt1
-      type(c_ptr), value :: dt2
-      type(c_ptr), value :: divisor
-      type(c_ptr), value :: quotient
-    end function my_divideDatetimeDifferenceInSeconds
-    !
-    function my_divideTwoDatetimeDiffsInSeconds(dt1_dividend, dt2_dividend, &
-         &                                      dt1_divisor,  dt2_divisor, denominator, quotient) &
-         &                                       result(ret_quotient) &
-         &                                       bind(c,name='divideTwoDatetimeDiffsInSeconds')
-      import :: c_ptr
-      type(c_ptr) :: ret_quotient
-      type(c_ptr), value :: dt1_dividend
-      type(c_ptr), value :: dt2_dividend
-      type(c_ptr), value :: dt1_divisor
-      type(c_ptr), value :: dt2_divisor
-      type(c_ptr), value :: denominator
-      type(c_ptr), value :: quotient
-    end function my_divideTwoDatetimeDiffsInSeconds
+    END FUNCTION my_divideTimeDelta
     !
   end interface
   !
diff --git a/src/mtime_t_timedelta.inc b/src/mtime_t_timedelta.inc
index 139d14766f2e574364fd06ff1999583781ad1456..83c6ef86b420adf8d75b73c97554ba8d609dfa94 100644
--- a/src/mtime_t_timedelta.inc
+++ b/src/mtime_t_timedelta.inc
@@ -217,7 +217,7 @@
     CALL my_deallocatetimedelta(c_pointer) !deallocate timedelta
   END FUNCTION t_timedelta_toString
 
-  FUNCTION t_timedelta_divideInSecondsBy (this, divisor, referenceDateTime) RESULT(quotient)
+  FUNCTION t_timedelta_divideBy (this, divisor, referenceDateTime) RESULT(quotient)
     CLASS(t_timedelta), INTENT(in), target  :: this
     TYPE(t_timedelta),  INTENT(in), target  :: divisor
     TYPE(t_datetime),   INTENT(IN), target  :: referenceDateTime
@@ -225,20 +225,12 @@
     TYPE(t_datetime), target                :: dt_tmp
     type(c_ptr)                             :: dummy_ptr
 
-    dummy_ptr = my_divideTimeDeltaInSeconds(c_loc(this%td), c_loc(divisor%td), &
+    dummy_ptr = my_divideTimeDelta(c_loc(this%td), c_loc(divisor%td), &
                       & c_loc(referenceDateTime%dt), c_loc(quotient))
 
-    IF (.NOT. c_associated(dummy_ptr)) THEN
-        dt_tmp = referenceDateTime + this
-        dummy_ptr = my_dividedatetimedifferenceinseconds(c_loc(dt_tmp%dt),    &
-          &                     c_loc(referenceDateTime%dt), &
-          &                     c_loc(divisor%td), c_loc(quotient))
-
-        CALL handle_errno(.not. c_associated(dummy_ptr), general_arithmetic_error, &
-          __FILE__, &
-          __LINE__)
-    END IF
-  END FUNCTION t_timedelta_divideInSecondsBy
+    CALL handle_errno(.not. c_associated(dummy_ptr), general_arithmetic_error, &
+                      __FILE__, __LINE__)
+  END FUNCTION t_timedelta_divideBy
 
   FUNCTION t_timedelta_toSeconds (this, td) RESULT(seconds)
     CLASS(t_timedelta), INTENT(in) :: this
diff --git a/src/mtime_timedelta.c b/src/mtime_timedelta.c
index e665d5f1d7ffbe1bcf7a77d3eb42c56e346b3327..e1a549df6f8ebd8232f310ec97439c0871fd1cac 100644
--- a/src/mtime_timedelta.c
+++ b/src/mtime_timedelta.c
@@ -1,11 +1,11 @@
 /**
  * @file mtime_timedelta.c
- * 
+ *
  * @brief TimeDelta and some operations supported on TimeDelta.
- * 
+ *
  * @author Luis Kornblueh, Rahul Sinha. MPIM.
  * @date March 2013
- * 
+ *
  */
 
 #include <stdio.h>
@@ -14,6 +14,7 @@
 #include <inttypes.h>
 #include <time.h>
 #include <math.h>
+#include "int_div.h"
 
 #include "mtime_timedelta.h"
 
@@ -36,7 +37,7 @@
  *         A pointer to char. The string should contain parameters with which TimeDelta is to be created.
  *
  * @return td
- *         A pointer to a filled TimeDelta. 
+ *         A pointer to a filled TimeDelta.
  */
 
 struct _timedelta*
@@ -52,7 +53,7 @@ newTimeDelta(const char* tds)
 
       if (
            ((duration_type_flag = verify_string_duration(tds, isoDuration)) != DURATION_MATCH_STD)
-           && 
+           &&
            (duration_type_flag != DURATION_MATCH_LONG)
          )
         {
@@ -83,7 +84,7 @@ newTimeDelta(const char* tds)
 
       //Cleanup.
       deallocate_iso8601_duration(isoDuration);
-      isoDuration = NULL;       
+      isoDuration = NULL;
 
       return td;
     }
@@ -96,7 +97,7 @@ newTimeDelta(const char* tds)
  * @brief Construct new TimeDelta using 'raw' numerical values.
  *
  * @param  _sign
- *         An char value denoting positive('+') or negative('-') TimeDelta. 
+ *         An char value denoting positive('+') or negative('-') TimeDelta.
  * @param  _year
  *         An "int64_t" value denoting the year part of TimeDelta.
  * @param  _month
@@ -113,14 +114,14 @@ newTimeDelta(const char* tds)
  *         An "int" value denoting the milli-second part of TimeDelta.
  *
  * @return td
- *         A pointer to a filled TimeDelta. 
+ *         A pointer to a filled TimeDelta.
  */
 
 struct _timedelta *
 newRawTimeDelta(char _sign, int64_t _year, int _month, int _day, int _hour, int _minute, int _second, int _ms)
 {
   struct _timedelta* td;
-  
+
   td = (struct _timedelta *) calloc(1, sizeof(struct _timedelta));
   if (td == NULL )
     {
@@ -135,7 +136,7 @@ newRawTimeDelta(char _sign, int64_t _year, int _month, int _day, int _hour, int
     {
       td->flag_std_form = false;
     }
-  
+
   td->sign   = _sign;
   td->year   = _year;
   td->month  = _month;
@@ -144,7 +145,7 @@ newRawTimeDelta(char _sign, int64_t _year, int _month, int _day, int _hour, int
   td->minute = _minute;
   td->second = _second;
   td->ms     = _ms;
-  
+
   return td;
 }
 
@@ -152,10 +153,10 @@ newRawTimeDelta(char _sign, int64_t _year, int _month, int _day, int _hour, int
  * @brief Copy the values and construct a new TimeDelta.
  *
  * @param  td
- *         A pointer to struct _timedelta. Values of td are used to initialize the new timedelta being created. 
+ *         A pointer to struct _timedelta. Values of td are used to initialize the new timedelta being created.
  *
  * @return _td
- *         A pointer to an initialized TimeDelta object. 
+ *         A pointer to an initialized TimeDelta object.
  */
 
 struct _timedelta*
@@ -215,7 +216,7 @@ deallocateTimeDelta(struct _timedelta* td)
  *         A pointer to struct _timedelta.
  *
  * @return boolean
- *         if (td1 > td2), return greater_than. If (td1 == td2), return equal_to. If (td1 < td2), return less_than. 
+ *         if (td1 > td2), return greater_than. If (td1 == td2), return equal_to. If (td1 < td2), return less_than.
  *         Return compare_error indicating error.
  */
 
@@ -234,12 +235,12 @@ compareTimeDelta(struct _timedelta* td1, struct _timedelta* td2)
 
       compare_return_val boolean = compare_error;
 
-      /* 
+      /*
        * add special case for seconds larger than 59 by creating
-       * temporary/-ies, otherwise proceed with normal processing 
+       * temporary/-ies, otherwise proceed with normal processing
        */
-      if ( (td1->second < 59) && td2->second < 59)  
-	{	 
+      if ( (td1->second < 59) && td2->second < 59)
+	{
 	  if(td1->year == td2->year)
 	    {
 	      if(td1->month == td2->month)
@@ -281,7 +282,7 @@ compareTimeDelta(struct _timedelta* td1, struct _timedelta* td2)
 			  else if(td1->minute < td2->minute)
 			    {
 			      boolean = less_than;
-			    }       
+			    }
 			}
 		      else if(td1->hour > td2->hour)
 			{
@@ -298,7 +299,7 @@ compareTimeDelta(struct _timedelta* td1, struct _timedelta* td2)
 		    }
 		  else  if(td1->day < td2->day)
 		    {
-		      boolean = less_than;                    
+		      boolean = less_than;
 		    }
 		}
 	      else if (td1->month > td2->month)
@@ -327,7 +328,7 @@ compareTimeDelta(struct _timedelta* td1, struct _timedelta* td2)
 	      if ( (td1->month == 0) && (td2->month == 0))
 		{
 		  int second;
-		  
+
 		  td1_tmp.day    = td1->day;
 		  td1_tmp.hour   = td1->hour;
 		  td1_tmp.minute = td1->minute;
@@ -389,7 +390,7 @@ compareTimeDelta(struct _timedelta* td1, struct _timedelta* td2)
 			  else if(td1_tmp.minute < td2_tmp.minute)
 			    {
 			      boolean = less_than;
-			    }       
+			    }
 			}
 		      else if(td1_tmp.hour > td2_tmp.hour)
 			{
@@ -406,7 +407,7 @@ compareTimeDelta(struct _timedelta* td1, struct _timedelta* td2)
 		    }
 		  else  if(td1_tmp.day < td2_tmp.day)
 		    {
-		      boolean = less_than;                    
+		      boolean = less_than;
 		    }
 		}
 	      else
@@ -416,16 +417,16 @@ compareTimeDelta(struct _timedelta* td1, struct _timedelta* td2)
 	    }
 	  else
 	    {
-	      return compare_error;	      
+	      return compare_error;
 	    }
 	}
-      
+
       if (td1->sign == '+' && td2->sign == '+' )
 	return boolean;
       else if (td1->sign == '-' && td2->sign == '-' )
 	return boolean*(-1);
     }
-  
+
   return compare_error;
 }
 
@@ -435,7 +436,7 @@ compareTimeDelta(struct _timedelta* td1, struct _timedelta* td2)
  *
  * Routine replaceTimeDelta copies the contents of source TimeDelta into a Destination TimeDelta object.
  *
- * @param  tdsrc         
+ * @param  tdsrc
  *         A pointer to struct _timedelta. Copy "FROM" TimeDelta object.
  *
  * @param  tddest
@@ -469,14 +470,14 @@ else
 }
 
 
-static 
+static
 struct _juliandelta*
 localTimeDeltaToJulianDelta_NonStandardTimeDelta(struct _timedelta* td, struct _datetime* base_dt, struct _juliandelta* jd_return)
   {
 
     if ( td->sign == '-' )
       {
-        /* Negative TimeDelta.  A negative juliandelta is represented in the following 
+        /* Negative TimeDelta.  A negative juliandelta is represented in the following
            way: -P01DT00.500S  = jd2->sign = '-', jd2->day = -30, jd2->ms = -500. */
 
         jd_return->sign = '-';
@@ -564,7 +565,7 @@ localTimeDeltaToJulianDelta_StandardTimeDelta_CalType360OR365(struct _timedelta*
 
   if ( td->sign == '-' )
     {
-      /* Negative TimeDelta.  A negative juliandelta is represented in the following 
+      /* Negative TimeDelta.  A negative juliandelta is represented in the following
          way: -P01DT00.500S  = jd2->sign = '-', jd2->day = -30, jd2->ms = -500. */
 
       jd_return->sign = '-';
@@ -628,21 +629,21 @@ localTimeDeltaToJulianDelta_StandardTimeDelta_CalTypeGREGORIAN(struct _timedelta
 
           /* In final year, the crucial point is the month of february. */
           if (
-                ( 
-                        (testYearIsLeapYear(base_dt->date.year + td->year + 1) ) 
-                        && 
-                        (base_dt->date.month >= 3) 
+                (
+                        (testYearIsLeapYear(base_dt->date.year + td->year + 1) )
+                        &&
+                        (base_dt->date.month >= 3)
                 )
                 ||
-                ( 
-                        (testYearIsLeapYear(base_dt->date.year + td->year)) 
-                        && 
-                        (base_dt->date.month < 3) 
+                (
+                        (testYearIsLeapYear(base_dt->date.year + td->year))
+                        &&
+                        (base_dt->date.month < 3)
                 )
             )
             {
-              /* If the (base year + delta year) is a year before a leap year and base month is >= 3 
-               OR 
+              /* If the (base year + delta year) is a year before a leap year and base month is >= 3
+               OR
                base year + delta year is a leap year and month is < 3
                => An addition of leap-year specific delta for each month.
               */
@@ -671,9 +672,9 @@ localTimeDeltaToJulianDelta_StandardTimeDelta_CalTypeGREGORIAN(struct _timedelta
                         ((testYearIsLeapYear(i)) && (base_dt->date.month < 3))
                  )
                 {
-                  /* If the next year is a leap year and month is >= 3 OR 
-                   this year is a leap year and month is less than 3 
-                   => delta of 1 year corresponds to 366 day julian delta. 
+                  /* If the next year is a leap year and month is >= 3 OR
+                   this year is a leap year and month is less than 3
+                   => delta of 1 year corresponds to 366 day julian delta.
                   */
                   jd_return->day = jd_return->day + NO_OF_DAYS_IN_A_LEAP_YEAR;
                 }
@@ -684,17 +685,17 @@ localTimeDeltaToJulianDelta_StandardTimeDelta_CalTypeGREGORIAN(struct _timedelta
                 }
             }
 
-          jd_return->day = jd_return->day + msdinm[base_dt->date.month - 1][td->month] 
+          jd_return->day = jd_return->day + msdinm[base_dt->date.month - 1][td->month]
                                           + td->day;
 
-          jd_return->ms =  (int64_t)td->hour   * NO_OF_MS_IN_A_HOUR 
-                        +  (int64_t)td->minute * NO_OF_MS_IN_A_MINUTE 
-                        +  (int64_t)td->second * NO_OF_MS_IN_A_SECOND 
+          jd_return->ms =  (int64_t)td->hour   * NO_OF_MS_IN_A_HOUR
+                        +  (int64_t)td->minute * NO_OF_MS_IN_A_MINUTE
+                        +  (int64_t)td->second * NO_OF_MS_IN_A_SECOND
                         + td->ms;
         }
       else if ( td->sign == '-' )
         {
-          /* A negative juliandelta is represented in the following way: 
+          /* A negative juliandelta is represented in the following way:
              -P01DT00.500S  = jd2->sign = '-', jd2->day = -30, jd2->ms = -500.  */
 
           jd_return->sign = '-';
@@ -705,8 +706,8 @@ localTimeDeltaToJulianDelta_StandardTimeDelta_CalTypeGREGORIAN(struct _timedelta
                 ||
                 ((testYearIsLeapYear(base_dt->date.year - td->year)) && (base_dt->date.month >= 3)))
             {
-              /* If the (base year - delta year) is a year after leap year and base month is < 3 
-               OR 
+              /* If the (base year - delta year) is a year after leap year and base month is < 3
+               OR
                base year - delta year is a leap year and month is >= 3
                => A substraction of leap-year specific delta for each month.
               */
@@ -735,9 +736,9 @@ localTimeDeltaToJulianDelta_StandardTimeDelta_CalTypeGREGORIAN(struct _timedelta
                         ((testYearIsLeapYear(i)) && (base_dt->date.month >= 3))
                  )
                 {
-                  /* If the previous year is a leap year and month is < 3 OR 
-                   this year is a leap year and month is >= 3 
-                   => delta of 1 year corresponds to 366 day julian delta. 
+                  /* If the previous year is a leap year and month is < 3 OR
+                   this year is a leap year and month is >= 3
+                   => delta of 1 year corresponds to 366 day julian delta.
                    */
                   jd_return->day = jd_return->day - NO_OF_DAYS_IN_A_LEAP_YEAR;
                 }
@@ -750,9 +751,9 @@ localTimeDeltaToJulianDelta_StandardTimeDelta_CalTypeGREGORIAN(struct _timedelta
           jd_return->day = jd_return->day + (-1) * (ndiny - msdinm[base_dt->date.month - 1][NO_OF_MONTHS_IN_A_YEAR - td->month])
                                           + (-1) * td->day;
 
-          jd_return->ms = (-1) *  (int64_t)td->hour * NO_OF_MS_IN_A_HOUR 
-                        + (-1) *  (int64_t)td->minute * NO_OF_MS_IN_A_MINUTE 
-                        + (-1) *  (int64_t)td->second * NO_OF_MS_IN_A_SECOND 
+          jd_return->ms = (-1) *  (int64_t)td->hour * NO_OF_MS_IN_A_HOUR
+                        + (-1) *  (int64_t)td->minute * NO_OF_MS_IN_A_MINUTE
+                        + (-1) *  (int64_t)td->second * NO_OF_MS_IN_A_SECOND
                         + (-1) * td->ms;
 
         }
@@ -766,26 +767,26 @@ localTimeDeltaToJulianDelta_StandardTimeDelta_CalTypeGREGORIAN(struct _timedelta
 
 
 /* Internal function. */
-/* Converts TimeDelta value to a lib defined juliandelta value. Juliadelta depends on the calendar type. 
- Notice that TimeDelta is not uniquely defined but depends on the definition of corresponding DateTime 
+/* Converts TimeDelta value to a lib defined juliandelta value. Juliadelta depends on the calendar type.
+ Notice that TimeDelta is not uniquely defined but depends on the definition of corresponding DateTime
  (Referred to as Anchor date in the following).
 
- The library assumes the following definition: Let A denote an anchor date and P a timedelta. For a 
- positive P, A + P = B where date B > A. Consequently, for a negative P, A + P = B where A > B. 
+ The library assumes the following definition: Let A denote an anchor date and P a timedelta. For a
+ positive P, A + P = B where date B > A. Consequently, for a negative P, A + P = B where A > B.
  Also, when P is positive, a delta of 1 month has as many days as in the month of anchor DateTime;
- a delta of 2 months corresponds to the number of days in the anchor date month and the next month and 
+ a delta of 2 months corresponds to the number of days in the anchor date month and the next month and
  so on. When P is negative, a delta of 1 month corresponds to as many days as in the month before the
  anchor date month; a delta of 2 month corresponds to as many days as in the month before the
- anchor date month and the month before that and so on. 
- 
+ anchor date month and the month before that and so on.
+
  For eg, TimeDelta of P01M, when the 'Anchor' DateTime is 2001-02-01T00:00:00.000 is equivalent to a
- Julian Delta of positive 28 days ( #days in February ) and 0 ms. Also, TimeDelta of P02M, when the 
+ Julian Delta of positive 28 days ( #days in February ) and 0 ms. Also, TimeDelta of P02M, when the
  'Anchor' DateTime is 2001-02-01T00:00:00.000 is equivalent to a Julian Delta of positive 28+31 days
- ( #days in February PLUS #days in March) and 0 ms. Similarly, a TimeDelta of -P01M with 'Anchor' 
- DateTime of 2001-02-01T00:00:00.000 is equivalent to a Julian Delta of negative 31 days 
- ( #days in January ) and 0 ms. Likewise, TimeDelta of -P02M with 'Anchor' DateTime of 
+ ( #days in February PLUS #days in March) and 0 ms. Similarly, a TimeDelta of -P01M with 'Anchor'
+ DateTime of 2001-02-01T00:00:00.000 is equivalent to a Julian Delta of negative 31 days
+ ( #days in January ) and 0 ms. Likewise, TimeDelta of -P02M with 'Anchor' DateTime of
  2001-02-01T00:00:00.000 is equivalent to a Julian Delta of negative 31+31 days ( #days in January
- PLUS #days in December) and 0 ms. The same logic is used for all calendar types (with corresponding 
+ PLUS #days in December) and 0 ms. The same logic is used for all calendar types (with corresponding
  days in months according to current calendar type.)
  */
 
@@ -804,7 +805,7 @@ timeDeltaToJulianDelta(struct _timedelta* td, struct _datetime* base_dt, struct
           {
             case YEAR_OF_365_DAYS:
             case YEAR_OF_360_DAYS:
-        
+
             localTimeDeltaToJulianDelta_StandardTimeDelta_CalType360OR365(td, base_dt, jd_return);
             break;
 
@@ -821,7 +822,7 @@ timeDeltaToJulianDelta(struct _timedelta* td, struct _datetime* base_dt, struct
           }
 
       }
-    return jd_return;  
+    return jd_return;
   }
   else
     return NULL;
@@ -829,26 +830,26 @@ timeDeltaToJulianDelta(struct _timedelta* td, struct _datetime* base_dt, struct
 
 
 /* Internal function. */
-/* Converts a lib defined Julian delta value to a TimeDelta value. TimeDelta depends on the calendar type. 
- Notice that TimeDelta is not uniquely defined but depends on the definition of corresponding DateTime 
+/* Converts a lib defined Julian delta value to a TimeDelta value. TimeDelta depends on the calendar type.
+ Notice that TimeDelta is not uniquely defined but depends on the definition of corresponding DateTime
  (Referred to as Anchor date in the following).
 
- The library assumes the following definition: Let A denote an anchor date and P a timedelta. For a 
- positive P, A + P = B where date B > A. Consequently, for a negative P, A + P = B where A > B. 
+ The library assumes the following definition: Let A denote an anchor date and P a timedelta. For a
+ positive P, A + P = B where date B > A. Consequently, for a negative P, A + P = B where A > B.
  Also, when P is positive, a delta of 1 month has as many days as in the month of anchor DateTime;
- a delta of 2 months corresponds to the number of days in the anchor date month and the next month and 
+ a delta of 2 months corresponds to the number of days in the anchor date month and the next month and
  so on. When P is negative, a delta of 1 month corresponds to as many days as in the month before the
  anchor date month; a delta of 2 month corresponds to as many days as in the month before the
- anchor date month and the month before that and so on. 
+ anchor date month and the month before that and so on.
 
  For eg, TimeDelta of P01M, when the 'Anchor' DateTime is 2001-02-01T00:00:00.000 is equivalent to a
- Julian Delta of positive 28 days ( #days in February ) and 0 ms. Also, TimeDelta of P02M, when the 
+ Julian Delta of positive 28 days ( #days in February ) and 0 ms. Also, TimeDelta of P02M, when the
  'Anchor' DateTime is 2001-02-01T00:00:00.000 is equivalent to a Julian Delta of positive 28+31 days
- ( #days in February PLUS #days in March) and 0 ms. Similarly, a TimeDelta of -P01M with 'Anchor' 
- DateTime of 2001-02-01T00:00:00.000 is equivalent to a Julian Delta of negative 31 days 
- ( #days in January ) and 0 ms. Likewise, TimeDelta of -P02M with 'Anchor' DateTime of 
+ ( #days in February PLUS #days in March) and 0 ms. Similarly, a TimeDelta of -P01M with 'Anchor'
+ DateTime of 2001-02-01T00:00:00.000 is equivalent to a Julian Delta of negative 31 days
+ ( #days in January ) and 0 ms. Likewise, TimeDelta of -P02M with 'Anchor' DateTime of
  2001-02-01T00:00:00.000 is equivalent to a Julian Delta of negative 31+31 days ( #days in January
- PLUS #days in December) and 0 ms. The same logic is used for all calendar types (with corresponding 
+ PLUS #days in December) and 0 ms. The same logic is used for all calendar types (with corresponding
  days in months according to current calendar type).
 */
 
@@ -887,8 +888,8 @@ julianDeltaToTimeDelta(struct _juliandelta* jd, struct _datetime* base_dt, struc
 
       if ( jd->sign == '-' )
         {
-          /* Negative TimeDelta.  A negative juliandelta is represented in the following 
-             way: -P01DT00.500S  = jd2->sign = '-', jd2->day = -30, jd2->ms = -500. */  
+          /* Negative TimeDelta.  A negative juliandelta is represented in the following
+             way: -P01DT00.500S  = jd2->sign = '-', jd2->day = -30, jd2->ms = -500. */
 
           td_return->sign = '-';
 
@@ -1132,7 +1133,7 @@ julianDeltaToTimeDelta(struct _juliandelta* jd, struct _datetime* base_dt, struc
     }
   else if (jd->sign == '-')
     {
-      /* Negative TimeDelta.  A negative juliandelta is represented in the following 
+      /* Negative TimeDelta.  A negative juliandelta is represented in the following
          way: -P01DT00.500S  = jd2->sign = '-', jd2->day = -30, jd2->ms = -500. */
 
       /* Negative delta. */
@@ -1160,7 +1161,7 @@ julianDeltaToTimeDelta(struct _juliandelta* jd, struct _datetime* base_dt, struc
       td_return->minute = ((-1)*(int) jd->ms - td_return->hour * NO_OF_MS_IN_A_HOUR) / NO_OF_MS_IN_A_MINUTE;
       td_return->second = ((-1)*(int) jd->ms - td_return->hour * NO_OF_MS_IN_A_HOUR - td_return->minute * NO_OF_MS_IN_A_MINUTE) / NO_OF_MS_IN_A_SECOND;
       td_return->ms     =  (-1)*(int) jd->ms - td_return->hour * NO_OF_MS_IN_A_HOUR - td_return->minute * NO_OF_MS_IN_A_MINUTE - td_return->second * NO_OF_MS_IN_A_SECOND;
-    } 
+    }
   else
   return NULL; /* ERROR: Sign of julian delta not defined. */
 
@@ -1172,168 +1173,92 @@ else
 /*! \endcond */
 
 
-// WE REALLY NEED SOMETHING BETTER HERE!
 struct _divisionquotienttimespan*
-divideTimeDeltaInSeconds(struct _timedelta* dividend, struct _timedelta* divisor, struct _datetime* base_dt, struct _divisionquotienttimespan* quo_ret)
+divideTimeDelta(struct _timedelta* dividend, struct _timedelta* divisor,
+		struct _datetime* base_dt, struct _divisionquotienttimespan* quo_ret)
 {
-  if ((dividend != NULL) && (divisor != NULL) && (quo_ret != NULL))
-    {
-      if ((dividend->year == 0) && (divisor->year == 0) && (divisor->month == 0))
-        {
-
-	  struct _juliandelta* jd = newJulianDelta('+', 0, 0);
-	  if ( jd == NULL )
-	    return 0;
-	  jd = timeDeltaToJulianDelta(dividend, base_dt, jd);
-	  if ( jd == NULL )
-	    {
-	      deallocateJulianDelta(jd);
-	      return 0;
-	    }
+  struct _divisionquotienttimespan* ret = quo_ret;
+  struct _juliandelta *jd_dividend = NULL, *jd_divisor = NULL;
+
+  if ((dividend == NULL) || (divisor == NULL) || (quo_ret == NULL))  ret = NULL;
+
+  /* convert dividend and divisor into _juliandelta data type. */
+  if (ret != NULL) {
+    jd_dividend = newJulianDelta('+', 0, 0);
+    if ( jd_dividend == NULL )
+      ret = NULL;
+    else {
+      jd_dividend = timeDeltaToJulianDelta(dividend, base_dt, jd_dividend);
+      if ( jd_dividend == NULL )  ret = NULL;
+    }
+  }
 
-	  intmax_t numerator   = (intmax_t) (((int64_t) jd->day * NO_OF_MS_IN_A_DAY + jd->ms));
-          intmax_t denominator = (intmax_t) (((int64_t) divisor->day * 86400 + 
-					      divisor->hour * 3600 + 
-					      divisor->minute * 60 + 
-					      divisor->second ) * 1000 + 
-					     divisor->ms);
-
-	  deallocateJulianDelta(jd);
-          if (denominator == 0) /* Division by zero is illegal. */
-	     return NULL;
-
-          imaxdiv_t div = imaxdiv(numerator, denominator);
-           
-          quo_ret->quotient        = (int64_t) div.quot;          
-          quo_ret->remainder_in_ms = (int64_t) div.rem;
-          return quo_ret;
-        }
+  if (ret != NULL) {
+    jd_divisor = newJulianDelta('+', 0, 0);
+    if ( jd_divisor == NULL )
+      ret = NULL;
+    else {
+      jd_divisor = timeDeltaToJulianDelta(divisor, base_dt, jd_divisor);
+      if (( jd_divisor == NULL )  || ( jd_divisor->sign == '-' )) ret = NULL;
     }
-  return NULL;
-}
+  }
 
-/**
- * @brief division by an interval given in of seconds.
- *
- * @param  refdt
- *         Reference date (a pointer to struct _datetime).
- * @param  dt
- *         A pointer to struct _datetime.
- * @param  intvlsec
- *         Interval given in seconds.
- *
- * @return result of division. NULL indicates error.
- */
-struct _divisionquotienttimespan*
-divideDatetimeDifferenceInSeconds(struct _datetime* dt1, struct _datetime* dt2, struct _timedelta* divisor, struct _divisionquotienttimespan* quo_ret)
-{
-  if ((dt1 != NULL) && (dt2 != NULL) && (divisor != NULL) && (quo_ret != NULL))
-    {
-      if ((divisor->year == 0) && (divisor->month == 0))
-        {
-          struct _julianday* jd1 = newJulianDay(0, 0);
-          if (jd1 != NULL) jd1 = date2julian(dt1, jd1);
-
-          struct _julianday* jd2 = newJulianDay(0, 0);
-          if (jd2 != NULL) jd2 = date2julian(dt2, jd2);
-
-          intmax_t denominator = (intmax_t) (((int64_t) divisor->day  * 86400 + divisor->hour * 3600 +divisor->minute  * 60 + divisor->second ) * 1000 + divisor->ms);
-
-          if (denominator == 0) /* Division by zero is illegal. */
-            return NULL;
-
-          struct _juliandelta* jd = newJulianDelta('+', 0, 0);
-          jd = substractJulianDay(jd1, jd2, jd);  
-          intmax_t numerator = (intmax_t) (jd->day * 86400000 + jd->ms);
-          
-          imaxdiv_t div = imaxdiv(numerator, denominator);
-          
-          quo_ret->quotient        = (int64_t) div.quot;          
-          quo_ret->remainder_in_ms = (int64_t) div.rem;
-          
-          if (jd1 != NULL)  deallocateJulianDay(jd1);
-          if (jd2 != NULL)  deallocateJulianDay(jd2);
-          if (jd  != NULL)  deallocateJulianDelta(jd);
-          return quo_ret;
-        }
+  if (ret != NULL) {
+    /* perform large integer division. */
+    t_int days1, day_fraction1, days2, day_fraction2;
+    t_int q_decimal  = 0, remainder_days  = 0, remainder_ms  = 0;
+
+    days1           = (t_int) jd_dividend->day;
+    day_fraction1   = (t_int) jd_dividend->ms;
+    days2           = (t_int) jd_divisor->day;
+    day_fraction2   = (t_int) jd_divisor->ms;
+
+    int iret = divide_timespan(days1, day_fraction1, days2, day_fraction2,
+			       &q_decimal, &remainder_days, &remainder_ms);
+
+    quo_ret->quotient        = (int64_t) 0;
+    quo_ret->remainder_days  = (int64_t) 0;
+    quo_ret->remainder_in_ms = (int64_t) 0;
+    if (iret != 0)
+      ret = NULL;
+    else {
+      quo_ret->quotient        = (int64_t) q_decimal;
+      quo_ret->remainder_days  = (int64_t) remainder_days;
+      quo_ret->remainder_in_ms = (int64_t) remainder_ms;
     }
-  return NULL;
-}
+  }
 
-/**
- * @brief division of two differences in datetimes.
- *
- * @param  dt1_dividend, dt2_dividend, dt1_divisor, dt2_divisor
- *         Reference date (a pointer to struct _datetime).
- *
- * @param  intvlsec
- *         Interval given in seconds.
- *
- * @return result of division. NULL indicates error.
- */
-struct _divisionquotienttimespan*
-divideTwoDatetimeDiffsInSeconds(struct _datetime* dt1_dividend, struct _datetime* dt2_dividend, struct _datetime* dt1_divisor, struct _datetime* dt2_divisor, int64_t* denominator_ret,  struct _divisionquotienttimespan* quo_ret)
-{
-  if ((dt1_dividend != NULL) && (dt2_dividend != NULL) && 
-      (dt1_divisor != NULL) && (dt2_divisor != NULL) && (quo_ret != NULL))
-    {
-          // dividend
-          struct _julianday* jd1_dividend = newJulianDay(0, 0);
-          if (jd1_dividend != NULL) jd1_dividend = date2julian(dt1_dividend, jd1_dividend);
-
-          struct _julianday* jd2_dividend = newJulianDay(0, 0);
-          if (jd2_dividend != NULL) jd2_dividend = date2julian(dt2_dividend, jd2_dividend);
-
-          struct _juliandelta* jd_dividend = newJulianDelta('+', 0, 0);
-          jd_dividend = substractJulianDay(jd1_dividend, jd2_dividend, jd_dividend);  
-          intmax_t numerator = (intmax_t) (jd_dividend->day * 86400000 + jd_dividend->ms);
-
-          // divisor
-          struct _julianday* jd1_divisor = newJulianDay(0, 0);
-          if (jd1_divisor != NULL) jd1_divisor = date2julian(dt1_divisor, jd1_divisor);
-
-          struct _julianday* jd2_divisor = newJulianDay(0, 0);
-          if (jd2_divisor != NULL) jd2_divisor = date2julian(dt2_divisor, jd2_divisor);
-
-          struct _juliandelta* jd_divisor = newJulianDelta('+', 0, 0);
-          jd_divisor = substractJulianDay(jd1_divisor, jd2_divisor, jd_divisor);  
-          intmax_t denominator = (intmax_t) (jd_divisor->day * 86400000 + jd_divisor->ms);
-          
-          imaxdiv_t div = imaxdiv(numerator, denominator);
-          
-          quo_ret->quotient        = (int64_t) div.quot;          
-          quo_ret->remainder_in_ms = (int64_t) div.rem;
-          *denominator_ret = (int64_t) (denominator / 1000);
-          
-          if (jd1_dividend != NULL)  deallocateJulianDay(jd1_dividend);
-          if (jd2_dividend != NULL)  deallocateJulianDay(jd2_dividend);
-          if (jd_dividend  != NULL)  deallocateJulianDelta(jd_dividend);
-          
-          if (jd1_divisor != NULL)  deallocateJulianDay(jd1_divisor);
-          if (jd2_divisor != NULL)  deallocateJulianDay(jd2_divisor);
-          if (jd_divisor  != NULL)  deallocateJulianDelta(jd_divisor);
-          
-          return quo_ret;          
+  if (ret != NULL) {
+    if ( jd_dividend->sign == '-' ) {
+      /* correct sign. */
+      quo_ret->quotient        *= -1;
+      quo_ret->remainder_days  *= -1;
+      quo_ret->remainder_in_ms *= -1;
     }
-  else 
-    return NULL;
+  }
+
+  if (jd_dividend != NULL)  deallocateJulianDelta(jd_dividend);
+  if (jd_divisor  != NULL)  deallocateJulianDelta(jd_divisor);
+
+  return ret;
 }
 
+
 /**
  * @brief Get the TimeDelta between two Dates d1 and d2 as (d1-d2).
  *
  * Routine getTimeDeltaFromDate 'substracts' two Dates and returns the TimeDelta between
- * them. Internally, Dates are converted to DateTimes and then delta is calculated using 
+ * them. Internally, Dates are converted to DateTimes and then delta is calculated using
  * getTimeDeltaFromDateTime().
- * 
- * This routine  handles all supported Calendar types; i.e. the translation from Calendar date 
- * to Julian date and conversion from Julian Delta to normal TimeDetla is Calendar-type dependent. 
- * For eg. for Calendar type Gregorian, the TimeDelta between 2001-02-01 and 2001-01-01 will be 1 month. 
- * Similarly, for Calendar of type 360-Day-Calendar, the TimeDelta will be 1 month. It must be noted 
+ *
+ * This routine  handles all supported Calendar types; i.e. the translation from Calendar date
+ * to Julian date and conversion from Julian Delta to normal TimeDetla is Calendar-type dependent.
+ * For eg. for Calendar type Gregorian, the TimeDelta between 2001-02-01 and 2001-01-01 will be 1 month.
+ * Similarly, for Calendar of type 360-Day-Calendar, the TimeDelta will be 1 month. It must be noted
  * however, that the two dates differ by 31 and 30 days respectively.
  *
  * @param  d1
- *         A pointer to struct _date. 
+ *         A pointer to struct _date.
  *
  * @param  d2
  *         A pointer to struct _date.
@@ -1383,17 +1308,17 @@ else
  *
  * Routine getTimeDeltaFromDateTime 'substracts' two DateTime's and returns the TimeDelta between
  * them. Each datetime is converted to an equivalent Julian Date. Substraction is then performed
- * on Julian axis. The "Julian delta" is finally converted back to normal calendar delta. 
- * 
- * This routine handles all supported Calendar types; i.e. the translation from Calendar date 
- * to Julian date and conversion from Julian Delta to normal TimeDetla is Calendar-type dependent. 
- * For eg. for Calendar type Gregorian, the TimeDelta between 2001-02-01T00:00:00.000 and 
- * 2001-01-01T00:00:00.000 will be 1 month. Similarly, for Calendar of type 360-Day-Calendar, 
- * the TimeDelta will be 1 month. It must be noted however, that the two dates differ by 31 and 
- * 30 days respectively. 
+ * on Julian axis. The "Julian delta" is finally converted back to normal calendar delta.
+ *
+ * This routine handles all supported Calendar types; i.e. the translation from Calendar date
+ * to Julian date and conversion from Julian Delta to normal TimeDetla is Calendar-type dependent.
+ * For eg. for Calendar type Gregorian, the TimeDelta between 2001-02-01T00:00:00.000 and
+ * 2001-01-01T00:00:00.000 will be 1 month. Similarly, for Calendar of type 360-Day-Calendar,
+ * the TimeDelta will be 1 month. It must be noted however, that the two dates differ by 31 and
+ * 30 days respectively.
  *
  * @param  dt1
- *         A pointer to struct _datetime. 
+ *         A pointer to struct _datetime.
  *
  * @param  dt2
  *         A pointer to struct _datetime.
@@ -1409,7 +1334,7 @@ struct _timedelta*
 getTimeDeltaFromDateTime(struct _datetime* dt1, struct _datetime* dt2, struct _timedelta* td_return)
 {
   if ((dt1 != NULL )&& (dt2 != NULL) && (td_return != NULL) ){
- 
+
   /* Convert dt1 to Julian. */
   struct _julianday* jd1 = newJulianDay(0, 0);
   if (jd1 == NULL)
@@ -1455,9 +1380,9 @@ else
 /**
  * @brief Get total number of milliseconds in timedelta.
  *
- * Routine getTotalMilliSecondsTimeDelta returns the total number of milliseconds in TimeDelta. 
- * Notice that TimeDelta is not uniquely defined but depends on the definition of corresponding 
- * DateTime. TimeDelta is first converted to corresponding delta on the Julian axis. Julian delta 
+ * Routine getTotalMilliSecondsTimeDelta returns the total number of milliseconds in TimeDelta.
+ * Notice that TimeDelta is not uniquely defined but depends on the definition of corresponding
+ * DateTime. TimeDelta is first converted to corresponding delta on the Julian axis. Julian delta
  * is finally converted to the correct millisecond value.
  *
  * @param  td
@@ -1501,9 +1426,9 @@ else
 /**
  * @brief Get total number of seconds in timedelta.
  *
- * Routine getTotalSecondsTimeDelta returns the total number of seconds in TimeDelta. Notice that TimeDelta 
+ * Routine getTotalSecondsTimeDelta returns the total number of seconds in TimeDelta. Notice that TimeDelta
  * is not uniquely defined but depends on the definition of corresponding DateTime. Internally, number of seconds
- * is calculated by calling the routine getTotalMilliSecondsTimeDelta() and then converting the millisecond value 
+ * is calculated by calling the routine getTotalMilliSecondsTimeDelta() and then converting the millisecond value
  * to seconds by dividing it by 1000.
  *
  * @param  td
@@ -1554,7 +1479,7 @@ timedeltaToString(struct _timedelta* td, char* toStr)
     strcpy (toStr,"P");
   else
     return NULL; /* ERROR: TD sign not set. Should never happen. */
-  
+
   if (td->year != 0)
     sprintf(&(toStr[strlen(toStr)]),"%" PRIi64 "Y",td->year);
 
@@ -1613,14 +1538,14 @@ else
 /**
 * @brief Add timedelta to Date.
 *
-* Routine addTimeDeltaToDate adds a timedelta to a Date and returns the new Date. Both Date 
+* Routine addTimeDeltaToDate adds a timedelta to a Date and returns the new Date. Both Date
 * and TimeDetla are first converted to corresponding values on the Julian axis. Addition is performed on
-* the julian axis and the resulting Julian Date is converted back to the corrsponding Date. 
+* the julian axis and the resulting Julian Date is converted back to the corrsponding Date.
 *
-* The library assumes the following definition: Let A denote an anchor date and P a timedelta. For a 
-* positive P, A + P = B where date B > A. Consequently, for a negative P, A + P = B where A > B. 
+* The library assumes the following definition: Let A denote an anchor date and P a timedelta. For a
+* positive P, A + P = B where date B > A. Consequently, for a negative P, A + P = B where A > B.
 * Also, when P is positive, a delta of 1 month has as many days as in the month of anchor DateTime;
-* a delta of 2 months corresponds to the number of days in the anchor date month and the next month and 
+* a delta of 2 months corresponds to the number of days in the anchor date month and the next month and
 * so on. When P is negative, a delta of 1 month corresponds to as many days as in the month before the
 * anchor date month; a delta of 2 month corresponds to as many days as in the month before the
 * anchor date month and the month before that and so on.
@@ -1676,14 +1601,14 @@ addTimeDeltaToDate(struct _date* d, struct _timedelta* td, struct _date* d_retur
 /**
  * @brief Add timedelta to DateTime.
  *
- * Routine addTimeDeltaToDateTime adds a timedelta to a DateTime and returns the new DateTime. Both DateTime 
+ * Routine addTimeDeltaToDateTime adds a timedelta to a DateTime and returns the new DateTime. Both DateTime
  * and TimeDetla are first converted to corresponding values on the Julian axis. Addition is performed on
  * the julian axis and the resulting Julian Date is converted back to the corrsponding DateTime.
  *
- * The library assumes the following definition: Let A denote an anchor date and P a timedelta. For a 
- * positive P, A + P = B where date B > A. Consequently, for a negative P, A + P = B where A > B. 
+ * The library assumes the following definition: Let A denote an anchor date and P a timedelta. For a
+ * positive P, A + P = B where date B > A. Consequently, for a negative P, A + P = B where A > B.
  * Also, when P is positive, a delta of 1 month has as many days as in the month of anchor DateTime;
- * a delta of 2 months corresponds to the number of days in the anchor date month and the next month and 
+ * a delta of 2 months corresponds to the number of days in the anchor date month and the next month and
  * so on. When P is negative, a delta of 1 month corresponds to as many days as in the month before the
  * anchor date month; a delta of 2 month corresponds to as many days as in the month before the
  * anchor date month and the month before that and so on.
@@ -1758,10 +1683,10 @@ addTimeDeltaToDateTime(struct _datetime* dt, struct _timedelta* td, struct _date
  * @brief Get the timedelta between current_dt and start_dt plus next integral-multiple-of-timestep (timedelta).
  *
  * Routine moduloTimeDeltaFromDateTime returns the timedelta between the current DateTime (current_dt) and the event's next-trigger time.
- * The next trigger time is defined as the the Anchor DateTime (start_dt) + N * TimeDelta(timestep) 
- * where N is the minimum positive integer for which this sum is >= Current DateTime. In case 
+ * The next trigger time is defined as the the Anchor DateTime (start_dt) + N * TimeDelta(timestep)
+ * where N is the minimum positive integer for which this sum is >= Current DateTime. In case
  * Anchor DateTime > Current DateTime, TimeDelta is calculated as start_dt - current_dt.
- * 
+ *
  * Notice that this TimeDelta will always be positive.
  *
  * @param  start_dt
@@ -1774,7 +1699,7 @@ addTimeDeltaToDateTime(struct _datetime* dt, struct _timedelta* td, struct _date
  *         A pointer to struct _datetime. The Current Date time.
  *
  * @param  modulo_td
- *         A pointer to struct _timedelta. The timedelta between 'current datetime' and 'Start Datetime plus next 
+ *         A pointer to struct _timedelta. The timedelta between 'current datetime' and 'Start Datetime plus next
  *         integral-multiple-of-timestep' is copied here.
  *
  * @return modulo_td
@@ -1832,7 +1757,7 @@ moduloTimeDeltaFromDateTime(struct _datetime* start_dt, struct _timedelta* times
  *         A pointer to struct _timedelta. The element-wise product of lambda and base_td is stored here.
  *
  * @return scaled_td
- *        A pointer to struct _timedelta. The filled structure containing the scaled timedelta values.    
+ *        A pointer to struct _timedelta. The filled structure containing the scaled timedelta values.
  */
 struct _timedelta*
 elementwiseScalarMultiplyTimeDelta(struct _timedelta* base_td, int64_t lambda, struct _timedelta* scaled_td)
@@ -1865,14 +1790,14 @@ elementwiseScalarMultiplyTimeDelta(struct _timedelta* base_td, int64_t lambda, s
           scaled_td->ms -= NO_OF_MS_IN_A_SECOND;
           scaled_td->second += 1;
         }
- 
+
       scaled_td->second += (int) lambda*base_td->second;
       while ( scaled_td->second >= 60 )
-        {       
+        {
           scaled_td->second -= 60;
           scaled_td->minute += 1;
         }
-   
+
       scaled_td->minute += (int) lambda*base_td->minute;
       while ( scaled_td->minute >= 60 )
         {
@@ -1894,7 +1819,7 @@ elementwiseScalarMultiplyTimeDelta(struct _timedelta* base_td, int64_t lambda, s
       scaled_td->year   = 0;
 
       return scaled_td;
-        
+
     }
   else
     return NULL;
@@ -1922,7 +1847,7 @@ elementwiseScalarMultiplyTimeDelta(struct _timedelta* base_td, int64_t lambda, s
  *         A pointer to struct _timedelta. The element-wise product of lambda and base_td is stored here.
  *
  * @return scaled_td
- *        A pointer to struct _timedelta. The filled structure containing the scaled timedelta values.    
+ *        A pointer to struct _timedelta. The filled structure containing the scaled timedelta values.
  */
 struct _timedelta*
 elementwiseScalarMultiplyTimeDeltaDP(struct _timedelta* base_td, double lambda, struct _timedelta* scaled_td)
@@ -1972,8 +1897,8 @@ elementwiseScalarMultiplyTimeDeltaDP(struct _timedelta* base_td, double lambda,
  * @brief Return the element-wise sum of two timedeltas.
  *
  * elementwiseAddTimeDeltatoTimeDelta adds two timedeltas elementwise and returns the result.
- * Timedeltas being added must be of the same sign; Substraction is not supported. 
- * The timedelta can not have days,months or years however: Only Timedeltas upto hours should call this routine. 
+ * Timedeltas being added must be of the same sign; Substraction is not supported.
+ * The timedelta can not have days,months or years however: Only Timedeltas upto hours should call this routine.
  * Also td_return->hour >= 24 will lead to an error.
  *
  *
@@ -1987,7 +1912,7 @@ elementwiseScalarMultiplyTimeDeltaDP(struct _timedelta* base_td, double lambda,
  *         A pointer to struct _timedelta. The element-wise sum of td1 and td2 is stored here.
  *
  * @return td_return
- *         A pointer to struct _timedelta. The filled structure containing the added timedelta values.   
+ *         A pointer to struct _timedelta. The filled structure containing the added timedelta values.
  */
 struct _timedelta*
 elementwiseAddTimeDeltatoTimeDelta(struct _timedelta* td1, struct _timedelta* td2,  struct _timedelta* td_return)
@@ -2000,13 +1925,13 @@ elementwiseAddTimeDeltatoTimeDelta(struct _timedelta* td1, struct _timedelta* td
 
       /*Reset td_return to 0.*/
       memset(td_return,0,sizeof(struct _timedelta));
-      
+
       if(td1->sign == td2->sign)
         {
           /* If signs match, do add. */
           td_return->sign = td1->sign;
 
-          td_return->ms += (td1->ms + td2->ms);  
+          td_return->ms += (td1->ms + td2->ms);
           if ( td_return->ms >= NO_OF_MS_IN_A_SECOND )
             {
               td_return->ms -= NO_OF_MS_IN_A_SECOND;
@@ -2037,22 +1962,22 @@ elementwiseAddTimeDeltatoTimeDelta(struct _timedelta* td1, struct _timedelta* td
           td_return->day        = 0;
           td_return->month      = 0;
           td_return->year       = 0;
-        } 
+        }
       else
         {
           /* Substraction not supported. */
-          return NULL; 
+          return NULL;
         }
 
       return td_return;
-    } 
+    }
   else
     return NULL;
 
 }
 
 /**
- * @brief Returns the remainder of timedelta a modulo timedelta p  
+ * @brief Returns the remainder of timedelta a modulo timedelta p
  *
  * moduloTimedelta(struct _timedelta *a, struct _timdelta *p) returns the remainder of
  * a modulo p. The function is restricted to input timedeltas less than 29 days.
@@ -2066,7 +1991,7 @@ elementwiseAddTimeDeltatoTimeDelta(struct _timedelta* td1, struct _timedelta* td
  * @param  p
  *         on return contains the quotient
  *
- * @return rem 
+ * @return rem
  *         returns remainder of the division of the numerator a by
  *         the denominator p. quotient * p + rem shall equal a.
  */
@@ -2075,15 +2000,15 @@ int64_t
 moduloTimedelta(struct _timedelta *a, struct _timedelta *p, int64_t *quot)
 {
   struct _datetime *dt = newDateTime("0001-01-01");
-    
+
   struct _juliandelta *jdd1 = newJulianDelta('+', (int64_t) 0, (int64_t) 0);
   struct _juliandelta *jdd2 = newJulianDelta('+', (int64_t) 0, (int64_t) 0);
-        
+
   jdd1 = timeDeltaToJulianDelta(a, dt, jdd1);
   jdd2 = timeDeltaToJulianDelta(p, dt, jdd2);
 
   int64_t d1, d2;
-  
+
   d1 = (int64_t) 86400000 * jdd1->day + jdd1->ms;
   d2 = (int64_t) 86400000 * jdd2->day + jdd2->ms;
 
@@ -2098,7 +2023,7 @@ moduloTimedelta(struct _timedelta *a, struct _timedelta *p, int64_t *quot)
 /**
  * @brief Return a PT String corresponding to arbitrary number of milliseconds.
  *
- * getPTStringFromMS() translates ms values to ISO 8601 compliant timedelta string. 
+ * getPTStringFromMS() translates ms values to ISO 8601 compliant timedelta string.
  * Conversion of ms >= 86400000 and  ms <= -86400000 not supported.
  *
  * @param  _ms
@@ -2108,16 +2033,16 @@ moduloTimedelta(struct _timedelta *a, struct _timedelta *p, int64_t *quot)
  *         A pointer to char. Translated string is written here.
  *
  * @return PTstr
- *         A pointer to char. The translated TimeDelta string.   
+ *         A pointer to char. The translated TimeDelta string.
  */
 
 char*
 getPTStringFromMS(int64_t _ms, char* PTstr)
 {
    /* Reuse the juliandelta to TimeDelta conversion routine. */
-   
+
    /* Create a _juliandelta object and copy the _ms to the jd->ms. */
-   struct _juliandelta* jd = NULL; 
+   struct _juliandelta* jd = NULL;
    /* if ((_ms >= NO_OF_MS_IN_A_DAY) || (_ms <= ((-1)*NO_OF_MS_IN_A_DAY))) */
    /*   { */
    /*     /\* ERROR: Conversion greater than 23:59:59:999 not supported. *\/ */
@@ -2125,17 +2050,17 @@ getPTStringFromMS(int64_t _ms, char* PTstr)
    /*   } */
 
    if (_ms >= 0)
-     jd = newJulianDelta('+', 0, _ms); 
+     jd = newJulianDelta('+', 0, _ms);
    else
      jd = newJulianDelta('-', 0, _ms);
 
    /* Create dummy variables for julianDeltaToTimeDelta() */
    struct _datetime* dumm_base_dt       = newDateTime(initDummyDTString);
    struct _timedelta* dummy_td_return   = newTimeDelta(initDummyTDString);
-   
+
    /* Get the translated TimeDelta and return the corresponding string. */
    PTstr = timedeltaToString(julianDeltaToTimeDelta(jd, dumm_base_dt, dummy_td_return), PTstr);
-     
+
    deallocateDateTime(dumm_base_dt);
    deallocateTimeDelta(dummy_td_return);
    if(jd)
@@ -2148,7 +2073,7 @@ getPTStringFromMS(int64_t _ms, char* PTstr)
 /**
  * @brief Return a PT String corresponding to arbitrary number of seconds.
  *
- * getPTStringFromSeconds() translates second values to ISO 8601 compliant timedelta string. 
+ * getPTStringFromSeconds() translates second values to ISO 8601 compliant timedelta string.
  * Conversion of s >= 86400 and  s <= -86400 not supported.
  *
  * @param  _s
@@ -2158,7 +2083,7 @@ getPTStringFromMS(int64_t _ms, char* PTstr)
  *         A pointer to char. Translated string is written here.
  *
  * @return PTstr
- *         A pointer to char. The translated TimeDelta string.   
+ *         A pointer to char. The translated TimeDelta string.
  */
 
 char*
@@ -2171,7 +2096,7 @@ getPTStringFromSeconds(int64_t _s, char* PTstr)
 /**
  * @brief Return a PT String corresponding to arbitrary number of seconds.
  *
- * getPTStringFromSecondsFloat() translates second values to ISO 8601 compliant timedelta string. 
+ * getPTStringFromSecondsFloat() translates second values to ISO 8601 compliant timedelta string.
  * Conversion of s >= 86400 and  s <= -86400 not supported. Returned PT string has a precision of 3
  * digits.
  *
@@ -2182,7 +2107,7 @@ getPTStringFromSeconds(int64_t _s, char* PTstr)
  *         A pointer to char. Translated string is written here.
  *
  * @return PTstr
- *         A pointer to char. The translated TimeDelta string.   
+ *         A pointer to char. The translated TimeDelta string.
  */
 
 char*
@@ -2195,7 +2120,7 @@ getPTStringFromSecondsFloat(float _s, char* PTstr)
 /**
  * @brief Return a PT String corresponding to arbitrary number of seconds.
  *
- * getPTStringFromSecondsDouble() translates second values to ISO 8601 compliant timedelta string. 
+ * getPTStringFromSecondsDouble() translates second values to ISO 8601 compliant timedelta string.
  * Conversion of s >= 86400 and  s <= -86400 not supported. Returned PT string has a precision of 3 digits.
  *
  * @param  _s
@@ -2205,7 +2130,7 @@ getPTStringFromSecondsFloat(float _s, char* PTstr)
  *         A pointer to char. Translated string is written here.
  *
  * @return PTstr
- *         A pointer to char. The translated TimeDelta string.   
+ *         A pointer to char. The translated TimeDelta string.
  */
 
 char*
@@ -2218,7 +2143,7 @@ getPTStringFromSecondsDouble(double _s, char* PTstr)
 /**
  * @brief Return a PT String corresponding to arbitrary number of minutes.
  *
- * getPTStringFromMinutes() translates minutes values to ISO 8601 compliant timedelta string. 
+ * getPTStringFromMinutes() translates minutes values to ISO 8601 compliant timedelta string.
  * Conversion of m >= 1440 and  m <= -1440 not supported.
  *
  * @param  _m
@@ -2228,7 +2153,7 @@ getPTStringFromSecondsDouble(double _s, char* PTstr)
  *         A pointer to char. Translated string is written here.
  *
  * @return PTstr
- *         A pointer to char. The translated TimeDelta string.   
+ *         A pointer to char. The translated TimeDelta string.
  */
 
 char*
@@ -2241,7 +2166,7 @@ getPTStringFromMinutes(int64_t _m, char* PTstr)
 /**
  * @brief Return a PT String corresponding to arbitrary number of Hours.
  *
- * getPTStringFromHours() translates hour values to ISO 8601 compliant timedelta string. 
+ * getPTStringFromHours() translates hour values to ISO 8601 compliant timedelta string.
  * Conversion of h >= 24 and  ms <= -24 not supported.
  *
  * @param  _h
@@ -2251,7 +2176,7 @@ getPTStringFromMinutes(int64_t _m, char* PTstr)
  *         A pointer to char. Translated string is written here.
  *
  * @return PTstr
- *         A pointer to char. The translated TimeDelta string.   
+ *         A pointer to char. The translated TimeDelta string.
  */
 
 char*