/* ecvt, fcvt -- convert floating-point numbers to ascii	-*- C -*-	*/

/* Copyright (c) 2009 Ian Piumarta
 * 
 * All rights reserved.
 * 
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the 'Software'), to
 * deal in the Software without restriction, including without limitation the
 * rights to use, copy, modify, merge, publish, distribute, and/or sell copies
 * of the Software, and to permit persons to whom the Software is furnished to
 * do so, provided that the above copyright notice(s) and this permission
 * notice appear in all copies of the Software.  Inclusion of the above
 * copyright notice(s) in supporting documentation would be appreciated, but is
 * not required.
 *
 * THE SOFTWARE IS PROVIDED 'AS IS'.  USE ENTIRELY AT YOUR OWN RISK.
 */

/* This file provides replacements for the functions ecvt() and fcvt() that
 * were recently deprecated in POSIX.  The interface and behaviour of these
 * functions is identical to the POSIX ones that they replace, except that:
 * 
 *   - the broken behaviour for -inf, -nan and -0.0 in glibc is not reproduced;
 * 
 *   - binary integers and fractions that can be represented precisely in the
 *     available number of sinificant binary digits (regardless of binary
 *     exponent) can are reproduced precisely in decimal, given sufficient
 *     significant decimal digits (glibc improperly rounds the fraction after
 *     15 decimal digits or so);
 * 
 *   - on my machines they run at about half the speed of the ones in libc.
 * 
 * They have been tested on 32- and 64-bit machines of both orientations.
 * 
 * For details on the use of these functions, see your ecvt(3) manual page.  If
 * you don't have one handy, there might still be one available here:
 * http://opengroup.org/onlinepubs/007908799/xsh/ecvt.html
 */

#include <stdint.h>
#include <sys/types.h>
#include <stdlib.h>
#include <string.h>

#define MP_NO_DECLS	1
#define MP_API		static
#define MP_INL		static inline

#include "mpint.h"
#include "_mpint.h"
#include "mp_clear.c"
#include "mp_normalise.c"
#include "mp_digitAtPut.c"
#include "mp_copy.c"
#include "mp_cmp.c"
#include "mp_shiftl.c"
#include "mp_shiftr.c"
#include "mp_shift.c"
#include "mp_mul_d.c"
#include "mp_div_d.c"

static inline void prepend(char **buf, int *len, int *pos, int c)
{
  if (--*pos < 0) {
    *buf= *buf ? realloc(*buf, *len += 32) : malloc(*len= 32);
    memmove(*buf + 32, *buf, *len - 32);
    *pos += 32;
  }
  (*buf)[*pos]= c;
}

static inline void append(char **buf, int *len, int *pos, int c)
{
  if (*pos >= *len) {
    *buf= *buf ? realloc(*buf, *len += 32) : malloc(*len= 32);
  }
  (*buf)[(*pos)++]= c;
}

static char *fconvert(double value, int ndigit, int *decpt, int *sign, int fflag)
{
  static char *buf= 0;
  static int buflen= 0;
  int pos= 0;
  union { uint32_t l;  float f; } x;
  x.f= value;
  int exp2= (0xff & (x.l >> 23)) - 127;
  uint32_t mant= x.l << 9;
  *sign= x.l >> 31;
  *decpt= 0;
  if (exp2 == 0x80) return mant ? "nan" : "inf";
  int norm= 1;
  if (exp2 == -127) {
    norm= 0;
    exp2= -126;
  }
  mp_t mp= MP_INITIALISER, tmp= MP_INITIALISER;
  mp_digitAtPut(&mp, 1, norm);
  mp_digitAtPut(&mp, 0, mant);
  mp_shift(&mp, exp2 - 32);
  int p= buflen;
  while (!mp_zerop(&mp)) prepend(&buf, &buflen, &p, '0' + mp_div_d(&mp, 10));
  *decpt= pos= buflen - p;
  memmove(buf, buf + p, pos);
  if (fflag) ndigit += pos;
  if (pos <= ndigit) {
    size_t index= _mp_max(0, 1 + ((23 - exp2) / 32));
    mp_clear(&mp);
    mp_digitAtPut(&mp, index, norm);
    if (index > 0) mp_digitAtPut(&mp, index - 1, mant);
    mp_shift(&mp, exp2);
    mp_digitAtPut(&mp, index, 0);
    if (!mp_zerop(&mp)) {
      mp_mul_d(&mp, 10);
      if (!pos)
	while (!mp_digitAt(&mp, index)) {
	  mp_mul_d(&mp, 10);
	  --*decpt;
	  if (fflag && !--ndigit) break;
	}
      while (pos <= ndigit && !mp_zerop(&mp)) {
	mp_digit_t r= mp_digitAt(&mp, index);
	append(&buf, &buflen, &pos, '0' + r);
	if (r) _mp_digitAtPut(&mp, index, 0);
	mp_mul_d(&mp, 10);
      }
    }
    if (!pos && !fflag) ++*decpt;
    while (pos <= ndigit) append(&buf, &buflen, &pos, '0');
  }
  append(&buf, &buflen, &pos, 0);
  append(&buf, &buflen, &pos, 0);
  if (buf[ndigit] >= '5') {
    for (p= ndigit;  p-- && ++buf[p] > '9';  buf[p]= '0');
    if (p < 0) {
      memmove(buf + 1, buf, ndigit);
      buf[0]= '1';
      ++*decpt;
    }
  }
  buf[ndigit]= 0;
  mp_clear(&tmp);
  mp_clear(&mp);
  return buf;
}

