diff --git a/CHANGES b/CHANGES
index c882168b3e..9fb37caa22 100644
--- a/CHANGES
+++ b/CHANGES
@@ -1,3 +1,16 @@
+
+4.1.0-dev.646 | 2021-05-18 11:47:25 -0700
+
+  * Omit unneeded decimal points in modp_dtoa2() scientific notation output (Jon Siwek, Corelight)
+
+    For example, "1e-13" is now used instead of "1.e-13".
+
+  * GH-1244: Change modp_dtoa2() to use scientific notation for small values (Jon Siwek, Corelight)
+
+    This fixes problems where printing floating point numbers less than
+    10^-6 output as "0.0".  Such numbers now use using scientific notation
+    and preserve the value's actual floating point representation.
+
 4.1.0-dev.643 | 2021-05-17 11:57:58 -0700
 
   * GH-1546: Make DictIterator() public, add copy/move operators (Tim Wojtulewicz, Corelight)
diff --git a/VERSION b/VERSION
index 4b012da258..c76b75482e 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-4.1.0-dev.643
+4.1.0-dev.646
diff --git a/src/modp_numtoa.c b/src/modp_numtoa.c
index d77b68580a..f6b7059d8a 100644
--- a/src/modp_numtoa.c
+++ b/src/modp_numtoa.c
@@ -23,6 +23,8 @@
  */
 static const double _pow10[] = {1, 10, 100, 1000, 10000, 100000, 1000000,
                                10000000, 100000000, 1000000000};
+static const double _pow10r[] = {1, .1, .01, .001, .0001, .00001, .000001,
+                                .0000001, .00000001, .000000001};
 
 static void strreverse(char* begin, char* end)
 {
@@ -53,6 +55,7 @@ static void sn_strip_trailing_zeros(char* str)
 	if ( ! frac )
 		return;
 
+	char* start_dec = frac;
 	char* exp = 0;
 	char* trailing_zeros = 0;
 
@@ -78,6 +81,9 @@ static void sn_strip_trailing_zeros(char* str)
 		++frac;
 		}
 
+	if ( trailing_zeros == start_dec )
+		--trailing_zeros;
+
 	if ( trailing_zeros && exp )
 		{
 		for ( ; ; )
@@ -287,6 +293,14 @@ void modp_dtoa2(double value, char* str, int prec)
         prec = 9;
     }
 
+    double smallest = _pow10r[prec];
+
+    if (value != 0.0 && value < smallest) {
+        sprintf(str, "%.*e", DBL_DECIMAL_DIG - 1, neg ? -value : value);
+        sn_strip_trailing_zeros(str);
+        return;
+    }
+
     int whole = (int) value;
     double tmp = (value - whole) * _pow10[prec];
     uint32_t frac = (uint32_t)(tmp);
diff --git a/src/modp_numtoa.h b/src/modp_numtoa.h
index 619a4ac767..c8c46827c8 100644
--- a/src/modp_numtoa.h
+++ b/src/modp_numtoa.h
@@ -91,6 +91,10 @@ void modp_dtoa(double value, char* buf, int precision);
  * will be switched exponential format and include as many precision digits
  * as needed to preserve information.
  *
+ * If a non-zero input value is less than 10^(-precision), the output format
+ * will be switched exponential format and include as many precision digits
+ * as needed to preserve information.
+ *
  * \param[in] value
  * \param[out] buf  The allocated output buffer.  Should be 32 chars or more.
  * \param[in] precision  Number of digits to the right of the decimal point.
diff --git a/testing/btest/Baseline/language.double/out b/testing/btest/Baseline/language.double/out
index 8ec12b40e1..143e1a7c81 100644
--- a/testing/btest/Baseline/language.double/out
+++ b/testing/btest/Baseline/language.double/out
@@ -27,3 +27,35 @@ relational operator (PASS)
 relational operator (PASS)
 division operator (PASS)
 max double value = 1.7976931348623157e+308 (PASS)
+
+4.9999999999999999e-13
+4.9999999999999997e-12
+5.0000000000000002e-11
+5.0000000000000003e-10
+5.0000000000000001e-09
+4.9999999999999998e-08
+4.9999999999999998e-07
+
+0.000005
+0.00005
+0.0005
+0.005
+0.05
+0.5
+5.0
+
+1e-13
+9.9999999999999998e-13
+9.9999999999999994e-12
+1e-10
+1.0000000000000001e-09
+1e-08
+9.9999999999999995e-08
+
+0.000001
+0.00001
+0.0001
+0.001
+0.01
+0.1
+1.0
diff --git a/testing/btest/language/double.zeek b/testing/btest/language/double.zeek
index 56ce711da2..f98abc839a 100644
--- a/testing/btest/language/double.zeek
+++ b/testing/btest/language/double.zeek
@@ -75,5 +75,40 @@ event zeek_init()
 	local str1 = fmt("max double value = %.16e", d19);
 	test_case( str1, str1 == "max double value = 1.7976931348623157e+308" );
 
+	# Printing small numbers: default precision is 6 with values smaller than
+	# 10^-6 rendered in scientific notation, preserving exact floating point
+	# representation.
+	print "";
+	print 0.0000000000005;
+	print 0.000000000005;
+	print 0.00000000005;
+	print 0.0000000005;
+	print 0.000000005;
+	print 0.00000005;
+	print 0.0000005;
+	print "";
+	print 0.000005;
+	print 0.00005;
+	print 0.0005;
+	print 0.005;
+	print 0.05;
+	print 0.5;
+	print 5.0;
+	print "";
+	print 0.0000000000001;
+	print 0.000000000001;
+	print 0.00000000001;
+	print 0.0000000001;
+	print 0.000000001;
+	print 0.00000001;
+	print 0.0000001;
+	print "";
+	print 0.000001;
+	print 0.00001;
+	print 0.0001;
+	print 0.001;
+	print 0.01;
+	print 0.1;
+	print 1.0;
 }