bid128_to_uint32.c 127 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853285428552856285728582859286028612862286328642865286628672868286928702871287228732874287528762877287828792880288128822883288428852886288728882889289028912892289328942895289628972898289929002901290229032904290529062907290829092910291129122913291429152916291729182919292029212922292329242925292629272928292929302931293229332934293529362937293829392940294129422943294429452946294729482949295029512952295329542955295629572958295929602961296229632964296529662967296829692970297129722973297429752976297729782979298029812982298329842985298629872988298929902991299229932994299529962997299829993000300130023003300430053006300730083009301030113012301330143015301630173018301930203021302230233024302530263027302830293030303130323033303430353036303730383039304030413042304330443045304630473048304930503051305230533054305530563057305830593060306130623063306430653066306730683069307030713072307330743075307630773078307930803081308230833084308530863087308830893090309130923093309430953096309730983099310031013102310331043105310631073108310931103111311231133114311531163117311831193120312131223123312431253126312731283129313031313132313331343135313631373138313931403141314231433144314531463147314831493150315131523153315431553156315731583159316031613162316331643165316631673168316931703171317231733174317531763177317831793180318131823183318431853186318731883189319031913192319331943195319631973198319932003201320232033204320532063207320832093210321132123213321432153216321732183219322032213222322332243225322632273228322932303231323232333234323532363237323832393240324132423243324432453246324732483249325032513252325332543255325632573258325932603261326232633264326532663267326832693270327132723273327432753276327732783279328032813282328332843285328632873288328932903291329232933294329532963297329832993300330133023303330433053306330733083309331033113312331333143315331633173318331933203321332233233324332533263327332833293330333133323333333433353336333733383339334033413342334333443345334633473348334933503351335233533354335533563357335833593360336133623363336433653366336733683369337033713372337333743375337633773378337933803381338233833384338533863387338833893390339133923393339433953396339733983399340034013402340334043405340634073408340934103411341234133414341534163417341834193420342134223423342434253426342734283429343034313432343334343435343634373438343934403441344234433444344534463447344834493450345134523453345434553456345734583459346034613462346334643465346634673468346934703471347234733474347534763477347834793480348134823483348434853486348734883489349034913492349334943495349634973498349935003501350235033504350535063507350835093510351135123513351435153516351735183519352035213522352335243525352635273528352935303531353235333534353535363537353835393540354135423543354435453546354735483549355035513552355335543555355635573558355935603561356235633564356535663567356835693570357135723573357435753576357735783579358035813582358335843585358635873588
  1. /* Copyright (C) 2007-2022 Free Software Foundation, Inc.
  2. This file is part of GCC.
  3. GCC is free software; you can redistribute it and/or modify it under
  4. the terms of the GNU General Public License as published by the Free
  5. Software Foundation; either version 3, or (at your option) any later
  6. version.
  7. GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  8. WARRANTY; without even the implied warranty of MERCHANTABILITY or
  9. FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
  10. for more details.
  11. Under Section 7 of GPL version 3, you are granted additional
  12. permissions described in the GCC Runtime Library Exception, version
  13. 3.1, as published by the Free Software Foundation.
  14. You should have received a copy of the GNU General Public License and
  15. a copy of the GCC Runtime Library Exception along with this program;
  16. see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
  17. <http://www.gnu.org/licenses/>. */
  18. #include "bid_internal.h"
  19. /*****************************************************************************
  20. * BID128_to_uint32_rnint
  21. ****************************************************************************/
  22. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
  23. bid128_to_uint32_rnint, x)
  24. unsigned int res;
  25. UINT64 x_sign;
  26. UINT64 x_exp;
  27. int exp; // unbiased exponent
  28. // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  29. UINT64 tmp64;
  30. BID_UI64DOUBLE tmp1;
  31. unsigned int x_nr_bits;
  32. int q, ind, shift;
  33. UINT128 C1, C;
  34. UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
  35. UINT256 fstar;
  36. UINT256 P256;
  37. // unpack x
  38. x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
  39. x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
  40. C1.w[1] = x.w[1] & MASK_COEFF;
  41. C1.w[0] = x.w[0];
  42. // check for NaN or Infinity
  43. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  44. // x is special
  45. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  46. if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
  47. // set invalid flag
  48. *pfpsf |= INVALID_EXCEPTION;
  49. // return Integer Indefinite
  50. res = 0x80000000;
  51. } else { // x is QNaN
  52. // set invalid flag
  53. *pfpsf |= INVALID_EXCEPTION;
  54. // return Integer Indefinite
  55. res = 0x80000000;
  56. }
  57. BID_RETURN (res);
  58. } else { // x is not a NaN, so it must be infinity
  59. if (!x_sign) { // x is +inf
  60. // set invalid flag
  61. *pfpsf |= INVALID_EXCEPTION;
  62. // return Integer Indefinite
  63. res = 0x80000000;
  64. } else { // x is -inf
  65. // set invalid flag
  66. *pfpsf |= INVALID_EXCEPTION;
  67. // return Integer Indefinite
  68. res = 0x80000000;
  69. }
  70. BID_RETURN (res);
  71. }
  72. }
  73. // check for non-canonical values (after the check for special values)
  74. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  75. || (C1.w[1] == 0x0001ed09bead87c0ull
  76. && (C1.w[0] > 0x378d8e63ffffffffull))
  77. || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  78. res = 0x00000000;
  79. BID_RETURN (res);
  80. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  81. // x is 0
  82. res = 0x00000000;
  83. BID_RETURN (res);
  84. } else { // x is not special and is not zero
  85. // q = nr. of decimal digits in x
  86. // determine first the nr. of bits in x
  87. if (C1.w[1] == 0) {
  88. if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
  89. // split the 64-bit value in two 32-bit halves to avoid rounding errors
  90. if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
  91. tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
  92. x_nr_bits =
  93. 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  94. } else { // x < 2^32
  95. tmp1.d = (double) (C1.w[0]); // exact conversion
  96. x_nr_bits =
  97. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  98. }
  99. } else { // if x < 2^53
  100. tmp1.d = (double) C1.w[0]; // exact conversion
  101. x_nr_bits =
  102. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  103. }
  104. } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  105. tmp1.d = (double) C1.w[1]; // exact conversion
  106. x_nr_bits =
  107. 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  108. }
  109. q = nr_digits[x_nr_bits - 1].digits;
  110. if (q == 0) {
  111. q = nr_digits[x_nr_bits - 1].digits1;
  112. if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  113. || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  114. && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  115. q++;
  116. }
  117. exp = (x_exp >> 49) - 6176;
  118. if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
  119. // set invalid flag
  120. *pfpsf |= INVALID_EXCEPTION;
  121. // return Integer Indefinite
  122. res = 0x80000000;
  123. BID_RETURN (res);
  124. } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
  125. // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
  126. // so x rounded to an integer may or may not fit in an unsigned 32-bit int
  127. // the cases that do not fit are identified here; the ones that fit
  128. // fall through and will be handled with other cases further,
  129. // under '1 <= q + exp <= 10'
  130. if (x_sign) { // if n < 0 and q + exp = 10
  131. // if n < -1/2 then n cannot be converted to uint32 with RN
  132. // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 1/2
  133. // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x05, 1<=q<=34
  134. if (q <= 11) {
  135. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  136. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  137. if (tmp64 > 0x05ull) {
  138. // set invalid flag
  139. *pfpsf |= INVALID_EXCEPTION;
  140. // return Integer Indefinite
  141. res = 0x80000000;
  142. BID_RETURN (res);
  143. }
  144. // else cases that can be rounded to a 32-bit int fall through
  145. // to '1 <= q + exp <= 10'
  146. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  147. // 0.c(0)c(1)...c(q-1) * 10^11 > 0x05 <=>
  148. // C > 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
  149. // (scale 1/2 up)
  150. tmp64 = 0x05ull;
  151. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  152. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  153. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  154. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  155. }
  156. if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
  157. // set invalid flag
  158. *pfpsf |= INVALID_EXCEPTION;
  159. // return Integer Indefinite
  160. res = 0x80000000;
  161. BID_RETURN (res);
  162. }
  163. // else cases that can be rounded to a 32-bit int fall through
  164. // to '1 <= q + exp <= 10'
  165. }
  166. } else { // if n > 0 and q + exp = 10
  167. // if n >= 2^32 - 1/2 then n is too large
  168. // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
  169. // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
  170. if (q <= 11) {
  171. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  172. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  173. if (tmp64 >= 0x9fffffffbull) {
  174. // set invalid flag
  175. *pfpsf |= INVALID_EXCEPTION;
  176. // return Integer Indefinite
  177. res = 0x80000000;
  178. BID_RETURN (res);
  179. }
  180. // else cases that can be rounded to a 32-bit int fall through
  181. // to '1 <= q + exp <= 10'
  182. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  183. // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
  184. // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
  185. // (scale 2^32-1/2 up)
  186. tmp64 = 0x9fffffffbull;
  187. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  188. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  189. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  190. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  191. }
  192. if (C1.w[1] > C.w[1]
  193. || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  194. // set invalid flag
  195. *pfpsf |= INVALID_EXCEPTION;
  196. // return Integer Indefinite
  197. res = 0x80000000;
  198. BID_RETURN (res);
  199. }
  200. // else cases that can be rounded to a 32-bit int fall through
  201. // to '1 <= q + exp <= 10'
  202. }
  203. }
  204. }
  205. // n is not too large to be converted to int32: -1/2 <= n < 2^32 - 1/2
  206. // Note: some of the cases tested for above fall through to this point
  207. if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
  208. // return 0
  209. res = 0x00000000;
  210. BID_RETURN (res);
  211. } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
  212. // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
  213. // res = 0
  214. // else if x > 0
  215. // res = +1
  216. // else // if x < 0
  217. // invalid exc
  218. ind = q - 1;
  219. if (ind <= 18) { // 0 <= ind <= 18
  220. if ((C1.w[1] == 0) && (C1.w[0] <= midpoint64[ind])) {
  221. res = 0x00000000; // return 0
  222. } else if (!x_sign) { // n > 0
  223. res = 0x00000001; // return +1
  224. } else {
  225. res = 0x80000000;
  226. *pfpsf |= INVALID_EXCEPTION;
  227. }
  228. } else { // 19 <= ind <= 33
  229. if ((C1.w[1] < midpoint128[ind - 19].w[1])
  230. || ((C1.w[1] == midpoint128[ind - 19].w[1])
  231. && (C1.w[0] <= midpoint128[ind - 19].w[0]))) {
  232. res = 0x00000000; // return 0
  233. } else if (!x_sign) { // n > 0
  234. res = 0x00000001; // return +1
  235. } else {
  236. res = 0x80000000;
  237. *pfpsf |= INVALID_EXCEPTION;
  238. }
  239. }
  240. } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
  241. if (x_sign) { // x <= -1
  242. // set invalid flag
  243. *pfpsf |= INVALID_EXCEPTION;
  244. // return Integer Indefinite
  245. res = 0x80000000;
  246. BID_RETURN (res);
  247. }
  248. // 1 <= x < 2^32-1/2 so x can be rounded
  249. // to nearest to a 32-bit unsigned integer
  250. if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
  251. ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
  252. // chop off ind digits from the lower part of C1
  253. // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
  254. tmp64 = C1.w[0];
  255. if (ind <= 19) {
  256. C1.w[0] = C1.w[0] + midpoint64[ind - 1];
  257. } else {
  258. C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
  259. C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
  260. }
  261. if (C1.w[0] < tmp64)
  262. C1.w[1]++;
  263. // calculate C* and f*
  264. // C* is actually floor(C*) in this case
  265. // C* and f* need shifting and masking, as shown by
  266. // shiftright128[] and maskhigh128[]
  267. // 1 <= x <= 33
  268. // kx = 10^(-x) = ten2mk128[ind - 1]
  269. // C* = (C1 + 1/2 * 10^x) * 10^(-x)
  270. // the approximation of 10^(-x) was rounded up to 118 bits
  271. __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  272. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  273. Cstar.w[1] = P256.w[3];
  274. Cstar.w[0] = P256.w[2];
  275. fstar.w[3] = 0;
  276. fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  277. fstar.w[1] = P256.w[1];
  278. fstar.w[0] = P256.w[0];
  279. } else { // 22 <= ind - 1 <= 33
  280. Cstar.w[1] = 0;
  281. Cstar.w[0] = P256.w[3];
  282. fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  283. fstar.w[2] = P256.w[2];
  284. fstar.w[1] = P256.w[1];
  285. fstar.w[0] = P256.w[0];
  286. }
  287. // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  288. // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  289. // if (0 < f* < 10^(-x)) then the result is a midpoint
  290. // if floor(C*) is even then C* = floor(C*) - logical right
  291. // shift; C* has p decimal digits, correct by Prop. 1)
  292. // else if floor(C*) is odd C* = floor(C*)-1 (logical right
  293. // shift; C* has p decimal digits, correct by Pr. 1)
  294. // else
  295. // C* = floor(C*) (logical right shift; C has p decimal digits,
  296. // correct by Property 1)
  297. // n = C* * 10^(e+x)
  298. // shift right C* by Ex-128 = shiftright128[ind]
  299. shift = shiftright128[ind - 1]; // 0 <= shift <= 102
  300. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  301. Cstar.w[0] =
  302. (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  303. // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  304. } else { // 22 <= ind - 1 <= 33
  305. Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
  306. }
  307. // if the result was a midpoint it was rounded away from zero, so
  308. // it will need a correction
  309. // check for midpoints
  310. if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
  311. && (fstar.w[1] || fstar.w[0])
  312. && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
  313. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  314. && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
  315. // the result is a midpoint; round to nearest
  316. if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
  317. // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
  318. Cstar.w[0]--; // Cstar.w[0] is now even
  319. } // else MP in [ODD, EVEN]
  320. }
  321. res = Cstar.w[0]; // the result is positive
  322. } else if (exp == 0) {
  323. // 1 <= q <= 10
  324. // res = C (exact)
  325. res = C1.w[0];
  326. } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
  327. // res = C * 10^exp (exact)
  328. res = C1.w[0] * ten2k64[exp];
  329. }
  330. }
  331. }
  332. BID_RETURN (res);
  333. }
  334. /*****************************************************************************
  335. * BID128_to_uint32_xrnint
  336. ****************************************************************************/
  337. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
  338. bid128_to_uint32_xrnint, x)
  339. unsigned int res;
  340. UINT64 x_sign;
  341. UINT64 x_exp;
  342. int exp; // unbiased exponent
  343. // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  344. UINT64 tmp64, tmp64A;
  345. BID_UI64DOUBLE tmp1;
  346. unsigned int x_nr_bits;
  347. int q, ind, shift;
  348. UINT128 C1, C;
  349. UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
  350. UINT256 fstar;
  351. UINT256 P256;
  352. unsigned int tmp_inexact = 0;
  353. // unpack x
  354. x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
  355. x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
  356. C1.w[1] = x.w[1] & MASK_COEFF;
  357. C1.w[0] = x.w[0];
  358. // check for NaN or Infinity
  359. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  360. // x is special
  361. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  362. if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
  363. // set invalid flag
  364. *pfpsf |= INVALID_EXCEPTION;
  365. // return Integer Indefinite
  366. res = 0x80000000;
  367. } else { // x is QNaN
  368. // set invalid flag
  369. *pfpsf |= INVALID_EXCEPTION;
  370. // return Integer Indefinite
  371. res = 0x80000000;
  372. }
  373. BID_RETURN (res);
  374. } else { // x is not a NaN, so it must be infinity
  375. if (!x_sign) { // x is +inf
  376. // set invalid flag
  377. *pfpsf |= INVALID_EXCEPTION;
  378. // return Integer Indefinite
  379. res = 0x80000000;
  380. } else { // x is -inf
  381. // set invalid flag
  382. *pfpsf |= INVALID_EXCEPTION;
  383. // return Integer Indefinite
  384. res = 0x80000000;
  385. }
  386. BID_RETURN (res);
  387. }
  388. }
  389. // check for non-canonical values (after the check for special values)
  390. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  391. || (C1.w[1] == 0x0001ed09bead87c0ull
  392. && (C1.w[0] > 0x378d8e63ffffffffull))
  393. || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  394. res = 0x00000000;
  395. BID_RETURN (res);
  396. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  397. // x is 0
  398. res = 0x00000000;
  399. BID_RETURN (res);
  400. } else { // x is not special and is not zero
  401. // q = nr. of decimal digits in x
  402. // determine first the nr. of bits in x
  403. if (C1.w[1] == 0) {
  404. if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
  405. // split the 64-bit value in two 32-bit halves to avoid rounding errors
  406. if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
  407. tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
  408. x_nr_bits =
  409. 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  410. } else { // x < 2^32
  411. tmp1.d = (double) (C1.w[0]); // exact conversion
  412. x_nr_bits =
  413. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  414. }
  415. } else { // if x < 2^53
  416. tmp1.d = (double) C1.w[0]; // exact conversion
  417. x_nr_bits =
  418. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  419. }
  420. } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  421. tmp1.d = (double) C1.w[1]; // exact conversion
  422. x_nr_bits =
  423. 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  424. }
  425. q = nr_digits[x_nr_bits - 1].digits;
  426. if (q == 0) {
  427. q = nr_digits[x_nr_bits - 1].digits1;
  428. if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  429. || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  430. && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  431. q++;
  432. }
  433. exp = (x_exp >> 49) - 6176;
  434. if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
  435. // set invalid flag
  436. *pfpsf |= INVALID_EXCEPTION;
  437. // return Integer Indefinite
  438. res = 0x80000000;
  439. BID_RETURN (res);
  440. } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
  441. // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
  442. // so x rounded to an integer may or may not fit in an unsigned 32-bit int
  443. // the cases that do not fit are identified here; the ones that fit
  444. // fall through and will be handled with other cases further,
  445. // under '1 <= q + exp <= 10'
  446. if (x_sign) { // if n < 0 and q + exp = 10
  447. // if n < -1/2 then n cannot be converted to uint32 with RN
  448. // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 1/2
  449. // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x05, 1<=q<=34
  450. if (q <= 11) {
  451. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  452. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  453. if (tmp64 > 0x05ull) {
  454. // set invalid flag
  455. *pfpsf |= INVALID_EXCEPTION;
  456. // return Integer Indefinite
  457. res = 0x80000000;
  458. BID_RETURN (res);
  459. }
  460. // else cases that can be rounded to a 32-bit int fall through
  461. // to '1 <= q + exp <= 10'
  462. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  463. // 0.c(0)c(1)...c(q-1) * 10^11 > 0x05 <=>
  464. // C > 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
  465. // (scale 1/2 up)
  466. tmp64 = 0x05ull;
  467. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  468. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  469. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  470. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  471. }
  472. if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
  473. // set invalid flag
  474. *pfpsf |= INVALID_EXCEPTION;
  475. // return Integer Indefinite
  476. res = 0x80000000;
  477. BID_RETURN (res);
  478. }
  479. // else cases that can be rounded to a 32-bit int fall through
  480. // to '1 <= q + exp <= 10'
  481. }
  482. } else { // if n > 0 and q + exp = 10
  483. // if n >= 2^32 - 1/2 then n is too large
  484. // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
  485. // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
  486. if (q <= 11) {
  487. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  488. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  489. if (tmp64 >= 0x9fffffffbull) {
  490. // set invalid flag
  491. *pfpsf |= INVALID_EXCEPTION;
  492. // return Integer Indefinite
  493. res = 0x80000000;
  494. BID_RETURN (res);
  495. }
  496. // else cases that can be rounded to a 32-bit int fall through
  497. // to '1 <= q + exp <= 10'
  498. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  499. // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
  500. // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
  501. // (scale 2^32-1/2 up)
  502. tmp64 = 0x9fffffffbull;
  503. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  504. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  505. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  506. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  507. }
  508. if (C1.w[1] > C.w[1]
  509. || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  510. // set invalid flag
  511. *pfpsf |= INVALID_EXCEPTION;
  512. // return Integer Indefinite
  513. res = 0x80000000;
  514. BID_RETURN (res);
  515. }
  516. // else cases that can be rounded to a 32-bit int fall through
  517. // to '1 <= q + exp <= 10'
  518. }
  519. }
  520. }
  521. // n is not too large to be converted to int32: -1/2 <= n < 2^32 - 1/2
  522. // Note: some of the cases tested for above fall through to this point
  523. if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
  524. // set inexact flag
  525. *pfpsf |= INEXACT_EXCEPTION;
  526. // return 0
  527. res = 0x00000000;
  528. BID_RETURN (res);
  529. } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
  530. // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
  531. // res = 0
  532. // else if x > 0
  533. // res = +1
  534. // else // if x < 0
  535. // invalid exc
  536. ind = q - 1;
  537. if (ind <= 18) { // 0 <= ind <= 18
  538. if ((C1.w[1] == 0) && (C1.w[0] <= midpoint64[ind])) {
  539. res = 0x00000000; // return 0
  540. } else if (!x_sign) { // n > 0
  541. res = 0x00000001; // return +1
  542. } else {
  543. res = 0x80000000;
  544. *pfpsf |= INVALID_EXCEPTION;
  545. BID_RETURN (res);
  546. }
  547. } else { // 19 <= ind <= 33
  548. if ((C1.w[1] < midpoint128[ind - 19].w[1])
  549. || ((C1.w[1] == midpoint128[ind - 19].w[1])
  550. && (C1.w[0] <= midpoint128[ind - 19].w[0]))) {
  551. res = 0x00000000; // return 0
  552. } else if (!x_sign) { // n > 0
  553. res = 0x00000001; // return +1
  554. } else {
  555. res = 0x80000000;
  556. *pfpsf |= INVALID_EXCEPTION;
  557. BID_RETURN (res);
  558. }
  559. }
  560. // set inexact flag
  561. *pfpsf |= INEXACT_EXCEPTION;
  562. } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
  563. if (x_sign) { // x <= -1
  564. // set invalid flag
  565. *pfpsf |= INVALID_EXCEPTION;
  566. // return Integer Indefinite
  567. res = 0x80000000;
  568. BID_RETURN (res);
  569. }
  570. // 1 <= x < 2^32-1/2 so x can be rounded
  571. // to nearest to a 32-bit unsigned integer
  572. if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
  573. ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
  574. // chop off ind digits from the lower part of C1
  575. // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
  576. tmp64 = C1.w[0];
  577. if (ind <= 19) {
  578. C1.w[0] = C1.w[0] + midpoint64[ind - 1];
  579. } else {
  580. C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
  581. C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
  582. }
  583. if (C1.w[0] < tmp64)
  584. C1.w[1]++;
  585. // calculate C* and f*
  586. // C* is actually floor(C*) in this case
  587. // C* and f* need shifting and masking, as shown by
  588. // shiftright128[] and maskhigh128[]
  589. // 1 <= x <= 33
  590. // kx = 10^(-x) = ten2mk128[ind - 1]
  591. // C* = (C1 + 1/2 * 10^x) * 10^(-x)
  592. // the approximation of 10^(-x) was rounded up to 118 bits
  593. __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  594. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  595. Cstar.w[1] = P256.w[3];
  596. Cstar.w[0] = P256.w[2];
  597. fstar.w[3] = 0;
  598. fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  599. fstar.w[1] = P256.w[1];
  600. fstar.w[0] = P256.w[0];
  601. } else { // 22 <= ind - 1 <= 33
  602. Cstar.w[1] = 0;
  603. Cstar.w[0] = P256.w[3];
  604. fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  605. fstar.w[2] = P256.w[2];
  606. fstar.w[1] = P256.w[1];
  607. fstar.w[0] = P256.w[0];
  608. }
  609. // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  610. // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  611. // if (0 < f* < 10^(-x)) then the result is a midpoint
  612. // if floor(C*) is even then C* = floor(C*) - logical right
  613. // shift; C* has p decimal digits, correct by Prop. 1)
  614. // else if floor(C*) is odd C* = floor(C*)-1 (logical right
  615. // shift; C* has p decimal digits, correct by Pr. 1)
  616. // else
  617. // C* = floor(C*) (logical right shift; C has p decimal digits,
  618. // correct by Property 1)
  619. // n = C* * 10^(e+x)
  620. // shift right C* by Ex-128 = shiftright128[ind]
  621. shift = shiftright128[ind - 1]; // 0 <= shift <= 102
  622. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  623. Cstar.w[0] =
  624. (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  625. // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  626. } else { // 22 <= ind - 1 <= 33
  627. Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
  628. }
  629. // determine inexactness of the rounding of C*
  630. // if (0 < f* - 1/2 < 10^(-x)) then
  631. // the result is exact
  632. // else // if (f* - 1/2 > T*) then
  633. // the result is inexact
  634. if (ind - 1 <= 2) {
  635. if (fstar.w[1] > 0x8000000000000000ull ||
  636. (fstar.w[1] == 0x8000000000000000ull
  637. && fstar.w[0] > 0x0ull)) {
  638. // f* > 1/2 and the result may be exact
  639. tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
  640. if (tmp64 > ten2mk128trunc[ind - 1].w[1]
  641. || (tmp64 == ten2mk128trunc[ind - 1].w[1]
  642. && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
  643. // set the inexact flag
  644. // *pfpsf |= INEXACT_EXCEPTION;
  645. tmp_inexact = 1;
  646. } // else the result is exact
  647. } else { // the result is inexact; f2* <= 1/2
  648. // set the inexact flag
  649. // *pfpsf |= INEXACT_EXCEPTION;
  650. tmp_inexact = 1;
  651. }
  652. } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
  653. if (fstar.w[3] > 0x0 ||
  654. (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
  655. (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
  656. (fstar.w[1] || fstar.w[0]))) {
  657. // f2* > 1/2 and the result may be exact
  658. // Calculate f2* - 1/2
  659. tmp64 = fstar.w[2] - onehalf128[ind - 1];
  660. tmp64A = fstar.w[3];
  661. if (tmp64 > fstar.w[2])
  662. tmp64A--;
  663. if (tmp64A || tmp64
  664. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  665. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  666. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  667. // set the inexact flag
  668. // *pfpsf |= INEXACT_EXCEPTION;
  669. tmp_inexact = 1;
  670. } // else the result is exact
  671. } else { // the result is inexact; f2* <= 1/2
  672. // set the inexact flag
  673. // *pfpsf |= INEXACT_EXCEPTION;
  674. tmp_inexact = 1;
  675. }
  676. } else { // if 22 <= ind <= 33
  677. if (fstar.w[3] > onehalf128[ind - 1] ||
  678. (fstar.w[3] == onehalf128[ind - 1] &&
  679. (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
  680. // f2* > 1/2 and the result may be exact
  681. // Calculate f2* - 1/2
  682. tmp64 = fstar.w[3] - onehalf128[ind - 1];
  683. if (tmp64 || fstar.w[2]
  684. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  685. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  686. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  687. // set the inexact flag
  688. // *pfpsf |= INEXACT_EXCEPTION;
  689. tmp_inexact = 1;
  690. } // else the result is exact
  691. } else { // the result is inexact; f2* <= 1/2
  692. // set the inexact flag
  693. // *pfpsf |= INEXACT_EXCEPTION;
  694. tmp_inexact = 1;
  695. }
  696. }
  697. // if the result was a midpoint it was rounded away from zero, so
  698. // it will need a correction
  699. // check for midpoints
  700. if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
  701. && (fstar.w[1] || fstar.w[0])
  702. && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
  703. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  704. && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
  705. // the result is a midpoint; round to nearest
  706. if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
  707. // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
  708. Cstar.w[0]--; // Cstar.w[0] is now even
  709. } // else MP in [ODD, EVEN]
  710. }
  711. res = Cstar.w[0]; // the result is positive
  712. if (tmp_inexact)
  713. *pfpsf |= INEXACT_EXCEPTION;
  714. } else if (exp == 0) {
  715. // 1 <= q <= 10
  716. // res = C (exact)
  717. res = C1.w[0];
  718. } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
  719. // res = C * 10^exp (exact)
  720. res = C1.w[0] * ten2k64[exp];
  721. }
  722. }
  723. }
  724. BID_RETURN (res);
  725. }
  726. /*****************************************************************************
  727. * BID128_to_uint32_floor
  728. ****************************************************************************/
  729. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
  730. bid128_to_uint32_floor, x)
  731. unsigned int res;
  732. UINT64 x_sign;
  733. UINT64 x_exp;
  734. int exp; // unbiased exponent
  735. // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  736. UINT64 tmp64, tmp64A;
  737. BID_UI64DOUBLE tmp1;
  738. unsigned int x_nr_bits;
  739. int q, ind, shift;
  740. UINT128 C1, C;
  741. UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
  742. UINT256 fstar;
  743. UINT256 P256;
  744. int is_inexact_gt_midpoint = 0;
  745. int is_midpoint_lt_even = 0;
  746. // unpack x
  747. x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
  748. x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
  749. C1.w[1] = x.w[1] & MASK_COEFF;
  750. C1.w[0] = x.w[0];
  751. // check for NaN or Infinity
  752. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  753. // x is special
  754. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  755. if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
  756. // set invalid flag
  757. *pfpsf |= INVALID_EXCEPTION;
  758. // return Integer Indefinite
  759. res = 0x80000000;
  760. } else { // x is QNaN
  761. // set invalid flag
  762. *pfpsf |= INVALID_EXCEPTION;
  763. // return Integer Indefinite
  764. res = 0x80000000;
  765. }
  766. BID_RETURN (res);
  767. } else { // x is not a NaN, so it must be infinity
  768. if (!x_sign) { // x is +inf
  769. // set invalid flag
  770. *pfpsf |= INVALID_EXCEPTION;
  771. // return Integer Indefinite
  772. res = 0x80000000;
  773. } else { // x is -inf
  774. // set invalid flag
  775. *pfpsf |= INVALID_EXCEPTION;
  776. // return Integer Indefinite
  777. res = 0x80000000;
  778. }
  779. BID_RETURN (res);
  780. }
  781. }
  782. // check for non-canonical values (after the check for special values)
  783. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  784. || (C1.w[1] == 0x0001ed09bead87c0ull
  785. && (C1.w[0] > 0x378d8e63ffffffffull))
  786. || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  787. res = 0x00000000;
  788. BID_RETURN (res);
  789. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  790. // x is 0
  791. res = 0x00000000;
  792. BID_RETURN (res);
  793. } else { // x is not special and is not zero
  794. // x < 0 is invalid
  795. if (x_sign) {
  796. // set invalid flag
  797. *pfpsf |= INVALID_EXCEPTION;
  798. // return Integer Indefinite
  799. res = 0x80000000;
  800. BID_RETURN (res);
  801. }
  802. // x > 0 from this point on
  803. // q = nr. of decimal digits in x
  804. // determine first the nr. of bits in x
  805. if (C1.w[1] == 0) {
  806. if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
  807. // split the 64-bit value in two 32-bit halves to avoid rounding errors
  808. if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
  809. tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
  810. x_nr_bits =
  811. 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  812. } else { // x < 2^32
  813. tmp1.d = (double) (C1.w[0]); // exact conversion
  814. x_nr_bits =
  815. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  816. }
  817. } else { // if x < 2^53
  818. tmp1.d = (double) C1.w[0]; // exact conversion
  819. x_nr_bits =
  820. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  821. }
  822. } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  823. tmp1.d = (double) C1.w[1]; // exact conversion
  824. x_nr_bits =
  825. 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  826. }
  827. q = nr_digits[x_nr_bits - 1].digits;
  828. if (q == 0) {
  829. q = nr_digits[x_nr_bits - 1].digits1;
  830. if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  831. || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  832. && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  833. q++;
  834. }
  835. exp = (x_exp >> 49) - 6176;
  836. if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
  837. // set invalid flag
  838. *pfpsf |= INVALID_EXCEPTION;
  839. // return Integer Indefinite
  840. res = 0x80000000;
  841. BID_RETURN (res);
  842. } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
  843. // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
  844. // so x rounded to an integer may or may not fit in a signed 32-bit int
  845. // the cases that do not fit are identified here; the ones that fit
  846. // fall through and will be handled with other cases further,
  847. // under '1 <= q + exp <= 10'
  848. // n > 0 and q + exp = 10
  849. // if n >= 2^32 then n is too large
  850. // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
  851. // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
  852. if (q <= 11) {
  853. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  854. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  855. if (tmp64 >= 0xa00000000ull) {
  856. // set invalid flag
  857. *pfpsf |= INVALID_EXCEPTION;
  858. // return Integer Indefinite
  859. res = 0x80000000;
  860. BID_RETURN (res);
  861. }
  862. // else cases that can be rounded to a 32-bit int fall through
  863. // to '1 <= q + exp <= 10'
  864. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  865. // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
  866. // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
  867. // (scale 2^32 up)
  868. tmp64 = 0xa00000000ull;
  869. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  870. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  871. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  872. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  873. }
  874. if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  875. // set invalid flag
  876. *pfpsf |= INVALID_EXCEPTION;
  877. // return Integer Indefinite
  878. res = 0x80000000;
  879. BID_RETURN (res);
  880. }
  881. // else cases that can be rounded to a 32-bit int fall through
  882. // to '1 <= q + exp <= 10'
  883. }
  884. }
  885. // n is not too large to be converted to int32: 0 <= n < 2^31
  886. // Note: some of the cases tested for above fall through to this point
  887. if ((q + exp) <= 0) {
  888. // n = +0.0...c(0)c(1)...c(q-1) or n = +0.c(0)c(1)...c(q-1)
  889. // return 0
  890. res = 0x00000000;
  891. BID_RETURN (res);
  892. } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
  893. // 1 <= x < 2^32 so x can be rounded down to a 32-bit unsigned integer
  894. if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
  895. ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
  896. // chop off ind digits from the lower part of C1
  897. // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
  898. tmp64 = C1.w[0];
  899. if (ind <= 19) {
  900. C1.w[0] = C1.w[0] + midpoint64[ind - 1];
  901. } else {
  902. C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
  903. C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
  904. }
  905. if (C1.w[0] < tmp64)
  906. C1.w[1]++;
  907. // calculate C* and f*
  908. // C* is actually floor(C*) in this case
  909. // C* and f* need shifting and masking, as shown by
  910. // shiftright128[] and maskhigh128[]
  911. // 1 <= x <= 33
  912. // kx = 10^(-x) = ten2mk128[ind - 1]
  913. // C* = (C1 + 1/2 * 10^x) * 10^(-x)
  914. // the approximation of 10^(-x) was rounded up to 118 bits
  915. __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  916. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  917. Cstar.w[1] = P256.w[3];
  918. Cstar.w[0] = P256.w[2];
  919. fstar.w[3] = 0;
  920. fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  921. fstar.w[1] = P256.w[1];
  922. fstar.w[0] = P256.w[0];
  923. } else { // 22 <= ind - 1 <= 33
  924. Cstar.w[1] = 0;
  925. Cstar.w[0] = P256.w[3];
  926. fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  927. fstar.w[2] = P256.w[2];
  928. fstar.w[1] = P256.w[1];
  929. fstar.w[0] = P256.w[0];
  930. }
  931. // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  932. // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  933. // if (0 < f* < 10^(-x)) then the result is a midpoint
  934. // if floor(C*) is even then C* = floor(C*) - logical right
  935. // shift; C* has p decimal digits, correct by Prop. 1)
  936. // else if floor(C*) is odd C* = floor(C*)-1 (logical right
  937. // shift; C* has p decimal digits, correct by Pr. 1)
  938. // else
  939. // C* = floor(C*) (logical right shift; C has p decimal digits,
  940. // correct by Property 1)
  941. // n = C* * 10^(e+x)
  942. // shift right C* by Ex-128 = shiftright128[ind]
  943. shift = shiftright128[ind - 1]; // 0 <= shift <= 102
  944. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  945. Cstar.w[0] =
  946. (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  947. // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  948. } else { // 22 <= ind - 1 <= 33
  949. Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
  950. }
  951. // determine inexactness of the rounding of C*
  952. // if (0 < f* - 1/2 < 10^(-x)) then
  953. // the result is exact
  954. // else // if (f* - 1/2 > T*) then
  955. // the result is inexact
  956. if (ind - 1 <= 2) {
  957. if (fstar.w[1] > 0x8000000000000000ull ||
  958. (fstar.w[1] == 0x8000000000000000ull
  959. && fstar.w[0] > 0x0ull)) {
  960. // f* > 1/2 and the result may be exact
  961. tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
  962. if (tmp64 > ten2mk128trunc[ind - 1].w[1]
  963. || (tmp64 == ten2mk128trunc[ind - 1].w[1]
  964. && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
  965. } // else the result is exact
  966. } else { // the result is inexact; f2* <= 1/2
  967. is_inexact_gt_midpoint = 1;
  968. }
  969. } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
  970. if (fstar.w[3] > 0x0 ||
  971. (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
  972. (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
  973. (fstar.w[1] || fstar.w[0]))) {
  974. // f2* > 1/2 and the result may be exact
  975. // Calculate f2* - 1/2
  976. tmp64 = fstar.w[2] - onehalf128[ind - 1];
  977. tmp64A = fstar.w[3];
  978. if (tmp64 > fstar.w[2])
  979. tmp64A--;
  980. if (tmp64A || tmp64
  981. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  982. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  983. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  984. } // else the result is exact
  985. } else { // the result is inexact; f2* <= 1/2
  986. is_inexact_gt_midpoint = 1;
  987. }
  988. } else { // if 22 <= ind <= 33
  989. if (fstar.w[3] > onehalf128[ind - 1] ||
  990. (fstar.w[3] == onehalf128[ind - 1] &&
  991. (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
  992. // f2* > 1/2 and the result may be exact
  993. // Calculate f2* - 1/2
  994. tmp64 = fstar.w[3] - onehalf128[ind - 1];
  995. if (tmp64 || fstar.w[2]
  996. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  997. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  998. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  999. } // else the result is exact
  1000. } else { // the result is inexact; f2* <= 1/2
  1001. is_inexact_gt_midpoint = 1;
  1002. }
  1003. }
  1004. // if the result was a midpoint it was rounded away from zero, so
  1005. // it will need a correction
  1006. // check for midpoints
  1007. if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
  1008. && (fstar.w[1] || fstar.w[0])
  1009. && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
  1010. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1011. && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
  1012. // the result is a midpoint; round to nearest
  1013. if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
  1014. // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
  1015. Cstar.w[0]--; // Cstar.w[0] is now even
  1016. is_inexact_gt_midpoint = 0;
  1017. } else { // else MP in [ODD, EVEN]
  1018. is_midpoint_lt_even = 1;
  1019. is_inexact_gt_midpoint = 0;
  1020. }
  1021. }
  1022. // general correction for RM
  1023. if (is_midpoint_lt_even || is_inexact_gt_midpoint) {
  1024. Cstar.w[0] = Cstar.w[0] - 1;
  1025. } else {
  1026. ; // the result is already correct
  1027. }
  1028. res = Cstar.w[0];
  1029. } else if (exp == 0) {
  1030. // 1 <= q <= 10
  1031. // res = +C (exact)
  1032. res = C1.w[0];
  1033. } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
  1034. // res = +C * 10^exp (exact)
  1035. res = C1.w[0] * ten2k64[exp];
  1036. }
  1037. }
  1038. }
  1039. BID_RETURN (res);
  1040. }
  1041. /*****************************************************************************
  1042. * BID128_to_uint32_xfloor
  1043. ****************************************************************************/
  1044. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
  1045. bid128_to_uint32_xfloor, x)
  1046. unsigned int res;
  1047. UINT64 x_sign;
  1048. UINT64 x_exp;
  1049. int exp; // unbiased exponent
  1050. // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  1051. UINT64 tmp64, tmp64A;
  1052. BID_UI64DOUBLE tmp1;
  1053. unsigned int x_nr_bits;
  1054. int q, ind, shift;
  1055. UINT128 C1, C;
  1056. UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
  1057. UINT256 fstar;
  1058. UINT256 P256;
  1059. int is_inexact_gt_midpoint = 0;
  1060. int is_midpoint_lt_even = 0;
  1061. // unpack x
  1062. x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
  1063. x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
  1064. C1.w[1] = x.w[1] & MASK_COEFF;
  1065. C1.w[0] = x.w[0];
  1066. // check for NaN or Infinity
  1067. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  1068. // x is special
  1069. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  1070. if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
  1071. // set invalid flag
  1072. *pfpsf |= INVALID_EXCEPTION;
  1073. // return Integer Indefinite
  1074. res = 0x80000000;
  1075. } else { // x is QNaN
  1076. // set invalid flag
  1077. *pfpsf |= INVALID_EXCEPTION;
  1078. // return Integer Indefinite
  1079. res = 0x80000000;
  1080. }
  1081. BID_RETURN (res);
  1082. } else { // x is not a NaN, so it must be infinity
  1083. if (!x_sign) { // x is +inf
  1084. // set invalid flag
  1085. *pfpsf |= INVALID_EXCEPTION;
  1086. // return Integer Indefinite
  1087. res = 0x80000000;
  1088. } else { // x is -inf
  1089. // set invalid flag
  1090. *pfpsf |= INVALID_EXCEPTION;
  1091. // return Integer Indefinite
  1092. res = 0x80000000;
  1093. }
  1094. BID_RETURN (res);
  1095. }
  1096. }
  1097. // check for non-canonical values (after the check for special values)
  1098. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  1099. || (C1.w[1] == 0x0001ed09bead87c0ull
  1100. && (C1.w[0] > 0x378d8e63ffffffffull))
  1101. || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  1102. res = 0x00000000;
  1103. BID_RETURN (res);
  1104. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  1105. // x is 0
  1106. res = 0x00000000;
  1107. BID_RETURN (res);
  1108. } else { // x is not special and is not zero
  1109. // x < 0 is invalid
  1110. if (x_sign) {
  1111. // set invalid flag
  1112. *pfpsf |= INVALID_EXCEPTION;
  1113. // return Integer Indefinite
  1114. res = 0x80000000;
  1115. BID_RETURN (res);
  1116. }
  1117. // x > 0 from this point on
  1118. // q = nr. of decimal digits in x
  1119. // determine first the nr. of bits in x
  1120. if (C1.w[1] == 0) {
  1121. if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
  1122. // split the 64-bit value in two 32-bit halves to avoid rounding errors
  1123. if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
  1124. tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
  1125. x_nr_bits =
  1126. 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1127. } else { // x < 2^32
  1128. tmp1.d = (double) (C1.w[0]); // exact conversion
  1129. x_nr_bits =
  1130. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1131. }
  1132. } else { // if x < 2^53
  1133. tmp1.d = (double) C1.w[0]; // exact conversion
  1134. x_nr_bits =
  1135. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1136. }
  1137. } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  1138. tmp1.d = (double) C1.w[1]; // exact conversion
  1139. x_nr_bits =
  1140. 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1141. }
  1142. q = nr_digits[x_nr_bits - 1].digits;
  1143. if (q == 0) {
  1144. q = nr_digits[x_nr_bits - 1].digits1;
  1145. if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  1146. || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  1147. && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  1148. q++;
  1149. }
  1150. exp = (x_exp >> 49) - 6176;
  1151. if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
  1152. // set invalid flag
  1153. *pfpsf |= INVALID_EXCEPTION;
  1154. // return Integer Indefinite
  1155. res = 0x80000000;
  1156. BID_RETURN (res);
  1157. } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
  1158. // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
  1159. // so x rounded to an integer may or may not fit in a signed 32-bit int
  1160. // the cases that do not fit are identified here; the ones that fit
  1161. // fall through and will be handled with other cases further,
  1162. // under '1 <= q + exp <= 10'
  1163. // n > 0 and q + exp = 10
  1164. // if n >= 2^32 then n is too large
  1165. // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
  1166. // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
  1167. if (q <= 11) {
  1168. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  1169. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  1170. if (tmp64 >= 0xa00000000ull) {
  1171. // set invalid flag
  1172. *pfpsf |= INVALID_EXCEPTION;
  1173. // return Integer Indefinite
  1174. res = 0x80000000;
  1175. BID_RETURN (res);
  1176. }
  1177. // else cases that can be rounded to a 32-bit int fall through
  1178. // to '1 <= q + exp <= 10'
  1179. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  1180. // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
  1181. // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
  1182. // (scale 2^32 up)
  1183. tmp64 = 0xa00000000ull;
  1184. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  1185. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  1186. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  1187. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  1188. }
  1189. if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  1190. // set invalid flag
  1191. *pfpsf |= INVALID_EXCEPTION;
  1192. // return Integer Indefinite
  1193. res = 0x80000000;
  1194. BID_RETURN (res);
  1195. }
  1196. // else cases that can be rounded to a 32-bit int fall through
  1197. // to '1 <= q + exp <= 10'
  1198. }
  1199. }
  1200. // n is not too large to be converted to int32: 0 <= n < 2^31
  1201. // Note: some of the cases tested for above fall through to this point
  1202. if ((q + exp) <= 0) {
  1203. // n = +0.0...c(0)c(1)...c(q-1) or n = +0.c(0)c(1)...c(q-1)
  1204. // set inexact flag
  1205. *pfpsf |= INEXACT_EXCEPTION;
  1206. // return 0
  1207. res = 0x00000000;
  1208. BID_RETURN (res);
  1209. } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
  1210. // 1 <= x < 2^32 so x can be rounded down to a 32-bit unsigned integer
  1211. if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
  1212. ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
  1213. // chop off ind digits from the lower part of C1
  1214. // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
  1215. tmp64 = C1.w[0];
  1216. if (ind <= 19) {
  1217. C1.w[0] = C1.w[0] + midpoint64[ind - 1];
  1218. } else {
  1219. C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
  1220. C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
  1221. }
  1222. if (C1.w[0] < tmp64)
  1223. C1.w[1]++;
  1224. // calculate C* and f*
  1225. // C* is actually floor(C*) in this case
  1226. // C* and f* need shifting and masking, as shown by
  1227. // shiftright128[] and maskhigh128[]
  1228. // 1 <= x <= 33
  1229. // kx = 10^(-x) = ten2mk128[ind - 1]
  1230. // C* = (C1 + 1/2 * 10^x) * 10^(-x)
  1231. // the approximation of 10^(-x) was rounded up to 118 bits
  1232. __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  1233. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  1234. Cstar.w[1] = P256.w[3];
  1235. Cstar.w[0] = P256.w[2];
  1236. fstar.w[3] = 0;
  1237. fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  1238. fstar.w[1] = P256.w[1];
  1239. fstar.w[0] = P256.w[0];
  1240. } else { // 22 <= ind - 1 <= 33
  1241. Cstar.w[1] = 0;
  1242. Cstar.w[0] = P256.w[3];
  1243. fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  1244. fstar.w[2] = P256.w[2];
  1245. fstar.w[1] = P256.w[1];
  1246. fstar.w[0] = P256.w[0];
  1247. }
  1248. // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  1249. // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  1250. // if (0 < f* < 10^(-x)) then the result is a midpoint
  1251. // if floor(C*) is even then C* = floor(C*) - logical right
  1252. // shift; C* has p decimal digits, correct by Prop. 1)
  1253. // else if floor(C*) is odd C* = floor(C*)-1 (logical right
  1254. // shift; C* has p decimal digits, correct by Pr. 1)
  1255. // else
  1256. // C* = floor(C*) (logical right shift; C has p decimal digits,
  1257. // correct by Property 1)
  1258. // n = C* * 10^(e+x)
  1259. // shift right C* by Ex-128 = shiftright128[ind]
  1260. shift = shiftright128[ind - 1]; // 0 <= shift <= 102
  1261. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  1262. Cstar.w[0] =
  1263. (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  1264. // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  1265. } else { // 22 <= ind - 1 <= 33
  1266. Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
  1267. }
  1268. // determine inexactness of the rounding of C*
  1269. // if (0 < f* - 1/2 < 10^(-x)) then
  1270. // the result is exact
  1271. // else // if (f* - 1/2 > T*) then
  1272. // the result is inexact
  1273. if (ind - 1 <= 2) {
  1274. if (fstar.w[1] > 0x8000000000000000ull ||
  1275. (fstar.w[1] == 0x8000000000000000ull
  1276. && fstar.w[0] > 0x0ull)) {
  1277. // f* > 1/2 and the result may be exact
  1278. tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
  1279. if (tmp64 > ten2mk128trunc[ind - 1].w[1]
  1280. || (tmp64 == ten2mk128trunc[ind - 1].w[1]
  1281. && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
  1282. // set the inexact flag
  1283. *pfpsf |= INEXACT_EXCEPTION;
  1284. } // else the result is exact
  1285. } else { // the result is inexact; f2* <= 1/2
  1286. // set the inexact flag
  1287. *pfpsf |= INEXACT_EXCEPTION;
  1288. is_inexact_gt_midpoint = 1;
  1289. }
  1290. } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
  1291. if (fstar.w[3] > 0x0 ||
  1292. (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
  1293. (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
  1294. (fstar.w[1] || fstar.w[0]))) {
  1295. // f2* > 1/2 and the result may be exact
  1296. // Calculate f2* - 1/2
  1297. tmp64 = fstar.w[2] - onehalf128[ind - 1];
  1298. tmp64A = fstar.w[3];
  1299. if (tmp64 > fstar.w[2])
  1300. tmp64A--;
  1301. if (tmp64A || tmp64
  1302. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  1303. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1304. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  1305. // set the inexact flag
  1306. *pfpsf |= INEXACT_EXCEPTION;
  1307. } // else the result is exact
  1308. } else { // the result is inexact; f2* <= 1/2
  1309. // set the inexact flag
  1310. *pfpsf |= INEXACT_EXCEPTION;
  1311. is_inexact_gt_midpoint = 1;
  1312. }
  1313. } else { // if 22 <= ind <= 33
  1314. if (fstar.w[3] > onehalf128[ind - 1] ||
  1315. (fstar.w[3] == onehalf128[ind - 1] &&
  1316. (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
  1317. // f2* > 1/2 and the result may be exact
  1318. // Calculate f2* - 1/2
  1319. tmp64 = fstar.w[3] - onehalf128[ind - 1];
  1320. if (tmp64 || fstar.w[2]
  1321. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  1322. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1323. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  1324. // set the inexact flag
  1325. *pfpsf |= INEXACT_EXCEPTION;
  1326. } // else the result is exact
  1327. } else { // the result is inexact; f2* <= 1/2
  1328. // set the inexact flag
  1329. *pfpsf |= INEXACT_EXCEPTION;
  1330. is_inexact_gt_midpoint = 1;
  1331. }
  1332. }
  1333. // if the result was a midpoint it was rounded away from zero, so
  1334. // it will need a correction
  1335. // check for midpoints
  1336. if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
  1337. && (fstar.w[1] || fstar.w[0])
  1338. && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
  1339. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1340. && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
  1341. // the result is a midpoint; round to nearest
  1342. if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
  1343. // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
  1344. Cstar.w[0]--; // Cstar.w[0] is now even
  1345. is_inexact_gt_midpoint = 0;
  1346. } else { // else MP in [ODD, EVEN]
  1347. is_midpoint_lt_even = 1;
  1348. is_inexact_gt_midpoint = 0;
  1349. }
  1350. }
  1351. // general correction for RM
  1352. if (is_midpoint_lt_even || is_inexact_gt_midpoint) {
  1353. Cstar.w[0] = Cstar.w[0] - 1;
  1354. } else {
  1355. ; // the result is already correct
  1356. }
  1357. res = Cstar.w[0];
  1358. } else if (exp == 0) {
  1359. // 1 <= q <= 10
  1360. // res = +C (exact)
  1361. res = C1.w[0];
  1362. } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
  1363. // res = +C * 10^exp (exact)
  1364. res = C1.w[0] * ten2k64[exp];
  1365. }
  1366. }
  1367. }
  1368. BID_RETURN (res);
  1369. }
  1370. /*****************************************************************************
  1371. * BID128_to_uint32_ceil
  1372. ****************************************************************************/
  1373. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
  1374. bid128_to_uint32_ceil, x)
  1375. unsigned int res;
  1376. UINT64 x_sign;
  1377. UINT64 x_exp;
  1378. int exp; // unbiased exponent
  1379. // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  1380. UINT64 tmp64, tmp64A;
  1381. BID_UI64DOUBLE tmp1;
  1382. unsigned int x_nr_bits;
  1383. int q, ind, shift;
  1384. UINT128 C1, C;
  1385. UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
  1386. UINT256 fstar;
  1387. UINT256 P256;
  1388. int is_inexact_lt_midpoint = 0;
  1389. int is_midpoint_gt_even = 0;
  1390. // unpack x
  1391. x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
  1392. x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
  1393. C1.w[1] = x.w[1] & MASK_COEFF;
  1394. C1.w[0] = x.w[0];
  1395. // check for NaN or Infinity
  1396. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  1397. // x is special
  1398. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  1399. if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
  1400. // set invalid flag
  1401. *pfpsf |= INVALID_EXCEPTION;
  1402. // return Integer Indefinite
  1403. res = 0x80000000;
  1404. } else { // x is QNaN
  1405. // set invalid flag
  1406. *pfpsf |= INVALID_EXCEPTION;
  1407. // return Integer Indefinite
  1408. res = 0x80000000;
  1409. }
  1410. BID_RETURN (res);
  1411. } else { // x is not a NaN, so it must be infinity
  1412. if (!x_sign) { // x is +inf
  1413. // set invalid flag
  1414. *pfpsf |= INVALID_EXCEPTION;
  1415. // return Integer Indefinite
  1416. res = 0x80000000;
  1417. } else { // x is -inf
  1418. // set invalid flag
  1419. *pfpsf |= INVALID_EXCEPTION;
  1420. // return Integer Indefinite
  1421. res = 0x80000000;
  1422. }
  1423. BID_RETURN (res);
  1424. }
  1425. }
  1426. // check for non-canonical values (after the check for special values)
  1427. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  1428. || (C1.w[1] == 0x0001ed09bead87c0ull
  1429. && (C1.w[0] > 0x378d8e63ffffffffull))
  1430. || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  1431. res = 0x00000000;
  1432. BID_RETURN (res);
  1433. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  1434. // x is 0
  1435. res = 0x00000000;
  1436. BID_RETURN (res);
  1437. } else { // x is not special and is not zero
  1438. // q = nr. of decimal digits in x
  1439. // determine first the nr. of bits in x
  1440. if (C1.w[1] == 0) {
  1441. if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
  1442. // split the 64-bit value in two 32-bit halves to avoid rounding errors
  1443. if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
  1444. tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
  1445. x_nr_bits =
  1446. 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1447. } else { // x < 2^32
  1448. tmp1.d = (double) (C1.w[0]); // exact conversion
  1449. x_nr_bits =
  1450. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1451. }
  1452. } else { // if x < 2^53
  1453. tmp1.d = (double) C1.w[0]; // exact conversion
  1454. x_nr_bits =
  1455. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1456. }
  1457. } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  1458. tmp1.d = (double) C1.w[1]; // exact conversion
  1459. x_nr_bits =
  1460. 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1461. }
  1462. q = nr_digits[x_nr_bits - 1].digits;
  1463. if (q == 0) {
  1464. q = nr_digits[x_nr_bits - 1].digits1;
  1465. if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  1466. || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  1467. && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  1468. q++;
  1469. }
  1470. exp = (x_exp >> 49) - 6176;
  1471. if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
  1472. // set invalid flag
  1473. *pfpsf |= INVALID_EXCEPTION;
  1474. // return Integer Indefinite
  1475. res = 0x80000000;
  1476. BID_RETURN (res);
  1477. } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
  1478. // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
  1479. // so x rounded to an integer may or may not fit in a signed 32-bit int
  1480. // the cases that do not fit are identified here; the ones that fit
  1481. // fall through and will be handled with other cases further,
  1482. // under '1 <= q + exp <= 10'
  1483. if (x_sign) { // if n < 0 and q + exp = 10
  1484. // if n <= -1 then n is too large
  1485. // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
  1486. // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34
  1487. if (q <= 11) {
  1488. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  1489. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  1490. if (tmp64 >= 0x0aull) {
  1491. // set invalid flag
  1492. *pfpsf |= INVALID_EXCEPTION;
  1493. // return Integer Indefinite
  1494. res = 0x80000000;
  1495. BID_RETURN (res);
  1496. }
  1497. // else cases that can be rounded to a 32-bit int fall through
  1498. // to '1 <= q + exp <= 10'
  1499. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  1500. // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
  1501. // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
  1502. // (scale 1 up)
  1503. tmp64 = 0x0aull;
  1504. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  1505. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  1506. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  1507. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  1508. }
  1509. if (C1.w[1] > C.w[1]
  1510. || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  1511. // set invalid flag
  1512. *pfpsf |= INVALID_EXCEPTION;
  1513. // return Integer Indefinite
  1514. res = 0x80000000;
  1515. BID_RETURN (res);
  1516. }
  1517. // else cases that can be rounded to a 32-bit int fall through
  1518. // to '1 <= q + exp <= 10'
  1519. }
  1520. } else { // if n > 0 and q + exp = 10
  1521. // if n > 2^32 - 1 then n is too large
  1522. // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^32 - 1
  1523. // too large if 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6, 1<=q<=34
  1524. if (q <= 11) {
  1525. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  1526. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  1527. if (tmp64 > 0x9fffffff6ull) {
  1528. // set invalid flag
  1529. *pfpsf |= INVALID_EXCEPTION;
  1530. // return Integer Indefinite
  1531. res = 0x80000000;
  1532. BID_RETURN (res);
  1533. }
  1534. // else cases that can be rounded to a 32-bit int fall through
  1535. // to '1 <= q + exp <= 10'
  1536. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  1537. // 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6 <=>
  1538. // C > 0x9fffffff6 * 10^(q-11) where 1 <= q - 11 <= 23
  1539. // (scale 2^32 up)
  1540. tmp64 = 0x9fffffff6ull;
  1541. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  1542. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  1543. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  1544. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  1545. }
  1546. if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
  1547. // set invalid flag
  1548. *pfpsf |= INVALID_EXCEPTION;
  1549. // return Integer Indefinite
  1550. res = 0x80000000;
  1551. BID_RETURN (res);
  1552. }
  1553. // else cases that can be rounded to a 32-bit int fall through
  1554. // to '1 <= q + exp <= 10'
  1555. }
  1556. }
  1557. }
  1558. // n is not too large to be converted to int32: -2^32-1 < n <= 2^32-1
  1559. // Note: some of the cases tested for above fall through to this point
  1560. if ((q + exp) <= 0) {
  1561. // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
  1562. // return 0
  1563. if (x_sign)
  1564. res = 0x00000000;
  1565. else
  1566. res = 0x00000001;
  1567. BID_RETURN (res);
  1568. } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
  1569. // -2^32-1 < x <= -1 or 1 <= x <= 2^32-1 so x can be rounded
  1570. // toward positive infinity to a 32-bit signed integer
  1571. if (x_sign) { // x <= -1 is invalid
  1572. // set invalid flag
  1573. *pfpsf |= INVALID_EXCEPTION;
  1574. // return Integer Indefinite
  1575. res = 0x80000000;
  1576. BID_RETURN (res);
  1577. }
  1578. // x > 0 from this point on
  1579. if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
  1580. ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
  1581. // chop off ind digits from the lower part of C1
  1582. // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
  1583. tmp64 = C1.w[0];
  1584. if (ind <= 19) {
  1585. C1.w[0] = C1.w[0] + midpoint64[ind - 1];
  1586. } else {
  1587. C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
  1588. C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
  1589. }
  1590. if (C1.w[0] < tmp64)
  1591. C1.w[1]++;
  1592. // calculate C* and f*
  1593. // C* is actually floor(C*) in this case
  1594. // C* and f* need shifting and masking, as shown by
  1595. // shiftright128[] and maskhigh128[]
  1596. // 1 <= x <= 33
  1597. // kx = 10^(-x) = ten2mk128[ind - 1]
  1598. // C* = (C1 + 1/2 * 10^x) * 10^(-x)
  1599. // the approximation of 10^(-x) was rounded up to 118 bits
  1600. __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  1601. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  1602. Cstar.w[1] = P256.w[3];
  1603. Cstar.w[0] = P256.w[2];
  1604. fstar.w[3] = 0;
  1605. fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  1606. fstar.w[1] = P256.w[1];
  1607. fstar.w[0] = P256.w[0];
  1608. } else { // 22 <= ind - 1 <= 33
  1609. Cstar.w[1] = 0;
  1610. Cstar.w[0] = P256.w[3];
  1611. fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  1612. fstar.w[2] = P256.w[2];
  1613. fstar.w[1] = P256.w[1];
  1614. fstar.w[0] = P256.w[0];
  1615. }
  1616. // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  1617. // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  1618. // if (0 < f* < 10^(-x)) then the result is a midpoint
  1619. // if floor(C*) is even then C* = floor(C*) - logical right
  1620. // shift; C* has p decimal digits, correct by Prop. 1)
  1621. // else if floor(C*) is odd C* = floor(C*)-1 (logical right
  1622. // shift; C* has p decimal digits, correct by Pr. 1)
  1623. // else
  1624. // C* = floor(C*) (logical right shift; C has p decimal digits,
  1625. // correct by Property 1)
  1626. // n = C* * 10^(e+x)
  1627. // shift right C* by Ex-128 = shiftright128[ind]
  1628. shift = shiftright128[ind - 1]; // 0 <= shift <= 102
  1629. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  1630. Cstar.w[0] =
  1631. (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  1632. // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  1633. } else { // 22 <= ind - 1 <= 33
  1634. Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
  1635. }
  1636. // determine inexactness of the rounding of C*
  1637. // if (0 < f* - 1/2 < 10^(-x)) then
  1638. // the result is exact
  1639. // else // if (f* - 1/2 > T*) then
  1640. // the result is inexact
  1641. if (ind - 1 <= 2) {
  1642. if (fstar.w[1] > 0x8000000000000000ull ||
  1643. (fstar.w[1] == 0x8000000000000000ull
  1644. && fstar.w[0] > 0x0ull)) {
  1645. // f* > 1/2 and the result may be exact
  1646. tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
  1647. if (tmp64 > ten2mk128trunc[ind - 1].w[1]
  1648. || (tmp64 == ten2mk128trunc[ind - 1].w[1]
  1649. && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
  1650. is_inexact_lt_midpoint = 1;
  1651. } // else the result is exact
  1652. } else { // the result is inexact; f2* <= 1/2
  1653. ;
  1654. }
  1655. } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
  1656. if (fstar.w[3] > 0x0 ||
  1657. (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
  1658. (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
  1659. (fstar.w[1] || fstar.w[0]))) {
  1660. // f2* > 1/2 and the result may be exact
  1661. // Calculate f2* - 1/2
  1662. tmp64 = fstar.w[2] - onehalf128[ind - 1];
  1663. tmp64A = fstar.w[3];
  1664. if (tmp64 > fstar.w[2])
  1665. tmp64A--;
  1666. if (tmp64A || tmp64
  1667. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  1668. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1669. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  1670. is_inexact_lt_midpoint = 1;
  1671. } // else the result is exact
  1672. } else { // the result is inexact; f2* <= 1/2
  1673. }
  1674. } else { // if 22 <= ind <= 33
  1675. if (fstar.w[3] > onehalf128[ind - 1] ||
  1676. (fstar.w[3] == onehalf128[ind - 1] &&
  1677. (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
  1678. // f2* > 1/2 and the result may be exact
  1679. // Calculate f2* - 1/2
  1680. tmp64 = fstar.w[3] - onehalf128[ind - 1];
  1681. if (tmp64 || fstar.w[2]
  1682. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  1683. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1684. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  1685. is_inexact_lt_midpoint = 1;
  1686. } // else the result is exact
  1687. } else { // the result is inexact; f2* <= 1/2
  1688. ;
  1689. }
  1690. }
  1691. // if the result was a midpoint it was rounded away from zero, so
  1692. // it will need a correction
  1693. // check for midpoints
  1694. if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
  1695. && (fstar.w[1] || fstar.w[0])
  1696. && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
  1697. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  1698. && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
  1699. // the result is a midpoint; round to nearest
  1700. if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
  1701. // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
  1702. Cstar.w[0]--; // Cstar.w[0] is now even
  1703. is_midpoint_gt_even = 1;
  1704. is_inexact_lt_midpoint = 0;
  1705. } else { // else MP in [ODD, EVEN]
  1706. is_inexact_lt_midpoint = 0;
  1707. }
  1708. }
  1709. // general correction for RM
  1710. if (is_midpoint_gt_even || is_inexact_lt_midpoint) {
  1711. Cstar.w[0] = Cstar.w[0] + 1;
  1712. } else {
  1713. ; // the result is already correct
  1714. }
  1715. res = Cstar.w[0];
  1716. } else if (exp == 0) {
  1717. // 1 <= q <= 10
  1718. // res = +C (exact)
  1719. res = C1.w[0];
  1720. } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
  1721. // res = +C * 10^exp (exact)
  1722. res = C1.w[0] * ten2k64[exp];
  1723. }
  1724. }
  1725. }
  1726. BID_RETURN (res);
  1727. }
  1728. /*****************************************************************************
  1729. * BID128_to_uint32_xceil
  1730. ****************************************************************************/
  1731. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
  1732. bid128_to_uint32_xceil, x)
  1733. unsigned int res;
  1734. UINT64 x_sign;
  1735. UINT64 x_exp;
  1736. int exp; // unbiased exponent
  1737. // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  1738. UINT64 tmp64, tmp64A;
  1739. BID_UI64DOUBLE tmp1;
  1740. unsigned int x_nr_bits;
  1741. int q, ind, shift;
  1742. UINT128 C1, C;
  1743. UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
  1744. UINT256 fstar;
  1745. UINT256 P256;
  1746. int is_inexact_lt_midpoint = 0;
  1747. int is_midpoint_gt_even = 0;
  1748. // unpack x
  1749. x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
  1750. x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
  1751. C1.w[1] = x.w[1] & MASK_COEFF;
  1752. C1.w[0] = x.w[0];
  1753. // check for NaN or Infinity
  1754. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  1755. // x is special
  1756. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  1757. if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
  1758. // set invalid flag
  1759. *pfpsf |= INVALID_EXCEPTION;
  1760. // return Integer Indefinite
  1761. res = 0x80000000;
  1762. } else { // x is QNaN
  1763. // set invalid flag
  1764. *pfpsf |= INVALID_EXCEPTION;
  1765. // return Integer Indefinite
  1766. res = 0x80000000;
  1767. }
  1768. BID_RETURN (res);
  1769. } else { // x is not a NaN, so it must be infinity
  1770. if (!x_sign) { // x is +inf
  1771. // set invalid flag
  1772. *pfpsf |= INVALID_EXCEPTION;
  1773. // return Integer Indefinite
  1774. res = 0x80000000;
  1775. } else { // x is -inf
  1776. // set invalid flag
  1777. *pfpsf |= INVALID_EXCEPTION;
  1778. // return Integer Indefinite
  1779. res = 0x80000000;
  1780. }
  1781. BID_RETURN (res);
  1782. }
  1783. }
  1784. // check for non-canonical values (after the check for special values)
  1785. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  1786. || (C1.w[1] == 0x0001ed09bead87c0ull
  1787. && (C1.w[0] > 0x378d8e63ffffffffull))
  1788. || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  1789. res = 0x00000000;
  1790. BID_RETURN (res);
  1791. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  1792. // x is 0
  1793. res = 0x00000000;
  1794. BID_RETURN (res);
  1795. } else { // x is not special and is not zero
  1796. // q = nr. of decimal digits in x
  1797. // determine first the nr. of bits in x
  1798. if (C1.w[1] == 0) {
  1799. if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
  1800. // split the 64-bit value in two 32-bit halves to avoid rounding errors
  1801. if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
  1802. tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
  1803. x_nr_bits =
  1804. 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1805. } else { // x < 2^32
  1806. tmp1.d = (double) (C1.w[0]); // exact conversion
  1807. x_nr_bits =
  1808. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1809. }
  1810. } else { // if x < 2^53
  1811. tmp1.d = (double) C1.w[0]; // exact conversion
  1812. x_nr_bits =
  1813. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1814. }
  1815. } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  1816. tmp1.d = (double) C1.w[1]; // exact conversion
  1817. x_nr_bits =
  1818. 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  1819. }
  1820. q = nr_digits[x_nr_bits - 1].digits;
  1821. if (q == 0) {
  1822. q = nr_digits[x_nr_bits - 1].digits1;
  1823. if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  1824. || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  1825. && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  1826. q++;
  1827. }
  1828. exp = (x_exp >> 49) - 6176;
  1829. if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
  1830. // set invalid flag
  1831. *pfpsf |= INVALID_EXCEPTION;
  1832. // return Integer Indefinite
  1833. res = 0x80000000;
  1834. BID_RETURN (res);
  1835. } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
  1836. // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
  1837. // so x rounded to an integer may or may not fit in a signed 32-bit int
  1838. // the cases that do not fit are identified here; the ones that fit
  1839. // fall through and will be handled with other cases further,
  1840. // under '1 <= q + exp <= 10'
  1841. if (x_sign) { // if n < 0 and q + exp = 10
  1842. // if n <= -1 then n is too large
  1843. // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
  1844. // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34
  1845. if (q <= 11) {
  1846. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  1847. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  1848. if (tmp64 >= 0x0aull) {
  1849. // set invalid flag
  1850. *pfpsf |= INVALID_EXCEPTION;
  1851. // return Integer Indefinite
  1852. res = 0x80000000;
  1853. BID_RETURN (res);
  1854. }
  1855. // else cases that can be rounded to a 32-bit int fall through
  1856. // to '1 <= q + exp <= 10'
  1857. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  1858. // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
  1859. // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
  1860. // (scale 1 up)
  1861. tmp64 = 0x0aull;
  1862. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  1863. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  1864. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  1865. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  1866. }
  1867. if (C1.w[1] > C.w[1]
  1868. || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  1869. // set invalid flag
  1870. *pfpsf |= INVALID_EXCEPTION;
  1871. // return Integer Indefinite
  1872. res = 0x80000000;
  1873. BID_RETURN (res);
  1874. }
  1875. // else cases that can be rounded to a 32-bit int fall through
  1876. // to '1 <= q + exp <= 10'
  1877. }
  1878. } else { // if n > 0 and q + exp = 10
  1879. // if n > 2^32 - 1 then n is too large
  1880. // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^32 - 1
  1881. // too large if 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6, 1<=q<=34
  1882. if (q <= 11) {
  1883. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  1884. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  1885. if (tmp64 > 0x9fffffff6ull) {
  1886. // set invalid flag
  1887. *pfpsf |= INVALID_EXCEPTION;
  1888. // return Integer Indefinite
  1889. res = 0x80000000;
  1890. BID_RETURN (res);
  1891. }
  1892. // else cases that can be rounded to a 32-bit int fall through
  1893. // to '1 <= q + exp <= 10'
  1894. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  1895. // 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6 <=>
  1896. // C > 0x9fffffff6 * 10^(q-11) where 1 <= q - 11 <= 23
  1897. // (scale 2^32 up)
  1898. tmp64 = 0x9fffffff6ull;
  1899. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  1900. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  1901. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  1902. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  1903. }
  1904. if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
  1905. // set invalid flag
  1906. *pfpsf |= INVALID_EXCEPTION;
  1907. // return Integer Indefinite
  1908. res = 0x80000000;
  1909. BID_RETURN (res);
  1910. }
  1911. // else cases that can be rounded to a 32-bit int fall through
  1912. // to '1 <= q + exp <= 10'
  1913. }
  1914. }
  1915. }
  1916. // n is not too large to be converted to int32: -2^32-1 < n <= 2^32-1
  1917. // Note: some of the cases tested for above fall through to this point
  1918. if ((q + exp) <= 0) {
  1919. // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
  1920. // set inexact flag
  1921. *pfpsf |= INEXACT_EXCEPTION;
  1922. // return 0
  1923. if (x_sign)
  1924. res = 0x00000000;
  1925. else
  1926. res = 0x00000001;
  1927. BID_RETURN (res);
  1928. } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
  1929. // -2^32-1 < x <= -1 or 1 <= x <= 2^32-1 so x can be rounded
  1930. // toward positive infinity to a 32-bit signed integer
  1931. if (x_sign) { // x <= -1 is invalid
  1932. // set invalid flag
  1933. *pfpsf |= INVALID_EXCEPTION;
  1934. // return Integer Indefinite
  1935. res = 0x80000000;
  1936. BID_RETURN (res);
  1937. }
  1938. // x > 0 from this point on
  1939. if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
  1940. ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
  1941. // chop off ind digits from the lower part of C1
  1942. // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
  1943. tmp64 = C1.w[0];
  1944. if (ind <= 19) {
  1945. C1.w[0] = C1.w[0] + midpoint64[ind - 1];
  1946. } else {
  1947. C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
  1948. C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
  1949. }
  1950. if (C1.w[0] < tmp64)
  1951. C1.w[1]++;
  1952. // calculate C* and f*
  1953. // C* is actually floor(C*) in this case
  1954. // C* and f* need shifting and masking, as shown by
  1955. // shiftright128[] and maskhigh128[]
  1956. // 1 <= x <= 33
  1957. // kx = 10^(-x) = ten2mk128[ind - 1]
  1958. // C* = (C1 + 1/2 * 10^x) * 10^(-x)
  1959. // the approximation of 10^(-x) was rounded up to 118 bits
  1960. __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  1961. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  1962. Cstar.w[1] = P256.w[3];
  1963. Cstar.w[0] = P256.w[2];
  1964. fstar.w[3] = 0;
  1965. fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  1966. fstar.w[1] = P256.w[1];
  1967. fstar.w[0] = P256.w[0];
  1968. } else { // 22 <= ind - 1 <= 33
  1969. Cstar.w[1] = 0;
  1970. Cstar.w[0] = P256.w[3];
  1971. fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  1972. fstar.w[2] = P256.w[2];
  1973. fstar.w[1] = P256.w[1];
  1974. fstar.w[0] = P256.w[0];
  1975. }
  1976. // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  1977. // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  1978. // if (0 < f* < 10^(-x)) then the result is a midpoint
  1979. // if floor(C*) is even then C* = floor(C*) - logical right
  1980. // shift; C* has p decimal digits, correct by Prop. 1)
  1981. // else if floor(C*) is odd C* = floor(C*)-1 (logical right
  1982. // shift; C* has p decimal digits, correct by Pr. 1)
  1983. // else
  1984. // C* = floor(C*) (logical right shift; C has p decimal digits,
  1985. // correct by Property 1)
  1986. // n = C* * 10^(e+x)
  1987. // shift right C* by Ex-128 = shiftright128[ind]
  1988. shift = shiftright128[ind - 1]; // 0 <= shift <= 102
  1989. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  1990. Cstar.w[0] =
  1991. (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  1992. // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  1993. } else { // 22 <= ind - 1 <= 33
  1994. Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
  1995. }
  1996. // determine inexactness of the rounding of C*
  1997. // if (0 < f* - 1/2 < 10^(-x)) then
  1998. // the result is exact
  1999. // else // if (f* - 1/2 > T*) then
  2000. // the result is inexact
  2001. if (ind - 1 <= 2) {
  2002. if (fstar.w[1] > 0x8000000000000000ull ||
  2003. (fstar.w[1] == 0x8000000000000000ull
  2004. && fstar.w[0] > 0x0ull)) {
  2005. // f* > 1/2 and the result may be exact
  2006. tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
  2007. if (tmp64 > ten2mk128trunc[ind - 1].w[1]
  2008. || (tmp64 == ten2mk128trunc[ind - 1].w[1]
  2009. && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
  2010. // set the inexact flag
  2011. *pfpsf |= INEXACT_EXCEPTION;
  2012. is_inexact_lt_midpoint = 1;
  2013. } // else the result is exact
  2014. } else { // the result is inexact; f2* <= 1/2
  2015. // set the inexact flag
  2016. *pfpsf |= INEXACT_EXCEPTION;
  2017. }
  2018. } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
  2019. if (fstar.w[3] > 0x0 ||
  2020. (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
  2021. (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
  2022. (fstar.w[1] || fstar.w[0]))) {
  2023. // f2* > 1/2 and the result may be exact
  2024. // Calculate f2* - 1/2
  2025. tmp64 = fstar.w[2] - onehalf128[ind - 1];
  2026. tmp64A = fstar.w[3];
  2027. if (tmp64 > fstar.w[2])
  2028. tmp64A--;
  2029. if (tmp64A || tmp64
  2030. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  2031. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  2032. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  2033. // set the inexact flag
  2034. *pfpsf |= INEXACT_EXCEPTION;
  2035. is_inexact_lt_midpoint = 1;
  2036. } // else the result is exact
  2037. } else { // the result is inexact; f2* <= 1/2
  2038. // set the inexact flag
  2039. *pfpsf |= INEXACT_EXCEPTION;
  2040. }
  2041. } else { // if 22 <= ind <= 33
  2042. if (fstar.w[3] > onehalf128[ind - 1] ||
  2043. (fstar.w[3] == onehalf128[ind - 1] &&
  2044. (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
  2045. // f2* > 1/2 and the result may be exact
  2046. // Calculate f2* - 1/2
  2047. tmp64 = fstar.w[3] - onehalf128[ind - 1];
  2048. if (tmp64 || fstar.w[2]
  2049. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  2050. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  2051. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  2052. // set the inexact flag
  2053. *pfpsf |= INEXACT_EXCEPTION;
  2054. is_inexact_lt_midpoint = 1;
  2055. } // else the result is exact
  2056. } else { // the result is inexact; f2* <= 1/2
  2057. // set the inexact flag
  2058. *pfpsf |= INEXACT_EXCEPTION;
  2059. }
  2060. }
  2061. // if the result was a midpoint it was rounded away from zero, so
  2062. // it will need a correction
  2063. // check for midpoints
  2064. if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
  2065. && (fstar.w[1] || fstar.w[0])
  2066. && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
  2067. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  2068. && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
  2069. // the result is a midpoint; round to nearest
  2070. if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
  2071. // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
  2072. Cstar.w[0]--; // Cstar.w[0] is now even
  2073. is_midpoint_gt_even = 1;
  2074. is_inexact_lt_midpoint = 0;
  2075. } else { // else MP in [ODD, EVEN]
  2076. is_inexact_lt_midpoint = 0;
  2077. }
  2078. }
  2079. // general correction for RM
  2080. if (is_midpoint_gt_even || is_inexact_lt_midpoint) {
  2081. Cstar.w[0] = Cstar.w[0] + 1;
  2082. } else {
  2083. ; // the result is already correct
  2084. }
  2085. res = Cstar.w[0];
  2086. } else if (exp == 0) {
  2087. // 1 <= q <= 10
  2088. // res = +C (exact)
  2089. res = C1.w[0];
  2090. } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
  2091. // res = +C * 10^exp (exact)
  2092. res = C1.w[0] * ten2k64[exp];
  2093. }
  2094. }
  2095. }
  2096. BID_RETURN (res);
  2097. }
  2098. /*****************************************************************************
  2099. * BID128_to_uint32_int
  2100. ****************************************************************************/
  2101. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
  2102. bid128_to_uint32_int, x)
  2103. int res;
  2104. UINT64 x_sign;
  2105. UINT64 x_exp;
  2106. int exp; // unbiased exponent
  2107. // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  2108. UINT64 tmp64, tmp64A;
  2109. BID_UI64DOUBLE tmp1;
  2110. unsigned int x_nr_bits;
  2111. int q, ind, shift;
  2112. UINT128 C1, C;
  2113. UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
  2114. UINT256 fstar;
  2115. UINT256 P256;
  2116. int is_inexact_gt_midpoint = 0;
  2117. int is_midpoint_lt_even = 0;
  2118. // unpack x
  2119. x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
  2120. x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
  2121. C1.w[1] = x.w[1] & MASK_COEFF;
  2122. C1.w[0] = x.w[0];
  2123. // check for NaN or Infinity
  2124. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  2125. // x is special
  2126. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  2127. if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
  2128. // set invalid flag
  2129. *pfpsf |= INVALID_EXCEPTION;
  2130. // return Integer Indefinite
  2131. res = 0x80000000;
  2132. } else { // x is QNaN
  2133. // set invalid flag
  2134. *pfpsf |= INVALID_EXCEPTION;
  2135. // return Integer Indefinite
  2136. res = 0x80000000;
  2137. }
  2138. BID_RETURN (res);
  2139. } else { // x is not a NaN, so it must be infinity
  2140. if (!x_sign) { // x is +inf
  2141. // set invalid flag
  2142. *pfpsf |= INVALID_EXCEPTION;
  2143. // return Integer Indefinite
  2144. res = 0x80000000;
  2145. } else { // x is -inf
  2146. // set invalid flag
  2147. *pfpsf |= INVALID_EXCEPTION;
  2148. // return Integer Indefinite
  2149. res = 0x80000000;
  2150. }
  2151. BID_RETURN (res);
  2152. }
  2153. }
  2154. // check for non-canonical values (after the check for special values)
  2155. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  2156. || (C1.w[1] == 0x0001ed09bead87c0ull
  2157. && (C1.w[0] > 0x378d8e63ffffffffull))
  2158. || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  2159. res = 0x00000000;
  2160. BID_RETURN (res);
  2161. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  2162. // x is 0
  2163. res = 0x00000000;
  2164. BID_RETURN (res);
  2165. } else { // x is not special and is not zero
  2166. // q = nr. of decimal digits in x
  2167. // determine first the nr. of bits in x
  2168. if (C1.w[1] == 0) {
  2169. if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
  2170. // split the 64-bit value in two 32-bit halves to avoid rounding errors
  2171. if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
  2172. tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
  2173. x_nr_bits =
  2174. 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2175. } else { // x < 2^32
  2176. tmp1.d = (double) (C1.w[0]); // exact conversion
  2177. x_nr_bits =
  2178. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2179. }
  2180. } else { // if x < 2^53
  2181. tmp1.d = (double) C1.w[0]; // exact conversion
  2182. x_nr_bits =
  2183. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2184. }
  2185. } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  2186. tmp1.d = (double) C1.w[1]; // exact conversion
  2187. x_nr_bits =
  2188. 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2189. }
  2190. q = nr_digits[x_nr_bits - 1].digits;
  2191. if (q == 0) {
  2192. q = nr_digits[x_nr_bits - 1].digits1;
  2193. if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  2194. || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  2195. && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  2196. q++;
  2197. }
  2198. exp = (x_exp >> 49) - 6176;
  2199. if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
  2200. // set invalid flag
  2201. *pfpsf |= INVALID_EXCEPTION;
  2202. // return Integer Indefinite
  2203. res = 0x80000000;
  2204. BID_RETURN (res);
  2205. } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
  2206. // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
  2207. // so x rounded to an integer may or may not fit in a signed 32-bit int
  2208. // the cases that do not fit are identified here; the ones that fit
  2209. // fall through and will be handled with other cases further,
  2210. // under '1 <= q + exp <= 10'
  2211. if (x_sign) { // if n < 0 and q + exp = 10
  2212. // if n <= -1 then n is too large
  2213. // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
  2214. // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a, 1<=q<=34
  2215. if (q <= 11) {
  2216. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  2217. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  2218. if (tmp64 >= 0x0aull) {
  2219. // set invalid flag
  2220. *pfpsf |= INVALID_EXCEPTION;
  2221. // return Integer Indefinite
  2222. res = 0x80000000;
  2223. BID_RETURN (res);
  2224. }
  2225. // else cases that can be rounded to a 32-bit uint fall through
  2226. // to '1 <= q + exp <= 10'
  2227. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  2228. // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
  2229. // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
  2230. // (scale 1 up)
  2231. tmp64 = 0x0aull;
  2232. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  2233. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  2234. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  2235. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  2236. }
  2237. if (C1.w[1] > C.w[1]
  2238. || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  2239. // set invalid flag
  2240. *pfpsf |= INVALID_EXCEPTION;
  2241. // return Integer Indefinite
  2242. res = 0x80000000;
  2243. BID_RETURN (res);
  2244. }
  2245. // else cases that can be rounded to a 32-bit int fall through
  2246. // to '1 <= q + exp <= 10'
  2247. }
  2248. } else { // if n > 0 and q + exp = 10
  2249. // if n >= 2^32 then n is too large
  2250. // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
  2251. // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
  2252. if (q <= 11) {
  2253. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  2254. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  2255. if (tmp64 >= 0xa00000000ull) {
  2256. // set invalid flag
  2257. *pfpsf |= INVALID_EXCEPTION;
  2258. // return Integer Indefinite
  2259. res = 0x80000000;
  2260. BID_RETURN (res);
  2261. }
  2262. // else cases that can be rounded to a 32-bit uint fall through
  2263. // to '1 <= q + exp <= 10'
  2264. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  2265. // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
  2266. // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
  2267. // (scale 2^32 up)
  2268. tmp64 = 0xa00000000ull;
  2269. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  2270. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  2271. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  2272. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  2273. }
  2274. if (C1.w[1] > C.w[1]
  2275. || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  2276. // set invalid flag
  2277. *pfpsf |= INVALID_EXCEPTION;
  2278. // return Integer Indefinite
  2279. res = 0x80000000;
  2280. BID_RETURN (res);
  2281. }
  2282. // else cases that can be rounded to a 32-bit int fall through
  2283. // to '1 <= q + exp <= 10'
  2284. }
  2285. }
  2286. }
  2287. // n is not too large to be converted to uint32: -2^32 < n < 2^32
  2288. // Note: some of the cases tested for above fall through to this point
  2289. if ((q + exp) <= 0) {
  2290. // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
  2291. // return 0
  2292. res = 0x00000000;
  2293. BID_RETURN (res);
  2294. } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
  2295. // x = d(0)...d(k).d(k+1)..., k >= 0, d(0) != 0
  2296. if (x_sign) { // x <= -1
  2297. // set invalid flag
  2298. *pfpsf |= INVALID_EXCEPTION;
  2299. // return Integer Indefinite
  2300. res = 0x80000000;
  2301. BID_RETURN (res);
  2302. }
  2303. // x > 0 from this point on
  2304. // 1 <= x < 2^32 so x can be rounded to zero to a 32-bit unsigned integer
  2305. if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
  2306. ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
  2307. // chop off ind digits from the lower part of C1
  2308. // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
  2309. tmp64 = C1.w[0];
  2310. if (ind <= 19) {
  2311. C1.w[0] = C1.w[0] + midpoint64[ind - 1];
  2312. } else {
  2313. C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
  2314. C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
  2315. }
  2316. if (C1.w[0] < tmp64)
  2317. C1.w[1]++;
  2318. // calculate C* and f*
  2319. // C* is actually floor(C*) in this case
  2320. // C* and f* need shifting and masking, as shown by
  2321. // shiftright128[] and maskhigh128[]
  2322. // 1 <= x <= 33
  2323. // kx = 10^(-x) = ten2mk128[ind - 1]
  2324. // C* = (C1 + 1/2 * 10^x) * 10^(-x)
  2325. // the approximation of 10^(-x) was rounded up to 118 bits
  2326. __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  2327. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  2328. Cstar.w[1] = P256.w[3];
  2329. Cstar.w[0] = P256.w[2];
  2330. fstar.w[3] = 0;
  2331. fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  2332. fstar.w[1] = P256.w[1];
  2333. fstar.w[0] = P256.w[0];
  2334. } else { // 22 <= ind - 1 <= 33
  2335. Cstar.w[1] = 0;
  2336. Cstar.w[0] = P256.w[3];
  2337. fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  2338. fstar.w[2] = P256.w[2];
  2339. fstar.w[1] = P256.w[1];
  2340. fstar.w[0] = P256.w[0];
  2341. }
  2342. // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  2343. // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  2344. // if (0 < f* < 10^(-x)) then the result is a midpoint
  2345. // if floor(C*) is even then C* = floor(C*) - logical right
  2346. // shift; C* has p decimal digits, correct by Prop. 1)
  2347. // else if floor(C*) is odd C* = floor(C*)-1 (logical right
  2348. // shift; C* has p decimal digits, correct by Pr. 1)
  2349. // else
  2350. // C* = floor(C*) (logical right shift; C has p decimal digits,
  2351. // correct by Property 1)
  2352. // n = C* * 10^(e+x)
  2353. // shift right C* by Ex-128 = shiftright128[ind]
  2354. shift = shiftright128[ind - 1]; // 0 <= shift <= 102
  2355. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  2356. Cstar.w[0] =
  2357. (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  2358. // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  2359. } else { // 22 <= ind - 1 <= 33
  2360. Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
  2361. }
  2362. // determine inexactness of the rounding of C*
  2363. // if (0 < f* - 1/2 < 10^(-x)) then
  2364. // the result is exact
  2365. // else // if (f* - 1/2 > T*) then
  2366. // the result is inexact
  2367. if (ind - 1 <= 2) {
  2368. if (fstar.w[1] > 0x8000000000000000ull ||
  2369. (fstar.w[1] == 0x8000000000000000ull
  2370. && fstar.w[0] > 0x0ull)) {
  2371. // f* > 1/2 and the result may be exact
  2372. tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
  2373. if (tmp64 > ten2mk128trunc[ind - 1].w[1]
  2374. || (tmp64 == ten2mk128trunc[ind - 1].w[1]
  2375. && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
  2376. } // else the result is exact
  2377. } else { // the result is inexact; f2* <= 1/2
  2378. is_inexact_gt_midpoint = 1;
  2379. }
  2380. } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
  2381. if (fstar.w[3] > 0x0 ||
  2382. (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
  2383. (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
  2384. (fstar.w[1] || fstar.w[0]))) {
  2385. // f2* > 1/2 and the result may be exact
  2386. // Calculate f2* - 1/2
  2387. tmp64 = fstar.w[2] - onehalf128[ind - 1];
  2388. tmp64A = fstar.w[3];
  2389. if (tmp64 > fstar.w[2])
  2390. tmp64A--;
  2391. if (tmp64A || tmp64
  2392. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  2393. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  2394. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  2395. } // else the result is exact
  2396. } else { // the result is inexact; f2* <= 1/2
  2397. is_inexact_gt_midpoint = 1;
  2398. }
  2399. } else { // if 22 <= ind <= 33
  2400. if (fstar.w[3] > onehalf128[ind - 1] ||
  2401. (fstar.w[3] == onehalf128[ind - 1] &&
  2402. (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
  2403. // f2* > 1/2 and the result may be exact
  2404. // Calculate f2* - 1/2
  2405. tmp64 = fstar.w[3] - onehalf128[ind - 1];
  2406. if (tmp64 || fstar.w[2]
  2407. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  2408. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  2409. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  2410. } // else the result is exact
  2411. } else { // the result is inexact; f2* <= 1/2
  2412. is_inexact_gt_midpoint = 1;
  2413. }
  2414. }
  2415. // if the result was a midpoint it was rounded away from zero, so
  2416. // it will need a correction
  2417. // check for midpoints
  2418. if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
  2419. && (fstar.w[1] || fstar.w[0])
  2420. && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
  2421. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  2422. && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
  2423. // the result is a midpoint; round to nearest
  2424. if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
  2425. // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
  2426. Cstar.w[0]--; // Cstar.w[0] is now even
  2427. is_inexact_gt_midpoint = 0;
  2428. } else { // else MP in [ODD, EVEN]
  2429. is_midpoint_lt_even = 1;
  2430. is_inexact_gt_midpoint = 0;
  2431. }
  2432. }
  2433. // general correction for RZ
  2434. if (is_midpoint_lt_even || is_inexact_gt_midpoint) {
  2435. Cstar.w[0] = Cstar.w[0] - 1;
  2436. } else {
  2437. ; // exact, the result is already correct
  2438. }
  2439. res = Cstar.w[0];
  2440. } else if (exp == 0) {
  2441. // 1 <= q <= 10
  2442. // res = +C (exact)
  2443. res = C1.w[0];
  2444. } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
  2445. // res = +C * 10^exp (exact)
  2446. res = C1.w[0] * ten2k64[exp];
  2447. }
  2448. }
  2449. }
  2450. BID_RETURN (res);
  2451. }
  2452. /*****************************************************************************
  2453. * BID128_to_uint32_xint
  2454. ****************************************************************************/
  2455. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
  2456. bid128_to_uint32_xint, x)
  2457. int res;
  2458. UINT64 x_sign;
  2459. UINT64 x_exp;
  2460. int exp; // unbiased exponent
  2461. // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  2462. UINT64 tmp64, tmp64A;
  2463. BID_UI64DOUBLE tmp1;
  2464. unsigned int x_nr_bits;
  2465. int q, ind, shift;
  2466. UINT128 C1, C;
  2467. UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
  2468. UINT256 fstar;
  2469. UINT256 P256;
  2470. int is_inexact_gt_midpoint = 0;
  2471. int is_midpoint_lt_even = 0;
  2472. // unpack x
  2473. x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
  2474. x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
  2475. C1.w[1] = x.w[1] & MASK_COEFF;
  2476. C1.w[0] = x.w[0];
  2477. // check for NaN or Infinity
  2478. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  2479. // x is special
  2480. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  2481. if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
  2482. // set invalid flag
  2483. *pfpsf |= INVALID_EXCEPTION;
  2484. // return Integer Indefinite
  2485. res = 0x80000000;
  2486. } else { // x is QNaN
  2487. // set invalid flag
  2488. *pfpsf |= INVALID_EXCEPTION;
  2489. // return Integer Indefinite
  2490. res = 0x80000000;
  2491. }
  2492. BID_RETURN (res);
  2493. } else { // x is not a NaN, so it must be infinity
  2494. if (!x_sign) { // x is +inf
  2495. // set invalid flag
  2496. *pfpsf |= INVALID_EXCEPTION;
  2497. // return Integer Indefinite
  2498. res = 0x80000000;
  2499. } else { // x is -inf
  2500. // set invalid flag
  2501. *pfpsf |= INVALID_EXCEPTION;
  2502. // return Integer Indefinite
  2503. res = 0x80000000;
  2504. }
  2505. BID_RETURN (res);
  2506. }
  2507. }
  2508. // check for non-canonical values (after the check for special values)
  2509. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  2510. || (C1.w[1] == 0x0001ed09bead87c0ull
  2511. && (C1.w[0] > 0x378d8e63ffffffffull))
  2512. || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  2513. res = 0x00000000;
  2514. BID_RETURN (res);
  2515. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  2516. // x is 0
  2517. res = 0x00000000;
  2518. BID_RETURN (res);
  2519. } else { // x is not special and is not zero
  2520. // q = nr. of decimal digits in x
  2521. // determine first the nr. of bits in x
  2522. if (C1.w[1] == 0) {
  2523. if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
  2524. // split the 64-bit value in two 32-bit halves to avoid rounding errors
  2525. if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
  2526. tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
  2527. x_nr_bits =
  2528. 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2529. } else { // x < 2^32
  2530. tmp1.d = (double) (C1.w[0]); // exact conversion
  2531. x_nr_bits =
  2532. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2533. }
  2534. } else { // if x < 2^53
  2535. tmp1.d = (double) C1.w[0]; // exact conversion
  2536. x_nr_bits =
  2537. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2538. }
  2539. } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  2540. tmp1.d = (double) C1.w[1]; // exact conversion
  2541. x_nr_bits =
  2542. 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2543. }
  2544. q = nr_digits[x_nr_bits - 1].digits;
  2545. if (q == 0) {
  2546. q = nr_digits[x_nr_bits - 1].digits1;
  2547. if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  2548. || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  2549. && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  2550. q++;
  2551. }
  2552. exp = (x_exp >> 49) - 6176;
  2553. if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
  2554. // set invalid flag
  2555. *pfpsf |= INVALID_EXCEPTION;
  2556. // return Integer Indefinite
  2557. res = 0x80000000;
  2558. BID_RETURN (res);
  2559. } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
  2560. // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
  2561. // so x rounded to an integer may or may not fit in a signed 32-bit int
  2562. // the cases that do not fit are identified here; the ones that fit
  2563. // fall through and will be handled with other cases further,
  2564. // under '1 <= q + exp <= 10'
  2565. if (x_sign) { // if n < 0 and q + exp = 10
  2566. // if n <= -1 then n is too large
  2567. // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
  2568. // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a, 1<=q<=34
  2569. if (q <= 11) {
  2570. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  2571. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  2572. if (tmp64 >= 0x0aull) {
  2573. // set invalid flag
  2574. *pfpsf |= INVALID_EXCEPTION;
  2575. // return Integer Indefinite
  2576. res = 0x80000000;
  2577. BID_RETURN (res);
  2578. }
  2579. // else cases that can be rounded to a 32-bit uint fall through
  2580. // to '1 <= q + exp <= 10'
  2581. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  2582. // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
  2583. // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
  2584. // (scale 1 up)
  2585. tmp64 = 0x0aull;
  2586. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  2587. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  2588. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  2589. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  2590. }
  2591. if (C1.w[1] > C.w[1]
  2592. || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  2593. // set invalid flag
  2594. *pfpsf |= INVALID_EXCEPTION;
  2595. // return Integer Indefinite
  2596. res = 0x80000000;
  2597. BID_RETURN (res);
  2598. }
  2599. // else cases that can be rounded to a 32-bit int fall through
  2600. // to '1 <= q + exp <= 10'
  2601. }
  2602. } else { // if n > 0 and q + exp = 10
  2603. // if n >= 2^32 then n is too large
  2604. // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
  2605. // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
  2606. if (q <= 11) {
  2607. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  2608. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  2609. if (tmp64 >= 0xa00000000ull) {
  2610. // set invalid flag
  2611. *pfpsf |= INVALID_EXCEPTION;
  2612. // return Integer Indefinite
  2613. res = 0x80000000;
  2614. BID_RETURN (res);
  2615. }
  2616. // else cases that can be rounded to a 32-bit uint fall through
  2617. // to '1 <= q + exp <= 10'
  2618. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  2619. // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
  2620. // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
  2621. // (scale 2^32 up)
  2622. tmp64 = 0xa00000000ull;
  2623. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  2624. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  2625. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  2626. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  2627. }
  2628. if (C1.w[1] > C.w[1]
  2629. || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  2630. // set invalid flag
  2631. *pfpsf |= INVALID_EXCEPTION;
  2632. // return Integer Indefinite
  2633. res = 0x80000000;
  2634. BID_RETURN (res);
  2635. }
  2636. // else cases that can be rounded to a 32-bit int fall through
  2637. // to '1 <= q + exp <= 10'
  2638. }
  2639. }
  2640. }
  2641. // n is not too large to be converted to uint32: -2^32 < n < 2^32
  2642. // Note: some of the cases tested for above fall through to this point
  2643. if ((q + exp) <= 0) {
  2644. // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
  2645. // set inexact flag
  2646. *pfpsf |= INEXACT_EXCEPTION;
  2647. // return 0
  2648. res = 0x00000000;
  2649. BID_RETURN (res);
  2650. } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
  2651. // x = d(0)...d(k).d(k+1)..., k >= 0, d(0) != 0
  2652. if (x_sign) { // x <= -1
  2653. // set invalid flag
  2654. *pfpsf |= INVALID_EXCEPTION;
  2655. // return Integer Indefinite
  2656. res = 0x80000000;
  2657. BID_RETURN (res);
  2658. }
  2659. // x > 0 from this point on
  2660. // 1 <= x < 2^32 so x can be rounded to zero to a 32-bit unsigned integer
  2661. if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
  2662. ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
  2663. // chop off ind digits from the lower part of C1
  2664. // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
  2665. tmp64 = C1.w[0];
  2666. if (ind <= 19) {
  2667. C1.w[0] = C1.w[0] + midpoint64[ind - 1];
  2668. } else {
  2669. C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
  2670. C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
  2671. }
  2672. if (C1.w[0] < tmp64)
  2673. C1.w[1]++;
  2674. // calculate C* and f*
  2675. // C* is actually floor(C*) in this case
  2676. // C* and f* need shifting and masking, as shown by
  2677. // shiftright128[] and maskhigh128[]
  2678. // 1 <= x <= 33
  2679. // kx = 10^(-x) = ten2mk128[ind - 1]
  2680. // C* = (C1 + 1/2 * 10^x) * 10^(-x)
  2681. // the approximation of 10^(-x) was rounded up to 118 bits
  2682. __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  2683. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  2684. Cstar.w[1] = P256.w[3];
  2685. Cstar.w[0] = P256.w[2];
  2686. fstar.w[3] = 0;
  2687. fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  2688. fstar.w[1] = P256.w[1];
  2689. fstar.w[0] = P256.w[0];
  2690. } else { // 22 <= ind - 1 <= 33
  2691. Cstar.w[1] = 0;
  2692. Cstar.w[0] = P256.w[3];
  2693. fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  2694. fstar.w[2] = P256.w[2];
  2695. fstar.w[1] = P256.w[1];
  2696. fstar.w[0] = P256.w[0];
  2697. }
  2698. // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  2699. // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  2700. // if (0 < f* < 10^(-x)) then the result is a midpoint
  2701. // if floor(C*) is even then C* = floor(C*) - logical right
  2702. // shift; C* has p decimal digits, correct by Prop. 1)
  2703. // else if floor(C*) is odd C* = floor(C*)-1 (logical right
  2704. // shift; C* has p decimal digits, correct by Pr. 1)
  2705. // else
  2706. // C* = floor(C*) (logical right shift; C has p decimal digits,
  2707. // correct by Property 1)
  2708. // n = C* * 10^(e+x)
  2709. // shift right C* by Ex-128 = shiftright128[ind]
  2710. shift = shiftright128[ind - 1]; // 0 <= shift <= 102
  2711. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  2712. Cstar.w[0] =
  2713. (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  2714. // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  2715. } else { // 22 <= ind - 1 <= 33
  2716. Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
  2717. }
  2718. // determine inexactness of the rounding of C*
  2719. // if (0 < f* - 1/2 < 10^(-x)) then
  2720. // the result is exact
  2721. // else // if (f* - 1/2 > T*) then
  2722. // the result is inexact
  2723. if (ind - 1 <= 2) {
  2724. if (fstar.w[1] > 0x8000000000000000ull ||
  2725. (fstar.w[1] == 0x8000000000000000ull
  2726. && fstar.w[0] > 0x0ull)) {
  2727. // f* > 1/2 and the result may be exact
  2728. tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
  2729. if (tmp64 > ten2mk128trunc[ind - 1].w[1]
  2730. || (tmp64 == ten2mk128trunc[ind - 1].w[1]
  2731. && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
  2732. // set the inexact flag
  2733. *pfpsf |= INEXACT_EXCEPTION;
  2734. } // else the result is exact
  2735. } else { // the result is inexact; f2* <= 1/2
  2736. // set the inexact flag
  2737. *pfpsf |= INEXACT_EXCEPTION;
  2738. is_inexact_gt_midpoint = 1;
  2739. }
  2740. } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
  2741. if (fstar.w[3] > 0x0 ||
  2742. (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
  2743. (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
  2744. (fstar.w[1] || fstar.w[0]))) {
  2745. // f2* > 1/2 and the result may be exact
  2746. // Calculate f2* - 1/2
  2747. tmp64 = fstar.w[2] - onehalf128[ind - 1];
  2748. tmp64A = fstar.w[3];
  2749. if (tmp64 > fstar.w[2])
  2750. tmp64A--;
  2751. if (tmp64A || tmp64
  2752. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  2753. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  2754. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  2755. // set the inexact flag
  2756. *pfpsf |= INEXACT_EXCEPTION;
  2757. } // else the result is exact
  2758. } else { // the result is inexact; f2* <= 1/2
  2759. // set the inexact flag
  2760. *pfpsf |= INEXACT_EXCEPTION;
  2761. is_inexact_gt_midpoint = 1;
  2762. }
  2763. } else { // if 22 <= ind <= 33
  2764. if (fstar.w[3] > onehalf128[ind - 1] ||
  2765. (fstar.w[3] == onehalf128[ind - 1] &&
  2766. (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
  2767. // f2* > 1/2 and the result may be exact
  2768. // Calculate f2* - 1/2
  2769. tmp64 = fstar.w[3] - onehalf128[ind - 1];
  2770. if (tmp64 || fstar.w[2]
  2771. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  2772. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  2773. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  2774. // set the inexact flag
  2775. *pfpsf |= INEXACT_EXCEPTION;
  2776. } // else the result is exact
  2777. } else { // the result is inexact; f2* <= 1/2
  2778. // set the inexact flag
  2779. *pfpsf |= INEXACT_EXCEPTION;
  2780. is_inexact_gt_midpoint = 1;
  2781. }
  2782. }
  2783. // if the result was a midpoint it was rounded away from zero, so
  2784. // it will need a correction
  2785. // check for midpoints
  2786. if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
  2787. && (fstar.w[1] || fstar.w[0])
  2788. && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
  2789. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  2790. && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
  2791. // the result is a midpoint; round to nearest
  2792. if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
  2793. // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
  2794. Cstar.w[0]--; // Cstar.w[0] is now even
  2795. is_inexact_gt_midpoint = 0;
  2796. } else { // else MP in [ODD, EVEN]
  2797. is_midpoint_lt_even = 1;
  2798. is_inexact_gt_midpoint = 0;
  2799. }
  2800. }
  2801. // general correction for RZ
  2802. if (is_midpoint_lt_even || is_inexact_gt_midpoint) {
  2803. Cstar.w[0] = Cstar.w[0] - 1;
  2804. } else {
  2805. ; // exact, the result is already correct
  2806. }
  2807. res = Cstar.w[0];
  2808. } else if (exp == 0) {
  2809. // 1 <= q <= 10
  2810. // res = +C (exact)
  2811. res = C1.w[0];
  2812. } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
  2813. // res = +C * 10^exp (exact)
  2814. res = C1.w[0] * ten2k64[exp];
  2815. }
  2816. }
  2817. }
  2818. BID_RETURN (res);
  2819. }
  2820. /*****************************************************************************
  2821. * BID128_to_uint32_rninta
  2822. ****************************************************************************/
  2823. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
  2824. bid128_to_uint32_rninta, x)
  2825. unsigned int res;
  2826. UINT64 x_sign;
  2827. UINT64 x_exp;
  2828. int exp; // unbiased exponent
  2829. // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  2830. UINT64 tmp64;
  2831. BID_UI64DOUBLE tmp1;
  2832. unsigned int x_nr_bits;
  2833. int q, ind, shift;
  2834. UINT128 C1, C;
  2835. UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
  2836. UINT256 P256;
  2837. // unpack x
  2838. x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
  2839. x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
  2840. C1.w[1] = x.w[1] & MASK_COEFF;
  2841. C1.w[0] = x.w[0];
  2842. // check for NaN or Infinity
  2843. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  2844. // x is special
  2845. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  2846. if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
  2847. // set invalid flag
  2848. *pfpsf |= INVALID_EXCEPTION;
  2849. // return Integer Indefinite
  2850. res = 0x80000000;
  2851. } else { // x is QNaN
  2852. // set invalid flag
  2853. *pfpsf |= INVALID_EXCEPTION;
  2854. // return Integer Indefinite
  2855. res = 0x80000000;
  2856. }
  2857. BID_RETURN (res);
  2858. } else { // x is not a NaN, so it must be infinity
  2859. if (!x_sign) { // x is +inf
  2860. // set invalid flag
  2861. *pfpsf |= INVALID_EXCEPTION;
  2862. // return Integer Indefinite
  2863. res = 0x80000000;
  2864. } else { // x is -inf
  2865. // set invalid flag
  2866. *pfpsf |= INVALID_EXCEPTION;
  2867. // return Integer Indefinite
  2868. res = 0x80000000;
  2869. }
  2870. BID_RETURN (res);
  2871. }
  2872. }
  2873. // check for non-canonical values (after the check for special values)
  2874. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  2875. || (C1.w[1] == 0x0001ed09bead87c0ull
  2876. && (C1.w[0] > 0x378d8e63ffffffffull))
  2877. || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  2878. res = 0x00000000;
  2879. BID_RETURN (res);
  2880. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  2881. // x is 0
  2882. res = 0x00000000;
  2883. BID_RETURN (res);
  2884. } else { // x is not special and is not zero
  2885. // q = nr. of decimal digits in x
  2886. // determine first the nr. of bits in x
  2887. if (C1.w[1] == 0) {
  2888. if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
  2889. // split the 64-bit value in two 32-bit halves to avoid rounding errors
  2890. if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
  2891. tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
  2892. x_nr_bits =
  2893. 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2894. } else { // x < 2^32
  2895. tmp1.d = (double) (C1.w[0]); // exact conversion
  2896. x_nr_bits =
  2897. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2898. }
  2899. } else { // if x < 2^53
  2900. tmp1.d = (double) C1.w[0]; // exact conversion
  2901. x_nr_bits =
  2902. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2903. }
  2904. } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  2905. tmp1.d = (double) C1.w[1]; // exact conversion
  2906. x_nr_bits =
  2907. 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  2908. }
  2909. q = nr_digits[x_nr_bits - 1].digits;
  2910. if (q == 0) {
  2911. q = nr_digits[x_nr_bits - 1].digits1;
  2912. if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  2913. || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  2914. && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  2915. q++;
  2916. }
  2917. exp = (x_exp >> 49) - 6176;
  2918. if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
  2919. // set invalid flag
  2920. *pfpsf |= INVALID_EXCEPTION;
  2921. // return Integer Indefinite
  2922. res = 0x80000000;
  2923. BID_RETURN (res);
  2924. } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
  2925. // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
  2926. // so x rounded to an integer may or may not fit in a signed 32-bit int
  2927. // the cases that do not fit are identified here; the ones that fit
  2928. // fall through and will be handled with other cases further,
  2929. // under '1 <= q + exp <= 10'
  2930. if (x_sign) { // if n < 0 and q + exp = 10
  2931. // if n <= -1/2 then n is too large
  2932. // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1/2
  2933. // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05, 1<=q<=34
  2934. if (q <= 11) {
  2935. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  2936. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  2937. if (tmp64 >= 0x05ull) {
  2938. // set invalid flag
  2939. *pfpsf |= INVALID_EXCEPTION;
  2940. // return Integer Indefinite
  2941. res = 0x80000000;
  2942. BID_RETURN (res);
  2943. }
  2944. // else cases that can be rounded to a 32-bit int fall through
  2945. // to '1 <= q + exp <= 10'
  2946. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  2947. // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05 <=>
  2948. // C >= 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
  2949. // (scale 1/2 up)
  2950. tmp64 = 0x05ull;
  2951. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  2952. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  2953. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  2954. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  2955. }
  2956. if (C1.w[1] > C.w[1]
  2957. || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  2958. // set invalid flag
  2959. *pfpsf |= INVALID_EXCEPTION;
  2960. // return Integer Indefinite
  2961. res = 0x80000000;
  2962. BID_RETURN (res);
  2963. }
  2964. // else cases that can be rounded to a 32-bit int fall through
  2965. // to '1 <= q + exp <= 10'
  2966. }
  2967. } else { // if n > 0 and q + exp = 10
  2968. // if n >= 2^32 - 1/2 then n is too large
  2969. // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
  2970. // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
  2971. if (q <= 11) {
  2972. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  2973. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  2974. if (tmp64 >= 0x9fffffffbull) {
  2975. // set invalid flag
  2976. *pfpsf |= INVALID_EXCEPTION;
  2977. // return Integer Indefinite
  2978. res = 0x80000000;
  2979. BID_RETURN (res);
  2980. }
  2981. // else cases that can be rounded to a 32-bit int fall through
  2982. // to '1 <= q + exp <= 10'
  2983. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  2984. // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
  2985. // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
  2986. // (scale 2^32-1/2 up)
  2987. tmp64 = 0x9fffffffbull;
  2988. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  2989. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  2990. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  2991. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  2992. }
  2993. if (C1.w[1] > C.w[1]
  2994. || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  2995. // set invalid flag
  2996. *pfpsf |= INVALID_EXCEPTION;
  2997. // return Integer Indefinite
  2998. res = 0x80000000;
  2999. BID_RETURN (res);
  3000. }
  3001. // else cases that can be rounded to a 32-bit int fall through
  3002. // to '1 <= q + exp <= 10'
  3003. }
  3004. }
  3005. }
  3006. // n is not too large to be converted to int32: -1/2 < n < 2^32 - 1/2
  3007. // Note: some of the cases tested for above fall through to this point
  3008. if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
  3009. // return 0
  3010. res = 0x00000000;
  3011. BID_RETURN (res);
  3012. } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
  3013. // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
  3014. // res = 0
  3015. // else if x > 0
  3016. // res = +1
  3017. // else // if x < 0
  3018. // invalid exc
  3019. ind = q - 1;
  3020. if (ind <= 18) { // 0 <= ind <= 18
  3021. if ((C1.w[1] == 0) && (C1.w[0] < midpoint64[ind])) {
  3022. res = 0x00000000; // return 0
  3023. } else if (!x_sign) { // n > 0
  3024. res = 0x00000001; // return +1
  3025. } else {
  3026. res = 0x80000000;
  3027. *pfpsf |= INVALID_EXCEPTION;
  3028. }
  3029. } else { // 19 <= ind <= 33
  3030. if ((C1.w[1] < midpoint128[ind - 19].w[1])
  3031. || ((C1.w[1] == midpoint128[ind - 19].w[1])
  3032. && (C1.w[0] < midpoint128[ind - 19].w[0]))) {
  3033. res = 0x00000000; // return 0
  3034. } else if (!x_sign) { // n > 0
  3035. res = 0x00000001; // return +1
  3036. } else {
  3037. res = 0x80000000;
  3038. *pfpsf |= INVALID_EXCEPTION;
  3039. }
  3040. }
  3041. } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
  3042. if (x_sign) { // x <= -1
  3043. // set invalid flag
  3044. *pfpsf |= INVALID_EXCEPTION;
  3045. // return Integer Indefinite
  3046. res = 0x80000000;
  3047. BID_RETURN (res);
  3048. }
  3049. // 1 <= x < 2^31-1/2 so x can be rounded
  3050. // to nearest-away to a 32-bit signed integer
  3051. if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
  3052. ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
  3053. // chop off ind digits from the lower part of C1
  3054. // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
  3055. tmp64 = C1.w[0];
  3056. if (ind <= 19) {
  3057. C1.w[0] = C1.w[0] + midpoint64[ind - 1];
  3058. } else {
  3059. C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
  3060. C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
  3061. }
  3062. if (C1.w[0] < tmp64)
  3063. C1.w[1]++;
  3064. // calculate C* and f*
  3065. // C* is actually floor(C*) in this case
  3066. // C* and f* need shifting and masking, as shown by
  3067. // shiftright128[] and maskhigh128[]
  3068. // 1 <= x <= 33
  3069. // kx = 10^(-x) = ten2mk128[ind - 1]
  3070. // C* = (C1 + 1/2 * 10^x) * 10^(-x)
  3071. // the approximation of 10^(-x) was rounded up to 118 bits
  3072. __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  3073. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  3074. Cstar.w[1] = P256.w[3];
  3075. Cstar.w[0] = P256.w[2];
  3076. } else { // 22 <= ind - 1 <= 33
  3077. Cstar.w[1] = 0;
  3078. Cstar.w[0] = P256.w[3];
  3079. }
  3080. // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  3081. // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  3082. // if (0 < f* < 10^(-x)) then the result is a midpoint
  3083. // if floor(C*) is even then C* = floor(C*) - logical right
  3084. // shift; C* has p decimal digits, correct by Prop. 1)
  3085. // else if floor(C*) is odd C* = floor(C*)-1 (logical right
  3086. // shift; C* has p decimal digits, correct by Pr. 1)
  3087. // else
  3088. // C* = floor(C*) (logical right shift; C has p decimal digits,
  3089. // correct by Property 1)
  3090. // n = C* * 10^(e+x)
  3091. // shift right C* by Ex-128 = shiftright128[ind]
  3092. shift = shiftright128[ind - 1]; // 0 <= shift <= 102
  3093. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  3094. Cstar.w[0] =
  3095. (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  3096. // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  3097. } else { // 22 <= ind - 1 <= 33
  3098. Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
  3099. }
  3100. // if the result was a midpoint, it was already rounded away from zero
  3101. res = Cstar.w[0]; // always positive
  3102. // no need to check for midpoints - already rounded away from zero!
  3103. } else if (exp == 0) {
  3104. // 1 <= q <= 10
  3105. // res = +C (exact)
  3106. res = C1.w[0];
  3107. } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
  3108. // res = +C * 10^exp (exact)
  3109. res = C1.w[0] * ten2k64[exp];
  3110. }
  3111. }
  3112. }
  3113. BID_RETURN (res);
  3114. }
  3115. /*****************************************************************************
  3116. * BID128_to_uint32_xrninta
  3117. ****************************************************************************/
  3118. BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
  3119. bid128_to_uint32_xrninta, x)
  3120. unsigned int res;
  3121. UINT64 x_sign;
  3122. UINT64 x_exp;
  3123. int exp; // unbiased exponent
  3124. // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
  3125. UINT64 tmp64, tmp64A;
  3126. BID_UI64DOUBLE tmp1;
  3127. unsigned int x_nr_bits;
  3128. int q, ind, shift;
  3129. UINT128 C1, C;
  3130. UINT128 Cstar; // C* represents up to 34 decimal digits ~ 113 bits
  3131. UINT256 fstar;
  3132. UINT256 P256;
  3133. unsigned int tmp_inexact = 0;
  3134. // unpack x
  3135. x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
  3136. x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bit positions
  3137. C1.w[1] = x.w[1] & MASK_COEFF;
  3138. C1.w[0] = x.w[0];
  3139. // check for NaN or Infinity
  3140. if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  3141. // x is special
  3142. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  3143. if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
  3144. // set invalid flag
  3145. *pfpsf |= INVALID_EXCEPTION;
  3146. // return Integer Indefinite
  3147. res = 0x80000000;
  3148. } else { // x is QNaN
  3149. // set invalid flag
  3150. *pfpsf |= INVALID_EXCEPTION;
  3151. // return Integer Indefinite
  3152. res = 0x80000000;
  3153. }
  3154. BID_RETURN (res);
  3155. } else { // x is not a NaN, so it must be infinity
  3156. if (!x_sign) { // x is +inf
  3157. // set invalid flag
  3158. *pfpsf |= INVALID_EXCEPTION;
  3159. // return Integer Indefinite
  3160. res = 0x80000000;
  3161. } else { // x is -inf
  3162. // set invalid flag
  3163. *pfpsf |= INVALID_EXCEPTION;
  3164. // return Integer Indefinite
  3165. res = 0x80000000;
  3166. }
  3167. BID_RETURN (res);
  3168. }
  3169. }
  3170. // check for non-canonical values (after the check for special values)
  3171. if ((C1.w[1] > 0x0001ed09bead87c0ull)
  3172. || (C1.w[1] == 0x0001ed09bead87c0ull
  3173. && (C1.w[0] > 0x378d8e63ffffffffull))
  3174. || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
  3175. res = 0x00000000;
  3176. BID_RETURN (res);
  3177. } else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
  3178. // x is 0
  3179. res = 0x00000000;
  3180. BID_RETURN (res);
  3181. } else { // x is not special and is not zero
  3182. // q = nr. of decimal digits in x
  3183. // determine first the nr. of bits in x
  3184. if (C1.w[1] == 0) {
  3185. if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
  3186. // split the 64-bit value in two 32-bit halves to avoid rounding errors
  3187. if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
  3188. tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
  3189. x_nr_bits =
  3190. 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  3191. } else { // x < 2^32
  3192. tmp1.d = (double) (C1.w[0]); // exact conversion
  3193. x_nr_bits =
  3194. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  3195. }
  3196. } else { // if x < 2^53
  3197. tmp1.d = (double) C1.w[0]; // exact conversion
  3198. x_nr_bits =
  3199. 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  3200. }
  3201. } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
  3202. tmp1.d = (double) C1.w[1]; // exact conversion
  3203. x_nr_bits =
  3204. 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
  3205. }
  3206. q = nr_digits[x_nr_bits - 1].digits;
  3207. if (q == 0) {
  3208. q = nr_digits[x_nr_bits - 1].digits1;
  3209. if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
  3210. || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
  3211. && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
  3212. q++;
  3213. }
  3214. exp = (x_exp >> 49) - 6176;
  3215. if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
  3216. // set invalid flag
  3217. *pfpsf |= INVALID_EXCEPTION;
  3218. // return Integer Indefinite
  3219. res = 0x80000000;
  3220. BID_RETURN (res);
  3221. } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
  3222. // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
  3223. // so x rounded to an integer may or may not fit in a signed 32-bit int
  3224. // the cases that do not fit are identified here; the ones that fit
  3225. // fall through and will be handled with other cases further,
  3226. // under '1 <= q + exp <= 10'
  3227. if (x_sign) { // if n < 0 and q + exp = 10
  3228. // if n <= -1/2 then n is too large
  3229. // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1/2
  3230. // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05, 1<=q<=34
  3231. if (q <= 11) {
  3232. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  3233. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  3234. if (tmp64 >= 0x05ull) {
  3235. // set invalid flag
  3236. *pfpsf |= INVALID_EXCEPTION;
  3237. // return Integer Indefinite
  3238. res = 0x80000000;
  3239. BID_RETURN (res);
  3240. }
  3241. // else cases that can be rounded to a 32-bit int fall through
  3242. // to '1 <= q + exp <= 10'
  3243. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  3244. // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05 <=>
  3245. // C >= 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
  3246. // (scale 1/2 up)
  3247. tmp64 = 0x05ull;
  3248. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  3249. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  3250. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  3251. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  3252. }
  3253. if (C1.w[1] > C.w[1]
  3254. || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  3255. // set invalid flag
  3256. *pfpsf |= INVALID_EXCEPTION;
  3257. // return Integer Indefinite
  3258. res = 0x80000000;
  3259. BID_RETURN (res);
  3260. }
  3261. // else cases that can be rounded to a 32-bit int fall through
  3262. // to '1 <= q + exp <= 10'
  3263. }
  3264. } else { // if n > 0 and q + exp = 10
  3265. // if n >= 2^32 - 1/2 then n is too large
  3266. // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
  3267. // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
  3268. if (q <= 11) {
  3269. tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
  3270. // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
  3271. if (tmp64 >= 0x9fffffffbull) {
  3272. // set invalid flag
  3273. *pfpsf |= INVALID_EXCEPTION;
  3274. // return Integer Indefinite
  3275. res = 0x80000000;
  3276. BID_RETURN (res);
  3277. }
  3278. // else cases that can be rounded to a 32-bit int fall through
  3279. // to '1 <= q + exp <= 10'
  3280. } else { // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
  3281. // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
  3282. // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
  3283. // (scale 2^32-1/2 up)
  3284. tmp64 = 0x9fffffffbull;
  3285. if (q - 11 <= 19) { // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
  3286. __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
  3287. } else { // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
  3288. __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
  3289. }
  3290. if (C1.w[1] > C.w[1]
  3291. || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
  3292. // set invalid flag
  3293. *pfpsf |= INVALID_EXCEPTION;
  3294. // return Integer Indefinite
  3295. res = 0x80000000;
  3296. BID_RETURN (res);
  3297. }
  3298. // else cases that can be rounded to a 32-bit int fall through
  3299. // to '1 <= q + exp <= 10'
  3300. }
  3301. }
  3302. }
  3303. // n is not too large to be converted to int32: -1/2 < n < 2^32 - 1/2
  3304. // Note: some of the cases tested for above fall through to this point
  3305. if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
  3306. // set inexact flag
  3307. *pfpsf |= INEXACT_EXCEPTION;
  3308. // return 0
  3309. res = 0x00000000;
  3310. BID_RETURN (res);
  3311. } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
  3312. // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
  3313. // res = 0
  3314. // else if x > 0
  3315. // res = +1
  3316. // else // if x < 0
  3317. // invalid exc
  3318. ind = q - 1;
  3319. if (ind <= 18) { // 0 <= ind <= 18
  3320. if ((C1.w[1] == 0) && (C1.w[0] < midpoint64[ind])) {
  3321. res = 0x00000000; // return 0
  3322. } else if (!x_sign) { // n > 0
  3323. res = 0x00000001; // return +1
  3324. } else {
  3325. res = 0x80000000;
  3326. *pfpsf |= INVALID_EXCEPTION;
  3327. BID_RETURN (res);
  3328. }
  3329. } else { // 19 <= ind <= 33
  3330. if ((C1.w[1] < midpoint128[ind - 19].w[1])
  3331. || ((C1.w[1] == midpoint128[ind - 19].w[1])
  3332. && (C1.w[0] < midpoint128[ind - 19].w[0]))) {
  3333. res = 0x00000000; // return 0
  3334. } else if (!x_sign) { // n > 0
  3335. res = 0x00000001; // return +1
  3336. } else {
  3337. res = 0x80000000;
  3338. *pfpsf |= INVALID_EXCEPTION;
  3339. BID_RETURN (res);
  3340. }
  3341. }
  3342. // set inexact flag
  3343. *pfpsf |= INEXACT_EXCEPTION;
  3344. } else { // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
  3345. if (x_sign) { // x <= -1
  3346. // set invalid flag
  3347. *pfpsf |= INVALID_EXCEPTION;
  3348. // return Integer Indefinite
  3349. res = 0x80000000;
  3350. BID_RETURN (res);
  3351. }
  3352. // 1 <= x < 2^31-1/2 so x can be rounded
  3353. // to nearest-away to a 32-bit signed integer
  3354. if (exp < 0) { // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
  3355. ind = -exp; // 1 <= ind <= 33; ind is a synonym for 'x'
  3356. // chop off ind digits from the lower part of C1
  3357. // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
  3358. tmp64 = C1.w[0];
  3359. if (ind <= 19) {
  3360. C1.w[0] = C1.w[0] + midpoint64[ind - 1];
  3361. } else {
  3362. C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
  3363. C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
  3364. }
  3365. if (C1.w[0] < tmp64)
  3366. C1.w[1]++;
  3367. // calculate C* and f*
  3368. // C* is actually floor(C*) in this case
  3369. // C* and f* need shifting and masking, as shown by
  3370. // shiftright128[] and maskhigh128[]
  3371. // 1 <= x <= 33
  3372. // kx = 10^(-x) = ten2mk128[ind - 1]
  3373. // C* = (C1 + 1/2 * 10^x) * 10^(-x)
  3374. // the approximation of 10^(-x) was rounded up to 118 bits
  3375. __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
  3376. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  3377. Cstar.w[1] = P256.w[3];
  3378. Cstar.w[0] = P256.w[2];
  3379. fstar.w[3] = 0;
  3380. fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
  3381. fstar.w[1] = P256.w[1];
  3382. fstar.w[0] = P256.w[0];
  3383. } else { // 22 <= ind - 1 <= 33
  3384. Cstar.w[1] = 0;
  3385. Cstar.w[0] = P256.w[3];
  3386. fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
  3387. fstar.w[2] = P256.w[2];
  3388. fstar.w[1] = P256.w[1];
  3389. fstar.w[0] = P256.w[0];
  3390. }
  3391. // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
  3392. // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
  3393. // if (0 < f* < 10^(-x)) then the result is a midpoint
  3394. // if floor(C*) is even then C* = floor(C*) - logical right
  3395. // shift; C* has p decimal digits, correct by Prop. 1)
  3396. // else if floor(C*) is odd C* = floor(C*)-1 (logical right
  3397. // shift; C* has p decimal digits, correct by Pr. 1)
  3398. // else
  3399. // C* = floor(C*) (logical right shift; C has p decimal digits,
  3400. // correct by Property 1)
  3401. // n = C* * 10^(e+x)
  3402. // shift right C* by Ex-128 = shiftright128[ind]
  3403. shift = shiftright128[ind - 1]; // 0 <= shift <= 102
  3404. if (ind - 1 <= 21) { // 0 <= ind - 1 <= 21
  3405. Cstar.w[0] =
  3406. (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
  3407. // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
  3408. } else { // 22 <= ind - 1 <= 33
  3409. Cstar.w[0] = (Cstar.w[0] >> (shift - 64)); // 2 <= shift - 64 <= 38
  3410. }
  3411. // if the result was a midpoint, it was already rounded away from zero
  3412. // determine inexactness of the rounding of C*
  3413. // if (0 < f* - 1/2 < 10^(-x)) then
  3414. // the result is exact
  3415. // else // if (f* - 1/2 > T*) then
  3416. // the result is inexact
  3417. if (ind - 1 <= 2) {
  3418. if (fstar.w[1] > 0x8000000000000000ull ||
  3419. (fstar.w[1] == 0x8000000000000000ull
  3420. && fstar.w[0] > 0x0ull)) {
  3421. // f* > 1/2 and the result may be exact
  3422. tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
  3423. if (tmp64 > ten2mk128trunc[ind - 1].w[1]
  3424. || (tmp64 == ten2mk128trunc[ind - 1].w[1]
  3425. && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
  3426. // set the inexact flag
  3427. // *pfpsf |= INEXACT_EXCEPTION;
  3428. tmp_inexact = 1;
  3429. } // else the result is exact
  3430. } else { // the result is inexact; f2* <= 1/2
  3431. // set the inexact flag
  3432. // *pfpsf |= INEXACT_EXCEPTION;
  3433. tmp_inexact = 1;
  3434. }
  3435. } else if (ind - 1 <= 21) { // if 3 <= ind <= 21
  3436. if (fstar.w[3] > 0x0 ||
  3437. (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
  3438. (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
  3439. (fstar.w[1] || fstar.w[0]))) {
  3440. // f2* > 1/2 and the result may be exact
  3441. // Calculate f2* - 1/2
  3442. tmp64 = fstar.w[2] - onehalf128[ind - 1];
  3443. tmp64A = fstar.w[3];
  3444. if (tmp64 > fstar.w[2])
  3445. tmp64A--;
  3446. if (tmp64A || tmp64
  3447. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  3448. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  3449. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  3450. // set the inexact flag
  3451. // *pfpsf |= INEXACT_EXCEPTION;
  3452. tmp_inexact = 1;
  3453. } // else the result is exact
  3454. } else { // the result is inexact; f2* <= 1/2
  3455. // set the inexact flag
  3456. // *pfpsf |= INEXACT_EXCEPTION;
  3457. tmp_inexact = 1;
  3458. }
  3459. } else { // if 22 <= ind <= 33
  3460. if (fstar.w[3] > onehalf128[ind - 1] ||
  3461. (fstar.w[3] == onehalf128[ind - 1] &&
  3462. (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
  3463. // f2* > 1/2 and the result may be exact
  3464. // Calculate f2* - 1/2
  3465. tmp64 = fstar.w[3] - onehalf128[ind - 1];
  3466. if (tmp64 || fstar.w[2]
  3467. || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
  3468. || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
  3469. && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
  3470. // set the inexact flag
  3471. // *pfpsf |= INEXACT_EXCEPTION;
  3472. tmp_inexact = 1;
  3473. } // else the result is exact
  3474. } else { // the result is inexact; f2* <= 1/2
  3475. // set the inexact flag
  3476. // *pfpsf |= INEXACT_EXCEPTION;
  3477. tmp_inexact = 1;
  3478. }
  3479. }
  3480. // no need to check for midpoints - already rounded away from zero!
  3481. res = Cstar.w[0]; // the result is positive
  3482. if (tmp_inexact)
  3483. *pfpsf |= INEXACT_EXCEPTION;
  3484. } else if (exp == 0) {
  3485. // 1 <= q <= 10
  3486. // res = +C (exact)
  3487. res = C1.w[0];
  3488. } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
  3489. // res = +C * 10^exp (exact)
  3490. res = C1.w[0] * ten2k64[exp];
  3491. }
  3492. }
  3493. }
  3494. BID_RETURN (res);
  3495. }