radix.odin 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480
  1. package math_big
  2. /*
  3. Copyright 2021 Jeroen van Rijn <[email protected]>.
  4. Made available under Odin's BSD-3 license.
  5. An arbitrary precision mathematics implementation in Odin.
  6. For the theoretical underpinnings, see Knuth's The Art of Computer Programming, Volume 2, section 4.3.
  7. The code started out as an idiomatic source port of libTomMath, which is in the public domain, with thanks.
  8. This file contains radix conversions, `string_to_int` (atoi) and `int_to_string` (itoa).
  9. TODO:
  10. - Use Barrett reduction for non-powers-of-two.
  11. - Also look at extracting and splatting several digits at once.
  12. */
  13. import "core:intrinsics"
  14. import "core:mem"
  15. /*
  16. This version of `itoa` allocates one behalf of the caller. The caller must free the string.
  17. */
  18. int_itoa_string :: proc(a: ^Int, radix := i8(-1), zero_terminate := false, allocator := context.allocator) -> (res: string, err: Error) {
  19. assert_if_nil(a);
  20. context.allocator = allocator;
  21. a := a; radix := radix;
  22. clear_if_uninitialized(a) or_return;
  23. /*
  24. Radix defaults to 10.
  25. */
  26. radix = radix if radix > 0 else 10;
  27. /*
  28. TODO: If we want to write a prefix for some of the radixes, we can oversize the buffer.
  29. Then after the digits are written and the string is reversed
  30. */
  31. /*
  32. Calculate the size of the buffer we need, and
  33. Exit if calculating the size returned an error.
  34. */
  35. size := radix_size(a, radix, zero_terminate) or_return;
  36. /*
  37. Allocate the buffer we need.
  38. */
  39. buffer := make([]u8, size);
  40. /*
  41. Write the digits out into the buffer.
  42. */
  43. written: int;
  44. written, err = int_itoa_raw(a, radix, buffer, size, zero_terminate);
  45. return string(buffer[:written]), err;
  46. }
  47. /*
  48. This version of `itoa` allocates one behalf of the caller. The caller must free the string.
  49. */
  50. int_itoa_cstring :: proc(a: ^Int, radix := i8(-1), allocator := context.allocator) -> (res: cstring, err: Error) {
  51. assert_if_nil(a);
  52. context.allocator = allocator;
  53. a := a; radix := radix;
  54. clear_if_uninitialized(a) or_return;
  55. /*
  56. Radix defaults to 10.
  57. */
  58. radix = radix if radix > 0 else 10;
  59. s: string;
  60. s, err = int_itoa_string(a, radix, true);
  61. return cstring(raw_data(s)), err;
  62. }
  63. /*
  64. A low-level `itoa` using a caller-provided buffer. `itoa_string` and `itoa_cstring` use this.
  65. You can use also use it if you want to pre-allocate a buffer and optionally reuse it.
  66. Use `radix_size` or `radix_size_estimate` to determine a buffer size big enough.
  67. You can pass the output of `radix_size` to `size` if you've previously called it to size
  68. the output buffer. If you haven't, this routine will call it. This way it knows if the buffer
  69. is the appropriate size, and we can write directly in place without a reverse step at the end.
  70. === === === IMPORTANT === === ===
  71. If you determined the buffer size using `radix_size_estimate`, or have a buffer
  72. that you reuse that you know is large enough, don't pass this size unless you know what you are doing,
  73. because we will always write backwards starting at last byte of the buffer.
  74. Keep in mind that if you set `size` yourself and it's smaller than the buffer,
  75. it'll result in buffer overflows, as we use it to avoid reversing at the end
  76. and having to perform a buffer overflow check each character.
  77. */
  78. int_itoa_raw :: proc(a: ^Int, radix: i8, buffer: []u8, size := int(-1), zero_terminate := false) -> (written: int, err: Error) {
  79. assert_if_nil(a);
  80. a := a; radix := radix; size := size;
  81. clear_if_uninitialized(a) or_return;
  82. /*
  83. Radix defaults to 10.
  84. */
  85. radix = radix if radix > 0 else 10;
  86. if radix < 2 || radix > 64 {
  87. return 0, .Invalid_Argument;
  88. }
  89. /*
  90. We weren't given a size. Let's compute it.
  91. */
  92. if size == -1 {
  93. size = radix_size(a, radix, zero_terminate) or_return;
  94. }
  95. /*
  96. Early exit if the buffer we were given is too small.
  97. */
  98. available := len(buffer);
  99. if available < size {
  100. return 0, .Buffer_Overflow;
  101. }
  102. /*
  103. Fast path for when `Int` == 0 or the entire `Int` fits in a single radix digit.
  104. */
  105. z, _ := is_zero(a);
  106. if z || (a.used == 1 && a.digit[0] < DIGIT(radix)) {
  107. if zero_terminate {
  108. available -= 1;
  109. buffer[available] = 0;
  110. }
  111. available -= 1;
  112. buffer[available] = RADIX_TABLE[a.digit[0]];
  113. if n, _ := is_neg(a); n {
  114. available -= 1;
  115. buffer[available] = '-';
  116. }
  117. /*
  118. If we overestimated the size, we need to move the buffer left.
  119. */
  120. written = len(buffer) - available;
  121. if written < size {
  122. diff := size - written;
  123. mem.copy(&buffer[0], &buffer[diff], written);
  124. }
  125. return written, nil;
  126. }
  127. /*
  128. Fast path for when `Int` fits within a `_WORD`.
  129. */
  130. if a.used == 1 || a.used == 2 {
  131. if zero_terminate {
  132. available -= 1;
  133. buffer[available] = 0;
  134. }
  135. val := _WORD(a.digit[1]) << _DIGIT_BITS + _WORD(a.digit[0]);
  136. for val > 0 {
  137. q := val / _WORD(radix);
  138. available -= 1;
  139. buffer[available] = RADIX_TABLE[val - (q * _WORD(radix))];
  140. val = q;
  141. }
  142. if n, _ := is_neg(a); n {
  143. available -= 1;
  144. buffer[available] = '-';
  145. }
  146. /*
  147. If we overestimated the size, we need to move the buffer left.
  148. */
  149. written = len(buffer) - available;
  150. if written < size {
  151. diff := size - written;
  152. mem.copy(&buffer[0], &buffer[diff], written);
  153. }
  154. return written, nil;
  155. }
  156. /*
  157. Fast path for radixes that are a power of two.
  158. */
  159. if is_power_of_two(int(radix)) {
  160. if zero_terminate {
  161. available -= 1;
  162. buffer[available] = 0;
  163. }
  164. shift, count: int;
  165. // mask := _WORD(radix - 1);
  166. shift, err = log(DIGIT(radix), 2);
  167. count, err = count_bits(a);
  168. digit: _WORD;
  169. for offset := 0; offset < count; offset += shift {
  170. bits_to_get := int(min(count - offset, shift));
  171. digit, err = int_bitfield_extract(a, offset, bits_to_get);
  172. if err != nil {
  173. return len(buffer) - available, .Invalid_Argument;
  174. }
  175. available -= 1;
  176. buffer[available] = RADIX_TABLE[digit];
  177. }
  178. if n, _ := is_neg(a); n {
  179. available -= 1;
  180. buffer[available] = '-';
  181. }
  182. /*
  183. If we overestimated the size, we need to move the buffer left.
  184. */
  185. written = len(buffer) - available;
  186. if written < size {
  187. diff := size - written;
  188. mem.copy(&buffer[0], &buffer[diff], written);
  189. }
  190. return written, nil;
  191. }
  192. return _itoa_raw_full(a, radix, buffer, zero_terminate);
  193. }
  194. itoa :: proc{int_itoa_string, int_itoa_raw};
  195. int_to_string :: int_itoa_string;
  196. int_to_cstring :: int_itoa_cstring;
  197. /*
  198. Read a string [ASCII] in a given radix.
  199. */
  200. int_atoi :: proc(res: ^Int, input: string, radix := i8(10), allocator := context.allocator) -> (err: Error) {
  201. assert_if_nil(res);
  202. input := input;
  203. context.allocator = allocator;
  204. /*
  205. Make sure the radix is ok.
  206. */
  207. if radix < 2 || radix > 64 { return .Invalid_Argument; }
  208. /*
  209. Set the integer to the default of zero.
  210. */
  211. internal_zero(res) or_return;
  212. /*
  213. We'll interpret an empty string as zero.
  214. */
  215. if len(input) == 0 {
  216. return nil;
  217. }
  218. /*
  219. If the leading digit is a minus set the sign to negative.
  220. Given the above early out, the length should be at least 1.
  221. */
  222. sign := Sign.Zero_or_Positive;
  223. if input[0] == '-' {
  224. input = input[1:];
  225. sign = .Negative;
  226. }
  227. /*
  228. Process each digit of the string.
  229. */
  230. ch: rune;
  231. for len(input) > 0 {
  232. /* if the radix <= 36 the conversion is case insensitive
  233. * this allows numbers like 1AB and 1ab to represent the same value
  234. * [e.g. in hex]
  235. */
  236. ch = rune(input[0]);
  237. if radix <= 36 && ch >= 'a' && ch <= 'z' {
  238. ch -= 32; // 'a' - 'A'
  239. }
  240. pos := ch - '+';
  241. if RADIX_TABLE_REVERSE_SIZE <= pos {
  242. break;
  243. }
  244. y := RADIX_TABLE_REVERSE[pos];
  245. /* if the char was found in the map
  246. * and is less than the given radix add it
  247. * to the number, otherwise exit the loop.
  248. */
  249. if y >= u8(radix) {
  250. break;
  251. }
  252. internal_mul(res, res, DIGIT(radix)) or_return;
  253. internal_add(res, res, DIGIT(y)) or_return;
  254. input = input[1:];
  255. }
  256. /*
  257. If an illegal character was found, fail.
  258. */
  259. if len(input) > 0 && ch != 0 && ch != '\r' && ch != '\n' {
  260. return .Invalid_Argument;
  261. }
  262. /*
  263. Set the sign only if res != 0.
  264. */
  265. if res.used > 0 {
  266. res.sign = sign;
  267. }
  268. return nil;
  269. }
  270. atoi :: proc { int_atoi, };
  271. /*
  272. We size for `string` by default.
  273. */
  274. radix_size :: proc(a: ^Int, radix: i8, zero_terminate := false, allocator := context.allocator) -> (size: int, err: Error) {
  275. a := a;
  276. assert_if_nil(a);
  277. if radix < 2 || radix > 64 { return -1, .Invalid_Argument; }
  278. clear_if_uninitialized(a) or_return;
  279. if internal_is_zero(a) {
  280. if zero_terminate {
  281. return 2, nil;
  282. }
  283. return 1, nil;
  284. }
  285. if internal_is_power_of_two(a) {
  286. /*
  287. Calculate `log` on a temporary "copy" with its sign set to positive.
  288. */
  289. t := &Int{
  290. used = a.used,
  291. sign = .Zero_or_Positive,
  292. digit = a.digit,
  293. };
  294. size = internal_log(t, DIGIT(radix)) or_return;
  295. } else {
  296. la, k := &Int{}, &Int{};
  297. defer internal_destroy(la, k);
  298. /* la = floor(log_2(a)) + 1 */
  299. bit_count := internal_count_bits(a);
  300. internal_set(la, bit_count) or_return;
  301. /* k = floor(2^29/log_2(radix)) + 1 */
  302. lb := _log_bases;
  303. internal_set(k, lb[radix]) or_return;
  304. /* n = floor((la * k) / 2^29) + 1 */
  305. internal_mul(k, la, k) or_return;
  306. internal_shr(k, k, _RADIX_SIZE_SCALE) or_return;
  307. /* The "+1" here is the "+1" in "floor((la * k) / 2^29) + 1" */
  308. /* n = n + 1 + EOS + sign */
  309. size_, _ := internal_get(k, u128);
  310. size = int(size_);
  311. }
  312. /*
  313. log truncates to zero, so we need to add one more, and one for `-` if negative.
  314. */
  315. size += 2 if a.sign == .Negative else 1;
  316. size += 1 if zero_terminate else 0;
  317. return size, nil;
  318. }
  319. /*
  320. Overestimate the size needed for the bigint to string conversion by a very small amount.
  321. The error is about 10^-8; it will overestimate the result by at most 11 elements for
  322. a number of the size 2^(2^31)-1 which is currently the largest possible in this library.
  323. Some short tests gave no results larger than 5 (plus 2 for sign and EOS).
  324. */
  325. /*
  326. Table of {0, INT(log_2([1..64])*2^p)+1 } where p is the scale
  327. factor defined in MP_RADIX_SIZE_SCALE and INT() extracts the integer part (truncating).
  328. Good for 32 bit "int". Set MP_RADIX_SIZE_SCALE = 61 and recompute values
  329. for 64 bit "int".
  330. */
  331. _RADIX_SIZE_SCALE :: 29;
  332. _log_bases :: [65]u32{
  333. 0, 0, 0x20000001, 0x14309399, 0x10000001,
  334. 0xdc81a35, 0xc611924, 0xb660c9e, 0xaaaaaab, 0xa1849cd,
  335. 0x9a209a9, 0x94004e1, 0x8ed19c2, 0x8a5ca7d, 0x867a000,
  336. 0x830cee3, 0x8000001, 0x7d42d60, 0x7ac8b32, 0x7887847,
  337. 0x7677349, 0x749131f, 0x72d0163, 0x712f657, 0x6fab5db,
  338. 0x6e40d1b, 0x6ced0d0, 0x6badbde, 0x6a80e3b, 0x6964c19,
  339. 0x6857d31, 0x6758c38, 0x6666667, 0x657fb21, 0x64a3b9f,
  340. 0x63d1ab4, 0x6308c92, 0x624869e, 0x618ff47, 0x60dedea,
  341. 0x6034ab0, 0x5f90e7b, 0x5ef32cb, 0x5e5b1b2, 0x5dc85c3,
  342. 0x5d3aa02, 0x5cb19d9, 0x5c2d10f, 0x5bacbbf, 0x5b3064f,
  343. 0x5ab7d68, 0x5a42df0, 0x59d1506, 0x5962ffe, 0x58f7c57,
  344. 0x588f7bc, 0x582a000, 0x57c7319, 0x5766f1d, 0x5709243,
  345. 0x56adad9, 0x565474d, 0x55fd61f, 0x55a85e8, 0x5555556,
  346. };
  347. /*
  348. Characters used in radix conversions.
  349. */
  350. RADIX_TABLE := "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
  351. RADIX_TABLE_REVERSE := [RADIX_TABLE_REVERSE_SIZE]u8{
  352. 0x3e, 0xff, 0xff, 0xff, 0x3f, 0x00, 0x01, 0x02, 0x03, 0x04, /* +,-./01234 */
  353. 0x05, 0x06, 0x07, 0x08, 0x09, 0xff, 0xff, 0xff, 0xff, 0xff, /* 56789:;<=> */
  354. 0xff, 0xff, 0x0a, 0x0b, 0x0c, 0x0d, 0x0e, 0x0f, 0x10, 0x11, /* ?@ABCDEFGH */
  355. 0x12, 0x13, 0x14, 0x15, 0x16, 0x17, 0x18, 0x19, 0x1a, 0x1b, /* IJKLMNOPQR */
  356. 0x1c, 0x1d, 0x1e, 0x1f, 0x20, 0x21, 0x22, 0x23, 0xff, 0xff, /* STUVWXYZ[\ */
  357. 0xff, 0xff, 0xff, 0xff, 0x24, 0x25, 0x26, 0x27, 0x28, 0x29, /* ]^_`abcdef */
  358. 0x2a, 0x2b, 0x2c, 0x2d, 0x2e, 0x2f, 0x30, 0x31, 0x32, 0x33, /* ghijklmnop */
  359. 0x34, 0x35, 0x36, 0x37, 0x38, 0x39, 0x3a, 0x3b, 0x3c, 0x3d, /* qrstuvwxyz */
  360. };
  361. RADIX_TABLE_REVERSE_SIZE :: 80;
  362. /*
  363. Stores a bignum as a ASCII string in a given radix (2..64)
  364. The buffer must be appropriately sized. This routine doesn't check.
  365. */
  366. _itoa_raw_full :: proc(a: ^Int, radix: i8, buffer: []u8, zero_terminate := false, allocator := context.allocator) -> (written: int, err: Error) {
  367. assert_if_nil(a);
  368. context.allocator = allocator;
  369. temp, denominator := &Int{}, &Int{};
  370. internal_copy(temp, a) or_return;
  371. internal_set(denominator, radix) or_return;
  372. available := len(buffer);
  373. if zero_terminate {
  374. available -= 1;
  375. buffer[available] = 0;
  376. }
  377. if a.sign == .Negative {
  378. temp.sign = .Zero_or_Positive;
  379. }
  380. remainder: DIGIT;
  381. for {
  382. if remainder, err = #force_inline internal_divmod(temp, temp, DIGIT(radix)); err != nil {
  383. internal_destroy(temp, denominator);
  384. return len(buffer) - available, err;
  385. }
  386. available -= 1;
  387. buffer[available] = RADIX_TABLE[remainder];
  388. if temp.used == 0 {
  389. break;
  390. }
  391. }
  392. if a.sign == .Negative {
  393. available -= 1;
  394. buffer[available] = '-';
  395. }
  396. internal_destroy(temp, denominator);
  397. /*
  398. If we overestimated the size, we need to move the buffer left.
  399. */
  400. written = len(buffer) - available;
  401. if written < len(buffer) {
  402. diff := len(buffer) - written;
  403. mem.copy(&buffer[0], &buffer[diff], written);
  404. }
  405. return written, nil;
  406. }