From nobody Tue Feb 10 08:31:26 2026 Delivered-To: importer@patchew.org Authentication-Results: mx.zohomail.com; dkim=pass; spf=pass (zohomail.com: domain of gnu.org designates 209.51.188.17 as permitted sender) smtp.mailfrom=qemu-devel-bounces+importer=patchew.org@nongnu.org; dmarc=pass(p=none dis=none) header.from=linaro.org ARC-Seal: i=1; a=rsa-sha256; t=1622757113; cv=none; d=zohomail.com; s=zohoarc; b=TGKVo8f/sRGn2tMZBYLIgeTNfB6Qq5lZ11LAFe6YG7zswVeRJBhAJuNFkdscysW+dTJ+Kv4YZiULE1XWaFuvcApWoxMDoFIgBWa+EO72UM01EG0GQbZ33OgBZkdtHNM789qpbEmu0NiQtbzRDSiOMdqJ+LaAuhGSC5y8V5PvnFo= ARC-Message-Signature: i=1; a=rsa-sha256; c=relaxed/relaxed; d=zohomail.com; s=zohoarc; t=1622757113; h=Content-Type:Content-Transfer-Encoding:Cc:Date:From:In-Reply-To:List-Subscribe:List-Post:List-Id:List-Archive:List-Help:List-Unsubscribe:MIME-Version:Message-ID:References:Sender:Subject:To; bh=PYPaf/j9AHT9fcW+1JhpwV4TzbljbGsjiwi3XwcqOTc=; b=QJPNTKkP0dlLhXo32pLYpsLYOzycy5FwqORD/FG2i7WBJYsX9LKMiFi1njRTWqJI/ua0moNrHeOTHHL2ACZkY3lgVwaaRpwu3ytRn+ZbAk4iHzmur0EHIEgyGqud0t+M51MZMCzbgut/yjm7iNSi1c9+rSiDIvrP1PK2RKw7410= ARC-Authentication-Results: i=1; mx.zohomail.com; dkim=pass; spf=pass (zohomail.com: domain of gnu.org designates 209.51.188.17 as permitted sender) smtp.mailfrom=qemu-devel-bounces+importer=patchew.org@nongnu.org; dmarc=pass header.from= (p=none dis=none) header.from= Return-Path: Received: from lists.gnu.org (lists.gnu.org [209.51.188.17]) by mx.zohomail.com with SMTPS id 1622757113172332.534079771944; Thu, 3 Jun 2021 14:51:53 -0700 (PDT) Received: from localhost ([::1]:59048 helo=lists1p.gnu.org) by lists.gnu.org with esmtp (Exim 4.90_1) (envelope-from ) id 1lovFn-00050T-UA for importer@patchew.org; Thu, 03 Jun 2021 17:51:51 -0400 Received: from eggs.gnu.org ([2001:470:142:3::10]:58582) by lists.gnu.org with esmtps (TLS1.2:ECDHE_RSA_AES_256_GCM_SHA384:256) (Exim 4.90_1) (envelope-from ) id 1lov6L-0003EO-5h for qemu-devel@nongnu.org; Thu, 03 Jun 2021 17:42:05 -0400 Received: from mail-pj1-x1033.google.com ([2607:f8b0:4864:20::1033]:44845) by eggs.gnu.org with esmtps (TLS1.2:ECDHE_RSA_AES_128_GCM_SHA256:128) (Exim 4.90_1) (envelope-from ) id 1lov68-0001BF-Se for qemu-devel@nongnu.org; Thu, 03 Jun 2021 17:42:04 -0400 Received: by mail-pj1-x1033.google.com with SMTP id h12-20020a17090aa88cb029016400fd8ad8so4683039pjq.3 for ; Thu, 03 Jun 2021 14:41:52 -0700 (PDT) Received: from localhost.localdomain (174-21-70-228.tukw.qwest.net. [174.21.70.228]) by smtp.gmail.com with ESMTPSA id p65sm40115pfb.62.2021.06.03.14.41.49 (version=TLS1_3 cipher=TLS_AES_256_GCM_SHA384 bits=256/256); Thu, 03 Jun 2021 14:41:50 -0700 (PDT) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=linaro.org; s=google; h=from:to:cc:subject:date:message-id:in-reply-to:references :mime-version:content-transfer-encoding; bh=PYPaf/j9AHT9fcW+1JhpwV4TzbljbGsjiwi3XwcqOTc=; b=ewFf7mWYNjU2Dbd5H05g+3/Wzxj6M6YmXqmCSDu15ym+wtOJPMvZs4ltepZFo2PKAM o2HxZJK6c09NSX6UXvDj8maJpjfA7PK0dhTwGFhuPgYwyPL/n1BSJvApWQJ6krF8/nCy 5N/ap885CSxv0kDPAQmHY6Ek1AIZ8jA8cl2qIYufubQpWc/5HH8mXBgPp3E3LMQp7S/K zqaBdwF4wv3e8rOgD+Xu5g6tnKQSTo93v9vNPQPP4EEsdKBFq2VCL8XoE+fjUUiJ8HV0 pElK5W2xabDxL9fBmXSGSrWp0OgdfbEHmiyS1T4clxqQqHu1q11FTk2EgRlsD7uciEQX ESwQ== X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20161025; h=x-gm-message-state:from:to:cc:subject:date:message-id:in-reply-to :references:mime-version:content-transfer-encoding; bh=PYPaf/j9AHT9fcW+1JhpwV4TzbljbGsjiwi3XwcqOTc=; b=DKN3EX0YfY4bgH62206+H4yKbaoG26Iy8IRyy4Q7pMM4ZRKTd7fsUdmVZR6+LOWW8i dl8fGu8gmy7Lxqy0NkKOIQHKJ9UJrwImxop2ClOOU2WTILbQzSsFtzu+NxJ3xogCFxfo Q1wr26fFzYPZ9wH/d2S92FUThzn1rwOWR72B4ezbKgUtIHvSaGALAcakbTaxogp1CWPc mOJsH5eUeLDteZnFzG4SIkYU2nI4/j0p1hkuc9rxfgqDMkJd4Vl6u1oFXAG9q0/WVc2V CmrTIVtEgO8YTi9IH2Kwde7buRVNPG4CVVh/b2COaEeJ842a1utH/3y7qBnc7sOwG/Um sM1g== X-Gm-Message-State: AOAM533qqkPKMacYOlhrPyCFROiK/jCXNAj1IqmFhO942QUxl+R8JKHR 7pkC2J0X+uU8h2tGGGep+BhlZQiswOoSLw== X-Google-Smtp-Source: ABdhPJxVP4Dk532zGdBPEtW4PNOPHEOkpmjEW80i/dz1E67khmFOdL98cGnRBUbcdnRu98Ke+xaPaA== X-Received: by 2002:a17:90a:8d82:: with SMTP id d2mr13359307pjo.106.1622756510446; Thu, 03 Jun 2021 14:41:50 -0700 (PDT) From: Richard Henderson To: qemu-devel@nongnu.org Subject: [PULL 27/29] softfloat: Convert modrem operations to FloatParts Date: Thu, 3 Jun 2021 14:41:29 -0700 Message-Id: <20210603214131.629841-28-richard.henderson@linaro.org> X-Mailer: git-send-email 2.25.1 In-Reply-To: <20210603214131.629841-1-richard.henderson@linaro.org> References: <20210603214131.629841-1-richard.henderson@linaro.org> MIME-Version: 1.0 Content-Type: text/plain; charset="utf-8" Content-Transfer-Encoding: quoted-printable Received-SPF: pass (zohomail.com: domain of gnu.org designates 209.51.188.17 as permitted sender) client-ip=209.51.188.17; envelope-from=qemu-devel-bounces+importer=patchew.org@nongnu.org; helo=lists.gnu.org; Received-SPF: pass client-ip=2607:f8b0:4864:20::1033; envelope-from=richard.henderson@linaro.org; helo=mail-pj1-x1033.google.com X-Spam_score_int: -20 X-Spam_score: -2.1 X-Spam_bar: -- X-Spam_report: (-2.1 / 5.0 requ) BAYES_00=-1.9, DKIM_SIGNED=0.1, DKIM_VALID=-0.1, DKIM_VALID_AU=-0.1, DKIM_VALID_EF=-0.1, RCVD_IN_DNSWL_NONE=-0.0001, SPF_HELO_NONE=0.001, SPF_PASS=-0.001 autolearn=ham autolearn_force=no X-Spam_action: no action X-BeenThere: qemu-devel@nongnu.org X-Mailman-Version: 2.1.23 Precedence: list List-Id: List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Cc: =?UTF-8?q?Alex=20Benn=C3=A9e?= Errors-To: qemu-devel-bounces+importer=patchew.org@nongnu.org Sender: "Qemu-devel" X-ZohoMail-DKIM: pass (identity @linaro.org) Rename to parts$N_modrem. This was the last use of a lot of the legacy infrastructure, so remove it as required. Reviewed-by: Alex Benn=C3=A9e Signed-off-by: Richard Henderson --- include/fpu/softfloat-macros.h | 34 + fpu/softfloat.c | 1339 +++++++------------------------- fpu/softfloat-parts.c.inc | 34 + fpu/softfloat-specialize.c.inc | 165 ---- 4 files changed, 329 insertions(+), 1243 deletions(-) diff --git a/include/fpu/softfloat-macros.h b/include/fpu/softfloat-macros.h index ec4e27a595..81c3fe8256 100644 --- a/include/fpu/softfloat-macros.h +++ b/include/fpu/softfloat-macros.h @@ -745,4 +745,38 @@ static inline bool ne128(uint64_t a0, uint64_t a1, uin= t64_t b0, uint64_t b1) return a0 !=3D b0 || a1 !=3D b1; } =20 +/* + * Similarly, comparisons of 192-bit values. + */ + +static inline bool eq192(uint64_t a0, uint64_t a1, uint64_t a2, + uint64_t b0, uint64_t b1, uint64_t b2) +{ + return ((a0 ^ b0) | (a1 ^ b1) | (a2 ^ b2)) =3D=3D 0; +} + +static inline bool le192(uint64_t a0, uint64_t a1, uint64_t a2, + uint64_t b0, uint64_t b1, uint64_t b2) +{ + if (a0 !=3D b0) { + return a0 < b0; + } + if (a1 !=3D b1) { + return a1 < b1; + } + return a2 <=3D b2; +} + +static inline bool lt192(uint64_t a0, uint64_t a1, uint64_t a2, + uint64_t b0, uint64_t b1, uint64_t b2) +{ + if (a0 !=3D b0) { + return a0 < b0; + } + if (a1 !=3D b1) { + return a1 < b1; + } + return a2 < b2; +} + #endif diff --git a/fpu/softfloat.c b/fpu/softfloat.c index c0fe191f4d..5026f518b0 100644 --- a/fpu/softfloat.c +++ b/fpu/softfloat.c @@ -401,60 +401,6 @@ float64_gen2(float64 xa, float64 xb, float_status *s, return soft(ua.s, ub.s, s); } =20 -/*------------------------------------------------------------------------= ---- -| Returns the fraction bits of the single-precision floating-point value `= a'. -*-------------------------------------------------------------------------= ---*/ - -static inline uint32_t extractFloat32Frac(float32 a) -{ - return float32_val(a) & 0x007FFFFF; -} - -/*------------------------------------------------------------------------= ---- -| Returns the exponent bits of the single-precision floating-point value `= a'. -*-------------------------------------------------------------------------= ---*/ - -static inline int extractFloat32Exp(float32 a) -{ - return (float32_val(a) >> 23) & 0xFF; -} - -/*------------------------------------------------------------------------= ---- -| Returns the sign bit of the single-precision floating-point value `a'. -*-------------------------------------------------------------------------= ---*/ - -static inline bool extractFloat32Sign(float32 a) -{ - return float32_val(a) >> 31; -} - -/*------------------------------------------------------------------------= ---- -| Returns the fraction bits of the double-precision floating-point value `= a'. -*-------------------------------------------------------------------------= ---*/ - -static inline uint64_t extractFloat64Frac(float64 a) -{ - return float64_val(a) & UINT64_C(0x000FFFFFFFFFFFFF); -} - -/*------------------------------------------------------------------------= ---- -| Returns the exponent bits of the double-precision floating-point value `= a'. -*-------------------------------------------------------------------------= ---*/ - -static inline int extractFloat64Exp(float64 a) -{ - return (float64_val(a) >> 52) & 0x7FF; -} - -/*------------------------------------------------------------------------= ---- -| Returns the sign bit of the double-precision floating-point value `a'. -*-------------------------------------------------------------------------= ---*/ - -static inline bool extractFloat64Sign(float64 a) -{ - return float64_val(a) >> 63; -} - /* * Classify a floating point number. Everything above float_class_qnan * is a NaN so cls >=3D float_class_qnan is any NaN. @@ -845,6 +791,14 @@ static FloatParts128 *parts128_div(FloatParts128 *a, F= loatParts128 *b, #define parts_div(A, B, S) \ PARTS_GENERIC_64_128(div, A)(A, B, S) =20 +static FloatParts64 *parts64_modrem(FloatParts64 *a, FloatParts64 *b, + uint64_t *mod_quot, float_status *s); +static FloatParts128 *parts128_modrem(FloatParts128 *a, FloatParts128 *b, + uint64_t *mod_quot, float_status *s); + +#define parts_modrem(A, B, Q, S) \ + PARTS_GENERIC_64_128(modrem, A)(A, B, Q, S) + static void parts64_sqrt(FloatParts64 *a, float_status *s, const FloatFmt = *f); static void parts128_sqrt(FloatParts128 *a, float_status *s, const FloatFm= t *f); =20 @@ -1229,6 +1183,186 @@ static int frac256_normalize(FloatParts256 *a) =20 #define frac_normalize(A) FRAC_GENERIC_64_128_256(normalize, A)(A) =20 +static void frac64_modrem(FloatParts64 *a, FloatParts64 *b, uint64_t *mod_= quot) +{ + uint64_t a0, a1, b0, t0, t1, q, quot; + int exp_diff =3D a->exp - b->exp; + int shift; + + a0 =3D a->frac; + a1 =3D 0; + + if (exp_diff < -1) { + if (mod_quot) { + *mod_quot =3D 0; + } + return; + } + if (exp_diff =3D=3D -1) { + a0 >>=3D 1; + exp_diff =3D 0; + } + + b0 =3D b->frac; + quot =3D q =3D b0 <=3D a0; + if (q) { + a0 -=3D b0; + } + + exp_diff -=3D 64; + while (exp_diff > 0) { + q =3D estimateDiv128To64(a0, a1, b0); + q =3D q > 2 ? q - 2 : 0; + mul64To128(b0, q, &t0, &t1); + sub128(a0, a1, t0, t1, &a0, &a1); + shortShift128Left(a0, a1, 62, &a0, &a1); + exp_diff -=3D 62; + quot =3D (quot << 62) + q; + } + + exp_diff +=3D 64; + if (exp_diff > 0) { + q =3D estimateDiv128To64(a0, a1, b0); + q =3D q > 2 ? (q - 2) >> (64 - exp_diff) : 0; + mul64To128(b0, q << (64 - exp_diff), &t0, &t1); + sub128(a0, a1, t0, t1, &a0, &a1); + shortShift128Left(0, b0, 64 - exp_diff, &t0, &t1); + while (le128(t0, t1, a0, a1)) { + ++q; + sub128(a0, a1, t0, t1, &a0, &a1); + } + quot =3D (exp_diff < 64 ? quot << exp_diff : 0) + q; + } else { + t0 =3D b0; + t1 =3D 0; + } + + if (mod_quot) { + *mod_quot =3D quot; + } else { + sub128(t0, t1, a0, a1, &t0, &t1); + if (lt128(t0, t1, a0, a1) || + (eq128(t0, t1, a0, a1) && (q & 1))) { + a0 =3D t0; + a1 =3D t1; + a->sign =3D !a->sign; + } + } + + if (likely(a0)) { + shift =3D clz64(a0); + shortShift128Left(a0, a1, shift, &a0, &a1); + } else if (likely(a1)) { + shift =3D clz64(a1); + a0 =3D a1 << shift; + a1 =3D 0; + shift +=3D 64; + } else { + a->cls =3D float_class_zero; + return; + } + + a->exp =3D b->exp + exp_diff - shift; + a->frac =3D a0 | (a1 !=3D 0); +} + +static void frac128_modrem(FloatParts128 *a, FloatParts128 *b, + uint64_t *mod_quot) +{ + uint64_t a0, a1, a2, b0, b1, t0, t1, t2, q, quot; + int exp_diff =3D a->exp - b->exp; + int shift; + + a0 =3D a->frac_hi; + a1 =3D a->frac_lo; + a2 =3D 0; + + if (exp_diff < -1) { + if (mod_quot) { + *mod_quot =3D 0; + } + return; + } + if (exp_diff =3D=3D -1) { + shift128Right(a0, a1, 1, &a0, &a1); + exp_diff =3D 0; + } + + b0 =3D b->frac_hi; + b1 =3D b->frac_lo; + + quot =3D q =3D le128(b0, b1, a0, a1); + if (q) { + sub128(a0, a1, b0, b1, &a0, &a1); + } + + exp_diff -=3D 64; + while (exp_diff > 0) { + q =3D estimateDiv128To64(a0, a1, b0); + q =3D q > 4 ? q - 4 : 0; + mul128By64To192(b0, b1, q, &t0, &t1, &t2); + sub192(a0, a1, a2, t0, t1, t2, &a0, &a1, &a2); + shortShift192Left(a0, a1, a2, 61, &a0, &a1, &a2); + exp_diff -=3D 61; + quot =3D (quot << 61) + q; + } + + exp_diff +=3D 64; + if (exp_diff > 0) { + q =3D estimateDiv128To64(a0, a1, b0); + q =3D q > 4 ? (q - 4) >> (64 - exp_diff) : 0; + mul128By64To192(b0, b1, q << (64 - exp_diff), &t0, &t1, &t2); + sub192(a0, a1, a2, t0, t1, t2, &a0, &a1, &a2); + shortShift192Left(0, b0, b1, 64 - exp_diff, &t0, &t1, &t2); + while (le192(t0, t1, t2, a0, a1, a2)) { + ++q; + sub192(a0, a1, a2, t0, t1, t2, &a0, &a1, &a2); + } + quot =3D (exp_diff < 64 ? quot << exp_diff : 0) + q; + } else { + t0 =3D b0; + t1 =3D b1; + t2 =3D 0; + } + + if (mod_quot) { + *mod_quot =3D quot; + } else { + sub192(t0, t1, t2, a0, a1, a2, &t0, &t1, &t2); + if (lt192(t0, t1, t2, a0, a1, a2) || + (eq192(t0, t1, t2, a0, a1, a2) && (q & 1))) { + a0 =3D t0; + a1 =3D t1; + a2 =3D t2; + a->sign =3D !a->sign; + } + } + + if (likely(a0)) { + shift =3D clz64(a0); + shortShift192Left(a0, a1, a2, shift, &a0, &a1, &a2); + } else if (likely(a1)) { + shift =3D clz64(a1); + shortShift128Left(a1, a2, shift, &a0, &a1); + a2 =3D 0; + shift +=3D 64; + } else if (likely(a2)) { + shift =3D clz64(a2); + a0 =3D a2 << shift; + a1 =3D a2 =3D 0; + shift +=3D 128; + } else { + a->cls =3D float_class_zero; + return; + } + + a->exp =3D b->exp + exp_diff - shift; + a->frac_hi =3D a0; + a->frac_lo =3D a1 | (a2 !=3D 0); +} + +#define frac_modrem(A, B, Q) FRAC_GENERIC_64_128(modrem, A)(A, B, Q) + static void frac64_shl(FloatParts64 *a, int c) { a->frac <<=3D c; @@ -2313,6 +2447,79 @@ floatx80 floatx80_div(floatx80 a, floatx80 b, float_= status *status) return floatx80_round_pack_canonical(pr, status); } =20 +/* + * Remainder + */ + +float32 float32_rem(float32 a, float32 b, float_status *status) +{ + FloatParts64 pa, pb, *pr; + + float32_unpack_canonical(&pa, a, status); + float32_unpack_canonical(&pb, b, status); + pr =3D parts_modrem(&pa, &pb, NULL, status); + + return float32_round_pack_canonical(pr, status); +} + +float64 float64_rem(float64 a, float64 b, float_status *status) +{ + FloatParts64 pa, pb, *pr; + + float64_unpack_canonical(&pa, a, status); + float64_unpack_canonical(&pb, b, status); + pr =3D parts_modrem(&pa, &pb, NULL, status); + + return float64_round_pack_canonical(pr, status); +} + +float128 float128_rem(float128 a, float128 b, float_status *status) +{ + FloatParts128 pa, pb, *pr; + + float128_unpack_canonical(&pa, a, status); + float128_unpack_canonical(&pb, b, status); + pr =3D parts_modrem(&pa, &pb, NULL, status); + + return float128_round_pack_canonical(pr, status); +} + +/* + * Returns the remainder of the extended double-precision floating-point v= alue + * `a' with respect to the corresponding value `b'. + * If 'mod' is false, the operation is performed according to the IEC/IEEE + * Standard for Binary Floating-Point Arithmetic. If 'mod' is true, return + * the remainder based on truncating the quotient toward zero instead and + * *quotient is set to the low 64 bits of the absolute value of the integer + * quotient. + */ +floatx80 floatx80_modrem(floatx80 a, floatx80 b, bool mod, + uint64_t *quotient, float_status *status) +{ + FloatParts128 pa, pb, *pr; + + *quotient =3D 0; + if (!floatx80_unpack_canonical(&pa, a, status) || + !floatx80_unpack_canonical(&pb, b, status)) { + return floatx80_default_nan(status); + } + pr =3D parts_modrem(&pa, &pb, mod ? quotient : NULL, status); + + return floatx80_round_pack_canonical(pr, status); +} + +floatx80 floatx80_rem(floatx80 a, floatx80 b, float_status *status) +{ + uint64_t quotient; + return floatx80_modrem(a, b, false, "ient, status); +} + +floatx80 floatx80_mod(floatx80 a, floatx80 b, float_status *status) +{ + uint64_t quotient; + return floatx80_modrem(a, b, true, "ient, status); +} + /* * Float to Float conversions * @@ -4262,300 +4469,6 @@ bfloat16 bfloat16_squash_input_denormal(bfloat16 a,= float_status *status) return a; } =20 -/*------------------------------------------------------------------------= ---- -| Normalizes the subnormal single-precision floating-point value represent= ed -| by the denormalized significand `aSig'. The normalized exponent and -| significand are stored at the locations pointed to by `zExpPtr' and -| `zSigPtr', respectively. -*-------------------------------------------------------------------------= ---*/ - -static void - normalizeFloat32Subnormal(uint32_t aSig, int *zExpPtr, uint32_t *zSigPtr) -{ - int8_t shiftCount; - - shiftCount =3D clz32(aSig) - 8; - *zSigPtr =3D aSig<float_rounding_mode; - roundNearestEven =3D ( roundingMode =3D=3D float_round_nearest_even ); - switch (roundingMode) { - case float_round_nearest_even: - case float_round_ties_away: - roundIncrement =3D 0x40; - break; - case float_round_to_zero: - roundIncrement =3D 0; - break; - case float_round_up: - roundIncrement =3D zSign ? 0 : 0x7f; - break; - case float_round_down: - roundIncrement =3D zSign ? 0x7f : 0; - break; - case float_round_to_odd: - roundIncrement =3D zSig & 0x80 ? 0 : 0x7f; - break; - default: - abort(); - break; - } - roundBits =3D zSig & 0x7F; - if ( 0xFD <=3D (uint16_t) zExp ) { - if ( ( 0xFD < zExp ) - || ( ( zExp =3D=3D 0xFD ) - && ( (int32_t) ( zSig + roundIncrement ) < 0 ) ) - ) { - bool overflow_to_inf =3D roundingMode !=3D float_round_to_odd = && - roundIncrement !=3D 0; - float_raise(float_flag_overflow | float_flag_inexact, status); - return packFloat32(zSign, 0xFF, -!overflow_to_inf); - } - if ( zExp < 0 ) { - if (status->flush_to_zero) { - float_raise(float_flag_output_denormal, status); - return packFloat32(zSign, 0, 0); - } - isTiny =3D status->tininess_before_rounding - || (zExp < -1) - || (zSig + roundIncrement < 0x80000000); - shift32RightJamming( zSig, - zExp, &zSig ); - zExp =3D 0; - roundBits =3D zSig & 0x7F; - if (isTiny && roundBits) { - float_raise(float_flag_underflow, status); - } - if (roundingMode =3D=3D float_round_to_odd) { - /* - * For round-to-odd case, the roundIncrement depends on - * zSig which just changed. - */ - roundIncrement =3D zSig & 0x80 ? 0 : 0x7f; - } - } - } - if (roundBits) { - float_raise(float_flag_inexact, status); - } - zSig =3D ( zSig + roundIncrement )>>7; - if (!(roundBits ^ 0x40) && roundNearestEven) { - zSig &=3D ~1; - } - if ( zSig =3D=3D 0 ) zExp =3D 0; - return packFloat32( zSign, zExp, zSig ); - -} - -/*------------------------------------------------------------------------= ---- -| Takes an abstract floating-point value having sign `zSign', exponent `zE= xp', -| and significand `zSig', and returns the proper single-precision floating- -| point value corresponding to the abstract input. This routine is just l= ike -| `roundAndPackFloat32' except that `zSig' does not have to be normalized. -| Bit 31 of `zSig' must be zero, and `zExp' must be 1 less than the ``true= '' -| floating-point exponent. -*-------------------------------------------------------------------------= ---*/ - -static float32 - normalizeRoundAndPackFloat32(bool zSign, int zExp, uint32_t zSig, - float_status *status) -{ - int8_t shiftCount; - - shiftCount =3D clz32(zSig) - 1; - return roundAndPackFloat32(zSign, zExp - shiftCount, zSig<float_rounding_mode; - roundNearestEven =3D ( roundingMode =3D=3D float_round_nearest_even ); - switch (roundingMode) { - case float_round_nearest_even: - case float_round_ties_away: - roundIncrement =3D 0x200; - break; - case float_round_to_zero: - roundIncrement =3D 0; - break; - case float_round_up: - roundIncrement =3D zSign ? 0 : 0x3ff; - break; - case float_round_down: - roundIncrement =3D zSign ? 0x3ff : 0; - break; - case float_round_to_odd: - roundIncrement =3D (zSig & 0x400) ? 0 : 0x3ff; - break; - default: - abort(); - } - roundBits =3D zSig & 0x3FF; - if ( 0x7FD <=3D (uint16_t) zExp ) { - if ( ( 0x7FD < zExp ) - || ( ( zExp =3D=3D 0x7FD ) - && ( (int64_t) ( zSig + roundIncrement ) < 0 ) ) - ) { - bool overflow_to_inf =3D roundingMode !=3D float_round_to_odd = && - roundIncrement !=3D 0; - float_raise(float_flag_overflow | float_flag_inexact, status); - return packFloat64(zSign, 0x7FF, -(!overflow_to_inf)); - } - if ( zExp < 0 ) { - if (status->flush_to_zero) { - float_raise(float_flag_output_denormal, status); - return packFloat64(zSign, 0, 0); - } - isTiny =3D status->tininess_before_rounding - || (zExp < -1) - || (zSig + roundIncrement < UINT64_C(0x8000000000000000)= ); - shift64RightJamming( zSig, - zExp, &zSig ); - zExp =3D 0; - roundBits =3D zSig & 0x3FF; - if (isTiny && roundBits) { - float_raise(float_flag_underflow, status); - } - if (roundingMode =3D=3D float_round_to_odd) { - /* - * For round-to-odd case, the roundIncrement depends on - * zSig which just changed. - */ - roundIncrement =3D (zSig & 0x400) ? 0 : 0x3ff; - } - } - } - if (roundBits) { - float_raise(float_flag_inexact, status); - } - zSig =3D ( zSig + roundIncrement )>>10; - if (!(roundBits ^ 0x200) && roundNearestEven) { - zSig &=3D ~1; - } - if ( zSig =3D=3D 0 ) zExp =3D 0; - return packFloat64( zSign, zExp, zSig ); - -} - -/*------------------------------------------------------------------------= ---- -| Takes an abstract floating-point value having sign `zSign', exponent `zE= xp', -| and significand `zSig', and returns the proper double-precision floating- -| point value corresponding to the abstract input. This routine is just l= ike -| `roundAndPackFloat64' except that `zSig' does not have to be normalized. -| Bit 63 of `zSig' must be zero, and `zExp' must be 1 less than the ``true= '' -| floating-point exponent. -*-------------------------------------------------------------------------= ---*/ - -static float64 - normalizeRoundAndPackFloat64(bool zSign, int zExp, uint64_t zSig, - float_status *status) -{ - int8_t shiftCount; - - shiftCount =3D clz64(zSig) - 1; - return roundAndPackFloat64(zSign, zExp - shiftCount, zSig<>48 ) & 0x7FFF; - -} - -/*------------------------------------------------------------------------= ---- -| Returns the sign bit of the quadruple-precision floating-point value `a'. -*-------------------------------------------------------------------------= ---*/ - -static inline bool extractFloat128Sign(float128 a) -{ - return a.high >> 63; -} - -/*------------------------------------------------------------------------= ---- -| Normalizes the subnormal quadruple-precision floating-point value -| represented by the denormalized significand formed by the concatenation = of -| `aSig0' and `aSig1'. The normalized exponent is stored at the location -| pointed to by `zExpPtr'. The most significant 49 bits of the normalized -| significand are stored at the location pointed to by `zSig0Ptr', and the -| least significant 64 bits of the normalized significand are stored at the -| location pointed to by `zSig1Ptr'. -*-------------------------------------------------------------------------= ---*/ - -static void - normalizeFloat128Subnormal( - uint64_t aSig0, - uint64_t aSig1, - int32_t *zExpPtr, - uint64_t *zSig0Ptr, - uint64_t *zSig1Ptr - ) -{ - int8_t shiftCount; - - if ( aSig0 =3D=3D 0 ) { - shiftCount =3D clz64(aSig1) - 15; - if ( shiftCount < 0 ) { - *zSig0Ptr =3D aSig1>>( - shiftCount ); - *zSig1Ptr =3D aSig1<<( shiftCount & 63 ); - } - else { - *zSig0Ptr =3D aSig1<float_rounding_mode; - roundNearestEven =3D ( roundingMode =3D=3D float_round_nearest_even ); - switch (roundingMode) { - case float_round_nearest_even: - case float_round_ties_away: - increment =3D ((int64_t)zSig2 < 0); - break; - case float_round_to_zero: - increment =3D 0; - break; - case float_round_up: - increment =3D !zSign && zSig2; - break; - case float_round_down: - increment =3D zSign && zSig2; - break; - case float_round_to_odd: - increment =3D !(zSig1 & 0x1) && zSig2; - break; - default: - abort(); - } - if ( 0x7FFD <=3D (uint32_t) zExp ) { - if ( ( 0x7FFD < zExp ) - || ( ( zExp =3D=3D 0x7FFD ) - && eq128( - UINT64_C(0x0001FFFFFFFFFFFF), - UINT64_C(0xFFFFFFFFFFFFFFFF), - zSig0, - zSig1 - ) - && increment - ) - ) { - float_raise(float_flag_overflow | float_flag_inexact, status); - if ( ( roundingMode =3D=3D float_round_to_zero ) - || ( zSign && ( roundingMode =3D=3D float_round_up ) ) - || ( ! zSign && ( roundingMode =3D=3D float_round_down ) ) - || (roundingMode =3D=3D float_round_to_odd) - ) { - return - packFloat128( - zSign, - 0x7FFE, - UINT64_C(0x0000FFFFFFFFFFFF), - UINT64_C(0xFFFFFFFFFFFFFFFF) - ); - } - return packFloat128( zSign, 0x7FFF, 0, 0 ); - } - if ( zExp < 0 ) { - if (status->flush_to_zero) { - float_raise(float_flag_output_denormal, status); - return packFloat128(zSign, 0, 0, 0); - } - isTiny =3D status->tininess_before_rounding - || (zExp < -1) - || !increment - || lt128(zSig0, zSig1, - UINT64_C(0x0001FFFFFFFFFFFF), - UINT64_C(0xFFFFFFFFFFFFFFFF)); - shift128ExtraRightJamming( - zSig0, zSig1, zSig2, - zExp, &zSig0, &zSig1, &zSig2 ); - zExp =3D 0; - if (isTiny && zSig2) { - float_raise(float_flag_underflow, status); - } - switch (roundingMode) { - case float_round_nearest_even: - case float_round_ties_away: - increment =3D ((int64_t)zSig2 < 0); - break; - case float_round_to_zero: - increment =3D 0; - break; - case float_round_up: - increment =3D !zSign && zSig2; - break; - case float_round_down: - increment =3D zSign && zSig2; - break; - case float_round_to_odd: - increment =3D !(zSig1 & 0x1) && zSig2; - break; - default: - abort(); - } - } - } - if (zSig2) { - float_raise(float_flag_inexact, status); - } - if ( increment ) { - add128( zSig0, zSig1, 0, 1, &zSig0, &zSig1 ); - if ((zSig2 + zSig2 =3D=3D 0) && roundNearestEven) { - zSig1 &=3D ~1; - } - } - else { - if ( ( zSig0 | zSig1 ) =3D=3D 0 ) zExp =3D 0; - } - return packFloat128( zSign, zExp, zSig0, zSig1 ); - -} - -/*------------------------------------------------------------------------= ---- -| Takes an abstract floating-point value having sign `zSign', exponent `zE= xp', -| and significand formed by the concatenation of `zSig0' and `zSig1', and -| returns the proper quadruple-precision floating-point value corresponding -| to the abstract input. This routine is just like `roundAndPackFloat128' -| except that the input significand has fewer bits and does not have to be -| normalized. In all cases, `zExp' must be 1 less than the ``true'' float= ing- -| point exponent. -*-------------------------------------------------------------------------= ---*/ - -static float128 normalizeRoundAndPackFloat128(bool zSign, int32_t zExp, - uint64_t zSig0, uint64_t zSi= g1, - float_status *status) -{ - int8_t shiftCount; - uint64_t zSig2; - - if ( zSig0 =3D=3D 0 ) { - zSig0 =3D zSig1; - zSig1 =3D 0; - zExp -=3D 64; - } - shiftCount =3D clz64(zSig0) - 15; - if ( 0 <=3D shiftCount ) { - zSig2 =3D 0; - shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 ); - } - else { - shift128ExtraRightJamming( - zSig0, zSig1, 0, - shiftCount, &zSig0, &zSig1, &zSig2 ); - } - zExp -=3D shiftCount; - return roundAndPackFloat128(zSign, zExp, zSig0, zSig1, zSig2, status); - -} - -/*------------------------------------------------------------------------= ---- -| Returns the remainder of the single-precision floating-point value `a' -| with respect to the corresponding value `b'. The operation is performed -| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. -*-------------------------------------------------------------------------= ---*/ - -float32 float32_rem(float32 a, float32 b, float_status *status) -{ - bool aSign, zSign; - int aExp, bExp, expDiff; - uint32_t aSig, bSig; - uint32_t q; - uint64_t aSig64, bSig64, q64; - uint32_t alternateASig; - int32_t sigMean; - a =3D float32_squash_input_denormal(a, status); - b =3D float32_squash_input_denormal(b, status); - - aSig =3D extractFloat32Frac( a ); - aExp =3D extractFloat32Exp( a ); - aSign =3D extractFloat32Sign( a ); - bSig =3D extractFloat32Frac( b ); - bExp =3D extractFloat32Exp( b ); - if ( aExp =3D=3D 0xFF ) { - if ( aSig || ( ( bExp =3D=3D 0xFF ) && bSig ) ) { - return propagateFloat32NaN(a, b, status); - } - float_raise(float_flag_invalid, status); - return float32_default_nan(status); - } - if ( bExp =3D=3D 0xFF ) { - if (bSig) { - return propagateFloat32NaN(a, b, status); - } - return a; - } - if ( bExp =3D=3D 0 ) { - if ( bSig =3D=3D 0 ) { - float_raise(float_flag_invalid, status); - return float32_default_nan(status); - } - normalizeFloat32Subnormal( bSig, &bExp, &bSig ); - } - if ( aExp =3D=3D 0 ) { - if ( aSig =3D=3D 0 ) return a; - normalizeFloat32Subnormal( aSig, &aExp, &aSig ); - } - expDiff =3D aExp - bExp; - aSig |=3D 0x00800000; - bSig |=3D 0x00800000; - if ( expDiff < 32 ) { - aSig <<=3D 8; - bSig <<=3D 8; - if ( expDiff < 0 ) { - if ( expDiff < -1 ) return a; - aSig >>=3D 1; - } - q =3D ( bSig <=3D aSig ); - if ( q ) aSig -=3D bSig; - if ( 0 < expDiff ) { - q =3D ( ( (uint64_t) aSig )<<32 ) / bSig; - q >>=3D 32 - expDiff; - bSig >>=3D 2; - aSig =3D ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q; - } - else { - aSig >>=3D 2; - bSig >>=3D 2; - } - } - else { - if ( bSig <=3D aSig ) aSig -=3D bSig; - aSig64 =3D ( (uint64_t) aSig )<<40; - bSig64 =3D ( (uint64_t) bSig )<<40; - expDiff -=3D 64; - while ( 0 < expDiff ) { - q64 =3D estimateDiv128To64( aSig64, 0, bSig64 ); - q64 =3D ( 2 < q64 ) ? q64 - 2 : 0; - aSig64 =3D - ( ( bSig * q64 )<<38 ); - expDiff -=3D 62; - } - expDiff +=3D 64; - q64 =3D estimateDiv128To64( aSig64, 0, bSig64 ); - q64 =3D ( 2 < q64 ) ? q64 - 2 : 0; - q =3D q64>>( 64 - expDiff ); - bSig <<=3D 6; - aSig =3D ( ( aSig64>>33 )<<( expDiff - 1 ) ) - bSig * q; - } - do { - alternateASig =3D aSig; - ++q; - aSig -=3D bSig; - } while ( 0 <=3D (int32_t) aSig ); - sigMean =3D aSig + alternateASig; - if ( ( sigMean < 0 ) || ( ( sigMean =3D=3D 0 ) && ( q & 1 ) ) ) { - aSig =3D alternateASig; - } - zSign =3D ( (int32_t) aSig < 0 ); - if ( zSign ) aSig =3D - aSig; - return normalizeRoundAndPackFloat32(aSign ^ zSign, bExp, aSig, status); -} - - - /*------------------------------------------------------------------------= ---- | Returns the binary exponential of the single-precision floating-point va= lue | `a'. The operation is performed according to the IEC/IEEE Standard for @@ -5273,94 +4804,6 @@ float32 float32_exp2(float32 a, float_status *status) return float32_round_pack_canonical(&rp, status); } =20 -/*------------------------------------------------------------------------= ---- -| Returns the remainder of the double-precision floating-point value `a' -| with respect to the corresponding value `b'. The operation is performed -| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. -*-------------------------------------------------------------------------= ---*/ - -float64 float64_rem(float64 a, float64 b, float_status *status) -{ - bool aSign, zSign; - int aExp, bExp, expDiff; - uint64_t aSig, bSig; - uint64_t q, alternateASig; - int64_t sigMean; - - a =3D float64_squash_input_denormal(a, status); - b =3D float64_squash_input_denormal(b, status); - aSig =3D extractFloat64Frac( a ); - aExp =3D extractFloat64Exp( a ); - aSign =3D extractFloat64Sign( a ); - bSig =3D extractFloat64Frac( b ); - bExp =3D extractFloat64Exp( b ); - if ( aExp =3D=3D 0x7FF ) { - if ( aSig || ( ( bExp =3D=3D 0x7FF ) && bSig ) ) { - return propagateFloat64NaN(a, b, status); - } - float_raise(float_flag_invalid, status); - return float64_default_nan(status); - } - if ( bExp =3D=3D 0x7FF ) { - if (bSig) { - return propagateFloat64NaN(a, b, status); - } - return a; - } - if ( bExp =3D=3D 0 ) { - if ( bSig =3D=3D 0 ) { - float_raise(float_flag_invalid, status); - return float64_default_nan(status); - } - normalizeFloat64Subnormal( bSig, &bExp, &bSig ); - } - if ( aExp =3D=3D 0 ) { - if ( aSig =3D=3D 0 ) return a; - normalizeFloat64Subnormal( aSig, &aExp, &aSig ); - } - expDiff =3D aExp - bExp; - aSig =3D (aSig | UINT64_C(0x0010000000000000)) << 11; - bSig =3D (bSig | UINT64_C(0x0010000000000000)) << 11; - if ( expDiff < 0 ) { - if ( expDiff < -1 ) return a; - aSig >>=3D 1; - } - q =3D ( bSig <=3D aSig ); - if ( q ) aSig -=3D bSig; - expDiff -=3D 64; - while ( 0 < expDiff ) { - q =3D estimateDiv128To64( aSig, 0, bSig ); - q =3D ( 2 < q ) ? q - 2 : 0; - aSig =3D - ( ( bSig>>2 ) * q ); - expDiff -=3D 62; - } - expDiff +=3D 64; - if ( 0 < expDiff ) { - q =3D estimateDiv128To64( aSig, 0, bSig ); - q =3D ( 2 < q ) ? q - 2 : 0; - q >>=3D 64 - expDiff; - bSig >>=3D 2; - aSig =3D ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q; - } - else { - aSig >>=3D 2; - bSig >>=3D 2; - } - do { - alternateASig =3D aSig; - ++q; - aSig -=3D bSig; - } while ( 0 <=3D (int64_t) aSig ); - sigMean =3D aSig + alternateASig; - if ( ( sigMean < 0 ) || ( ( sigMean =3D=3D 0 ) && ( q & 1 ) ) ) { - aSig =3D alternateASig; - } - zSign =3D ( (int64_t) aSig < 0 ); - if ( zSign ) aSig =3D - aSig; - return normalizeRoundAndPackFloat64(aSign ^ zSign, bExp, aSig, status); - -} - /*------------------------------------------------------------------------= ---- | Rounds the extended double-precision floating-point value `a' | to the precision provided by floatx80_rounding_precision and returns the @@ -5379,266 +4822,6 @@ floatx80 floatx80_round(floatx80 a, float_status *s= tatus) return floatx80_round_pack_canonical(&p, status); } =20 -/*------------------------------------------------------------------------= ---- -| Returns the remainder of the extended double-precision floating-point va= lue -| `a' with respect to the corresponding value `b'. The operation is perfo= rmed -| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic, -| if 'mod' is false; if 'mod' is true, return the remainder based on trunc= ating -| the quotient toward zero instead. '*quotient' is set to the low 64 bits= of -| the absolute value of the integer quotient. -*-------------------------------------------------------------------------= ---*/ - -floatx80 floatx80_modrem(floatx80 a, floatx80 b, bool mod, uint64_t *quoti= ent, - float_status *status) -{ - bool aSign, zSign; - int32_t aExp, bExp, expDiff, aExpOrig; - uint64_t aSig0, aSig1, bSig; - uint64_t q, term0, term1, alternateASig0, alternateASig1; - - *quotient =3D 0; - if (floatx80_invalid_encoding(a) || floatx80_invalid_encoding(b)) { - float_raise(float_flag_invalid, status); - return floatx80_default_nan(status); - } - aSig0 =3D extractFloatx80Frac( a ); - aExpOrig =3D aExp =3D extractFloatx80Exp( a ); - aSign =3D extractFloatx80Sign( a ); - bSig =3D extractFloatx80Frac( b ); - bExp =3D extractFloatx80Exp( b ); - if ( aExp =3D=3D 0x7FFF ) { - if ( (uint64_t) ( aSig0<<1 ) - || ( ( bExp =3D=3D 0x7FFF ) && (uint64_t) ( bSig<<1 ) ) ) { - return propagateFloatx80NaN(a, b, status); - } - goto invalid; - } - if ( bExp =3D=3D 0x7FFF ) { - if ((uint64_t)(bSig << 1)) { - return propagateFloatx80NaN(a, b, status); - } - if (aExp =3D=3D 0 && aSig0 >> 63) { - /* - * Pseudo-denormal argument must be returned in normalized - * form. - */ - return packFloatx80(aSign, 1, aSig0); - } - return a; - } - if ( bExp =3D=3D 0 ) { - if ( bSig =3D=3D 0 ) { - invalid: - float_raise(float_flag_invalid, status); - return floatx80_default_nan(status); - } - normalizeFloatx80Subnormal( bSig, &bExp, &bSig ); - } - if ( aExp =3D=3D 0 ) { - if ( aSig0 =3D=3D 0 ) return a; - normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 ); - } - zSign =3D aSign; - expDiff =3D aExp - bExp; - aSig1 =3D 0; - if ( expDiff < 0 ) { - if ( mod || expDiff < -1 ) { - if (aExp =3D=3D 1 && aExpOrig =3D=3D 0) { - /* - * Pseudo-denormal argument must be returned in - * normalized form. - */ - return packFloatx80(aSign, aExp, aSig0); - } - return a; - } - shift128Right( aSig0, 0, 1, &aSig0, &aSig1 ); - expDiff =3D 0; - } - *quotient =3D q =3D ( bSig <=3D aSig0 ); - if ( q ) aSig0 -=3D bSig; - expDiff -=3D 64; - while ( 0 < expDiff ) { - q =3D estimateDiv128To64( aSig0, aSig1, bSig ); - q =3D ( 2 < q ) ? q - 2 : 0; - mul64To128( bSig, q, &term0, &term1 ); - sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 ); - shortShift128Left( aSig0, aSig1, 62, &aSig0, &aSig1 ); - expDiff -=3D 62; - *quotient <<=3D 62; - *quotient +=3D q; - } - expDiff +=3D 64; - if ( 0 < expDiff ) { - q =3D estimateDiv128To64( aSig0, aSig1, bSig ); - q =3D ( 2 < q ) ? q - 2 : 0; - q >>=3D 64 - expDiff; - mul64To128( bSig, q<<( 64 - expDiff ), &term0, &term1 ); - sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 ); - shortShift128Left( 0, bSig, 64 - expDiff, &term0, &term1 ); - while ( le128( term0, term1, aSig0, aSig1 ) ) { - ++q; - sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 ); - } - if (expDiff < 64) { - *quotient <<=3D expDiff; - } else { - *quotient =3D 0; - } - *quotient +=3D q; - } - else { - term1 =3D 0; - term0 =3D bSig; - } - if (!mod) { - sub128( term0, term1, aSig0, aSig1, &alternateASig0, &alternateASi= g1 ); - if ( lt128( alternateASig0, alternateASig1, aSig0, aSig1 ) - || ( eq128( alternateASig0, alternateASig1, aSig0, aSig= 1 ) - && ( q & 1 ) ) - ) { - aSig0 =3D alternateASig0; - aSig1 =3D alternateASig1; - zSign =3D ! zSign; - ++*quotient; - } - } - return - normalizeRoundAndPackFloatx80( - floatx80_precision_x, zSign, bExp + expDiff, aSig0, aSig1, sta= tus); - -} - -/*------------------------------------------------------------------------= ---- -| Returns the remainder of the extended double-precision floating-point va= lue -| `a' with respect to the corresponding value `b'. The operation is perfo= rmed -| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. -*-------------------------------------------------------------------------= ---*/ - -floatx80 floatx80_rem(floatx80 a, floatx80 b, float_status *status) -{ - uint64_t quotient; - return floatx80_modrem(a, b, false, "ient, status); -} - -/*------------------------------------------------------------------------= ---- -| Returns the remainder of the extended double-precision floating-point va= lue -| `a' with respect to the corresponding value `b', with the quotient trunc= ated -| toward zero. -*-------------------------------------------------------------------------= ---*/ - -floatx80 floatx80_mod(floatx80 a, floatx80 b, float_status *status) -{ - uint64_t quotient; - return floatx80_modrem(a, b, true, "ient, status); -} - -/*------------------------------------------------------------------------= ---- -| Returns the remainder of the quadruple-precision floating-point value `a' -| with respect to the corresponding value `b'. The operation is performed -| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic. -*-------------------------------------------------------------------------= ---*/ - -float128 float128_rem(float128 a, float128 b, float_status *status) -{ - bool aSign, zSign; - int32_t aExp, bExp, expDiff; - uint64_t aSig0, aSig1, bSig0, bSig1, q, term0, term1, term2; - uint64_t allZero, alternateASig0, alternateASig1, sigMean1; - int64_t sigMean0; - - aSig1 =3D extractFloat128Frac1( a ); - aSig0 =3D extractFloat128Frac0( a ); - aExp =3D extractFloat128Exp( a ); - aSign =3D extractFloat128Sign( a ); - bSig1 =3D extractFloat128Frac1( b ); - bSig0 =3D extractFloat128Frac0( b ); - bExp =3D extractFloat128Exp( b ); - if ( aExp =3D=3D 0x7FFF ) { - if ( ( aSig0 | aSig1 ) - || ( ( bExp =3D=3D 0x7FFF ) && ( bSig0 | bSig1 ) ) ) { - return propagateFloat128NaN(a, b, status); - } - goto invalid; - } - if ( bExp =3D=3D 0x7FFF ) { - if (bSig0 | bSig1) { - return propagateFloat128NaN(a, b, status); - } - return a; - } - if ( bExp =3D=3D 0 ) { - if ( ( bSig0 | bSig1 ) =3D=3D 0 ) { - invalid: - float_raise(float_flag_invalid, status); - return float128_default_nan(status); - } - normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 ); - } - if ( aExp =3D=3D 0 ) { - if ( ( aSig0 | aSig1 ) =3D=3D 0 ) return a; - normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 ); - } - expDiff =3D aExp - bExp; - if ( expDiff < -1 ) return a; - shortShift128Left( - aSig0 | UINT64_C(0x0001000000000000), - aSig1, - 15 - ( expDiff < 0 ), - &aSig0, - &aSig1 - ); - shortShift128Left( - bSig0 | UINT64_C(0x0001000000000000), bSig1, 15, &bSig0, &bSig1 ); - q =3D le128( bSig0, bSig1, aSig0, aSig1 ); - if ( q ) sub128( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 ); - expDiff -=3D 64; - while ( 0 < expDiff ) { - q =3D estimateDiv128To64( aSig0, aSig1, bSig0 ); - q =3D ( 4 < q ) ? q - 4 : 0; - mul128By64To192( bSig0, bSig1, q, &term0, &term1, &term2 ); - shortShift192Left( term0, term1, term2, 61, &term1, &term2, &allZe= ro ); - shortShift128Left( aSig0, aSig1, 61, &aSig0, &allZero ); - sub128( aSig0, 0, term1, term2, &aSig0, &aSig1 ); - expDiff -=3D 61; - } - if ( -64 < expDiff ) { - q =3D estimateDiv128To64( aSig0, aSig1, bSig0 ); - q =3D ( 4 < q ) ? q - 4 : 0; - q >>=3D - expDiff; - shift128Right( bSig0, bSig1, 12, &bSig0, &bSig1 ); - expDiff +=3D 52; - if ( expDiff < 0 ) { - shift128Right( aSig0, aSig1, - expDiff, &aSig0, &aSig1 ); - } - else { - shortShift128Left( aSig0, aSig1, expDiff, &aSig0, &aSig1 ); - } - mul128By64To192( bSig0, bSig1, q, &term0, &term1, &term2 ); - sub128( aSig0, aSig1, term1, term2, &aSig0, &aSig1 ); - } - else { - shift128Right( aSig0, aSig1, 12, &aSig0, &aSig1 ); - shift128Right( bSig0, bSig1, 12, &bSig0, &bSig1 ); - } - do { - alternateASig0 =3D aSig0; - alternateASig1 =3D aSig1; - ++q; - sub128( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 ); - } while ( 0 <=3D (int64_t) aSig0 ); - add128( - aSig0, aSig1, alternateASig0, alternateASig1, (uint64_t *)&sigMean= 0, &sigMean1 ); - if ( ( sigMean0 < 0 ) - || ( ( ( sigMean0 | sigMean1 ) =3D=3D 0 ) && ( q & 1 ) ) ) { - aSig0 =3D alternateASig0; - aSig1 =3D alternateASig1; - } - zSign =3D ( (int64_t) aSig0 < 0 ); - if ( zSign ) sub128( 0, 0, aSig0, aSig1, &aSig0, &aSig1 ); - return normalizeRoundAndPackFloat128(aSign ^ zSign, bExp - 4, aSig0, a= Sig1, - status); -} static void __attribute__((constructor)) softfloat_init(void) { union_float64 ua, ub, uc, ur; diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc index d1bd5c6edf..dddee92d6e 100644 --- a/fpu/softfloat-parts.c.inc +++ b/fpu/softfloat-parts.c.inc @@ -626,6 +626,40 @@ static FloatPartsN *partsN(div)(FloatPartsN *a, FloatP= artsN *b, return a; } =20 +/* + * Floating point remainder, per IEC/IEEE, or modulus. + */ +static FloatPartsN *partsN(modrem)(FloatPartsN *a, FloatPartsN *b, + uint64_t *mod_quot, float_status *s) +{ + int ab_mask =3D float_cmask(a->cls) | float_cmask(b->cls); + + if (likely(ab_mask =3D=3D float_cmask_normal)) { + frac_modrem(a, b, mod_quot); + return a; + } + + if (mod_quot) { + *mod_quot =3D 0; + } + + /* All the NaN cases */ + if (unlikely(ab_mask & float_cmask_anynan)) { + return parts_pick_nan(a, b, s); + } + + /* Inf % N; N % 0 */ + if (a->cls =3D=3D float_class_inf || b->cls =3D=3D float_class_zero) { + float_raise(float_flag_invalid, s); + parts_default_nan(a, s); + return a; + } + + /* N % Inf; 0 % N */ + g_assert(b->cls =3D=3D float_class_inf || a->cls =3D=3D float_class_ze= ro); + return a; +} + /* * Square Root * diff --git a/fpu/softfloat-specialize.c.inc b/fpu/softfloat-specialize.c.inc index 95e5325f67..12467bb9bb 100644 --- a/fpu/softfloat-specialize.c.inc +++ b/fpu/softfloat-specialize.c.inc @@ -641,62 +641,6 @@ static int pickNaNMulAdd(FloatClass a_cls, FloatClass = b_cls, FloatClass c_cls, #endif } =20 -/*------------------------------------------------------------------------= ---- -| Takes two single-precision floating-point values `a' and `b', one of whi= ch -| is a NaN, and returns the appropriate NaN result. If either `a' or `b' = is a -| signaling NaN, the invalid exception is raised. -*-------------------------------------------------------------------------= ---*/ - -static float32 propagateFloat32NaN(float32 a, float32 b, float_status *sta= tus) -{ - bool aIsLargerSignificand; - uint32_t av, bv; - FloatClass a_cls, b_cls; - - /* This is not complete, but is good enough for pickNaN. */ - a_cls =3D (!float32_is_any_nan(a) - ? float_class_normal - : float32_is_signaling_nan(a, status) - ? float_class_snan - : float_class_qnan); - b_cls =3D (!float32_is_any_nan(b) - ? float_class_normal - : float32_is_signaling_nan(b, status) - ? float_class_snan - : float_class_qnan); - - av =3D float32_val(a); - bv =3D float32_val(b); - - if (is_snan(a_cls) || is_snan(b_cls)) { - float_raise(float_flag_invalid, status); - } - - if (status->default_nan_mode) { - return float32_default_nan(status); - } - - if ((uint32_t)(av << 1) < (uint32_t)(bv << 1)) { - aIsLargerSignificand =3D 0; - } else if ((uint32_t)(bv << 1) < (uint32_t)(av << 1)) { - aIsLargerSignificand =3D 1; - } else { - aIsLargerSignificand =3D (av < bv) ? 1 : 0; - } - - if (pickNaN(a_cls, b_cls, aIsLargerSignificand, status)) { - if (is_snan(b_cls)) { - return float32_silence_nan(b, status); - } - return b; - } else { - if (is_snan(a_cls)) { - return float32_silence_nan(a, status); - } - return a; - } -} - /*------------------------------------------------------------------------= ---- | Returns 1 if the double-precision floating-point value `a' is a quiet | NaN; otherwise returns 0. @@ -737,62 +681,6 @@ bool float64_is_signaling_nan(float64 a_, float_status= *status) } } =20 -/*------------------------------------------------------------------------= ---- -| Takes two double-precision floating-point values `a' and `b', one of whi= ch -| is a NaN, and returns the appropriate NaN result. If either `a' or `b' = is a -| signaling NaN, the invalid exception is raised. -*-------------------------------------------------------------------------= ---*/ - -static float64 propagateFloat64NaN(float64 a, float64 b, float_status *sta= tus) -{ - bool aIsLargerSignificand; - uint64_t av, bv; - FloatClass a_cls, b_cls; - - /* This is not complete, but is good enough for pickNaN. */ - a_cls =3D (!float64_is_any_nan(a) - ? float_class_normal - : float64_is_signaling_nan(a, status) - ? float_class_snan - : float_class_qnan); - b_cls =3D (!float64_is_any_nan(b) - ? float_class_normal - : float64_is_signaling_nan(b, status) - ? float_class_snan - : float_class_qnan); - - av =3D float64_val(a); - bv =3D float64_val(b); - - if (is_snan(a_cls) || is_snan(b_cls)) { - float_raise(float_flag_invalid, status); - } - - if (status->default_nan_mode) { - return float64_default_nan(status); - } - - if ((uint64_t)(av << 1) < (uint64_t)(bv << 1)) { - aIsLargerSignificand =3D 0; - } else if ((uint64_t)(bv << 1) < (uint64_t)(av << 1)) { - aIsLargerSignificand =3D 1; - } else { - aIsLargerSignificand =3D (av < bv) ? 1 : 0; - } - - if (pickNaN(a_cls, b_cls, aIsLargerSignificand, status)) { - if (is_snan(b_cls)) { - return float64_silence_nan(b, status); - } - return b; - } else { - if (is_snan(a_cls)) { - return float64_silence_nan(a, status); - } - return a; - } -} - /*------------------------------------------------------------------------= ---- | Returns 1 if the extended double-precision floating-point value `a' is a | quiet NaN; otherwise returns 0. This slightly differs from the same @@ -947,56 +835,3 @@ bool float128_is_signaling_nan(float128 a, float_statu= s *status) } } } - -/*------------------------------------------------------------------------= ---- -| Takes two quadruple-precision floating-point values `a' and `b', one of -| which is a NaN, and returns the appropriate NaN result. If either `a' or -| `b' is a signaling NaN, the invalid exception is raised. -*-------------------------------------------------------------------------= ---*/ - -static float128 propagateFloat128NaN(float128 a, float128 b, - float_status *status) -{ - bool aIsLargerSignificand; - FloatClass a_cls, b_cls; - - /* This is not complete, but is good enough for pickNaN. */ - a_cls =3D (!float128_is_any_nan(a) - ? float_class_normal - : float128_is_signaling_nan(a, status) - ? float_class_snan - : float_class_qnan); - b_cls =3D (!float128_is_any_nan(b) - ? float_class_normal - : float128_is_signaling_nan(b, status) - ? float_class_snan - : float_class_qnan); - - if (is_snan(a_cls) || is_snan(b_cls)) { - float_raise(float_flag_invalid, status); - } - - if (status->default_nan_mode) { - return float128_default_nan(status); - } - - if (lt128(a.high << 1, a.low, b.high << 1, b.low)) { - aIsLargerSignificand =3D 0; - } else if (lt128(b.high << 1, b.low, a.high << 1, a.low)) { - aIsLargerSignificand =3D 1; - } else { - aIsLargerSignificand =3D (a.high < b.high) ? 1 : 0; - } - - if (pickNaN(a_cls, b_cls, aIsLargerSignificand, status)) { - if (is_snan(b_cls)) { - return float128_silence_nan(b, status); - } - return b; - } else { - if (is_snan(a_cls)) { - return float128_silence_nan(a, status); - } - return a; - } -} --=20 2.25.1