/*
Copyright (c) 2000-2010, Dirk Krause
All rights reserved.

Redistribution and use in source and binary forms,
with or without modification, are permitted provided
that the following conditions are met:

* Redistributions of source code must retain the above
  copyright notice, this list of conditions and the
  following disclaimer.
* Redistributions in binary form must reproduce the above 
  opyright notice, this list of conditions and the following
  disclaimer in the documentation and/or other materials
  provided with the distribution.
* Neither the name of the Dirk Krause nor the names of
  contributors may be used to endorse or promote
  products derived from this software without specific
  prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED.
IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGE.
*/



/**	@file	dkma.c	Mathematical operations.
The operations do checks for range overflows and zero divisions.
*/


/**	The code for negli et al was modified, we use the new code.
*/
#define OLD_VERSION_NEGLI 0

#include "dk.h"
#include "dkstr.h"
#include "dkma.h"
#include "dkerror.h"

#if DK_HAVE_STRING_H
#include <string.h>
#endif
#if DK_HAVE_STRINGS_H
#include <strings.h>
#endif
#if DK_HAVE_MATH_H
#include <math.h>
#endif
#if DK_HAVE_STDLIB_H
#include <stdlib.h>
#endif




#line 69 "dkma.ctr"




#if DK_HAVE_LONG_LONG_INT

/**	MSB of unsigned long long.
*/
#define MSBULLI	0x8000000000000000ULL
/**	Maximum unsigned long long number, expressed as double.
*/
#define MSBULLD 9223372036854775808.0
/**	Maximum unsigned long long number.
*/
#define BITSULL 64
/**	Maximum unsigned long long number.
*/
#define MAXULLI 0xFFFFFFFFFFFFFFFFULL

#else
/* if DK_HAVE_LONG_LONG_INT */

/**	MSB of unsigned long long.
*/
#define MSBULLI 0x80000000UL
/**	Maximum unsigned long long number, expressed as double.
*/
#define MSBULLD 2147483648.0
/**	Number of bits for unsigned long long.
*/
#define BITSULL 32
/**	Maximum unsigned long long number.
*/
#define MAXULLI 0xFFFFFFFFUL

#endif
/* if DK_HAVE_LONG_LONG_INT */



#if SIZEOF_LONG == 8
/**	MSB for unsigned long in hex notation.
*/
#define MSBULI 0x8000000000000000UL
/**	MSB for unsigned long in double notation.
*/
#define MSBULD 9223372036854775808.0
/** 0xFFFFFFFFFFFFFFFFUL in double notation.
*/
#define MAXULD 18446744073709551615.0
/**	MSB for long in hex notation.
*/
#define MSBLI  0x4000000000000000L
/**	MSB for long in double notation.
*/
#define MSBLD  4611686018427387904.0
/**	Maximum value for unsigned long, expressed as double.
*/
#define MAXLD  9223372036854775807.0
/**	Maximum value for signed int variable, expressed as long.
*/
#define MAXLI  0x7FFFFFFFFFFFFFFFL
/**	Number of bits in unsigned long.
*/
#define BITSUL 64
/**	Number of bits in long.
*/
#define BITSL  63
#if OLD_VERSION_NEGLI
/**     Negativest int value, expressed as long.
*/
#define NEGLI  (-9223372036854775808L)
#else
/**	Negativest int value, expressed as long.
*/
#define NEGLI  (-9223372036854775807L)
#endif
/** Negativest long value, expressed as double.
*/
#define NEGLD  (-9223372036854775807.0)
#else
/**	leftmost bit of unsigned long.
*/
#define MSBULI 0x80000000UL
/**	Leftmost bit of unsigned long, expressed as double.
*/
#define MSBULD 2147483648.0
/**	Maximum unsigned long value, expressed as double.
*/
#define MAXULD 4294967295.0
/**	Leftmost bit of signed int, expressed as long.
*/
#define MSBLI  0x40000000L
/**	Leftmost bit of signed long value, expressed as double.
*/
#define MSBLD  1073741824.0
/**	Maximum long value, expressed as double.
*/
#define MAXLD  2147483647.0
/**	Maximum long value.
*/
#define MAXLI  0x7FFFFFFFL
/**	Number of bits in an unsigned long variable.
*/
#define BITSUL 32
/**	Number of bits in a long variable.
*/
#define BITSL  31
#if OLD_VERSION_NEGLI
/**	Negativest int value, expressed as long.
*/
#define NEGLI  (-2147483648L)
#else
/* define NEGLI  0x80000000L */
/**	Negativest int value, expressed as long.
*/
#define NEGLI  (-2147483647L)
#endif
/** Negativest long value, expressed as double.
*/
#define NEGLD  (-2147483647.0)
#endif



