12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964196519661967196819691970197119721973197419751976197719781979198019811982198319841985198619871988198919901991199219931994199519961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018201920202021202220232024202520262027202820292030203120322033203420352036203720382039204020412042204320442045204620472048204920502051205220532054205520562057205820592060206120622063206420652066206720682069207020712072207320742075207620772078207920802081208220832084208520862087208820892090209120922093209420952096209720982099210021012102210321042105210621072108210921102111211221132114211521162117211821192120212121222123212421252126212721282129213021312132213321342135213621372138213921402141214221432144214521462147214821492150215121522153215421552156215721582159216021612162216321642165216621672168216921702171217221732174217521762177217821792180218121822183218421852186218721882189219021912192219321942195219621972198219922002201220222032204220522062207220822092210221122122213221422152216221722182219222022212222222322242225222622272228222922302231223222332234223522362237223822392240224122422243224422452246224722482249225022512252225322542255225622572258225922602261226222632264226522662267226822692270227122722273 |
- /* Copyright (C) 2007-2022 Free Software Foundation, Inc.
- This file is part of GCC.
- GCC is free software; you can redistribute it and/or modify it under
- the terms of the GNU General Public License as published by the Free
- Software Foundation; either version 3, or (at your option) any later
- version.
- GCC is distributed in the hope that it will be useful, but WITHOUT ANY
- WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
- Under Section 7 of GPL version 3, you are granted additional
- permissions described in the GCC Runtime Library Exception, version
- 3.1, as published by the Free Software Foundation.
- You should have received a copy of the GNU General Public License and
- a copy of the GCC Runtime Library Exception along with this program;
- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
- <http://www.gnu.org/licenses/>. */
- #include "bid_internal.h"
- /*****************************************************************************
- * BID64_to_uint64_rnint
- ****************************************************************************/
- #if DECIMAL_CALL_BY_REFERENCE
- void
- bid64_to_uint64_rnint (UINT64 * pres, UINT64 * px
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- UINT64 x = *px;
- #else
- UINT64
- bid64_to_uint64_rnint (UINT64 x
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- #endif
- UINT64 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1 represents x_significand (UINT64)
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT64 C1;
- UINT128 C;
- UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
- UINT128 fstar;
- UINT128 P128;
- // check for NaN or Infinity
- if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // unpack x
- x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
- if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
- x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
- C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
- if (C1 > 9999999999999999ull) { // non-canonical
- x_exp = 0;
- C1 = 0;
- }
- } else {
- x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
- C1 = x & MASK_BINARY_SIG1;
- }
- // check for zeros (possibly from non-canonical values)
- if (C1 == 0x0ull) {
- // x is 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- // q = nr. of decimal digits in x (1 <= q <= 54)
- // determine first the nr. of bits in x
- if (C1 >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1 >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1 >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- q = nr_digits[x_nr_bits - 1].digits;
- if (q == 0) {
- q = nr_digits[x_nr_bits - 1].digits1;
- if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
- q++;
- }
- exp = x_exp - 398; // unbiased exponent
- if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
- // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
- // so x rounded to an integer may or may not fit in an unsigned 64-bit int
- // the cases that do not fit are identified here; the ones that fit
- // fall through and will be handled with other cases further,
- // under '1 <= q + exp <= 20'
- if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1/2
- // => set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else { // if n > 0 and q + exp = 20
- // if n >= 2^64 - 1/2 then n is too large
- // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
- // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=16
- if (q == 1) {
- // C * 10^20 >= 0x9fffffffffffffffb
- __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C
- if (C.w[1] > 0x09 ||
- (C.w[1] == 0x09 && C.w[0] >= 0xfffffffffffffffbull)) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
- // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffffb
- // has 21 digits
- __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]);
- if (C.w[1] > 0x09 ||
- (C.w[1] == 0x09 && C.w[0] >= 0xfffffffffffffffbull)) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- }
- }
- }
- // n is not too large to be converted to int64 if -1/2 <= n < 2^64 - 1/2
- // Note: some of the cases tested for above fall through to this point
- if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
- // return 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
- // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
- // res = 0
- // else if x > 0
- // res = +1
- // else // if x < 0
- // invalid exc
- ind = q - 1; // 0 <= ind <= 15
- if (C1 <= midpoint64[ind]) {
- res = 0x0000000000000000ull; // return 0
- } else if (!x_sign) { // n > 0
- res = 0x0000000000000001ull; // return +1
- } else { // if n < 0
- res = 0x8000000000000000ull;
- *pfpsf |= INVALID_EXCEPTION;
- BID_RETURN (res);
- }
- } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
- // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
- // to nearest to a 64-bit unsigned signed integer
- if (x_sign) { // x <= -1
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // 1 <= x < 2^64-1/2 so x can be rounded
- // to nearest to a 64-bit unsigned integer
- if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
- ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
- // chop off ind digits from the lower part of C1
- // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
- C1 = C1 + midpoint64[ind - 1];
- // calculate C* and f*
- // C* is actually floor(C*) in this case
- // C* and f* need shifting and masking, as shown by
- // shiftright128[] and maskhigh128[]
- // 1 <= x <= 15
- // kx = 10^(-x) = ten2mk64[ind - 1]
- // C* = (C1 + 1/2 * 10^x) * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 54 bits
- __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
- Cstar = P128.w[1];
- fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
- fstar.w[0] = P128.w[0];
- // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
- // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
- // if (0 < f* < 10^(-x)) then the result is a midpoint
- // if floor(C*) is even then C* = floor(C*) - logical right
- // shift; C* has p decimal digits, correct by Prop. 1)
- // else if floor(C*) is odd C* = floor(C*)-1 (logical right
- // shift; C* has p decimal digits, correct by Pr. 1)
- // else
- // C* = floor(C*) (logical right shift; C has p decimal digits,
- // correct by Property 1)
- // n = C* * 10^(e+x)
- // shift right C* by Ex-64 = shiftright128[ind]
- shift = shiftright128[ind - 1]; // 0 <= shift <= 39
- Cstar = Cstar >> shift;
- // if the result was a midpoint it was rounded away from zero, so
- // it will need a correction
- // check for midpoints
- if ((fstar.w[1] == 0) && fstar.w[0] &&
- (fstar.w[0] <= ten2mk128trunc[ind - 1].w[1])) {
- // ten2mk128trunc[ind -1].w[1] is identical to
- // ten2mk128[ind -1].w[1]
- // the result is a midpoint; round to nearest
- if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD]
- // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
- Cstar--; // Cstar is now even
- } // else MP in [ODD, EVEN]
- }
- res = Cstar; // the result is positive
- } else if (exp == 0) {
- // 1 <= q <= 10
- // res = +C (exact)
- res = C1; // the result is positive
- } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
- // res = +C * 10^exp (exact)
- res = C1 * ten2k64[exp]; // the result is positive
- }
- }
- BID_RETURN (res);
- }
- /*****************************************************************************
- * BID64_to_uint64_xrnint
- ****************************************************************************/
- #if DECIMAL_CALL_BY_REFERENCE
- void
- bid64_to_uint64_xrnint (UINT64 * pres, UINT64 * px
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- UINT64 x = *px;
- #else
- UINT64
- bid64_to_uint64_xrnint (UINT64 x
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- #endif
- UINT64 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1 represents x_significand (UINT64)
- UINT64 tmp64;
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT64 C1;
- UINT128 C;
- UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
- UINT128 fstar;
- UINT128 P128;
- // check for NaN or Infinity
- if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // unpack x
- x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
- if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
- x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
- C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
- if (C1 > 9999999999999999ull) { // non-canonical
- x_exp = 0;
- C1 = 0;
- }
- } else {
- x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
- C1 = x & MASK_BINARY_SIG1;
- }
- // check for zeros (possibly from non-canonical values)
- if (C1 == 0x0ull) {
- // x is 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- // q = nr. of decimal digits in x (1 <= q <= 54)
- // determine first the nr. of bits in x
- if (C1 >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1 >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1 >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- q = nr_digits[x_nr_bits - 1].digits;
- if (q == 0) {
- q = nr_digits[x_nr_bits - 1].digits1;
- if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
- q++;
- }
- exp = x_exp - 398; // unbiased exponent
- if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
- // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
- // so x rounded to an integer may or may not fit in an unsigned 64-bit int
- // the cases that do not fit are identified here; the ones that fit
- // fall through and will be handled with other cases further,
- // under '1 <= q + exp <= 20'
- if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1/2
- // => set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else { // if n > 0 and q + exp = 20
- // if n >= 2^64 - 1/2 then n is too large
- // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
- // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=16
- if (q == 1) {
- // C * 10^20 >= 0x9fffffffffffffffb
- __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C
- if (C.w[1] > 0x09 ||
- (C.w[1] == 0x09 && C.w[0] >= 0xfffffffffffffffbull)) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
- // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffffb
- // has 21 digits
- __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]);
- if (C.w[1] > 0x09 ||
- (C.w[1] == 0x09 && C.w[0] >= 0xfffffffffffffffbull)) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- }
- }
- }
- // n is not too large to be converted to int64 if -1/2 <= n < 2^64 - 1/2
- // Note: some of the cases tested for above fall through to this point
- if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
- // set inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- // return 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
- // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
- // res = 0
- // else if x > 0
- // res = +1
- // else // if x < 0
- // invalid exc
- ind = q - 1; // 0 <= ind <= 15
- if (C1 <= midpoint64[ind]) {
- res = 0x0000000000000000ull; // return 0
- } else if (!x_sign) { // n > 0
- res = 0x0000000000000001ull; // return +1
- } else { // if n < 0
- res = 0x8000000000000000ull;
- *pfpsf |= INVALID_EXCEPTION;
- BID_RETURN (res);
- }
- // set inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
- // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
- // to nearest to a 64-bit unsigned signed integer
- if (x_sign) { // x <= -1
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // 1 <= x < 2^64-1/2 so x can be rounded
- // to nearest to a 64-bit unsigned integer
- if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
- ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
- // chop off ind digits from the lower part of C1
- // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
- C1 = C1 + midpoint64[ind - 1];
- // calculate C* and f*
- // C* is actually floor(C*) in this case
- // C* and f* need shifting and masking, as shown by
- // shiftright128[] and maskhigh128[]
- // 1 <= x <= 15
- // kx = 10^(-x) = ten2mk64[ind - 1]
- // C* = (C1 + 1/2 * 10^x) * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 54 bits
- __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
- Cstar = P128.w[1];
- fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
- fstar.w[0] = P128.w[0];
- // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
- // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
- // if (0 < f* < 10^(-x)) then the result is a midpoint
- // if floor(C*) is even then C* = floor(C*) - logical right
- // shift; C* has p decimal digits, correct by Prop. 1)
- // else if floor(C*) is odd C* = floor(C*)-1 (logical right
- // shift; C* has p decimal digits, correct by Pr. 1)
- // else
- // C* = floor(C*) (logical right shift; C has p decimal digits,
- // correct by Property 1)
- // n = C* * 10^(e+x)
- // shift right C* by Ex-64 = shiftright128[ind]
- shift = shiftright128[ind - 1]; // 0 <= shift <= 39
- Cstar = Cstar >> shift;
- // determine inexactness of the rounding of C*
- // if (0 < f* - 1/2 < 10^(-x)) then
- // the result is exact
- // else // if (f* - 1/2 > T*) then
- // the result is inexact
- if (ind - 1 <= 2) { // fstar.w[1] is 0
- if (fstar.w[0] > 0x8000000000000000ull) {
- // f* > 1/2 and the result may be exact
- tmp64 = fstar.w[0] - 0x8000000000000000ull; // f* - 1/2
- if ((tmp64 > ten2mk128trunc[ind - 1].w[1])) {
- // ten2mk128trunc[ind -1].w[1] is identical to
- // ten2mk128[ind -1].w[1]
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- } // else the result is exact
- } else { // the result is inexact; f2* <= 1/2
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- }
- } else { // if 3 <= ind - 1 <= 14
- if (fstar.w[1] > onehalf128[ind - 1] ||
- (fstar.w[1] == onehalf128[ind - 1] && fstar.w[0])) {
- // f2* > 1/2 and the result may be exact
- // Calculate f2* - 1/2
- tmp64 = fstar.w[1] - onehalf128[ind - 1];
- if (tmp64 || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
- // ten2mk128trunc[ind -1].w[1] is identical to
- // ten2mk128[ind -1].w[1]
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- } // else the result is exact
- } else { // the result is inexact; f2* <= 1/2
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- }
- }
- // if the result was a midpoint it was rounded away from zero, so
- // it will need a correction
- // check for midpoints
- if ((fstar.w[1] == 0) && fstar.w[0] &&
- (fstar.w[0] <= ten2mk128trunc[ind - 1].w[1])) {
- // ten2mk128trunc[ind -1].w[1] is identical to
- // ten2mk128[ind -1].w[1]
- // the result is a midpoint; round to nearest
- if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD]
- // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
- Cstar--; // Cstar is now even
- } // else MP in [ODD, EVEN]
- }
- res = Cstar; // the result is positive
- } else if (exp == 0) {
- // 1 <= q <= 10
- // res = +C (exact)
- res = C1; // the result is positive
- } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
- // res = +C * 10^exp (exact)
- res = C1 * ten2k64[exp]; // the result is positive
- }
- }
- BID_RETURN (res);
- }
- /*****************************************************************************
- * BID64_to_uint64_floor
- ****************************************************************************/
- #if DECIMAL_CALL_BY_REFERENCE
- void
- bid64_to_uint64_floor (UINT64 * pres, UINT64 * px
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- UINT64 x = *px;
- #else
- UINT64
- bid64_to_uint64_floor (UINT64 x
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- #endif
- UINT64 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1 represents x_significand (UINT64)
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT64 C1;
- UINT128 C;
- UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
- UINT128 P128;
- // check for NaN or Infinity
- if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // unpack x
- x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
- if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
- x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
- C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
- if (C1 > 9999999999999999ull) { // non-canonical
- x_exp = 0;
- C1 = 0;
- }
- } else {
- x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
- C1 = x & MASK_BINARY_SIG1;
- }
- // check for zeros (possibly from non-canonical values)
- if (C1 == 0x0ull) {
- // x is 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- if (x_sign) { // if n < 0 the conversion is invalid
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // q = nr. of decimal digits in x (1 <= q <= 54)
- // determine first the nr. of bits in x
- if (C1 >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1 >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1 >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- q = nr_digits[x_nr_bits - 1].digits;
- if (q == 0) {
- q = nr_digits[x_nr_bits - 1].digits1;
- if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
- q++;
- }
- exp = x_exp - 398; // unbiased exponent
- if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
- // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
- // so x rounded to an integer may or may not fit in an unsigned 64-bit int
- // the cases that do not fit are identified here; the ones that fit
- // fall through and will be handled with other cases further,
- // under '1 <= q + exp <= 20'
- // n > 0 and q + exp = 20
- // if n >= 2^64 then n is too large
- // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65)
- // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=16
- if (q == 1) {
- // C * 10^20 >= 0xa0000000000000000
- __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C
- if (C.w[1] >= 0x0a) {
- // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
- // Note: C * 10^(21-q) has 20 or 21 digits; 0xa0000000000000000
- // has 21 digits
- __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]);
- if (C.w[1] >= 0x0a) {
- // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- }
- }
- // n is not too large to be converted to int64 if -1 < n < 2^64
- // Note: some of the cases tested for above fall through to this point
- if ((q + exp) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1)
- // return 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
- // 1 <= x < 2^64 so x can be rounded
- // to nearest to a 64-bit unsigned signed integer
- if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
- ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
- // chop off ind digits from the lower part of C1
- // C1 fits in 64 bits
- // calculate C* and f*
- // C* is actually floor(C*) in this case
- // C* and f* need shifting and masking, as shown by
- // shiftright128[] and maskhigh128[]
- // 1 <= x <= 15
- // kx = 10^(-x) = ten2mk64[ind - 1]
- // C* = C1 * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 54 bits
- __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
- Cstar = P128.w[1];
- // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
- // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
- // C* = floor(C*) (logical right shift; C has p decimal digits,
- // correct by Property 1)
- // n = C* * 10^(e+x)
- // shift right C* by Ex-64 = shiftright128[ind]
- shift = shiftright128[ind - 1]; // 0 <= shift <= 39
- Cstar = Cstar >> shift;
- res = Cstar; // the result is positive
- } else if (exp == 0) {
- // 1 <= q <= 10
- // res = +C (exact)
- res = C1; // the result is positive
- } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
- // res = +C * 10^exp (exact)
- res = C1 * ten2k64[exp]; // the result is positive
- }
- }
- BID_RETURN (res);
- }
- /*****************************************************************************
- * BID64_to_uint64_xfloor
- ****************************************************************************/
- #if DECIMAL_CALL_BY_REFERENCE
- void
- bid64_to_uint64_xfloor (UINT64 * pres, UINT64 * px
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- UINT64 x = *px;
- #else
- UINT64
- bid64_to_uint64_xfloor (UINT64 x
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- #endif
- UINT64 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1 represents x_significand (UINT64)
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT64 C1;
- UINT128 C;
- UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
- UINT128 fstar;
- UINT128 P128;
- // check for NaN or Infinity
- if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // unpack x
- x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
- if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
- x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
- C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
- if (C1 > 9999999999999999ull) { // non-canonical
- x_exp = 0;
- C1 = 0;
- }
- } else {
- x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
- C1 = x & MASK_BINARY_SIG1;
- }
- // check for zeros (possibly from non-canonical values)
- if (C1 == 0x0ull) {
- // x is 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- if (x_sign) { // if n < 0 the conversion is invalid
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // q = nr. of decimal digits in x (1 <= q <= 54)
- // determine first the nr. of bits in x
- if (C1 >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1 >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1 >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- q = nr_digits[x_nr_bits - 1].digits;
- if (q == 0) {
- q = nr_digits[x_nr_bits - 1].digits1;
- if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
- q++;
- }
- exp = x_exp - 398; // unbiased exponent
- if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
- // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
- // so x rounded to an integer may or may not fit in an unsigned 64-bit int
- // the cases that do not fit are identified here; the ones that fit
- // fall through and will be handled with other cases further,
- // under '1 <= q + exp <= 20'
- // n > 0 and q + exp = 20
- // if n >= 2^64 then n is too large
- // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65)
- // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=16
- if (q == 1) {
- // C * 10^20 >= 0xa0000000000000000
- __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C
- if (C.w[1] >= 0x0a) {
- // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
- // Note: C * 10^(21-q) has 20 or 21 digits; 0xa0000000000000000
- // has 21 digits
- __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]);
- if (C.w[1] >= 0x0a) {
- // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- }
- }
- // n is not too large to be converted to int64 if -1 < n < 2^64
- // Note: some of the cases tested for above fall through to this point
- if ((q + exp) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1)
- // set inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- // return 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
- // 1 <= x < 2^64 so x can be rounded
- // to nearest to a 64-bit unsigned signed integer
- if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
- ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
- // chop off ind digits from the lower part of C1
- // C1 fits in 64 bits
- // calculate C* and f*
- // C* is actually floor(C*) in this case
- // C* and f* need shifting and masking, as shown by
- // shiftright128[] and maskhigh128[]
- // 1 <= x <= 15
- // kx = 10^(-x) = ten2mk64[ind - 1]
- // C* = C1 * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 54 bits
- __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
- Cstar = P128.w[1];
- fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
- fstar.w[0] = P128.w[0];
- // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
- // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
- // C* = floor(C*) (logical right shift; C has p decimal digits,
- // correct by Property 1)
- // n = C* * 10^(e+x)
- // shift right C* by Ex-64 = shiftright128[ind]
- shift = shiftright128[ind - 1]; // 0 <= shift <= 39
- Cstar = Cstar >> shift;
- // determine inexactness of the rounding of C*
- // if (0 < f* < 10^(-x)) then
- // the result is exact
- // else // if (f* > T*) then
- // the result is inexact
- if (ind - 1 <= 2) { // fstar.w[1] is 0
- if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
- // ten2mk128trunc[ind -1].w[1] is identical to
- // ten2mk128[ind -1].w[1]
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- } // else the result is exact
- } else { // if 3 <= ind - 1 <= 14
- if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
- // ten2mk128trunc[ind -1].w[1] is identical to
- // ten2mk128[ind -1].w[1]
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- } // else the result is exact
- }
- res = Cstar; // the result is positive
- } else if (exp == 0) {
- // 1 <= q <= 10
- // res = +C (exact)
- res = C1; // the result is positive
- } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
- // res = +C * 10^exp (exact)
- res = C1 * ten2k64[exp]; // the result is positive
- }
- }
- BID_RETURN (res);
- }
- /*****************************************************************************
- * BID64_to_uint64_ceil
- ****************************************************************************/
- #if DECIMAL_CALL_BY_REFERENCE
- void
- bid64_to_uint64_ceil (UINT64 * pres, UINT64 * px
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- UINT64 x = *px;
- #else
- UINT64
- bid64_to_uint64_ceil (UINT64 x
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- #endif
- UINT64 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1 represents x_significand (UINT64)
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT64 C1;
- UINT128 C;
- UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
- UINT128 fstar;
- UINT128 P128;
- // check for NaN or Infinity
- if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // unpack x
- x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
- if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
- x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
- C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
- if (C1 > 9999999999999999ull) { // non-canonical
- x_exp = 0;
- C1 = 0;
- }
- } else {
- x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
- C1 = x & MASK_BINARY_SIG1;
- }
- // check for zeros (possibly from non-canonical values)
- if (C1 == 0x0ull) {
- // x is 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- // q = nr. of decimal digits in x (1 <= q <= 54)
- // determine first the nr. of bits in x
- if (C1 >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1 >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1 >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- q = nr_digits[x_nr_bits - 1].digits;
- if (q == 0) {
- q = nr_digits[x_nr_bits - 1].digits1;
- if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
- q++;
- }
- exp = x_exp - 398; // unbiased exponent
- if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
- // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
- // so x rounded to an integer may or may not fit in an unsigned 64-bit int
- // the cases that do not fit are identified here; the ones that fit
- // fall through and will be handled with other cases further,
- // under '1 <= q + exp <= 20'
- if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1
- // => set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else { // if n > 0 and q + exp = 20
- // if n > 2^64 - 1 then n is too large
- // <=> c(0)c(1)...c(19).c(20)...c(q-1) > 2^64 - 1
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 > 2^64 - 1
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65 - 2)
- // <=> C * 10^(21-q) > 0x9fffffffffffffff6, 1<=q<=16
- if (q == 1) {
- // C * 10^20 > 0x9fffffffffffffff6
- __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C
- if (C.w[1] > 0x09 ||
- (C.w[1] == 0x09 && C.w[0] > 0xfffffffffffffff6ull)) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
- // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffff6
- // has 21 digits
- __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]);
- if (C.w[1] > 0x09 ||
- (C.w[1] == 0x09 && C.w[0] > 0xfffffffffffffff6ull)) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- }
- }
- }
- // n is not too large to be converted to int64 if -1 < n < 2^64
- // Note: some of the cases tested for above fall through to this point
- if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
- // return 0 or 1
- if (x_sign)
- res = 0x0000000000000000ull;
- else
- res = 0x0000000000000001ull;
- BID_RETURN (res);
- } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
- // x <= -1 or 1 <= x <= 2^64 - 1 so if positive x can be rounded
- // to nearest to a 64-bit unsigned signed integer
- if (x_sign) { // x <= -1
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // 1 <= x <= 2^64 - 1 so x can be rounded
- // to nearest to a 64-bit unsigned integer
- if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
- ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
- // chop off ind digits from the lower part of C1
- // C1 fits in 64 bits
- // calculate C* and f*
- // C* is actually floor(C*) in this case
- // C* and f* need shifting and masking, as shown by
- // shiftright128[] and maskhigh128[]
- // 1 <= x <= 15
- // kx = 10^(-x) = ten2mk64[ind - 1]
- // C* = C1 * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 54 bits
- __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
- Cstar = P128.w[1];
- fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
- fstar.w[0] = P128.w[0];
- // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
- // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
- // C* = floor(C*) (logical right shift; C has p decimal digits,
- // correct by Property 1)
- // n = C* * 10^(e+x)
- // shift right C* by Ex-64 = shiftright128[ind]
- shift = shiftright128[ind - 1]; // 0 <= shift <= 39
- Cstar = Cstar >> shift;
- // determine inexactness of the rounding of C*
- // if (0 < f* < 10^(-x)) then
- // the result is exact
- // else // if (f* > T*) then
- // the result is inexact
- if (ind - 1 <= 2) { // fstar.w[1] is 0
- if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
- // ten2mk128trunc[ind -1].w[1] is identical to
- // ten2mk128[ind -1].w[1]
- Cstar++;
- } // else the result is exact
- } else { // if 3 <= ind - 1 <= 14
- if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
- // ten2mk128trunc[ind -1].w[1] is identical to
- // ten2mk128[ind -1].w[1]
- Cstar++;
- } // else the result is exact
- }
- res = Cstar; // the result is positive
- } else if (exp == 0) {
- // 1 <= q <= 10
- // res = +C (exact)
- res = C1; // the result is positive
- } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
- // res = +C * 10^exp (exact)
- res = C1 * ten2k64[exp]; // the result is positive
- }
- }
- BID_RETURN (res);
- }
- /*****************************************************************************
- * BID64_to_uint64_xceil
- ****************************************************************************/
- #if DECIMAL_CALL_BY_REFERENCE
- void
- bid64_to_uint64_xceil (UINT64 * pres, UINT64 * px
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- UINT64 x = *px;
- #else
- UINT64
- bid64_to_uint64_xceil (UINT64 x
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- #endif
- UINT64 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1 represents x_significand (UINT64)
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT64 C1;
- UINT128 C;
- UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
- UINT128 fstar;
- UINT128 P128;
- // check for NaN or Infinity
- if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // unpack x
- x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
- if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
- x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
- C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
- if (C1 > 9999999999999999ull) { // non-canonical
- x_exp = 0;
- C1 = 0;
- }
- } else {
- x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
- C1 = x & MASK_BINARY_SIG1;
- }
- // check for zeros (possibly from non-canonical values)
- if (C1 == 0x0ull) {
- // x is 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- // q = nr. of decimal digits in x (1 <= q <= 54)
- // determine first the nr. of bits in x
- if (C1 >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1 >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1 >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- q = nr_digits[x_nr_bits - 1].digits;
- if (q == 0) {
- q = nr_digits[x_nr_bits - 1].digits1;
- if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
- q++;
- }
- exp = x_exp - 398; // unbiased exponent
- if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
- // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
- // so x rounded to an integer may or may not fit in an unsigned 64-bit int
- // the cases that do not fit are identified here; the ones that fit
- // fall through and will be handled with other cases further,
- // under '1 <= q + exp <= 20'
- if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1
- // => set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else { // if n > 0 and q + exp = 20
- // if n > 2^64 - 1 then n is too large
- // <=> c(0)c(1)...c(19).c(20)...c(q-1) > 2^64 - 1
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 > 2^64 - 1
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65 - 2)
- // <=> C * 10^(21-q) > 0x9fffffffffffffff6, 1<=q<=16
- if (q == 1) {
- // C * 10^20 > 0x9fffffffffffffff6
- __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C
- if (C.w[1] > 0x09 ||
- (C.w[1] == 0x09 && C.w[0] > 0xfffffffffffffff6ull)) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
- // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffff6
- // has 21 digits
- __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]);
- if (C.w[1] > 0x09 ||
- (C.w[1] == 0x09 && C.w[0] > 0xfffffffffffffff6ull)) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- }
- }
- }
- // n is not too large to be converted to int64 if -1 < n < 2^64
- // Note: some of the cases tested for above fall through to this point
- if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
- // set inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- // return 0 or 1
- if (x_sign)
- res = 0x0000000000000000ull;
- else
- res = 0x0000000000000001ull;
- BID_RETURN (res);
- } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
- // x <= -1 or 1 <= x <= 2^64 - 1 so if positive x can be rounded
- // to nearest to a 64-bit unsigned signed integer
- if (x_sign) { // x <= -1
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // 1 <= x <= 2^64 - 1 so x can be rounded
- // to nearest to a 64-bit unsigned integer
- if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
- ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
- // chop off ind digits from the lower part of C1
- // C1 fits in 64 bits
- // calculate C* and f*
- // C* is actually floor(C*) in this case
- // C* and f* need shifting and masking, as shown by
- // shiftright128[] and maskhigh128[]
- // 1 <= x <= 15
- // kx = 10^(-x) = ten2mk64[ind - 1]
- // C* = C1 * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 54 bits
- __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
- Cstar = P128.w[1];
- fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
- fstar.w[0] = P128.w[0];
- // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
- // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
- // C* = floor(C*) (logical right shift; C has p decimal digits,
- // correct by Property 1)
- // n = C* * 10^(e+x)
- // shift right C* by Ex-64 = shiftright128[ind]
- shift = shiftright128[ind - 1]; // 0 <= shift <= 39
- Cstar = Cstar >> shift;
- // determine inexactness of the rounding of C*
- // if (0 < f* < 10^(-x)) then
- // the result is exact
- // else // if (f* > T*) then
- // the result is inexact
- if (ind - 1 <= 2) { // fstar.w[1] is 0
- if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
- // ten2mk128trunc[ind -1].w[1] is identical to
- // ten2mk128[ind -1].w[1]
- Cstar++;
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- } // else the result is exact
- } else { // if 3 <= ind - 1 <= 14
- if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
- // ten2mk128trunc[ind -1].w[1] is identical to
- // ten2mk128[ind -1].w[1]
- Cstar++;
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- } // else the result is exact
- }
- res = Cstar; // the result is positive
- } else if (exp == 0) {
- // 1 <= q <= 10
- // res = +C (exact)
- res = C1; // the result is positive
- } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
- // res = +C * 10^exp (exact)
- res = C1 * ten2k64[exp]; // the result is positive
- }
- }
- BID_RETURN (res);
- }
- /*****************************************************************************
- * BID64_to_uint64_int
- ****************************************************************************/
- #if DECIMAL_CALL_BY_REFERENCE
- void
- bid64_to_uint64_int (UINT64 * pres, UINT64 * px
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM)
- {
- UINT64 x = *px;
- #else
- UINT64
- bid64_to_uint64_int (UINT64 x
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM)
- {
- #endif
- UINT64 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1 represents x_significand (UINT64)
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT64 C1;
- UINT128 C;
- UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
- UINT128 P128;
- // check for NaN or Infinity
- if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // unpack x
- x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
- if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
- x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
- C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
- if (C1 > 9999999999999999ull) { // non-canonical
- x_exp = 0;
- C1 = 0;
- }
- } else {
- x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
- C1 = x & MASK_BINARY_SIG1;
- }
- // check for zeros (possibly from non-canonical values)
- if (C1 == 0x0ull) {
- // x is 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- // q = nr. of decimal digits in x (1 <= q <= 54)
- // determine first the nr. of bits in x
- if (C1 >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1 >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1 >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- q = nr_digits[x_nr_bits - 1].digits;
- if (q == 0) {
- q = nr_digits[x_nr_bits - 1].digits1;
- if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
- q++;
- }
- exp = x_exp - 398; // unbiased exponent
- if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
- // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
- // so x rounded to an integer may or may not fit in an unsigned 64-bit int
- // the cases that do not fit are identified here; the ones that fit
- // fall through and will be handled with other cases further,
- // under '1 <= q + exp <= 20'
- if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1
- // => set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else { // if n > 0 and q + exp = 20
- // if n >= 2^64 then n is too large
- // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65)
- // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=16
- if (q == 1) {
- // C * 10^20 >= 0xa0000000000000000
- __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C
- if (C.w[1] >= 0x0a) {
- // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
- // Note: C * 10^(21-q) has 20 or 21 digits; 0xa0000000000000000
- // has 21 digits
- __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]);
- if (C.w[1] >= 0x0a) {
- // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- }
- }
- }
- // n is not too large to be converted to int64 if -1 < n < 2^64
- // Note: some of the cases tested for above fall through to this point
- if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
- // return 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
- // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
- // to nearest to a 64-bit unsigned signed integer
- if (x_sign) { // x <= -1
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // 1 <= x < 2^64 so x can be rounded
- // to nearest to a 64-bit unsigned integer
- if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
- ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
- // chop off ind digits from the lower part of C1
- // C1 fits in 64 bits
- // calculate C* and f*
- // C* is actually floor(C*) in this case
- // C* and f* need shifting and masking, as shown by
- // shiftright128[] and maskhigh128[]
- // 1 <= x <= 15
- // kx = 10^(-x) = ten2mk64[ind - 1]
- // C* = C1 * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 54 bits
- __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
- Cstar = P128.w[1];
- // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
- // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
- // C* = floor(C*) (logical right shift; C has p decimal digits,
- // correct by Property 1)
- // n = C* * 10^(e+x)
- // shift right C* by Ex-64 = shiftright128[ind]
- shift = shiftright128[ind - 1]; // 0 <= shift <= 39
- Cstar = Cstar >> shift;
- res = Cstar; // the result is positive
- } else if (exp == 0) {
- // 1 <= q <= 10
- // res = +C (exact)
- res = C1; // the result is positive
- } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
- // res = +C * 10^exp (exact)
- res = C1 * ten2k64[exp]; // the result is positive
- }
- }
- BID_RETURN (res);
- }
- /*****************************************************************************
- * BID64_to_uint64_xint
- ****************************************************************************/
- #if DECIMAL_CALL_BY_REFERENCE
- void
- bid64_to_uint64_xint (UINT64 * pres, UINT64 * px
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- UINT64 x = *px;
- #else
- UINT64
- bid64_to_uint64_xint (UINT64 x
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- #endif
- UINT64 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1 represents x_significand (UINT64)
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT64 C1;
- UINT128 C;
- UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
- UINT128 fstar;
- UINT128 P128;
- // check for NaN or Infinity
- if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // unpack x
- x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
- if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
- x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
- C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
- if (C1 > 9999999999999999ull) { // non-canonical
- x_exp = 0;
- C1 = 0;
- }
- } else {
- x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
- C1 = x & MASK_BINARY_SIG1;
- }
- // check for zeros (possibly from non-canonical values)
- if (C1 == 0x0ull) {
- // x is 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- // q = nr. of decimal digits in x (1 <= q <= 54)
- // determine first the nr. of bits in x
- if (C1 >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1 >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1 >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- q = nr_digits[x_nr_bits - 1].digits;
- if (q == 0) {
- q = nr_digits[x_nr_bits - 1].digits1;
- if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
- q++;
- }
- exp = x_exp - 398; // unbiased exponent
- if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
- // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
- // so x rounded to an integer may or may not fit in an unsigned 64-bit int
- // the cases that do not fit are identified here; the ones that fit
- // fall through and will be handled with other cases further,
- // under '1 <= q + exp <= 20'
- if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1
- // => set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else { // if n > 0 and q + exp = 20
- // if n >= 2^64 then n is too large
- // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65)
- // <=> C * 10^(21-q) >= 0xa0000000000000000, 1<=q<=16
- if (q == 1) {
- // C * 10^20 >= 0xa0000000000000000
- __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C
- if (C.w[1] >= 0x0a) {
- // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
- // Note: C * 10^(21-q) has 20 or 21 digits; 0xa0000000000000000
- // has 21 digits
- __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]);
- if (C.w[1] >= 0x0a) {
- // actually C.w[1] == 0x0a && C.w[0] >= 0x0000000000000000ull) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- }
- }
- }
- // n is not too large to be converted to int64 if -1 < n < 2^64
- // Note: some of the cases tested for above fall through to this point
- if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
- // set inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- // return 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
- // x <= -1 or 1 <= x < 2^64 so if positive x can be rounded
- // to nearest to a 64-bit unsigned signed integer
- if (x_sign) { // x <= -1
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // 1 <= x < 2^64 so x can be rounded
- // to nearest to a 64-bit unsigned integer
- if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
- ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
- // chop off ind digits from the lower part of C1
- // C1 fits in 64 bits
- // calculate C* and f*
- // C* is actually floor(C*) in this case
- // C* and f* need shifting and masking, as shown by
- // shiftright128[] and maskhigh128[]
- // 1 <= x <= 15
- // kx = 10^(-x) = ten2mk64[ind - 1]
- // C* = C1 * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 54 bits
- __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
- Cstar = P128.w[1];
- fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
- fstar.w[0] = P128.w[0];
- // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
- // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
- // C* = floor(C*) (logical right shift; C has p decimal digits,
- // correct by Property 1)
- // n = C* * 10^(e+x)
- // shift right C* by Ex-64 = shiftright128[ind]
- shift = shiftright128[ind - 1]; // 0 <= shift <= 39
- Cstar = Cstar >> shift;
- // determine inexactness of the rounding of C*
- // if (0 < f* < 10^(-x)) then
- // the result is exact
- // else // if (f* > T*) then
- // the result is inexact
- if (ind - 1 <= 2) { // fstar.w[1] is 0
- if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
- // ten2mk128trunc[ind -1].w[1] is identical to
- // ten2mk128[ind -1].w[1]
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- } // else the result is exact
- } else { // if 3 <= ind - 1 <= 14
- if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
- // ten2mk128trunc[ind -1].w[1] is identical to
- // ten2mk128[ind -1].w[1]
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- } // else the result is exact
- }
- res = Cstar; // the result is positive
- } else if (exp == 0) {
- // 1 <= q <= 10
- // res = +C (exact)
- res = C1; // the result is positive
- } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
- // res = +C * 10^exp (exact)
- res = C1 * ten2k64[exp]; // the result is positive
- }
- }
- BID_RETURN (res);
- }
- /*****************************************************************************
- * BID64_to_uint64_rninta
- ****************************************************************************/
- #if DECIMAL_CALL_BY_REFERENCE
- void
- bid64_to_uint64_rninta (UINT64 * pres, UINT64 * px
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- UINT64 x = *px;
- #else
- UINT64
- bid64_to_uint64_rninta (UINT64 x
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- #endif
- UINT64 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1 represents x_significand (UINT64)
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT64 C1;
- UINT128 C;
- UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
- UINT128 P128;
- // check for NaN or Infinity
- if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // unpack x
- x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
- if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
- x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
- C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
- if (C1 > 9999999999999999ull) { // non-canonical
- x_exp = 0;
- C1 = 0;
- }
- } else {
- x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
- C1 = x & MASK_BINARY_SIG1;
- }
- // check for zeros (possibly from non-canonical values)
- if (C1 == 0x0ull) {
- // x is 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- // q = nr. of decimal digits in x (1 <= q <= 54)
- // determine first the nr. of bits in x
- if (C1 >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1 >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1 >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- q = nr_digits[x_nr_bits - 1].digits;
- if (q == 0) {
- q = nr_digits[x_nr_bits - 1].digits1;
- if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
- q++;
- }
- exp = x_exp - 398; // unbiased exponent
- if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
- // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
- // so x rounded to an integer may or may not fit in an unsigned 64-bit int
- // the cases that do not fit are identified here; the ones that fit
- // fall through and will be handled with other cases further,
- // under '1 <= q + exp <= 20'
- if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1/2
- // => set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else { // if n > 0 and q + exp = 20
- // if n >= 2^64 - 1/2 then n is too large
- // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
- // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=16
- if (q == 1) {
- // C * 10^20 >= 0x9fffffffffffffffb
- __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C
- if (C.w[1] > 0x09 ||
- (C.w[1] == 0x09 && C.w[0] >= 0xfffffffffffffffbull)) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
- // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffffb
- // has 21 digits
- __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]);
- if (C.w[1] > 0x09 ||
- (C.w[1] == 0x09 && C.w[0] >= 0xfffffffffffffffbull)) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- }
- }
- }
- // n is not too large to be converted to int64 if -1/2 <= n < 2^64 - 1/2
- // Note: some of the cases tested for above fall through to this point
- if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
- // return 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
- // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
- // res = 0
- // else if x > 0
- // res = +1
- // else // if x < 0
- // invalid exc
- ind = q - 1; // 0 <= ind <= 15
- if (C1 < midpoint64[ind]) {
- res = 0x0000000000000000ull; // return 0
- } else if (!x_sign) { // n > 0
- res = 0x0000000000000001ull; // return +1
- } else { // if n < 0
- res = 0x8000000000000000ull;
- *pfpsf |= INVALID_EXCEPTION;
- BID_RETURN (res);
- }
- } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
- // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
- // to nearest to a 64-bit unsigned signed integer
- if (x_sign) { // x <= -1
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // 1 <= x < 2^64-1/2 so x can be rounded
- // to nearest to a 64-bit unsigned integer
- if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
- ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
- // chop off ind digits from the lower part of C1
- // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
- C1 = C1 + midpoint64[ind - 1];
- // calculate C* and f*
- // C* is actually floor(C*) in this case
- // C* and f* need shifting and masking, as shown by
- // shiftright128[] and maskhigh128[]
- // 1 <= x <= 15
- // kx = 10^(-x) = ten2mk64[ind - 1]
- // C* = (C1 + 1/2 * 10^x) * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 54 bits
- __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
- Cstar = P128.w[1];
- // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
- // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
- // if (0 < f* < 10^(-x)) then the result is a midpoint
- // if floor(C*) is even then C* = floor(C*) - logical right
- // shift; C* has p decimal digits, correct by Prop. 1)
- // else if floor(C*) is odd C* = floor(C*)-1 (logical right
- // shift; C* has p decimal digits, correct by Pr. 1)
- // else
- // C* = floor(C*) (logical right shift; C has p decimal digits,
- // correct by Property 1)
- // n = C* * 10^(e+x)
- // shift right C* by Ex-64 = shiftright128[ind]
- shift = shiftright128[ind - 1]; // 0 <= shift <= 39
- Cstar = Cstar >> shift;
- // if the result was a midpoint it was rounded away from zero
- res = Cstar; // the result is positive
- } else if (exp == 0) {
- // 1 <= q <= 10
- // res = +C (exact)
- res = C1; // the result is positive
- } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
- // res = +C * 10^exp (exact)
- res = C1 * ten2k64[exp]; // the result is positive
- }
- }
- BID_RETURN (res);
- }
- /*****************************************************************************
- * BID64_to_uint64_xrninta
- ****************************************************************************/
- #if DECIMAL_CALL_BY_REFERENCE
- void
- bid64_to_uint64_xrninta (UINT64 * pres, UINT64 * px
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- UINT64 x = *px;
- #else
- UINT64
- bid64_to_uint64_xrninta (UINT64 x
- _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
- _EXC_INFO_PARAM) {
- #endif
- UINT64 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1 represents x_significand (UINT64)
- UINT64 tmp64;
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT64 C1;
- UINT128 C;
- UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
- UINT128 fstar;
- UINT128 P128;
- // check for NaN or Infinity
- if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // unpack x
- x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
- if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
- x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
- C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
- if (C1 > 9999999999999999ull) { // non-canonical
- x_exp = 0;
- C1 = 0;
- }
- } else {
- x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
- C1 = x & MASK_BINARY_SIG1;
- }
- // check for zeros (possibly from non-canonical values)
- if (C1 == 0x0ull) {
- // x is 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- // q = nr. of decimal digits in x (1 <= q <= 54)
- // determine first the nr. of bits in x
- if (C1 >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1 >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1 >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- q = nr_digits[x_nr_bits - 1].digits;
- if (q == 0) {
- q = nr_digits[x_nr_bits - 1].digits1;
- if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
- q++;
- }
- exp = x_exp - 398; // unbiased exponent
- if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits)
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1)
- // in this case 2^63.11... ~= 10^19 <= x < 10^20 ~= 2^66.43...
- // so x rounded to an integer may or may not fit in an unsigned 64-bit int
- // the cases that do not fit are identified here; the ones that fit
- // fall through and will be handled with other cases further,
- // under '1 <= q + exp <= 20'
- if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1/2
- // => set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- } else { // if n > 0 and q + exp = 20
- // if n >= 2^64 - 1/2 then n is too large
- // <=> c(0)c(1)...c(19).c(20)...c(q-1) >= 2^64-1/2
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^20 >= 2^64-1/2
- // <=> 0.c(0)c(1)...c(19)c(20)...c(q-1) * 10^21 >= 5*(2^65-1)
- // <=> C * 10^(21-q) >= 0x9fffffffffffffffb, 1<=q<=16
- if (q == 1) {
- // C * 10^20 >= 0x9fffffffffffffffb
- __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C
- if (C.w[1] > 0x09 ||
- (C.w[1] == 0x09 && C.w[0] >= 0xfffffffffffffffbull)) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19
- // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffffb
- // has 21 digits
- __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]);
- if (C.w[1] > 0x09 ||
- (C.w[1] == 0x09 && C.w[0] >= 0xfffffffffffffffbull)) {
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // else cases that can be rounded to a 64-bit int fall through
- // to '1 <= q + exp <= 20'
- }
- }
- }
- // n is not too large to be converted to int64 if -1/2 <= n < 2^64 - 1/2
- // Note: some of the cases tested for above fall through to this point
- if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
- // set inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- // return 0
- res = 0x0000000000000000ull;
- BID_RETURN (res);
- } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
- // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
- // res = 0
- // else if x > 0
- // res = +1
- // else // if x < 0
- // invalid exc
- ind = q - 1; // 0 <= ind <= 15
- if (C1 < midpoint64[ind]) {
- res = 0x0000000000000000ull; // return 0
- } else if (!x_sign) { // n > 0
- res = 0x0000000000000001ull; // return +1
- } else { // if n < 0
- res = 0x8000000000000000ull;
- *pfpsf |= INVALID_EXCEPTION;
- BID_RETURN (res);
- }
- // set inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19)
- // x <= -1 or 1 <= x < 2^64-1/2 so if positive x can be rounded
- // to nearest to a 64-bit unsigned signed integer
- if (x_sign) { // x <= -1
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return Integer Indefinite
- res = 0x8000000000000000ull;
- BID_RETURN (res);
- }
- // 1 <= x < 2^64-1/2 so x can be rounded
- // to nearest to a 64-bit unsigned integer
- if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 20
- ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
- // chop off ind digits from the lower part of C1
- // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
- C1 = C1 + midpoint64[ind - 1];
- // calculate C* and f*
- // C* is actually floor(C*) in this case
- // C* and f* need shifting and masking, as shown by
- // shiftright128[] and maskhigh128[]
- // 1 <= x <= 15
- // kx = 10^(-x) = ten2mk64[ind - 1]
- // C* = (C1 + 1/2 * 10^x) * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 54 bits
- __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
- Cstar = P128.w[1];
- fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
- fstar.w[0] = P128.w[0];
- // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
- // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
- // if (0 < f* < 10^(-x)) then the result is a midpoint
- // if floor(C*) is even then C* = floor(C*) - logical right
- // shift; C* has p decimal digits, correct by Prop. 1)
- // else if floor(C*) is odd C* = floor(C*)-1 (logical right
- // shift; C* has p decimal digits, correct by Pr. 1)
- // else
- // C* = floor(C*) (logical right shift; C has p decimal digits,
- // correct by Property 1)
- // n = C* * 10^(e+x)
- // shift right C* by Ex-64 = shiftright128[ind]
- shift = shiftright128[ind - 1]; // 0 <= shift <= 39
- Cstar = Cstar >> shift;
- // determine inexactness of the rounding of C*
- // if (0 < f* - 1/2 < 10^(-x)) then
- // the result is exact
- // else // if (f* - 1/2 > T*) then
- // the result is inexact
- if (ind - 1 <= 2) { // fstar.w[1] is 0
- if (fstar.w[0] > 0x8000000000000000ull) {
- // f* > 1/2 and the result may be exact
- tmp64 = fstar.w[0] - 0x8000000000000000ull; // f* - 1/2
- if ((tmp64 > ten2mk128trunc[ind - 1].w[1])) {
- // ten2mk128trunc[ind -1].w[1] is identical to
- // ten2mk128[ind -1].w[1]
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- } // else the result is exact
- } else { // the result is inexact; f2* <= 1/2
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- }
- } else { // if 3 <= ind - 1 <= 14
- if (fstar.w[1] > onehalf128[ind - 1] ||
- (fstar.w[1] == onehalf128[ind - 1] && fstar.w[0])) {
- // f2* > 1/2 and the result may be exact
- // Calculate f2* - 1/2
- tmp64 = fstar.w[1] - onehalf128[ind - 1];
- if (tmp64 || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
- // ten2mk128trunc[ind -1].w[1] is identical to
- // ten2mk128[ind -1].w[1]
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- } // else the result is exact
- } else { // the result is inexact; f2* <= 1/2
- // set the inexact flag
- *pfpsf |= INEXACT_EXCEPTION;
- }
- }
- // if the result was a midpoint it was rounded away from zero
- res = Cstar; // the result is positive
- } else if (exp == 0) {
- // 1 <= q <= 10
- // res = +C (exact)
- res = C1; // the result is positive
- } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
- // res = +C * 10^exp (exact)
- res = C1 * ten2k64[exp]; // the result is positive
- }
- }
- BID_RETURN (res);
- }
|