123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853285428552856285728582859286028612862286328642865286628672868286928702871287228732874287528762877287828792880288128822883288428852886288728882889289028912892289328942895289628972898289929002901290229032904290529062907290829092910291129122913291429152916291729182919292029212922292329242925292629272928292929302931293229332934293529362937293829392940294129422943294429452946294729482949295029512952295329542955295629572958295929602961296229632964296529662967296829692970 |
- //+ignore
- /*
- Copyright 2021 Jeroen van Rijn <[email protected]>.
- Made available under Odin's BSD-3 license.
- A BigInt implementation in Odin.
- For the theoretical underpinnings, see Knuth's The Art of Computer Programming, Volume 2, section 4.3.
- The code started out as an idiomatic source port of libTomMath, which is in the public domain, with thanks.
- ========================== Low-level routines ==========================
- IMPORTANT: `internal_*` procedures make certain assumptions about their input.
- The public functions that call them are expected to satisfy their sanity check requirements.
- This allows `internal_*` call `internal_*` without paying this overhead multiple times.
- Where errors can occur, they are of course still checked and returned as appropriate.
- When importing `math:core/big` to implement an involved algorithm of your own, you are welcome
- to use these procedures instead of their public counterparts.
- Most inputs and outputs are expected to be passed an initialized `Int`, for example.
- Exceptions include `quotient` and `remainder`, which are allowed to be `nil` when the calling code doesn't need them.
- Check the comments above each `internal_*` implementation to see what constraints it expects to have met.
- We pass the custom allocator to procedures by default using the pattern `context.allocator = allocator`.
- This way we don't have to add `, allocator` at the end of each call.
- TODO: Handle +/- Infinity and NaN.
- */
- package math_big
- import "core:mem"
- import "core:intrinsics"
- import rnd "core:math/rand"
- /*
- Low-level addition, unsigned. Handbook of Applied Cryptography, algorithm 14.7.
- Assumptions:
- `dest`, `a` and `b` != `nil` and have been initalized.
- */
- internal_int_add_unsigned :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- dest := dest; x := a; y := b
- context.allocator = allocator
- old_used, min_used, max_used, i: int
- if x.used < y.used {
- x, y = y, x
- }
- min_used = y.used
- max_used = x.used
- old_used = dest.used
- internal_grow(dest, max(max_used + 1, _DEFAULT_DIGIT_COUNT)) or_return
- dest.used = max_used + 1
- /*
- All parameters have been initialized.
- */
- /* Zero the carry */
- carry := DIGIT(0)
- #no_bounds_check for i = 0; i < min_used; i += 1 {
- /*
- Compute the sum one _DIGIT at a time.
- dest[i] = a[i] + b[i] + carry;
- */
- dest.digit[i] = x.digit[i] + y.digit[i] + carry
- /*
- Compute carry
- */
- carry = dest.digit[i] >> _DIGIT_BITS
- /*
- Mask away carry from result digit.
- */
- dest.digit[i] &= _MASK
- }
- if min_used != max_used {
- /*
- Now copy higher words, if any, in A+B.
- If A or B has more digits, add those in.
- */
- #no_bounds_check for ; i < max_used; i += 1 {
- dest.digit[i] = x.digit[i] + carry
- /*
- Compute carry
- */
- carry = dest.digit[i] >> _DIGIT_BITS
- /*
- Mask away carry from result digit.
- */
- dest.digit[i] &= _MASK
- }
- }
- /*
- Add remaining carry.
- */
- dest.digit[i] = carry
- /*
- Zero remainder.
- */
- internal_zero_unused(dest, old_used)
- /*
- Adjust dest.used based on leading zeroes.
- */
- return internal_clamp(dest)
- }
- internal_add_unsigned :: proc { internal_int_add_unsigned, }
- /*
- Low-level addition, signed. Handbook of Applied Cryptography, algorithm 14.7.
- Assumptions:
- `dest`, `a` and `b` != `nil` and have been initalized.
- */
- internal_int_add_signed :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- x := a; y := b
- context.allocator = allocator
- /*
- Handle both negative or both positive.
- */
- if x.sign == y.sign {
- dest.sign = x.sign
- return #force_inline internal_int_add_unsigned(dest, x, y)
- }
- /*
- One positive, the other negative.
- Subtract the one with the greater magnitude from the other.
- The result gets the sign of the one with the greater magnitude.
- */
- if #force_inline internal_lt_abs(a, b) {
- x, y = y, x
- }
- dest.sign = x.sign
- return #force_inline internal_int_sub_unsigned(dest, x, y)
- }
- internal_add_signed :: proc { internal_int_add_signed, }
- /*
- Low-level addition Int+DIGIT, signed. Handbook of Applied Cryptography, algorithm 14.7.
- Assumptions:
- `dest` and `a` != `nil` and have been initalized.
- `dest` is large enough (a.used + 1) to fit result.
- */
- internal_int_add_digit :: proc(dest, a: ^Int, digit: DIGIT, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- internal_grow(dest, a.used + 1) or_return
- /*
- Fast paths for destination and input Int being the same.
- */
- if dest == a {
- /*
- Fast path for dest.digit[0] + digit fits in dest.digit[0] without overflow.
- */
- if dest.sign == .Zero_or_Positive && (dest.digit[0] + digit < _DIGIT_MAX) {
- dest.digit[0] += digit
- dest.used += 1
- return internal_clamp(dest)
- }
- /*
- Can be subtracted from dest.digit[0] without underflow.
- */
- if a.sign == .Negative && (dest.digit[0] > digit) {
- dest.digit[0] -= digit
- dest.used += 1
- return internal_clamp(dest)
- }
- }
- /*
- If `a` is negative and `|a|` >= `digit`, call `dest = |a| - digit`
- */
- if a.sign == .Negative && (a.used > 1 || a.digit[0] >= digit) {
- /*
- Temporarily fix `a`'s sign.
- */
- a.sign = .Zero_or_Positive
- /*
- dest = |a| - digit
- */
- if err = #force_inline internal_int_add_digit(dest, a, digit); err != nil {
- /*
- Restore a's sign.
- */
- a.sign = .Negative
- return err
- }
- /*
- Restore sign and set `dest` sign.
- */
- a.sign = .Negative
- dest.sign = .Negative
- return internal_clamp(dest)
- }
- /*
- Remember the currently used number of digits in `dest`.
- */
- old_used := dest.used
- /*
- If `a` is positive
- */
- if a.sign == .Zero_or_Positive {
- /*
- Add digits, use `carry`.
- */
- i: int
- carry := digit
- #no_bounds_check for i = 0; i < a.used; i += 1 {
- dest.digit[i] = a.digit[i] + carry
- carry = dest.digit[i] >> _DIGIT_BITS
- dest.digit[i] &= _MASK
- }
- /*
- Set final carry.
- */
- dest.digit[i] = carry
- /*
- Set `dest` size.
- */
- dest.used = a.used + 1
- } else {
- /*
- `a` was negative and |a| < digit.
- */
- dest.used = 1
- /*
- The result is a single DIGIT.
- */
- dest.digit[0] = digit - a.digit[0] if a.used == 1 else digit
- }
- /*
- Sign is always positive.
- */
- dest.sign = .Zero_or_Positive
- /*
- Zero remainder.
- */
- internal_zero_unused(dest, old_used)
- /*
- Adjust dest.used based on leading zeroes.
- */
- return internal_clamp(dest)
- }
- internal_add :: proc { internal_int_add_signed, internal_int_add_digit, }
- internal_int_incr :: proc(dest: ^Int, allocator := context.allocator) -> (err: Error) {
- return #force_inline internal_add(dest, dest, 1)
- }
- internal_incr :: proc { internal_int_incr, }
- /*
- Low-level subtraction, dest = number - decrease. Assumes |number| > |decrease|.
- Handbook of Applied Cryptography, algorithm 14.9.
- Assumptions:
- `dest`, `number` and `decrease` != `nil` and have been initalized.
- */
- internal_int_sub_unsigned :: proc(dest, number, decrease: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- dest := dest; x := number; y := decrease
- old_used := dest.used
- min_used := y.used
- max_used := x.used
- i: int
- grow(dest, max(max_used, _DEFAULT_DIGIT_COUNT)) or_return
- dest.used = max_used
- /*
- All parameters have been initialized.
- */
- borrow := DIGIT(0)
- #no_bounds_check for i = 0; i < min_used; i += 1 {
- dest.digit[i] = (x.digit[i] - y.digit[i] - borrow)
- /*
- borrow = carry bit of dest[i]
- Note this saves performing an AND operation since if a carry does occur,
- it will propagate all the way to the MSB.
- As a result a single shift is enough to get the carry.
- */
- borrow = dest.digit[i] >> ((size_of(DIGIT) * 8) - 1)
- /*
- Clear borrow from dest[i].
- */
- dest.digit[i] &= _MASK
- }
- /*
- Now copy higher words if any, e.g. if A has more digits than B
- */
- #no_bounds_check for ; i < max_used; i += 1 {
- dest.digit[i] = x.digit[i] - borrow
- /*
- borrow = carry bit of dest[i]
- Note this saves performing an AND operation since if a carry does occur,
- it will propagate all the way to the MSB.
- As a result a single shift is enough to get the carry.
- */
- borrow = dest.digit[i] >> ((size_of(DIGIT) * 8) - 1)
- /*
- Clear borrow from dest[i].
- */
- dest.digit[i] &= _MASK
- }
- /*
- Zero remainder.
- */
- internal_zero_unused(dest, old_used)
- /*
- Adjust dest.used based on leading zeroes.
- */
- return internal_clamp(dest)
- }
- internal_sub_unsigned :: proc { internal_int_sub_unsigned, }
- /*
- Low-level subtraction, signed. Handbook of Applied Cryptography, algorithm 14.9.
- dest = number - decrease. Assumes |number| > |decrease|.
- Assumptions:
- `dest`, `number` and `decrease` != `nil` and have been initalized.
- */
- internal_int_sub_signed :: proc(dest, number, decrease: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- number := number; decrease := decrease
- if number.sign != decrease.sign {
- /*
- Subtract a negative from a positive, OR subtract a positive from a negative.
- In either case, ADD their magnitudes and use the sign of the first number.
- */
- dest.sign = number.sign
- return #force_inline internal_int_add_unsigned(dest, number, decrease)
- }
- /*
- Subtract a positive from a positive, OR negative from a negative.
- First, take the difference between their magnitudes, then...
- */
- if #force_inline internal_lt_abs(number, decrease) {
- /*
- The second has a larger magnitude.
- The result has the *opposite* sign from the first number.
- */
- dest.sign = .Negative if number.sign == .Zero_or_Positive else .Zero_or_Positive
- number, decrease = decrease, number
- } else {
- /*
- The first has a larger or equal magnitude.
- Copy the sign from the first.
- */
- dest.sign = number.sign
- }
- return #force_inline internal_int_sub_unsigned(dest, number, decrease)
- }
- /*
- Low-level subtraction, signed. Handbook of Applied Cryptography, algorithm 14.9.
- dest = number - decrease. Assumes |number| > |decrease|.
- Assumptions:
- `dest`, `number` != `nil` and have been initalized.
- `dest` is large enough (number.used + 1) to fit result.
- */
- internal_int_sub_digit :: proc(dest, number: ^Int, digit: DIGIT, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- internal_grow(dest, number.used + 1) or_return
- dest := dest; digit := digit
- /*
- All parameters have been initialized.
- Fast paths for destination and input Int being the same.
- */
- if dest == number {
- /*
- Fast path for `dest` is negative and unsigned addition doesn't overflow the lowest digit.
- */
- if dest.sign == .Negative && (dest.digit[0] + digit < _DIGIT_MAX) {
- dest.digit[0] += digit
- return nil
- }
- /*
- Can be subtracted from dest.digit[0] without underflow.
- */
- if number.sign == .Zero_or_Positive && (dest.digit[0] > digit) {
- dest.digit[0] -= digit
- return nil
- }
- }
- /*
- If `a` is negative, just do an unsigned addition (with fudged signs).
- */
- if number.sign == .Negative {
- t := number
- t.sign = .Zero_or_Positive
- err = #force_inline internal_int_add_digit(dest, t, digit)
- dest.sign = .Negative
- internal_clamp(dest)
- return err
- }
- old_used := dest.used
- /*
- if `a`<= digit, simply fix the single digit.
- */
- if number.used == 1 && (number.digit[0] <= digit) || number.used == 0 {
- dest.digit[0] = digit - number.digit[0] if number.used == 1 else digit
- dest.sign = .Negative
- dest.used = 1
- } else {
- dest.sign = .Zero_or_Positive
- dest.used = number.used
- /*
- Subtract with carry.
- */
- carry := digit
- #no_bounds_check for i := 0; i < number.used; i += 1 {
- dest.digit[i] = number.digit[i] - carry
- carry = dest.digit[i] >> (_DIGIT_TYPE_BITS - 1)
- dest.digit[i] &= _MASK
- }
- }
- /*
- Zero remainder.
- */
- internal_zero_unused(dest, old_used)
- /*
- Adjust dest.used based on leading zeroes.
- */
- return internal_clamp(dest)
- }
- internal_sub :: proc { internal_int_sub_signed, internal_int_sub_digit, }
- internal_int_decr :: proc(dest: ^Int, allocator := context.allocator) -> (err: Error) {
- return #force_inline internal_sub(dest, dest, 1)
- }
- internal_decr :: proc { internal_int_decr, }
- /*
- dest = src / 2
- dest = src >> 1
- Assumes `dest` and `src` not to be `nil` and have been initialized.
- We make no allocations here.
- */
- internal_int_shr1 :: proc(dest, src: ^Int) -> (err: Error) {
- old_used := dest.used; dest.used = src.used
- /*
- Carry
- */
- fwd_carry := DIGIT(0)
- #no_bounds_check for x := dest.used - 1; x >= 0; x -= 1 {
- /*
- Get the carry for the next iteration.
- */
- src_digit := src.digit[x]
- carry := src_digit & 1
- /*
- Shift the current digit, add in carry and store.
- */
- dest.digit[x] = (src_digit >> 1) | (fwd_carry << (_DIGIT_BITS - 1))
- /*
- Forward carry to next iteration.
- */
- fwd_carry = carry
- }
- /*
- Zero remainder.
- */
- internal_zero_unused(dest, old_used)
- /*
- Adjust dest.used based on leading zeroes.
- */
- dest.sign = src.sign
- return internal_clamp(dest)
- }
- /*
- dest = src * 2
- dest = src << 1
- */
- internal_int_shl1 :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- internal_copy(dest, src) or_return
- /*
- Grow `dest` to accommodate the additional bits.
- */
- digits_needed := dest.used + 1
- internal_grow(dest, digits_needed) or_return
- dest.used = digits_needed
- mask := (DIGIT(1) << uint(1)) - DIGIT(1)
- shift := DIGIT(_DIGIT_BITS - 1)
- carry := DIGIT(0)
- #no_bounds_check for x:= 0; x < dest.used; x+= 1 {
- fwd_carry := (dest.digit[x] >> shift) & mask
- dest.digit[x] = (dest.digit[x] << uint(1) | carry) & _MASK
- carry = fwd_carry
- }
- /*
- Use final carry.
- */
- if carry != 0 {
- dest.digit[dest.used] = carry
- dest.used += 1
- }
- return internal_clamp(dest)
- }
- /*
- Multiply bigint `a` with int `d` and put the result in `dest`.
- Like `internal_int_mul_digit` but with an integer as the small input.
- */
- internal_int_mul_integer :: proc(dest, a: ^Int, b: $T, allocator := context.allocator) -> (err: Error)
- where intrinsics.type_is_integer(T) && T != DIGIT {
- context.allocator = allocator
- t := &Int{}
- defer internal_destroy(t)
- /*
- DIGIT might be smaller than a long, which excludes the use of `internal_int_mul_digit` here.
- */
- internal_set(t, b) or_return
- internal_mul(dest, a, t) or_return
- return
- }
- /*
- Multiply by a DIGIT.
- */
- internal_int_mul_digit :: proc(dest, src: ^Int, multiplier: DIGIT, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- assert_if_nil(dest, src)
- if multiplier == 0 {
- return internal_zero(dest)
- }
- if multiplier == 1 {
- return internal_copy(dest, src)
- }
- /*
- Power of two?
- */
- if multiplier == 2 {
- return #force_inline internal_int_shl1(dest, src)
- }
- if #force_inline platform_int_is_power_of_two(int(multiplier)) {
- ix := internal_log(multiplier, 2) or_return
- return internal_shl(dest, src, ix)
- }
- /*
- Ensure `dest` is big enough to hold `src` * `multiplier`.
- */
- grow(dest, max(src.used + 1, _DEFAULT_DIGIT_COUNT)) or_return
- /*
- Save the original used count.
- */
- old_used := dest.used
- /*
- Set the sign.
- */
- dest.sign = src.sign
- /*
- Set up carry.
- */
- carry := _WORD(0)
- /*
- Compute columns.
- */
- ix := 0
- #no_bounds_check for ; ix < src.used; ix += 1 {
- /*
- Compute product and carry sum for this term
- */
- product := carry + _WORD(src.digit[ix]) * _WORD(multiplier)
- /*
- Mask off higher bits to get a single DIGIT.
- */
- dest.digit[ix] = DIGIT(product & _WORD(_MASK))
- /*
- Send carry into next iteration
- */
- carry = product >> _DIGIT_BITS
- }
- /*
- Store final carry [if any] and increment used.
- */
- dest.digit[ix] = DIGIT(carry)
- dest.used = src.used + 1
- /*
- Zero remainder.
- */
- internal_zero_unused(dest, old_used)
- return internal_clamp(dest)
- }
- /*
- High level multiplication (handles sign).
- */
- internal_int_mul :: proc(dest, src, multiplier: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- /*
- Early out for `multiplier` is zero; Set `dest` to zero.
- */
- if multiplier.used == 0 || src.used == 0 { return internal_zero(dest) }
- neg := src.sign != multiplier.sign
- if src == multiplier {
- /*
- Do we need to square?
- */
- if src.used >= SQR_TOOM_CUTOFF {
- /*
- Use Toom-Cook?
- */
- err = #force_inline _private_int_sqr_toom(dest, src)
- } else if src.used >= SQR_KARATSUBA_CUTOFF {
- /*
- Karatsuba?
- */
- err = #force_inline _private_int_sqr_karatsuba(dest, src)
- } else if ((src.used * 2) + 1) < _WARRAY && src.used < (_MAX_COMBA / 2) {
- /*
- Fast comba?
- */
- err = #force_inline _private_int_sqr_comba(dest, src)
- } else {
- err = #force_inline _private_int_sqr(dest, src)
- }
- } else {
- /*
- Can we use the balance method? Check sizes.
- * The smaller one needs to be larger than the Karatsuba cut-off.
- * The bigger one needs to be at least about one `_MUL_KARATSUBA_CUTOFF` bigger
- * to make some sense, but it depends on architecture, OS, position of the stars... so YMMV.
- * Using it to cut the input into slices small enough for _mul_comba
- * was actually slower on the author's machine, but YMMV.
- */
- min_used := min(src.used, multiplier.used)
- max_used := max(src.used, multiplier.used)
- digits := src.used + multiplier.used + 1
- if min_used >= MUL_KARATSUBA_CUTOFF && (max_used / 2) >= MUL_KARATSUBA_CUTOFF && max_used >= (2 * min_used) {
- /*
- Not much effect was observed below a ratio of 1:2, but again: YMMV.
- */
- err = _private_int_mul_balance(dest, src, multiplier)
- } else if min_used >= MUL_TOOM_CUTOFF {
- /*
- Toom path commented out until it no longer fails Factorial 10k or 100k,
- as reveaved in the long test.
- */
- err = #force_inline _private_int_mul_toom(dest, src, multiplier)
- } else if min_used >= MUL_KARATSUBA_CUTOFF {
- err = #force_inline _private_int_mul_karatsuba(dest, src, multiplier)
- } else if digits < _WARRAY && min_used <= _MAX_COMBA {
- /*
- Can we use the fast multiplier?
- * The fast multiplier can be used if the output will
- * have less than MP_WARRAY digits and the number of
- * digits won't affect carry propagation
- */
- err = #force_inline _private_int_mul_comba(dest, src, multiplier, digits)
- } else {
- err = #force_inline _private_int_mul(dest, src, multiplier, digits)
- }
- }
- dest.sign = .Negative if dest.used > 0 && neg else .Zero_or_Positive
- return err
- }
- internal_mul :: proc { internal_int_mul, internal_int_mul_digit, internal_int_mul_integer }
- internal_sqr :: proc (dest, src: ^Int, allocator := context.allocator) -> (res: Error) {
- /*
- We call `internal_mul` and not e.g. `_private_int_sqr` because the former
- will dispatch to the optimal implementation depending on the source.
- */
- return #force_inline internal_mul(dest, src, src, allocator)
- }
- /*
- divmod.
- Both the quotient and remainder are optional and may be passed a nil.
- `numerator` and `denominator` are expected not to be `nil` and have been initialized.
- */
- internal_int_divmod :: proc(quotient, remainder, numerator, denominator: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- if denominator.used == 0 { return .Division_by_Zero }
- /*
- If numerator < denominator then quotient = 0, remainder = numerator.
- */
- if #force_inline internal_lt_abs(numerator, denominator) {
- if remainder != nil {
- internal_copy(remainder, numerator) or_return
- }
- if quotient != nil {
- internal_zero(quotient)
- }
- return nil
- }
- if (denominator.used > 2 * MUL_KARATSUBA_CUTOFF) && (denominator.used <= (numerator.used / 3) * 2) {
- assert(denominator.used >= 160 && numerator.used >= 240, "MUL_KARATSUBA_CUTOFF global not properly set.")
- err = _private_int_div_recursive(quotient, remainder, numerator, denominator)
- } else {
- when true {
- err = #force_inline _private_int_div_school(quotient, remainder, numerator, denominator)
- } else {
- /*
- NOTE(Jeroen): We no longer need or use `_private_int_div_small`.
- We'll keep it around for a bit until we're reasonably certain div_school is bug free.
- */
- err = _private_int_div_small(quotient, remainder, numerator, denominator)
- }
- }
- return
- }
- /*
- Single digit division (based on routine from MPI).
- The quotient is optional and may be passed a nil.
- */
- internal_int_divmod_digit :: proc(quotient, numerator: ^Int, denominator: DIGIT, allocator := context.allocator) -> (remainder: DIGIT, err: Error) {
- context.allocator = allocator
- /*
- Cannot divide by zero.
- */
- if denominator == 0 { return 0, .Division_by_Zero }
- /*
- Quick outs.
- */
- if denominator == 1 || numerator.used == 0 {
- if quotient != nil {
- return 0, internal_copy(quotient, numerator)
- }
- return 0, err
- }
- /*
- Power of two?
- */
- if denominator == 2 {
- if numerator.used > 0 && numerator.digit[0] & 1 != 0 {
- // Remainder is 1 if numerator is odd.
- remainder = 1
- }
- if quotient == nil {
- return remainder, nil
- }
- return remainder, internal_shr(quotient, numerator, 1)
- }
- ix: int
- if platform_int_is_power_of_two(int(denominator)) {
- ix = 1
- for ix < _DIGIT_BITS && denominator != (1 << uint(ix)) {
- ix += 1
- }
- remainder = numerator.digit[0] & ((1 << uint(ix)) - 1)
- if quotient == nil {
- return remainder, nil
- }
- return remainder, internal_shr(quotient, numerator, int(ix))
- }
- /*
- Three?
- */
- if denominator == 3 {
- return _private_int_div_3(quotient, numerator)
- }
- /*
- No easy answer [c'est la vie]. Just division.
- */
- q := &Int{}
- internal_grow(q, numerator.used) or_return
- q.used = numerator.used
- q.sign = numerator.sign
- w := _WORD(0)
- for ix = numerator.used - 1; ix >= 0; ix -= 1 {
- t := DIGIT(0)
- w = (w << _WORD(_DIGIT_BITS) | _WORD(numerator.digit[ix]))
- if w >= _WORD(denominator) {
- t = DIGIT(w / _WORD(denominator))
- w -= _WORD(t) * _WORD(denominator)
- }
- q.digit[ix] = t
- }
- remainder = DIGIT(w)
- if quotient != nil {
- internal_clamp(q)
- internal_swap(q, quotient)
- }
- internal_destroy(q)
- return remainder, nil
- }
- internal_divmod :: proc { internal_int_divmod, internal_int_divmod_digit, }
- /*
- Asssumes quotient, numerator and denominator to have been initialized and not to be nil.
- */
- internal_int_div :: proc(quotient, numerator, denominator: ^Int, allocator := context.allocator) -> (err: Error) {
- return #force_inline internal_int_divmod(quotient, nil, numerator, denominator, allocator)
- }
- internal_div :: proc { internal_int_div, }
- /*
- remainder = numerator % denominator.
- 0 <= remainder < denominator if denominator > 0
- denominator < remainder <= 0 if denominator < 0
- Asssumes quotient, numerator and denominator to have been initialized and not to be nil.
- */
- internal_int_mod :: proc(remainder, numerator, denominator: ^Int, allocator := context.allocator) -> (err: Error) {
- #force_inline internal_int_divmod(nil, remainder, numerator, denominator, allocator) or_return
- if remainder.used == 0 || denominator.sign == remainder.sign { return nil }
- return #force_inline internal_add(remainder, remainder, denominator, allocator)
- }
- internal_int_mod_digit :: proc(numerator: ^Int, denominator: DIGIT, allocator := context.allocator) -> (remainder: DIGIT, err: Error) {
- return internal_int_divmod_digit(nil, numerator, denominator, allocator)
- }
- internal_mod :: proc{ internal_int_mod, internal_int_mod_digit, }
- /*
- remainder = (number + addend) % modulus.
- */
- internal_int_addmod :: proc(remainder, number, addend, modulus: ^Int, allocator := context.allocator) -> (err: Error) {
- #force_inline internal_add(remainder, number, addend, allocator) or_return
- return #force_inline internal_mod(remainder, remainder, modulus, allocator)
- }
- internal_addmod :: proc { internal_int_addmod, }
- /*
- remainder = (number - decrease) % modulus.
- */
- internal_int_submod :: proc(remainder, number, decrease, modulus: ^Int, allocator := context.allocator) -> (err: Error) {
- #force_inline internal_sub(remainder, number, decrease, allocator) or_return
- return #force_inline internal_mod(remainder, remainder, modulus, allocator)
- }
- internal_submod :: proc { internal_int_submod, }
- /*
- remainder = (number * multiplicand) % modulus.
- */
- internal_int_mulmod :: proc(remainder, number, multiplicand, modulus: ^Int, allocator := context.allocator) -> (err: Error) {
- #force_inline internal_mul(remainder, number, multiplicand, allocator) or_return
- return #force_inline internal_mod(remainder, remainder, modulus, allocator)
- }
- internal_mulmod :: proc { internal_int_mulmod, }
- /*
- remainder = (number * number) % modulus.
- */
- internal_int_sqrmod :: proc(remainder, number, modulus: ^Int, allocator := context.allocator) -> (err: Error) {
- #force_inline internal_sqr(remainder, number, allocator) or_return
- return #force_inline internal_mod(remainder, remainder, modulus, allocator)
- }
- internal_sqrmod :: proc { internal_int_sqrmod, }
- /*
- TODO: Use Sterling's Approximation to estimate log2(N!) to size the result.
- This way we'll have to reallocate less, possibly not at all.
- */
- internal_int_factorial :: proc(res: ^Int, n: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- if n >= FACTORIAL_BINARY_SPLIT_CUTOFF {
- return _private_int_factorial_binary_split(res, n)
- }
- i := len(_factorial_table)
- if n < i {
- return #force_inline internal_set(res, _factorial_table[n])
- }
- #force_inline internal_set(res, _factorial_table[i - 1]) or_return
- for {
- if err = #force_inline internal_mul(res, res, DIGIT(i)); err != nil || i == n {
- return err
- }
- i += 1
- }
- return nil
- }
- /*
- Returns GCD, LCM or both.
- Assumes `a` and `b` to have been initialized.
- `res_gcd` and `res_lcm` can be nil or ^Int depending on which results are desired.
- */
- internal_int_gcd_lcm :: proc(res_gcd, res_lcm, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- if res_gcd == nil && res_lcm == nil { return nil }
- return #force_inline _private_int_gcd_lcm(res_gcd, res_lcm, a, b, allocator)
- }
- internal_int_gcd :: proc(res_gcd, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- return #force_inline _private_int_gcd_lcm(res_gcd, nil, a, b, allocator)
- }
- internal_int_lcm :: proc(res_lcm, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- return #force_inline _private_int_gcd_lcm(nil, res_lcm, a, b, allocator)
- }
- /*
- remainder = numerator % (1 << bits)
- Assumes `remainder` and `numerator` both not to be `nil` and `bits` to be >= 0.
- */
- internal_int_mod_bits :: proc(remainder, numerator: ^Int, bits: int, allocator := context.allocator) -> (err: Error) {
- /*
- Everything is divisible by 1 << 0 == 1, so this returns 0.
- */
- if bits == 0 { return internal_zero(remainder) }
- /*
- If the modulus is larger than the value, return the value.
- */
- internal_copy(remainder, numerator) or_return
- if bits >= (numerator.used * _DIGIT_BITS) {
- return
- }
- /*
- Zero digits above the last digit of the modulus.
- */
- zero_count := (bits / _DIGIT_BITS)
- zero_count += 0 if (bits % _DIGIT_BITS == 0) else 1
- /*
- Zero remainder. Special case, can't use `internal_zero_unused`.
- */
- if zero_count > 0 {
- mem.zero_slice(remainder.digit[zero_count:])
- }
- /*
- Clear the digit that is not completely outside/inside the modulus.
- */
- remainder.digit[bits / _DIGIT_BITS] &= DIGIT(1 << DIGIT(bits % _DIGIT_BITS)) - DIGIT(1)
- return internal_clamp(remainder)
- }
- /*
- ============================= Low-level helpers =============================
- `internal_*` helpers don't return an `Error` like their public counterparts do,
- because they expect not to be passed `nil` or uninitialized inputs.
- This makes them more suitable for `internal_*` functions and some of the
- public ones that have already satisfied these constraints.
- */
- /*
- This procedure returns the allocated capacity of an Int.
- Assumes `a` not to be `nil`.
- */
- internal_int_allocated_cap :: #force_inline proc(a: ^Int) -> (cap: int) {
- raw := transmute(mem.Raw_Dynamic_Array)a.digit
- return raw.cap
- }
- /*
- This procedure will return `true` if the `Int` is initialized, `false` if not.
- Assumes `a` not to be `nil`.
- */
- internal_int_is_initialized :: #force_inline proc(a: ^Int) -> (initialized: bool) {
- return internal_int_allocated_cap(a) >= _MIN_DIGIT_COUNT
- }
- internal_is_initialized :: proc { internal_int_is_initialized, }
- /*
- This procedure will return `true` if the `Int` is zero, `false` if not.
- Assumes `a` not to be `nil`.
- */
- internal_int_is_zero :: #force_inline proc(a: ^Int) -> (zero: bool) {
- return a.used == 0
- }
- internal_is_zero :: proc { internal_int_is_zero }
- /*
- This procedure will return `true` if the `Int` is positive, `false` if not.
- Assumes `a` not to be `nil`.
- */
- internal_int_is_positive :: #force_inline proc(a: ^Int) -> (positive: bool) {
- return a.sign == .Zero_or_Positive
- }
- internal_is_positive :: proc { internal_int_is_positive, }
- /*
- This procedure will return `true` if the `Int` is negative, `false` if not.
- Assumes `a` not to be `nil`.
- */
- internal_int_is_negative :: #force_inline proc(a: ^Int) -> (negative: bool) {
- return a.sign == .Negative
- }
- internal_is_negative :: proc { internal_int_is_negative, }
- /*
- This procedure will return `true` if the `Int` is even, `false` if not.
- Assumes `a` not to be `nil`.
- */
- internal_int_is_even :: #force_inline proc(a: ^Int) -> (even: bool) {
- if internal_is_zero(a) { return true }
- /*
- `a.used` > 0 here, because the above handled `is_zero`.
- We don't need to explicitly test it.
- */
- return a.digit[0] & 1 == 0
- }
- internal_is_even :: proc { internal_int_is_even, }
- /*
- This procedure will return `true` if the `Int` is even, `false` if not.
- Assumes `a` not to be `nil`.
- */
- internal_int_is_odd :: #force_inline proc(a: ^Int) -> (odd: bool) {
- return !internal_int_is_even(a)
- }
- internal_is_odd :: proc { internal_int_is_odd, }
- /*
- This procedure will return `true` if the `Int` is a power of two, `false` if not.
- Assumes `a` not to be `nil`.
- */
- internal_int_is_power_of_two :: #force_inline proc(a: ^Int) -> (power_of_two: bool) {
- /*
- Early out for Int == 0.
- */
- if #force_inline internal_is_zero(a) { return true }
- /*
- For an `Int` to be a power of two, its bottom limb has to be a power of two.
- */
- if ! #force_inline platform_int_is_power_of_two(int(a.digit[a.used - 1])) { return false }
- /*
- We've established that the bottom limb is a power of two.
- If it's the only limb, that makes the entire Int a power of two.
- */
- if a.used == 1 { return true }
- /*
- For an `Int` to be a power of two, all limbs except the top one have to be zero.
- */
- for i := 1; i < a.used && a.digit[i - 1] != 0; i += 1 { return false }
- return true
- }
- internal_is_power_of_two :: proc { internal_int_is_power_of_two, }
- /*
- Compare two `Int`s, signed.
- Returns -1 if `a` < `b`, 0 if `a` == `b` and 1 if `b` > `a`.
- Expects `a` and `b` both to be valid `Int`s, i.e. initialized and not `nil`.
- */
- internal_int_compare :: #force_inline proc(a, b: ^Int) -> (comparison: int) {
- assert_if_nil(a, b)
- a_is_negative := #force_inline internal_is_negative(a)
- /*
- Compare based on sign.
- */
- if a.sign != b.sign { return -1 if a_is_negative else +1 }
- /*
- If `a` is negative, compare in the opposite direction */
- if a_is_negative { return #force_inline internal_compare_magnitude(b, a) }
- return #force_inline internal_compare_magnitude(a, b)
- }
- internal_compare :: proc { internal_int_compare, internal_int_compare_digit, }
- internal_cmp :: internal_compare
- /*
- Compare an `Int` to an unsigned number upto `DIGIT & _MASK`.
- Returns -1 if `a` < `b`, 0 if `a` == `b` and 1 if `b` > `a`.
- Expects: `a` and `b` both to be valid `Int`s, i.e. initialized and not `nil`.
- */
- internal_int_compare_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (comparison: int) {
- assert_if_nil(a)
- a_is_negative := #force_inline internal_is_negative(a)
- switch {
- /*
- Compare based on sign first.
- */
- case a_is_negative: return -1
- /*
- Then compare on magnitude.
- */
- case a.used > 1: return +1
- /*
- We have only one digit. Compare it against `b`.
- */
- case a.digit[0] < b: return -1
- case a.digit[0] == b: return 0
- case a.digit[0] > b: return +1
- /*
- Unreachable.
- Just here because Odin complains about a missing return value at the bottom of the proc otherwise.
- */
- case: return
- }
- }
- internal_compare_digit :: proc { internal_int_compare_digit, }
- internal_cmp_digit :: internal_compare_digit
- /*
- Compare the magnitude of two `Int`s, unsigned.
- */
- internal_int_compare_magnitude :: #force_inline proc(a, b: ^Int) -> (comparison: int) {
- assert_if_nil(a, b)
- /*
- Compare based on used digits.
- */
- if a.used != b.used {
- if a.used > b.used {
- return +1
- }
- return -1
- }
- /*
- Same number of used digits, compare based on their value.
- */
- #no_bounds_check for n := a.used - 1; n >= 0; n -= 1 {
- if a.digit[n] != b.digit[n] {
- if a.digit[n] > b.digit[n] {
- return +1
- }
- return -1
- }
- }
- return 0
- }
- internal_compare_magnitude :: proc { internal_int_compare_magnitude, }
- internal_cmp_mag :: internal_compare_magnitude
- /*
- bool := a < b
- */
- internal_int_less_than :: #force_inline proc(a, b: ^Int) -> (less_than: bool) {
- return internal_cmp(a, b) == -1
- }
- /*
- bool := a < b
- */
- internal_int_less_than_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (less_than: bool) {
- return internal_cmp_digit(a, b) == -1
- }
- /*
- bool := |a| < |b|
- Compares the magnitudes only, ignores the sign.
- */
- internal_int_less_than_abs :: #force_inline proc(a, b: ^Int) -> (less_than: bool) {
- return internal_cmp_mag(a, b) == -1
- }
- internal_less_than :: proc {
- internal_int_less_than,
- internal_int_less_than_digit,
- }
- internal_lt :: internal_less_than
- internal_less_than_abs :: proc {
- internal_int_less_than_abs,
- }
- internal_lt_abs :: internal_less_than_abs
- /*
- bool := a <= b
- */
- internal_int_less_than_or_equal :: #force_inline proc(a, b: ^Int) -> (less_than_or_equal: bool) {
- return internal_cmp(a, b) <= 0
- }
- /*
- bool := a <= b
- */
- internal_int_less_than_or_equal_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (less_than_or_equal: bool) {
- return internal_cmp_digit(a, b) <= 0
- }
- /*
- bool := |a| <= |b|
- Compares the magnitudes only, ignores the sign.
- */
- internal_int_less_than_or_equal_abs :: #force_inline proc(a, b: ^Int) -> (less_than_or_equal: bool) {
- return internal_cmp_mag(a, b) <= 0
- }
- internal_less_than_or_equal :: proc {
- internal_int_less_than_or_equal,
- internal_int_less_than_or_equal_digit,
- }
- internal_lte :: internal_less_than_or_equal
- internal_less_than_or_equal_abs :: proc {
- internal_int_less_than_or_equal_abs,
- }
- internal_lte_abs :: internal_less_than_or_equal_abs
- /*
- bool := a == b
- */
- internal_int_equals :: #force_inline proc(a, b: ^Int) -> (equals: bool) {
- return internal_cmp(a, b) == 0
- }
- /*
- bool := a == b
- */
- internal_int_equals_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (equals: bool) {
- return internal_cmp_digit(a, b) == 0
- }
- /*
- bool := |a| == |b|
- Compares the magnitudes only, ignores the sign.
- */
- internal_int_equals_abs :: #force_inline proc(a, b: ^Int) -> (equals: bool) {
- return internal_cmp_mag(a, b) == 0
- }
- internal_equals :: proc {
- internal_int_equals,
- internal_int_equals_digit,
- }
- internal_eq :: internal_equals
- internal_equals_abs :: proc {
- internal_int_equals_abs,
- }
- internal_eq_abs :: internal_equals_abs
- /*
- bool := a >= b
- */
- internal_int_greater_than_or_equal :: #force_inline proc(a, b: ^Int) -> (greater_than_or_equal: bool) {
- return internal_cmp(a, b) >= 0
- }
- /*
- bool := a >= b
- */
- internal_int_greater_than_or_equal_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (greater_than_or_equal: bool) {
- return internal_cmp_digit(a, b) >= 0
- }
- /*
- bool := |a| >= |b|
- Compares the magnitudes only, ignores the sign.
- */
- internal_int_greater_than_or_equal_abs :: #force_inline proc(a, b: ^Int) -> (greater_than_or_equal: bool) {
- return internal_cmp_mag(a, b) >= 0
- }
- internal_greater_than_or_equal :: proc {
- internal_int_greater_than_or_equal,
- internal_int_greater_than_or_equal_digit,
- }
- internal_gte :: internal_greater_than_or_equal
- internal_greater_than_or_equal_abs :: proc {
- internal_int_greater_than_or_equal_abs,
- }
- internal_gte_abs :: internal_greater_than_or_equal_abs
- /*
- bool := a > b
- */
- internal_int_greater_than :: #force_inline proc(a, b: ^Int) -> (greater_than: bool) {
- return internal_cmp(a, b) == 1
- }
- /*
- bool := a > b
- */
- internal_int_greater_than_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (greater_than: bool) {
- return internal_cmp_digit(a, b) == 1
- }
- /*
- bool := |a| > |b|
- Compares the magnitudes only, ignores the sign.
- */
- internal_int_greater_than_abs :: #force_inline proc(a, b: ^Int) -> (greater_than: bool) {
- return internal_cmp_mag(a, b) == 1
- }
- internal_greater_than :: proc {
- internal_int_greater_than,
- internal_int_greater_than_digit,
- }
- internal_gt :: internal_greater_than
- internal_greater_than_abs :: proc {
- internal_int_greater_than_abs,
- }
- internal_gt_abs :: internal_greater_than_abs
- /*
- Check if remainders are possible squares - fast exclude non-squares.
- Returns `true` if `a` is a square, `false` if not.
- Assumes `a` not to be `nil` and to have been initialized.
- */
- internal_int_is_square :: proc(a: ^Int, allocator := context.allocator) -> (square: bool, err: Error) {
- context.allocator = allocator
- /*
- Default to Non-square :)
- */
- square = false
- if internal_is_negative(a) { return }
- if internal_is_zero(a) { return }
- /*
- First check mod 128 (suppose that _DIGIT_BITS is at least 7).
- */
- if _private_int_rem_128[127 & a.digit[0]] == 1 { return }
- /*
- Next check mod 105 (3*5*7).
- */
- c: DIGIT
- c, err = internal_mod(a, 105)
- if _private_int_rem_105[c] == 1 { return }
- t := &Int{}
- defer destroy(t)
- set(t, 11 * 13 * 17 * 19 * 23 * 29 * 31) or_return
- internal_mod(t, a, t) or_return
- r: u64
- r, err = internal_int_get(t, u64)
- /*
- Check for other prime modules, note it's not an ERROR but we must
- free "t" so the easiest way is to goto LBL_ERR. We know that err
- is already equal to MP_OKAY from the mp_mod call
- */
- if (1 << (r % 11) & 0x5C4) != 0 { return }
- if (1 << (r % 13) & 0x9E4) != 0 { return }
- if (1 << (r % 17) & 0x5CE8) != 0 { return }
- if (1 << (r % 19) & 0x4F50C) != 0 { return }
- if (1 << (r % 23) & 0x7ACCA0) != 0 { return }
- if (1 << (r % 29) & 0xC2EDD0C) != 0 { return }
- if (1 << (r % 31) & 0x6DE2B848) != 0 { return }
- /*
- Final check - is sqr(sqrt(arg)) == arg?
- */
- sqrt(t, a) or_return
- sqr(t, t) or_return
- square = internal_eq_abs(t, a)
- return
- }
- /*
- ========================= Logs, powers and roots ============================
- */
- /*
- Returns log_base(a).
- Assumes `a` to not be `nil` and have been iniialized.
- */
- internal_int_log :: proc(a: ^Int, base: DIGIT) -> (res: int, err: Error) {
- if base < 2 || DIGIT(base) > _DIGIT_MAX { return -1, .Invalid_Argument }
- if internal_is_negative(a) { return -1, .Math_Domain_Error }
- if internal_is_zero(a) { return -1, .Math_Domain_Error }
- /*
- Fast path for bases that are a power of two.
- */
- if platform_int_is_power_of_two(int(base)) { return _private_log_power_of_two(a, base) }
- /*
- Fast path for `Int`s that fit within a single `DIGIT`.
- */
- if a.used == 1 { return internal_log(a.digit[0], DIGIT(base)) }
- return _private_int_log(a, base)
- }
- /*
- Returns log_base(a), where `a` is a DIGIT.
- */
- internal_digit_log :: proc(a: DIGIT, base: DIGIT) -> (log: int, err: Error) {
- /*
- If the number is smaller than the base, it fits within a fraction.
- Therefore, we return 0.
- */
- if a < base { return 0, nil }
- /*
- If a number equals the base, the log is 1.
- */
- if a == base { return 1, nil }
- N := _WORD(a)
- bracket_low := _WORD(1)
- bracket_high := _WORD(base)
- high := 1
- low := 0
- for bracket_high < N {
- low = high
- bracket_low = bracket_high
- high <<= 1
- bracket_high *= bracket_high
- }
- for high - low > 1 {
- mid := (low + high) >> 1
- bracket_mid := bracket_low * #force_inline internal_small_pow(_WORD(base), _WORD(mid - low))
- if N < bracket_mid {
- high = mid
- bracket_high = bracket_mid
- }
- if N > bracket_mid {
- low = mid
- bracket_low = bracket_mid
- }
- if N == bracket_mid {
- return mid, nil
- }
- }
- if bracket_high == N {
- return high, nil
- } else {
- return low, nil
- }
- }
- internal_log :: proc { internal_int_log, internal_digit_log, }
- /*
- Calculate dest = base^power using a square-multiply algorithm.
- Assumes `dest` and `base` not to be `nil` and to have been initialized.
- */
- internal_int_pow :: proc(dest, base: ^Int, power: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- power := power
- /*
- Early outs.
- */
- if #force_inline internal_is_zero(base) {
- /*
- A zero base is a special case.
- */
- if power < 0 {
- internal_zero(dest) or_return
- return .Math_Domain_Error
- }
- if power == 0 { return internal_one(dest) }
- if power > 0 { return internal_zero(dest) }
- }
- if power < 0 {
- /*
- Fraction, so we'll return zero.
- */
- return internal_zero(dest)
- }
- switch(power) {
- case 0:
- /*
- Any base to the power zero is one.
- */
- return #force_inline internal_one(dest)
- case 1:
- /*
- Any base to the power one is itself.
- */
- return copy(dest, base)
- case 2:
- return #force_inline internal_sqr(dest, base)
- }
- g := &Int{}
- internal_copy(g, base) or_return
- /*
- Set initial result.
- */
- internal_one(dest) or_return
- defer internal_destroy(g)
- for power > 0 {
- /*
- If the bit is set, multiply.
- */
- if power & 1 != 0 {
- internal_mul(dest, g, dest) or_return
- }
- /*
- Square.
- */
- if power > 1 {
- internal_sqr(g, g) or_return
- }
- /* shift to next bit */
- power >>= 1
- }
- return
- }
- /*
- Calculate `dest = base^power`.
- Assumes `dest` not to be `nil` and to have been initialized.
- */
- internal_int_pow_int :: proc(dest: ^Int, base, power: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- base_t := &Int{}
- defer internal_destroy(base_t)
- internal_set(base_t, base) or_return
- return #force_inline internal_int_pow(dest, base_t, power)
- }
- internal_pow :: proc { internal_int_pow, internal_int_pow_int, }
- internal_exp :: pow
- /*
- */
- internal_small_pow :: proc(base: _WORD, exponent: _WORD) -> (result: _WORD) {
- exponent := exponent; base := base
- result = _WORD(1)
- for exponent != 0 {
- if exponent & 1 == 1 {
- result *= base
- }
- exponent >>= 1
- base *= base
- }
- return result
- }
- /*
- This function is less generic than `root_n`, simpler and faster.
- Assumes `dest` and `src` not to be `nil` and to have been initialized.
- */
- internal_int_sqrt :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- /*
- Must be positive.
- */
- if #force_inline internal_is_negative(src) { return .Invalid_Argument }
- /*
- Easy out. If src is zero, so is dest.
- */
- if #force_inline internal_is_zero(src) { return internal_zero(dest) }
- /*
- Set up temporaries.
- */
- x, y, t1, t2 := &Int{}, &Int{}, &Int{}, &Int{}
- defer internal_destroy(x, y, t1, t2)
- count := #force_inline internal_count_bits(src)
- a, b := count >> 1, count & 1
- internal_int_power_of_two(x, a+b, allocator) or_return
- for {
- /*
- y = (x + n // x) // 2
- */
- internal_div(t1, src, x) or_return
- internal_add(t2, t1, x) or_return
- internal_shr(y, t2, 1) or_return
- if internal_gte(y, x) {
- internal_swap(dest, x)
- return nil
- }
- internal_swap(x, y)
- }
- internal_swap(dest, x)
- return err
- }
- internal_sqrt :: proc { internal_int_sqrt, }
- /*
- Find the nth root of an Integer.
- Result found such that `(dest)**n <= src` and `(dest+1)**n > src`
- This algorithm uses Newton's approximation `x[i+1] = x[i] - f(x[i])/f'(x[i])`,
- which will find the root in `log(n)` time where each step involves a fair bit.
- Assumes `dest` and `src` not to be `nil` and have been initialized.
- */
- internal_int_root_n :: proc(dest, src: ^Int, n: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- /*
- Fast path for n == 2
- */
- if n == 2 { return #force_inline internal_sqrt(dest, src) }
- if n < 0 || n > int(_DIGIT_MAX) { return .Invalid_Argument }
- if n & 1 == 0 && #force_inline internal_is_negative(src) { return .Invalid_Argument }
- /*
- Set up temporaries.
- */
- t1, t2, t3, a := &Int{}, &Int{}, &Int{}, &Int{}
- defer internal_destroy(t1, t2, t3)
- /*
- If `src` is negative fudge the sign but keep track.
- */
- a.sign = .Zero_or_Positive
- a.used = src.used
- a.digit = src.digit
- /*
- If "n" is larger than INT_MAX it is also larger than
- log_2(src) because the bit-length of the "src" is measured
- with an int and hence the root is always < 2 (two).
- */
- if n > max(int) / 2 {
- err = set(dest, 1)
- dest.sign = a.sign
- return err
- }
- /*
- Compute seed: 2^(log_2(src)/n + 2)
- */
- ilog2 := internal_count_bits(src)
- /*
- "src" is smaller than max(int), we can cast safely.
- */
- if ilog2 < n {
- err = internal_one(dest)
- dest.sign = a.sign
- return err
- }
- ilog2 /= n
- if ilog2 == 0 {
- err = internal_one(dest)
- dest.sign = a.sign
- return err
- }
- /*
- Start value must be larger than root.
- */
- ilog2 += 2
- internal_int_power_of_two(t2, ilog2) or_return
- c: int
- iterations := 0
- for {
- /* t1 = t2 */
- internal_copy(t1, t2) or_return
- /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */
- /* t3 = t1**(b-1) */
- internal_pow(t3, t1, n-1) or_return
- /* numerator */
- /* t2 = t1**b */
- internal_mul(t2, t1, t3) or_return
- /* t2 = t1**b - a */
- internal_sub(t2, t2, a) or_return
- /* denominator */
- /* t3 = t1**(b-1) * b */
- internal_mul(t3, t3, DIGIT(n)) or_return
- /* t3 = (t1**b - a)/(b * t1**(b-1)) */
- internal_div(t3, t2, t3) or_return
- internal_sub(t2, t1, t3) or_return
- /*
- Number of rounds is at most log_2(root). If it is more it
- got stuck, so break out of the loop and do the rest manually.
- */
- if ilog2 -= 1; ilog2 == 0 { break }
- if internal_eq(t1, t2) { break }
- iterations += 1
- if iterations == MAX_ITERATIONS_ROOT_N {
- return .Max_Iterations_Reached
- }
- }
- /* Result can be off by a few so check. */
- /* Loop beneath can overshoot by one if found root is smaller than actual root. */
- iterations = 0
- for {
- internal_pow(t2, t1, n) or_return
- c = internal_cmp(t2, a)
- if c == 0 {
- swap(dest, t1)
- return nil
- } else if c == -1 {
- internal_add(t1, t1, DIGIT(1)) or_return
- } else {
- break
- }
- iterations += 1
- if iterations == MAX_ITERATIONS_ROOT_N {
- return .Max_Iterations_Reached
- }
- }
- iterations = 0
- /*
- Correct overshoot from above or from recurrence.
- */
- for {
- internal_pow(t2, t1, n) or_return
-
- if internal_lt(t2, a) { break }
-
- internal_sub(t1, t1, DIGIT(1)) or_return
- iterations += 1
- if iterations == MAX_ITERATIONS_ROOT_N {
- return .Max_Iterations_Reached
- }
- }
- /*
- Set the result.
- */
- internal_swap(dest, t1)
- /*
- Set the sign of the result.
- */
- dest.sign = src.sign
- return err
- }
- internal_root_n :: proc { internal_int_root_n, }
- /*
- Other internal helpers
- */
- /*
- Deallocates the backing memory of one or more `Int`s.
- Asssumes none of the `integers` to be a `nil`.
- */
- internal_int_destroy :: proc(integers: ..^Int) {
- integers := integers
- for a in &integers {
- if internal_int_allocated_cap(a) > 0 {
- mem.zero_slice(a.digit[:])
- free(&a.digit[0])
- }
- a = &Int{}
- }
- }
- internal_destroy :: proc{ internal_int_destroy, }
- /*
- Helpers to set an `Int` to a specific value.
- */
- internal_int_set_from_integer :: proc(dest: ^Int, src: $T, minimize := false, allocator := context.allocator) -> (err: Error)
- where intrinsics.type_is_integer(T) {
- context.allocator = allocator
- src := src
- internal_error_if_immutable(dest) or_return
- /*
- Most internal procs asssume an Int to have already been initialize,
- but as this is one of the procs that initializes, we have to check the following.
- */
- internal_clear_if_uninitialized_single(dest) or_return
- dest.flags = {} // We're not -Inf, Inf, NaN or Immutable.
- dest.used = 0
- dest.sign = .Zero_or_Positive if src >= 0 else .Negative
- src = internal_abs(src)
- #no_bounds_check for src != 0 {
- dest.digit[dest.used] = DIGIT(src) & _MASK
- dest.used += 1
- src >>= _DIGIT_BITS
- }
- internal_zero_unused(dest)
- return nil
- }
- internal_set :: proc { internal_int_set_from_integer, internal_int_copy, int_atoi }
- internal_copy_digits :: #force_inline proc(dest, src: ^Int, digits: int, offset := int(0)) -> (err: Error) {
- #force_inline internal_error_if_immutable(dest) or_return
- /*
- If dest == src, do nothing
- */
- return #force_inline _private_copy_digits(dest, src, digits, offset)
- }
- /*
- Copy one `Int` to another.
- */
- internal_int_copy :: proc(dest, src: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- /*
- If dest == src, do nothing
- */
- if (dest == src) { return nil }
- internal_error_if_immutable(dest) or_return
- /*
- Grow `dest` to fit `src`.
- If `dest` is not yet initialized, it will be using `allocator`.
- */
- needed := src.used if minimize else max(src.used, _DEFAULT_DIGIT_COUNT)
- internal_grow(dest, needed, minimize) or_return
- /*
- Copy everything over and zero high digits.
- */
- internal_copy_digits(dest, src, src.used)
- dest.used = src.used
- dest.sign = src.sign
- dest.flags = src.flags &~ {.Immutable}
- internal_zero_unused(dest)
- return nil
- }
- internal_copy :: proc { internal_int_copy, }
- /*
- In normal code, you can also write `a, b = b, a`.
- However, that only swaps within the current scope.
- This helper swaps completely.
- */
- internal_int_swap :: #force_inline proc(a, b: ^Int) {
- a := a; b := b
- a.used, b.used = b.used, a.used
- a.sign, b.sign = b.sign, a.sign
- a.digit, b.digit = b.digit, a.digit
- }
- internal_swap :: proc { internal_int_swap, }
- /*
- Set `dest` to |`src`|.
- */
- internal_int_abs :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- /*
- If `dest == src`, just fix `dest`'s sign.
- */
- if (dest == src) {
- dest.sign = .Zero_or_Positive
- return nil
- }
- /*
- Copy `src` to `dest`
- */
- internal_copy(dest, src) or_return
- /*
- Fix sign.
- */
- dest.sign = .Zero_or_Positive
- return nil
- }
- internal_platform_abs :: proc(n: $T) -> T where intrinsics.type_is_integer(T) {
- return n if n >= 0 else -n
- }
- internal_abs :: proc{ internal_int_abs, internal_platform_abs, }
- /*
- Set `dest` to `-src`.
- */
- internal_int_neg :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- /*
- If `dest == src`, just fix `dest`'s sign.
- */
- sign := Sign.Negative
- if #force_inline internal_is_zero(src) || #force_inline internal_is_negative(src) {
- sign = .Zero_or_Positive
- }
- if dest == src {
- dest.sign = sign
- return nil
- }
- /*
- Copy `src` to `dest`
- */
- internal_copy(dest, src) or_return
- /*
- Fix sign.
- */
- dest.sign = sign
- return nil
- }
- internal_neg :: proc { internal_int_neg, }
- /*
- hac 14.61, pp608.
- */
- internal_int_inverse_modulo :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- /*
- For all n in N and n > 0, n = 0 mod 1.
- */
- if internal_is_positive(a) && internal_eq(b, 1) { return internal_zero(dest) }
- /*
- `b` cannot be negative and has to be > 1
- */
- if internal_is_negative(b) || internal_gt(b, 1) { return .Invalid_Argument }
- /*
- If the modulus is odd we can use a faster routine instead.
- */
- if internal_is_odd(b) { return _private_inverse_modulo_odd(dest, a, b) }
- return _private_inverse_modulo(dest, a, b)
- }
- internal_invmod :: proc{ internal_int_inverse_modulo, }
- /*
- Helpers to extract values from the `Int`.
- Offset is zero indexed.
- */
- internal_int_bitfield_extract_bool :: proc(a: ^Int, offset: int) -> (val: bool, err: Error) {
- limb := offset / _DIGIT_BITS
- if limb < 0 || limb >= a.used { return false, .Invalid_Argument }
- i := _WORD(1 << _WORD((offset % _DIGIT_BITS)))
- return bool(_WORD(a.digit[limb]) & i), nil
- }
- internal_int_bitfield_extract_single :: proc(a: ^Int, offset: int) -> (bit: _WORD, err: Error) {
- limb := offset / _DIGIT_BITS
- if limb < 0 || limb >= a.used { return 0, .Invalid_Argument }
- i := _WORD(1 << _WORD((offset % _DIGIT_BITS)))
- return 1 if ((_WORD(a.digit[limb]) & i) != 0) else 0, nil
- }
- internal_int_bitfield_extract :: proc(a: ^Int, offset, count: int) -> (res: _WORD, err: Error) #no_bounds_check {
- /*
- Early out for single bit.
- */
- if count == 1 {
- limb := offset / _DIGIT_BITS
- if limb < 0 || limb >= a.used { return 0, .Invalid_Argument }
- i := _WORD(1 << _WORD((offset % _DIGIT_BITS)))
- return 1 if ((_WORD(a.digit[limb]) & i) != 0) else 0, nil
- }
- if count > _WORD_BITS || count < 1 { return 0, .Invalid_Argument }
- /*
- There are 3 possible cases.
- - [offset:][:count] covers 1 DIGIT,
- e.g. offset: 0, count: 60 = bits 0..59
- - [offset:][:count] covers 2 DIGITS,
- e.g. offset: 5, count: 60 = bits 5..59, 0..4
- e.g. offset: 0, count: 120 = bits 0..59, 60..119
- - [offset:][:count] covers 3 DIGITS,
- e.g. offset: 40, count: 100 = bits 40..59, 0..59, 0..19
- e.g. offset: 40, count: 120 = bits 40..59, 0..59, 0..39
- */
- limb := offset / _DIGIT_BITS
- bits_left := count
- bits_offset := offset % _DIGIT_BITS
- num_bits := min(bits_left, _DIGIT_BITS - bits_offset)
- shift := offset % _DIGIT_BITS
- mask := (_WORD(1) << uint(num_bits)) - 1
- res = (_WORD(a.digit[limb]) >> uint(shift)) & mask
- bits_left -= num_bits
- if bits_left == 0 { return res, nil }
- res_shift := num_bits
- num_bits = min(bits_left, _DIGIT_BITS)
- mask = (1 << uint(num_bits)) - 1
- res |= (_WORD(a.digit[limb + 1]) & mask) << uint(res_shift)
- bits_left -= num_bits
- if bits_left == 0 { return res, nil }
- mask = (1 << uint(bits_left)) - 1
- res_shift += _DIGIT_BITS
- res |= (_WORD(a.digit[limb + 2]) & mask) << uint(res_shift)
- return res, nil
- }
- /*
- Helpers to (un)set a bit in an Int.
- Offset is zero indexed.
- */
- internal_int_bitfield_set_single :: proc(a: ^Int, offset: int) -> (err: Error) {
- limb := offset / _DIGIT_BITS
- if limb < 0 || limb >= a.used { return .Invalid_Argument }
- i := DIGIT(1 << uint((offset % _DIGIT_BITS)))
- a.digit[limb] |= i
- return
- }
- internal_int_bitfield_unset_single :: proc(a: ^Int, offset: int) -> (err: Error) {
- limb := offset / _DIGIT_BITS
- if limb < 0 || limb >= a.used { return .Invalid_Argument }
- i := DIGIT(1 << uint((offset % _DIGIT_BITS)))
- a.digit[limb] &= _MASK - i
- return
- }
- internal_int_bitfield_toggle_single :: proc(a: ^Int, offset: int) -> (err: Error) {
- limb := offset / _DIGIT_BITS
- if limb < 0 || limb >= a.used { return .Invalid_Argument }
- i := DIGIT(1 << uint((offset % _DIGIT_BITS)))
- a.digit[limb] ~= i
- return
- }
- /*
- Resize backing store.
- We don't need to pass the allocator, because the storage itself stores it.
- Assumes `a` not to be `nil`, and to have already been initialized.
- */
- internal_int_shrink :: proc(a: ^Int) -> (err: Error) {
- needed := max(_MIN_DIGIT_COUNT, a.used)
- if a.used != needed { return internal_grow(a, needed, true) }
- return nil
- }
- internal_shrink :: proc { internal_int_shrink, }
- internal_int_grow :: proc(a: ^Int, digits: int, allow_shrink := false, allocator := context.allocator) -> (err: Error) {
- /*
- We need at least _MIN_DIGIT_COUNT or a.used digits, whichever is bigger.
- The caller is asking for `digits`. Let's be accomodating.
- */
- cap := internal_int_allocated_cap(a)
- needed := max(_MIN_DIGIT_COUNT, a.used, digits)
- if !allow_shrink {
- needed = max(needed, cap)
- }
- /*
- If not yet iniialized, initialize the `digit` backing with the allocator we were passed.
- */
- if cap == 0 {
- a.digit = make([dynamic]DIGIT, needed, allocator)
- } else if cap != needed {
- /*
- `[dynamic]DIGIT` already knows what allocator was used for it, so resize will do the right thing.
- */
- resize(&a.digit, needed)
- }
- /*
- Let's see if the allocation/resize worked as expected.
- */
- if len(a.digit) != needed {
- return .Out_Of_Memory
- }
- return nil
- }
- internal_grow :: proc { internal_int_grow, }
- /*
- Clear `Int` and resize it to the default size.
- Assumes `a` not to be `nil`.
- */
- internal_int_clear :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
- raw := transmute(mem.Raw_Dynamic_Array)a.digit
- if raw.cap != 0 {
- mem.zero_slice(a.digit[:a.used])
- }
- a.sign = .Zero_or_Positive
- a.used = 0
- return #force_inline internal_grow(a, a.used, minimize, allocator)
- }
- internal_clear :: proc { internal_int_clear, }
- internal_zero :: internal_clear
- /*
- Set the `Int` to 1 and optionally shrink it to the minimum backing size.
- */
- internal_int_one :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
- return internal_copy(a, INT_ONE, minimize, allocator)
- }
- internal_one :: proc { internal_int_one, }
- /*
- Set the `Int` to -1 and optionally shrink it to the minimum backing size.
- */
- internal_int_minus_one :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
- return internal_copy(a, INT_MINUS_ONE, minimize, allocator)
- }
- internal_minus_one :: proc { internal_int_minus_one, }
- /*
- Set the `Int` to Inf and optionally shrink it to the minimum backing size.
- */
- internal_int_inf :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
- return internal_copy(a, INT_INF, minimize, allocator)
- }
- internal_inf :: proc { internal_int_inf, }
- /*
- Set the `Int` to -Inf and optionally shrink it to the minimum backing size.
- */
- internal_int_minus_inf :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
- return internal_copy(a, INT_MINUS_INF, minimize, allocator)
- }
- internal_minus_inf :: proc { internal_int_inf, }
- /*
- Set the `Int` to NaN and optionally shrink it to the minimum backing size.
- */
- internal_int_nan :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
- return internal_copy(a, INT_NAN, minimize, allocator)
- }
- internal_nan :: proc { internal_int_nan, }
- internal_int_power_of_two :: proc(a: ^Int, power: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- if power < 0 || power > _MAX_BIT_COUNT { return .Invalid_Argument }
- /*
- Grow to accomodate the single bit.
- */
- a.used = (power / _DIGIT_BITS) + 1
- internal_grow(a, a.used) or_return
- /*
- Zero the entirety.
- */
- mem.zero_slice(a.digit[:])
- /*
- Set the bit.
- */
- a.digit[power / _DIGIT_BITS] = 1 << uint((power % _DIGIT_BITS))
- return nil
- }
- internal_int_get_u128 :: proc(a: ^Int) -> (res: u128, err: Error) {
- return internal_int_get(a, u128)
- }
- internal_get_u128 :: proc { internal_int_get_u128, }
- internal_int_get_i128 :: proc(a: ^Int) -> (res: i128, err: Error) {
- return internal_int_get(a, i128)
- }
- internal_get_i128 :: proc { internal_int_get_i128, }
- internal_int_get_u64 :: proc(a: ^Int) -> (res: u64, err: Error) {
- return internal_int_get(a, u64)
- }
- internal_get_u64 :: proc { internal_int_get_u64, }
- internal_int_get_i64 :: proc(a: ^Int) -> (res: i64, err: Error) {
- return internal_int_get(a, i64)
- }
- internal_get_i64 :: proc { internal_int_get_i64, }
- internal_int_get_u32 :: proc(a: ^Int) -> (res: u32, err: Error) {
- return internal_int_get(a, u32)
- }
- internal_get_u32 :: proc { internal_int_get_u32, }
- internal_int_get_i32 :: proc(a: ^Int) -> (res: i32, err: Error) {
- return internal_int_get(a, i32)
- }
- internal_get_i32 :: proc { internal_int_get_i32, }
- /*
- TODO: Think about using `count_bits` to check if the value could be returned completely,
- and maybe return max(T), .Integer_Overflow if not?
- */
- internal_int_get :: proc(a: ^Int, $T: typeid) -> (res: T, err: Error) where intrinsics.type_is_integer(T) {
- size_in_bits := int(size_of(T) * 8)
- i := int((size_in_bits + _DIGIT_BITS - 1) / _DIGIT_BITS)
- i = min(int(a.used), i)
- #no_bounds_check for ; i >= 0; i -= 1 {
- res <<= uint(0) if size_in_bits <= _DIGIT_BITS else _DIGIT_BITS
- res |= T(a.digit[i])
- if size_in_bits <= _DIGIT_BITS {
- break
- }
- }
- when !intrinsics.type_is_unsigned(T) {
- /*
- Mask off sign bit.
- */
- res ~= 1 << uint(size_in_bits - 1)
- /*
- Set the sign.
- */
- if a.sign == .Negative { res = -res }
- }
- return
- }
- internal_get :: proc { internal_int_get, }
- internal_int_get_float :: proc(a: ^Int) -> (res: f64, err: Error) {
- /*
- log2(max(f64)) is approximately 1020, or 17 legs with the 64-bit storage.
- */
- legs :: 1020 / _DIGIT_BITS
- l := min(a.used, legs)
- fac := f64(1 << _DIGIT_BITS)
- d := 0.0
- #no_bounds_check for i := l; i >= 0; i -= 1 {
- d = (d * fac) + f64(a.digit[i])
- }
- res = -d if a.sign == .Negative else d
- return
- }
- /*
- The `and`, `or` and `xor` binops differ in two lines only.
- We could handle those with a switch, but that adds overhead.
- TODO: Implement versions that take a DIGIT immediate.
- */
- /*
- 2's complement `and`, returns `dest = a & b;`
- */
- internal_int_and :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- used := max(a.used, b.used) + 1
- /*
- Grow the destination to accomodate the result.
- */
- internal_grow(dest, used) or_return
- neg_a := #force_inline internal_is_negative(a)
- neg_b := #force_inline internal_is_negative(b)
- neg := neg_a && neg_b
- ac, bc, cc := DIGIT(1), DIGIT(1), DIGIT(1)
- #no_bounds_check for i := 0; i < used; i += 1 {
- x, y: DIGIT
- /*
- Convert to 2's complement if negative.
- */
- if neg_a {
- ac += _MASK if i >= a.used else (~a.digit[i] & _MASK)
- x = ac & _MASK
- ac >>= _DIGIT_BITS
- } else {
- x = 0 if i >= a.used else a.digit[i]
- }
- /*
- Convert to 2's complement if negative.
- */
- if neg_b {
- bc += _MASK if i >= b.used else (~b.digit[i] & _MASK)
- y = bc & _MASK
- bc >>= _DIGIT_BITS
- } else {
- y = 0 if i >= b.used else b.digit[i]
- }
- dest.digit[i] = x & y
- /*
- Convert to to sign-magnitude if negative.
- */
- if neg {
- cc += ~dest.digit[i] & _MASK
- dest.digit[i] = cc & _MASK
- cc >>= _DIGIT_BITS
- }
- }
- dest.used = used
- dest.sign = .Negative if neg else .Zero_or_Positive
- return internal_clamp(dest)
- }
- internal_and :: proc { internal_int_and, }
- /*
- 2's complement `or`, returns `dest = a | b;`
- */
- internal_int_or :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- used := max(a.used, b.used) + 1
- /*
- Grow the destination to accomodate the result.
- */
- internal_grow(dest, used) or_return
- neg_a := #force_inline internal_is_negative(a)
- neg_b := #force_inline internal_is_negative(b)
- neg := neg_a || neg_b
- ac, bc, cc := DIGIT(1), DIGIT(1), DIGIT(1)
- #no_bounds_check for i := 0; i < used; i += 1 {
- x, y: DIGIT
- /*
- Convert to 2's complement if negative.
- */
- if neg_a {
- ac += _MASK if i >= a.used else (~a.digit[i] & _MASK)
- x = ac & _MASK
- ac >>= _DIGIT_BITS
- } else {
- x = 0 if i >= a.used else a.digit[i]
- }
- /*
- Convert to 2's complement if negative.
- */
- if neg_b {
- bc += _MASK if i >= b.used else (~b.digit[i] & _MASK)
- y = bc & _MASK
- bc >>= _DIGIT_BITS
- } else {
- y = 0 if i >= b.used else b.digit[i]
- }
- dest.digit[i] = x | y
- /*
- Convert to to sign-magnitude if negative.
- */
- if neg {
- cc += ~dest.digit[i] & _MASK
- dest.digit[i] = cc & _MASK
- cc >>= _DIGIT_BITS
- }
- }
- dest.used = used
- dest.sign = .Negative if neg else .Zero_or_Positive
- return internal_clamp(dest)
- }
- internal_or :: proc { internal_int_or, }
- /*
- 2's complement `xor`, returns `dest = a ~ b;`
- */
- internal_int_xor :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- used := max(a.used, b.used) + 1
- /*
- Grow the destination to accomodate the result.
- */
- internal_grow(dest, used) or_return
- neg_a := #force_inline internal_is_negative(a)
- neg_b := #force_inline internal_is_negative(b)
- neg := neg_a != neg_b
- ac, bc, cc := DIGIT(1), DIGIT(1), DIGIT(1)
- #no_bounds_check for i := 0; i < used; i += 1 {
- x, y: DIGIT
- /*
- Convert to 2's complement if negative.
- */
- if neg_a {
- ac += _MASK if i >= a.used else (~a.digit[i] & _MASK)
- x = ac & _MASK
- ac >>= _DIGIT_BITS
- } else {
- x = 0 if i >= a.used else a.digit[i]
- }
- /*
- Convert to 2's complement if negative.
- */
- if neg_b {
- bc += _MASK if i >= b.used else (~b.digit[i] & _MASK)
- y = bc & _MASK
- bc >>= _DIGIT_BITS
- } else {
- y = 0 if i >= b.used else b.digit[i]
- }
- dest.digit[i] = x ~ y
- /*
- Convert to to sign-magnitude if negative.
- */
- if neg {
- cc += ~dest.digit[i] & _MASK
- dest.digit[i] = cc & _MASK
- cc >>= _DIGIT_BITS
- }
- }
- dest.used = used
- dest.sign = .Negative if neg else .Zero_or_Positive
- return internal_clamp(dest)
- }
- internal_xor :: proc { internal_int_xor, }
- /*
- dest = ~src
- */
- internal_int_complement :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- /*
- Temporarily fix sign.
- */
- old_sign := src.sign
- neg := #force_inline internal_is_zero(src) || #force_inline internal_is_positive(src)
- src.sign = .Negative if neg else .Zero_or_Positive
- err = #force_inline internal_sub(dest, src, 1)
- /*
- Restore sign.
- */
- src.sign = old_sign
- return err
- }
- internal_complement :: proc { internal_int_complement, }
- /*
- quotient, remainder := numerator >> bits;
- `remainder` is allowed to be passed a `nil`, in which case `mod` won't be computed.
- */
- internal_int_shrmod :: proc(quotient, remainder, numerator: ^Int, bits: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- bits := bits
- if bits < 0 { return .Invalid_Argument }
- internal_copy(quotient, numerator) or_return
- /*
- Shift right by a certain bit count (store quotient and optional remainder.)
- `numerator` should not be used after this.
- */
- if remainder != nil {
- internal_int_mod_bits(remainder, numerator, bits) or_return
- }
- /*
- Shift by as many digits in the bit count.
- */
- if bits >= _DIGIT_BITS {
- internal_shr_digit(quotient, bits / _DIGIT_BITS) or_return
- }
- /*
- Shift any bit count < _DIGIT_BITS.
- */
- bits %= _DIGIT_BITS
- if bits != 0 {
- mask := DIGIT(1 << uint(bits)) - 1
- shift := DIGIT(_DIGIT_BITS - bits)
- carry := DIGIT(0)
- #no_bounds_check for x := quotient.used - 1; x >= 0; x -= 1 {
- /*
- Get the lower bits of this word in a temp.
- */
- fwd_carry := quotient.digit[x] & mask
- /*
- Shift the current word and mix in the carry bits from the previous word.
- */
- quotient.digit[x] = (quotient.digit[x] >> uint(bits)) | (carry << shift)
- /*
- Update carry from forward carry.
- */
- carry = fwd_carry
- }
- }
- return internal_clamp(numerator)
- }
- internal_shrmod :: proc { internal_int_shrmod, }
- internal_int_shr :: proc(dest, source: ^Int, bits: int, allocator := context.allocator) -> (err: Error) {
- return #force_inline internal_shrmod(dest, nil, source, bits, allocator)
- }
- internal_shr :: proc { internal_int_shr, }
- /*
- Shift right by `digits` * _DIGIT_BITS bits.
- */
- internal_int_shr_digit :: proc(quotient: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- if digits <= 0 { return nil }
- /*
- If digits > used simply zero and return.
- */
- if digits > quotient.used { return internal_zero(quotient) }
- /*
- Much like `int_shl_digit`, this is implemented using a sliding window,
- except the window goes the other way around.
- b-2 | b-1 | b0 | b1 | b2 | ... | bb | ---->
- /\ | ---->
- \-------------------/ ---->
- */
- #no_bounds_check for x := 0; x < (quotient.used - digits); x += 1 {
- quotient.digit[x] = quotient.digit[x + digits]
- }
- quotient.used -= digits
- internal_zero_unused(quotient)
- return internal_clamp(quotient)
- }
- internal_shr_digit :: proc { internal_int_shr_digit, }
- /*
- Shift right by a certain bit count with sign extension.
- */
- internal_int_shr_signed :: proc(dest, src: ^Int, bits: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- if src.sign == .Zero_or_Positive {
- return internal_shr(dest, src, bits)
- }
- internal_int_add_digit(dest, src, DIGIT(1)) or_return
- internal_shr(dest, dest, bits) or_return
- return internal_sub(dest, src, DIGIT(1))
- }
- internal_shr_signed :: proc { internal_int_shr_signed, }
- /*
- Shift left by a certain bit count.
- */
- internal_int_shl :: proc(dest, src: ^Int, bits: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- bits := bits
- if bits < 0 { return .Invalid_Argument }
- internal_copy(dest, src) or_return
- /*
- Grow `dest` to accommodate the additional bits.
- */
- digits_needed := dest.used + (bits / _DIGIT_BITS) + 1
- internal_grow(dest, digits_needed) or_return
- dest.used = digits_needed
- /*
- Shift by as many digits in the bit count as we have.
- */
- if bits >= _DIGIT_BITS {
- internal_shl_digit(dest, bits / _DIGIT_BITS) or_return
- }
- /*
- Shift any remaining bit count < _DIGIT_BITS
- */
- bits %= _DIGIT_BITS
- if bits != 0 {
- mask := (DIGIT(1) << uint(bits)) - DIGIT(1)
- shift := DIGIT(_DIGIT_BITS - bits)
- carry := DIGIT(0)
- #no_bounds_check for x:= 0; x < dest.used; x+= 1 {
- fwd_carry := (dest.digit[x] >> shift) & mask
- dest.digit[x] = (dest.digit[x] << uint(bits) | carry) & _MASK
- carry = fwd_carry
- }
- /*
- Use final carry.
- */
- if carry != 0 {
- dest.digit[dest.used] = carry
- dest.used += 1
- }
- }
- return internal_clamp(dest)
- }
- internal_shl :: proc { internal_int_shl, }
- /*
- Shift left by `digits` * _DIGIT_BITS bits.
- */
- internal_int_shl_digit :: proc(quotient: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- if digits <= 0 { return nil }
- /*
- No need to shift a zero.
- */
- if #force_inline internal_is_zero(quotient) {
- return nil
- }
- /*
- Resize `quotient` to accomodate extra digits.
- */
- #force_inline internal_grow(quotient, quotient.used + digits) or_return
- /*
- Increment the used by the shift amount then copy upwards.
- */
- /*
- Much like `int_shr_digit`, this is implemented using a sliding window,
- except the window goes the other way around.
- */
- #no_bounds_check for x := quotient.used; x > 0; x -= 1 {
- quotient.digit[x+digits-1] = quotient.digit[x-1]
- }
- quotient.used += digits
- mem.zero_slice(quotient.digit[:digits])
- return nil
- }
- internal_shl_digit :: proc { internal_int_shl_digit, }
- /*
- Count bits in an `Int`.
- Assumes `a` not to be `nil` and to have been initialized.
- */
- internal_count_bits :: proc(a: ^Int) -> (count: int) {
- /*
- Fast path for zero.
- */
- if #force_inline internal_is_zero(a) { return {} }
- /*
- Get the number of DIGITs and use it.
- */
- count = (a.used - 1) * _DIGIT_BITS
- /*
- Take the last DIGIT and count the bits in it.
- */
- clz := int(intrinsics.count_leading_zeros(a.digit[a.used - 1]))
- count += (_DIGIT_TYPE_BITS - clz)
- return
- }
- /*
- Returns the number of trailing zeroes before the first one.
- Differs from regular `ctz` in that 0 returns 0.
- Assumes `a` not to be `nil` and have been initialized.
- */
- internal_int_count_lsb :: proc(a: ^Int) -> (count: int, err: Error) {
- /*
- Easy out.
- */
- if #force_inline internal_is_zero(a) { return {}, nil }
- /*
- Scan lower digits until non-zero.
- */
- x: int
- #no_bounds_check for x = 0; x < a.used && a.digit[x] == 0; x += 1 {}
- when true {
- q := a.digit[x]
- x *= _DIGIT_BITS
- x += internal_count_lsb(q)
- } else {
- lnz := []int{
- 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
- }
- q := a.digit[x]
- x *= _DIGIT_BITS
- if q & 1 == 0 {
- p: DIGIT
- for {
- p = q & 15
- x += lnz[p]
- q >>= 4
- if p != 0 { break }
- }
- }
- }
- return x, nil
- }
- internal_platform_count_lsb :: #force_inline proc(a: $T) -> (count: int)
- where intrinsics.type_is_integer(T) && intrinsics.type_is_unsigned(T) {
- return int(intrinsics.count_trailing_zeros(a)) if a > 0 else 0
- }
- internal_count_lsb :: proc { internal_int_count_lsb, internal_platform_count_lsb, }
- internal_int_random_digit :: proc(r: ^rnd.Rand = nil) -> (res: DIGIT) {
- when _DIGIT_BITS == 60 { // DIGIT = u64
- return DIGIT(rnd.uint64(r)) & _MASK
- } else when _DIGIT_BITS == 28 { // DIGIT = u32
- return DIGIT(rnd.uint32(r)) & _MASK
- } else {
- panic("Unsupported DIGIT size.")
- }
- return 0 // We shouldn't get here.
- }
- internal_int_random :: proc(dest: ^Int, bits: int, r: ^rnd.Rand = nil, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- bits := bits
- if bits <= 0 { return .Invalid_Argument }
- digits := bits / _DIGIT_BITS
- bits %= _DIGIT_BITS
- if bits > 0 {
- digits += 1
- }
- #force_inline internal_grow(dest, digits) or_return
- for i := 0; i < digits; i += 1 {
- dest.digit[i] = int_random_digit(r) & _MASK
- }
- if bits > 0 {
- dest.digit[digits - 1] &= ((1 << uint(bits)) - 1)
- }
- dest.used = digits
- return nil
- }
- internal_random :: proc { internal_int_random, }
- /*
- Internal helpers.
- */
- internal_assert_initialized :: proc(a: ^Int, loc := #caller_location) {
- assert(internal_is_initialized(a), "`Int` was not properly initialized.", loc)
- }
- internal_clear_if_uninitialized_single :: proc(arg: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- if ! #force_inline internal_is_initialized(arg) {
- return #force_inline internal_grow(arg, _DEFAULT_DIGIT_COUNT)
- }
- return err
- }
- internal_clear_if_uninitialized_multi :: proc(args: ..^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- for i in args {
- if ! #force_inline internal_is_initialized(i) {
- e := #force_inline internal_grow(i, _DEFAULT_DIGIT_COUNT)
- if e != nil { err = e }
- }
- }
- return err
- }
- internal_clear_if_uninitialized :: proc {internal_clear_if_uninitialized_single, internal_clear_if_uninitialized_multi, }
- internal_error_if_immutable_single :: proc(arg: ^Int) -> (err: Error) {
- if arg != nil && .Immutable in arg.flags { return .Assignment_To_Immutable }
- return nil
- }
- internal_error_if_immutable_multi :: proc(args: ..^Int) -> (err: Error) {
- for i in args {
- if i != nil && .Immutable in i.flags { return .Assignment_To_Immutable }
- }
- return nil
- }
- internal_error_if_immutable :: proc {internal_error_if_immutable_single, internal_error_if_immutable_multi, }
- /*
- Allocates several `Int`s at once.
- */
- internal_int_init_multi :: proc(integers: ..^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator
- integers := integers
- for a in &integers {
- internal_clear(a) or_return
- }
- return nil
- }
- internal_init_multi :: proc { internal_int_init_multi, }
- /*
- Trim unused digits.
- This is used to ensure that leading zero digits are trimmed and the leading "used" digit will be non-zero.
- Typically very fast. Also fixes the sign if there are no more leading digits.
- */
- internal_clamp :: proc(a: ^Int) -> (err: Error) {
- for a.used > 0 && a.digit[a.used - 1] == 0 { a.used -= 1 }
- if #force_inline internal_is_zero(a) { a.sign = .Zero_or_Positive }
- return nil
- }
- internal_int_zero_unused :: #force_inline proc(dest: ^Int, old_used := -1) {
- /*
- If we don't pass the number of previously used DIGITs, we zero all remaining ones.
- */
- zero_count: int
- if old_used == -1 {
- zero_count = len(dest.digit) - dest.used
- } else {
- zero_count = old_used - dest.used
- }
- /*
- Zero remainder.
- */
- if zero_count > 0 && dest.used < len(dest.digit) {
- mem.zero_slice(dest.digit[dest.used:][:zero_count])
- }
- }
- internal_zero_unused :: proc { internal_int_zero_unused, }
- /*
- ========================== End of low-level routines ==========================
- */
|