/**	This error code variable is set if no pointer is provided
	to another error code variable.
*/
static int error_code = DK_ERR_NONE;



int dkma_get_error DK_P1(int, del)
{
  int back;
  back = error_code;
  if(del) { error_code = 0; }
  
  return back;
}



/**	Addition of 2 unsigned long values.
	@param	u1	First operand.
	@param	u2	Second operand.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
unsigned long
st_add_ulong_ok DK_P3(unsigned long, u1, unsigned long, u2, int *, ok)
{
  unsigned long back;
  
  if((DK_MAX_ULONG - u1) < u2) {
    if(ok) { *ok = DK_ERR_MATH_OOR; }
  }
  back = u1 + u2;
  
  return back;
}



/**	Addition of 2 unsigned long long values.
	@param	o1	First operand.
	@param	o2	Second operand.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
dk_long_long_unsigned_t
st_add_ull_ok DK_P3(dk_long_long_unsigned_t,o1,dk_long_long_unsigned_t,o2,int *,ok)
{
  dk_long_long_unsigned_t back;
  if((DK_MAX_LONG_LONG_UNSIGNED - o1) < o2) {
    if(ok) { *ok = DK_ERR_MATH_OOR; }
  }
  back = o1 + o2;
  return back;
}



/**	Substraction of 2 unsigned long values.
	@param	u1	First operand.
	@param	u2	Second operand.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
unsigned long
st_sub_ulong_ok DK_P3(unsigned long, u1, unsigned long, u2, int *, ok)
{
  unsigned long back;
  
  if(u2 > u1) {
    if(ok) { *ok = DK_ERR_MATH_OOR; }
  }
  back = u1 - u2;
  
  return back;
}



/**	Substraction of 2 unsigned long long values.
	@param	o1	First operand.
	@param	o2	Second operand.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
dk_long_long_unsigned_t
st_sub_ull_ok DK_P3(\
dk_long_long_unsigned_t,o1,dk_long_long_unsigned_t,o2,int *,ok)
{
  dk_long_long_unsigned_t back;
  if(o2 > o1) {
    if(ok) { *ok = DK_ERR_MATH_OOR; }
  }
  back = o1 - o2;
  return back;
}



/**	Multiplication of 2 unsigned long values.
	@param	u1	First operand.
	@param	u2	Second operand.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
unsigned long
st_mul_ulong_ok DK_P3(unsigned long, u1, unsigned long, u2, int *, ok)
{
  unsigned long back;
  
  if(u1 > (DK_MAX_ULONG / u2)) {
    if(ok) { *ok = DK_ERR_MATH_OOR; }
  }
  back = u1 * u2;
  
  return back;
}



/**	Multiplication of 2 unsigned long long values.
	@param	o1	First operand.
	@param	o2	Second operand.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
dk_long_long_unsigned_t
st_mul_ull_ok DK_P3(\
dk_long_long_unsigned_t,o1,dk_long_long_unsigned_t,o2,int *,ok)
{
  dk_long_long_unsigned_t back;
  if(o1 > (DK_MAX_LONG_LONG_UNSIGNED / o2)) {
    if(ok) { *ok = DK_ERR_MATH_OOR; }
  }
  back = o1 * o2;
  return back;
}




/**	Division of 2 unsigned long values.
	@param	u1	First operand.
	@param	u2	Second operand.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
unsigned long
st_div_ulong_ok DK_P3(unsigned long, u1, unsigned long, u2, int *, ok)
{
  unsigned long back;
  
  if(u2) {
    back = u1 / u2;
  } else {
    back = 0UL;
    if(ok) { *ok = DK_ERR_DIV_ZERO; }
  }
  
  return back;
}



/**	Division of 2 unsigned long long values.
	@param	o1	First operand.
	@param	o2	Second operand.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
dk_long_long_unsigned_t
st_div_ull_ok DK_P3(\
dk_long_long_unsigned_t,o1,dk_long_long_unsigned_t,o2,int *,ok)
{
  dk_long_long_unsigned_t back;
  if(o2) {
    back = o1 / o2;
  } else {
    back = DK_ZERO_LONG_LONG_UNSIGNED;
    if(ok) { *ok = DK_ERR_DIV_ZERO; }
  }
  return back;
}



unsigned long
dkma_add_ulong_ok DK_P3(unsigned long, u1, unsigned long, u2, int *, ok)
{
  unsigned long back;
  back = st_add_ulong_ok(u1,u2,(ok ? ok : &error_code));
  return back;
}



dk_long_long_unsigned_t
dkma_add_ull_ok DK_P3(dk_long_long_unsigned_t,o1,dk_long_long_unsigned_t,o2,int *,ok)
{
  dk_long_long_unsigned_t back;
  back = st_add_ull_ok(o1,o2,(ok ? ok : &error_code));
  return back;
}



unsigned long
dkma_sub_ulong_ok DK_P3(unsigned long, u1, unsigned long, u2, int *, ok)
{
  unsigned long back;
  back = st_sub_ulong_ok(u1,u2,(ok ? ok : &error_code));
  return back;
}



dk_long_long_unsigned_t
dkma_sub_ull_ok DK_P3(dk_long_long_unsigned_t,o1,dk_long_long_unsigned_t,o2,int *,ok)
{
  dk_long_long_unsigned_t back;
  back = st_sub_ull_ok(o1,o2,(ok ? ok : &error_code));
  return back;
}



unsigned long
dkma_mul_ulong_ok DK_P3(unsigned long, u1, unsigned long, u2, int *, ok)
{
  unsigned long back;
  back = st_mul_ulong_ok(u1,u2,(ok ? ok : &error_code));
  return back;
}



dk_long_long_unsigned_t
dkma_mul_ull_ok DK_P3(dk_long_long_unsigned_t,o1,dk_long_long_unsigned_t,o2,int *,ok)
{
  dk_long_long_unsigned_t back;
  back = st_mul_ull_ok(o1,o2,(ok ? ok : &error_code));
  return back;
}



unsigned long
dkma_div_ulong_ok DK_P3(unsigned long, u1, unsigned long, u2, int *, ok)
{
  unsigned long back;
  back = st_div_ulong_ok(u1,u2,(ok ? ok : &error_code));
  return back;
}



dk_long_long_unsigned_t
dkma_div_ull_ok DK_P3(dk_long_long_unsigned_t,o1,dk_long_long_unsigned_t,o2,int *,ok)
{
  dk_long_long_unsigned_t back;
  back = st_div_ull_ok(o1,o2,(ok ? ok : &error_code));
  return back;
}



unsigned long
dkma_add_ulong DK_P2(unsigned long, u1, unsigned long, u2)
{
  unsigned long back;
  back = st_add_ulong_ok(u1,u2,&error_code);
  return back;
}



unsigned long
dkma_sub_ulong DK_P2(unsigned long, u1, unsigned long, u2)
{
  unsigned long back;
  back = st_sub_ulong_ok(u1,u2,&error_code);
  return back;
}



unsigned long
dkma_mul_ulong DK_P2(unsigned long, u1, unsigned long, u2)
{
  unsigned long back;
  back = st_mul_ulong_ok(u1,u2,&error_code);
  return back;
}



unsigned long
dkma_div_ulong DK_P2(unsigned long, u1, unsigned long, u2)
{
  unsigned long back;
  back = st_div_ulong_ok(u1,u2,&error_code);
  return back;
}



dk_long_long_unsigned_t
dkma_add_ull DK_P2(dk_long_long_unsigned_t,o1,dk_long_long_unsigned_t,o2)
{
  dk_long_long_unsigned_t back;
  back = st_add_ull_ok(o1,o2,&error_code);
  return back;
}



dk_long_long_unsigned_t
dkma_sub_ull DK_P2(dk_long_long_unsigned_t,o1,dk_long_long_unsigned_t,o2)
{
  dk_long_long_unsigned_t back;
  back = st_sub_ull_ok(o1,o2,&error_code);
  return back;
}



dk_long_long_unsigned_t
dkma_mul_ull DK_P2(dk_long_long_unsigned_t,o1,dk_long_long_unsigned_t,o2)
{
  dk_long_long_unsigned_t back;
  back = st_mul_ull_ok(o1,o2,&error_code);
  return back;
}



dk_long_long_unsigned_t
dkma_div_ull DK_P2(dk_long_long_unsigned_t,o1,dk_long_long_unsigned_t,o2)
{
  dk_long_long_unsigned_t back;
  back = st_div_ull_ok(o1,o2,&error_code);
  return back;
}



/**	Addition of 2 long numbers.
	@param	l1	First operand.
	@param	l2	Second operand.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
long
st_add_long_ok DK_P3(long, l1, long, l2, int *, ok)
{
  long back;
  
  back = l1 + l2;
  if(l1 >= 0L) {
    if(l2 >= 0L) {
      if((DK_MAX_LONG - l1) < l2) {
        if(ok) { *ok = DK_ERR_MATH_OOR; }
      }
    }
  } else {
    if(l2 < 0L) {
      if((DK_MAX_LONG + l1) < (0L - l2)) {
        if(ok) { *ok = DK_ERR_MATH_OOR; }
      }
    }
  }
  
  return back;
}



/**	Substraction of 2 long numbers.
	@param	l1	First operand.
	@param	l2	Second operand.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
long
st_sub_long_ok DK_P3(long, l1, long, l2, int *, ok)
{
  long back;
  
  back = l1 - l2;
  if(l1 >= 0) {
    if(l2 < 0L) {
      if((DK_MAX_LONG + l2) < l1) {
        if(ok) { *ok = DK_ERR_MATH_OOR; }
      }
    }
  } else {
    if(l2 > 0) {
      if((DK_MAX_LONG + l1) < l2) {
        if(ok) { *ok = DK_ERR_MATH_OOR; }
      }
    }
  }
  
  return back;
}



/**	Multiplication of 2 long numbers.
	@param	l1	First operand.
	@param	l2	Second operand.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
long
st_mul_long_ok DK_P3(long, l1, long, l2, int *, ok)
{
  long back, a1, a2;
  
  back = l1 * l2;
  a1 = labs(l1); if(a1 > 1L) {
    a2 = labs(l2); if(a2 > 1L) {
      if((DK_MAX_LONG / a1) < a2) {
        if(ok) { *ok = DK_ERR_MATH_OOR; }
      }
    }
  }
  
  return back;
}



/**	Division of 2 long numbers.
	@param	l1	First operand.
	@param	l2	Second operand.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
long
st_div_long_ok DK_P3(long, l1, long, l2, int *, ok)
{
  long back;
  
  if(l2) {
    back = l1 / l2;
  } else {
    back = 0L;
    if(ok) { *ok = DK_ERR_DIV_ZERO; }
  } 
  return back;
}



long dkma_add_long_ok DK_P3(long, l1, long, l2, int *, ok)
{
  long back;
  back = st_add_long_ok(l1,l2,(ok ? ok : &error_code));
  return back;
}



long dkma_sub_long_ok DK_P3(long, l1, long, l2, int *, ok)
{
  long back;
  back = st_sub_long_ok(l1,l2,(ok ? ok : &error_code));
  return back;
}



long dkma_mul_long_ok DK_P3(long, l1, long, l2, int *, ok)
{
  long back;
  back = st_mul_long_ok(l1,l2,(ok ? ok : &error_code));
  return back;
}



long dkma_div_long_ok DK_P3(long, l1, long, l2, int *, ok)
{
  long back;
  back = st_div_long_ok(l1,l2,(ok ? ok : &error_code));
  return back;
}



long dkma_add_long DK_P2(long, l1, long, l2)
{
  long back;
  back = st_add_long_ok(l1,l2,&error_code);
  return back;
}



long dkma_sub_long DK_P2(long, l1, long, l2)
{
  long back;
  back = st_sub_long_ok(l1,l2,&error_code);
  return back;
}



long dkma_mul_long DK_P2(long, l1, long, l2)
{
  long back;
  back = st_mul_long_ok(l1,l2,&error_code);
  return back;
}



long dkma_div_long DK_P2(long, l1, long, l2)
{
  long back;
  back = st_div_long_ok(l1,l2,&error_code);
  return back;
}



/**	Addition of 2 double values.
	@param	d1	First operand.
	@param	d2	Second operand.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
double
st_add_double_ok DK_P3(double, d1, double, d2, int *, ok)
{
  double back;
  
  back = d1 + d2;
  if(d1 >= 0.0) {
    if(d2 >= 0.0) {
      if((DK_MAX_DOUBLE - d1) < d2) {
        if(ok) { *ok = DK_ERR_MATH_OOR; }
      }
    }
  } else {
    if(d2 < 0.0) {
      if((DK_MAX_DOUBLE + d1) < (0.0 - d2)) {
        if(ok) { *ok = DK_ERR_MATH_OOR; }
      }
    }
  }
  
  return back;
}



/**	Subtraction of 2 double values.
	@param	d1	First operand.
	@param	d2	Second operand.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
double
st_sub_double_ok DK_P3(double, d1, double, d2, int *, ok)
{
  double back;
  
  back = d1 - d2;
  if(d1 >= 0.0) {
    if(d2 < 0.0) {
      if((DK_MAX_DOUBLE + d2) < d1) {
        if(ok) { *ok = DK_ERR_MATH_OOR; }
      }
    }
  } else {
    if(d2 >= 0.0) {
      if((DK_MAX_DOUBLE + d1) < d2) {
        if(ok) { *ok = DK_ERR_MATH_OOR; }
      }
    }
  }
  
  return back;
}



/**	Multiplication of 2 double values.
	@param	d1	Factor.
	@param	d2	Factor.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
double
st_mul_double_ok DK_P3(double, d1, double, d2, int *, ok)
{
  double back, a1, a2;
  
  back = d1 * d2;
  a1 = fabs(d1);
  if(a1 > 1.0) {
    a2 = fabs(d2);
    if(a2 > 1.0) {
      if((DK_MAX_DOUBLE / a1) < a2) {
        if(ok) { *ok = DK_ERR_MATH_OOR; }
      }
    }
  }
  
  return back;
}



/**	Division of 2 double values.
	@param	d1	Dividend.
	@param	d2	Divisor.
	@param	ok	Pointer to error code variable.
	@return	Result.
*/
static
double
st_div_double_ok DK_P3(double, d1, double, d2, int *, ok)
{
  double back, a1, a2;
  
  a2 = fabs(d2);
  if(a2 < 1.0) {
    a1 = fabs(d1);
    if((DK_MAX_DOUBLE * a2) > a1) {
      back = d1 / d2;
    } else {
      back = 0.0;
      if(ok) { *ok = DK_ERR_MATH_OOR; }
    }
  } else {
    back = d1 / d2;
  }
  
  return back;
}



