bigint.lua 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512
  1. #!/usr/bin/env lua
  2. -- If this variable is true, then strict type checking is performed for all
  3. -- operations. This may result in slower code, but it will allow you to catch
  4. -- errors and bugs earlier.
  5. local strict = true
  6. --------------------------------------------------------------------------------
  7. local bigint = {}
  8. local magnitudes = require("ordersofmagnitude")
  9. -- Create a new bigint or convert a number or string into a big
  10. -- Returns an empty, positive bigint if no number or string is given
  11. function bigint.new(num)
  12. local self = {
  13. sign = "+",
  14. digits = {}
  15. }
  16. -- Return a new bigint with the same sign and digits
  17. function self:clone()
  18. local newint = bigint.new()
  19. newint.sign = self.sign
  20. for _, digit in pairs(self.digits) do
  21. newint.digits[#newint.digits + 1] = digit
  22. end
  23. return newint
  24. end
  25. if (num) then
  26. local num_string = tostring(num)
  27. for digit in string.gmatch(num_string, "[0-9]") do
  28. table.insert(self.digits, tonumber(digit))
  29. end
  30. if string.sub(num_string, 1, 1) == "-" then
  31. self.sign = "-"
  32. end
  33. end
  34. return self
  35. end
  36. -- Check the type of a big
  37. -- Normally only runs when global variable "strict" == true, but checking can be
  38. -- forced by supplying "true" as the second argument.
  39. function bigint.check(big, force)
  40. if (strict or force) then
  41. assert(#big.digits > 0, "bigint is empty")
  42. assert(type(big.sign) == "string", "bigint is unsigned")
  43. for _, digit in pairs(big.digits) do
  44. assert(type(digit) == "number", digit .. " is not a number")
  45. assert(digit < 10, digit .. " is greater than or equal to 10")
  46. end
  47. end
  48. return true
  49. end
  50. -- Return a new big with the same digits but with a positive sign (absolute
  51. -- value)
  52. function bigint.abs(big)
  53. bigint.check(big)
  54. local result = big:clone()
  55. result.sign = "+"
  56. return result
  57. end
  58. -- Convert a big to a number or string
  59. function bigint.unserialize(big, output_type, precision)
  60. bigint.check(big)
  61. local num = ""
  62. if big.sign == "-" then
  63. num = "-"
  64. end
  65. if ((output_type == nil)
  66. or (output_type == "number")
  67. or (output_type == "n")
  68. or (output_type == "string")
  69. or (output_type == "s")) then
  70. -- Unserialization to a string or number requires reconstructing the
  71. -- entire number
  72. for _, digit in pairs(big.digits) do
  73. num = num .. math.floor(digit) -- lazy way of getting rid of .0$
  74. end
  75. if ((output_type == nil)
  76. or (output_type == "number")
  77. or (output_type == "n")) then
  78. return tonumber(num)
  79. else
  80. return num
  81. end
  82. else
  83. -- Unserialization to human-readable form or scientific notation only
  84. -- requires reading the first few digits
  85. if (precision == nil) then
  86. precision = 3
  87. else
  88. assert(precision > 0, "Precision cannot be less than 1")
  89. assert(math.floor(precision) == precision,
  90. "Precision must be a positive integer")
  91. end
  92. num = num .. big.digits[1]
  93. if (precision > 1) then
  94. num = num .. "."
  95. for i = 1, (precision - 1) do
  96. num = num .. big.digits[i + 1]
  97. end
  98. end
  99. if ((output_type == "human-readable")
  100. or (output_type == "human")
  101. or (output_type == "h")) then
  102. local original_number = bigint.unserialize(big,"n") or 0 -- Calling function from within itself lool spooky
  103. local mag = ""
  104. local numb = 0
  105. for i = #big.digits-1 ,#big.digits-6,-1 do
  106. v = magnitudes[i]
  107. if v then
  108. numb = original_number/(10^(i))
  109. mag = magnitudes[i]
  110. break
  111. end
  112. end
  113. local formated_num = string.format("%0.3f",tostring(numb)).." "..mag
  114. return formated_num
  115. else
  116. return num .. "*10^" .. (#big.digits - 1)
  117. end
  118. end
  119. end
  120. -- Basic comparisons
  121. -- Accepts symbols (<, >=, ~=) and Unix shell-like options (lt, ge, ne)
  122. function bigint.compare(big1, big2, comparison)
  123. bigint.check(big1)
  124. bigint.check(big2)
  125. local greater = false -- If big1.digits > big2.digits
  126. local equal = false
  127. if (big1.sign == "-") and (big2.sign == "+") then
  128. greater = false
  129. elseif (#big1.digits > #big2.digits)
  130. or ((big1.sign == "+") and (big2.sign == "-")) then
  131. greater = true
  132. elseif (#big1.digits == #big2.digits) then
  133. -- Walk left to right, comparing digits
  134. for digit = 1, #big1.digits do
  135. if (big1.digits[digit] > big2.digits[digit]) then
  136. greater = true
  137. break
  138. elseif (big2.digits[digit] > big1.digits[digit]) then
  139. break
  140. elseif (digit == #big1.digits)
  141. and (big1.digits[digit] == big2.digits[digit]) then
  142. equal = true
  143. end
  144. end
  145. end
  146. -- If both numbers are negative, then the requirements for greater are
  147. -- reversed
  148. if (not equal) and (big1.sign == "-") and (big2.sign == "-") then
  149. greater = not greater
  150. end
  151. return (((comparison == "<") or (comparison == "lt"))
  152. and ((not greater) and (not equal)) and true)
  153. or (((comparison == ">") or (comparison == "gt"))
  154. and ((greater) and (not equal)) and true)
  155. or (((comparison == "==") or (comparison == "eq"))
  156. and (equal) and true)
  157. or (((comparison == ">=") or (comparison == "ge"))
  158. and (equal or greater) and true)
  159. or (((comparison == "<=") or (comparison == "le"))
  160. and (equal or not greater) and true)
  161. or (((comparison == "~=") or (comparison == "!=") or (comparison == "ne"))
  162. and (not equal) and true)
  163. or false
  164. end
  165. -- BACKEND: Add big1 and big2, ignoring signs
  166. function bigint.add_raw(big1, big2)
  167. bigint.check(big1)
  168. bigint.check(big2)
  169. local result = bigint.new()
  170. local max_digits = 0
  171. local carry = 0
  172. if (#big1.digits >= #big2.digits) then
  173. max_digits = #big1.digits
  174. else
  175. max_digits = #big2.digits
  176. end
  177. -- Walk backwards right to left, like in long addition
  178. for digit = 0, max_digits - 1 do
  179. local sum = (big1.digits[#big1.digits - digit] or 0)
  180. + (big2.digits[#big2.digits - digit] or 0)
  181. + carry
  182. if (sum >= 10) then
  183. carry = 1
  184. sum = sum - 10
  185. else
  186. carry = 0
  187. end
  188. result.digits[max_digits - digit] = sum
  189. end
  190. -- Leftover carry in cases when #big1.digits == #big2.digits and sum > 10, ex. 7 + 9
  191. if (carry == 1) then
  192. table.insert(result.digits, 1, 1)
  193. end
  194. return result
  195. end
  196. -- BACKEND: Subtract big2 from big1, ignoring signs
  197. function bigint.subtract_raw(big1, big2)
  198. -- Type checking is done by bigint.compare
  199. assert(bigint.compare(bigint.abs(big1), bigint.abs(big2), ">="),
  200. "Size of " .. bigint.unserialize(big1, "string") .. " is less than "
  201. .. bigint.unserialize(big2, "string"))
  202. local result = big1:clone()
  203. local max_digits = #big1.digits
  204. local borrow = 0
  205. -- Logic mostly copied from bigint.add_raw ---------------------------------
  206. -- Walk backwards right to left, like in long subtraction
  207. for digit = 0, max_digits - 1 do
  208. local diff = (big1.digits[#big1.digits - digit] or 0)
  209. - (big2.digits[#big2.digits - digit] or 0)
  210. - borrow
  211. if (diff < 0) then
  212. borrow = 1
  213. diff = diff + 10
  214. else
  215. borrow = 0
  216. end
  217. result.digits[max_digits - digit] = diff
  218. end
  219. ----------------------------------------------------------------------------
  220. -- Strip leading zeroes if any, but not if 0 is the only digit
  221. while (#result.digits > 1) and (result.digits[1] == 0) do
  222. table.remove(result.digits, 1)
  223. end
  224. return result
  225. end
  226. -- FRONTEND: Addition and subtraction operations, accounting for signs
  227. function bigint.add(big1, big2)
  228. -- Type checking is done by bigint.compare
  229. local result
  230. -- If adding numbers of different sign, subtract the smaller sized one from
  231. -- the bigger sized one and take the sign of the bigger sized one
  232. if (big1.sign ~= big2.sign) then
  233. if (bigint.compare(bigint.abs(big1), bigint.abs(big2), ">")) then
  234. result = bigint.subtract_raw(big1, big2)
  235. result.sign = big1.sign
  236. else
  237. result = bigint.subtract_raw(big2, big1)
  238. result.sign = big2.sign
  239. end
  240. elseif (big1.sign == "+") and (big2.sign == "+") then
  241. result = bigint.add_raw(big1, big2)
  242. elseif (big1.sign == "-") and (big2.sign == "-") then
  243. result = bigint.add_raw(big1, big2)
  244. result.sign = "-"
  245. end
  246. return result
  247. end
  248. function bigint.subtract(big1, big2)
  249. -- Type checking is done by bigint.compare in bigint.add
  250. -- Subtracting is like adding a negative
  251. local big2_local = big2:clone()
  252. if (big2.sign == "+") then
  253. big2_local.sign = "-"
  254. else
  255. big2_local.sign = "+"
  256. end
  257. return bigint.add(big1, big2_local)
  258. end
  259. -- BACKEND: Multiply a big by a single digit big, ignoring signs
  260. function bigint.multiply_single(big1, big2)
  261. bigint.check(big1)
  262. bigint.check(big2)
  263. assert(#big2.digits == 1, bigint.unserialize(big2, "string")
  264. .. " has more than one digit")
  265. local result = bigint.new()
  266. local carry = 0
  267. -- Logic mostly copied from bigint.add_raw ---------------------------------
  268. -- Walk backwards right to left, like in long multiplication
  269. for digit = 0, #big1.digits - 1 do
  270. local this_digit = big1.digits[#big1.digits - digit]
  271. * big2.digits[1]
  272. + carry
  273. if (this_digit >= 10) then
  274. carry = math.floor(this_digit / 10)
  275. this_digit = this_digit - (carry * 10)
  276. else
  277. carry = 0
  278. end
  279. result.digits[#big1.digits - digit] = this_digit
  280. end
  281. -- Leftover carry in cases when big1.digits[1] * big2.digits[1] > 0
  282. if (carry > 0) then
  283. table.insert(result.digits, 1, carry)
  284. end
  285. ----------------------------------------------------------------------------
  286. return result
  287. end
  288. -- FRONTEND: Multiply two bigs, accounting for signs
  289. function bigint.multiply(big1, big2)
  290. -- Type checking done by bigint.multiply_single
  291. local result = bigint.new(0)
  292. local larger, smaller -- Larger and smaller in terms of digits, not size
  293. if (bigint.unserialize(big1) == 0) or (bigint.unserialize(big2) == 0) then
  294. return result
  295. end
  296. if (#big1.digits >= #big2.digits) then
  297. larger = big1
  298. smaller = big2
  299. else
  300. larger = big2
  301. smaller = big1
  302. end
  303. -- Walk backwards right to left, like in long multiplication
  304. for digit = 0, #smaller.digits - 1 do
  305. -- Sorry for going over column 80! There's lots of big names here
  306. local this_digit_product = bigint.multiply_single(larger,
  307. bigint.new(smaller.digits[#smaller.digits - digit]))
  308. -- "Placeholding zeroes"
  309. if (digit > 0) then
  310. for placeholder = 1, digit do
  311. table.insert(this_digit_product.digits, 0)
  312. end
  313. end
  314. result = bigint.add(result, this_digit_product)
  315. end
  316. if (larger.sign == smaller.sign) then
  317. result.sign = "+"
  318. else
  319. result.sign = "-"
  320. end
  321. return result
  322. end
  323. -- Raise a big to a positive integer or big power (TODO: negative integer power)
  324. function bigint.exponentiate(big, power)
  325. -- Type checking for big done by bigint.multiply
  326. assert(bigint.compare(power, bigint.new(0), ">"),
  327. " negative powers are not supported")
  328. local exp = power:clone()
  329. if (bigint.compare(exp, bigint.new(0), "==")) then
  330. return bigint.new(1)
  331. elseif (bigint.compare(exp, bigint.new(1), "==")) then
  332. return big
  333. else
  334. local result = big:clone()
  335. while (bigint.compare(exp, bigint.new(1), ">")) do
  336. result = bigint.multiply(result, big)
  337. exp = bigint.subtract(exp, bigint.new(1))
  338. end
  339. return result
  340. end
  341. end
  342. -- BACKEND: Divide two bigs (decimals not supported), returning big result and
  343. -- big remainder
  344. -- WARNING: Only supports positive integers
  345. function bigint.divide_raw(big1, big2)
  346. -- Type checking done by bigint.compare
  347. if (bigint.compare(big1, big2, "==")) then
  348. return bigint.new(1), bigint.new(0)
  349. elseif (bigint.compare(big1, big2, "<")) then
  350. return bigint.new(0), bigint.new(0)
  351. else
  352. assert(bigint.compare(big2, bigint.new(0), "!="), "error: divide by zero")
  353. assert(big1.sign == "+", "error: big1 is not positive")
  354. assert(big2.sign == "+", "error: big2 is not positive")
  355. local result = bigint.new()
  356. local dividend = bigint.new() -- Dividend of a single operation, not the
  357. -- dividend of the overall function
  358. local divisor = big2:clone()
  359. local factor = 1
  360. -- Walk left to right among digits in the dividend, like in long
  361. -- division
  362. for _, digit in pairs(big1.digits) do
  363. dividend.digits[#dividend.digits + 1] = digit
  364. -- The dividend is smaller than the divisor, so a zero is appended
  365. -- to the result and the loop ends
  366. if (bigint.compare(dividend, divisor, "<")) then
  367. if (#result.digits > 0) then -- Don't add leading zeroes
  368. result.digits[#result.digits + 1] = 0
  369. end
  370. else
  371. -- Find the maximum number of divisors that fit into the
  372. -- dividend
  373. factor = 0
  374. while (bigint.compare(divisor, dividend, "<")) do
  375. divisor = bigint.add(divisor, big2)
  376. factor = factor + 1
  377. end
  378. -- Append the factor to the result
  379. if (factor == 10) then
  380. -- Fixes a weird bug that introduces a new bug if fixed by
  381. -- changing the comparison in the while loop to "<="
  382. result.digits[#result.digits] = 1
  383. result.digits[#result.digits + 1] = 0
  384. else
  385. result.digits[#result.digits + 1] = factor
  386. end
  387. -- Subtract the divisor from the dividend to obtain the
  388. -- remainder, which is the new dividend for the next loop
  389. dividend = bigint.subtract(dividend,
  390. bigint.subtract(divisor, big2))
  391. -- Reset the divisor
  392. divisor = big2:clone()
  393. end
  394. end
  395. -- The remainder of the final loop is returned as the function's
  396. -- overall remainder
  397. return result, dividend
  398. end
  399. end
  400. -- FRONTEND: Divide two bigs (decimals not supported), returning big result and
  401. -- big remainder, accounting for signs
  402. function bigint.divide(big1, big2)
  403. local result, remainder = bigint.divide_raw(bigint.abs(big1),
  404. bigint.abs(big2))
  405. if (big1.sign == big2.sign) then
  406. result.sign = "+"
  407. else
  408. result.sign = "-"
  409. end
  410. return result, remainder
  411. end
  412. -- FRONTEND: Return only the remainder from bigint.divide
  413. function bigint.modulus(big1, big2)
  414. local result, remainder = bigint.divide(big1, big2)
  415. -- Remainder will always have the same sign as the dividend per C standard
  416. -- https://en.wikipedia.org/wiki/Modulo_operation#Remainder_calculation_for_the_modulo_operation
  417. remainder.sign = big1.sign
  418. return remainder
  419. end
  420. return bigint