internal.odin 73 KB

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