double
dkma_add_double_ok DK_P3(double, d1, double, d2, int *, ok)
{
  double back;
  back = st_add_double_ok(d1,d2,(ok ? ok : &error_code));
  return back;
}



double
dkma_sub_double_ok DK_P3(double, d1, double, d2, int *, ok)
{
  double back;
  back = st_sub_double_ok(d1,d2,(ok ? ok : &error_code));
  return back;
}



double
dkma_mul_double_ok DK_P3(double, d1, double, d2, int *, ok)
{
  double back;
  back = st_mul_double_ok(d1,d2,(ok ? ok : &error_code));
  return back;
}



double
dkma_div_double_ok DK_P3(double, d1, double, d2, int *, ok)
{
  double back;
  back = st_div_double_ok(d1,d2,(ok ? ok : &error_code));
  return back;
}



double dkma_add_double DK_P2(double, d1, double, d2)
{
  double back;
  back = st_add_double_ok(d1,d2,&error_code);
  
  return back;
}



double dkma_sub_double DK_P2(double, d1, double, d2)
{
  double back;
  back = st_sub_double_ok(d1,d2,&error_code);
  
  return back;
}



double dkma_mul_double DK_P2(double, d1, double, d2)
{
  double back;
  back = st_mul_double_ok(d1,d2,&error_code);
  
  return back;
}



