internal.odin 74 KB

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