addition sum of two positiv double numbers

src/addition.c

@@ -145,3 +145,32 @@ unsigned long long addition_precision_double(unsigned long long p1, unsigned lon
     unsigned int vornep1 = p1 >> 32, hintenp1 = p1 & hinten, vornep2 = p2 >> 32, hintenp2 = p2 & hinten;
     return ((unsigned long long) addition_uint(vornep1, vornep2) << 32) + (unsigned long long) addition_uint(hintenp1, hintenp2);
 }
+
+// addition of two double numbers
+
+double addition_double(double number1, double number2) {
+    if (number2 == 0.0) return number1;
+    else if (number1 == 0.0) return number2;
+    else if (number2 > number1) return addition_double(number2, number1);
+
+    unsigned long long e = 0x0010000000000000;
+    struct datad num1, num2, num3;
+    num1.number.dnum = number1, num2.number.dnum = number2;
+    num1.sign = sign_double(num1.number.lnum), num2.sign = sign_double(num2.number.lnum);
+    num1.exponent = exponent_double(num1.number.lnum), num2.exponent = exponent_double(num2.number.lnum);
+    num1.precision = precision_double(num1.number.lnum) | e, num2.precision = precision_double(num2.number.lnum) | e;
+
+    unsigned int k = (num1.exponent - num2.exponent > 52 ? 53 : num1.exponent - num2.exponent);
+    num3.precision = addition_precision_double(num1.precision, num2.precision >> k);
+    if (num3.precision > 0x001fffffffffffff) {
+        num3.exponent = addition_uint((unsigned int) num1.exponent, 1);
+        num3.precision = (num3.precision >> 1);
+    }
+    else {
+        num3.exponent = num1.exponent;
+    }
+    num3.precision ^= e;
+    num3.number.lnum = output_double(num1.sign, num3.exponent, num3.precision);
+
+    return num3.number.dnum;
+}

src/addition.h

@@ -33,4 +33,6 @@ unsigned long long output_double(unsigned long long sign, unsigned long long exp
 
 unsigned long long addition_precision_double(unsigned long long p1, unsigned long long p2);
 
+double addition_double(double number1, double number2);
+
 #endif // ADDITION_H