double dkma_div_double DK_P2(double, d1, double, d2)
{
  double back;
  back = st_div_double_ok(d1,d2,&error_code);
  
  return back;
}



double dkma_rint DK_P1(double, d)
{
  double back = 0.0;
#if DK_HAVE_RINT
  back = rint(d);
#else
  if(d >= 0.0) {
    if(d < (DK_MAX_DOUBLE - 0.5)) {
      back = floor(d + 0.5);
    } else {
      back = DK_MAX_DOUBLE;
    }
  } else {
    back = 0.0 - dkma_rint(0.0 - d);
  }
#endif
  
  return back;
}



/**	Convert double value to unsigned long.
	@param	d	Value to convert.
	@param	ok	Pointer to error code variable.
	@return	Conversion result.
*/
static unsigned long
st_double_to_ul_ok DK_P2(double, d, int *, ok)
{
  unsigned long back = 0UL;
  double xval, testd; unsigned long testval; int trials;
  
  if(d >= 0.0) {
    xval = dkma_rint(d); 
    if(xval > MAXULD) {
      *ok = DK_ERR_MATH_OOR;
    } else {
      testd = MSBULD; testval = MSBULI; trials = BITSUL;
      while(trials--) {
	if(xval >= testd) {
	  back |= testval;
	  xval -= testd;
	}
	testd = testd / 2.0; testval = testval / 2UL;
      }
    }
  } else {
    *ok = DK_ERR_MATH_OOR;
  } 
  return back;
}



