--- a/QDouble.st	Sun Jun 18 23:20:39 2017 +0200
+++ b/QDouble.st	Mon Jun 19 00:16:48 2017 +0200
@@ -1,5 +1,3 @@
-"{ Encoding: utf8 }"
-
 "
  COPYRIGHT (c) 2017 by eXept Software AG
 	      All Rights Reserved
@@ -24,590 +22,14 @@
 !
 
 !QDouble primitiveDefinitions!
-%{
-#include <stdio.h>
-#include <errno.h>
-
-#define __USE_ISOC9X 1
-#define __USE_ISOC99 1
-#include <math.h>
-
-/*
- * on some systems, errno is a macro ... check for it here
- */
-#ifndef errno
- extern errno;
-#endif
-
-#if !defined (__win32__)
-# include <locale.h>
-#endif
-
-#if defined (__aix__)
-# include <float.h>
-#endif
-
-#if defined(__irix__) || defined(__solaris__) || defined(__sunos__)
-# include <nan.h>
-#endif
-
-#if defined(__linux__)
-# ifndef NAN
-#  include <bits/nan.h>
-# endif
-#endif
-
-#ifdef __win32__
-# ifndef isinf
-#  define isinf(x) \
-	((((unsigned int *)(&x))[0] == 0x00000000) && \
-	 ((((unsigned int *)(&x))[1] & 0x7FF00000) == 0x7FF00000))
-# endif
-#endif
-
-#if defined(__x86__) || defined(__x86_64__)
-
-# ifndef _FPU_EXTENDED
-#  define _FPU_EXTENDED 0x0300
-# endif
-
-# ifndef _FPU_DOUBLE
-#  define _FPU_DOUBLE 0x0200
-# endif
-
-# if defined( __win32__ ) && (defined( __BORLANDC__ ) || defined( __VISUALC__ ))
-
-#  define fpu_fix_start(old_cw_ptr)\
-    {\
-	*old_cw_ptr = _control87(0, 0); \
-	_control87(_FPU_DOUBLE, _FPU_EXTENDED);\
-    }
-
-#  define fpu_fix_end(old_cw_ptr)\
-    {\
-	_control87(*old_cw_ptr, _FPU_EXTENDED);\
-    }
-
-# else // assume MINGW, GCC or CLANG
-
-#  ifndef _FPU_GETCW
-#   define _FPU_GETCW(x) asm volatile ("fnstcw %0":"=m" (x));
-#  endif
-#  ifndef _FPU_SETCW
-#   define _FPU_SETCW(x) asm volatile ("fldcw %0": :"m" (x));
-#  endif
-
-#  define fpu_fix_start(old_cw_ptr)\
-    {\
-	volatile unsigned short cw, new_cw;\
-	_FPU_GETCW(cw);\
-	new_cw = (cw & ~_FPU_EXTENDED) | _FPU_DOUBLE;\
-	_FPU_SETCW(new_cw);\
-	*old_cw_ptr = cw;\
-    }
-
-#  define fpu_fix_end(old_cw_ptr)\
-    {\
-	volatile unsigned short cw = *old_cw_ptr;\
-	_FPU_SETCW(cw);\
-    }
-
-# endif
-
-#endif
-
-
-struct qd_real {
-    double x[4];    /* The Components. */
-};
-
-struct __quadDoubleStruct {
-	STX_OBJ_HEADER
-#ifdef __NEED_DOUBLE_ALIGN
-	__FILLTYPE_DOUBLE f_filler;
-#endif
-	double d_quadDoubleValue[4];
-};
-
-#define __QuadDoubleInstPtr(obj)      ((struct __quadDoubleStruct *)(__objPtr(obj)))
-
-#ifndef __isQuadDouble
-# define __isQuadDouble(o) \
-	(__Class(o) == @global(QuadDouble))
-#endif
-
-#ifndef __qIsQuadDouble
-# define __qIsQuadDouble(o) \
-	(__qClass(o) == @global(QuadDouble))
-#endif
-
-#define __qNew_qdReal(newQD, d0,d1,d2,d3) { \
-    __qNew(newQD, sizeof(struct __quadDoubleStruct)); \
-    __stx_setClass(newQD, QDouble);                \
-    __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[0] = d0;   \
-    __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[1] = d1;   \
-    __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[2] = d2;   \
-    __QuadDoubleInstPtr(newQD)->d_quadDoubleValue[3] = d3;   \
-}
-
-// sigh: not all compilers (borland) support inline functions;
-// therefore we have to use macros...
-// sigh2: c-macros are unhygienic - to avoid catching/hiding variable bindings,
-// use different names in each macros (i.e. a_xxx)
-
-// qd_real(c0, c1, c2, c3);
-
-#define _QD_SPLITTER 134217729.0               // = 2^27 + 1
-#define _QD_SPLIT_THRESH 6.69692879491417e+299 // = 2^996
-
-#define m_quick_two_sum(rslt_1, a_1, b_1, err_1)\
-{\
-    double s_1 = (a_1) + (b_1);\
-    (err_1) = (b_1) - (s_1 - (a_1));\
-    (rslt_1) = s_1; \
-}
-
-#define m_quick_two_diff(rslt_2, a_2, b_2, err_2)\
-{\
-    double s_2 = (a_2) - (b_2);\
-    (err_2) = ((a_2) - s_2) - (b_2);\
-    (rslt_2) = s_2;\
-}
-
-#define m_two_sum(rslt_3, a_3, b_3, err_3)\
-{\
-    double s_3 = a_3 + b_3;\
-    double bb_3 = s_3 - a_3;\
-    err_3 = (a_3 - (s_3 - bb_3)) + (b_3 - bb_3);\
-    rslt_3 = s_3;\
-}
-
-/* Computes fl(a-b) and err(a-b).  */
-#define m_two_diff(rslt_4, a_4, b_4, err_4)\
-{\
-    double s_4 = a_4 - b_4;\
-    double bb_4 = s_4 - a_4;\
-    err_4 = (a_4 - (s_4 - bb_4)) - (b_4 + bb_4);\
-    rslt_4 = s_4;\
-}
-
-#define m_three_sum(a_5, b_5, c_5)\
-{ \
-    double t1_5, t2_5, t3_5; \
-    m_two_sum(t1_5, a_5, b_5, t2_5); \
-    m_two_sum(a_5, c_5, t1_5, t3_5); \
-    m_two_sum(b_5, t2_5, t3_5, c_5); \
-}
-
-#define m_three_sum2(a_6, b_6, c_6)\
-{\
-    double t1_6, t2_6, t3_6;\
-    m_two_sum(t1_6, a_6, b_6, t2_6);\
-    m_two_sum(a_6, c_6, t1_6, t3_6);\
-    b_6 = t2_6 + t3_6;\
-}
-
-#ifndef QD_FMS
-
-/* Computes high word and lo word of a */
-#define m_split(a_7, hi_7, lo_7)\
-{\
-    double temp_7;\
-    if (a_7 > _QD_SPLIT_THRESH || a_7 < -_QD_SPLIT_THRESH) {\
-	a_7 *= 3.7252902984619140625e-09;\
-	temp_7 = _QD_SPLITTER * a_7;\
-	hi_7 = temp_7 - (temp_7 - a_7);\
-	lo_7 = a_7 - hi_7;\
-	hi_7 *= 268435456.0;\
-	lo_7 *= 268435456.0;\
-    } else {\
-	temp_7 = _QD_SPLITTER * a_7;\
-	hi_7 = temp_7 - (temp_7 - a_7);\
-	lo_7 = a_7 - hi_7;\
-    }\
-}
-
-#endif
-
-
-#ifdef QD_FMS
-
-/* Computes fl(a*b) and err(a*b). */
-
-#define m_two_prod(rslt_8, a_8, b_8, err_8)\
-{\
-    double p_8 = a_8 * b_8;\
-    err_8 = QD_FMS(a_8, b_8, p_8);\
-    rslt_8 = p_8; \
-}
-
-#else
-
-#define m_two_prod(rslt_8, a_8, b_8, err_8)\
-{\
-    double a_hi_8, a_lo_8, b_hi_8, b_lo_8;\
-    double p_8 = a_8 * b_8;\
-    m_split(a_8, a_hi_8, a_lo_8);\
-    m_split(b_8, b_hi_8, b_lo_8);\
-    err_8 = ((a_hi_8 * b_hi_8 - p_8) + a_hi_8 * b_lo_8 + a_lo_8 * b_hi_8) + a_lo_8 * b_lo_8;\
-    rslt_8 = p_8; \
-}
-
-#endif
-
-#ifdef QD_FMS
-
-#define m_two_sqr(rslt_9, a_9, err_9)\
-{\
-    double p_9 = a_9 * a_9;\
-    err_9 = QD_FMS(a_9, a_9, p_9);\
-    rslt_9 = p_9;\
-}
-
-#else
-
-#define m_two_sqr(rslt_9, a_9, err_9)\
-{\
-    double hi_9, lo_9;\
-    double q_9 = a_9 * a_9;\
-    m_split(a_9, hi_9, lo_9);\
-    err_9 = ((hi_9 * hi_9 - q_9) + 2.0 * hi_9 * lo_9) + lo_9 * lo_9;\
-    rslt_9 = q_9;\
-}
-
-#endif
-
-#define m_renorm4(c0_10, c1_10, c2_10, c3_10)\
-{\
-    double s0_10, s1_10, s2_10 = 0.0, s3_10 = 0.0;\
-\
-    if (! isinf(c0_10)) { \
-\
-	m_quick_two_sum(s0_10, c2_10, c3_10, c3_10);\
-	m_quick_two_sum(s0_10, c1_10, s0_10, c2_10);\
-	m_quick_two_sum(c0_10, c0_10, s0_10, c1_10);\
-\
-	s0_10 = c0_10;\
-	s1_10 = c1_10;\
-	if (s1_10 != 0.0) {\
-	     m_quick_two_sum(s1_10, s1_10, c2_10, s2_10);\
-	    if (s2_10 != 0.0) {\
-		 m_quick_two_sum(s2_10, s2_10, c3_10, s3_10);\
-	    } else {\
-		 m_quick_two_sum(s1_10, s1_10, c3_10, s2_10);\
-	    }\
-	} else {\
-	    m_quick_two_sum(s0_10, s0_10, c2_10, s1_10);\
-	    if (s1_10 != 0.0) {\
-		 m_quick_two_sum(s1_10, s1_10, c3_10, s2_10);\
-	    } else {\
-		 m_quick_two_sum(s0_10, s0_10, c3_10, s1_10);\
-	    }\
-	}\
-\
-	c0_10 = s0_10;\
-	c1_10 = s1_10;\
-	c2_10 = s2_10;\
-	c3_10 = s3_10;\
-    }\
-}
-
-#define m_renorm5(c0_11, c1_11, c2_11, c3_11, c4_11)\
-{\
-    double s0_11, s1_11, s2_11 = 0.0, s3_11 = 0.0; \
-\
-    if (! isinf(c0_11)) { \
-	m_quick_two_sum(s0_11, c3_11, c4_11, c4_11); \
-	m_quick_two_sum(s0_11, c2_11, s0_11, c3_11); \
-	m_quick_two_sum(s0_11, c1_11, s0_11, c2_11); \
-	m_quick_two_sum(c0_11, c0_11, s0_11, c1_11); \
-\
-	s0_11 = c0_11; \
-	s1_11 = c1_11; \
-\
-	m_quick_two_sum(s0_11, c0_11, c1_11, s1_11); \
-	if (s1_11 != 0.0) { \
-	    m_quick_two_sum(s1_11, s1_11, c2_11, s2_11); \
-	    if (s2_11 != 0.0) { \
-		m_quick_two_sum(s2_11 ,s2_11, c3_11, s3_11); \
-		if (s3_11 != 0.0) {\
-		    s3_11 += c4_11; \
-		} else {\
-		    s2_11 += c4_11;\
-		}\
-	    } else { \
-		m_quick_two_sum(s1_11, s1_11, c3_11, s2_11); \
-		if (s2_11 != 0.0) {\
-		    m_quick_two_sum(s2_11, s2_11, c4_11, s3_11); \
-		} else { \
-		    m_quick_two_sum(s1_11, s1_11, c4_11, s2_11); \
-		} \
-	    } \
-	} else { \
-	    m_quick_two_sum(s0_11,s0_11, c2_11, s1_11); \
-	    if (s1_11 != 0.0) { \
-		m_quick_two_sum(s1_11,s1_11, c3_11, s2_11); \
-		if (s2_11 != 0.0) {\
-		    m_quick_two_sum(s2_11,s2_11, c4_11, s3_11); \
-		} else { \
-		    m_quick_two_sum(s1_11 ,s1_11, c4_11, s2_11); \
-		} \
-	    } else { \
-		m_quick_two_sum(s0_11,s0_11, c3_11, s1_11); \
-		if (s1_11 != 0.0) { \
-		    m_quick_two_sum(s1_11,s1_11, c4_11, s2_11); \
-		} else { \
-		    m_quick_two_sum(s0_11,s0_11, c4_11, s1_11); \
-		} \
-	    } \
-	} \
- \
-	c0_11 = s0_11; \
-	c1_11 = s1_11; \
-	c2_11 = s2_11; \
-	c3_11 = s3_11; \
-    } \
-}
-
-%}
+ry:'Magnitude-Numbers'
 ! !
 
 !QDouble primitiveFunctions!
