math.h (3293B)
1 #define SWAP(x, y) do { x ^= y; y ^= x; x ^= y; } while (0) 2 #define MIN(x, y) ((x) < (y) ? (x) : (y)) 3 #define MAX(x, y) ((x) > (y) ? (x) : (y)) 4 #define DIV_ROUND_UP(n, d) (((n) + (d) - 1) / (d)) 5 #define POSITIVE_MOD(x, y) (((x) % (y) + (y)) % (y)) 6 7 STATIC uint64_t permtoindex(size_t, const uint8_t *); 8 STATIC void indextoperm(uint64_t, size_t, uint8_t *); 9 STATIC int permsign(size_t, const uint8_t *); 10 STATIC uint64_t digitstosumzero(size_t, const uint8_t *, uint8_t); 11 STATIC void sumzerotodigits(uint64_t, size_t, uint8_t, uint8_t *); 12 STATIC double intpow(double, uint64_t); 13 14 /* This code is only used for assertions in debug mode */ 15 #ifdef DEBUG 16 STATIC bool isperm(size_t, const uint8_t *); 17 STATIC bool 18 isperm(size_t n, const uint8_t *a) 19 { 20 size_t i; 21 bool aux[FACTORIAL_MAX+1]; 22 23 memset(aux, false, n); 24 25 for (i = 0; i < n; i++) { 26 if (a[i] >= n) 27 return false; 28 else 29 aux[a[i]] = true; 30 } 31 32 for (i = 0; i < n; i++) 33 if (!aux[i]) 34 return false; 35 36 return true; 37 } 38 #endif 39 40 STATIC uint64_t 41 permtoindex(size_t n, const uint8_t *a) 42 { 43 size_t i, j; 44 uint64_t c, ret; 45 46 DBG_ASSERT(n <= FACTORIAL_MAX, "Error: cannot compute permtoindex() " 47 "for set of size %zu > %" PRIu64 "\n", n, FACTORIAL_MAX); 48 DBG_ASSERT(isperm(n, a), "Error: cannot compute permtoindex() for " 49 "invalid permutation\n"); 50 51 for (i = 0, ret = 0; i < n; i++) { 52 for (j = i+1, c = 0; j < n; j++) 53 c += a[i] > a[j]; 54 ret += factorial[n-i-1] * c; 55 } 56 57 return ret; 58 } 59 60 STATIC void 61 indextoperm(uint64_t p, size_t n, uint8_t *r) 62 { 63 uint64_t c, k; 64 size_t i, j, used; 65 66 DBG_ASSERT(n <= FACTORIAL_MAX, "Error: cannot compute indextoperm() " 67 "for set of size %zu > %" PRIu64 "\n", n, FACTORIAL_MAX); 68 DBG_ASSERT(p < factorial[n], "Error: invalid permutation index %" 69 PRIu64 " for set of size %zu\n", p, n); 70 71 for (i = 0, used = 0; i < n; i++) { 72 k = p / factorial[n-i-1]; 73 74 /* Find k-th unused number */ 75 for (j = 0, c = 0; c <= k; j++) 76 c += 1 - ((used & (1<<j)) >> j); 77 78 r[i] = j-1; 79 used |= 1 << (j-1); 80 p %= factorial[n-i-1]; 81 } 82 83 return; 84 } 85 86 STATIC int 87 permsign(size_t n, const uint8_t *a) 88 { 89 size_t i, j, ret; 90 91 for (i = 0, ret = 0; i < n; i++) 92 for (j = i+1; j < n; j++) 93 ret += a[i] > a[j]; 94 95 return ret % 2; 96 } 97 98 STATIC uint64_t 99 digitstosumzero(size_t n, const uint8_t *a, uint8_t b) 100 { 101 uint64_t ret, p; 102 size_t i, sum; 103 104 DBG_ASSERT((n == 8 && b == 3 ) || (n == 12 && b == 2), 105 "Error: digitstosumzero() called with n=%zu and b=%" PRIu8 106 " (use n=8 b=3 or n=12 b=2)\n", n, b); 107 108 for (i = 1, ret = 0, p = 1, sum = 0; i < n; i++, p *= (uint64_t)b) { 109 DBG_ASSERT(a[i] < b, "Error: digit %" PRIu8 110 " > %" PRIu8 "in digitstosumzero()\n", a[i], b); 111 sum += a[i]; 112 ret += p * (uint64_t)a[i]; 113 } 114 115 return ret; 116 } 117 118 STATIC void 119 sumzerotodigits(uint64_t d, size_t n, uint8_t b, uint8_t *a) 120 { 121 uint8_t sum; 122 size_t i; 123 124 DBG_ASSERT((n == 8 && b == 3 ) || (n == 12 && b == 2), 125 "Error: sumzerotodigits() called with n=%zu and b=%" PRIu8 126 " (use n=8 b=3 or n=12 b=2)\n", n, b); 127 128 for (i = 1, sum = 0; i < n; i++, d /= (uint64_t)b) { 129 a[i] = (uint8_t)(d % (uint64_t)b); 130 sum += a[i]; 131 } 132 a[0] = (b - (sum % b)) % b; 133 134 return; 135 } 136 137 STATIC double 138 intpow(double b, uint64_t e) 139 { 140 double r; 141 142 if (e == 0) 143 return 1.0; 144 145 r = intpow(b, e/2); 146 147 return e % 2 == 0 ? r * r : b * r * r; 148 }