/**	Convert double value to long
	@param	d	Value to convert.
	@param	ok	Pointer to error code variable.
	@return	Conversion result.
*/
static long
st_double_to_l_ok  DK_P2(double, d, int *, ok)
{
  long back = 0L; long x;
  long testval;
  double xval, testd;
  int trials;
  
  if(d < 0.0) {
    x = st_double_to_l_ok((0.0 - d), ok);
    if(x <= DK_MAX_LONG) {
      back = 0L - x;
    } else {
      *ok = DK_ERR_MATH_OOR;
    }
  } else {
    xval = dkma_rint(d);
    if(xval > MAXLD) {
      *ok = DK_ERR_MATH_OOR;
    } else {
      testval = MSBLI; testd = MSBLD; trials = BITSL;
      while(trials--) {
        if(xval >= testd) {
	  xval -= testd;
	  back |= testval;
        }
        testval = testval / 2L; testd = testd / 2.0;
      }
    }
  } 
  return back;
}



/**	Convert unsigned long value to double.
	@param	l	Value to convert.
	@param	ok	Pointer to error code variable.
	@return	Conversion result.
*/
static double
st_ul_to_double_ok DK_P2(unsigned long, l, int *, ok)
{
  double back = 0.0;
  double testd; unsigned long testval; int trials;
  
  testd = MSBULD; testval = MSBULI; trials = BITSUL;
  while(trials--) {
    if(l & testval) { back += testd; }
    testd = testd / 2.0; testval = testval / 2UL;
  } 
  return back;
}



