misc.c 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664
  1. /*
  2. * misc.c: Miscellaneous helpful functions.
  3. */
  4. #include <assert.h>
  5. #include <ctype.h>
  6. #ifdef NO_TGMATH_H
  7. # include <math.h>
  8. #else
  9. # include <tgmath.h>
  10. #endif
  11. #include <stdlib.h>
  12. #include <string.h>
  13. #include <stdio.h>
  14. #include "puzzles.h"
  15. char MOVE_UI_UPDATE[] = "";
  16. char MOVE_NO_EFFECT[] = "";
  17. char MOVE_UNUSED[] = "";
  18. void free_cfg(config_item *cfg)
  19. {
  20. config_item *i;
  21. for (i = cfg; i->type != C_END; i++)
  22. if (i->type == C_STRING)
  23. sfree(i->u.string.sval);
  24. sfree(cfg);
  25. }
  26. void free_keys(key_label *keys, int nkeys)
  27. {
  28. int i;
  29. for(i = 0; i < nkeys; i++)
  30. sfree(keys[i].label);
  31. sfree(keys);
  32. }
  33. /*
  34. * The Mines (among others) game descriptions contain the location of every
  35. * mine, and can therefore be used to cheat.
  36. *
  37. * It would be pointless to attempt to _prevent_ this form of
  38. * cheating by encrypting the description, since Mines is
  39. * open-source so anyone can find out the encryption key. However,
  40. * I think it is worth doing a bit of gentle obfuscation to prevent
  41. * _accidental_ spoilers: if you happened to note that the game ID
  42. * starts with an F, for example, you might be unable to put the
  43. * knowledge of those mines out of your mind while playing. So,
  44. * just as discussions of film endings are rot13ed to avoid
  45. * spoiling it for people who don't want to be told, we apply a
  46. * keyless, reversible, but visually completely obfuscatory masking
  47. * function to the mine bitmap.
  48. */
  49. void obfuscate_bitmap(unsigned char *bmp, int bits, bool decode)
  50. {
  51. int bytes, firsthalf, secondhalf;
  52. struct step {
  53. unsigned char *seedstart;
  54. int seedlen;
  55. unsigned char *targetstart;
  56. int targetlen;
  57. } steps[2];
  58. int i, j;
  59. /*
  60. * My obfuscation algorithm is similar in concept to the OAEP
  61. * encoding used in some forms of RSA. Here's a specification
  62. * of it:
  63. *
  64. * + We have a `masking function' which constructs a stream of
  65. * pseudorandom bytes from a seed of some number of input
  66. * bytes.
  67. *
  68. * + We pad out our input bit stream to a whole number of
  69. * bytes by adding up to 7 zero bits on the end. (In fact
  70. * the bitmap passed as input to this function will already
  71. * have had this done in practice.)
  72. *
  73. * + We divide the _byte_ stream exactly in half, rounding the
  74. * half-way position _down_. So an 81-bit input string, for
  75. * example, rounds up to 88 bits or 11 bytes, and then
  76. * dividing by two gives 5 bytes in the first half and 6 in
  77. * the second half.
  78. *
  79. * + We generate a mask from the second half of the bytes, and
  80. * XOR it over the first half.
  81. *
  82. * + We generate a mask from the (encoded) first half of the
  83. * bytes, and XOR it over the second half. Any null bits at
  84. * the end which were added as padding are cleared back to
  85. * zero even if this operation would have made them nonzero.
  86. *
  87. * To de-obfuscate, the steps are precisely the same except
  88. * that the final two are reversed.
  89. *
  90. * Finally, our masking function. Given an input seed string of
  91. * bytes, the output mask consists of concatenating the SHA-1
  92. * hashes of the seed string and successive decimal integers,
  93. * starting from 0.
  94. */
  95. bytes = (bits + 7) / 8;
  96. firsthalf = bytes / 2;
  97. secondhalf = bytes - firsthalf;
  98. steps[decode ? 1 : 0].seedstart = bmp + firsthalf;
  99. steps[decode ? 1 : 0].seedlen = secondhalf;
  100. steps[decode ? 1 : 0].targetstart = bmp;
  101. steps[decode ? 1 : 0].targetlen = firsthalf;
  102. steps[decode ? 0 : 1].seedstart = bmp;
  103. steps[decode ? 0 : 1].seedlen = firsthalf;
  104. steps[decode ? 0 : 1].targetstart = bmp + firsthalf;
  105. steps[decode ? 0 : 1].targetlen = secondhalf;
  106. for (i = 0; i < 2; i++) {
  107. SHA_State base, final;
  108. unsigned char digest[20];
  109. char numberbuf[80];
  110. int digestpos = 20, counter = 0;
  111. SHA_Init(&base);
  112. SHA_Bytes(&base, steps[i].seedstart, steps[i].seedlen);
  113. for (j = 0; j < steps[i].targetlen; j++) {
  114. if (digestpos >= 20) {
  115. sprintf(numberbuf, "%d", counter++);
  116. final = base;
  117. SHA_Bytes(&final, numberbuf, strlen(numberbuf));
  118. SHA_Final(&final, digest);
  119. digestpos = 0;
  120. }
  121. steps[i].targetstart[j] ^= digest[digestpos++];
  122. }
  123. /*
  124. * Mask off the pad bits in the final byte after both steps.
  125. */
  126. if (bits % 8)
  127. bmp[bits / 8] &= 0xFF & (0xFF00 >> (bits % 8));
  128. }
  129. }
  130. /* err, yeah, these two pretty much rely on unsigned char being 8 bits.
  131. * Platforms where this is not the case probably have bigger problems
  132. * than just making these two work, though... */
  133. char *bin2hex(const unsigned char *in, int inlen)
  134. {
  135. char *ret = snewn(inlen*2 + 1, char), *p = ret;
  136. int i;
  137. for (i = 0; i < inlen*2; i++) {
  138. int v = in[i/2];
  139. if (i % 2 == 0) v >>= 4;
  140. *p++ = "0123456789abcdef"[v & 0xF];
  141. }
  142. *p = '\0';
  143. return ret;
  144. }
  145. unsigned char *hex2bin(const char *in, int outlen)
  146. {
  147. unsigned char *ret = snewn(outlen, unsigned char);
  148. int i;
  149. memset(ret, 0, outlen*sizeof(unsigned char));
  150. for (i = 0; i < outlen*2; i++) {
  151. int c = in[i];
  152. int v;
  153. assert(c != 0);
  154. if (c >= '0' && c <= '9')
  155. v = c - '0';
  156. else if (c >= 'a' && c <= 'f')
  157. v = c - 'a' + 10;
  158. else if (c >= 'A' && c <= 'F')
  159. v = c - 'A' + 10;
  160. else
  161. v = 0;
  162. ret[i / 2] |= v << (4 * (1 - (i % 2)));
  163. }
  164. return ret;
  165. }
  166. char *fgetline(FILE *fp)
  167. {
  168. char *ret = snewn(512, char);
  169. int size = 512, len = 0;
  170. while (fgets(ret + len, size - len, fp)) {
  171. len += strlen(ret + len);
  172. if (ret[len-1] == '\n')
  173. break; /* got a newline, we're done */
  174. size = len + 512;
  175. ret = sresize(ret, size, char);
  176. }
  177. if (len == 0) { /* first fgets returned NULL */
  178. sfree(ret);
  179. return NULL;
  180. }
  181. ret[len] = '\0';
  182. return ret;
  183. }
  184. int getenv_bool(const char *name, int dflt)
  185. {
  186. char *env = getenv(name);
  187. if (env == NULL) return dflt;
  188. if (strchr("yYtT", env[0])) return true;
  189. return false;
  190. }
  191. /* Utility functions for colour manipulation. */
  192. static float colour_distance(const float a[3], const float b[3])
  193. {
  194. return (float)sqrt((a[0]-b[0]) * (a[0]-b[0]) +
  195. (a[1]-b[1]) * (a[1]-b[1]) +
  196. (a[2]-b[2]) * (a[2]-b[2]));
  197. }
  198. void colour_mix(const float src1[3], const float src2[3], float p, float dst[3])
  199. {
  200. int i;
  201. for (i = 0; i < 3; i++)
  202. dst[i] = src1[i] * (1.0F - p) + src2[i] * p;
  203. }
  204. void game_mkhighlight_specific(frontend *fe, float *ret,
  205. int background, int highlight, int lowlight)
  206. {
  207. static const float black[3] = { 0.0F, 0.0F, 0.0F };
  208. static const float white[3] = { 1.0F, 1.0F, 1.0F };
  209. float db, dw;
  210. int i;
  211. /*
  212. * New geometric highlight-generation algorithm: Draw a line from
  213. * the base colour to white. The point K distance along this line
  214. * from the base colour is the highlight colour. Similarly, draw
  215. * a line from the base colour to black. The point on this line
  216. * at a distance K from the base colour is the shadow. If either
  217. * of these colours is imaginary (for reasonable K at most one
  218. * will be), _extrapolate_ the base colour along the same line
  219. * until it's a distance K from white (or black) and start again
  220. * with that as the base colour.
  221. *
  222. * This preserves the hue of the base colour, ensures that of the
  223. * three the base colour is the most saturated, and only ever
  224. * flattens the highlight and shadow to pure white or pure black.
  225. *
  226. * K must be at most sqrt(3)/2, or mid grey would be too close to
  227. * both white and black. Here K is set to sqrt(3)/6 so that this
  228. * code produces the same results as the former code in the common
  229. * case where the background is grey and the highlight saturates
  230. * to white.
  231. */
  232. const float k = sqrt(3)/6.0F;
  233. if (lowlight >= 0) {
  234. db = colour_distance(&ret[background*3], black);
  235. if (db < k) {
  236. for (i = 0; i < 3; i++) ret[lowlight*3+i] = black[i];
  237. if (db == 0.0F)
  238. colour_mix(black, white, k/sqrt(3), &ret[background*3]);
  239. else
  240. colour_mix(black, &ret[background*3], k/db, &ret[background*3]);
  241. } else {
  242. colour_mix(&ret[background*3], black, k/db, &ret[lowlight*3]);
  243. }
  244. }
  245. if (highlight >= 0) {
  246. dw = colour_distance(&ret[background*3], white);
  247. if (dw < k) {
  248. for (i = 0; i < 3; i++) ret[highlight*3+i] = white[i];
  249. if (dw == 0.0F)
  250. colour_mix(white, black, k/sqrt(3), &ret[background*3]);
  251. else
  252. colour_mix(white, &ret[background*3], k/dw, &ret[background*3]);
  253. /* Background has changed; recalculate lowlight. */
  254. if (lowlight >= 0)
  255. colour_mix(&ret[background*3], black, k/db, &ret[lowlight*3]);
  256. } else {
  257. colour_mix(&ret[background*3], white, k/dw, &ret[highlight*3]);
  258. }
  259. }
  260. }
  261. void game_mkhighlight(frontend *fe, float *ret,
  262. int background, int highlight, int lowlight)
  263. {
  264. frontend_default_colour(fe, &ret[background * 3]);
  265. game_mkhighlight_specific(fe, ret, background, highlight, lowlight);
  266. }
  267. static void memswap(void *av, void *bv, int size)
  268. {
  269. char tmpbuf[512];
  270. char *a = av, *b = bv;
  271. while (size > 0) {
  272. int thislen = min(size, sizeof(tmpbuf));
  273. memcpy(tmpbuf, a, thislen);
  274. memcpy(a, b, thislen);
  275. memcpy(b, tmpbuf, thislen);
  276. a += thislen;
  277. b += thislen;
  278. size -= thislen;
  279. }
  280. }
  281. void shuffle(void *array, int nelts, int eltsize, random_state *rs)
  282. {
  283. char *carray = (char *)array;
  284. int i;
  285. for (i = nelts; i-- > 1 ;) {
  286. int j = random_upto(rs, i+1);
  287. if (j != i)
  288. memswap(carray + eltsize * i, carray + eltsize * j, eltsize);
  289. }
  290. }
  291. void draw_rect_outline(drawing *dr, int x, int y, int w, int h, int colour)
  292. {
  293. int x0 = x, x1 = x+w-1, y0 = y, y1 = y+h-1;
  294. int coords[8];
  295. coords[0] = x0;
  296. coords[1] = y0;
  297. coords[2] = x0;
  298. coords[3] = y1;
  299. coords[4] = x1;
  300. coords[5] = y1;
  301. coords[6] = x1;
  302. coords[7] = y0;
  303. draw_polygon(dr, coords, 4, -1, colour);
  304. }
  305. void draw_rect_corners(drawing *dr, int cx, int cy, int r, int col)
  306. {
  307. draw_line(dr, cx - r, cy - r, cx - r, cy - r/2, col);
  308. draw_line(dr, cx - r, cy - r, cx - r/2, cy - r, col);
  309. draw_line(dr, cx - r, cy + r, cx - r, cy + r/2, col);
  310. draw_line(dr, cx - r, cy + r, cx - r/2, cy + r, col);
  311. draw_line(dr, cx + r, cy - r, cx + r, cy - r/2, col);
  312. draw_line(dr, cx + r, cy - r, cx + r/2, cy - r, col);
  313. draw_line(dr, cx + r, cy + r, cx + r, cy + r/2, col);
  314. draw_line(dr, cx + r, cy + r, cx + r/2, cy + r, col);
  315. }
  316. void move_cursor(int button, int *x, int *y, int maxw, int maxh, bool wrap)
  317. {
  318. int dx = 0, dy = 0;
  319. switch (button) {
  320. case CURSOR_UP: dy = -1; break;
  321. case CURSOR_DOWN: dy = 1; break;
  322. case CURSOR_RIGHT: dx = 1; break;
  323. case CURSOR_LEFT: dx = -1; break;
  324. default: return;
  325. }
  326. if (wrap) {
  327. *x = (*x + dx + maxw) % maxw;
  328. *y = (*y + dy + maxh) % maxh;
  329. } else {
  330. *x = min(max(*x+dx, 0), maxw - 1);
  331. *y = min(max(*y+dy, 0), maxh - 1);
  332. }
  333. }
  334. /* Used in netslide.c and sixteen.c for cursor movement around edge. */
  335. int c2pos(int w, int h, int cx, int cy)
  336. {
  337. if (cy == -1)
  338. return cx; /* top row, 0 .. w-1 (->) */
  339. else if (cx == w)
  340. return w + cy; /* R col, w .. w+h -1 (v) */
  341. else if (cy == h)
  342. return w + h + (w-cx-1); /* bottom row, w+h .. w+h+w-1 (<-) */
  343. else if (cx == -1)
  344. return w + h + w + (h-cy-1); /* L col, w+h+w .. w+h+w+h-1 (^) */
  345. assert(!"invalid cursor pos!");
  346. return -1; /* not reached */
  347. }
  348. int c2diff(int w, int h, int cx, int cy, int button)
  349. {
  350. int diff = 0;
  351. assert(IS_CURSOR_MOVE(button));
  352. /* Obvious moves around edge. */
  353. if (cy == -1)
  354. diff = (button == CURSOR_RIGHT) ? +1 : (button == CURSOR_LEFT) ? -1 : diff;
  355. if (cy == h)
  356. diff = (button == CURSOR_RIGHT) ? -1 : (button == CURSOR_LEFT) ? +1 : diff;
  357. if (cx == -1)
  358. diff = (button == CURSOR_UP) ? +1 : (button == CURSOR_DOWN) ? -1 : diff;
  359. if (cx == w)
  360. diff = (button == CURSOR_UP) ? -1 : (button == CURSOR_DOWN) ? +1 : diff;
  361. if (button == CURSOR_LEFT && cx == w && (cy == 0 || cy == h-1))
  362. diff = (cy == 0) ? -1 : +1;
  363. if (button == CURSOR_RIGHT && cx == -1 && (cy == 0 || cy == h-1))
  364. diff = (cy == 0) ? +1 : -1;
  365. if (button == CURSOR_DOWN && cy == -1 && (cx == 0 || cx == w-1))
  366. diff = (cx == 0) ? -1 : +1;
  367. if (button == CURSOR_UP && cy == h && (cx == 0 || cx == w-1))
  368. diff = (cx == 0) ? +1 : -1;
  369. debug(("cx,cy = %d,%d; w%d h%d, diff = %d", cx, cy, w, h, diff));
  370. return diff;
  371. }
  372. void pos2c(int w, int h, int pos, int *cx, int *cy)
  373. {
  374. int max = w+h+w+h;
  375. pos = (pos + max) % max;
  376. if (pos < w) {
  377. *cx = pos; *cy = -1; return;
  378. }
  379. pos -= w;
  380. if (pos < h) {
  381. *cx = w; *cy = pos; return;
  382. }
  383. pos -= h;
  384. if (pos < w) {
  385. *cx = w-pos-1; *cy = h; return;
  386. }
  387. pos -= w;
  388. if (pos < h) {
  389. *cx = -1; *cy = h-pos-1; return;
  390. }
  391. assert(!"invalid pos, huh?"); /* limited by % above! */
  392. }
  393. void draw_text_outline(drawing *dr, int x, int y, int fonttype,
  394. int fontsize, int align,
  395. int text_colour, int outline_colour, const char *text)
  396. {
  397. if (outline_colour > -1) {
  398. draw_text(dr, x-1, y, fonttype, fontsize, align, outline_colour, text);
  399. draw_text(dr, x+1, y, fonttype, fontsize, align, outline_colour, text);
  400. draw_text(dr, x, y-1, fonttype, fontsize, align, outline_colour, text);
  401. draw_text(dr, x, y+1, fonttype, fontsize, align, outline_colour, text);
  402. }
  403. draw_text(dr, x, y, fonttype, fontsize, align, text_colour, text);
  404. }
  405. /* kludge for sprintf() in Rockbox not supporting "%-8.8s" */
  406. void copy_left_justified(char *buf, size_t sz, const char *str)
  407. {
  408. size_t len = strlen(str);
  409. assert(sz > 0);
  410. memset(buf, ' ', sz - 1);
  411. assert(len <= sz - 1);
  412. memcpy(buf, str, len);
  413. buf[sz - 1] = 0;
  414. }
  415. /* Returns a dynamically allocated label for a generic button.
  416. * Game-specific buttons should go into the `label' field of key_label
  417. * instead. */
  418. char *button2label(int button)
  419. {
  420. /* check if it's a keyboard button */
  421. if(('A' <= button && button <= 'Z') ||
  422. ('a' <= button && button <= 'z') ||
  423. ('0' <= button && button <= '9') )
  424. {
  425. char str[2];
  426. str[0] = button;
  427. str[1] = '\0';
  428. return dupstr(str);
  429. }
  430. switch(button)
  431. {
  432. case CURSOR_UP:
  433. return dupstr("Up");
  434. case CURSOR_DOWN:
  435. return dupstr("Down");
  436. case CURSOR_LEFT:
  437. return dupstr("Left");
  438. case CURSOR_RIGHT:
  439. return dupstr("Right");
  440. case CURSOR_SELECT:
  441. return dupstr("Select");
  442. case '\b':
  443. return dupstr("Clear");
  444. default:
  445. fatal("unknown generic key");
  446. }
  447. /* should never get here */
  448. return NULL;
  449. }
  450. char *make_prefs_path(const char *dir, const char *sep,
  451. const game *game, const char *suffix)
  452. {
  453. size_t dirlen = strlen(dir);
  454. size_t seplen = strlen(sep);
  455. size_t gamelen = strlen(game->name);
  456. size_t suffixlen = strlen(suffix);
  457. char *path, *p;
  458. const char *q;
  459. if (!dir || !sep || !game || !suffix)
  460. return NULL;
  461. path = snewn(dirlen + seplen + gamelen + suffixlen + 1, char);
  462. p = path;
  463. memcpy(p, dir, dirlen);
  464. p += dirlen;
  465. memcpy(p, sep, seplen);
  466. p += seplen;
  467. for (q = game->name; *q; q++)
  468. if (*q != ' ')
  469. *p++ = tolower((unsigned char)*q);
  470. memcpy(p, suffix, suffixlen);
  471. p += suffixlen;
  472. *p = '\0';
  473. return path;
  474. }
  475. /*
  476. * Calculate the nearest integer to n*sqrt(k), via a bitwise algorithm
  477. * that avoids floating point.
  478. *
  479. * (It would probably be OK in practice to use floating point, but I
  480. * felt like overengineering it for fun. With FP, there's at least a
  481. * theoretical risk of rounding the wrong way, due to the three
  482. * successive roundings involved - rounding sqrt(k), rounding its
  483. * product with n, and then rounding to the nearest integer. This
  484. * approach avoids that: it's exact.)
  485. */
  486. int n_times_root_k(int n_signed, int k)
  487. {
  488. unsigned x, r, m;
  489. int sign = n_signed < 0 ? -1 : +1;
  490. unsigned n = n_signed * sign;
  491. unsigned bitpos;
  492. /*
  493. * Method:
  494. *
  495. * We transform m gradually from zero into n, by multiplying it by
  496. * 2 in each step and optionally adding 1, so that it's always
  497. * floor(n/2^something).
  498. *
  499. * At the start of each step, x is the largest integer less than
  500. * or equal to m*sqrt(k). We transform m to 2m+bit, and therefore
  501. * we must transform x to 2x+something to match. The 'something'
  502. * we add to 2x is at most floor(sqrt(k))+2. (Worst case is if m
  503. * sqrt(k) was equal to x + 1-eps for some tiny eps, and then the
  504. * incoming bit of m is 1, so that (2m+1)sqrt(k) =
  505. * 2x+2+sqrt(k)-2eps.)
  506. *
  507. * To compute this, we also track the residual value r such that
  508. * x^2+r = km^2.
  509. *
  510. * The algorithm below is very similar to the usual approach for
  511. * taking the square root of an integer in binary. The wrinkle is
  512. * that we have an integer multiplier, i.e. we're computing
  513. * n*sqrt(k) rather than just sqrt(k). Of course in principle we
  514. * could just take sqrt(n^2k), but we'd need an integer twice the
  515. * width to hold n^2. Pulling out n and treating it specially
  516. * makes overflow less likely.
  517. */
  518. x = r = m = 0;
  519. for (bitpos = UINT_MAX & ~(UINT_MAX >> 1); bitpos; bitpos >>= 1) {
  520. unsigned a, b = (n & bitpos) ? 1 : 0;
  521. /*
  522. * Check invariants. We expect that x^2 + r = km^2 (i.e. our
  523. * residual term is correct), and also that r < 2x+1 (because
  524. * if not, then we could replace x with x+1 and still get a
  525. * value that made r non-negative, i.e. x would not be the
  526. * _largest_ integer less than m sqrt(k)).
  527. */
  528. assert(x*x + r == k*m*m);
  529. assert(r < 2*x+1);
  530. /*
  531. * We're going to replace m with 2m+b, and x with 2x+a for
  532. * some a we haven't decided on yet.
  533. *
  534. * The new value of the residual will therefore be
  535. *
  536. * k (2m+b)^2 - (2x+a)^2
  537. * = (4km^2 + 4kmb + kb^2) - (4x^2 + 4xa + a^2)
  538. * = 4 (km^2 - x^2) + 4kmb + kb^2 - 4xa - a^2
  539. * = 4r + 4kmb + kb^2 - 4xa - a^2 (because r = km^2 - x^2)
  540. * = 4r + (4m + 1)kb - 4xa - a^2 (b is 0 or 1, so b = b^2)
  541. */
  542. for (a = 0;; a++) {
  543. /* If we made this routine handle square roots of numbers
  544. * significantly bigger than 3 or 5 then it would be
  545. * sensible to make this a binary search. Here, it hardly
  546. * seems important. */
  547. unsigned pos = 4*r + k*b*(4*m + 1);
  548. unsigned neg = 4*a*x + a*a;
  549. if (pos < neg)
  550. break; /* this value of a is too big */
  551. }
  552. /* The above loop will have terminated with a one too big. So
  553. * now decrementing a will give us the right value to add. */
  554. a--;
  555. r = 4*r + b*k*(4*m + 1) - (4*a*x + a*a);
  556. m = 2*m+b;
  557. x = 2*x+a;
  558. }
  559. /*
  560. * Finally, round to the nearest integer. At present, x is the
  561. * largest integer that is _at most_ m sqrt(k). But we want the
  562. * _nearest_ integer, whether that's rounded up or down. So check
  563. * whether (x + 1/2) is still less than m sqrt(k), i.e. whether
  564. * (x + 1/2)^2 < km^2; if it is, then we increment x.
  565. *
  566. * We have km^2 - (x + 1/2)^2 = km^2 - x^2 - x - 1/4
  567. * = r - x - 1/4
  568. *
  569. * and since r and x are integers, this is greater than 0 if and
  570. * only if r > x.
  571. *
  572. * (There's no need to worry about tie-breaking exact halfway
  573. * rounding cases. sqrt(k) is irrational, so none such exist.)
  574. */
  575. if (r > x)
  576. x++;
  577. /*
  578. * Put the sign back on, and convert back from unsigned to int.
  579. */
  580. if (sign == +1) {
  581. return x;
  582. } else {
  583. /* Be a little careful to avoid compilers deciding I've just
  584. * perpetrated signed-integer overflow. This should optimise
  585. * down to no actual code. */
  586. return INT_MIN + (int)(-x - (unsigned)INT_MIN);
  587. }
  588. }
  589. /* vim: set shiftwidth=4 tabstop=8: */