1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951 |
- /* 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/>. */
- #define BID_128RES
- #include "bid_internal.h"
- /*****************************************************************************
- * BID128_round_integral_exact
- ****************************************************************************/
- BID128_FUNCTION_ARG1 (bid128_round_integral_exact, x)
- UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
- };
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
- UINT64 tmp64;
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT128 C1;
- UINT256 fstar;
- UINT256 P256;
- // check for NaN or Infinity
- if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
- // x is special
- if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
- // if x = NaN, then res = Q (x)
- // check first for non-canonical NaN payload
- if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
- (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
- (x.w[0] > 0x38c15b09ffffffffull))) {
- x.w[1] = x.w[1] & 0xffffc00000000000ull;
- x.w[0] = 0x0ull;
- }
- if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return quiet (x)
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
- res.w[0] = x.w[0];
- } else { // x is QNaN
- // return x
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
- res.w[0] = x.w[0];
- }
- BID_RETURN (res)
- } else { // x is not a NaN, so it must be infinity
- if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
- // return +inf
- res.w[1] = 0x7800000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- } else { // x is -inf
- // return -inf
- res.w[1] = 0xf800000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- }
- BID_RETURN (res);
- }
- }
- // unpack x
- x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- C1.w[1] = x.w[1] & MASK_COEFF;
- C1.w[0] = x.w[0];
- // check for non-canonical values (treated as zero)
- if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
- // non-canonical
- x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
- C1.w[1] = 0; // significand high
- C1.w[0] = 0; // significand low
- } else { // G0_G1 != 11
- x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
- if (C1.w[1] > 0x0001ed09bead87c0ull ||
- (C1.w[1] == 0x0001ed09bead87c0ull
- && C1.w[0] > 0x378d8e63ffffffffull)) {
- // x is non-canonical if coefficient is larger than 10^34 -1
- C1.w[1] = 0;
- C1.w[0] = 0;
- } else { // canonical
- ;
- }
- }
- // test for input equal to zero
- if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
- // x is 0
- // return 0 preserving the sign bit and the preferred exponent
- // of MAX(Q(x), 0)
- if (x_exp <= (0x1820ull << 49)) {
- res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
- } else {
- res.w[1] = x_sign | x_exp;
- }
- res.w[0] = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- switch (rnd_mode) {
- case ROUNDING_TO_NEAREST:
- case ROUNDING_TIES_AWAY:
- // if (exp <= -(p+1)) return 0.0
- if (x_exp <= 0x2ffa000000000000ull) { // 0x2ffa000000000000ull == -35
- res.w[1] = x_sign | 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- *pfpsf |= INEXACT_EXCEPTION;
- BID_RETURN (res);
- }
- break;
- case ROUNDING_DOWN:
- // if (exp <= -p) return -1.0 or +0.0
- if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffa000000000000ull == -34
- if (x_sign) {
- // if negative, return negative 1, because we know coefficient
- // is non-zero (would have been caught above)
- res.w[1] = 0xb040000000000000ull;
- res.w[0] = 0x0000000000000001ull;
- } else {
- // if positive, return positive 0, because we know coefficient is
- // non-zero (would have been caught above)
- res.w[1] = 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- }
- *pfpsf |= INEXACT_EXCEPTION;
- BID_RETURN (res);
- }
- break;
- case ROUNDING_UP:
- // if (exp <= -p) return -0.0 or +1.0
- if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34
- if (x_sign) {
- // if negative, return negative 0, because we know the coefficient
- // is non-zero (would have been caught above)
- res.w[1] = 0xb040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- } else {
- // if positive, return positive 1, because we know coefficient is
- // non-zero (would have been caught above)
- res.w[1] = 0x3040000000000000ull;
- res.w[0] = 0x0000000000000001ull;
- }
- *pfpsf |= INEXACT_EXCEPTION;
- BID_RETURN (res);
- }
- break;
- case ROUNDING_TO_ZERO:
- // if (exp <= -p) return -0.0 or +0.0
- if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34
- res.w[1] = x_sign | 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- *pfpsf |= INEXACT_EXCEPTION;
- BID_RETURN (res);
- }
- break;
- }
- // q = nr. of decimal digits in x
- // determine first the nr. of bits in x
- if (C1.w[1] == 0) {
- if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) (C1.w[0]); // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1.w[0]; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
- tmp1.d = (double) C1.w[1]; // exact conversion
- x_nr_bits =
- 65 + ((((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.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
- (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
- C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
- q++;
- }
- exp = (x_exp >> 49) - 6176;
- if (exp >= 0) { // -exp <= 0
- // the argument is an integer already
- res.w[1] = x.w[1];
- res.w[0] = x.w[0];
- BID_RETURN (res);
- }
- // exp < 0
- switch (rnd_mode) {
- case ROUNDING_TO_NEAREST:
- if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
- // need to shift right -exp digits from the coefficient; exp will be 0
- ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
- // chop off ind digits from the lower part of C1
- // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
- tmp64 = C1.w[0];
- if (ind <= 19) {
- C1.w[0] = C1.w[0] + midpoint64[ind - 1];
- } else {
- C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
- C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
- }
- if (C1.w[0] < tmp64)
- C1.w[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 <= 34
- // kx = 10^(-x) = ten2mk128[ind - 1]
- // C* = (C1 + 1/2 * 10^x) * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 118 bits
- __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
- // determine the value of res and fstar
- // 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
- // Note: we are going to use ten2mk128[] instead of ten2mk128trunc[]
- if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
- // redundant shift = shiftright128[ind - 1]; // shift = 0
- res.w[1] = P256.w[3];
- res.w[0] = P256.w[2];
- // redundant fstar.w[3] = 0;
- // redundant fstar.w[2] = 0;
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* < 10^(-x) <=> midpoint
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- // if 0 < fstar < 10^(-x), subtract 1 if odd (for rounding to even)
- if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
- ((fstar.w[1] < (ten2mk128[ind - 1].w[1]))
- || ((fstar.w[1] == ten2mk128[ind - 1].w[1])
- && (fstar.w[0] < ten2mk128[ind - 1].w[0])))) {
- // subract 1 to make even
- if (res.w[0]-- == 0) {
- res.w[1]--;
- }
- }
- if (fstar.w[1] > 0x8000000000000000ull ||
- (fstar.w[1] == 0x8000000000000000ull
- && fstar.w[0] > 0x0ull)) {
- // f* > 1/2 and the result may be exact
- tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
- if (tmp64 > ten2mk128[ind - 1].w[1] ||
- (tmp64 == ten2mk128[ind - 1].w[1] &&
- fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- // 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 (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
- shift = shiftright128[ind - 1]; // 3 <= shift <= 63
- res.w[1] = (P256.w[3] >> shift);
- res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
- // redundant fstar.w[3] = 0;
- fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* < 10^(-x) <=> midpoint
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
- fstar.w[2] == 0 && (fstar.w[1] < ten2mk128[ind - 1].w[1] ||
- (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
- fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
- // subract 1 to make even
- if (res.w[0]-- == 0) {
- res.w[1]--;
- }
- }
- if (fstar.w[2] > onehalf128[ind - 1] ||
- (fstar.w[2] == onehalf128[ind - 1]
- && (fstar.w[1] || fstar.w[0]))) {
- // f2* > 1/2 and the result may be exact
- // Calculate f2* - 1/2
- tmp64 = fstar.w[2] - onehalf128[ind - 1];
- if (tmp64 || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
- (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
- fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- // 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 { // 22 <= ind - 1 <= 33
- shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
- res.w[1] = 0;
- res.w[0] = P256.w[3] >> shift;
- fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
- fstar.w[2] = P256.w[2];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* < 10^(-x) <=> midpoint
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
- fstar.w[3] == 0 && fstar.w[2] == 0 &&
- (fstar.w[1] < ten2mk128[ind - 1].w[1] ||
- (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
- fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
- // subract 1 to make even
- if (res.w[0]-- == 0) {
- res.w[1]--;
- }
- }
- if (fstar.w[3] > onehalf128[ind - 1] ||
- (fstar.w[3] == onehalf128[ind - 1] &&
- (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
- // f2* > 1/2 and the result may be exact
- // Calculate f2* - 1/2
- tmp64 = fstar.w[3] - onehalf128[ind - 1];
- if (tmp64 || fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1]
- || (fstar.w[1] == ten2mk128[ind - 1].w[1]
- && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- // 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;
- }
- }
- res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
- BID_RETURN (res);
- } else { // if ((q + exp) < 0) <=> q < -exp
- // the result is +0 or -0
- res.w[1] = x_sign | 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- *pfpsf |= INEXACT_EXCEPTION;
- BID_RETURN (res);
- }
- break;
- case ROUNDING_TIES_AWAY:
- if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
- // need to shift right -exp digits from the coefficient; exp will be 0
- ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
- // chop off ind digits from the lower part of C1
- // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
- tmp64 = C1.w[0];
- if (ind <= 19) {
- C1.w[0] = C1.w[0] + midpoint64[ind - 1];
- } else {
- C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
- C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
- }
- if (C1.w[0] < tmp64)
- C1.w[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 <= 34
- // kx = 10^(-x) = ten2mk128[ind - 1]
- // C* = (C1 + 1/2 * 10^x) * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 118 bits
- __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
- // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
- // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
- // 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)
- // determine also the 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
- // Note: we are going to use ten2mk128[] instead of ten2mk128trunc[]
- // shift right C* by Ex-128 = shiftright128[ind]
- if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
- // redundant shift = shiftright128[ind - 1]; // shift = 0
- res.w[1] = P256.w[3];
- res.w[0] = P256.w[2];
- // redundant fstar.w[3] = 0;
- // redundant fstar.w[2] = 0;
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- if (fstar.w[1] > 0x8000000000000000ull ||
- (fstar.w[1] == 0x8000000000000000ull
- && fstar.w[0] > 0x0ull)) {
- // f* > 1/2 and the result may be exact
- tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
- if ((tmp64 > ten2mk128[ind - 1].w[1] ||
- (tmp64 == ten2mk128[ind - 1].w[1] &&
- fstar.w[0] >= ten2mk128[ind - 1].w[0]))) {
- // 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 (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
- shift = shiftright128[ind - 1]; // 3 <= shift <= 63
- res.w[1] = (P256.w[3] >> shift);
- res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
- // redundant fstar.w[3] = 0;
- fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- if (fstar.w[2] > onehalf128[ind - 1] ||
- (fstar.w[2] == onehalf128[ind - 1]
- && (fstar.w[1] || fstar.w[0]))) {
- // f2* > 1/2 and the result may be exact
- // Calculate f2* - 1/2
- tmp64 = fstar.w[2] - onehalf128[ind - 1];
- if (tmp64 || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
- (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
- fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- // 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 { // 22 <= ind - 1 <= 33
- shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
- res.w[1] = 0;
- res.w[0] = P256.w[3] >> shift;
- fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
- fstar.w[2] = P256.w[2];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- if (fstar.w[3] > onehalf128[ind - 1] ||
- (fstar.w[3] == onehalf128[ind - 1] &&
- (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
- // f2* > 1/2 and the result may be exact
- // Calculate f2* - 1/2
- tmp64 = fstar.w[3] - onehalf128[ind - 1];
- if (tmp64 || fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1]
- || (fstar.w[1] == ten2mk128[ind - 1].w[1]
- && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- // 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 already rounded away from zero
- res.w[1] |= x_sign | 0x3040000000000000ull;
- BID_RETURN (res);
- } else { // if ((q + exp) < 0) <=> q < -exp
- // the result is +0 or -0
- res.w[1] = x_sign | 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- *pfpsf |= INEXACT_EXCEPTION;
- BID_RETURN (res);
- }
- break;
- case ROUNDING_DOWN:
- if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
- // need to shift right -exp digits from the coefficient; exp will be 0
- ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
- // (number of digits to be chopped off)
- // chop off ind digits from the lower part of C1
- // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
- // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
- // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
- // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
- // tmp64 = C1.w[0];
- // if (ind <= 19) {
- // C1.w[0] = C1.w[0] + midpoint64[ind - 1];
- // } else {
- // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
- // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
- // }
- // if (C1.w[0] < tmp64) C1.w[1]++;
- // if carry-out from C1.w[0], increment C1.w[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 <= 34
- // kx = 10^(-x) = ten2mk128[ind - 1]
- // C* = (C1 + 1/2 * 10^x) * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 118 bits
- __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
- if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
- res.w[1] = P256.w[3];
- res.w[0] = P256.w[2];
- // redundant fstar.w[3] = 0;
- // redundant fstar.w[2] = 0;
- // redundant fstar.w[1] = P256.w[1];
- // redundant fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if ((P256.w[1] > ten2mk128[ind - 1].w[1])
- || (P256.w[1] == ten2mk128[ind - 1].w[1]
- && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
- *pfpsf |= INEXACT_EXCEPTION;
- // if positive, the truncated value is already the correct result
- if (x_sign) { // if negative
- if (++res.w[0] == 0) {
- res.w[1]++;
- }
- }
- }
- } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
- shift = shiftright128[ind - 1]; // 0 <= shift <= 102
- res.w[1] = (P256.w[3] >> shift);
- res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
- // redundant fstar.w[3] = 0;
- fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
- (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
- fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- *pfpsf |= INEXACT_EXCEPTION;
- // if positive, the truncated value is already the correct result
- if (x_sign) { // if negative
- if (++res.w[0] == 0) {
- res.w[1]++;
- }
- }
- }
- } else { // 22 <= ind - 1 <= 33
- shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
- res.w[1] = 0;
- res.w[0] = P256.w[3] >> shift;
- fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
- fstar.w[2] = P256.w[2];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if (fstar.w[3] || fstar.w[2]
- || fstar.w[1] > ten2mk128[ind - 1].w[1]
- || (fstar.w[1] == ten2mk128[ind - 1].w[1]
- && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- *pfpsf |= INEXACT_EXCEPTION;
- // if positive, the truncated value is already the correct result
- if (x_sign) { // if negative
- if (++res.w[0] == 0) {
- res.w[1]++;
- }
- }
- }
- }
- res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
- BID_RETURN (res);
- } else { // if exp < 0 and q + exp <= 0
- if (x_sign) { // negative rounds down to -1.0
- res.w[1] = 0xb040000000000000ull;
- res.w[0] = 0x0000000000000001ull;
- } else { // positive rpunds down to +0.0
- res.w[1] = 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- }
- *pfpsf |= INEXACT_EXCEPTION;
- BID_RETURN (res);
- }
- break;
- case ROUNDING_UP:
- if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
- // need to shift right -exp digits from the coefficient; exp will be 0
- ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
- // (number of digits to be chopped off)
- // chop off ind digits from the lower part of C1
- // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
- // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
- // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
- // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
- // tmp64 = C1.w[0];
- // if (ind <= 19) {
- // C1.w[0] = C1.w[0] + midpoint64[ind - 1];
- // } else {
- // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
- // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
- // }
- // if (C1.w[0] < tmp64) C1.w[1]++;
- // if carry-out from C1.w[0], increment C1.w[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 <= 34
- // kx = 10^(-x) = ten2mk128[ind - 1]
- // C* = C1 * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 118 bits
- __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
- if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
- res.w[1] = P256.w[3];
- res.w[0] = P256.w[2];
- // redundant fstar.w[3] = 0;
- // redundant fstar.w[2] = 0;
- // redundant fstar.w[1] = P256.w[1];
- // redundant fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if ((P256.w[1] > ten2mk128[ind - 1].w[1])
- || (P256.w[1] == ten2mk128[ind - 1].w[1]
- && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
- *pfpsf |= INEXACT_EXCEPTION;
- // if negative, the truncated value is already the correct result
- if (!x_sign) { // if positive
- if (++res.w[0] == 0) {
- res.w[1]++;
- }
- }
- }
- } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
- shift = shiftright128[ind - 1]; // 3 <= shift <= 63
- res.w[1] = (P256.w[3] >> shift);
- res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
- // redundant fstar.w[3] = 0;
- fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
- (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
- fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- *pfpsf |= INEXACT_EXCEPTION;
- // if negative, the truncated value is already the correct result
- if (!x_sign) { // if positive
- if (++res.w[0] == 0) {
- res.w[1]++;
- }
- }
- }
- } else { // 22 <= ind - 1 <= 33
- shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
- res.w[1] = 0;
- res.w[0] = P256.w[3] >> shift;
- fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
- fstar.w[2] = P256.w[2];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if (fstar.w[3] || fstar.w[2]
- || fstar.w[1] > ten2mk128[ind - 1].w[1]
- || (fstar.w[1] == ten2mk128[ind - 1].w[1]
- && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- *pfpsf |= INEXACT_EXCEPTION;
- // if negative, the truncated value is already the correct result
- if (!x_sign) { // if positive
- if (++res.w[0] == 0) {
- res.w[1]++;
- }
- }
- }
- }
- res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
- BID_RETURN (res);
- } else { // if exp < 0 and q + exp <= 0
- if (x_sign) { // negative rounds up to -0.0
- res.w[1] = 0xb040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- } else { // positive rpunds up to +1.0
- res.w[1] = 0x3040000000000000ull;
- res.w[0] = 0x0000000000000001ull;
- }
- *pfpsf |= INEXACT_EXCEPTION;
- BID_RETURN (res);
- }
- break;
- case ROUNDING_TO_ZERO:
- if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
- // need to shift right -exp digits from the coefficient; exp will be 0
- ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
- // (number of digits to be chopped off)
- // chop off ind digits from the lower part of C1
- // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
- // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
- // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
- // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
- //tmp64 = C1.w[0];
- // if (ind <= 19) {
- // C1.w[0] = C1.w[0] + midpoint64[ind - 1];
- // } else {
- // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
- // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
- // }
- // if (C1.w[0] < tmp64) C1.w[1]++;
- // if carry-out from C1.w[0], increment C1.w[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 <= 34
- // kx = 10^(-x) = ten2mk128[ind - 1]
- // C* = (C1 + 1/2 * 10^x) * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 118 bits
- __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
- if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
- res.w[1] = P256.w[3];
- res.w[0] = P256.w[2];
- // redundant fstar.w[3] = 0;
- // redundant fstar.w[2] = 0;
- // redundant fstar.w[1] = P256.w[1];
- // redundant fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if ((P256.w[1] > ten2mk128[ind - 1].w[1])
- || (P256.w[1] == ten2mk128[ind - 1].w[1]
- && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
- *pfpsf |= INEXACT_EXCEPTION;
- }
- } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
- shift = shiftright128[ind - 1]; // 3 <= shift <= 63
- res.w[1] = (P256.w[3] >> shift);
- res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
- // redundant fstar.w[3] = 0;
- fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
- (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
- fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- *pfpsf |= INEXACT_EXCEPTION;
- }
- } else { // 22 <= ind - 1 <= 33
- shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
- res.w[1] = 0;
- res.w[0] = P256.w[3] >> shift;
- fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
- fstar.w[2] = P256.w[2];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if (fstar.w[3] || fstar.w[2]
- || fstar.w[1] > ten2mk128[ind - 1].w[1]
- || (fstar.w[1] == ten2mk128[ind - 1].w[1]
- && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- *pfpsf |= INEXACT_EXCEPTION;
- }
- }
- res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
- BID_RETURN (res);
- } else { // if exp < 0 and q + exp <= 0 the result is +0 or -0
- res.w[1] = x_sign | 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- *pfpsf |= INEXACT_EXCEPTION;
- BID_RETURN (res);
- }
- break;
- }
- BID_RETURN (res);
- }
- /*****************************************************************************
- * BID128_round_integral_nearest_even
- ****************************************************************************/
- BID128_FUNCTION_ARG1_NORND (bid128_round_integral_nearest_even, x)
- UINT128 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
- UINT64 tmp64;
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT128 C1;
- // UINT128 res is C* at first - represents up to 34 decimal digits ~ 113 bits
- UINT256 fstar;
- UINT256 P256;
- // check for NaN or Infinity
- if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
- // x is special
- if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
- // if x = NaN, then res = Q (x)
- // check first for non-canonical NaN payload
- if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
- (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
- (x.w[0] > 0x38c15b09ffffffffull))) {
- x.w[1] = x.w[1] & 0xffffc00000000000ull;
- x.w[0] = 0x0ull;
- }
- if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return quiet (x)
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
- res.w[0] = x.w[0];
- } else { // x is QNaN
- // return x
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
- res.w[0] = x.w[0];
- }
- BID_RETURN (res)
- } else { // x is not a NaN, so it must be infinity
- if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
- // return +inf
- res.w[1] = 0x7800000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- } else { // x is -inf
- // return -inf
- res.w[1] = 0xf800000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- }
- BID_RETURN (res);
- }
- }
- // unpack x
- x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- C1.w[1] = x.w[1] & MASK_COEFF;
- C1.w[0] = x.w[0];
- // check for non-canonical values (treated as zero)
- if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
- // non-canonical
- x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
- C1.w[1] = 0; // significand high
- C1.w[0] = 0; // significand low
- } else { // G0_G1 != 11
- x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
- if (C1.w[1] > 0x0001ed09bead87c0ull ||
- (C1.w[1] == 0x0001ed09bead87c0ull
- && C1.w[0] > 0x378d8e63ffffffffull)) {
- // x is non-canonical if coefficient is larger than 10^34 -1
- C1.w[1] = 0;
- C1.w[0] = 0;
- } else { // canonical
- ;
- }
- }
- // test for input equal to zero
- if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
- // x is 0
- // return 0 preserving the sign bit and the preferred exponent
- // of MAX(Q(x), 0)
- if (x_exp <= (0x1820ull << 49)) {
- res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
- } else {
- res.w[1] = x_sign | x_exp;
- }
- res.w[0] = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- // if (exp <= -(p+1)) return 0
- if (x_exp <= 0x2ffa000000000000ull) { // 0x2ffa000000000000ull == -35
- res.w[1] = x_sign | 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // q = nr. of decimal digits in x
- // determine first the nr. of bits in x
- if (C1.w[1] == 0) {
- if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) (C1.w[0]); // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1.w[0]; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
- tmp1.d = (double) C1.w[1]; // exact conversion
- x_nr_bits =
- 65 + ((((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.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
- || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
- C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
- q++;
- }
- exp = (x_exp >> 49) - 6176;
- if (exp >= 0) { // -exp <= 0
- // the argument is an integer already
- res.w[1] = x.w[1];
- res.w[0] = x.w[0];
- BID_RETURN (res);
- } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
- // need to shift right -exp digits from the coefficient; the exp will be 0
- ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
- // chop off ind digits from the lower part of C1
- // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
- tmp64 = C1.w[0];
- if (ind <= 19) {
- C1.w[0] = C1.w[0] + midpoint64[ind - 1];
- } else {
- C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
- C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
- }
- if (C1.w[0] < tmp64)
- C1.w[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 <= 34
- // kx = 10^(-x) = ten2mk128[ind - 1]
- // C* = (C1 + 1/2 * 10^x) * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 118 bits
- __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
- // determine the value of res and fstar
- if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
- // redundant shift = shiftright128[ind - 1]; // shift = 0
- res.w[1] = P256.w[3];
- res.w[0] = P256.w[2];
- // redundant fstar.w[3] = 0;
- // redundant fstar.w[2] = 0;
- // redundant fstar.w[1] = P256.w[1];
- // redundant fstar.w[0] = P256.w[0];
- // fraction f* < 10^(-x) <=> midpoint
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- // if 0 < fstar < 10^(-x), subtract 1 if odd (for rounding to even)
- if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
- ((P256.w[1] < (ten2mk128[ind - 1].w[1]))
- || ((P256.w[1] == ten2mk128[ind - 1].w[1])
- && (P256.w[0] < ten2mk128[ind - 1].w[0])))) {
- // subract 1 to make even
- if (res.w[0]-- == 0) {
- res.w[1]--;
- }
- }
- } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
- shift = shiftright128[ind - 1]; // 3 <= shift <= 63
- res.w[1] = (P256.w[3] >> shift);
- res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
- // redundant fstar.w[3] = 0;
- fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* < 10^(-x) <=> midpoint
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
- fstar.w[2] == 0 && (fstar.w[1] < ten2mk128[ind - 1].w[1] ||
- (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
- fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
- // subract 1 to make even
- if (res.w[0]-- == 0) {
- res.w[1]--;
- }
- }
- } else { // 22 <= ind - 1 <= 33
- shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
- res.w[1] = 0;
- res.w[0] = P256.w[3] >> shift;
- fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
- fstar.w[2] = P256.w[2];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* < 10^(-x) <=> midpoint
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
- fstar.w[3] == 0 && fstar.w[2] == 0
- && (fstar.w[1] < ten2mk128[ind - 1].w[1]
- || (fstar.w[1] == ten2mk128[ind - 1].w[1]
- && fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
- // subract 1 to make even
- if (res.w[0]-- == 0) {
- res.w[1]--;
- }
- }
- }
- res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
- BID_RETURN (res);
- } else { // if ((q + exp) < 0) <=> q < -exp
- // the result is +0 or -0
- res.w[1] = x_sign | 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- }
- /*****************************************************************************
- * BID128_round_integral_negative
- ****************************************************************************/
- BID128_FUNCTION_ARG1_NORND (bid128_round_integral_negative, x)
- UINT128 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo
- // (all are UINT64)
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT128 C1;
- // UINT128 res is C* at first - represents up to 34 decimal digits ~
- // 113 bits
- UINT256 fstar;
- UINT256 P256;
- // check for NaN or Infinity
- if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
- // x is special
- if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
- // if x = NaN, then res = Q (x)
- // check first for non-canonical NaN payload
- if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
- (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
- (x.w[0] > 0x38c15b09ffffffffull))) {
- x.w[1] = x.w[1] & 0xffffc00000000000ull;
- x.w[0] = 0x0ull;
- }
- if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return quiet (x)
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
- res.w[0] = x.w[0];
- } else { // x is QNaN
- // return x
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
- res.w[0] = x.w[0];
- }
- BID_RETURN (res)
- } else { // x is not a NaN, so it must be infinity
- if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
- // return +inf
- res.w[1] = 0x7800000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- } else { // x is -inf
- // return -inf
- res.w[1] = 0xf800000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- }
- BID_RETURN (res);
- }
- }
- // unpack x
- x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- C1.w[1] = x.w[1] & MASK_COEFF;
- C1.w[0] = x.w[0];
- // check for non-canonical values (treated as zero)
- if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
- // non-canonical
- x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
- C1.w[1] = 0; // significand high
- C1.w[0] = 0; // significand low
- } else { // G0_G1 != 11
- x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
- if (C1.w[1] > 0x0001ed09bead87c0ull ||
- (C1.w[1] == 0x0001ed09bead87c0ull
- && C1.w[0] > 0x378d8e63ffffffffull)) {
- // x is non-canonical if coefficient is larger than 10^34 -1
- C1.w[1] = 0;
- C1.w[0] = 0;
- } else { // canonical
- ;
- }
- }
- // test for input equal to zero
- if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
- // x is 0
- // return 0 preserving the sign bit and the preferred exponent
- // of MAX(Q(x), 0)
- if (x_exp <= (0x1820ull << 49)) {
- res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
- } else {
- res.w[1] = x_sign | x_exp;
- }
- res.w[0] = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- // if (exp <= -p) return -1.0 or +0.0
- if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34
- if (x_sign) {
- // if negative, return negative 1, because we know the coefficient
- // is non-zero (would have been caught above)
- res.w[1] = 0xb040000000000000ull;
- res.w[0] = 0x0000000000000001ull;
- } else {
- // if positive, return positive 0, because we know coefficient is
- // non-zero (would have been caught above)
- res.w[1] = 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- }
- BID_RETURN (res);
- }
- // q = nr. of decimal digits in x
- // determine first the nr. of bits in x
- if (C1.w[1] == 0) {
- if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) (C1.w[0]); // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1.w[0]; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
- tmp1.d = (double) C1.w[1]; // exact conversion
- x_nr_bits =
- 65 + ((((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.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
- (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
- C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
- q++;
- }
- exp = (x_exp >> 49) - 6176;
- if (exp >= 0) { // -exp <= 0
- // the argument is an integer already
- res.w[1] = x.w[1];
- res.w[0] = x.w[0];
- BID_RETURN (res);
- } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
- // need to shift right -exp digits from the coefficient; the exp will be 0
- ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
- // (number of digits to be chopped off)
- // chop off ind digits from the lower part of C1
- // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
- // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
- // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
- // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
- //tmp64 = C1.w[0];
- // if (ind <= 19) {
- // C1.w[0] = C1.w[0] + midpoint64[ind - 1];
- // } else {
- // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
- // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
- // }
- // if (C1.w[0] < tmp64) C1.w[1]++;
- // if carry-out from C1.w[0], increment C1.w[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 <= 34
- // kx = 10^(-x) = ten2mk128[ind - 1]
- // C* = (C1 + 1/2 * 10^x) * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 118 bits
- __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
- if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
- res.w[1] = P256.w[3];
- res.w[0] = P256.w[2];
- // if positive, the truncated value is already the correct result
- if (x_sign) { // if negative
- // redundant fstar.w[3] = 0;
- // redundant fstar.w[2] = 0;
- // redundant fstar.w[1] = P256.w[1];
- // redundant fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if ((P256.w[1] > ten2mk128[ind - 1].w[1])
- || (P256.w[1] == ten2mk128[ind - 1].w[1]
- && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
- if (++res.w[0] == 0) {
- res.w[1]++;
- }
- }
- }
- } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
- shift = shiftright128[ind - 1]; // 0 <= shift <= 102
- res.w[1] = (P256.w[3] >> shift);
- res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
- // if positive, the truncated value is already the correct result
- if (x_sign) { // if negative
- // redundant fstar.w[3] = 0;
- fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
- (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
- fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- if (++res.w[0] == 0) {
- res.w[1]++;
- }
- }
- }
- } else { // 22 <= ind - 1 <= 33
- shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
- res.w[1] = 0;
- res.w[0] = P256.w[3] >> shift;
- // if positive, the truncated value is already the correct result
- if (x_sign) { // if negative
- fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
- fstar.w[2] = P256.w[2];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if (fstar.w[3] || fstar.w[2]
- || fstar.w[1] > ten2mk128[ind - 1].w[1]
- || (fstar.w[1] == ten2mk128[ind - 1].w[1]
- && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- if (++res.w[0] == 0) {
- res.w[1]++;
- }
- }
- }
- }
- res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
- BID_RETURN (res);
- } else { // if exp < 0 and q + exp <= 0
- if (x_sign) { // negative rounds down to -1.0
- res.w[1] = 0xb040000000000000ull;
- res.w[0] = 0x0000000000000001ull;
- } else { // positive rpunds down to +0.0
- res.w[1] = 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- }
- BID_RETURN (res);
- }
- }
- /*****************************************************************************
- * BID128_round_integral_positive
- ****************************************************************************/
- BID128_FUNCTION_ARG1_NORND (bid128_round_integral_positive, x)
- UINT128 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo
- // (all are UINT64)
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT128 C1;
- // UINT128 res is C* at first - represents up to 34 decimal digits ~
- // 113 bits
- UINT256 fstar;
- UINT256 P256;
- // check for NaN or Infinity
- if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
- // x is special
- if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
- // if x = NaN, then res = Q (x)
- // check first for non-canonical NaN payload
- if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
- (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
- (x.w[0] > 0x38c15b09ffffffffull))) {
- x.w[1] = x.w[1] & 0xffffc00000000000ull;
- x.w[0] = 0x0ull;
- }
- if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return quiet (x)
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
- res.w[0] = x.w[0];
- } else { // x is QNaN
- // return x
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
- res.w[0] = x.w[0];
- }
- BID_RETURN (res)
- } else { // x is not a NaN, so it must be infinity
- if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
- // return +inf
- res.w[1] = 0x7800000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- } else { // x is -inf
- // return -inf
- res.w[1] = 0xf800000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- }
- BID_RETURN (res);
- }
- }
- // unpack x
- x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- C1.w[1] = x.w[1] & MASK_COEFF;
- C1.w[0] = x.w[0];
- // check for non-canonical values (treated as zero)
- if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
- // non-canonical
- x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
- C1.w[1] = 0; // significand high
- C1.w[0] = 0; // significand low
- } else { // G0_G1 != 11
- x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
- if (C1.w[1] > 0x0001ed09bead87c0ull ||
- (C1.w[1] == 0x0001ed09bead87c0ull
- && C1.w[0] > 0x378d8e63ffffffffull)) {
- // x is non-canonical if coefficient is larger than 10^34 -1
- C1.w[1] = 0;
- C1.w[0] = 0;
- } else { // canonical
- ;
- }
- }
- // test for input equal to zero
- if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
- // x is 0
- // return 0 preserving the sign bit and the preferred exponent
- // of MAX(Q(x), 0)
- if (x_exp <= (0x1820ull << 49)) {
- res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
- } else {
- res.w[1] = x_sign | x_exp;
- }
- res.w[0] = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- // if (exp <= -p) return -0.0 or +1.0
- if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34
- if (x_sign) {
- // if negative, return negative 0, because we know the coefficient
- // is non-zero (would have been caught above)
- res.w[1] = 0xb040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- } else {
- // if positive, return positive 1, because we know coefficient is
- // non-zero (would have been caught above)
- res.w[1] = 0x3040000000000000ull;
- res.w[0] = 0x0000000000000001ull;
- }
- BID_RETURN (res);
- }
- // q = nr. of decimal digits in x
- // determine first the nr. of bits in x
- if (C1.w[1] == 0) {
- if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
- // split 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) (C1.w[0]); // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1.w[0]; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
- tmp1.d = (double) C1.w[1]; // exact conversion
- x_nr_bits =
- 65 + ((((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.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
- (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
- C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
- q++;
- }
- exp = (x_exp >> 49) - 6176;
- if (exp >= 0) { // -exp <= 0
- // the argument is an integer already
- res.w[1] = x.w[1];
- res.w[0] = x.w[0];
- BID_RETURN (res);
- } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
- // need to shift right -exp digits from the coefficient; exp will be 0
- ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
- // (number of digits to be chopped off)
- // chop off ind digits from the lower part of C1
- // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
- // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
- // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
- // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
- // tmp64 = C1.w[0];
- // if (ind <= 19) {
- // C1.w[0] = C1.w[0] + midpoint64[ind - 1];
- // } else {
- // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
- // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
- // }
- // if (C1.w[0] < tmp64) C1.w[1]++;
- // if carry-out from C1.w[0], increment C1.w[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 <= 34
- // kx = 10^(-x) = ten2mk128[ind - 1]
- // C* = C1 * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 118 bits
- __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
- if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
- res.w[1] = P256.w[3];
- res.w[0] = P256.w[2];
- // if negative, the truncated value is already the correct result
- if (!x_sign) { // if positive
- // redundant fstar.w[3] = 0;
- // redundant fstar.w[2] = 0;
- // redundant fstar.w[1] = P256.w[1];
- // redundant fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if ((P256.w[1] > ten2mk128[ind - 1].w[1])
- || (P256.w[1] == ten2mk128[ind - 1].w[1]
- && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
- if (++res.w[0] == 0) {
- res.w[1]++;
- }
- }
- }
- } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
- shift = shiftright128[ind - 1]; // 3 <= shift <= 63
- res.w[1] = (P256.w[3] >> shift);
- res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
- // if negative, the truncated value is already the correct result
- if (!x_sign) { // if positive
- // redundant fstar.w[3] = 0;
- fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
- (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
- fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- if (++res.w[0] == 0) {
- res.w[1]++;
- }
- }
- }
- } else { // 22 <= ind - 1 <= 33
- shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
- res.w[1] = 0;
- res.w[0] = P256.w[3] >> shift;
- // if negative, the truncated value is already the correct result
- if (!x_sign) { // if positive
- fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
- fstar.w[2] = P256.w[2];
- fstar.w[1] = P256.w[1];
- fstar.w[0] = P256.w[0];
- // fraction f* > 10^(-x) <=> inexact
- // f* is in the right position to be compared with
- // 10^(-x) from ten2mk128[]
- if (fstar.w[3] || fstar.w[2]
- || fstar.w[1] > ten2mk128[ind - 1].w[1]
- || (fstar.w[1] == ten2mk128[ind - 1].w[1]
- && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
- if (++res.w[0] == 0) {
- res.w[1]++;
- }
- }
- }
- }
- res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
- BID_RETURN (res);
- } else { // if exp < 0 and q + exp <= 0
- if (x_sign) { // negative rounds up to -0.0
- res.w[1] = 0xb040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- } else { // positive rpunds up to +1.0
- res.w[1] = 0x3040000000000000ull;
- res.w[0] = 0x0000000000000001ull;
- }
- BID_RETURN (res);
- }
- }
- /*****************************************************************************
- * BID128_round_integral_zero
- ****************************************************************************/
- BID128_FUNCTION_ARG1_NORND (bid128_round_integral_zero, x)
- UINT128 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo
- // (all are UINT64)
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT128 C1;
- // UINT128 res is C* at first - represents up to 34 decimal digits ~
- // 113 bits
- UINT256 P256;
- // check for NaN or Infinity
- if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
- // x is special
- if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
- // if x = NaN, then res = Q (x)
- // check first for non-canonical NaN payload
- if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
- (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
- (x.w[0] > 0x38c15b09ffffffffull))) {
- x.w[1] = x.w[1] & 0xffffc00000000000ull;
- x.w[0] = 0x0ull;
- }
- if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return quiet (x)
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
- res.w[0] = x.w[0];
- } else { // x is QNaN
- // return x
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
- res.w[0] = x.w[0];
- }
- BID_RETURN (res)
- } else { // x is not a NaN, so it must be infinity
- if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
- // return +inf
- res.w[1] = 0x7800000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- } else { // x is -inf
- // return -inf
- res.w[1] = 0xf800000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- }
- BID_RETURN (res);
- }
- }
- // unpack x
- x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- C1.w[1] = x.w[1] & MASK_COEFF;
- C1.w[0] = x.w[0];
- // check for non-canonical values (treated as zero)
- if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
- // non-canonical
- x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
- C1.w[1] = 0; // significand high
- C1.w[0] = 0; // significand low
- } else { // G0_G1 != 11
- x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
- if (C1.w[1] > 0x0001ed09bead87c0ull ||
- (C1.w[1] == 0x0001ed09bead87c0ull
- && C1.w[0] > 0x378d8e63ffffffffull)) {
- // x is non-canonical if coefficient is larger than 10^34 -1
- C1.w[1] = 0;
- C1.w[0] = 0;
- } else { // canonical
- ;
- }
- }
- // test for input equal to zero
- if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
- // x is 0
- // return 0 preserving the sign bit and the preferred exponent
- // of MAX(Q(x), 0)
- if (x_exp <= (0x1820ull << 49)) {
- res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
- } else {
- res.w[1] = x_sign | x_exp;
- }
- res.w[0] = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- // if (exp <= -p) return -0.0 or +0.0
- if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34
- res.w[1] = x_sign | 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // q = nr. of decimal digits in x
- // determine first the nr. of bits in x
- if (C1.w[1] == 0) {
- if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) (C1.w[0]); // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1.w[0]; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
- tmp1.d = (double) C1.w[1]; // exact conversion
- x_nr_bits =
- 65 + ((((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.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
- (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
- C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
- q++;
- }
- exp = (x_exp >> 49) - 6176;
- if (exp >= 0) { // -exp <= 0
- // the argument is an integer already
- res.w[1] = x.w[1];
- res.w[0] = x.w[0];
- BID_RETURN (res);
- } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
- // need to shift right -exp digits from the coefficient; the exp will be 0
- ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
- // (number of digits to be chopped off)
- // chop off ind digits from the lower part of C1
- // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
- // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
- // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
- // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
- //tmp64 = C1.w[0];
- // if (ind <= 19) {
- // C1.w[0] = C1.w[0] + midpoint64[ind - 1];
- // } else {
- // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
- // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
- // }
- // if (C1.w[0] < tmp64) C1.w[1]++;
- // if carry-out from C1.w[0], increment C1.w[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 <= 34
- // kx = 10^(-x) = ten2mk128[ind - 1]
- // C* = (C1 + 1/2 * 10^x) * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 118 bits
- __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
- if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
- res.w[1] = P256.w[3];
- res.w[0] = P256.w[2];
- } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
- shift = shiftright128[ind - 1]; // 3 <= shift <= 63
- res.w[1] = (P256.w[3] >> shift);
- res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
- } else { // 22 <= ind - 1 <= 33
- shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
- res.w[1] = 0;
- res.w[0] = P256.w[3] >> shift;
- }
- res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
- BID_RETURN (res);
- } else { // if exp < 0 and q + exp <= 0 the result is +0 or -0
- res.w[1] = x_sign | 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- }
- /*****************************************************************************
- * BID128_round_integral_nearest_away
- ****************************************************************************/
- BID128_FUNCTION_ARG1_NORND (bid128_round_integral_nearest_away, x)
- UINT128 res;
- UINT64 x_sign;
- UINT64 x_exp;
- int exp; // unbiased exponent
- // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo
- // (all are UINT64)
- UINT64 tmp64;
- BID_UI64DOUBLE tmp1;
- unsigned int x_nr_bits;
- int q, ind, shift;
- UINT128 C1;
- // UINT128 res is C* at first - represents up to 34 decimal digits ~
- // 113 bits
- // UINT256 fstar;
- UINT256 P256;
- // check for NaN or Infinity
- if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
- // x is special
- if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
- // if x = NaN, then res = Q (x)
- // check first for non-canonical NaN payload
- if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
- (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
- (x.w[0] > 0x38c15b09ffffffffull))) {
- x.w[1] = x.w[1] & 0xffffc00000000000ull;
- x.w[0] = 0x0ull;
- }
- if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
- // set invalid flag
- *pfpsf |= INVALID_EXCEPTION;
- // return quiet (x)
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
- res.w[0] = x.w[0];
- } else { // x is QNaN
- // return x
- res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
- res.w[0] = x.w[0];
- }
- BID_RETURN (res)
- } else { // x is not a NaN, so it must be infinity
- if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
- // return +inf
- res.w[1] = 0x7800000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- } else { // x is -inf
- // return -inf
- res.w[1] = 0xf800000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- }
- BID_RETURN (res);
- }
- }
- // unpack x
- x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
- C1.w[1] = x.w[1] & MASK_COEFF;
- C1.w[0] = x.w[0];
- // check for non-canonical values (treated as zero)
- if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
- // non-canonical
- x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
- C1.w[1] = 0; // significand high
- C1.w[0] = 0; // significand low
- } else { // G0_G1 != 11
- x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
- if (C1.w[1] > 0x0001ed09bead87c0ull ||
- (C1.w[1] == 0x0001ed09bead87c0ull
- && C1.w[0] > 0x378d8e63ffffffffull)) {
- // x is non-canonical if coefficient is larger than 10^34 -1
- C1.w[1] = 0;
- C1.w[0] = 0;
- } else { // canonical
- ;
- }
- }
- // test for input equal to zero
- if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
- // x is 0
- // return 0 preserving the sign bit and the preferred exponent
- // of MAX(Q(x), 0)
- if (x_exp <= (0x1820ull << 49)) {
- res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
- } else {
- res.w[1] = x_sign | x_exp;
- }
- res.w[0] = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // x is not special and is not zero
- // if (exp <= -(p+1)) return 0.0
- if (x_exp <= 0x2ffa000000000000ull) { // 0x2ffa000000000000ull == -35
- res.w[1] = x_sign | 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- // q = nr. of decimal digits in x
- // determine first the nr. of bits in x
- if (C1.w[1] == 0) {
- if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
- // split the 64-bit value in two 32-bit halves to avoid rounding errors
- if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
- tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
- x_nr_bits =
- 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- } else { // x < 2^32
- tmp1.d = (double) (C1.w[0]); // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // if x < 2^53
- tmp1.d = (double) C1.w[0]; // exact conversion
- x_nr_bits =
- 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
- }
- } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
- tmp1.d = (double) C1.w[1]; // exact conversion
- x_nr_bits =
- 65 + ((((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.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
- (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
- C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
- q++;
- }
- exp = (x_exp >> 49) - 6176;
- if (exp >= 0) { // -exp <= 0
- // the argument is an integer already
- res.w[1] = x.w[1];
- res.w[0] = x.w[0];
- BID_RETURN (res);
- } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
- // need to shift right -exp digits from the coefficient; the exp will be 0
- ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
- // chop off ind digits from the lower part of C1
- // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
- tmp64 = C1.w[0];
- if (ind <= 19) {
- C1.w[0] = C1.w[0] + midpoint64[ind - 1];
- } else {
- C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
- C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
- }
- if (C1.w[0] < tmp64)
- C1.w[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 <= 34
- // kx = 10^(-x) = ten2mk128[ind - 1]
- // C* = (C1 + 1/2 * 10^x) * 10^(-x)
- // the approximation of 10^(-x) was rounded up to 118 bits
- __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
- // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
- // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
- // 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-128 = shiftright128[ind]
- if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
- res.w[1] = P256.w[3];
- res.w[0] = P256.w[2];
- } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
- shift = shiftright128[ind - 1]; // 3 <= shift <= 63
- res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
- res.w[1] = (P256.w[3] >> shift);
- } else { // 22 <= ind - 1 <= 33
- shift = shiftright128[ind - 1]; // 2 <= shift <= 38
- res.w[1] = 0;
- res.w[0] = (P256.w[3] >> (shift - 64)); // 2 <= shift - 64 <= 38
- }
- // if the result was a midpoint, it was already rounded away from zero
- res.w[1] |= x_sign | 0x3040000000000000ull;
- BID_RETURN (res);
- } else { // if ((q + exp) < 0) <=> q < -exp
- // the result is +0 or -0
- res.w[1] = x_sign | 0x3040000000000000ull;
- res.w[0] = 0x0000000000000000ull;
- BID_RETURN (res);
- }
- }
|