/**	Convert unsigned long long long value to double.
	@param	l	Value to convert.
	@param	ok	Pointer to error code variable.
	@return	Conversion result.
*/
static double
st_ull_to_double_ok DK_P2(dk_long_long_unsigned_t,l,int *,ok)
{
  double back = 0.0;
  double testd;
  dk_long_long_unsigned_t testval; int trials;
#if DK_HAVE_LONG_LONG_INT
  
#else
  
#endif
  testd = MSBULLD;
  testval = MSBULLI;
  trials = BITSULL; 
  while(trials--) {
#if DK_HAVE_LONG_LONG_INT
    
#else
    
#endif
    
    if(l & testval) {
      back += testd;
      
    }
#if DK_HAVE_LONG_LONG_INT
    testd = testd / 2.0; testval = testval / 2ULL;
#else
    testd = testd / 2.0; testval = testval / 2UL;
#endif
  } 
  return back;
}



/**	Convert long value to double.
	@param	l	Value to convert.
	@param	ok	Pointer to error code variable.
	@return	Conversion result.
*/
static double
st_l_to_double_ok  DK_P2(long, l, int *, ok)
{
  double back = 0.0;
  double testd;
  long testval;
  int trials;
  
  if(l < 0) {
    if(l == (long)NEGLI) {
      back = NEGLD; 
    } else {
      back = 0.0 - st_l_to_double_ok((0L -l), ok);
    }
  } else {
    testd = MSBLD; testval = MSBLI; trials = BITSL;
    while(trials--) {
      if(l & testval) { back += testd; }
      testd = testd / 2.0; testval = testval / 2UL;
    }
  } 
  return back;
}



unsigned long
dkma_double_to_ul_ok DK_P2(double, d, int *, ok)
{
  unsigned long back;
  back = st_double_to_ul_ok(d,(ok ? ok : &error_code));
  return back;
}



long
dkma_double_to_l_ok  DK_P2(double, d, int *, ok)
{
  long back;
  back = st_double_to_l_ok(d,(ok ? ok : &error_code));
  return back;
}



double
dkma_ul_to_double_ok DK_P2(unsigned long, l, int *, ok)
{
  double back;
  back = st_ul_to_double_ok(l,(ok ? ok : &error_code));
  return back;
}



double
dkma_ull_to_double_ok DK_P2(dk_long_long_unsigned_t,l,int *,ok)
{
  double back;
  back = st_ull_to_double_ok(l,(ok ? ok : &error_code));
  return back;
}



double
dkma_l_to_double_ok  DK_P2(long, l, int *, ok)
{
  double back;
  back = st_l_to_double_ok(l, (ok ? ok : &error_code));
  return back;
}



unsigned long dkma_double_to_ul DK_P1(double, d)
{ unsigned long back; back = st_double_to_ul_ok(d, &error_code); return back; }



long dkma_double_to_l DK_P1(double, d)
{ long back; back = st_double_to_l_ok(d, &error_code); return back; }



double dkma_ul_to_double DK_P1(unsigned long, ul)
{ double back; back = st_ul_to_double_ok(ul, &error_code); return back; }



