internal.odin 72 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653265426552656265726582659266026612662266326642665266626672668266926702671267226732674267526762677267826792680268126822683268426852686268726882689269026912692269326942695269626972698269927002701270227032704270527062707270827092710271127122713271427152716271727182719272027212722272327242725272627272728272927302731273227332734273527362737273827392740274127422743274427452746274727482749275027512752275327542755275627572758275927602761276227632764276527662767276827692770277127722773277427752776277727782779278027812782278327842785278627872788278927902791279227932794279527962797279827992800280128022803280428052806280728082809281028112812281328142815281628172818281928202821282228232824282528262827282828292830283128322833283428352836283728382839284028412842284328442845284628472848284928502851285228532854285528562857285828592860286128622863286428652866286728682869287028712872287328742875287628772878287928802881288228832884288528862887288828892890289128922893289428952896289728982899290029012902290329042905290629072908290929102911291229132914291529162917291829192920292129222923292429252926292729282929293029312932293329342935293629372938293929402941294229432944294529462947294829492950295129522953295429552956295729582959
  1. /*
  2. Copyright 2021 Jeroen van Rijn <[email protected]>.
  3. Made available under Odin's BSD-3 license.
  4. ========================== Low-level routines ==========================
  5. IMPORTANT: `internal_*` procedures make certain assumptions about their input.
  6. The public functions that call them are expected to satisfy their sanity check requirements.
  7. This allows `internal_*` call `internal_*` without paying this overhead multiple times.
  8. Where errors can occur, they are of course still checked and returned as appropriate.
  9. When importing `math:core/big` to implement an involved algorithm of your own, you are welcome
  10. to use these procedures instead of their public counterparts.
  11. Most inputs and outputs are expected to be passed an initialized `Int`, for example.
  12. Exceptions include `quotient` and `remainder`, which are allowed to be `nil` when the calling code doesn't need them.
  13. Check the comments above each `internal_*` implementation to see what constraints it expects to have met.
  14. We pass the custom allocator to procedures by default using the pattern `context.allocator = allocator`.
  15. This way we don't have to add `, allocator` at the end of each call.
  16. TODO: Handle +/- Infinity and NaN.
  17. */
  18. //+ignore
  19. package math_big
  20. import "core:mem"
  21. import "core:intrinsics"
  22. import rnd "core:math/rand"
  23. import "core:builtin"
  24. /*
  25. Low-level addition, unsigned. Handbook of Applied Cryptography, algorithm 14.7.
  26. Assumptions:
  27. `dest`, `a` and `b` != `nil` and have been initalized.
  28. */
  29. internal_int_add_unsigned :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
  30. dest := dest; x := a; y := b
  31. context.allocator = allocator
  32. old_used, min_used, max_used, i: int
  33. if x.used < y.used {
  34. x, y = y, x
  35. }
  36. min_used = y.used
  37. max_used = x.used
  38. old_used = dest.used
  39. internal_grow(dest, max(max_used + 1, _DEFAULT_DIGIT_COUNT)) or_return
  40. dest.used = max_used + 1
  41. /*
  42. All parameters have been initialized.
  43. */
  44. /* Zero the carry */
  45. carry := DIGIT(0)
  46. #no_bounds_check for i = 0; i < min_used; i += 1 {
  47. /*
  48. Compute the sum one _DIGIT at a time.
  49. dest[i] = a[i] + b[i] + carry;
  50. */
  51. dest.digit[i] = x.digit[i] + y.digit[i] + carry
  52. /*
  53. Compute carry
  54. */
  55. carry = dest.digit[i] >> _DIGIT_BITS
  56. /*
  57. Mask away carry from result digit.
  58. */
  59. dest.digit[i] &= _MASK
  60. }
  61. if min_used != max_used {
  62. /*
  63. Now copy higher words, if any, in A+B.
  64. If A or B has more digits, add those in.
  65. */
  66. #no_bounds_check for ; i < max_used; i += 1 {
  67. dest.digit[i] = x.digit[i] + carry
  68. /*
  69. Compute carry
  70. */
  71. carry = dest.digit[i] >> _DIGIT_BITS
  72. /*
  73. Mask away carry from result digit.
  74. */
  75. dest.digit[i] &= _MASK
  76. }
  77. }
  78. /*
  79. Add remaining carry.
  80. */
  81. dest.digit[i] = carry
  82. /*
  83. Zero remainder.
  84. */
  85. internal_zero_unused(dest, old_used)
  86. /*
  87. Adjust dest.used based on leading zeroes.
  88. */
  89. return internal_clamp(dest)
  90. }
  91. internal_add_unsigned :: proc { internal_int_add_unsigned, }
  92. /*
  93. Low-level addition, signed. Handbook of Applied Cryptography, algorithm 14.7.
  94. Assumptions:
  95. `dest`, `a` and `b` != `nil` and have been initalized.
  96. */
  97. internal_int_add_signed :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
  98. x := a; y := b
  99. context.allocator = allocator
  100. /*
  101. Handle both negative or both positive.
  102. */
  103. if x.sign == y.sign {
  104. dest.sign = x.sign
  105. return #force_inline internal_int_add_unsigned(dest, x, y)
  106. }
  107. /*
  108. One positive, the other negative.
  109. Subtract the one with the greater magnitude from the other.
  110. The result gets the sign of the one with the greater magnitude.
  111. */
  112. if #force_inline internal_lt_abs(a, b) {
  113. x, y = y, x
  114. }
  115. dest.sign = x.sign
  116. return #force_inline internal_int_sub_unsigned(dest, x, y)
  117. }
  118. internal_add_signed :: proc { internal_int_add_signed, }
  119. /*
  120. Low-level addition Int+DIGIT, signed. Handbook of Applied Cryptography, algorithm 14.7.
  121. Assumptions:
  122. `dest` and `a` != `nil` and have been initalized.
  123. `dest` is large enough (a.used + 1) to fit result.
  124. */
  125. internal_int_add_digit :: proc(dest, a: ^Int, digit: DIGIT, allocator := context.allocator) -> (err: Error) {
  126. context.allocator = allocator
  127. internal_grow(dest, a.used + 1) or_return
  128. /*
  129. Fast paths for destination and input Int being the same.
  130. */
  131. if dest == a {
  132. /*
  133. Fast path for dest.digit[0] + digit fits in dest.digit[0] without overflow.
  134. */
  135. if dest.sign == .Zero_or_Positive && (dest.digit[0] + digit < _DIGIT_MAX) {
  136. dest.digit[0] += digit
  137. dest.used += 1
  138. return internal_clamp(dest)
  139. }
  140. /*
  141. Can be subtracted from dest.digit[0] without underflow.
  142. */
  143. if a.sign == .Negative && (dest.digit[0] > digit) {
  144. dest.digit[0] -= digit
  145. dest.used += 1
  146. return internal_clamp(dest)
  147. }
  148. }
  149. /*
  150. If `a` is negative and `|a|` >= `digit`, call `dest = |a| - digit`
  151. */
  152. if a.sign == .Negative && (a.used > 1 || a.digit[0] >= digit) {
  153. /*
  154. Temporarily fix `a`'s sign.
  155. */
  156. a.sign = .Zero_or_Positive
  157. /*
  158. dest = |a| - digit
  159. */
  160. if err = #force_inline internal_int_add_digit(dest, a, digit); err != nil {
  161. /*
  162. Restore a's sign.
  163. */
  164. a.sign = .Negative
  165. return err
  166. }
  167. /*
  168. Restore sign and set `dest` sign.
  169. */
  170. a.sign = .Negative
  171. dest.sign = .Negative
  172. return internal_clamp(dest)
  173. }
  174. /*
  175. Remember the currently used number of digits in `dest`.
  176. */
  177. old_used := dest.used
  178. /*
  179. If `a` is positive
  180. */
  181. if a.sign == .Zero_or_Positive {
  182. /*
  183. Add digits, use `carry`.
  184. */
  185. i: int
  186. carry := digit
  187. #no_bounds_check for i = 0; i < a.used; i += 1 {
  188. dest.digit[i] = a.digit[i] + carry
  189. carry = dest.digit[i] >> _DIGIT_BITS
  190. dest.digit[i] &= _MASK
  191. }
  192. /*
  193. Set final carry.
  194. */
  195. dest.digit[i] = carry
  196. /*
  197. Set `dest` size.
  198. */
  199. dest.used = a.used + 1
  200. } else {
  201. /*
  202. `a` was negative and |a| < digit.
  203. */
  204. dest.used = 1
  205. /*
  206. The result is a single DIGIT.
  207. */
  208. dest.digit[0] = digit - a.digit[0] if a.used == 1 else digit
  209. }
  210. /*
  211. Sign is always positive.
  212. */
  213. dest.sign = .Zero_or_Positive
  214. /*
  215. Zero remainder.
  216. */
  217. internal_zero_unused(dest, old_used)
  218. /*
  219. Adjust dest.used based on leading zeroes.
  220. */
  221. return internal_clamp(dest)
  222. }
  223. internal_add :: proc { internal_int_add_signed, internal_int_add_digit, }
  224. internal_int_incr :: proc(dest: ^Int, allocator := context.allocator) -> (err: Error) {
  225. return #force_inline internal_add(dest, dest, 1)
  226. }
  227. internal_incr :: proc { internal_int_incr, }
  228. /*
  229. Low-level subtraction, dest = number - decrease. Assumes |number| > |decrease|.
  230. Handbook of Applied Cryptography, algorithm 14.9.
  231. Assumptions:
  232. `dest`, `number` and `decrease` != `nil` and have been initalized.
  233. */
  234. internal_int_sub_unsigned :: proc(dest, number, decrease: ^Int, allocator := context.allocator) -> (err: Error) {
  235. context.allocator = allocator
  236. dest := dest; x := number; y := decrease
  237. old_used := dest.used
  238. min_used := y.used
  239. max_used := x.used
  240. i: int
  241. grow(dest, max(max_used, _DEFAULT_DIGIT_COUNT)) or_return
  242. dest.used = max_used
  243. /*
  244. All parameters have been initialized.
  245. */
  246. borrow := DIGIT(0)
  247. #no_bounds_check for i = 0; i < min_used; i += 1 {
  248. dest.digit[i] = (x.digit[i] - y.digit[i] - borrow)
  249. /*
  250. borrow = carry bit of dest[i]
  251. Note this saves performing an AND operation since if a carry does occur,
  252. it will propagate all the way to the MSB.
  253. As a result a single shift is enough to get the carry.
  254. */
  255. borrow = dest.digit[i] >> ((size_of(DIGIT) * 8) - 1)
  256. /*
  257. Clear borrow from dest[i].
  258. */
  259. dest.digit[i] &= _MASK
  260. }
  261. /*
  262. Now copy higher words if any, e.g. if A has more digits than B
  263. */
  264. #no_bounds_check for ; i < max_used; i += 1 {
  265. dest.digit[i] = x.digit[i] - borrow
  266. /*
  267. borrow = carry bit of dest[i]
  268. Note this saves performing an AND operation since if a carry does occur,
  269. it will propagate all the way to the MSB.
  270. As a result a single shift is enough to get the carry.
  271. */
  272. borrow = dest.digit[i] >> ((size_of(DIGIT) * 8) - 1)
  273. /*
  274. Clear borrow from dest[i].
  275. */
  276. dest.digit[i] &= _MASK
  277. }
  278. /*
  279. Zero remainder.
  280. */
  281. internal_zero_unused(dest, old_used)
  282. /*
  283. Adjust dest.used based on leading zeroes.
  284. */
  285. return internal_clamp(dest)
  286. }
  287. internal_sub_unsigned :: proc { internal_int_sub_unsigned, }
  288. /*
  289. Low-level subtraction, signed. Handbook of Applied Cryptography, algorithm 14.9.
  290. dest = number - decrease. Assumes |number| > |decrease|.
  291. Assumptions:
  292. `dest`, `number` and `decrease` != `nil` and have been initalized.
  293. */
  294. internal_int_sub_signed :: proc(dest, number, decrease: ^Int, allocator := context.allocator) -> (err: Error) {
  295. context.allocator = allocator
  296. number := number; decrease := decrease
  297. if number.sign != decrease.sign {
  298. /*
  299. Subtract a negative from a positive, OR subtract a positive from a negative.
  300. In either case, ADD their magnitudes and use the sign of the first number.
  301. */
  302. dest.sign = number.sign
  303. return #force_inline internal_int_add_unsigned(dest, number, decrease)
  304. }
  305. /*
  306. Subtract a positive from a positive, OR negative from a negative.
  307. First, take the difference between their magnitudes, then...
  308. */
  309. if #force_inline internal_lt_abs(number, decrease) {
  310. /*
  311. The second has a larger magnitude.
  312. The result has the *opposite* sign from the first number.
  313. */
  314. dest.sign = .Negative if number.sign == .Zero_or_Positive else .Zero_or_Positive
  315. number, decrease = decrease, number
  316. } else {
  317. /*
  318. The first has a larger or equal magnitude.
  319. Copy the sign from the first.
  320. */
  321. dest.sign = number.sign
  322. }
  323. return #force_inline internal_int_sub_unsigned(dest, number, decrease)
  324. }
  325. /*
  326. Low-level subtraction, signed. Handbook of Applied Cryptography, algorithm 14.9.
  327. dest = number - decrease. Assumes |number| > |decrease|.
  328. Assumptions:
  329. `dest`, `number` != `nil` and have been initalized.
  330. `dest` is large enough (number.used + 1) to fit result.
  331. */
  332. internal_int_sub_digit :: proc(dest, number: ^Int, digit: DIGIT, allocator := context.allocator) -> (err: Error) {
  333. context.allocator = allocator
  334. internal_grow(dest, number.used + 1) or_return
  335. dest := dest; digit := digit
  336. /*
  337. All parameters have been initialized.
  338. Fast paths for destination and input Int being the same.
  339. */
  340. if dest == number {
  341. /*
  342. Fast path for `dest` is negative and unsigned addition doesn't overflow the lowest digit.
  343. */
  344. if dest.sign == .Negative && (dest.digit[0] + digit < _DIGIT_MAX) {
  345. dest.digit[0] += digit
  346. return nil
  347. }
  348. /*
  349. Can be subtracted from dest.digit[0] without underflow.
  350. */
  351. if number.sign == .Zero_or_Positive && (dest.digit[0] > digit) {
  352. dest.digit[0] -= digit
  353. return nil
  354. }
  355. }
  356. /*
  357. If `a` is negative, just do an unsigned addition (with fudged signs).
  358. */
  359. if number.sign == .Negative {
  360. t := number
  361. t.sign = .Zero_or_Positive
  362. err = #force_inline internal_int_add_digit(dest, t, digit)
  363. dest.sign = .Negative
  364. internal_clamp(dest)
  365. return err
  366. }
  367. old_used := dest.used
  368. /*
  369. if `a`<= digit, simply fix the single digit.
  370. */
  371. if number.used == 1 && (number.digit[0] <= digit) || number.used == 0 {
  372. dest.digit[0] = digit - number.digit[0] if number.used == 1 else digit
  373. dest.sign = .Negative
  374. dest.used = 1
  375. } else {
  376. dest.sign = .Zero_or_Positive
  377. dest.used = number.used
  378. /*
  379. Subtract with carry.
  380. */
  381. carry := digit
  382. #no_bounds_check for i := 0; i < number.used; i += 1 {
  383. dest.digit[i] = number.digit[i] - carry
  384. carry = dest.digit[i] >> (_DIGIT_TYPE_BITS - 1)
  385. dest.digit[i] &= _MASK
  386. }
  387. }
  388. /*
  389. Zero remainder.
  390. */
  391. internal_zero_unused(dest, old_used)
  392. /*
  393. Adjust dest.used based on leading zeroes.
  394. */
  395. return internal_clamp(dest)
  396. }
  397. internal_sub :: proc { internal_int_sub_signed, internal_int_sub_digit, }
  398. internal_int_decr :: proc(dest: ^Int, allocator := context.allocator) -> (err: Error) {
  399. return #force_inline internal_sub(dest, dest, 1)
  400. }
  401. internal_decr :: proc { internal_int_decr, }
  402. /*
  403. dest = src / 2
  404. dest = src >> 1
  405. Assumes `dest` and `src` not to be `nil` and have been initialized.
  406. We make no allocations here.
  407. */
  408. internal_int_shr1 :: proc(dest, src: ^Int) -> (err: Error) {
  409. old_used := dest.used; dest.used = src.used
  410. /*
  411. Carry
  412. */
  413. fwd_carry := DIGIT(0)
  414. #no_bounds_check for x := dest.used - 1; x >= 0; x -= 1 {
  415. /*
  416. Get the carry for the next iteration.
  417. */
  418. src_digit := src.digit[x]
  419. carry := src_digit & 1
  420. /*
  421. Shift the current digit, add in carry and store.
  422. */
  423. dest.digit[x] = (src_digit >> 1) | (fwd_carry << (_DIGIT_BITS - 1))
  424. /*
  425. Forward carry to next iteration.
  426. */
  427. fwd_carry = carry
  428. }
  429. /*
  430. Zero remainder.
  431. */
  432. internal_zero_unused(dest, old_used)
  433. /*
  434. Adjust dest.used based on leading zeroes.
  435. */
  436. dest.sign = src.sign
  437. return internal_clamp(dest)
  438. }
  439. /*
  440. dest = src * 2
  441. dest = src << 1
  442. */
  443. internal_int_shl1 :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
  444. context.allocator = allocator
  445. internal_copy(dest, src) or_return
  446. /*
  447. Grow `dest` to accommodate the additional bits.
  448. */
  449. digits_needed := dest.used + 1
  450. internal_grow(dest, digits_needed) or_return
  451. dest.used = digits_needed
  452. mask := (DIGIT(1) << uint(1)) - DIGIT(1)
  453. shift := DIGIT(_DIGIT_BITS - 1)
  454. carry := DIGIT(0)
  455. #no_bounds_check for x:= 0; x < dest.used; x+= 1 {
  456. fwd_carry := (dest.digit[x] >> shift) & mask
  457. dest.digit[x] = (dest.digit[x] << uint(1) | carry) & _MASK
  458. carry = fwd_carry
  459. }
  460. /*
  461. Use final carry.
  462. */
  463. if carry != 0 {
  464. dest.digit[dest.used] = carry
  465. dest.used += 1
  466. }
  467. return internal_clamp(dest)
  468. }
  469. /*
  470. Multiply bigint `a` with int `d` and put the result in `dest`.
  471. Like `internal_int_mul_digit` but with an integer as the small input.
  472. */
  473. internal_int_mul_integer :: proc(dest, a: ^Int, b: $T, allocator := context.allocator) -> (err: Error)
  474. where intrinsics.type_is_integer(T) && T != DIGIT {
  475. context.allocator = allocator
  476. t := &Int{}
  477. defer internal_destroy(t)
  478. /*
  479. DIGIT might be smaller than a long, which excludes the use of `internal_int_mul_digit` here.
  480. */
  481. internal_set(t, b) or_return
  482. internal_mul(dest, a, t) or_return
  483. return
  484. }
  485. /*
  486. Multiply by a DIGIT.
  487. */
  488. internal_int_mul_digit :: proc(dest, src: ^Int, multiplier: DIGIT, allocator := context.allocator) -> (err: Error) {
  489. context.allocator = allocator
  490. assert_if_nil(dest, src)
  491. if multiplier == 0 {
  492. return internal_zero(dest)
  493. }
  494. if multiplier == 1 {
  495. return internal_copy(dest, src)
  496. }
  497. /*
  498. Power of two?
  499. */
  500. if multiplier == 2 {
  501. return #force_inline internal_int_shl1(dest, src)
  502. }
  503. if #force_inline platform_int_is_power_of_two(int(multiplier)) {
  504. ix := internal_log(multiplier, 2) or_return
  505. return internal_shl(dest, src, ix)
  506. }
  507. /*
  508. Ensure `dest` is big enough to hold `src` * `multiplier`.
  509. */
  510. grow(dest, max(src.used + 1, _DEFAULT_DIGIT_COUNT)) or_return
  511. /*
  512. Save the original used count.
  513. */
  514. old_used := dest.used
  515. /*
  516. Set the sign.
  517. */
  518. dest.sign = src.sign
  519. /*
  520. Set up carry.
  521. */
  522. carry := _WORD(0)
  523. /*
  524. Compute columns.
  525. */
  526. ix := 0
  527. #no_bounds_check for ; ix < src.used; ix += 1 {
  528. /*
  529. Compute product and carry sum for this term
  530. */
  531. product := carry + _WORD(src.digit[ix]) * _WORD(multiplier)
  532. /*
  533. Mask off higher bits to get a single DIGIT.
  534. */
  535. dest.digit[ix] = DIGIT(product & _WORD(_MASK))
  536. /*
  537. Send carry into next iteration
  538. */
  539. carry = product >> _DIGIT_BITS
  540. }
  541. /*
  542. Store final carry [if any] and increment used.
  543. */
  544. dest.digit[ix] = DIGIT(carry)
  545. dest.used = src.used + 1
  546. /*
  547. Zero remainder.
  548. */
  549. internal_zero_unused(dest, old_used)
  550. return internal_clamp(dest)
  551. }
  552. /*
  553. High level multiplication (handles sign).
  554. */
  555. internal_int_mul :: proc(dest, src, multiplier: ^Int, allocator := context.allocator) -> (err: Error) {
  556. context.allocator = allocator
  557. /*
  558. Early out for `multiplier` is zero; Set `dest` to zero.
  559. */
  560. if multiplier.used == 0 || src.used == 0 { return internal_zero(dest) }
  561. neg := src.sign != multiplier.sign
  562. if src == multiplier {
  563. /*
  564. Do we need to square?
  565. */
  566. if src.used >= SQR_TOOM_CUTOFF {
  567. /*
  568. Use Toom-Cook?
  569. */
  570. err = #force_inline _private_int_sqr_toom(dest, src)
  571. } else if src.used >= SQR_KARATSUBA_CUTOFF {
  572. /*
  573. Karatsuba?
  574. */
  575. err = #force_inline _private_int_sqr_karatsuba(dest, src)
  576. } else if ((src.used * 2) + 1) < _WARRAY && src.used < (_MAX_COMBA / 2) {
  577. /*
  578. Fast comba?
  579. */
  580. err = #force_inline _private_int_sqr_comba(dest, src)
  581. } else {
  582. err = #force_inline _private_int_sqr(dest, src)
  583. }
  584. } else {
  585. /*
  586. Can we use the balance method? Check sizes.
  587. * The smaller one needs to be larger than the Karatsuba cut-off.
  588. * The bigger one needs to be at least about one `_MUL_KARATSUBA_CUTOFF` bigger
  589. * to make some sense, but it depends on architecture, OS, position of the stars... so YMMV.
  590. * Using it to cut the input into slices small enough for _mul_comba
  591. * was actually slower on the author's machine, but YMMV.
  592. */
  593. min_used := min(src.used, multiplier.used)
  594. max_used := max(src.used, multiplier.used)
  595. digits := src.used + multiplier.used + 1
  596. if min_used >= MUL_KARATSUBA_CUTOFF && (max_used / 2) >= MUL_KARATSUBA_CUTOFF && max_used >= (2 * min_used) {
  597. /*
  598. Not much effect was observed below a ratio of 1:2, but again: YMMV.
  599. */
  600. err = _private_int_mul_balance(dest, src, multiplier)
  601. } else if min_used >= MUL_TOOM_CUTOFF {
  602. /*
  603. Toom path commented out until it no longer fails Factorial 10k or 100k,
  604. as reveaved in the long test.
  605. */
  606. err = #force_inline _private_int_mul_toom(dest, src, multiplier)
  607. } else if min_used >= MUL_KARATSUBA_CUTOFF {
  608. err = #force_inline _private_int_mul_karatsuba(dest, src, multiplier)
  609. } else if digits < _WARRAY && min_used <= _MAX_COMBA {
  610. /*
  611. Can we use the fast multiplier?
  612. * The fast multiplier can be used if the output will
  613. * have less than MP_WARRAY digits and the number of
  614. * digits won't affect carry propagation
  615. */
  616. err = #force_inline _private_int_mul_comba(dest, src, multiplier, digits)
  617. } else {
  618. err = #force_inline _private_int_mul(dest, src, multiplier, digits)
  619. }
  620. }
  621. dest.sign = .Negative if dest.used > 0 && neg else .Zero_or_Positive
  622. return err
  623. }
  624. internal_mul :: proc { internal_int_mul, internal_int_mul_digit, internal_int_mul_integer }
  625. internal_sqr :: proc (dest, src: ^Int, allocator := context.allocator) -> (res: Error) {
  626. /*
  627. We call `internal_mul` and not e.g. `_private_int_sqr` because the former
  628. will dispatch to the optimal implementation depending on the source.
  629. */
  630. return #force_inline internal_mul(dest, src, src, allocator)
  631. }
  632. /*
  633. divmod.
  634. Both the quotient and remainder are optional and may be passed a nil.
  635. `numerator` and `denominator` are expected not to be `nil` and have been initialized.
  636. */
  637. internal_int_divmod :: proc(quotient, remainder, numerator, denominator: ^Int, allocator := context.allocator) -> (err: Error) {
  638. context.allocator = allocator
  639. if denominator.used == 0 { return .Division_by_Zero }
  640. /*
  641. If numerator < denominator then quotient = 0, remainder = numerator.
  642. */
  643. if #force_inline internal_lt_abs(numerator, denominator) {
  644. if remainder != nil {
  645. internal_copy(remainder, numerator) or_return
  646. }
  647. if quotient != nil {
  648. internal_zero(quotient)
  649. }
  650. return nil
  651. }
  652. if (denominator.used > 2 * MUL_KARATSUBA_CUTOFF) && (denominator.used <= (numerator.used / 3) * 2) {
  653. assert(denominator.used >= 160 && numerator.used >= 240, "MUL_KARATSUBA_CUTOFF global not properly set.")
  654. err = _private_int_div_recursive(quotient, remainder, numerator, denominator)
  655. } else {
  656. when true {
  657. err = #force_inline _private_int_div_school(quotient, remainder, numerator, denominator)
  658. } else {
  659. /*
  660. NOTE(Jeroen): We no longer need or use `_private_int_div_small`.
  661. We'll keep it around for a bit until we're reasonably certain div_school is bug free.
  662. */
  663. err = _private_int_div_small(quotient, remainder, numerator, denominator)
  664. }
  665. }
  666. return
  667. }
  668. /*
  669. Single digit division (based on routine from MPI).
  670. The quotient is optional and may be passed a nil.
  671. */
  672. internal_int_divmod_digit :: proc(quotient, numerator: ^Int, denominator: DIGIT, allocator := context.allocator) -> (remainder: DIGIT, err: Error) {
  673. context.allocator = allocator
  674. /*
  675. Cannot divide by zero.
  676. */
  677. if denominator == 0 { return 0, .Division_by_Zero }
  678. /*
  679. Quick outs.
  680. */
  681. if denominator == 1 || numerator.used == 0 {
  682. if quotient != nil {
  683. return 0, internal_copy(quotient, numerator)
  684. }
  685. return 0, err
  686. }
  687. /*
  688. Power of two?
  689. */
  690. if denominator == 2 {
  691. if numerator.used > 0 && numerator.digit[0] & 1 != 0 {
  692. // Remainder is 1 if numerator is odd.
  693. remainder = 1
  694. }
  695. if quotient == nil {
  696. return remainder, nil
  697. }
  698. return remainder, internal_shr(quotient, numerator, 1)
  699. }
  700. ix: int
  701. if platform_int_is_power_of_two(int(denominator)) {
  702. ix = 1
  703. for ix < _DIGIT_BITS && denominator != (1 << uint(ix)) {
  704. ix += 1
  705. }
  706. remainder = numerator.digit[0] & ((1 << uint(ix)) - 1)
  707. if quotient == nil {
  708. return remainder, nil
  709. }
  710. return remainder, internal_shr(quotient, numerator, int(ix))
  711. }
  712. /*
  713. Three?
  714. */
  715. if denominator == 3 {
  716. return _private_int_div_3(quotient, numerator)
  717. }
  718. /*
  719. No easy answer [c'est la vie]. Just division.
  720. */
  721. q := &Int{}
  722. internal_grow(q, numerator.used) or_return
  723. q.used = numerator.used
  724. q.sign = numerator.sign
  725. w := _WORD(0)
  726. for ix = numerator.used - 1; ix >= 0; ix -= 1 {
  727. t := DIGIT(0)
  728. w = (w << _WORD(_DIGIT_BITS) | _WORD(numerator.digit[ix]))
  729. if w >= _WORD(denominator) {
  730. t = DIGIT(w / _WORD(denominator))
  731. w -= _WORD(t) * _WORD(denominator)
  732. }
  733. q.digit[ix] = t
  734. }
  735. remainder = DIGIT(w)
  736. if quotient != nil {
  737. internal_clamp(q)
  738. internal_swap(q, quotient)
  739. }
  740. internal_destroy(q)
  741. return remainder, nil
  742. }
  743. internal_divmod :: proc { internal_int_divmod, internal_int_divmod_digit, }
  744. /*
  745. Asssumes quotient, numerator and denominator to have been initialized and not to be nil.
  746. */
  747. internal_int_div :: proc(quotient, numerator, denominator: ^Int, allocator := context.allocator) -> (err: Error) {
  748. return #force_inline internal_int_divmod(quotient, nil, numerator, denominator, allocator)
  749. }
  750. internal_div :: proc { internal_int_div, }
  751. /*
  752. remainder = numerator % denominator.
  753. 0 <= remainder < denominator if denominator > 0
  754. denominator < remainder <= 0 if denominator < 0
  755. Asssumes quotient, numerator and denominator to have been initialized and not to be nil.
  756. */
  757. internal_int_mod :: proc(remainder, numerator, denominator: ^Int, allocator := context.allocator) -> (err: Error) {
  758. #force_inline internal_int_divmod(nil, remainder, numerator, denominator, allocator) or_return
  759. if remainder.used == 0 || denominator.sign == remainder.sign { return nil }
  760. return #force_inline internal_add(remainder, remainder, denominator, allocator)
  761. }
  762. internal_int_mod_digit :: proc(numerator: ^Int, denominator: DIGIT, allocator := context.allocator) -> (remainder: DIGIT, err: Error) {
  763. return internal_int_divmod_digit(nil, numerator, denominator, allocator)
  764. }
  765. internal_mod :: proc{ internal_int_mod, internal_int_mod_digit, }
  766. /*
  767. remainder = (number + addend) % modulus.
  768. */
  769. internal_int_addmod :: proc(remainder, number, addend, modulus: ^Int, allocator := context.allocator) -> (err: Error) {
  770. #force_inline internal_add(remainder, number, addend, allocator) or_return
  771. return #force_inline internal_mod(remainder, remainder, modulus, allocator)
  772. }
  773. internal_addmod :: proc { internal_int_addmod, }
  774. /*
  775. remainder = (number - decrease) % modulus.
  776. */
  777. internal_int_submod :: proc(remainder, number, decrease, modulus: ^Int, allocator := context.allocator) -> (err: Error) {
  778. #force_inline internal_sub(remainder, number, decrease, allocator) or_return
  779. return #force_inline internal_mod(remainder, remainder, modulus, allocator)
  780. }
  781. internal_submod :: proc { internal_int_submod, }
  782. /*
  783. remainder = (number * multiplicand) % modulus.
  784. */
  785. internal_int_mulmod :: proc(remainder, number, multiplicand, modulus: ^Int, allocator := context.allocator) -> (err: Error) {
  786. #force_inline internal_mul(remainder, number, multiplicand, allocator) or_return
  787. return #force_inline internal_mod(remainder, remainder, modulus, allocator)
  788. }
  789. internal_mulmod :: proc { internal_int_mulmod, }
  790. /*
  791. remainder = (number * number) % modulus.
  792. */
  793. internal_int_sqrmod :: proc(remainder, number, modulus: ^Int, allocator := context.allocator) -> (err: Error) {
  794. #force_inline internal_sqr(remainder, number, allocator) or_return
  795. return #force_inline internal_mod(remainder, remainder, modulus, allocator)
  796. }
  797. internal_sqrmod :: proc { internal_int_sqrmod, }
  798. /*
  799. TODO: Use Sterling's Approximation to estimate log2(N!) to size the result.
  800. This way we'll have to reallocate less, possibly not at all.
  801. */
  802. internal_int_factorial :: proc(res: ^Int, n: int, allocator := context.allocator) -> (err: Error) {
  803. context.allocator = allocator
  804. if n >= FACTORIAL_BINARY_SPLIT_CUTOFF {
  805. return _private_int_factorial_binary_split(res, n)
  806. }
  807. i := len(_factorial_table)
  808. if n < i {
  809. return #force_inline internal_set(res, _factorial_table[n])
  810. }
  811. #force_inline internal_set(res, _factorial_table[i - 1]) or_return
  812. for {
  813. if err = #force_inline internal_mul(res, res, DIGIT(i)); err != nil || i == n {
  814. return err
  815. }
  816. i += 1
  817. }
  818. return nil
  819. }
  820. /*
  821. Returns GCD, LCM or both.
  822. Assumes `a` and `b` to have been initialized.
  823. `res_gcd` and `res_lcm` can be nil or ^Int depending on which results are desired.
  824. */
  825. internal_int_gcd_lcm :: proc(res_gcd, res_lcm, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
  826. if res_gcd == nil && res_lcm == nil { return nil }
  827. return #force_inline _private_int_gcd_lcm(res_gcd, res_lcm, a, b, allocator)
  828. }
  829. internal_int_gcd :: proc(res_gcd, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
  830. return #force_inline _private_int_gcd_lcm(res_gcd, nil, a, b, allocator)
  831. }
  832. internal_int_lcm :: proc(res_lcm, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
  833. return #force_inline _private_int_gcd_lcm(nil, res_lcm, a, b, allocator)
  834. }
  835. /*
  836. remainder = numerator % (1 << bits)
  837. Assumes `remainder` and `numerator` both not to be `nil` and `bits` to be >= 0.
  838. */
  839. internal_int_mod_bits :: proc(remainder, numerator: ^Int, bits: int, allocator := context.allocator) -> (err: Error) {
  840. /*
  841. Everything is divisible by 1 << 0 == 1, so this returns 0.
  842. */
  843. if bits == 0 { return internal_zero(remainder) }
  844. /*
  845. If the modulus is larger than the value, return the value.
  846. */
  847. internal_copy(remainder, numerator) or_return
  848. if bits >= (numerator.used * _DIGIT_BITS) {
  849. return
  850. }
  851. /*
  852. Zero digits above the last digit of the modulus.
  853. */
  854. zero_count := (bits / _DIGIT_BITS)
  855. zero_count += 0 if (bits % _DIGIT_BITS == 0) else 1
  856. /*
  857. Zero remainder. Special case, can't use `internal_zero_unused`.
  858. */
  859. if zero_count > 0 {
  860. mem.zero_slice(remainder.digit[zero_count:])
  861. }
  862. /*
  863. Clear the digit that is not completely outside/inside the modulus.
  864. */
  865. remainder.digit[bits / _DIGIT_BITS] &= DIGIT(1 << DIGIT(bits % _DIGIT_BITS)) - DIGIT(1)
  866. return internal_clamp(remainder)
  867. }
  868. /*
  869. ============================= Low-level helpers =============================
  870. `internal_*` helpers don't return an `Error` like their public counterparts do,
  871. because they expect not to be passed `nil` or uninitialized inputs.
  872. This makes them more suitable for `internal_*` functions and some of the
  873. public ones that have already satisfied these constraints.
  874. */
  875. /*
  876. This procedure returns the allocated capacity of an Int.
  877. Assumes `a` not to be `nil`.
  878. */
  879. internal_int_allocated_cap :: #force_inline proc(a: ^Int) -> (cap: int) {
  880. raw := transmute(mem.Raw_Dynamic_Array)a.digit
  881. return raw.cap
  882. }
  883. /*
  884. This procedure will return `true` if the `Int` is initialized, `false` if not.
  885. Assumes `a` not to be `nil`.
  886. */
  887. internal_int_is_initialized :: #force_inline proc(a: ^Int) -> (initialized: bool) {
  888. return internal_int_allocated_cap(a) >= _MIN_DIGIT_COUNT
  889. }
  890. internal_is_initialized :: proc { internal_int_is_initialized, }
  891. /*
  892. This procedure will return `true` if the `Int` is zero, `false` if not.
  893. Assumes `a` not to be `nil`.
  894. */
  895. internal_int_is_zero :: #force_inline proc(a: ^Int) -> (zero: bool) {
  896. return a.used == 0
  897. }
  898. internal_is_zero :: proc {
  899. internal_rat_is_zero,
  900. internal_int_is_zero,
  901. }
  902. /*
  903. This procedure will return `true` if the `Int` is positive, `false` if not.
  904. Assumes `a` not to be `nil`.
  905. */
  906. internal_int_is_positive :: #force_inline proc(a: ^Int) -> (positive: bool) {
  907. return a.sign == .Zero_or_Positive
  908. }
  909. internal_is_positive :: proc { internal_int_is_positive, }
  910. /*
  911. This procedure will return `true` if the `Int` is negative, `false` if not.
  912. Assumes `a` not to be `nil`.
  913. */
  914. internal_int_is_negative :: #force_inline proc(a: ^Int) -> (negative: bool) {
  915. return a.sign == .Negative
  916. }
  917. internal_is_negative :: proc { internal_int_is_negative, }
  918. /*
  919. This procedure will return `true` if the `Int` is even, `false` if not.
  920. Assumes `a` not to be `nil`.
  921. */
  922. internal_int_is_even :: #force_inline proc(a: ^Int) -> (even: bool) {
  923. if internal_is_zero(a) { return true }
  924. /*
  925. `a.used` > 0 here, because the above handled `is_zero`.
  926. We don't need to explicitly test it.
  927. */
  928. return a.digit[0] & 1 == 0
  929. }
  930. internal_is_even :: proc { internal_int_is_even, }
  931. /*
  932. This procedure will return `true` if the `Int` is even, `false` if not.
  933. Assumes `a` not to be `nil`.
  934. */
  935. internal_int_is_odd :: #force_inline proc(a: ^Int) -> (odd: bool) {
  936. return !internal_int_is_even(a)
  937. }
  938. internal_is_odd :: proc { internal_int_is_odd, }
  939. /*
  940. This procedure will return `true` if the `Int` is a power of two, `false` if not.
  941. Assumes `a` not to be `nil`.
  942. */
  943. internal_int_is_power_of_two :: #force_inline proc(a: ^Int) -> (power_of_two: bool) {
  944. /*
  945. Early out for Int == 0.
  946. */
  947. if #force_inline internal_is_zero(a) { return true }
  948. /*
  949. For an `Int` to be a power of two, its bottom limb has to be a power of two.
  950. */
  951. if ! #force_inline platform_int_is_power_of_two(int(a.digit[a.used - 1])) { return false }
  952. /*
  953. We've established that the bottom limb is a power of two.
  954. If it's the only limb, that makes the entire Int a power of two.
  955. */
  956. if a.used == 1 { return true }
  957. /*
  958. For an `Int` to be a power of two, all limbs except the top one have to be zero.
  959. */
  960. for i := 1; i < a.used && a.digit[i - 1] != 0; i += 1 { return false }
  961. return true
  962. }
  963. internal_is_power_of_two :: proc { internal_int_is_power_of_two, }
  964. /*
  965. Compare two `Int`s, signed.
  966. Returns -1 if `a` < `b`, 0 if `a` == `b` and 1 if `b` > `a`.
  967. Expects `a` and `b` both to be valid `Int`s, i.e. initialized and not `nil`.
  968. */
  969. internal_int_compare :: #force_inline proc(a, b: ^Int) -> (comparison: int) {
  970. assert_if_nil(a, b)
  971. a_is_negative := #force_inline internal_is_negative(a)
  972. /*
  973. Compare based on sign.
  974. */
  975. if a.sign != b.sign { return -1 if a_is_negative else +1 }
  976. /*
  977. If `a` is negative, compare in the opposite direction */
  978. if a_is_negative { return #force_inline internal_compare_magnitude(b, a) }
  979. return #force_inline internal_compare_magnitude(a, b)
  980. }
  981. internal_compare :: proc { internal_int_compare, internal_int_compare_digit, }
  982. internal_cmp :: internal_compare
  983. /*
  984. Compare an `Int` to an unsigned number upto `DIGIT & _MASK`.
  985. Returns -1 if `a` < `b`, 0 if `a` == `b` and 1 if `b` > `a`.
  986. Expects: `a` and `b` both to be valid `Int`s, i.e. initialized and not `nil`.
  987. */
  988. internal_int_compare_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (comparison: int) {
  989. assert_if_nil(a)
  990. a_is_negative := #force_inline internal_is_negative(a)
  991. switch {
  992. /*
  993. Compare based on sign first.
  994. */
  995. case a_is_negative: return -1
  996. /*
  997. Then compare on magnitude.
  998. */
  999. case a.used > 1: return +1
  1000. /*
  1001. We have only one digit. Compare it against `b`.
  1002. */
  1003. case a.digit[0] < b: return -1
  1004. case a.digit[0] == b: return 0
  1005. case a.digit[0] > b: return +1
  1006. /*
  1007. Unreachable.
  1008. Just here because Odin complains about a missing return value at the bottom of the proc otherwise.
  1009. */
  1010. case: return
  1011. }
  1012. }
  1013. internal_compare_digit :: proc { internal_int_compare_digit, }
  1014. internal_cmp_digit :: internal_compare_digit
  1015. /*
  1016. Compare the magnitude of two `Int`s, unsigned.
  1017. */
  1018. internal_int_compare_magnitude :: #force_inline proc(a, b: ^Int) -> (comparison: int) {
  1019. assert_if_nil(a, b)
  1020. /*
  1021. Compare based on used digits.
  1022. */
  1023. if a.used != b.used {
  1024. if a.used > b.used {
  1025. return +1
  1026. }
  1027. return -1
  1028. }
  1029. /*
  1030. Same number of used digits, compare based on their value.
  1031. */
  1032. #no_bounds_check for n := a.used - 1; n >= 0; n -= 1 {
  1033. if a.digit[n] != b.digit[n] {
  1034. if a.digit[n] > b.digit[n] {
  1035. return +1
  1036. }
  1037. return -1
  1038. }
  1039. }
  1040. return 0
  1041. }
  1042. internal_compare_magnitude :: proc { internal_int_compare_magnitude, }
  1043. internal_cmp_mag :: internal_compare_magnitude
  1044. /*
  1045. bool := a < b
  1046. */
  1047. internal_int_less_than :: #force_inline proc(a, b: ^Int) -> (less_than: bool) {
  1048. return internal_cmp(a, b) == -1
  1049. }
  1050. /*
  1051. bool := a < b
  1052. */
  1053. internal_int_less_than_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (less_than: bool) {
  1054. return internal_cmp_digit(a, b) == -1
  1055. }
  1056. /*
  1057. bool := |a| < |b|
  1058. Compares the magnitudes only, ignores the sign.
  1059. */
  1060. internal_int_less_than_abs :: #force_inline proc(a, b: ^Int) -> (less_than: bool) {
  1061. return internal_cmp_mag(a, b) == -1
  1062. }
  1063. internal_less_than :: proc {
  1064. internal_int_less_than,
  1065. internal_int_less_than_digit,
  1066. }
  1067. internal_lt :: internal_less_than
  1068. internal_less_than_abs :: proc {
  1069. internal_int_less_than_abs,
  1070. }
  1071. internal_lt_abs :: internal_less_than_abs
  1072. /*
  1073. bool := a <= b
  1074. */
  1075. internal_int_less_than_or_equal :: #force_inline proc(a, b: ^Int) -> (less_than_or_equal: bool) {
  1076. return internal_cmp(a, b) <= 0
  1077. }
  1078. /*
  1079. bool := a <= b
  1080. */
  1081. internal_int_less_than_or_equal_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (less_than_or_equal: bool) {
  1082. return internal_cmp_digit(a, b) <= 0
  1083. }
  1084. /*
  1085. bool := |a| <= |b|
  1086. Compares the magnitudes only, ignores the sign.
  1087. */
  1088. internal_int_less_than_or_equal_abs :: #force_inline proc(a, b: ^Int) -> (less_than_or_equal: bool) {
  1089. return internal_cmp_mag(a, b) <= 0
  1090. }
  1091. internal_less_than_or_equal :: proc {
  1092. internal_int_less_than_or_equal,
  1093. internal_int_less_than_or_equal_digit,
  1094. }
  1095. internal_lte :: internal_less_than_or_equal
  1096. internal_less_than_or_equal_abs :: proc {
  1097. internal_int_less_than_or_equal_abs,
  1098. }
  1099. internal_lte_abs :: internal_less_than_or_equal_abs
  1100. /*
  1101. bool := a == b
  1102. */
  1103. internal_int_equals :: #force_inline proc(a, b: ^Int) -> (equals: bool) {
  1104. return internal_cmp(a, b) == 0
  1105. }
  1106. /*
  1107. bool := a == b
  1108. */
  1109. internal_int_equals_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (equals: bool) {
  1110. return internal_cmp_digit(a, b) == 0
  1111. }
  1112. /*
  1113. bool := |a| == |b|
  1114. Compares the magnitudes only, ignores the sign.
  1115. */
  1116. internal_int_equals_abs :: #force_inline proc(a, b: ^Int) -> (equals: bool) {
  1117. return internal_cmp_mag(a, b) == 0
  1118. }
  1119. internal_equals :: proc {
  1120. internal_int_equals,
  1121. internal_int_equals_digit,
  1122. }
  1123. internal_eq :: internal_equals
  1124. internal_equals_abs :: proc {
  1125. internal_int_equals_abs,
  1126. }
  1127. internal_eq_abs :: internal_equals_abs
  1128. /*
  1129. bool := a >= b
  1130. */
  1131. internal_int_greater_than_or_equal :: #force_inline proc(a, b: ^Int) -> (greater_than_or_equal: bool) {
  1132. return internal_cmp(a, b) >= 0
  1133. }
  1134. /*
  1135. bool := a >= b
  1136. */
  1137. internal_int_greater_than_or_equal_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (greater_than_or_equal: bool) {
  1138. return internal_cmp_digit(a, b) >= 0
  1139. }
  1140. /*
  1141. bool := |a| >= |b|
  1142. Compares the magnitudes only, ignores the sign.
  1143. */
  1144. internal_int_greater_than_or_equal_abs :: #force_inline proc(a, b: ^Int) -> (greater_than_or_equal: bool) {
  1145. return internal_cmp_mag(a, b) >= 0
  1146. }
  1147. internal_greater_than_or_equal :: proc {
  1148. internal_int_greater_than_or_equal,
  1149. internal_int_greater_than_or_equal_digit,
  1150. }
  1151. internal_gte :: internal_greater_than_or_equal
  1152. internal_greater_than_or_equal_abs :: proc {
  1153. internal_int_greater_than_or_equal_abs,
  1154. }
  1155. internal_gte_abs :: internal_greater_than_or_equal_abs
  1156. /*
  1157. bool := a > b
  1158. */
  1159. internal_int_greater_than :: #force_inline proc(a, b: ^Int) -> (greater_than: bool) {
  1160. return internal_cmp(a, b) == 1
  1161. }
  1162. /*
  1163. bool := a > b
  1164. */
  1165. internal_int_greater_than_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (greater_than: bool) {
  1166. return internal_cmp_digit(a, b) == 1
  1167. }
  1168. /*
  1169. bool := |a| > |b|
  1170. Compares the magnitudes only, ignores the sign.
  1171. */
  1172. internal_int_greater_than_abs :: #force_inline proc(a, b: ^Int) -> (greater_than: bool) {
  1173. return internal_cmp_mag(a, b) == 1
  1174. }
  1175. internal_greater_than :: proc {
  1176. internal_int_greater_than,
  1177. internal_int_greater_than_digit,
  1178. }
  1179. internal_gt :: internal_greater_than
  1180. internal_greater_than_abs :: proc {
  1181. internal_int_greater_than_abs,
  1182. }
  1183. internal_gt_abs :: internal_greater_than_abs
  1184. /*
  1185. Check if remainders are possible squares - fast exclude non-squares.
  1186. Returns `true` if `a` is a square, `false` if not.
  1187. Assumes `a` not to be `nil` and to have been initialized.
  1188. */
  1189. internal_int_is_square :: proc(a: ^Int, allocator := context.allocator) -> (square: bool, err: Error) {
  1190. context.allocator = allocator
  1191. /*
  1192. Default to Non-square :)
  1193. */
  1194. square = false
  1195. if internal_is_negative(a) { return }
  1196. if internal_is_zero(a) { return }
  1197. /*
  1198. First check mod 128 (suppose that _DIGIT_BITS is at least 7).
  1199. */
  1200. if _private_int_rem_128[127 & a.digit[0]] == 1 { return }
  1201. /*
  1202. Next check mod 105 (3*5*7).
  1203. */
  1204. c: DIGIT
  1205. c, err = internal_mod(a, 105)
  1206. if _private_int_rem_105[c] == 1 { return }
  1207. t := &Int{}
  1208. defer destroy(t)
  1209. set(t, 11 * 13 * 17 * 19 * 23 * 29 * 31) or_return
  1210. internal_mod(t, a, t) or_return
  1211. r: u64
  1212. r, err = internal_int_get(t, u64)
  1213. /*
  1214. Check for other prime modules, note it's not an ERROR but we must
  1215. free "t" so the easiest way is to goto LBL_ERR. We know that err
  1216. is already equal to MP_OKAY from the mp_mod call
  1217. */
  1218. if (1 << (r % 11) & 0x5C4) != 0 { return }
  1219. if (1 << (r % 13) & 0x9E4) != 0 { return }
  1220. if (1 << (r % 17) & 0x5CE8) != 0 { return }
  1221. if (1 << (r % 19) & 0x4F50C) != 0 { return }
  1222. if (1 << (r % 23) & 0x7ACCA0) != 0 { return }
  1223. if (1 << (r % 29) & 0xC2EDD0C) != 0 { return }
  1224. if (1 << (r % 31) & 0x6DE2B848) != 0 { return }
  1225. /*
  1226. Final check - is sqr(sqrt(arg)) == arg?
  1227. */
  1228. sqrt(t, a) or_return
  1229. sqr(t, t) or_return
  1230. square = internal_eq_abs(t, a)
  1231. return
  1232. }
  1233. /*
  1234. ========================= Logs, powers and roots ============================
  1235. */
  1236. /*
  1237. Returns log_base(a).
  1238. Assumes `a` to not be `nil` and have been iniialized.
  1239. */
  1240. internal_int_log :: proc(a: ^Int, base: DIGIT) -> (res: int, err: Error) {
  1241. if base < 2 || DIGIT(base) > _DIGIT_MAX { return -1, .Invalid_Argument }
  1242. if internal_is_negative(a) { return -1, .Math_Domain_Error }
  1243. if internal_is_zero(a) { return -1, .Math_Domain_Error }
  1244. /*
  1245. Fast path for bases that are a power of two.
  1246. */
  1247. if platform_int_is_power_of_two(int(base)) { return _private_log_power_of_two(a, base) }
  1248. /*
  1249. Fast path for `Int`s that fit within a single `DIGIT`.
  1250. */
  1251. if a.used == 1 { return internal_log(a.digit[0], DIGIT(base)) }
  1252. return _private_int_log(a, base)
  1253. }
  1254. /*
  1255. Returns log_base(a), where `a` is a DIGIT.
  1256. */
  1257. internal_digit_log :: proc(a: DIGIT, base: DIGIT) -> (log: int, err: Error) {
  1258. /*
  1259. If the number is smaller than the base, it fits within a fraction.
  1260. Therefore, we return 0.
  1261. */
  1262. if a < base { return 0, nil }
  1263. /*
  1264. If a number equals the base, the log is 1.
  1265. */
  1266. if a == base { return 1, nil }
  1267. N := _WORD(a)
  1268. bracket_low := _WORD(1)
  1269. bracket_high := _WORD(base)
  1270. high := 1
  1271. low := 0
  1272. for bracket_high < N {
  1273. low = high
  1274. bracket_low = bracket_high
  1275. high <<= 1
  1276. bracket_high *= bracket_high
  1277. }
  1278. for high - low > 1 {
  1279. mid := (low + high) >> 1
  1280. bracket_mid := bracket_low * #force_inline internal_small_pow(_WORD(base), _WORD(mid - low))
  1281. if N < bracket_mid {
  1282. high = mid
  1283. bracket_high = bracket_mid
  1284. }
  1285. if N > bracket_mid {
  1286. low = mid
  1287. bracket_low = bracket_mid
  1288. }
  1289. if N == bracket_mid {
  1290. return mid, nil
  1291. }
  1292. }
  1293. if bracket_high == N {
  1294. return high, nil
  1295. } else {
  1296. return low, nil
  1297. }
  1298. }
  1299. internal_log :: proc { internal_int_log, internal_digit_log, }
  1300. /*
  1301. Calculate dest = base^power using a square-multiply algorithm.
  1302. Assumes `dest` and `base` not to be `nil` and to have been initialized.
  1303. */
  1304. internal_int_pow :: proc(dest, base: ^Int, power: int, allocator := context.allocator) -> (err: Error) {
  1305. context.allocator = allocator
  1306. power := power
  1307. /*
  1308. Early outs.
  1309. */
  1310. if #force_inline internal_is_zero(base) {
  1311. /*
  1312. A zero base is a special case.
  1313. */
  1314. if power < 0 {
  1315. internal_zero(dest) or_return
  1316. return .Math_Domain_Error
  1317. }
  1318. if power == 0 { return internal_one(dest) }
  1319. if power > 0 { return internal_zero(dest) }
  1320. }
  1321. if power < 0 {
  1322. /*
  1323. Fraction, so we'll return zero.
  1324. */
  1325. return internal_zero(dest)
  1326. }
  1327. switch(power) {
  1328. case 0:
  1329. /*
  1330. Any base to the power zero is one.
  1331. */
  1332. return #force_inline internal_one(dest)
  1333. case 1:
  1334. /*
  1335. Any base to the power one is itself.
  1336. */
  1337. return copy(dest, base)
  1338. case 2:
  1339. return #force_inline internal_sqr(dest, base)
  1340. }
  1341. g := &Int{}
  1342. internal_copy(g, base) or_return
  1343. /*
  1344. Set initial result.
  1345. */
  1346. internal_one(dest) or_return
  1347. defer internal_destroy(g)
  1348. for power > 0 {
  1349. /*
  1350. If the bit is set, multiply.
  1351. */
  1352. if power & 1 != 0 {
  1353. internal_mul(dest, g, dest) or_return
  1354. }
  1355. /*
  1356. Square.
  1357. */
  1358. if power > 1 {
  1359. internal_sqr(g, g) or_return
  1360. }
  1361. /* shift to next bit */
  1362. power >>= 1
  1363. }
  1364. return
  1365. }
  1366. /*
  1367. Calculate `dest = base^power`.
  1368. Assumes `dest` not to be `nil` and to have been initialized.
  1369. */
  1370. internal_int_pow_int :: proc(dest: ^Int, base, power: int, allocator := context.allocator) -> (err: Error) {
  1371. context.allocator = allocator
  1372. base_t := &Int{}
  1373. defer internal_destroy(base_t)
  1374. internal_set(base_t, base) or_return
  1375. return #force_inline internal_int_pow(dest, base_t, power)
  1376. }
  1377. internal_pow :: proc { internal_int_pow, internal_int_pow_int, }
  1378. internal_exp :: pow
  1379. /*
  1380. */
  1381. internal_small_pow :: proc(base: _WORD, exponent: _WORD) -> (result: _WORD) {
  1382. exponent := exponent; base := base
  1383. result = _WORD(1)
  1384. for exponent != 0 {
  1385. if exponent & 1 == 1 {
  1386. result *= base
  1387. }
  1388. exponent >>= 1
  1389. base *= base
  1390. }
  1391. return result
  1392. }
  1393. /*
  1394. This function is less generic than `root_n`, simpler and faster.
  1395. Assumes `dest` and `src` not to be `nil` and to have been initialized.
  1396. */
  1397. internal_int_sqrt :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
  1398. context.allocator = allocator
  1399. /*
  1400. Must be positive.
  1401. */
  1402. if #force_inline internal_is_negative(src) { return .Invalid_Argument }
  1403. /*
  1404. Easy out. If src is zero, so is dest.
  1405. */
  1406. if #force_inline internal_is_zero(src) { return internal_zero(dest) }
  1407. /*
  1408. Set up temporaries.
  1409. */
  1410. x, y, t1, t2 := &Int{}, &Int{}, &Int{}, &Int{}
  1411. defer internal_destroy(x, y, t1, t2)
  1412. count := #force_inline internal_count_bits(src)
  1413. a, b := count >> 1, count & 1
  1414. internal_int_power_of_two(x, a+b, allocator) or_return
  1415. for {
  1416. /*
  1417. y = (x + n // x) // 2
  1418. */
  1419. internal_div(t1, src, x) or_return
  1420. internal_add(t2, t1, x) or_return
  1421. internal_shr(y, t2, 1) or_return
  1422. if internal_gte(y, x) {
  1423. internal_swap(dest, x)
  1424. return nil
  1425. }
  1426. internal_swap(x, y)
  1427. }
  1428. internal_swap(dest, x)
  1429. return err
  1430. }
  1431. internal_sqrt :: proc { internal_int_sqrt, }
  1432. /*
  1433. Find the nth root of an Integer.
  1434. Result found such that `(dest)**n <= src` and `(dest+1)**n > src`
  1435. This algorithm uses Newton's approximation `x[i+1] = x[i] - f(x[i])/f'(x[i])`,
  1436. which will find the root in `log(n)` time where each step involves a fair bit.
  1437. Assumes `dest` and `src` not to be `nil` and have been initialized.
  1438. */
  1439. internal_int_root_n :: proc(dest, src: ^Int, n: int, allocator := context.allocator) -> (err: Error) {
  1440. context.allocator = allocator
  1441. /*
  1442. Fast path for n == 2
  1443. */
  1444. if n == 2 { return #force_inline internal_sqrt(dest, src) }
  1445. if n < 0 || n > int(_DIGIT_MAX) { return .Invalid_Argument }
  1446. if n & 1 == 0 && #force_inline internal_is_negative(src) { return .Invalid_Argument }
  1447. /*
  1448. Set up temporaries.
  1449. */
  1450. t1, t2, t3, a := &Int{}, &Int{}, &Int{}, &Int{}
  1451. defer internal_destroy(t1, t2, t3)
  1452. /*
  1453. If `src` is negative fudge the sign but keep track.
  1454. */
  1455. a.sign = .Zero_or_Positive
  1456. a.used = src.used
  1457. a.digit = src.digit
  1458. /*
  1459. If "n" is larger than INT_MAX it is also larger than
  1460. log_2(src) because the bit-length of the "src" is measured
  1461. with an int and hence the root is always < 2 (two).
  1462. */
  1463. if n > max(int) / 2 {
  1464. err = set(dest, 1)
  1465. dest.sign = a.sign
  1466. return err
  1467. }
  1468. /*
  1469. Compute seed: 2^(log_2(src)/n + 2)
  1470. */
  1471. ilog2 := internal_count_bits(src)
  1472. /*
  1473. "src" is smaller than max(int), we can cast safely.
  1474. */
  1475. if ilog2 < n {
  1476. err = internal_one(dest)
  1477. dest.sign = a.sign
  1478. return err
  1479. }
  1480. ilog2 /= n
  1481. if ilog2 == 0 {
  1482. err = internal_one(dest)
  1483. dest.sign = a.sign
  1484. return err
  1485. }
  1486. /*
  1487. Start value must be larger than root.
  1488. */
  1489. ilog2 += 2
  1490. internal_int_power_of_two(t2, ilog2) or_return
  1491. c: int
  1492. iterations := 0
  1493. for {
  1494. /* t1 = t2 */
  1495. internal_copy(t1, t2) or_return
  1496. /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */
  1497. /* t3 = t1**(b-1) */
  1498. internal_pow(t3, t1, n-1) or_return
  1499. /* numerator */
  1500. /* t2 = t1**b */
  1501. internal_mul(t2, t1, t3) or_return
  1502. /* t2 = t1**b - a */
  1503. internal_sub(t2, t2, a) or_return
  1504. /* denominator */
  1505. /* t3 = t1**(b-1) * b */
  1506. internal_mul(t3, t3, DIGIT(n)) or_return
  1507. /* t3 = (t1**b - a)/(b * t1**(b-1)) */
  1508. internal_div(t3, t2, t3) or_return
  1509. internal_sub(t2, t1, t3) or_return
  1510. /*
  1511. Number of rounds is at most log_2(root). If it is more it
  1512. got stuck, so break out of the loop and do the rest manually.
  1513. */
  1514. if ilog2 -= 1; ilog2 == 0 { break }
  1515. if internal_eq(t1, t2) { break }
  1516. iterations += 1
  1517. if iterations == MAX_ITERATIONS_ROOT_N {
  1518. return .Max_Iterations_Reached
  1519. }
  1520. }
  1521. /* Result can be off by a few so check. */
  1522. /* Loop beneath can overshoot by one if found root is smaller than actual root. */
  1523. iterations = 0
  1524. for {
  1525. internal_pow(t2, t1, n) or_return
  1526. c = internal_cmp(t2, a)
  1527. if c == 0 {
  1528. swap(dest, t1)
  1529. return nil
  1530. } else if c == -1 {
  1531. internal_add(t1, t1, DIGIT(1)) or_return
  1532. } else {
  1533. break
  1534. }
  1535. iterations += 1
  1536. if iterations == MAX_ITERATIONS_ROOT_N {
  1537. return .Max_Iterations_Reached
  1538. }
  1539. }
  1540. iterations = 0
  1541. /*
  1542. Correct overshoot from above or from recurrence.
  1543. */
  1544. for {
  1545. internal_pow(t2, t1, n) or_return
  1546. if internal_lt(t2, a) { break }
  1547. internal_sub(t1, t1, DIGIT(1)) or_return
  1548. iterations += 1
  1549. if iterations == MAX_ITERATIONS_ROOT_N {
  1550. return .Max_Iterations_Reached
  1551. }
  1552. }
  1553. /*
  1554. Set the result.
  1555. */
  1556. internal_swap(dest, t1)
  1557. /*
  1558. Set the sign of the result.
  1559. */
  1560. dest.sign = src.sign
  1561. return err
  1562. }
  1563. internal_root_n :: proc { internal_int_root_n, }
  1564. /*
  1565. Other internal helpers
  1566. */
  1567. /*
  1568. Deallocates the backing memory of one or more `Int`s.
  1569. Asssumes none of the `integers` to be a `nil`.
  1570. */
  1571. internal_int_destroy :: proc(integers: ..^Int) {
  1572. integers := integers
  1573. for a in &integers {
  1574. if internal_int_allocated_cap(a) > 0 {
  1575. mem.zero_slice(a.digit[:])
  1576. free(&a.digit[0])
  1577. }
  1578. a = &Int{}
  1579. }
  1580. }
  1581. internal_destroy :: proc{
  1582. internal_int_destroy,
  1583. internal_rat_destroy,
  1584. }
  1585. /*
  1586. Helpers to set an `Int` to a specific value.
  1587. */
  1588. internal_int_set_from_integer :: proc(dest: ^Int, src: $T, minimize := false, allocator := context.allocator) -> (err: Error)
  1589. where intrinsics.type_is_integer(T) {
  1590. context.allocator = allocator
  1591. internal_error_if_immutable(dest) or_return
  1592. /*
  1593. Most internal procs asssume an Int to have already been initialize,
  1594. but as this is one of the procs that initializes, we have to check the following.
  1595. */
  1596. internal_clear_if_uninitialized_single(dest) or_return
  1597. dest.flags = {} // We're not -Inf, Inf, NaN or Immutable.
  1598. dest.used = 0
  1599. dest.sign = .Negative if src < 0 else .Zero_or_Positive
  1600. temp := src
  1601. is_maximally_negative := src == min(T)
  1602. if is_maximally_negative {
  1603. /*
  1604. Prevent overflow on abs()
  1605. */
  1606. temp += 1
  1607. }
  1608. temp = -temp if temp < 0 else temp
  1609. #no_bounds_check for temp != 0 {
  1610. dest.digit[dest.used] = DIGIT(temp) & _MASK
  1611. dest.used += 1
  1612. temp >>= _DIGIT_BITS
  1613. }
  1614. if is_maximally_negative {
  1615. return internal_sub(dest, dest, 1)
  1616. }
  1617. internal_zero_unused(dest)
  1618. return nil
  1619. }
  1620. internal_set :: proc { internal_int_set_from_integer, internal_int_copy, int_atoi }
  1621. internal_copy_digits :: #force_inline proc(dest, src: ^Int, digits: int, offset := int(0)) -> (err: Error) {
  1622. #force_inline internal_error_if_immutable(dest) or_return
  1623. /*
  1624. If dest == src, do nothing
  1625. */
  1626. return #force_inline _private_copy_digits(dest, src, digits, offset)
  1627. }
  1628. /*
  1629. Copy one `Int` to another.
  1630. */
  1631. internal_int_copy :: proc(dest, src: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
  1632. context.allocator = allocator
  1633. /*
  1634. If dest == src, do nothing
  1635. */
  1636. if (dest == src) { return nil }
  1637. internal_error_if_immutable(dest) or_return
  1638. /*
  1639. Grow `dest` to fit `src`.
  1640. If `dest` is not yet initialized, it will be using `allocator`.
  1641. */
  1642. needed := src.used if minimize else max(src.used, _DEFAULT_DIGIT_COUNT)
  1643. internal_grow(dest, needed, minimize) or_return
  1644. /*
  1645. Copy everything over and zero high digits.
  1646. */
  1647. internal_copy_digits(dest, src, src.used)
  1648. dest.used = src.used
  1649. dest.sign = src.sign
  1650. dest.flags = src.flags &~ {.Immutable}
  1651. internal_zero_unused(dest)
  1652. return nil
  1653. }
  1654. internal_copy :: proc { internal_int_copy, }
  1655. /*
  1656. In normal code, you can also write `a, b = b, a`.
  1657. However, that only swaps within the current scope.
  1658. This helper swaps completely.
  1659. */
  1660. internal_int_swap :: #force_inline proc(a, b: ^Int) {
  1661. a.used, b.used = b.used, a.used
  1662. a.sign, b.sign = b.sign, a.sign
  1663. a.digit, b.digit = b.digit, a.digit
  1664. }
  1665. internal_swap :: proc {
  1666. internal_int_swap,
  1667. internal_rat_swap,
  1668. }
  1669. /*
  1670. Set `dest` to |`src`|.
  1671. */
  1672. internal_int_abs :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
  1673. context.allocator = allocator
  1674. /*
  1675. If `dest == src`, just fix `dest`'s sign.
  1676. */
  1677. if (dest == src) {
  1678. dest.sign = .Zero_or_Positive
  1679. return nil
  1680. }
  1681. /*
  1682. Copy `src` to `dest`
  1683. */
  1684. internal_copy(dest, src) or_return
  1685. /*
  1686. Fix sign.
  1687. */
  1688. dest.sign = .Zero_or_Positive
  1689. return nil
  1690. }
  1691. internal_platform_abs :: proc(n: $T) -> T where intrinsics.type_is_integer(T) {
  1692. return n if n >= 0 else -n
  1693. }
  1694. internal_abs :: proc{ internal_int_abs, internal_platform_abs, }
  1695. /*
  1696. Set `dest` to `-src`.
  1697. */
  1698. internal_int_neg :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
  1699. context.allocator = allocator
  1700. /*
  1701. If `dest == src`, just fix `dest`'s sign.
  1702. */
  1703. sign := Sign.Negative
  1704. if #force_inline internal_is_zero(src) || #force_inline internal_is_negative(src) {
  1705. sign = .Zero_or_Positive
  1706. }
  1707. if dest == src {
  1708. dest.sign = sign
  1709. return nil
  1710. }
  1711. /*
  1712. Copy `src` to `dest`
  1713. */
  1714. internal_copy(dest, src) or_return
  1715. /*
  1716. Fix sign.
  1717. */
  1718. dest.sign = sign
  1719. return nil
  1720. }
  1721. internal_neg :: proc { internal_int_neg, }
  1722. /*
  1723. hac 14.61, pp608.
  1724. */
  1725. internal_int_inverse_modulo :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
  1726. context.allocator = allocator
  1727. /*
  1728. For all n in N and n > 0, n = 0 mod 1.
  1729. */
  1730. if internal_is_positive(a) && internal_eq(b, 1) { return internal_zero(dest) }
  1731. /*
  1732. `b` cannot be negative and has to be > 1
  1733. */
  1734. if internal_is_negative(b) || internal_gt(b, 1) { return .Invalid_Argument }
  1735. /*
  1736. If the modulus is odd we can use a faster routine instead.
  1737. */
  1738. if internal_is_odd(b) { return _private_inverse_modulo_odd(dest, a, b) }
  1739. return _private_inverse_modulo(dest, a, b)
  1740. }
  1741. internal_invmod :: proc{ internal_int_inverse_modulo, }
  1742. /*
  1743. Helpers to extract values from the `Int`.
  1744. Offset is zero indexed.
  1745. */
  1746. internal_int_bitfield_extract_bool :: proc(a: ^Int, offset: int) -> (val: bool, err: Error) {
  1747. limb := offset / _DIGIT_BITS
  1748. if limb < 0 || limb >= a.used { return false, .Invalid_Argument }
  1749. i := _WORD(1 << _WORD((offset % _DIGIT_BITS)))
  1750. return bool(_WORD(a.digit[limb]) & i), nil
  1751. }
  1752. internal_int_bitfield_extract_single :: proc(a: ^Int, offset: int) -> (bit: _WORD, err: Error) {
  1753. limb := offset / _DIGIT_BITS
  1754. if limb < 0 || limb >= a.used { return 0, .Invalid_Argument }
  1755. i := _WORD(1 << _WORD((offset % _DIGIT_BITS)))
  1756. return 1 if ((_WORD(a.digit[limb]) & i) != 0) else 0, nil
  1757. }
  1758. internal_int_bitfield_extract :: proc(a: ^Int, offset, count: int) -> (res: _WORD, err: Error) #no_bounds_check {
  1759. /*
  1760. Early out for single bit.
  1761. */
  1762. if count == 1 {
  1763. limb := offset / _DIGIT_BITS
  1764. if limb < 0 || limb >= a.used { return 0, .Invalid_Argument }
  1765. i := _WORD(1 << _WORD((offset % _DIGIT_BITS)))
  1766. return 1 if ((_WORD(a.digit[limb]) & i) != 0) else 0, nil
  1767. }
  1768. if count > _WORD_BITS || count < 1 { return 0, .Invalid_Argument }
  1769. /*
  1770. There are 3 possible cases.
  1771. - [offset:][:count] covers 1 DIGIT,
  1772. e.g. offset: 0, count: 60 = bits 0..59
  1773. - [offset:][:count] covers 2 DIGITS,
  1774. e.g. offset: 5, count: 60 = bits 5..59, 0..4
  1775. e.g. offset: 0, count: 120 = bits 0..59, 60..119
  1776. - [offset:][:count] covers 3 DIGITS,
  1777. e.g. offset: 40, count: 100 = bits 40..59, 0..59, 0..19
  1778. e.g. offset: 40, count: 120 = bits 40..59, 0..59, 0..39
  1779. */
  1780. limb := offset / _DIGIT_BITS
  1781. bits_left := count
  1782. bits_offset := offset % _DIGIT_BITS
  1783. num_bits := min(bits_left, _DIGIT_BITS - bits_offset)
  1784. shift := offset % _DIGIT_BITS
  1785. mask := (_WORD(1) << uint(num_bits)) - 1
  1786. res = (_WORD(a.digit[limb]) >> uint(shift)) & mask
  1787. bits_left -= num_bits
  1788. if bits_left == 0 { return res, nil }
  1789. res_shift := num_bits
  1790. num_bits = min(bits_left, _DIGIT_BITS)
  1791. mask = (1 << uint(num_bits)) - 1
  1792. res |= (_WORD(a.digit[limb + 1]) & mask) << uint(res_shift)
  1793. bits_left -= num_bits
  1794. if bits_left == 0 { return res, nil }
  1795. mask = (1 << uint(bits_left)) - 1
  1796. res_shift += _DIGIT_BITS
  1797. res |= (_WORD(a.digit[limb + 2]) & mask) << uint(res_shift)
  1798. return res, nil
  1799. }
  1800. /*
  1801. Helpers to (un)set a bit in an Int.
  1802. Offset is zero indexed.
  1803. */
  1804. internal_int_bitfield_set_single :: proc(a: ^Int, offset: int) -> (err: Error) {
  1805. limb := offset / _DIGIT_BITS
  1806. if limb < 0 || limb >= a.used { return .Invalid_Argument }
  1807. i := DIGIT(1 << uint((offset % _DIGIT_BITS)))
  1808. a.digit[limb] |= i
  1809. return
  1810. }
  1811. internal_int_bitfield_unset_single :: proc(a: ^Int, offset: int) -> (err: Error) {
  1812. limb := offset / _DIGIT_BITS
  1813. if limb < 0 || limb >= a.used { return .Invalid_Argument }
  1814. i := DIGIT(1 << uint((offset % _DIGIT_BITS)))
  1815. a.digit[limb] &= _MASK - i
  1816. return
  1817. }
  1818. internal_int_bitfield_toggle_single :: proc(a: ^Int, offset: int) -> (err: Error) {
  1819. limb := offset / _DIGIT_BITS
  1820. if limb < 0 || limb >= a.used { return .Invalid_Argument }
  1821. i := DIGIT(1 << uint((offset % _DIGIT_BITS)))
  1822. a.digit[limb] ~= i
  1823. return
  1824. }
  1825. /*
  1826. Resize backing store.
  1827. We don't need to pass the allocator, because the storage itself stores it.
  1828. Assumes `a` not to be `nil`, and to have already been initialized.
  1829. */
  1830. internal_int_shrink :: proc(a: ^Int) -> (err: Error) {
  1831. needed := max(_MIN_DIGIT_COUNT, a.used)
  1832. if a.used != needed { return internal_grow(a, needed, true) }
  1833. return nil
  1834. }
  1835. internal_shrink :: proc { internal_int_shrink, }
  1836. internal_int_grow :: proc(a: ^Int, digits: int, allow_shrink := false, allocator := context.allocator) -> (err: Error) {
  1837. /*
  1838. We need at least _MIN_DIGIT_COUNT or a.used digits, whichever is bigger.
  1839. The caller is asking for `digits`. Let's be accomodating.
  1840. */
  1841. cap := internal_int_allocated_cap(a)
  1842. needed := max(_MIN_DIGIT_COUNT, a.used, digits)
  1843. if !allow_shrink {
  1844. needed = max(needed, cap)
  1845. }
  1846. /*
  1847. If not yet iniialized, initialize the `digit` backing with the allocator we were passed.
  1848. */
  1849. if cap == 0 {
  1850. a.digit = make([dynamic]DIGIT, needed, allocator)
  1851. } else if cap != needed {
  1852. /*
  1853. `[dynamic]DIGIT` already knows what allocator was used for it, so resize will do the right thing.
  1854. */
  1855. resize(&a.digit, needed)
  1856. }
  1857. /*
  1858. Let's see if the allocation/resize worked as expected.
  1859. */
  1860. if len(a.digit) != needed {
  1861. return .Out_Of_Memory
  1862. }
  1863. return nil
  1864. }
  1865. internal_grow :: proc { internal_int_grow, }
  1866. /*
  1867. Clear `Int` and resize it to the default size.
  1868. Assumes `a` not to be `nil`.
  1869. */
  1870. internal_int_clear :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
  1871. raw := transmute(mem.Raw_Dynamic_Array)a.digit
  1872. if raw.cap != 0 {
  1873. mem.zero_slice(a.digit[:a.used])
  1874. }
  1875. a.sign = .Zero_or_Positive
  1876. a.used = 0
  1877. return #force_inline internal_grow(a, a.used, minimize, allocator)
  1878. }
  1879. internal_clear :: proc { internal_int_clear, }
  1880. internal_zero :: internal_clear
  1881. /*
  1882. Set the `Int` to 1 and optionally shrink it to the minimum backing size.
  1883. */
  1884. internal_int_one :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
  1885. return internal_copy(a, INT_ONE, minimize, allocator)
  1886. }
  1887. internal_one :: proc { internal_int_one, }
  1888. /*
  1889. Set the `Int` to -1 and optionally shrink it to the minimum backing size.
  1890. */
  1891. internal_int_minus_one :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
  1892. return internal_copy(a, INT_MINUS_ONE, minimize, allocator)
  1893. }
  1894. internal_minus_one :: proc { internal_int_minus_one, }
  1895. /*
  1896. Set the `Int` to Inf and optionally shrink it to the minimum backing size.
  1897. */
  1898. internal_int_inf :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
  1899. return internal_copy(a, INT_INF, minimize, allocator)
  1900. }
  1901. internal_inf :: proc { internal_int_inf, }
  1902. /*
  1903. Set the `Int` to -Inf and optionally shrink it to the minimum backing size.
  1904. */
  1905. internal_int_minus_inf :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
  1906. return internal_copy(a, INT_MINUS_INF, minimize, allocator)
  1907. }
  1908. internal_minus_inf :: proc { internal_int_inf, }
  1909. /*
  1910. Set the `Int` to NaN and optionally shrink it to the minimum backing size.
  1911. */
  1912. internal_int_nan :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
  1913. return internal_copy(a, INT_NAN, minimize, allocator)
  1914. }
  1915. internal_nan :: proc { internal_int_nan, }
  1916. internal_int_power_of_two :: proc(a: ^Int, power: int, allocator := context.allocator) -> (err: Error) {
  1917. context.allocator = allocator
  1918. if power < 0 || power > _MAX_BIT_COUNT { return .Invalid_Argument }
  1919. /*
  1920. Grow to accomodate the single bit.
  1921. */
  1922. a.used = (power / _DIGIT_BITS) + 1
  1923. internal_grow(a, a.used) or_return
  1924. /*
  1925. Zero the entirety.
  1926. */
  1927. mem.zero_slice(a.digit[:])
  1928. /*
  1929. Set the bit.
  1930. */
  1931. a.digit[power / _DIGIT_BITS] = 1 << uint((power % _DIGIT_BITS))
  1932. return nil
  1933. }
  1934. internal_int_get_u128 :: proc(a: ^Int) -> (res: u128, err: Error) {
  1935. return internal_int_get(a, u128)
  1936. }
  1937. internal_get_u128 :: proc { internal_int_get_u128, }
  1938. internal_int_get_i128 :: proc(a: ^Int) -> (res: i128, err: Error) {
  1939. return internal_int_get(a, i128)
  1940. }
  1941. internal_get_i128 :: proc { internal_int_get_i128, }
  1942. internal_int_get_u64 :: proc(a: ^Int) -> (res: u64, err: Error) {
  1943. return internal_int_get(a, u64)
  1944. }
  1945. internal_get_u64 :: proc { internal_int_get_u64, }
  1946. internal_int_get_i64 :: proc(a: ^Int) -> (res: i64, err: Error) {
  1947. return internal_int_get(a, i64)
  1948. }
  1949. internal_get_i64 :: proc { internal_int_get_i64, }
  1950. internal_int_get_u32 :: proc(a: ^Int) -> (res: u32, err: Error) {
  1951. return internal_int_get(a, u32)
  1952. }
  1953. internal_get_u32 :: proc { internal_int_get_u32, }
  1954. internal_int_get_i32 :: proc(a: ^Int) -> (res: i32, err: Error) {
  1955. return internal_int_get(a, i32)
  1956. }
  1957. internal_get_i32 :: proc { internal_int_get_i32, }
  1958. internal_get_low_u32 :: proc(a: ^Int) -> u32 #no_bounds_check {
  1959. if a == nil {
  1960. return 0
  1961. }
  1962. if a.used == 0 {
  1963. return 0
  1964. }
  1965. return u32(a.digit[0])
  1966. }
  1967. internal_get_low_u64 :: proc(a: ^Int) -> u64 #no_bounds_check {
  1968. if a == nil {
  1969. return 0
  1970. }
  1971. if a.used == 0 {
  1972. return 0
  1973. }
  1974. v := u64(a.digit[0])
  1975. when size_of(DIGIT) == 4 {
  1976. if a.used > 1 {
  1977. return u64(a.digit[1])<<32 | v
  1978. }
  1979. }
  1980. return v
  1981. }
  1982. /*
  1983. TODO: Think about using `count_bits` to check if the value could be returned completely,
  1984. and maybe return max(T), .Integer_Overflow if not?
  1985. */
  1986. internal_int_get :: proc(a: ^Int, $T: typeid) -> (res: T, err: Error) where intrinsics.type_is_integer(T) {
  1987. /*
  1988. Calculate target bit size.
  1989. */
  1990. target_bit_size := int(size_of(T) * 8)
  1991. when !intrinsics.type_is_unsigned(T) {
  1992. if a.sign == .Zero_or_Positive {
  1993. target_bit_size -= 1
  1994. }
  1995. } else {
  1996. if a.sign == .Negative {
  1997. return 0, .Integer_Underflow
  1998. }
  1999. }
  2000. bits_used := internal_count_bits(a)
  2001. if bits_used > target_bit_size {
  2002. if a.sign == .Negative {
  2003. return min(T), .Integer_Underflow
  2004. }
  2005. return max(T), .Integer_Overflow
  2006. }
  2007. for i := a.used; i > 0; i -= 1 {
  2008. res <<= _DIGIT_BITS
  2009. res |= T(a.digit[i - 1])
  2010. }
  2011. when !intrinsics.type_is_unsigned(T) {
  2012. /*
  2013. Set the sign.
  2014. */
  2015. if a.sign == .Negative { res = -res }
  2016. }
  2017. return
  2018. }
  2019. internal_get :: proc { internal_int_get, }
  2020. internal_int_get_float :: proc(a: ^Int) -> (res: f64, err: Error) {
  2021. /*
  2022. log2(max(f64)) is approximately 1020, or 17 legs with the 64-bit storage.
  2023. */
  2024. legs :: 1020 / _DIGIT_BITS
  2025. l := min(a.used, legs)
  2026. fac := f64(1 << _DIGIT_BITS)
  2027. d := 0.0
  2028. #no_bounds_check for i := l; i >= 0; i -= 1 {
  2029. d = (d * fac) + f64(a.digit[i])
  2030. }
  2031. res = -d if a.sign == .Negative else d
  2032. return
  2033. }
  2034. /*
  2035. The `and`, `or` and `xor` binops differ in two lines only.
  2036. We could handle those with a switch, but that adds overhead.
  2037. TODO: Implement versions that take a DIGIT immediate.
  2038. */
  2039. /*
  2040. 2's complement `and`, returns `dest = a & b;`
  2041. */
  2042. internal_int_and :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
  2043. context.allocator = allocator
  2044. used := max(a.used, b.used) + 1
  2045. /*
  2046. Grow the destination to accomodate the result.
  2047. */
  2048. internal_grow(dest, used) or_return
  2049. neg_a := #force_inline internal_is_negative(a)
  2050. neg_b := #force_inline internal_is_negative(b)
  2051. neg := neg_a && neg_b
  2052. ac, bc, cc := DIGIT(1), DIGIT(1), DIGIT(1)
  2053. #no_bounds_check for i := 0; i < used; i += 1 {
  2054. x, y: DIGIT
  2055. /*
  2056. Convert to 2's complement if negative.
  2057. */
  2058. if neg_a {
  2059. ac += _MASK if i >= a.used else (~a.digit[i] & _MASK)
  2060. x = ac & _MASK
  2061. ac >>= _DIGIT_BITS
  2062. } else {
  2063. x = 0 if i >= a.used else a.digit[i]
  2064. }
  2065. /*
  2066. Convert to 2's complement if negative.
  2067. */
  2068. if neg_b {
  2069. bc += _MASK if i >= b.used else (~b.digit[i] & _MASK)
  2070. y = bc & _MASK
  2071. bc >>= _DIGIT_BITS
  2072. } else {
  2073. y = 0 if i >= b.used else b.digit[i]
  2074. }
  2075. dest.digit[i] = x & y
  2076. /*
  2077. Convert to to sign-magnitude if negative.
  2078. */
  2079. if neg {
  2080. cc += ~dest.digit[i] & _MASK
  2081. dest.digit[i] = cc & _MASK
  2082. cc >>= _DIGIT_BITS
  2083. }
  2084. }
  2085. dest.used = used
  2086. dest.sign = .Negative if neg else .Zero_or_Positive
  2087. return internal_clamp(dest)
  2088. }
  2089. internal_and :: proc { internal_int_and, }
  2090. /*
  2091. 2's complement `or`, returns `dest = a | b;`
  2092. */
  2093. internal_int_or :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
  2094. context.allocator = allocator
  2095. used := max(a.used, b.used) + 1
  2096. /*
  2097. Grow the destination to accomodate the result.
  2098. */
  2099. internal_grow(dest, used) or_return
  2100. neg_a := #force_inline internal_is_negative(a)
  2101. neg_b := #force_inline internal_is_negative(b)
  2102. neg := neg_a || neg_b
  2103. ac, bc, cc := DIGIT(1), DIGIT(1), DIGIT(1)
  2104. #no_bounds_check for i := 0; i < used; i += 1 {
  2105. x, y: DIGIT
  2106. /*
  2107. Convert to 2's complement if negative.
  2108. */
  2109. if neg_a {
  2110. ac += _MASK if i >= a.used else (~a.digit[i] & _MASK)
  2111. x = ac & _MASK
  2112. ac >>= _DIGIT_BITS
  2113. } else {
  2114. x = 0 if i >= a.used else a.digit[i]
  2115. }
  2116. /*
  2117. Convert to 2's complement if negative.
  2118. */
  2119. if neg_b {
  2120. bc += _MASK if i >= b.used else (~b.digit[i] & _MASK)
  2121. y = bc & _MASK
  2122. bc >>= _DIGIT_BITS
  2123. } else {
  2124. y = 0 if i >= b.used else b.digit[i]
  2125. }
  2126. dest.digit[i] = x | y
  2127. /*
  2128. Convert to to sign-magnitude if negative.
  2129. */
  2130. if neg {
  2131. cc += ~dest.digit[i] & _MASK
  2132. dest.digit[i] = cc & _MASK
  2133. cc >>= _DIGIT_BITS
  2134. }
  2135. }
  2136. dest.used = used
  2137. dest.sign = .Negative if neg else .Zero_or_Positive
  2138. return internal_clamp(dest)
  2139. }
  2140. internal_or :: proc { internal_int_or, }
  2141. /*
  2142. 2's complement `xor`, returns `dest = a ~ b;`
  2143. */
  2144. internal_int_xor :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
  2145. context.allocator = allocator
  2146. used := max(a.used, b.used) + 1
  2147. /*
  2148. Grow the destination to accomodate the result.
  2149. */
  2150. internal_grow(dest, used) or_return
  2151. neg_a := #force_inline internal_is_negative(a)
  2152. neg_b := #force_inline internal_is_negative(b)
  2153. neg := neg_a != neg_b
  2154. ac, bc, cc := DIGIT(1), DIGIT(1), DIGIT(1)
  2155. #no_bounds_check for i := 0; i < used; i += 1 {
  2156. x, y: DIGIT
  2157. /*
  2158. Convert to 2's complement if negative.
  2159. */
  2160. if neg_a {
  2161. ac += _MASK if i >= a.used else (~a.digit[i] & _MASK)
  2162. x = ac & _MASK
  2163. ac >>= _DIGIT_BITS
  2164. } else {
  2165. x = 0 if i >= a.used else a.digit[i]
  2166. }
  2167. /*
  2168. Convert to 2's complement if negative.
  2169. */
  2170. if neg_b {
  2171. bc += _MASK if i >= b.used else (~b.digit[i] & _MASK)
  2172. y = bc & _MASK
  2173. bc >>= _DIGIT_BITS
  2174. } else {
  2175. y = 0 if i >= b.used else b.digit[i]
  2176. }
  2177. dest.digit[i] = x ~ y
  2178. /*
  2179. Convert to to sign-magnitude if negative.
  2180. */
  2181. if neg {
  2182. cc += ~dest.digit[i] & _MASK
  2183. dest.digit[i] = cc & _MASK
  2184. cc >>= _DIGIT_BITS
  2185. }
  2186. }
  2187. dest.used = used
  2188. dest.sign = .Negative if neg else .Zero_or_Positive
  2189. return internal_clamp(dest)
  2190. }
  2191. internal_xor :: proc { internal_int_xor, }
  2192. /*
  2193. dest = ~src
  2194. */
  2195. internal_int_complement :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
  2196. context.allocator = allocator
  2197. /*
  2198. Temporarily fix sign.
  2199. */
  2200. old_sign := src.sign
  2201. neg := #force_inline internal_is_zero(src) || #force_inline internal_is_positive(src)
  2202. src.sign = .Negative if neg else .Zero_or_Positive
  2203. err = #force_inline internal_sub(dest, src, 1)
  2204. /*
  2205. Restore sign.
  2206. */
  2207. src.sign = old_sign
  2208. return err
  2209. }
  2210. internal_complement :: proc { internal_int_complement, }
  2211. /*
  2212. quotient, remainder := numerator >> bits;
  2213. `remainder` is allowed to be passed a `nil`, in which case `mod` won't be computed.
  2214. */
  2215. internal_int_shrmod :: proc(quotient, remainder, numerator: ^Int, bits: int, allocator := context.allocator) -> (err: Error) {
  2216. context.allocator = allocator
  2217. bits := bits
  2218. if bits < 0 { return .Invalid_Argument }
  2219. internal_copy(quotient, numerator) or_return
  2220. /*
  2221. Shift right by a certain bit count (store quotient and optional remainder.)
  2222. `numerator` should not be used after this.
  2223. */
  2224. if remainder != nil {
  2225. internal_int_mod_bits(remainder, numerator, bits) or_return
  2226. }
  2227. /*
  2228. Shift by as many digits in the bit count.
  2229. */
  2230. if bits >= _DIGIT_BITS {
  2231. _private_int_shr_leg(quotient, bits / _DIGIT_BITS) or_return
  2232. }
  2233. /*
  2234. Shift any bit count < _DIGIT_BITS.
  2235. */
  2236. bits %= _DIGIT_BITS
  2237. if bits != 0 {
  2238. mask := DIGIT(1 << uint(bits)) - 1
  2239. shift := DIGIT(_DIGIT_BITS - bits)
  2240. carry := DIGIT(0)
  2241. #no_bounds_check for x := quotient.used - 1; x >= 0; x -= 1 {
  2242. /*
  2243. Get the lower bits of this word in a temp.
  2244. */
  2245. fwd_carry := quotient.digit[x] & mask
  2246. /*
  2247. Shift the current word and mix in the carry bits from the previous word.
  2248. */
  2249. quotient.digit[x] = (quotient.digit[x] >> uint(bits)) | (carry << shift)
  2250. /*
  2251. Update carry from forward carry.
  2252. */
  2253. carry = fwd_carry
  2254. }
  2255. }
  2256. return internal_clamp(numerator)
  2257. }
  2258. internal_shrmod :: proc { internal_int_shrmod, }
  2259. internal_int_shr :: proc(dest, source: ^Int, bits: int, allocator := context.allocator) -> (err: Error) {
  2260. return #force_inline internal_shrmod(dest, nil, source, bits, allocator)
  2261. }
  2262. internal_shr :: proc { internal_int_shr, }
  2263. /*
  2264. Shift right by a certain bit count with sign extension.
  2265. */
  2266. internal_int_shr_signed :: proc(dest, src: ^Int, bits: int, allocator := context.allocator) -> (err: Error) {
  2267. context.allocator = allocator
  2268. if src.sign == .Zero_or_Positive {
  2269. return internal_shr(dest, src, bits)
  2270. }
  2271. internal_int_add_digit(dest, src, DIGIT(1)) or_return
  2272. internal_shr(dest, dest, bits) or_return
  2273. return internal_sub(dest, src, DIGIT(1))
  2274. }
  2275. internal_shr_signed :: proc { internal_int_shr_signed, }
  2276. /*
  2277. Shift left by a certain bit count.
  2278. */
  2279. internal_int_shl :: proc(dest, src: ^Int, bits: int, allocator := context.allocator) -> (err: Error) {
  2280. context.allocator = allocator
  2281. bits := bits
  2282. if bits < 0 { return .Invalid_Argument }
  2283. internal_copy(dest, src) or_return
  2284. /*
  2285. Grow `dest` to accommodate the additional bits.
  2286. */
  2287. digits_needed := dest.used + (bits / _DIGIT_BITS) + 1
  2288. internal_grow(dest, digits_needed) or_return
  2289. dest.used = digits_needed
  2290. /*
  2291. Shift by as many digits in the bit count as we have.
  2292. */
  2293. if bits >= _DIGIT_BITS {
  2294. _private_int_shl_leg(dest, bits / _DIGIT_BITS) or_return
  2295. }
  2296. /*
  2297. Shift any remaining bit count < _DIGIT_BITS
  2298. */
  2299. bits %= _DIGIT_BITS
  2300. if bits != 0 {
  2301. mask := (DIGIT(1) << uint(bits)) - DIGIT(1)
  2302. shift := DIGIT(_DIGIT_BITS - bits)
  2303. carry := DIGIT(0)
  2304. #no_bounds_check for x:= 0; x < dest.used; x+= 1 {
  2305. fwd_carry := (dest.digit[x] >> shift) & mask
  2306. dest.digit[x] = (dest.digit[x] << uint(bits) | carry) & _MASK
  2307. carry = fwd_carry
  2308. }
  2309. /*
  2310. Use final carry.
  2311. */
  2312. if carry != 0 {
  2313. dest.digit[dest.used] = carry
  2314. dest.used += 1
  2315. }
  2316. }
  2317. return internal_clamp(dest)
  2318. }
  2319. internal_shl :: proc { internal_int_shl, }
  2320. /*
  2321. Count bits in an `Int`.
  2322. Assumes `a` not to be `nil` and to have been initialized.
  2323. */
  2324. internal_count_bits :: proc(a: ^Int) -> (count: int) {
  2325. /*
  2326. Fast path for zero.
  2327. */
  2328. if #force_inline internal_is_zero(a) { return {} }
  2329. /*
  2330. Get the number of DIGITs and use it.
  2331. */
  2332. count = (a.used - 1) * _DIGIT_BITS
  2333. /*
  2334. Take the last DIGIT and count the bits in it.
  2335. */
  2336. clz := int(intrinsics.count_leading_zeros(a.digit[a.used - 1]))
  2337. count += (_DIGIT_TYPE_BITS - clz)
  2338. return
  2339. }
  2340. /*
  2341. Returns the number of trailing zeroes before the first one.
  2342. Differs from regular `ctz` in that 0 returns 0.
  2343. Assumes `a` not to be `nil` and have been initialized.
  2344. */
  2345. internal_int_count_lsb :: proc(a: ^Int) -> (count: int, err: Error) {
  2346. /*
  2347. Easy out.
  2348. */
  2349. if #force_inline internal_is_zero(a) { return {}, nil }
  2350. /*
  2351. Scan lower digits until non-zero.
  2352. */
  2353. x: int
  2354. #no_bounds_check for x = 0; x < a.used && a.digit[x] == 0; x += 1 {}
  2355. when true {
  2356. q := a.digit[x]
  2357. x *= _DIGIT_BITS
  2358. x += internal_count_lsb(q)
  2359. } else {
  2360. lnz := []int{
  2361. 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
  2362. }
  2363. q := a.digit[x]
  2364. x *= _DIGIT_BITS
  2365. if q & 1 == 0 {
  2366. p: DIGIT
  2367. for {
  2368. p = q & 15
  2369. x += lnz[p]
  2370. q >>= 4
  2371. if p != 0 { break }
  2372. }
  2373. }
  2374. }
  2375. return x, nil
  2376. }
  2377. internal_platform_count_lsb :: #force_inline proc(a: $T) -> (count: int)
  2378. where intrinsics.type_is_integer(T) && intrinsics.type_is_unsigned(T) {
  2379. return int(intrinsics.count_trailing_zeros(a)) if a > 0 else 0
  2380. }
  2381. internal_count_lsb :: proc { internal_int_count_lsb, internal_platform_count_lsb, }
  2382. internal_int_random_digit :: proc(r: ^rnd.Rand = nil) -> (res: DIGIT) {
  2383. when _DIGIT_BITS == 60 { // DIGIT = u64
  2384. return DIGIT(rnd.uint64(r)) & _MASK
  2385. } else when _DIGIT_BITS == 28 { // DIGIT = u32
  2386. return DIGIT(rnd.uint32(r)) & _MASK
  2387. } else {
  2388. panic("Unsupported DIGIT size.")
  2389. }
  2390. return 0 // We shouldn't get here.
  2391. }
  2392. internal_int_random :: proc(dest: ^Int, bits: int, r: ^rnd.Rand = nil, allocator := context.allocator) -> (err: Error) {
  2393. context.allocator = allocator
  2394. bits := bits
  2395. if bits <= 0 { return .Invalid_Argument }
  2396. digits := bits / _DIGIT_BITS
  2397. bits %= _DIGIT_BITS
  2398. if bits > 0 {
  2399. digits += 1
  2400. }
  2401. #force_inline internal_grow(dest, digits) or_return
  2402. for i := 0; i < digits; i += 1 {
  2403. dest.digit[i] = int_random_digit(r) & _MASK
  2404. }
  2405. if bits > 0 {
  2406. dest.digit[digits - 1] &= ((1 << uint(bits)) - 1)
  2407. }
  2408. dest.used = digits
  2409. return nil
  2410. }
  2411. internal_random :: proc { internal_int_random, }
  2412. /*
  2413. Internal helpers.
  2414. */
  2415. internal_assert_initialized :: proc(a: ^Int, loc := #caller_location) {
  2416. assert(internal_is_initialized(a), "`Int` was not properly initialized.", loc)
  2417. }
  2418. internal_clear_if_uninitialized_single :: proc(arg: ^Int, allocator := context.allocator) -> (err: Error) {
  2419. context.allocator = allocator
  2420. if ! #force_inline internal_is_initialized(arg) {
  2421. return #force_inline internal_grow(arg, _DEFAULT_DIGIT_COUNT)
  2422. }
  2423. return err
  2424. }
  2425. internal_clear_if_uninitialized_multi :: proc(args: ..^Int, allocator := context.allocator) -> (err: Error) {
  2426. context.allocator = allocator
  2427. for i in args {
  2428. if ! #force_inline internal_is_initialized(i) {
  2429. e := #force_inline internal_grow(i, _DEFAULT_DIGIT_COUNT)
  2430. if e != nil { err = e }
  2431. }
  2432. }
  2433. return err
  2434. }
  2435. internal_clear_if_uninitialized :: proc {internal_clear_if_uninitialized_single, internal_clear_if_uninitialized_multi, }
  2436. internal_error_if_immutable_single :: proc(arg: ^Int) -> (err: Error) {
  2437. if arg != nil && .Immutable in arg.flags { return .Assignment_To_Immutable }
  2438. return nil
  2439. }
  2440. internal_error_if_immutable_multi :: proc(args: ..^Int) -> (err: Error) {
  2441. for i in args {
  2442. if i != nil && .Immutable in i.flags { return .Assignment_To_Immutable }
  2443. }
  2444. return nil
  2445. }
  2446. internal_error_if_immutable :: proc {internal_error_if_immutable_single, internal_error_if_immutable_multi, }
  2447. /*
  2448. Allocates several `Int`s at once.
  2449. */
  2450. internal_int_init_multi :: proc(integers: ..^Int, allocator := context.allocator) -> (err: Error) {
  2451. context.allocator = allocator
  2452. integers := integers
  2453. for a in &integers {
  2454. internal_clear(a) or_return
  2455. }
  2456. return nil
  2457. }
  2458. internal_init_multi :: proc { internal_int_init_multi, }
  2459. /*
  2460. Trim unused digits.
  2461. This is used to ensure that leading zero digits are trimmed and the leading "used" digit will be non-zero.
  2462. Typically very fast. Also fixes the sign if there are no more leading digits.
  2463. */
  2464. internal_clamp :: proc(a: ^Int) -> (err: Error) {
  2465. for a.used > 0 && a.digit[a.used - 1] == 0 { a.used -= 1 }
  2466. if #force_inline internal_is_zero(a) { a.sign = .Zero_or_Positive }
  2467. return nil
  2468. }
  2469. internal_int_zero_unused :: #force_inline proc(dest: ^Int, old_used := -1) {
  2470. /*
  2471. If we don't pass the number of previously used DIGITs, we zero all remaining ones.
  2472. */
  2473. zero_count: int
  2474. if old_used == -1 {
  2475. zero_count = len(dest.digit) - dest.used
  2476. } else {
  2477. zero_count = old_used - dest.used
  2478. }
  2479. /*
  2480. Zero remainder.
  2481. */
  2482. if zero_count > 0 && dest.used < len(dest.digit) {
  2483. mem.zero_slice(dest.digit[dest.used:][:zero_count])
  2484. }
  2485. }
  2486. internal_zero_unused :: proc { internal_int_zero_unused, }
  2487. /*
  2488. ========================== End of low-level routines ==========================
  2489. */