-%{
-
-#if 0
-
-/*********** Basic Functions ************/
-/* Computes fl(a+b) and err(a+b).  Assumes |a| >= |b|. */
-inline double
-quick_two_sum(double a, double b, double *errPtr) {
-  double s = a + b;
-  *errPtr = b - (s - a);
-  return s;
-}
-
-/* Computes fl(a-b) and err(a-b).  Assumes |a| >= |b| */
-inline double
-quick_two_diff(double a, double b, double *errPtr) {
-  double s = a - b;
-  *errPtr = (a - s) - b;
-  return s;
-}
-
-/* Computes fl(a+b) and err(a+b).  */
-inline double
-two_sum(double a, double b, double *errPtr) {
-  double s = a + b;
-  double bb = s - a;
-  *errPtr = (a - (s - bb)) + (b - bb);
-  return s;
-}
-
-/* Computes fl(a-b) and err(a-b).  */
-inline double
-two_diff(double a, double b, double *errPtr) {
-  double s = a - b;
-  double bb = s - a;
-  *errPtr = (a - (s - bb)) - (b + bb);
-  return s;
-}
-
-#ifndef QD_FMS
-/* Computes high word and lo word of a */
-inline void
-split(double a, double *hiPtr, double *loPtr) {
-  double temp;
-  if (a > _QD_SPLIT_THRESH || a < -_QD_SPLIT_THRESH) {
-    a *= 3.7252902984619140625e-09;  // 2^-28
-    temp = _QD_SPLITTER * a;
-    *hiPtr = temp - (temp - a);
-    *loPtr = a - *hiPtr;
-    *hiPtr *= 268435456.0;          // 2^28
-    *loPtr *= 268435456.0;          // 2^28
-  } else {
-    temp = _QD_SPLITTER * a;
-    *hiPtr = temp - (temp - a);
-    *loPtr = a - *hiPtr;
-  }
-}
-#endif
-
-/* Computes fl(a*b) and err(a*b). */
-inline double
-two_prod(double a, double b, double *errPtr) {
-#ifdef QD_FMS
-  double p = a * b;
-  *errPtr = QD_FMS(a, b, p);
-  return p;
-#else
-  double a_hi, a_lo, b_hi, b_lo;
-  double p = a * b;
-  split(a, &a_hi, &a_lo);
-  split(b, &b_hi, &b_lo);
-  *errPtr = ((a_hi * b_hi - p) + a_hi * b_lo + a_lo * b_hi) + a_lo * b_lo;
-  return p;
-#endif
-}
-
-/* Computes fl(a*a) and err(a*a).  Faster than the above method. */
-inline double
-two_sqr(double a, double *errPtr) {
-#ifdef QD_FMS
-  double p = a * a;
-  *errPtr = QD_FMS(a, a, p);
-  return p;
-#else
-  double hi, lo;
-  double q = a * a;
-  split(a, &hi, &lo);
-  *errPtr = ((hi * hi - q) + 2.0 * hi * lo) + lo * lo;
-  return q;
-#endif
-}
-
-/* Computes the nearest integer to d. */
-inline double
-nint(double d) {
-  if (d == floor(d))
-    return d;
-  return floor(d + 0.5);
-}
-
-/* Computes the truncated integer. */
-inline double
-aint(double d) {
-  return (d >= 0.0) ? floor(d) : ceil(d);
-}
-
-inline void
-renorm4(double *c0Ptr, double *c1Ptr, double *c2Ptr, double *c3Ptr) {
-  double s0, s1, s2 = 0.0, s3 = 0.0;
-  double c0 = *c0Ptr;
-
-  if (isinf(c0)) return;
-
-  s0 = quick_two_sum(*c2Ptr, *c3Ptr, c3Ptr);
-  s0 = quick_two_sum(*c1Ptr, s0, c2Ptr);
-  c0 = quick_two_sum(c0, s0, c1Ptr);
-
-  s0 = c0;
-  s1 = *c1Ptr;
-  if (s1 != 0.0) {
-    s1 = quick_two_sum(s1, *c2Ptr, &s2);
-    if (s2 != 0.0)
-      s2 = quick_two_sum(s2, *c3Ptr, &s3);
-    else
-      s1 = quick_two_sum(s1, *c3Ptr, &s2);
+= quick_two_sum(s1, *c3Ptr, &s2);
   } else {
     s0 = quick_two_sum(s0, *c2Ptr, &s1);
-    if (s1 != 0.0)
-      s1 = quick_two_sum(s1, *c3Ptr, &s2);
-    else
-      s0 = quick_two_sum(s0, *c3Ptr, &s1);
-  }
-
-  *c0Ptr = s0;
-  *c1Ptr = s1;
-  *c2Ptr = s2;
-  *c3Ptr = s3;
-}
-
-inline void
-renorm5(double *c0Ptr, double *c1Ptr, double *c2Ptr, double *c3Ptr, double *c4Ptr) {
-  double s0, s1, s2 = 0.0, s3 = 0.0;
-
-  if (isinf(*c0Ptr)) return;
-
-  s0 = quick_two_sum(*c3Ptr, *c4Ptr, c4Ptr);
-  s0 = quick_two_sum(*c2Ptr, s0, c3Ptr);
-  s0 = quick_two_sum(*c1Ptr, s0, c2Ptr);
-  *c0Ptr = quick_two_sum(*c0Ptr, s0, c1Ptr);
-
-  s0 = *c0Ptr;
-  s1 = *c1Ptr;
-
-  s0 = quick_two_sum(*c0Ptr, *c1Ptr, &s1);
-  if (s1 != 0.0) {
-    s1 = quick_two_sum(s1, *c2Ptr, &s2);
-    if (s2 != 0.0) {
-      s2 =quick_two_sum(s2, *c3Ptr, &s3);
-      if (s3 != 0.0)
-	s3 += *c4Ptr;
-      else
-	s2 += *c4Ptr;
-    } else {
-      s1 = quick_two_sum(s1, *c3Ptr, &s2);
-      if (s2 != 0.0)
-	s2 = quick_two_sum(s2, *c4Ptr, &s3);
-      else
-	s1 = quick_two_sum(s1, *c4Ptr, &s2);
-    }
-  } else {
-    s0 = quick_two_sum(s0, *c2Ptr, &s1);
-    if (s1 != 0.0) {
-      s1 = quick_two_sum(s1, *c3Ptr, &s2);
-      if (s2 != 0.0)
-	s2 = quick_two_sum(s2, *c4Ptr, &s3);
-      else
-	s1 = quick_two_sum(s1, *c4Ptr, &s2);
-    } else {
-      s0 = quick_two_sum(s0, *c3Ptr, &s1);
-      if (s1 != 0.0)
-	s1 = quick_two_sum(s1, *c4Ptr, &s2);
-      else
-	s0 = quick_two_sum(s0, *c4Ptr, &s1);
-    }
-  }
-
-  *c0Ptr = s0;
-  *c1Ptr = s1;
-  *c2Ptr = s2;
-  *c3Ptr = s3;
-}
-
-inline void
-three_sum(double *aPtr, double *bPtr, double *cPtr) {
-  double t1, t2, t3;
-  t1 = two_sum(*aPtr, *bPtr, &t2);
-  *aPtr  = two_sum(*cPtr, t1, &t3);
-  *bPtr  = two_sum(t2, t3, cPtr);
-}
-
-inline void three_sum2(double *aPtr, double *bPtr, double *cPtr) {
-  double t1, t2, t3;
-  t1 = two_sum(*aPtr, *bPtr, &t2);
-  *aPtr  = two_sum(*cPtr, t1, &t3);
-  *bPtr = t2 + t3;
-}
-#endif
-
-#if 0
-/* These are provided to give consistent
-   interface for double with double-double and quad-double. */
-inline void
-sincosh(double t, double &sinh_t, double &cosh_t) {
-  sinh_t = sinh(t);
-  cosh_t = cosh(t);
-}
-
-inline double
-sqr(double t) {
-  return t * t;
-}
-
-inline double
-to_double(double a) {
-    return a;
-}
-
-inline int
-to_int(double a)    {
-    return static_cast<int>(a);
-}
-#endif
-
-%}
-! !
+    if (s1 ! !
 
 !QDouble class methodsFor:'documentation'!
 
@@ -886,6 +308,12 @@
     "Created: / 14-06-2017 / 19:14:49 / cg"
 !
 
+infinity
+    ^ Infinity positive
+
+    "Created: / 18-06-2017 / 23:41:07 / cg"
+!
+
 ln10
     "return the constant e as quad precision double.
      (returns approx. 212 bits of precision)"
@@ -926,6 +354,12 @@
     "Created: / 12-06-2017 / 18:31:34 / cg"
 !
 
+negativeInfinity
+    ^ Infinity negative
+
+    "Created: / 18-06-2017 / 23:40:47 / cg"
+!
+
 pi
     "return the constant pi as quad precision double.
      (returns approx. 212 bits of precision)"
@@ -1816,6 +1250,16 @@
     "Modified (comment): / 13-06-2017 / 17:33:19 / cg"
 !
 
+ln
+    self isOne ifTrue:[ ^ self class zero ].
+
+    "
+     1.0 asQDouble ln
+    "
+
+    "Created: / 18-06-2017 / 23:32:54 / cg"
+!
+
 negated
     ^ self class
 	d0:(self d0) negated
@@ -2397,6 +1841,21 @@
     ^ self d0 isNaN
 
     "Created: / 15-06-2017 / 01:57:35 / cg"
+!
+
+isOne
+    ^ self d0 = 1.0
+    and:[ self d1 = 0.0
+    and:[ self d2 = 0.0
+    and:[ self d3 = 0.0 ]]]
+
+    "Created: / 18-06-2017 / 23:29:07 / cg"
+!
+
+isZero
+    ^ self d0 = 0.0
+
+    "Created: / 18-06-2017 / 23:29:25 / cg"
 ! !
 
 !QDouble class methodsFor:'documentation'!
@@ -2408,3 +1867,4 @@
 version_CVS
     ^ '$Header$'
 ! !
+