round.c

// Rounding methods implementation in C
// 2014 kyontan No Rights Reserved (CC0)
// https://creativecommons.org/publicdomain/zero/1.0/deed.ja

/* Supporting methods
 - ceil (round up)
 - floor (round down)
 - truncate (round towards zero)
 - round off (round half up)
 - JIS Z8401:1999 round (round to the nearest even)
 */

#include <stdio.h>
#include <stdlib.h>

double abs_(double n);
int sign_(double n);
double pow_(double n, int k);
double mod_(double n, double k);

double round_(double n, double precision, double (*func)(double, double));
double round_dk(double n, int d, int k, double (*func)(double, double));

// Rounding test
// flag (bitfield): 1=ceil, 2=floor, 4 = truncate, 8 = round_off, 16 = JISRound
void test(double n, double precision, short flag);
void test_dk(double n, int d, int k, short flag);

// round width: d * 10^k (0 < d && d < 10)

double Ceil(double n, double precision);
double Floor(double n, double precision);
double Truncate(double n, double precision);
double Round_off(double n, double precision);
double JISround(double n, double precision); // JIS Z8401:1999

int main(int argc, char **argv) {

  double n, p;
  int d = 1, k;

  switch(argc) {
    case 2:
      {
        // Rounding test
        double t[] = { 1.35, 1.4, 1.44, 1.45, 1.46, 1.5, 1.55 };

        for(int f=0; f<5; f++) {
          for(int i=0; i<7; i++)
            test(t[6-i], 0.1, 1 << f);
          for(int i=0; i<7; i++)
            test(-t[i], 0.1, 1 << f);
        }
      }
      break;

    case 3:
      n = atof(argv[1]);
      p = atof(argv[2]);
      break;

    case 4:
      n = atof(argv[1]);
      d = atof(argv[2]);
      k = atof(argv[3]);
      p = d * pow_(10, k);
      break;

    default:
      printf("This program test's round methods.\n"
             "Round n by d * 10^k (0 < d < 10)\n");

      printf("Input n, d, k: ");
      scanf("%lf, %d, %d", &n, &d, &k);
      p = d * pow_(10, k);
      break;
  }

  if (d < 1 || 9 < d) {
    printf("d must be range between 1..9\n");
    return -1;
  }

  test(n, p, 31);

  return 0;
}

double abs_(double n) {
  return sign_(n) * n;
}

int sign_(double n) {
  return (n < 0) ? -1 : 1;
}

double pow_(double n, int k) {
  double x = 1.0;

  int times = abs_(k);
  double mul = (0 < sign_(k)) ? n : 1.0/n;

  for(int i=0; i<times; i++)
    x *= mul;

  return x;
}

double mod_(double n, double k) {
  if (k < 1)
    return k*(abs_(n/k) - (int)abs_(n/k));
  else
    return (abs_(n) - k*(int)(abs_(n)/k));
}

double round_(double n, double precision, double (*func)(double, double)) {
  if(precision < 0)
    precision *= -1;

  return func(n, precision);
}

double round_dk(double n, int d, int k, double (*func)(double, double)) {
  return round_(n, d * pow_(10, k), func);
}

void test(double n, double precision, short flag) {
  if (flag & 1 << 0)
    printf("     ceil(% lf, %lf) = % lf\n", n, precision, round_(n, precision, Ceil));

  if (flag & 1 << 1)
    printf("    floor(% lf, %lf) = % lf\n", n, precision, round_(n, precision, Floor));

  if (flag & 1 << 2)
    printf(" truncate(% lf, %lf) = % lf\n", n, precision, round_(n, precision, Truncate));

  if (flag & 1 << 3)
    printf("round_off(% lf, %lf) = % lf\n", n, precision, round_(n, precision, Round_off));

  if (flag & 1 << 4)
    printf(" JISround(% lf, %lf) = % lf\n", n, precision, round_(n, precision, JISround));
}

void test_dk(double n, int d, int k, short flag) {
  test(n, d * pow_(10, k), flag);
}


double Ceil(double n, double precision) {
  double m = mod_(n, precision);

  if (n < 0 || m == 0.0)
    return n + m;
  else
    return n + (precision - m);
}

double Floor(double n, double precision) {
  double m = mod_(n, precision);

  if (0 < n || m == 0.0)
    return n - m;
  else
    return n - (precision - mod_(n, precision));
}

double Truncate(double n, double precision) {
  return sign_(n) * Floor(abs_(n), precision);
}

double Round_off(double n, double precision) {
  return sign_(n) * Floor(abs_(n) + precision * 0.5, precision);
}

double JISround(double n, double precision) {
  double  c = Ceil(n, precision),
          f = Floor(n, precision),
          na = abs_(n),
          ca = abs_(c),
          fa = abs_(f);

  // fa <= na <= ca
  // abs_((ca-na)-(na-fa))
  if (abs_(ca+fa-2*na) > 1e-8)
    return Round_off(n, precision);

  if ((int)(ca*(1.0/precision)) % 2 == 1)
    return f;
  else
    return c;
}