static char *dconvert(double value, int ndigit, int *decpt, int *sign, int fflag)
{
  static char *buf= 0;
  static int buflen= 0;
  int pos= 0;
  union { uint64_t l;  double f; } x;
  x.f= value;
  int exp2= (0x7ff & (x.l >> 52)) - 1023;
  uint64_t mant= x.l << 12;
  *sign= x.l >> 63;
  *decpt= 0;
  if (exp2 == 0x400) return mant ? "nan" : "inf";
  int norm= 1;
  if (exp2 == -1023) {
    norm= 0;
    exp2= -1022;
  }
  mp_t mp= MP_INITIALISER, tmp= MP_INITIALISER;
  mp_digitAtPut(&mp, 2, norm);
  mp_digitAtPut(&mp, 1, mant >> 32);
  mp_digitAtPut(&mp, 0, mant & 0xffffffff);
  mp_shift(&mp, exp2 - 64);
  int p= buflen;
  while (!mp_zerop(&mp)) prepend(&buf, &buflen, &p, '0' + mp_div_d(&mp, 10));
  *decpt= pos= buflen - p;
  memmove(buf, buf + p, pos);
  if (fflag) ndigit += pos;
  if (pos <= ndigit) {
    size_t index= _mp_max(2, 2 - ((exp2 - 20) / 32));
    mp_clear(&mp);
    mp_digitAtPut(&mp, index,     norm);
    mp_digitAtPut(&mp, index - 1, mant >> 32);
    mp_digitAtPut(&mp, index - 2, mant & 0xffffffff);
    mp_shift(&mp, exp2);
    mp_digitAtPut(&mp, index, 0);
    if (!mp_zerop(&mp)) {
      mp_mul_d(&mp, 10);
      if (!pos)
	while (!mp_digitAt(&mp, index)) {
	  mp_mul_d(&mp, 10);
	  --*decpt;
	  if (fflag && !--ndigit) break;
	}
      while (pos <= ndigit && !mp_zerop(&mp)) {
	mp_digit_t r= mp_digitAt(&mp, index);
	append(&buf, &buflen, &pos, '0' + r);
	if (r) _mp_digitAtPut(&mp, index, 0);
	mp_mul_d(&mp, 10);
      }
    }
    if (!pos &&!fflag) {
      ++*decpt;
    }
    while (pos <= ndigit)
      append(&buf, &buflen, &pos, '0');
  }
  append(&buf, &buflen, &pos, 0);
  append(&buf, &buflen, &pos, 0);
  if (buf[ndigit] >= '5') {
    for (p= ndigit;  p-- && ++buf[p] > '9';  buf[p]= '0');
    if (p < 0) {
      memmove(buf + 1, buf, ndigit);
      buf[0]= '1';
      ++*decpt;
    }
  }
  buf[ndigit]= 0;
  mp_clear(&tmp);
  mp_clear(&mp);
  return buf;
}

char *fcvtf(float in, int ndigit, int *decpt, int *sign)	{ return fconvert(in, ndigit, decpt, sign, 1); }
char *ecvtf(float in, int ndigit, int *decpt, int *sign)	{ return fconvert(in, ndigit, decpt, sign, 0); }