double dkma_ull_to_double DK_P1(dk_long_long_unsigned_t, ul)
{ double back; back = st_ull_to_double_ok(ul, &error_code); return back; }



double dkma_l_to_double DK_P1(long, l)
{ double back; back = st_l_to_double_ok(l, &error_code); return back; }



double
dkma_double_restrict_digits DK_P2(double,x,unsigned,dig)
{
  double back, newval;
  unsigned cc, didcc;
  int errc;
  
  newval = back = x; cc = dig; didcc = 0; errc = 0;
  while (cc--) {
    newval = dkma_mul_double_ok(back, 10.0, &errc);
    if (errc) {
      cc = 0;
    } else {
      back = newval; didcc++;
    }
  }
  back = dkma_rint(back);
  while (didcc--) {
    back = back / 10.0;
  }
  
  return back;
}



double
dkma_double_restrict_downwards DK_P2(double,x,unsigned,dig)
{
  double back, newval;
  unsigned cc, didcc;
  int errc;
  
  newval = back = x; cc = dig; didcc = 0; errc = 0;
  while (cc--) {
    newval = dkma_mul_double_ok(back, 10.0, &errc);
    if (errc) {
      cc = 0;
    } else {
      back = newval; didcc++;
    }
  }
  back = floor(back);
  while (didcc--) {
    back = back / 10.0;
  }
  
  return back;
}



double
dkma_atan2 DK_P2(double,y, double,x)
{
  double back = 0.0;
#if DK_HAVE_ATAN2
  
  back = atan2(y, x);
#else
  int mathok = 0;
  
  back = dkma_div_double_ok(y, x, &mathok);
  if(mathok) {	
    if(y < 0.0) {	
      back = 0.0 - M_PI_2;
    } else {		
      back = M_PI_2;
    }
  } else {	
    back = atan(back);	
    if(x < 0.0) {	
      back += M_PI;
    }
    while(back < (0.0 - M_PI)) {	
      back += 2.0 * M_PI;
    }
    while(back > M_PI) {		
      back -= 2.0 * M_PI;
    }
  }
#endif
  
  return back;
}



dk_long_long_unsigned_t
dkma_ull_gcd DK_P2(dk_long_long_unsigned_t,u1, dk_long_long_unsigned_t,u2) {
  dk_long_long_unsigned_t	back	= DK_ZERO_LONG_LONG_UNSIGNED,
  				a	= DK_ZERO_LONG_LONG_UNSIGNED,
				b	= DK_ZERO_LONG_LONG_UNSIGNED,
				h	= DK_ZERO_LONG_LONG_UNSIGNED;
  while(b > DK_ZERO_LONG_LONG_UNSIGNED) {
    h = a % b; a = b; b = h;
  } back = a;
  return back;
}



dk_long_long_unsigned_t
dkma_ull_lcm_ok DK_P3(\
dk_long_long_unsigned_t,u1, dk_long_long_unsigned_t,u2, int *,ec)
{
  dk_long_long_unsigned_t back;
  back = dkma_div_ull_ok(u1, dkma_ull_gcd(u1, u2), ec);
  back = dkma_mul_ull_ok(back, u2, ec);
  return back;
}



dk_long_long_unsigned_t
dkma_ull_lcm DK_P2(dk_long_long_unsigned_t,u1, dk_long_long_unsigned_t,u2) {
  dk_long_long_unsigned_t back;
  back = dkma_div_ull(u1, dkma_ull_gcd(u1, u2));
  back = dkma_mul_ull(back, u2);
  return back;
}



unsigned long
dkma_ul_gcd DK_P2(unsigned long,ul1, unsigned long,ul2) {
  unsigned long back = 0UL, a = 0UL, b = 0UL, h = 0UL;
  a = ul1; b = ul2;
  while(b > 0UL) {
    h = a % b;
    a = b;
    b = h;
  } back = a;
  return back;
}



unsigned long
dkma_ul_lcm_ok DK_P3(unsigned long,ul1, unsigned long,ul2, int *,ec) {
  unsigned long back;
  back = dkma_div_ulong_ok(ul1, dkma_ul_gcd(ul1, ul2), ec);
  back = dkma_mul_ulong_ok(back, ul2, ec);
  return back;
}



unsigned long
dkma_ul_lcm DK_P2(unsigned long,ul1, unsigned long,ul2) {
  unsigned long back;
  back = dkma_div_ulong(ul1, dkma_ul_gcd(ul1, ul2));
  back = dkma_mul_ulong(back, ul2);
  return back;
}



