stb/stb_divide.h

   1 // stb_divide.h - v0.91 - public domain - Sean Barrett, Feb 2010
   2 // Three kinds of divide/modulus of signed integers.
   3 //
   4 // HISTORY
   5 //
   6 //   v0.91  2010-02-27  Fix euclidean division by INT_MIN for non-truncating C
   7 //                      Check result with 64-bit math to catch such cases
   8 //   v0.90  2010-02-24  First public release
   9 //
  10 // USAGE
  11 //
  12 // In *ONE* source file, put:
  13 //
  14 //    #define STB_DIVIDE_IMPLEMENTATION
  15 //    // #define C_INTEGER_DIVISION_TRUNCATES  // see Note 1
  16 //    // #define C_INTEGER_DIVISION_FLOORS     // see Note 2
  17 //    #include "stb_divide.h"
  18 //
  19 // Other source files should just include stb_divide.h
  20 //
  21 // Note 1: On platforms/compilers that you know signed C division
  22 // truncates, you can #define C_INTEGER_DIVISION_TRUNCATES.
  23 //
  24 // Note 2: On platforms/compilers that you know signed C division
  25 // floors (rounds to negative infinity), you can #define
  26 // C_INTEGER_DIVISION_FLOORS.
  27 //
  28 // You can #define STB_DIVIDE_TEST in which case the implementation
  29 // will generate a main() and compiling the result will create a
  30 // program that tests the implementation. Run it with no arguments
  31 // and any output indicates an error; run it with any argument and
  32 // it will also print the test results. Define STB_DIVIDE_TEST_64
  33 // to a 64-bit integer type to avoid overflows in the result-checking
  34 // which give false negatives.
  35 //
  36 // ABOUT
  37 //
  38 // This file provides three different consistent divide/mod pairs
  39 // implemented on top of arbitrary C/C++ division, including correct
  40 // handling of overflow of intermediate calculations:
  41 //
  42 //     trunc:   a/b truncates to 0,           a%b has same sign as a
  43 //     floor:   a/b truncates to -inf,        a%b has same sign as b
  44 //     eucl:    a/b truncates to sign(b)*inf, a%b is non-negative
  45 //
  46 // Not necessarily optimal; I tried to keep it generally efficient,
  47 // but there may be better ways.
  48 //
  49 // Briefly, for those who are not familiar with the problem, we note
  50 // the reason these divides exist and are interesting:
  51 //
  52 //     'trunc' is easy to implement in hardware (strip the signs,
  53 //          compute, reapply the signs), thus is commonly defined
  54 //          by many languages (including C99)
  55 //
  56 //     'floor' is simple to define and better behaved than trunc;
  57 //          for example it divides integers into fixed-size buckets
  58 //          without an extra-wide bucket at 0, and for a fixed
  59 //          divisor N there are only |N| possible moduli.
  60 //
  61 //     'eucl' guarantees fixed-sized buckets *and* a non-negative
  62 //          modulus and defines division to be whatever is needed
  63 //          to achieve that result.
  64 //
  65 // See "The Euclidean definition of the functions div and mod"
  66 // by Raymond Boute (1992), or "Division and Modulus for Computer
  67 // Scientists" by Daan Leijen (2001)
  68 //
  69 // We assume of the built-in C division:
  70 //     (a) modulus is the remainder for the corresponding division
  71 //     (b) a/b truncates if a and b are the same sign
  72 //
  73 // Property (a) requires (a/b)*b + (a%b)==a, and is required by C.
  74 // Property (b) seems to be true of all hardware but is *not* satisfied
  75 // by the euclidean division operator we define, so it's possibly not
  76 // always true. If any such platform turns up, we can add more cases.
  77 // (Possibly only stb_div_trunc currently relies on property (b).)
  78 //
  79 // LICENSE
  80 //
  81 //   This software is dual-licensed to the public domain and under the following
  82 //   license: you are granted a perpetual, irrevocable license to copy, modify,
  83 //   publish, and distribute this file as you see fit.
  84
  85
  86 #ifndef INCLUDE_STB_DIVIDE_H
  87 #define INCLUDE_STB_DIVIDE_H
  88
  89 #ifdef __cplusplus
  90 extern "C" {
  91 #endif
  92
  93 extern int stb_div_trunc(int value_to_be_divided, int value_to_divide_by);
  94 extern int stb_div_floor(int value_to_be_divided, int value_to_divide_by);
  95 extern int stb_div_eucl (int value_to_be_divided, int value_to_divide_by);
  96 extern int stb_mod_trunc(int value_to_be_divided, int value_to_divide_by);
  97 extern int stb_mod_floor(int value_to_be_divided, int value_to_divide_by);
  98 extern int stb_mod_eucl (int value_to_be_divided, int value_to_divide_by);
  99
 100 #ifdef __cplusplus
 101 }
 102 #endif
 103
 104 #ifdef STB_DIVIDE_IMPLEMENTATION
 105
 106 #if defined(__STDC_VERSION) && __STDC_VERSION__ >= 19901
 107    #ifndef C_INTEGER_DIVISION_TRUNCATES
 108       #define C_INTEGER_DIVISION_TRUNCATES
 109    #endif
 110 #endif
 111
 112 #ifndef INT_MIN
 113 #include <limits.h> // if you have no limits.h, #define INT_MIN yourself
 114 #endif
 115
 116 // the following macros are designed to allow testing
 117 // other platforms by simulating them
 118 #ifndef STB_DIVIDE_TEST_FLOOR
 119    #define stb__div(a,b)  ((a)/(b))
 120    #define stb__mod(a,b)  ((a)%(b))
 121 #else
 122    // implement floor-style divide on trunc platform
 123    #ifndef C_INTEGER_DIVISION_TRUNCATES
 124    #error "floor test requires truncating division"
 125    #endif
 126    #undef C_INTEGER_DIVISION_TRUNCATES
 127    int stb__div(int v1, int v2)
 128    {
 129       int q = v1/v2, r = v1%v2;
 130       if ((r > 0 && v2 < 0) || (r < 0 && v2 > 0))
 131          return q-1;
 132       else
 133          return q;
 134    }
 135
 136    int stb__mod(int v1, int v2)
 137    {
 138       int r = v1%v2;
 139       if ((r > 0 && v2 < 0) || (r < 0 && v2 > 0))
 140          return r+v2;
 141       else
 142          return r;
 143    }
 144 #endif
 145
 146 int stb_div_trunc(int v1, int v2)
 147 {
 148    #ifdef C_INTEGER_DIVISION_TRUNCATES
 149    return v1/v2;
 150    #else
 151    if (v1 >= 0 && v2 <= 0)
 152       return -stb__div(-v1,v2);  // both negative to avoid overflow
 153    if (v1 <= 0 && v2 >= 0)
 154       if (v1 != INT_MIN)
 155          return -stb__div(v1,-v2);    // both negative to avoid overflow
 156       else
 157          return -stb__div(v1+v2,-v2)-1; // push v1 away from wrap point
 158    else
 159       return v1/v2;            // same sign, so expect truncation
 160    #endif
 161 }
 162
 163 int stb_div_floor(int v1, int v2)
 164 {
 165    #ifdef C_INTEGER_DIVISION_FLOORS
 166    return v1/v2;
 167    #else
 168    if (v1 >= 0 && v2 < 0)
 169       if ((-v1)+v2+1 < 0) // check if increasing v1's magnitude overflows
 170          return -stb__div(-v1+v2+1,v2); // nope, so just compute it
 171       else
 172          return -stb__div(-v1,v2) + ((-v1)%v2 ? -1 : 0);
 173    if (v1 < 0 && v2 >= 0)
 174       if (v1 != INT_MIN)
 175          if (v1-v2+1 < 0) // check if increasing v1's magnitude overflows
 176             return -stb__div(v1-v2+1,-v2); // nope, so just compute it
 177          else
 178             return -stb__div(-v1,v2) + (stb__mod(v1,-v2) ? -1 : 0);
 179       else // it must be possible to compute -(v1+v2) without overflowing
 180          return -stb__div(-(v1+v2),v2) + (stb__mod(-(v1+v2),v2) ? -2 : -1);
 181    else
 182       return v1/v2;           // same sign, so expect truncation
 183    #endif
 184 }
 185
 186 int stb_div_eucl(int v1, int v2)
 187 {
 188    int q,r;
 189    #ifdef C_INTEGER_DIVISION_TRUNCATES
 190    q = v1/v2;
 191    r = v1%v2;
 192    #else
 193    // handle every quadrant separately, since we can't rely on q and r flor
 194    if (v1 >= 0)
 195       if (v2 >= 0)
 196          return stb__div(v1,v2);
 197       else if (v2 != INT_MIN)
 198          q = -stb__div(v1,-v2), r = stb__mod(v1,-v2);
 199       else
 200          q = 0, r = v1;
 201    else if (v1 != INT_MIN)
 202       if (v2 >= 0)
 203          q = -stb__div(-v1,v2), r = -stb__mod(-v1,v2);
 204       else if (v2 != INT_MIN)
 205          q = stb__div(-v1,-v2), r = -stb__mod(-v1,-v2);
 206       else // if v2 is INT_MIN, then we can't use -v2, but we can't divide by v2
 207          q = 1, r = v1-q*v2;
 208    else // if v1 is INT_MIN, we have to move away from overflow place
 209       if (v2 >= 0)
 210          q = -stb__div(-(v1+v2),v2)-1, r = -stb__mod(-(v1+v2),v2);
 211       else
 212          q = stb__div(-(v1-v2),-v2)+1, r = -stb__mod(-(v1-v2),-v2);
 213    #endif
 214    if (r >= 0)
 215       return q;
 216    else
 217       return q + (v2 > 0 ? -1 : 1);
 218 }
 219
 220 int stb_mod_trunc(int v1, int v2)
 221 {
 222    #ifdef C_INTEGER_DIVISION_TRUNCATES
 223    return v1%v2;
 224    #else
 225    if (v1 >= 0) { // modulus result should always be positive
 226       int r = stb__mod(v1,v2);
 227       if (r >= 0)
 228          return r;
 229       else
 230          return r + (v2 > 0 ? v2 : -v2);
 231    } else {    // modulus result should always be negative
 232       int r = stb__mod(v1,v2);
 233       if (r <= 0)
 234          return r;
 235       else
 236          return r - (v2 > 0 ? v2 : -v2);
 237    }
 238    #endif
 239 }
 240
 241 int stb_mod_floor(int v1, int v2)
 242 {
 243    #ifdef C_INTEGER_DIVISION_FLOORS
 244    return v1%v2;
 245    #else
 246    if (v2 >= 0) { // result should always be positive
 247       int r = stb__mod(v1,v2);
 248       if (r >= 0)
 249          return r;
 250       else
 251          return r + v2;
 252    } else { // result should always be negative
 253       int r = stb__mod(v1,v2);
 254       if (r <= 0)
 255          return r;
 256       else
 257          return r + v2;
 258    }
 259    #endif
 260 }
 261
 262 int stb_mod_eucl(int v1, int v2)
 263 {
 264    int r = stb__mod(v1,v2);
 265
 266    if (r >= 0)
 267       return r;
 268    else
 269       return r + (v2 > 0 ? v2 : -v2); // abs()
 270 }
 271
 272 #ifdef STB_DIVIDE_TEST
 273 #include <stdio.h>
 274 #include <math.h>
 275 #include <limits.h>
 276
 277 int show=0;
 278
 279 void stbdiv_check(int q, int r, int a, int b, char *type, int dir)
 280 {
 281    if ((dir > 0 && r < 0) || (dir < 0 && r > 0))
 282       fprintf(stderr, "FAILED: %s(%d,%d) remainder %d in wrong direction\n", type,a,b,r);
 283    else
 284       if (b != INT_MIN) // can't compute abs(), but if b==INT_MIN all remainders are valid
 285          if (r <= -abs(b) || r >= abs(b))
 286             fprintf(stderr, "FAILED: %s(%d,%d) remainder %d out of range\n", type,a,b,r);
 287    #ifdef STB_DIVIDE_TEST_64
 288    {
 289       STB_DIVIDE_TEST_64 q64 = q, r64=r, a64=a, b64=b;
 290       if (q64*b64+r64 != a64)
 291          fprintf(stderr, "FAILED: %s(%d,%d) remainder %d doesn't match quotient %d\n", type,a,b,r,q);
 292    }
 293    #else
 294    if (q*b+r != a)
 295       fprintf(stderr, "FAILED: %s(%d,%d) remainder %d doesn't match quotient %d\n", type,a,b,r,q);
 296    #endif
 297 }
 298
 299 void test(int a, int b)
 300 {
 301    int q,r;
 302    if (show) printf("(%+11d,%+d) |  ", a,b);
 303    q = stb_div_trunc(a,b), r = stb_mod_trunc(a,b);
 304    if (show) printf("(%+11d,%+2d)  ", q,r); stbdiv_check(q,r,a,b, "trunc",a);
 305    q = stb_div_floor(a,b), r = stb_mod_floor(a,b);
 306    if (show) printf("(%+11d,%+2d)  ", q,r); stbdiv_check(q,r,a,b, "floor",b);
 307    q = stb_div_eucl (a,b), r = stb_mod_eucl (a,b);
 308    if (show) printf("(%+11d,%+2d)\n", q,r); stbdiv_check(q,r,a,b, "euclidean",1);
 309 }
 310
 311 void testh(int a, int b)
 312 {
 313    int q,r;
 314    if (show) printf("(%08x,%08x) |\n", a,b);
 315    q = stb_div_trunc(a,b), r = stb_mod_trunc(a,b); stbdiv_check(q,r,a,b, "trunc",a);
 316    if (show) printf("             (%08x,%08x)", q,r);
 317    q = stb_div_floor(a,b), r = stb_mod_floor(a,b); stbdiv_check(q,r,a,b, "floor",b);
 318    if (show) printf("   (%08x,%08x)", q,r);
 319    q = stb_div_eucl (a,b), r = stb_mod_eucl (a,b); stbdiv_check(q,r,a,b, "euclidean",1);
 320    if (show) printf("   (%08x,%08x)\n ", q,r);
 321 }
 322
 323 int main(int argc, char **argv)
 324 {
 325    if (argc > 1) show=1;
 326
 327    test(8,3);
 328    test(8,-3);
 329    test(-8,3);
 330    test(-8,-3);
 331    test(1,2);
 332    test(1,-2);
 333    test(-1,2);
 334    test(-1,-2);
 335    test(8,4);
 336    test(8,-4);
 337    test(-8,4);
 338    test(-8,-4);
 339
 340    test(INT_MAX,1);
 341    test(INT_MIN,1);
 342    test(INT_MIN+1,1);
 343    test(INT_MAX,-1);
 344    //test(INT_MIN,-1); // this traps in MSVC, so we leave it untested
 345    test(INT_MIN+1,-1);
 346    test(INT_MIN,-2);
 347    test(INT_MIN+1,2);
 348    test(INT_MIN+1,-2);
 349    test(INT_MAX,2);
 350    test(INT_MAX,-2);
 351    test(INT_MIN+1,2);
 352    test(INT_MIN+1,-2);
 353    test(INT_MIN,2);
 354    test(INT_MIN,-2);
 355    test(INT_MIN,7);
 356    test(INT_MIN,-7);
 357    test(INT_MIN+1,4);
 358    test(INT_MIN+1,-4);
 359
 360    testh(-7, INT_MIN);
 361    testh(-1, INT_MIN);
 362    testh(1, INT_MIN);
 363    testh(7, INT_MIN);
 364
 365    testh(INT_MAX-1, INT_MIN);
 366    testh(INT_MAX,   INT_MIN);
 367    testh(INT_MIN,   INT_MIN);
 368    testh(INT_MIN+1, INT_MIN);
 369
 370    testh(INT_MAX-1, INT_MAX);
 371    testh(INT_MAX  , INT_MAX);
 372    testh(INT_MIN  , INT_MAX);
 373    testh(INT_MIN+1, INT_MAX);
 374
 375    return 0;
 376 }
 377 #endif // STB_DIVIDE_TEST
 378 #endif // STB_DIVIDE_IMPLEMENTATION
 379 #endif // INCLUDE_STB_DIVIDE_H