Annotation of brogue-ce/src/brogue/Math.c, Revision 1.1
1.1 ! rubenllo 1: /*
! 2: * RogueMain.c
! 3: * Brogue
! 4: *
! 5: * Created by Brian Walker on 12/26/08.
! 6: * Copyright 2012. All rights reserved.
! 7: *
! 8: * This file is part of Brogue.
! 9: *
! 10: * This program is free software: you can redistribute it and/or modify
! 11: * it under the terms of the GNU Affero General Public License as
! 12: * published by the Free Software Foundation, either version 3 of the
! 13: * License, or (at your option) any later version.
! 14: *
! 15: * This program is distributed in the hope that it will be useful,
! 16: * but WITHOUT ANY WARRANTY; without even the implied warranty of
! 17: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! 18: * GNU Affero General Public License for more details.
! 19: *
! 20: * You should have received a copy of the GNU Affero General Public License
! 21: * along with this program. If not, see <http://www.gnu.org/licenses/>.
! 22: */
! 23:
! 24: #include <time.h>
! 25: #include <limits.h>
! 26: #include <stdint.h> // C99
! 27:
! 28: #include "Rogue.h"
! 29: #include "IncludeGlobals.h"
! 30:
! 31:
! 32: // Random number generation
! 33:
! 34: short randClump(randomRange theRange) {
! 35: return randClumpedRange(theRange.lowerBound, theRange.upperBound, theRange.clumpFactor);
! 36: }
! 37:
! 38: // Get a random int between lowerBound and upperBound, inclusive, with probability distribution
! 39: // affected by clumpFactor.
! 40: short randClumpedRange(short lowerBound, short upperBound, short clumpFactor) {
! 41: if (upperBound <= lowerBound) {
! 42: return lowerBound;
! 43: }
! 44: if (clumpFactor <= 1) {
! 45: return rand_range(lowerBound, upperBound);
! 46: }
! 47:
! 48: short i, total = 0, numSides = (upperBound - lowerBound) / clumpFactor;
! 49:
! 50: for(i=0; i < (upperBound - lowerBound) % clumpFactor; i++) {
! 51: total += rand_range(0, numSides + 1);
! 52: }
! 53:
! 54: for(; i< clumpFactor; i++) {
! 55: total += rand_range(0, numSides);
! 56: }
! 57:
! 58: return (total + lowerBound);
! 59: }
! 60:
! 61: // Get a random int between lowerBound and upperBound, inclusive
! 62: boolean rand_percent(short percent) {
! 63: return (rand_range(0, 99) < clamp(percent, 0, 100));
! 64: }
! 65:
! 66: void shuffleList(short *list, short listLength) {
! 67: short i, r, buf;
! 68: for (i=0; i<listLength; i++) {
! 69: r = rand_range(0, listLength-1);
! 70: if (i != r) {
! 71: buf = list[r];
! 72: list[r] = list[i];
! 73: list[i] = buf;
! 74: }
! 75: }
! 76: }
! 77:
! 78: void fillSequentialList(short *list, short listLength) {
! 79: short i;
! 80: for (i=0; i<listLength; i++) {
! 81: list[i] = i;
! 82: }
! 83: }
! 84:
! 85: //typedef unsigned long int u4;
! 86: typedef uint32_t u4;
! 87: typedef struct ranctx { u4 a; u4 b; u4 c; u4 d; } ranctx;
! 88:
! 89: static ranctx RNGState[2];
! 90:
! 91: #define rot(x,k) (((x)<<(k))|((x)>>(32-(k))))
! 92: u4 ranval( ranctx *x ) {
! 93: u4 e = x->a - rot(x->b, 27);
! 94: x->a = x->b ^ rot(x->c, 17);
! 95: x->b = x->c + x->d;
! 96: x->c = x->d + e;
! 97: x->d = e + x->a;
! 98: return x->d;
! 99: }
! 100:
! 101: void raninit( ranctx *x, u4 seed ) {
! 102: u4 i;
! 103: x->a = 0xf1ea5eed, x->b = x->c = x->d = seed;
! 104: for (i=0; i<20; ++i) {
! 105: (void)ranval(x);
! 106: }
! 107: }
! 108:
! 109: /* ----------------------------------------------------------------------
! 110: range
! 111:
! 112: returns a number between 0 and N-1
! 113: without any bias.
! 114:
! 115: */
! 116:
! 117: #define RAND_MAX_COMBO ((unsigned long) UINT32_MAX)
! 118:
! 119: long range(long n, short RNG) {
! 120: unsigned long div;
! 121: long r;
! 122:
! 123: div = RAND_MAX_COMBO/n;
! 124:
! 125: do {
! 126: r = ranval(&(RNGState[RNG])) / div;
! 127: } while (r >= n);
! 128:
! 129: return r;
! 130: }
! 131:
! 132: // Get a random int between lowerBound and upperBound, inclusive, with uniform probability distribution
! 133:
! 134: #ifdef AUDIT_RNG // debug version
! 135: long rand_range(long lowerBound, long upperBound) {
! 136: int retval;
! 137: char RNGMessage[100];
! 138: if (upperBound <= lowerBound) {
! 139: return lowerBound;
! 140: }
! 141: long interval = upperBound - lowerBound + 1;
! 142: brogueAssert(interval > 1); // to verify that we didn't wrap around
! 143: retval = lowerBound + range(interval, rogue.RNG);
! 144: if (rogue.RNG == RNG_SUBSTANTIVE) {
! 145: randomNumbersGenerated++;
! 146: if (1) { //randomNumbersGenerated >= 1128397) {
! 147: sprintf(RNGMessage, "\n#%lu, %ld to %ld: %ld", randomNumbersGenerated, lowerBound, upperBound, retval);
! 148: RNGLog(RNGMessage);
! 149: }
! 150: }
! 151: return retval;
! 152: }
! 153: #else // normal version
! 154: long rand_range(long lowerBound, long upperBound) {
! 155: if (upperBound <= lowerBound) {
! 156: return lowerBound;
! 157: }
! 158: if (rogue.RNG == RNG_SUBSTANTIVE) {
! 159: randomNumbersGenerated++;
! 160: }
! 161: long interval = upperBound - lowerBound + 1;
! 162: brogueAssert(interval > 1); // to verify that we didn't wrap around
! 163: return lowerBound + range(interval, rogue.RNG);
! 164: }
! 165: #endif
! 166:
! 167: // seeds with the time if called with a parameter of 0; returns the seed regardless.
! 168: // All RNGs are seeded simultaneously and identically.
! 169: unsigned long seedRandomGenerator(unsigned long seed) {
! 170: if (seed == 0) {
! 171: seed = (unsigned long) time(NULL) - 1352700000;
! 172: }
! 173: raninit(&(RNGState[RNG_SUBSTANTIVE]), seed);
! 174: raninit(&(RNGState[RNG_COSMETIC]), seed);
! 175: return seed;
! 176: }
! 177:
! 178:
! 179: // Fixed-point arithmetic
! 180:
! 181: fixpt fp_round(fixpt x) {
! 182: long long div = x / FP_FACTOR, rem = x % FP_FACTOR;
! 183: int sign = (x >= 0) - (x < 0);
! 184:
! 185: if (rem >= FP_FACTOR / 2 || rem <= -FP_FACTOR / 2) {
! 186: return div + sign;
! 187: } else {
! 188: return div;
! 189: }
! 190: }
! 191:
! 192: // Returns the bit position of the most significant bit of x, where the unit
! 193: // bit has position 1. Returns 0 if x=0.
! 194: static int msbpos(unsigned long long x) {
! 195: if (x == 0) return 0;
! 196: int n = 0;
! 197: do {
! 198: n += 1;
! 199: } while (x >>= 1);
! 200: return n;
! 201: }
! 202:
! 203: static fixpt fp_exp2(int n) {
! 204: return (n >= 0 ? FP_FACTOR << n : FP_FACTOR >> -n);
! 205: }
! 206:
! 207: // Calculates sqrt(u) using the bisection method to find the root of
! 208: // f(x) = x^2 - u.
! 209: fixpt fp_sqrt(fixpt u) {
! 210:
! 211: static const fixpt SQUARE_ROOTS[] = { // values were computed by the code that follows
! 212: 0, 65536, 92682, 113511, 131073, 146543, 160529, 173392, 185363, 196608, 207243, 217359, 227023, 236293, 245213, 253819,
! 213: 262145, 270211, 278045, 285665, 293086, 300323, 307391, 314299, 321059, 327680, 334169, 340535, 346784, 352923, 358955, 364889,
! 214: 370727, 376475, 382137, 387717, 393216, 398640, 403991, 409273, 414487, 419635, 424721, 429749, 434717, 439629, 444487, 449293,
! 215: 454047, 458752, 463409, 468021, 472587, 477109, 481589, 486028, 490427, 494786, 499107, 503391, 507639, 511853, 516031, 520175,
! 216: 524289, 528369, 532417, 536435, 540423, 544383, 548313, 552217, 556091, 559939, 563762, 567559, 571329, 575077, 578799, 582497,
! 217: 586171, 589824, 593453, 597061, 600647, 604213, 607755, 611279, 614783, 618265, 621729, 625173, 628599, 632007, 635395, 638765,
! 218: 642119, 645455, 648773, 652075, 655360, 658629, 661881, 665117, 668339, 671545, 674735, 677909, 681071, 684215, 687347, 690465,
! 219: 693567, 696657, 699733, 702795, 705845, 708881, 711903, 714913, 717911, 720896, 723869, 726829, 729779, 732715, 735639, 738553
! 220: };
! 221:
! 222: if (u < 0) return -fp_sqrt(-u);
! 223:
! 224: if ((u & (127LL << FP_BASE)) == u) {
! 225: // u is an integer between 0 and 127
! 226: return SQUARE_ROOTS[u >> FP_BASE];
! 227: }
! 228:
! 229: // Find the unique k such that 2^(k-1) <= u < 2^k
! 230: // FP_BASE is the msbpos-1 of FP_FACTOR ("one")
! 231: int k = msbpos(u) - FP_BASE;
! 232:
! 233: fixpt x = 0, fx, upper, lower;
! 234: // Since 2^(k-1) <= u < 2^k, we have 2^(ceil(k/2)-1) <= sqrt(u) < 2^ceil(k/2).
! 235: // First ineq. from sqrt(u) >= 2^[(k-1)/2] = 2^[k/2 + 1/2 - 1] >= 2^(ceil(k/2) - 1)
! 236: // To calculate ceil(k/2), do k/2 but add 1 to k if positive.
! 237: upper = fp_exp2((k + (k > 0))/2);
! 238: lower = upper / 2;
! 239:
! 240: while (upper != lower + 1) {
! 241: x = (upper + lower) / 2;
! 242: fx = FP_MUL(x, x) - u;
! 243:
! 244: if (fx == 0) {
! 245: break;
! 246: } else if (fx > 0) {
! 247: upper = x;
! 248: } else {
! 249: lower = x;
! 250: }
! 251: }
! 252:
! 253: return x;
! 254: }
! 255:
! 256: // Returns base to the power of expn
! 257: fixpt fp_pow(fixpt base, int expn) {
! 258: if (base == 0) return 0;
! 259:
! 260: if (expn < 0) {
! 261: base = FP_DIV(FP_FACTOR, base);
! 262: expn = -expn;
! 263: }
! 264:
! 265: fixpt res = FP_FACTOR, err = 0;
! 266: while (expn--) {
! 267: res = res * base + (err * base) / FP_FACTOR;
! 268: err = res % FP_FACTOR;
! 269: res /= FP_FACTOR;
! 270: }
! 271:
! 272: return res + fp_round(err);
! 273: }
CVSweb
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to Math.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [contributed] / brogue-ce / src / brogue