internal.odin 72 KB

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