char *fcvtd(double in, int ndigit, int *decpt, int *sign)	{ return dconvert(in, ndigit, decpt, sign, 1); }
char *ecvtd(double in, int ndigit, int *decpt, int *sign)	{ return dconvert(in, ndigit, decpt, sign, 0); }

#if defined(TEST_FCVT)

#include <stdio.h>
#include <math.h>

int cmp(char *a, char *b)
{
  int i= 0;
  while (a[i] == b[i])
    if (!a[i])	return 1;
    else	++i;
  if      (a[i] == b[i] + 1) while (b[i + 1] == '9' || b[i + 1] == '.') ++i;
  else if (b[i] == a[i] + 1) while (a[i + 1] == '9' || a[i + 1] == '.') ++i;
  if (!a[i + 1] && !b[i + 1]) return 1;
  printf("%d <%s>\n", i, a + i);
  printf("%d <%s>\n", i, b + i);
  return 0;
}

int fcmp(double value, int ossign, int osdecpt, char *osbuf, int mysign, int mydecpt, char *mybuf)
{
  if (!ossign && mysign && !strcmp(osbuf + 1, mybuf))	// glibc is wrong for -inf
    return isinf(value);
  if (!ossign && mysign && !strcmp(osbuf, mybuf))	// glibc is wrong for -0.0 and -nan
    return isnan(value) || !value;
  if (osdecpt == mydecpt + 1)				// glibc has redundant 0. for +/-0.0
    return !value && !strcmp(osbuf + 1, mybuf);
  return ossign == mysign && osdecpt == mydecpt && cmp(osbuf, mybuf);
}

#define METHOD 1	/* 0 = os fcvt;  1 = my fcvt;  2 = none, for timing baseline */

void test(double f, int verbose)
{
  int   ossign, osdecpt, mysign, mydecpt;
  char *osbuf, *mybuf;
  osbuf= strdup(ecvt(f, 15, &osdecpt, &ossign));
  switch (METHOD) {
    case 0:	mybuf= ecvt (f, 15, &mydecpt, &mysign);			break;
    case 1:	mybuf= ecvtd(f, 15, &mydecpt, &mysign);			break;
    default:	mybuf= osbuf;  mydecpt= osdecpt;  mysign= ossign;	break;
  }
  int ok= 0;
  ok= fcmp(f, ossign, osdecpt, osbuf, mysign, mydecpt, mybuf);
  if (verbose || !ok) {
    if (!ok) printf("FAIL ECVT %.50f\n", f);
    printf("e %d\t%d\t<%s>\n", ossign, osdecpt, osbuf);
    printf("e %d\t%d\t<%s>\n", mysign, mydecpt, mybuf);
  }
  if (!ok) abort();
  free(osbuf);
  osbuf= strdup(fcvt(f, 15, &osdecpt, &ossign));
  switch(METHOD) {
    case 0:	mybuf= fcvt (f, 15, &mydecpt, &mysign);			break;
    case 1:	mybuf= fcvtd(f, 15, &mydecpt, &mysign);			break;
    default:	mybuf= osbuf;  mydecpt= osdecpt;  mysign= ossign;	break;
  }
  ok= fcmp(f, ossign, osdecpt, osbuf, mysign, mydecpt, mybuf);
  if (verbose || !ok) {
    if (!ok) printf("FAIL FCVT %.25f\n", f);
    printf("f %d\t%d\t<%s>\n", ossign, osdecpt, osbuf);
    printf("f %d\t%d\t<%s>\n", mysign, mydecpt, mybuf);
  }
  if (!ok) abort();
  free(osbuf);
}

#include <unistd.h>

int main(int argc, char **argv)
{
  if (argc > 1)
    while (argc > 1) {
      --argc;
      ++argv;
      test(atof(*argv), 1);
    }
  else {
    srandom(getpid());
    int i;
    for (i= 0;  i < 100000;  ++i) {
      if (0 == i % 1000) { printf("  %i\r", i);  fflush(stdout); }
      union { uint32_t l[4];  float f;  double d; } x;
      x.l[0]= (random() << 16) ^ random();
      x.l[1]= (random() << 16) ^ random();
      x.l[2]= (random() << 16) ^ random();
      x.l[3]= (random() << 16) ^ random();
      test(x.f, 0);
      test(x.d, 0);
    }
    printf("  %i\nok\n", i);
  }

  return 0;
}

#endif /* defined(TEST_FCVT) */