long
dkma_l_gcd DK_P2(long,l1, long,l2) {
  long back = 0L, a = 0L, b = 0L, h = 0L;
  if(l1 < 0L) {
    back = dkma_l_gcd(0L - l1, l2);
  } else {
    if(l2 < 0L) {
      back = dkma_l_gcd(l1, 0L - l2);
    } else {
      a = l1; b = l2;
      while(b > 0L) {
        h = a % b;
	a = b;
	b = h;
      } back = a;
    }
  }
  return back;
}



long
dkma_l_lcm_ok DK_P3(long,l1, long,l2, int *,ec) {
  long back;
  back = dkma_div_long_ok(l1, dkma_l_gcd(l1, l2), ec);
  back = dkma_mul_long_ok(back, l2, ec);
  return back;
}



long
dkma_l_lcm DK_P2(long,l1, long,l2) {
  long back;
  back = dkma_div_long(l1, dkma_l_gcd(l1, l2));
  back = dkma_mul_long(back, l2);
  return back;
}



/**	Check character whether it is a digit and not zero.
	@param	c	Character to check.
	@return	1 on success, 0 on error.
*/
static int
is_non_zero DK_P1(char,c) {
  int back = 0;
  switch(c) {
    case '1': case '2': case '3': case '4': case '5':
    case '6': case '7': case '8': case '9': {
      back = 1;
    } break;
  }
  return back;
}



/**	Find first non-zero digit in string.
	@param	s	String to search.
	@return	Index of first digit not zero.
*/
static int
first_non_zero DK_P1(char *,s) {
  int back = -1, i = 0;
  char *p;
  p = s; i = 0;
  while((back == -1) && (*p)) {
    if(is_non_zero(*p)) { back = i; }
    p++; i++;
  }
  return back;
}



void
dkma_fputs_double_str_no_exp DK_P2(FILE *,f, char *,s) {
  char buffer[64];
  char *p1, *p2, *pe, *sptr, *dptr;
  int evalue, itest, dotpos, fnz, max, i, sl;
  size_t sz;
  p1 = dkstr_start(s, NULL);
  if(p1) {
    if(*p1 == '-') { fputc(*(p1++), f); }
    pe = strchr(p1, 'e');
    if(!pe) { pe = strchr(p1, 'E'); }
    if(pe) {
      *(pe++) = '\0';
      if(sscanf(pe, "%d", &itest) == 1) {
        evalue = itest;
	if(evalue != 0) {
	  p2 = strchr(p1, '.');
	  if(p2) { *(p2++) = '\0'; }
	  sz = 0;
	  if(p1) { sz += strlen(p1); }
	  if(p2) { sz += strlen(p2); }
	  if(sz < sizeof(buffer)) {
	    buffer[0] = '\0';
	    if(p1) { strcat(buffer, p1); }
	    if(p2) { strcat(buffer, p2); }
	    dotpos = 0;
	    if(p1) { dotpos = strlen(p1); }
	    fnz = first_non_zero(buffer);
	    if(fnz > -1) {
	      if(fnz > 0) {
	        dptr = buffer; sptr = &(buffer[fnz]);
		while(*sptr) { *(dptr++) = *(sptr++); }
		*dptr = '\0';
		dotpos -= fnz;
	      }
	      dotpos += evalue;
	      if(dotpos <= 0) {	
	        fputc('0', f); fputc('.', f);
		max = 0 - dotpos;
		for(i = 0; i < max; i++) { fputc('0', f); }
		fputs(buffer, f);
	      } else {		
	        sl = max = strlen(buffer);
		if(dotpos > max) { max = dotpos; }
		for(i = 0; i < max; i++) {
		  if(i < sl) {
		    fputc(buffer[i], f);
		  } else {
		    fputc('0', f);
		  }
		  if(((i + 1) == dotpos) && (dotpos < sl)) {
		    fputc('.', f);
		  }
		}
		/* if(dotpos >= sl) { fputc('0', f); } */
	      }
	    } else {		
	      fputc('0', f);
	    }
	  } else {		
	    if(p1) { fputs(p1, f); }
	    if(p2) { fputc('.', f); fputs(p2, f); }
	    fprintf(f, "e%d", evalue);
	  }
	} else {		
	  fputs(p1, f);
	}
      } else {			
        fputs(p1, f);
      }
    } else {			
      fputs(p1, f);
    }
  } else {			
    fputc('0', f);
  }
}


void
dkma_fputs_double_no_exp DK_P2(FILE *,f, double,d) {
  char buffer[64];
  sprintf(buffer, "%lg", d);
  dkma_fputs_double_str_no_exp(f, buffer);
}


