rounded-hw 0.4.0.2 → 0.4.0.3
raw patch · 19 files changed
+1688/−107 lines, 19 filesdep ~QuickCheckdep ~basedep ~doctestPVP ok
version bump matches the API change (PVP)
Dependency ranges changed: QuickCheck, base, doctest, tagged
API changes (from Hackage documentation)
Files
- ChangeLog.md +13/−0
- LICENSE +1/−1
- README.md +1/−0
- benchmark/Conversion.hs +7/−0
- cbits/interval-prim-x86_64-avx512.S +37/−0
- cbits/interval-prim-x86_64-sse2.S +156/−66
- cbits/rounded-riscv.inl +1263/−0
- cbits/rounded.c +87/−11
- rounded-hw.cabal +17/−10
- src/Numeric/Rounded/Hardware/Backend/Default.hs +9/−6
- src/Numeric/Rounded/Hardware/Backend/ViaRational.hs +1/−1
- src/Numeric/Rounded/Hardware/Internal/FloatUtil.hs +3/−2
- src/Numeric/Rounded/Hardware/Interval.hs +17/−3
- src/Numeric/Rounded/Hardware/Interval/ElementaryFunctions.hs +2/−2
- src/Numeric/Rounded/Hardware/Interval/NonEmpty.hs +17/−3
- test/IntervalArithmeticNESpec.hs +44/−0
- test/IntervalArithmeticSpec.hs +6/−1
- test/RoundedArithmeticSpec.hs +5/−1
- test/Spec.hs +2/−0
ChangeLog.md view
@@ -1,5 +1,18 @@ # Changelog for rounded-hw +## 0.4.0.3 (2026-09-28)++* Fix range reduction of `cosI`/`sinI` (#9, by @sheaf).+* Fix `roundedSqrt Infinity` of `ViaRational`.+* Fix Interval `powInt` (correctness issue).+* Fix `distanceUlp` on `[-maxFinite, maxFinite]`.+* `powInt i 0` now return `I 1 1`.+* Make use of static rounding modes on RISC-V (#6).+* SSE2 backend now uses a fixed MXCSR value, ignoring ambient floating-point environment.+* Support GHC 10.0.+* Support no-TNTC GHC (#8).+* Support unregisterised GHC (#8).+ ## 0.4.0.2 (2025-12-29) * Support GHC 9.14.
LICENSE view
@@ -1,4 +1,4 @@-Copyright ARATA Mizuki (c) 2020-2024+Copyright ARATA Mizuki (c) 2020-2026 All rights reserved.
README.md view
@@ -62,6 +62,7 @@ * AVX512 EVEX encoding (`_mm_*_round_*`) * x87 Control Word (for x87 long double) * AArch64 FPCR+ * RISC-V static rounding mode * On x86_64, `foreign import prim` is used to provide faster interval addition/subtraction. By default, C FFI is used and an appropriate technology is detected.
benchmark/Conversion.hs view
@@ -1,3 +1,4 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE HexFloatLiterals #-} {-# LANGUAGE NumericUnderscores #-}@@ -11,7 +12,9 @@ import Numeric.Floating.IEEE import qualified Numeric.Floating.IEEE.Internal as IEEE.Internal import Numeric.Rounded.Hardware+#if defined(USE_FFI) import qualified Numeric.Rounded.Hardware.Backend.C as C+#endif import Numeric.Rounded.Hardware.Class import Numeric.Rounded.Hardware.Interval import Test.Tasty.Bench@@ -78,7 +81,9 @@ , bench "Interval/fromIntegralR" . nf (\n -> case IEEE.Internal.fromIntegralR n of Pair (IEEE.Internal.RoundTowardNegative x) (IEEE.Internal.RoundTowardPositive y) -> (x, y) :: (Double, Double) )+#if defined(USE_FFI) , bench "Interval/individual/C" . nf (\n -> (C.roundedDoubleFromInt64 TowardNegInf n, C.roundedDoubleFromInt64 TowardInf n))+#endif ] | (name, value) <- [ ("small", -2^50 + 2^13 + 127) , ("medium", -2^60 + 42 * 2^53 - 137 * 2^24 + 3)@@ -100,7 +105,9 @@ , bench "Interval/fromIntegralR" . nf (\n -> case IEEE.Internal.fromIntegralR n of Pair (IEEE.Internal.RoundTowardNegative x) (IEEE.Internal.RoundTowardPositive y) -> (x, y) :: (Double, Double) )+#if defined(USE_FFI) , bench "Interval/individual/C" . nf (\n -> (C.roundedDoubleFromWord64 TowardNegInf n, C.roundedDoubleFromWord64 TowardInf n))+#endif ] | (name, value) <- [ ("small", 2^50 + 2^13 + 127) , ("medium", 2^63 + 42 * 2^53 - 137 * 2^24 + 3)
cbits/interval-prim-x86_64-avx512.S view
@@ -21,11 +21,17 @@ # , Double# -- upper (%xmm2) # #) #+ .p2align 3 .globl SYMBOL(rounded_hw_interval_add) SYMBOL(rounded_hw_interval_add): vaddsd {rd-sae}, %xmm3, %xmm1, %xmm1 # xmm1 = xmm1[0] + xmm3[0], xmm1[1] (downward) vaddsd {ru-sae}, %xmm4, %xmm2, %xmm2 # xmm2 = xmm2[0] + xmm4[0], xmm2[1] (upward)+#if defined(TABLES_NEXT_TO_CODE) jmp *(%rbp)+#else+ movq (%rbp), %rax+ jmp *(%rax)+#endif # # rounded_hw_interval_sub@@ -37,11 +43,17 @@ # , Double# -- upper (%xmm2) # #) #+ .p2align 3 .globl SYMBOL(rounded_hw_interval_sub) SYMBOL(rounded_hw_interval_sub): vsubsd {rd-sae}, %xmm4, %xmm1, %xmm1 # xmm1 = xmm1[0] - xmm4[0], xmm1[1] (downward) vsubsd {ru-sae}, %xmm3, %xmm2, %xmm2 # xmm2 = xmm2[0] - xmm3[0], xmm2[1] (upward)+#if defined(TABLES_NEXT_TO_CODE) jmp *(%rbp)+#else+ movq (%rbp), %rax+ jmp *(%rax)+#endif # # rounded_hw_interval_recip@@ -51,13 +63,20 @@ # , Double# -- upper (%xmm2) # #) #+ .p2align 3 .globl SYMBOL(rounded_hw_interval_recip) SYMBOL(rounded_hw_interval_recip): vmovsd LC0(%rip), %xmm4 # xmm4 = 1.0, zero vdivsd {rd-sae}, %xmm2, %xmm4, %xmm3 # xmm3 = xmm4[0] / xmm2[0], xmm4[1] (downward) vdivsd {ru-sae}, %xmm1, %xmm4, %xmm2 # xmm2 = xmm4[0] / xmm1[0], xmm4[1] (upward) vmovapd %xmm3, %xmm1 # xmm1 = xmm3+#if defined(TABLES_NEXT_TO_CODE) jmp *(%rbp)+#else+ movq (%rbp), %rax+ jmp *(%rax)+#endif+ .p2align 2 LC0: .quad 0x3FF0000000000000 # 1.0 in binary64 # 0b0011_1111_1111_0000_..._0000@@ -74,13 +93,19 @@ # , Double# -- upper (%xmm2) # #) #+ .p2align 3 .globl SYMBOL(rounded_hw_interval_sqrt) SYMBOL(rounded_hw_interval_sqrt): vmovq %xmm1, %xmm1 # xmm1 = xmm1[0], zero vsqrtsd {rd-sae}, %xmm1, %xmm1, %xmm1 # xmm1 = sqrt(xmm1[0]), xmm1[1] (downward) vmovq %xmm2, %xmm2 # xmm2 = xmm2[0], zero vsqrtsd {ru-sae}, %xmm2, %xmm2, %xmm2 # xmm2 = sqrt(xmm2[0]), xmm2[1] (upward)+#if defined(TABLES_NEXT_TO_CODE) jmp *(%rbp)+#else+ movq (%rbp), %rax+ jmp *(%rax)+#endif # # rounded_hw_interval_from_int64@@ -89,12 +114,18 @@ # , Double# -- upper (%xmm2) # #) #+ .p2align 3 .globl SYMBOL(rounded_hw_interval_from_int64) SYMBOL(rounded_hw_interval_from_int64): vxorps %xmm2, %xmm2, %xmm2 # xmm2 = zero vcvtsi2sdq %rbx, {rd-sae}, %xmm2, %xmm1 # xmm1 = (double)(int64)rbx, xmm2[1] (downward) vcvtsi2sdq %rbx, {ru-sae}, %xmm2, %xmm2 # xmm2 = (double)(int64)rbx, xmm2[1] (upward)+#if defined(TABLES_NEXT_TO_CODE) jmp *(%rbp)+#else+ movq (%rbp), %rax+ jmp *(%rax)+#endif # # rounded_hw_interval_from_word64@@ -103,9 +134,15 @@ # , Double# -- upper (%xmm2) # #) #+ .p2align 3 .globl SYMBOL(rounded_hw_interval_from_word64) SYMBOL(rounded_hw_interval_from_word64): vxorps %xmm2, %xmm2, %xmm2 # xmm2 = zero vcvtusi2sdq %rbx, {rd-sae}, %xmm2, %xmm1 # xmm1 = (double)(uint64)rbx, xmm2[1] (downward) vcvtusi2sdq %rbx, {ru-sae}, %xmm2, %xmm2 # xmm2 = (double)(uint64)rbx, xmm2[1] (upward)+#if defined(TABLES_NEXT_TO_CODE) jmp *(%rbp)+#else+ movq (%rbp), %rax+ jmp *(%rax)+#endif
cbits/interval-prim-x86_64-sse2.S view
@@ -7,10 +7,30 @@ #define SYMBOL(name) name #endif +#define USE_STATIC_MXCSR++ # MXCSR value:+ # 0b0_01_111111_0_000000 = 0x3f80+ # ^ ^^ ^^^^^^ ^ ^^^^^^+ # | | | | +--- Flags: PE,UE,OE,ZE,DE,IE+ # | | | +-------- DAZ (Denormals are Zeros)+ # | | +------------ Masks: PM,UM,OM,ZM,DM,IM+ # | +----------------- Rounding Control+ # +-------------------- Flush to Zero+ .p2align 3+L_MXCSR_DOWNWARD:+ .long 0x00003f80+L_MXCSR_UPWARD:+ # 0b0_10_111111_0_000000+ .long 0x00005f80+ .globl SYMBOL(rounded_hw_interval_backend_name) SYMBOL(rounded_hw_interval_backend_name): .string "SSE2" +#define PREV_MXCSR 8(%rsp)+#define TMP_MXCSR 12(%rsp)+ # # rounded_hw_interval_add # :: Double# -- lower 1 (%xmm1)@@ -21,21 +41,35 @@ # , Double# -- upper (%xmm2) # #) #+ .p2align 4 .globl SYMBOL(rounded_hw_interval_add) SYMBOL(rounded_hw_interval_add):- stmxcsr -8(%rbp) # *(int32*)(rbp-8) = MXCSR- movl -8(%rbp), %ecx # ecx = *(int32*)(rbp-8)- andl $0x9FFF, %ecx # ecx = ecx & 0x9FFF; clear Rounding Control field- orl $0x2000, %ecx # ecx = ecx | 0x2000; set RC = downward- movl %ecx, -4(%rbp) # *(int32*)(rbp-4) = ecx- ldmxcsr -4(%rbp) # MXCSR = *(int32*)(rbp-4)- addsd %xmm3, %xmm1 # xmm1 = xmm1[0] + xmm3[0], xmm1[1]- xorl $0x6000, %ecx # ecx = ecx ^ 0x6000; downward -> upward- movl %ecx, -4(%rbp) # *(int32*)(rbp-4) = ecx- ldmxcsr -4(%rbp) # MXCSR = *(int32*)(rbp-4)- addsd %xmm4, %xmm2 # xmm2 = xmm2[0] + xmm4[0], xmm2[1]- ldmxcsr -8(%rbp) # MXCSR = *(int32*)(rbp-8)+ stmxcsr PREV_MXCSR # PREV_MXCSR = MXCSR+#if defined(USE_STATIC_MXCSR)+ ldmxcsr L_MXCSR_DOWNWARD(%rip)+#else+ movl PREV_MXCSR, %ecx # ecx = PREV_MXCSR+ andl $0x9FFF, %ecx # ecx = ecx & 0x9FFF; clear Rounding Control field+ orl $0x2000, %ecx # ecx = ecx | 0x2000; set RC = downward+ movl %ecx, TMP_MXCSR # TMP_MXCSR = ecx+ ldmxcsr TMP_MXCSR # MXCSR = TMP_MXCSR+#endif+ addsd %xmm3, %xmm1 # xmm1 = xmm1[0] + xmm3[0], xmm1[1]+#if defined(USE_STATIC_MXCSR)+ ldmxcsr L_MXCSR_UPWARD(%rip)+#else+ xorl $0x6000, %ecx # ecx = ecx ^ 0x6000; downward -> upward+ movl %ecx, TMP_MXCSR # TMP_MXCSR = ecx+ ldmxcsr TMP_MXCSR # MXCSR = TMP_MXCSR+#endif+ addsd %xmm4, %xmm2 # xmm2 = xmm2[0] + xmm4[0], xmm2[1]+ ldmxcsr PREV_MXCSR # MXCSR = PREV_MXCSR+#if defined(TABLES_NEXT_TO_CODE) jmp *(%rbp)+#else+ movq (%rbp), %rax+ jmp *(%rax)+#endif # # rounded_hw_interval_sub@@ -47,21 +81,35 @@ # , Double# -- upper (%xmm2) # #) #+ .p2align 4 .globl SYMBOL(rounded_hw_interval_sub) SYMBOL(rounded_hw_interval_sub):- stmxcsr -8(%rbp) # *(int32*)(rbp-8) = MXCSR- movl -8(%rbp), %ecx # ecx = *(int32*)(rbp-8)- andl $0x9FFF, %ecx # ecx = ecx & 0x9FFF; clear Rounding Control field- orl $0x2000, %ecx # ecx = ecx | 0x2000; set RC = downward- movl %ecx, -4(%rbp) # *(int32*)(rbp-4) = ecx- ldmxcsr -4(%rbp) # MXCSR = *(int32*)(rbp-4)- subsd %xmm4, %xmm1 # xmm1 = xmm1[0] - xmm4[0], xmm1[1]- xorl $0x6000, %ecx # ecx = ecx ^ 0x6000; downward -> upward- movl %ecx, -4(%rbp) # *(int32*)(rbp-4) = ecx- ldmxcsr -4(%rbp) # MXCSR = *(int32*)(rbp-4)- subsd %xmm3, %xmm2 # xmm2 = xmm2[0] - xmm3[0], xmm2[1]- ldmxcsr -8(%rbp) # MXCSR = *(int32*)(rbp-4)+ stmxcsr PREV_MXCSR # PREV_MXCSR = MXCSR+#if defined(USE_STATIC_MXCSR)+ ldmxcsr L_MXCSR_DOWNWARD(%rip)+#else+ movl PREV_MXCSR, %ecx # ecx = PREV_MXCSR+ andl $0x9FFF, %ecx # ecx = ecx & 0x9FFF; clear Rounding Control field+ orl $0x2000, %ecx # ecx = ecx | 0x2000; set RC = downward+ movl %ecx, TMP_MXCSR # TMP_MXCSR = ecx+ ldmxcsr TMP_MXCSR # MXCSR = TMP_MXCSR+#endif+ subsd %xmm4, %xmm1 # xmm1 = xmm1[0] - xmm4[0], xmm1[1]+#if defined(USE_STATIC_MXCSR)+ ldmxcsr L_MXCSR_UPWARD(%rip)+#else+ xorl $0x6000, %ecx # ecx = ecx ^ 0x6000; downward -> upward+ movl %ecx, TMP_MXCSR # TMP_MXCSR = ecx+ ldmxcsr TMP_MXCSR # MXCSR = TMP_MXCSR+#endif+ subsd %xmm3, %xmm2 # xmm2 = xmm2[0] - xmm3[0], xmm2[1]+ ldmxcsr PREV_MXCSR # MXCSR = PREV_MXCSR+#if defined(TABLES_NEXT_TO_CODE) jmp *(%rbp)+#else+ movq (%rbp), %rax+ jmp *(%rax)+#endif # # rounded_hw_interval_recip@@ -71,25 +119,39 @@ # , Double# -- upper (%xmm2) # #) #+ .p2align 4 .globl SYMBOL(rounded_hw_interval_recip) SYMBOL(rounded_hw_interval_recip):- stmxcsr -8(%rbp) # *(int32*)(rbp-8) = MXCSR- movl -8(%rbp), %ecx # ecx = *(int32*)(rbp-8)- andl $0x9FFF, %ecx # ecx = ecx & 0x9FFF; clear Rounding Control field- orl $0x2000, %ecx # ecx = ecx | 0x2000; set RC = downward- movl %ecx, -4(%rbp) # *(int32*)(rbp-4) = ecx- ldmxcsr -4(%rbp) # MXCSR = *(int32*)(rbp-4); set downward- movsd LC0(%rip), %xmm3 # xmm3 = (double)1.0,zero- movapd %xmm3, %xmm4 # xmm4 = xmm3- divsd %xmm2, %xmm3 # xmm3 = xmm3[0] / xmm2[0], xmm3[1]- xorl $0x6000, %ecx # ecx = ecx ^ 0x6000; downward -> upward- movl %ecx, -4(%rbp) # *(int32*)(rbp-4) = ecx- ldmxcsr -4(%rbp) # MXCSR = *(int32*)(rbp-4); set upward- divsd %xmm1, %xmm4 # xmm4 = xmm4[0] / xmm1[0], xmm4[1]- ldmxcsr -8(%rbp) # MXCSR = *(int32*)(rbp-8); restore- movapd %xmm3, %xmm1 # xmm1 = xmm3- movapd %xmm4, %xmm2 # xmm2 = xmm4+ stmxcsr PREV_MXCSR # PREV_MXCSR = MXCSR+#if defined(USE_STATIC_MXCSR)+ ldmxcsr L_MXCSR_DOWNWARD(%rip)+#else+ movl PREV_MXCSR, %ecx # ecx = PREV_MXCSR+ andl $0x9FFF, %ecx # ecx = ecx & 0x9FFF; clear Rounding Control field+ orl $0x2000, %ecx # ecx = ecx | 0x2000; set RC = downward+ movl %ecx, TMP_MXCSR # TMP_MXCSR = ecx+ ldmxcsr TMP_MXCSR # MXCSR = TMP_MXCSR; set downward+#endif+ movsd LC0(%rip), %xmm3 # xmm3 = (double)1.0,zero+ movapd %xmm3, %xmm4 # xmm4 = xmm3+ divsd %xmm2, %xmm3 # xmm3 = xmm3[0] / xmm2[0], xmm3[1]+#if defined(USE_STATIC_MXCSR)+ ldmxcsr L_MXCSR_UPWARD(%rip)+#else+ xorl $0x6000, %ecx # ecx = ecx ^ 0x6000; downward -> upward+ movl %ecx, TMP_MXCSR # TMP_MXCSR = ecx+ ldmxcsr TMP_MXCSR # MXCSR = TMP_MXCSR; set upward+#endif+ divsd %xmm1, %xmm4 # xmm4 = xmm4[0] / xmm1[0], xmm4[1]+ ldmxcsr PREV_MXCSR # MXCSR = PREV_MXCSR; restore+ movapd %xmm3, %xmm1 # xmm1 = xmm3+ movapd %xmm4, %xmm2 # xmm2 = xmm4+#if defined(TABLES_NEXT_TO_CODE) jmp *(%rbp)+#else+ movq (%rbp), %rax+ jmp *(%rax)+#endif LC0: .quad 0x3FF0000000000000 # 1.0 in binary64 # 0b0011_1111_1111_0000_..._0000@@ -106,21 +168,35 @@ # , Double# -- upper (%xmm2) # #) #+ .p2align 4 .globl SYMBOL(rounded_hw_interval_sqrt) SYMBOL(rounded_hw_interval_sqrt):- stmxcsr -8(%rbp) # *(int32*)(rbp-8) = MXCSR- movl -8(%rbp), %ecx # ecx = *(int32*)(rbp-8)- andl $0x9FFF, %ecx # ecx = ecx & 0x9FFF; clear Rounding Control field- orl $0x2000, %ecx # ecx = ecx | 0x2000; set RC = downward- movl %ecx, -4(%rbp) # *(int32*)(rbp-4) = ecx- ldmxcsr -4(%rbp) # MXCSR = *(int32*)(rbp-4); set downward- sqrtsd %xmm1, %xmm1 # xmm1 = sqrt(xmm1[0]), xmm1[1]- xorl $0x6000, %ecx # ecx = ecx ^ 0x6000; downward -> upward- movl %ecx, -4(%rbp) # *(int32*)(rbp-4) = ecx- ldmxcsr -4(%rbp) # MXCSR = *(int32*)(rbp-4); set upward- sqrtsd %xmm2, %xmm2 # xmm2 = sqrt(xmm2[0]), xmm2[1]- ldmxcsr -8(%rbp) # MXCSR = *(int32*)(rbp-8); restore+ stmxcsr PREV_MXCSR # PREV_MXCSR = MXCSR+#if defined(USE_STATIC_MXCSR)+ ldmxcsr L_MXCSR_DOWNWARD(%rip)+#else+ movl PREV_MXCSR, %ecx # ecx = PREV_MXCSR+ andl $0x9FFF, %ecx # ecx = ecx & 0x9FFF; clear Rounding Control field+ orl $0x2000, %ecx # ecx = ecx | 0x2000; set RC = downward+ movl %ecx, TMP_MXCSR # TMP_MXCSR = ecx+ ldmxcsr TMP_MXCSR # MXCSR = TMP_MXCSR; set downward+#endif+ sqrtsd %xmm1, %xmm1 # xmm1 = sqrt(xmm1[0]), xmm1[1]+#if defined(USE_STATIC_MXCSR)+ ldmxcsr L_MXCSR_UPWARD(%rip)+#else+ xorl $0x6000, %ecx # ecx = ecx ^ 0x6000; downward -> upward+ movl %ecx, TMP_MXCSR # TMP_MXCSR = ecx+ ldmxcsr TMP_MXCSR # MXCSR = TMP_MXCSR; set upward+#endif+ sqrtsd %xmm2, %xmm2 # xmm2 = sqrt(xmm2[0]), xmm2[1]+ ldmxcsr PREV_MXCSR # MXCSR = PREV_MXCSR; restore+#if defined(TABLES_NEXT_TO_CODE) jmp *(%rbp)+#else+ movq (%rbp), %rax+ jmp *(%rax)+#endif # # rounded_hw_interval_from_int64@@ -129,20 +205,34 @@ # , Double# -- upper (%xmm2) # #) #+ .p2align 4 .globl SYMBOL(rounded_hw_interval_from_int64) SYMBOL(rounded_hw_interval_from_int64):- stmxcsr -8(%rbp) # *(int32*)(rbp-8) = MXCSR- movl -8(%rbp), %ecx # ecx = *(int32*)(rbp-8)- andl $0x9FFF, %ecx # ecx = ecx & 0x9FFF; clear Rounding Control field- orl $0x2000, %ecx # ecx = ecx | 0x2000; set RC = downward- movl %ecx, -4(%rbp) # *(int32*)(rbp-4) = ecx- ldmxcsr -4(%rbp) # MXCSR = *(int32*)(rbp-4); set downward- pxor %xmm1, %xmm1 # xmm1 = zero- cvtsi2sdq %rbx, %xmm1 # xmm1 = (double)(int64)rbx, xmm1[1]- xorl $0x6000, %ecx # ecx = ecx ^ 0x6000; downward -> upward- movl %ecx, -4(%rbp) # *(int32*)(rbp-4) = ecx- ldmxcsr -4(%rbp) # MXCSR = *(int32*)(rbp-4); set upward- pxor %xmm2, %xmm2 # xmm2 = zero- cvtsi2sdq %rbx, %xmm2 # xmm2 = (double)(int64)rbx, xmm2[1]- ldmxcsr -8(%rbp) # MXCSR = *(int32*)(rbp-8); restore+ stmxcsr PREV_MXCSR # PREV_MXCSR = MXCSR+#if defined(USE_STATIC_MXCSR)+ ldmxcsr L_MXCSR_DOWNWARD(%rip)+#else+ movl PREV_MXCSR, %ecx # ecx = PREV_MXCSR+ andl $0x9FFF, %ecx # ecx = ecx & 0x9FFF; clear Rounding Control field+ orl $0x2000, %ecx # ecx = ecx | 0x2000; set RC = downward+ movl %ecx, TMP_MXCSR # TMP_MXCSR = ecx+ ldmxcsr TMP_MXCSR # MXCSR = TMP_MXCSR; set downward+#endif+ pxor %xmm1, %xmm1 # xmm1 = zero+ cvtsi2sdq %rbx, %xmm1 # xmm1 = (double)(int64)rbx, xmm1[1]+#if defined(USE_STATIC_MXCSR)+ ldmxcsr L_MXCSR_UPWARD(%rip)+#else+ xorl $0x6000, %ecx # ecx = ecx ^ 0x6000; downward -> upward+ movl %ecx, TMP_MXCSR # TMP_MXCSR = ecx+ ldmxcsr TMP_MXCSR # MXCSR = TMP_MXCSR; set upward+#endif+ pxor %xmm2, %xmm2 # xmm2 = zero+ cvtsi2sdq %rbx, %xmm2 # xmm2 = (double)(int64)rbx, xmm2[1]+ ldmxcsr PREV_MXCSR # MXCSR = PREV_MXCSR; restore+#if defined(TABLES_NEXT_TO_CODE) jmp *(%rbp)+#else+ movq (%rbp), %rax+ jmp *(%rax)+#endif
+ cbits/rounded-riscv.inl view
@@ -0,0 +1,1263 @@+/* This file was generated by etc/gen-rounded-riscv.sh. */++//+// double+//++static inline ALWAYS_INLINE+double rounded_add_impl_double(native_rounding_mode mode, double a, double b)+{+ double result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fadd.d %0, %1, %2, rne" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_DOWNWARD:+ __asm__("fadd.d %0, %1, %2, rdn" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_UPWARD:+ __asm__("fadd.d %0, %1, %2, rup" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fadd.d %0, %1, %2, rtz" : "=f"(result) : "f"(a), "f"(b));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern double rounded_hw_add_double(HsInt mode, double a, double b)+{ return rounded_add_impl_double(hs_rounding_mode_to_native(mode), a, b); }+extern double rounded_hw_add_double_up(double a, double b)+{ return rounded_add_impl_double(ROUND_UPWARD, a, b); }+extern double rounded_hw_add_double_down(double a, double b)+{ return rounded_add_impl_double(ROUND_DOWNWARD, a, b); }+extern double rounded_hw_add_double_zero(double a, double b)+{ return rounded_add_impl_double(ROUND_TOWARDZERO, a, b); }++static inline ALWAYS_INLINE+double rounded_sub_impl_double(native_rounding_mode mode, double a, double b)+{+ double result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fsub.d %0, %1, %2, rne" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_DOWNWARD:+ __asm__("fsub.d %0, %1, %2, rdn" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_UPWARD:+ __asm__("fsub.d %0, %1, %2, rup" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fsub.d %0, %1, %2, rtz" : "=f"(result) : "f"(a), "f"(b));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern double rounded_hw_sub_double(HsInt mode, double a, double b)+{ return rounded_sub_impl_double(hs_rounding_mode_to_native(mode), a, b); }+extern double rounded_hw_sub_double_up(double a, double b)+{ return rounded_sub_impl_double(ROUND_UPWARD, a, b); }+extern double rounded_hw_sub_double_down(double a, double b)+{ return rounded_sub_impl_double(ROUND_DOWNWARD, a, b); }+extern double rounded_hw_sub_double_zero(double a, double b)+{ return rounded_sub_impl_double(ROUND_TOWARDZERO, a, b); }++static inline ALWAYS_INLINE+double rounded_mul_impl_double(native_rounding_mode mode, double a, double b)+{+ double result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fmul.d %0, %1, %2, rne" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_DOWNWARD:+ __asm__("fmul.d %0, %1, %2, rdn" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_UPWARD:+ __asm__("fmul.d %0, %1, %2, rup" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fmul.d %0, %1, %2, rtz" : "=f"(result) : "f"(a), "f"(b));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern double rounded_hw_mul_double(HsInt mode, double a, double b)+{ return rounded_mul_impl_double(hs_rounding_mode_to_native(mode), a, b); }+extern double rounded_hw_mul_double_up(double a, double b)+{ return rounded_mul_impl_double(ROUND_UPWARD, a, b); }+extern double rounded_hw_mul_double_down(double a, double b)+{ return rounded_mul_impl_double(ROUND_DOWNWARD, a, b); }+extern double rounded_hw_mul_double_zero(double a, double b)+{ return rounded_mul_impl_double(ROUND_TOWARDZERO, a, b); }++static inline ALWAYS_INLINE+double rounded_div_impl_double(native_rounding_mode mode, double a, double b)+{+ double result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fdiv.d %0, %1, %2, rne" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_DOWNWARD:+ __asm__("fdiv.d %0, %1, %2, rdn" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_UPWARD:+ __asm__("fdiv.d %0, %1, %2, rup" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fdiv.d %0, %1, %2, rtz" : "=f"(result) : "f"(a), "f"(b));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern double rounded_hw_div_double(HsInt mode, double a, double b)+{ return rounded_div_impl_double(hs_rounding_mode_to_native(mode), a, b); }+extern double rounded_hw_div_double_up(double a, double b)+{ return rounded_div_impl_double(ROUND_UPWARD, a, b); }+extern double rounded_hw_div_double_down(double a, double b)+{ return rounded_div_impl_double(ROUND_DOWNWARD, a, b); }+extern double rounded_hw_div_double_zero(double a, double b)+{ return rounded_div_impl_double(ROUND_TOWARDZERO, a, b); }++static inline ALWAYS_INLINE+double rounded_sqrt_impl_double(native_rounding_mode mode, double a)+{+ double result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fsqrt.d %0, %1, rne" : "=f"(result) : "f"(a));+ break;+ case ROUND_DOWNWARD:+ __asm__("fsqrt.d %0, %1, rdn" : "=f"(result) : "f"(a));+ break;+ case ROUND_UPWARD:+ __asm__("fsqrt.d %0, %1, rup" : "=f"(result) : "f"(a));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fsqrt.d %0, %1, rtz" : "=f"(result) : "f"(a));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern double rounded_hw_sqrt_double(HsInt mode, double a)+{ return rounded_sqrt_impl_double(hs_rounding_mode_to_native(mode), a); }+extern double rounded_hw_sqrt_double_up(double a)+{ return rounded_sqrt_impl_double(ROUND_UPWARD, a); }+extern double rounded_hw_sqrt_double_down(double a)+{ return rounded_sqrt_impl_double(ROUND_DOWNWARD, a); }+extern double rounded_hw_sqrt_double_zero(double a)+{ return rounded_sqrt_impl_double(ROUND_TOWARDZERO, a); }++static inline ALWAYS_INLINE+double rounded_fma_impl_double(native_rounding_mode mode, double a, double b, double c)+{+ double result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fmadd.d %0, %1, %2, %3, rne" : "=f"(result) : "f"(a), "f"(b), "f"(c));+ break;+ case ROUND_DOWNWARD:+ __asm__("fmadd.d %0, %1, %2, %3, rdn" : "=f"(result) : "f"(a), "f"(b), "f"(c));+ break;+ case ROUND_UPWARD:+ __asm__("fmadd.d %0, %1, %2, %3, rup" : "=f"(result) : "f"(a), "f"(b), "f"(c));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fmadd.d %0, %1, %2, %3, rtz" : "=f"(result) : "f"(a), "f"(b), "f"(c));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern double rounded_hw_fma_double(HsInt mode, double a, double b, double c)+{ return rounded_fma_impl_double(hs_rounding_mode_to_native(mode), a, b, c); }+extern double rounded_hw_fma_double_up(double a, double b, double c)+{ return rounded_fma_impl_double(ROUND_UPWARD, a, b, c); }+extern double rounded_hw_fma_double_down(double a, double b, double c)+{ return rounded_fma_impl_double(ROUND_DOWNWARD, a, b, c); }+extern double rounded_hw_fma_double_zero(double a, double b, double c)+{ return rounded_fma_impl_double(ROUND_TOWARDZERO, a, b, c); }++extern double rounded_hw_fma_if_fast_double(HsInt mode, double a, double b, double c)+{ return rounded_fma_impl_double(hs_rounding_mode_to_native(mode), a, b, c); }+extern double rounded_hw_fma_if_fast_double_up(double a, double b, double c)+{ return rounded_fma_impl_double(ROUND_UPWARD, a, b, c); }+extern double rounded_hw_fma_if_fast_double_down(double a, double b, double c)+{ return rounded_fma_impl_double(ROUND_DOWNWARD, a, b, c); }+extern double rounded_hw_fma_if_fast_double_zero(double a, double b, double c)+{ return rounded_fma_impl_double(ROUND_TOWARDZERO, a, b, c); }++//+// Conversion+//++static inline double rounded_int64_to_double_impl(native_rounding_mode mode, int64_t x)+{+ double result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fcvt.d.l %0, %1, rne" : "=f"(result) : "r"(x));+ break;+ case ROUND_DOWNWARD:+ __asm__("fcvt.d.l %0, %1, rdn" : "=f"(result) : "r"(x));+ break;+ case ROUND_UPWARD:+ __asm__("fcvt.d.l %0, %1, rup" : "=f"(result) : "r"(x));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fcvt.d.l %0, %1, rtz" : "=f"(result) : "r"(x));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern double rounded_hw_int64_to_double(HsInt mode, int64_t x)+{ return rounded_int64_to_double_impl(hs_rounding_mode_to_native(mode), x); }+extern double rounded_hw_int64_to_double_up(int64_t x)+{ return rounded_int64_to_double_impl(ROUND_UPWARD, x); }+extern double rounded_hw_int64_to_double_down(int64_t x)+{ return rounded_int64_to_double_impl(ROUND_DOWNWARD, x); }+extern double rounded_hw_int64_to_double_zero(int64_t x)+{ return rounded_int64_to_double_impl(ROUND_TOWARDZERO, x); }++static inline double rounded_word64_to_double_impl(native_rounding_mode mode, uint64_t x)+{+ double result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fcvt.d.lu %0, %1, rne" : "=f"(result) : "r"(x));+ break;+ case ROUND_DOWNWARD:+ __asm__("fcvt.d.lu %0, %1, rdn" : "=f"(result) : "r"(x));+ break;+ case ROUND_UPWARD:+ __asm__("fcvt.d.lu %0, %1, rup" : "=f"(result) : "r"(x));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fcvt.d.lu %0, %1, rtz" : "=f"(result) : "r"(x));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern double rounded_hw_word64_to_double(HsInt mode, uint64_t x)+{ return rounded_word64_to_double_impl(hs_rounding_mode_to_native(mode), x); }+extern double rounded_hw_word64_to_double_up(uint64_t x)+{ return rounded_word64_to_double_impl(ROUND_UPWARD, x); }+extern double rounded_hw_word64_to_double_down(uint64_t x)+{ return rounded_word64_to_double_impl(ROUND_DOWNWARD, x); }+extern double rounded_hw_word64_to_double_zero(uint64_t x)+{ return rounded_word64_to_double_impl(ROUND_TOWARDZERO, x); }++//+// Interval arithmetic+//++static inline double fast_fmax_double(double x, double y)+{+ // IEEE 754-2019 maximumNumber+ double result;+ __asm__("fmax.d %0, %1, %2" : "=f"(result) : "f"(x), "f"(y));+ return result;+}+static inline double fast_fmax4_double(double x, double y, double z, double w)+{+ return fast_fmax_double(fast_fmax_double(x, y), fast_fmax_double(z, w));+}+static inline double fast_fmin_double(double x, double y)+{+ // IEEE 754-2019 minimumNumber+ double result;+ __asm__("fmin.d %0, %1, %2" : "=f"(result) : "f"(x), "f"(y));+ return result;+}+static inline double fast_fmin4_double(double x, double y, double z, double w)+{+ return fast_fmin_double(fast_fmin_double(x, y), fast_fmin_double(z, w));+}++extern double rounded_hw_interval_mul_double_up(double lo1, double hi1, double lo2, double hi2)+{+ double x = rounded_mul_impl_double(ROUND_UPWARD, lo1, lo2);+ double y = rounded_mul_impl_double(ROUND_UPWARD, lo1, hi2);+ double z = rounded_mul_impl_double(ROUND_UPWARD, hi1, lo2);+ double w = rounded_mul_impl_double(ROUND_UPWARD, hi1, hi2);+ if (isnan(x)) x = 0.0; /* 0 * inf -> 0 */+ if (isnan(y)) y = 0.0; /* 0 * inf -> 0 */+ if (isnan(z)) z = 0.0; /* 0 * inf -> 0 */+ if (isnan(w)) w = 0.0; /* 0 * inf -> 0 */+ return fast_fmax4_double(x, y, z, w);+}++extern double rounded_hw_interval_mul_double_down(double lo1, double hi1, double lo2, double hi2)+{+ double x = rounded_mul_impl_double(ROUND_DOWNWARD, lo1, lo2);+ double y = rounded_mul_impl_double(ROUND_DOWNWARD, lo1, hi2);+ double z = rounded_mul_impl_double(ROUND_DOWNWARD, hi1, lo2);+ double w = rounded_mul_impl_double(ROUND_DOWNWARD, hi1, hi2);+ if (isnan(x)) x = 0.0; /* 0 * inf -> 0 */+ if (isnan(y)) y = 0.0; /* 0 * inf -> 0 */+ if (isnan(z)) z = 0.0; /* 0 * inf -> 0 */+ if (isnan(w)) w = 0.0; /* 0 * inf -> 0 */+ return fast_fmin4_double(x, y, z, w);+}++extern double rounded_hw_interval_mul_add_double_up(double lo1, double hi1, double lo2, double hi2, double hi3)+{+ double x = rounded_mul_impl_double(ROUND_UPWARD, lo1, lo2);+ double y = rounded_mul_impl_double(ROUND_UPWARD, lo1, hi2);+ double z = rounded_mul_impl_double(ROUND_UPWARD, hi1, lo2);+ double w = rounded_mul_impl_double(ROUND_UPWARD, hi1, hi2);+ if (isnan(x)) x = 0.0; /* 0 * inf -> 0 */+ if (isnan(y)) y = 0.0; /* 0 * inf -> 0 */+ if (isnan(z)) z = 0.0; /* 0 * inf -> 0 */+ if (isnan(w)) w = 0.0; /* 0 * inf -> 0 */+ double p = fast_fmax4_double(x, y, z, w);+ return rounded_add_impl_double(ROUND_UPWARD, p, hi3);+}++extern double rounded_hw_interval_mul_add_double_down(double lo1, double hi1, double lo2, double hi2, double lo3)+{+ double x = rounded_mul_impl_double(ROUND_DOWNWARD, lo1, lo2);+ double y = rounded_mul_impl_double(ROUND_DOWNWARD, lo1, hi2);+ double z = rounded_mul_impl_double(ROUND_DOWNWARD, hi1, lo2);+ double w = rounded_mul_impl_double(ROUND_DOWNWARD, hi1, hi2);+ if (isnan(x)) x = 0.0; /* 0 * inf -> 0 */+ if (isnan(y)) y = 0.0; /* 0 * inf -> 0 */+ if (isnan(z)) z = 0.0; /* 0 * inf -> 0 */+ if (isnan(w)) w = 0.0; /* 0 * inf -> 0 */+ double p = fast_fmin4_double(x, y, z, w);+ return rounded_add_impl_double(ROUND_DOWNWARD, p, lo3);+}++extern double rounded_hw_interval_div_double_up(double lo1, double hi1, double lo2, double hi2)+{+ double x = rounded_div_impl_double(ROUND_UPWARD, lo1, lo2);+ double y = rounded_div_impl_double(ROUND_UPWARD, lo1, hi2);+ double z = rounded_div_impl_double(ROUND_UPWARD, hi1, lo2);+ double w = rounded_div_impl_double(ROUND_UPWARD, hi1, hi2);+ if (isnan(x)) x = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(y)) y = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(z)) z = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(w)) w = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ return fast_fmax4_double(x, y, z, w);+}++extern double rounded_hw_interval_div_double_down(double lo1, double hi1, double lo2, double hi2)+{+ double x = rounded_div_impl_double(ROUND_DOWNWARD, lo1, lo2);+ double y = rounded_div_impl_double(ROUND_DOWNWARD, lo1, hi2);+ double z = rounded_div_impl_double(ROUND_DOWNWARD, hi1, lo2);+ double w = rounded_div_impl_double(ROUND_DOWNWARD, hi1, hi2);+ if (isnan(x)) x = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(y)) y = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(z)) z = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(w)) w = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ return fast_fmin4_double(x, y, z, w);+}++extern double rounded_hw_interval_div_add_double_up(double lo1, double hi1, double lo2, double hi2, double hi3)+{+ double x = rounded_div_impl_double(ROUND_UPWARD, lo1, lo2);+ double y = rounded_div_impl_double(ROUND_UPWARD, lo1, hi2);+ double z = rounded_div_impl_double(ROUND_UPWARD, hi1, lo2);+ double w = rounded_div_impl_double(ROUND_UPWARD, hi1, hi2);+ if (isnan(x)) x = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(y)) y = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(z)) z = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(w)) w = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ double p = fast_fmax4_double(x, y, z, w);+ return rounded_add_impl_double(ROUND_UPWARD, p, hi3);+}++extern double rounded_hw_interval_div_add_double_down(double lo1, double hi1, double lo2, double hi2, double lo3)+{+ double x = rounded_div_impl_double(ROUND_DOWNWARD, lo1, lo2);+ double y = rounded_div_impl_double(ROUND_DOWNWARD, lo1, hi2);+ double z = rounded_div_impl_double(ROUND_DOWNWARD, hi1, lo2);+ double w = rounded_div_impl_double(ROUND_DOWNWARD, hi1, hi2);+ if (isnan(x)) x = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(y)) y = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(z)) z = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(w)) w = 0.0; /* 0 / 0, +-inf / +-inf -> 0 */+ double p = fast_fmin4_double(x, y, z, w);+ return rounded_add_impl_double(ROUND_DOWNWARD, p, lo3);+}++//+// Vector Operations+//++extern double rounded_hw_vector_sum_double(HsInt mode, HsInt length, HsInt offset, const double *a)+{+ switch (hs_rounding_mode_to_native(mode)) {+ case ROUND_TONEAREST:+ {+ double s = 0.0;+ for (HsInt i = 0; i < length; ++i) {+ s = rounded_add_impl_double(ROUND_TONEAREST, s, a[offset + i]);+ }+ return s;+ }+ case ROUND_DOWNWARD:+ {+ double s = 0.0;+ for (HsInt i = 0; i < length; ++i) {+ s = rounded_add_impl_double(ROUND_DOWNWARD, s, a[offset + i]);+ }+ return s;+ }+ case ROUND_UPWARD:+ {+ double s = 0.0;+ for (HsInt i = 0; i < length; ++i) {+ s = rounded_add_impl_double(ROUND_UPWARD, s, a[offset + i]);+ }+ return s;+ }+ case ROUND_TOWARDZERO:+ {+ double s = 0.0;+ for (HsInt i = 0; i < length; ++i) {+ s = rounded_add_impl_double(ROUND_TOWARDZERO, s, a[offset + i]);+ }+ return s;+ }+ default:+ UNREACHABLE();+ abort();+ }+}++extern void rounded_hw_vector_add_double(HsInt mode, HsInt length, HsInt offsetR, double * restrict result, HsInt offsetA, const double * restrict a, HsInt offsetB, const double * restrict b)+{+ // TODO: Use SIMD+ switch (hs_rounding_mode_to_native(mode)) {+ case ROUND_TONEAREST:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_add_impl_double(ROUND_TONEAREST, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_DOWNWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_add_impl_double(ROUND_DOWNWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_UPWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_add_impl_double(ROUND_UPWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_TOWARDZERO:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_add_impl_double(ROUND_TOWARDZERO, a[offsetA + i], b[offsetB + i]);+ }+ break;+ default:+ UNREACHABLE();+ abort();+ }+}++extern void rounded_hw_vector_sub_double(HsInt mode, HsInt length, HsInt offsetR, double * restrict result, HsInt offsetA, const double * restrict a, HsInt offsetB, const double * restrict b)+{+ // TODO: Use SIMD+ switch (hs_rounding_mode_to_native(mode)) {+ case ROUND_TONEAREST:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sub_impl_double(ROUND_TONEAREST, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_DOWNWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sub_impl_double(ROUND_DOWNWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_UPWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sub_impl_double(ROUND_UPWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_TOWARDZERO:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sub_impl_double(ROUND_TOWARDZERO, a[offsetA + i], b[offsetB + i]);+ }+ break;+ default:+ UNREACHABLE();+ abort();+ }+}++extern void rounded_hw_vector_mul_double(HsInt mode, HsInt length, HsInt offsetR, double * restrict result, HsInt offsetA, const double * restrict a, HsInt offsetB, const double * restrict b)+{+ // TODO: Use SIMD+ switch (hs_rounding_mode_to_native(mode)) {+ case ROUND_TONEAREST:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_mul_impl_double(ROUND_TONEAREST, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_DOWNWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_mul_impl_double(ROUND_DOWNWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_UPWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_mul_impl_double(ROUND_UPWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_TOWARDZERO:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_mul_impl_double(ROUND_TOWARDZERO, a[offsetA + i], b[offsetB + i]);+ }+ break;+ default:+ UNREACHABLE();+ abort();+ }+}++extern void rounded_hw_vector_fma_double(HsInt mode, HsInt length, HsInt offsetR, double * restrict result, HsInt offsetA, const double * restrict a, HsInt offsetB, const double * restrict b, HsInt offsetC, const double * restrict c)+{+ // TODO: Use SIMD+ switch (hs_rounding_mode_to_native(mode)) {+ case ROUND_TONEAREST:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_fma_impl_double(ROUND_TONEAREST, a[offsetA + i], b[offsetB + i], c[offsetC + i]);+ }+ break;+ case ROUND_DOWNWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_fma_impl_double(ROUND_DOWNWARD, a[offsetA + i], b[offsetB + i], c[offsetC + i]);+ }+ break;+ case ROUND_UPWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_fma_impl_double(ROUND_UPWARD, a[offsetA + i], b[offsetB + i], c[offsetC + i]);+ }+ break;+ case ROUND_TOWARDZERO:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_fma_impl_double(ROUND_TOWARDZERO, a[offsetA + i], b[offsetB + i], c[offsetC + i]);+ }+ break;+ default:+ UNREACHABLE();+ abort();+ }+}++extern void rounded_hw_vector_div_double(HsInt mode, HsInt length, HsInt offsetR, double * restrict result, HsInt offsetA, const double * restrict a, HsInt offsetB, const double * restrict b)+{+ // TODO: Use SIMD+ switch (hs_rounding_mode_to_native(mode)) {+ case ROUND_TONEAREST:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_div_impl_double(ROUND_TONEAREST, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_DOWNWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_div_impl_double(ROUND_DOWNWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_UPWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_div_impl_double(ROUND_UPWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_TOWARDZERO:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_div_impl_double(ROUND_TOWARDZERO, a[offsetA + i], b[offsetB + i]);+ }+ break;+ default:+ UNREACHABLE();+ abort();+ }+}++extern void rounded_hw_vector_sqrt_double(HsInt mode, HsInt length, HsInt offsetR, double * restrict result, HsInt offsetA, const double * restrict a)+{+ // TODO: Use SIMD+ switch (hs_rounding_mode_to_native(mode)) {+ case ROUND_TONEAREST:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sqrt_impl_double(ROUND_TONEAREST, a[offsetA + i]);+ }+ break;+ case ROUND_DOWNWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sqrt_impl_double(ROUND_DOWNWARD, a[offsetA + i]);+ }+ break;+ case ROUND_UPWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sqrt_impl_double(ROUND_UPWARD, a[offsetA + i]);+ }+ break;+ case ROUND_TOWARDZERO:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sqrt_impl_double(ROUND_TOWARDZERO, a[offsetA + i]);+ }+ break;+ default:+ UNREACHABLE();+ abort();+ }+}++//+// float+//++static inline ALWAYS_INLINE+float rounded_add_impl_float(native_rounding_mode mode, float a, float b)+{+ float result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fadd.s %0, %1, %2, rne" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_DOWNWARD:+ __asm__("fadd.s %0, %1, %2, rdn" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_UPWARD:+ __asm__("fadd.s %0, %1, %2, rup" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fadd.s %0, %1, %2, rtz" : "=f"(result) : "f"(a), "f"(b));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern float rounded_hw_add_float(HsInt mode, float a, float b)+{ return rounded_add_impl_float(hs_rounding_mode_to_native(mode), a, b); }+extern float rounded_hw_add_float_up(float a, float b)+{ return rounded_add_impl_float(ROUND_UPWARD, a, b); }+extern float rounded_hw_add_float_down(float a, float b)+{ return rounded_add_impl_float(ROUND_DOWNWARD, a, b); }+extern float rounded_hw_add_float_zero(float a, float b)+{ return rounded_add_impl_float(ROUND_TOWARDZERO, a, b); }++static inline ALWAYS_INLINE+float rounded_sub_impl_float(native_rounding_mode mode, float a, float b)+{+ float result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fsub.s %0, %1, %2, rne" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_DOWNWARD:+ __asm__("fsub.s %0, %1, %2, rdn" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_UPWARD:+ __asm__("fsub.s %0, %1, %2, rup" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fsub.s %0, %1, %2, rtz" : "=f"(result) : "f"(a), "f"(b));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern float rounded_hw_sub_float(HsInt mode, float a, float b)+{ return rounded_sub_impl_float(hs_rounding_mode_to_native(mode), a, b); }+extern float rounded_hw_sub_float_up(float a, float b)+{ return rounded_sub_impl_float(ROUND_UPWARD, a, b); }+extern float rounded_hw_sub_float_down(float a, float b)+{ return rounded_sub_impl_float(ROUND_DOWNWARD, a, b); }+extern float rounded_hw_sub_float_zero(float a, float b)+{ return rounded_sub_impl_float(ROUND_TOWARDZERO, a, b); }++static inline ALWAYS_INLINE+float rounded_mul_impl_float(native_rounding_mode mode, float a, float b)+{+ float result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fmul.s %0, %1, %2, rne" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_DOWNWARD:+ __asm__("fmul.s %0, %1, %2, rdn" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_UPWARD:+ __asm__("fmul.s %0, %1, %2, rup" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fmul.s %0, %1, %2, rtz" : "=f"(result) : "f"(a), "f"(b));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern float rounded_hw_mul_float(HsInt mode, float a, float b)+{ return rounded_mul_impl_float(hs_rounding_mode_to_native(mode), a, b); }+extern float rounded_hw_mul_float_up(float a, float b)+{ return rounded_mul_impl_float(ROUND_UPWARD, a, b); }+extern float rounded_hw_mul_float_down(float a, float b)+{ return rounded_mul_impl_float(ROUND_DOWNWARD, a, b); }+extern float rounded_hw_mul_float_zero(float a, float b)+{ return rounded_mul_impl_float(ROUND_TOWARDZERO, a, b); }++static inline ALWAYS_INLINE+float rounded_div_impl_float(native_rounding_mode mode, float a, float b)+{+ float result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fdiv.s %0, %1, %2, rne" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_DOWNWARD:+ __asm__("fdiv.s %0, %1, %2, rdn" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_UPWARD:+ __asm__("fdiv.s %0, %1, %2, rup" : "=f"(result) : "f"(a), "f"(b));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fdiv.s %0, %1, %2, rtz" : "=f"(result) : "f"(a), "f"(b));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern float rounded_hw_div_float(HsInt mode, float a, float b)+{ return rounded_div_impl_float(hs_rounding_mode_to_native(mode), a, b); }+extern float rounded_hw_div_float_up(float a, float b)+{ return rounded_div_impl_float(ROUND_UPWARD, a, b); }+extern float rounded_hw_div_float_down(float a, float b)+{ return rounded_div_impl_float(ROUND_DOWNWARD, a, b); }+extern float rounded_hw_div_float_zero(float a, float b)+{ return rounded_div_impl_float(ROUND_TOWARDZERO, a, b); }++static inline ALWAYS_INLINE+float rounded_sqrt_impl_float(native_rounding_mode mode, float a)+{+ float result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fsqrt.s %0, %1, rne" : "=f"(result) : "f"(a));+ break;+ case ROUND_DOWNWARD:+ __asm__("fsqrt.s %0, %1, rdn" : "=f"(result) : "f"(a));+ break;+ case ROUND_UPWARD:+ __asm__("fsqrt.s %0, %1, rup" : "=f"(result) : "f"(a));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fsqrt.s %0, %1, rtz" : "=f"(result) : "f"(a));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern float rounded_hw_sqrt_float(HsInt mode, float a)+{ return rounded_sqrt_impl_float(hs_rounding_mode_to_native(mode), a); }+extern float rounded_hw_sqrt_float_up(float a)+{ return rounded_sqrt_impl_float(ROUND_UPWARD, a); }+extern float rounded_hw_sqrt_float_down(float a)+{ return rounded_sqrt_impl_float(ROUND_DOWNWARD, a); }+extern float rounded_hw_sqrt_float_zero(float a)+{ return rounded_sqrt_impl_float(ROUND_TOWARDZERO, a); }++static inline ALWAYS_INLINE+float rounded_fma_impl_float(native_rounding_mode mode, float a, float b, float c)+{+ float result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fmadd.s %0, %1, %2, %3, rne" : "=f"(result) : "f"(a), "f"(b), "f"(c));+ break;+ case ROUND_DOWNWARD:+ __asm__("fmadd.s %0, %1, %2, %3, rdn" : "=f"(result) : "f"(a), "f"(b), "f"(c));+ break;+ case ROUND_UPWARD:+ __asm__("fmadd.s %0, %1, %2, %3, rup" : "=f"(result) : "f"(a), "f"(b), "f"(c));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fmadd.s %0, %1, %2, %3, rtz" : "=f"(result) : "f"(a), "f"(b), "f"(c));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern float rounded_hw_fma_float(HsInt mode, float a, float b, float c)+{ return rounded_fma_impl_float(hs_rounding_mode_to_native(mode), a, b, c); }+extern float rounded_hw_fma_float_up(float a, float b, float c)+{ return rounded_fma_impl_float(ROUND_UPWARD, a, b, c); }+extern float rounded_hw_fma_float_down(float a, float b, float c)+{ return rounded_fma_impl_float(ROUND_DOWNWARD, a, b, c); }+extern float rounded_hw_fma_float_zero(float a, float b, float c)+{ return rounded_fma_impl_float(ROUND_TOWARDZERO, a, b, c); }++extern float rounded_hw_fma_if_fast_float(HsInt mode, float a, float b, float c)+{ return rounded_fma_impl_float(hs_rounding_mode_to_native(mode), a, b, c); }+extern float rounded_hw_fma_if_fast_float_up(float a, float b, float c)+{ return rounded_fma_impl_float(ROUND_UPWARD, a, b, c); }+extern float rounded_hw_fma_if_fast_float_down(float a, float b, float c)+{ return rounded_fma_impl_float(ROUND_DOWNWARD, a, b, c); }+extern float rounded_hw_fma_if_fast_float_zero(float a, float b, float c)+{ return rounded_fma_impl_float(ROUND_TOWARDZERO, a, b, c); }++//+// Conversion+//++static inline float rounded_int64_to_float_impl(native_rounding_mode mode, int64_t x)+{+ float result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fcvt.s.l %0, %1, rne" : "=f"(result) : "r"(x));+ break;+ case ROUND_DOWNWARD:+ __asm__("fcvt.s.l %0, %1, rdn" : "=f"(result) : "r"(x));+ break;+ case ROUND_UPWARD:+ __asm__("fcvt.s.l %0, %1, rup" : "=f"(result) : "r"(x));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fcvt.s.l %0, %1, rtz" : "=f"(result) : "r"(x));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern float rounded_hw_int64_to_float(HsInt mode, int64_t x)+{ return rounded_int64_to_float_impl(hs_rounding_mode_to_native(mode), x); }+extern float rounded_hw_int64_to_float_up(int64_t x)+{ return rounded_int64_to_float_impl(ROUND_UPWARD, x); }+extern float rounded_hw_int64_to_float_down(int64_t x)+{ return rounded_int64_to_float_impl(ROUND_DOWNWARD, x); }+extern float rounded_hw_int64_to_float_zero(int64_t x)+{ return rounded_int64_to_float_impl(ROUND_TOWARDZERO, x); }++static inline float rounded_word64_to_float_impl(native_rounding_mode mode, uint64_t x)+{+ float result;+ switch (mode) {+ case ROUND_TONEAREST:+ __asm__("fcvt.s.lu %0, %1, rne" : "=f"(result) : "r"(x));+ break;+ case ROUND_DOWNWARD:+ __asm__("fcvt.s.lu %0, %1, rdn" : "=f"(result) : "r"(x));+ break;+ case ROUND_UPWARD:+ __asm__("fcvt.s.lu %0, %1, rup" : "=f"(result) : "r"(x));+ break;+ case ROUND_TOWARDZERO:+ __asm__("fcvt.s.lu %0, %1, rtz" : "=f"(result) : "r"(x));+ break;+ default:+ UNREACHABLE();+ abort();+ }+ return result;+}+extern float rounded_hw_word64_to_float(HsInt mode, uint64_t x)+{ return rounded_word64_to_float_impl(hs_rounding_mode_to_native(mode), x); }+extern float rounded_hw_word64_to_float_up(uint64_t x)+{ return rounded_word64_to_float_impl(ROUND_UPWARD, x); }+extern float rounded_hw_word64_to_float_down(uint64_t x)+{ return rounded_word64_to_float_impl(ROUND_DOWNWARD, x); }+extern float rounded_hw_word64_to_float_zero(uint64_t x)+{ return rounded_word64_to_float_impl(ROUND_TOWARDZERO, x); }++//+// Interval arithmetic+//++static inline float fast_fmax_float(float x, float y)+{+ // IEEE 754-2019 maximumNumber+ float result;+ __asm__("fmax.s %0, %1, %2" : "=f"(result) : "f"(x), "f"(y));+ return result;+}+static inline float fast_fmax4_float(float x, float y, float z, float w)+{+ return fast_fmax_float(fast_fmax_float(x, y), fast_fmax_float(z, w));+}+static inline float fast_fmin_float(float x, float y)+{+ // IEEE 754-2019 minimumNumber+ float result;+ __asm__("fmin.s %0, %1, %2" : "=f"(result) : "f"(x), "f"(y));+ return result;+}+static inline float fast_fmin4_float(float x, float y, float z, float w)+{+ return fast_fmin_float(fast_fmin_float(x, y), fast_fmin_float(z, w));+}++extern float rounded_hw_interval_mul_float_up(float lo1, float hi1, float lo2, float hi2)+{+ float x = rounded_mul_impl_float(ROUND_UPWARD, lo1, lo2);+ float y = rounded_mul_impl_float(ROUND_UPWARD, lo1, hi2);+ float z = rounded_mul_impl_float(ROUND_UPWARD, hi1, lo2);+ float w = rounded_mul_impl_float(ROUND_UPWARD, hi1, hi2);+ if (isnan(x)) x = 0.0f; /* 0 * inf -> 0 */+ if (isnan(y)) y = 0.0f; /* 0 * inf -> 0 */+ if (isnan(z)) z = 0.0f; /* 0 * inf -> 0 */+ if (isnan(w)) w = 0.0f; /* 0 * inf -> 0 */+ return fast_fmax4_float(x, y, z, w);+}++extern float rounded_hw_interval_mul_float_down(float lo1, float hi1, float lo2, float hi2)+{+ float x = rounded_mul_impl_float(ROUND_DOWNWARD, lo1, lo2);+ float y = rounded_mul_impl_float(ROUND_DOWNWARD, lo1, hi2);+ float z = rounded_mul_impl_float(ROUND_DOWNWARD, hi1, lo2);+ float w = rounded_mul_impl_float(ROUND_DOWNWARD, hi1, hi2);+ if (isnan(x)) x = 0.0f; /* 0 * inf -> 0 */+ if (isnan(y)) y = 0.0f; /* 0 * inf -> 0 */+ if (isnan(z)) z = 0.0f; /* 0 * inf -> 0 */+ if (isnan(w)) w = 0.0f; /* 0 * inf -> 0 */+ return fast_fmin4_float(x, y, z, w);+}++extern float rounded_hw_interval_mul_add_float_up(float lo1, float hi1, float lo2, float hi2, float hi3)+{+ float x = rounded_mul_impl_float(ROUND_UPWARD, lo1, lo2);+ float y = rounded_mul_impl_float(ROUND_UPWARD, lo1, hi2);+ float z = rounded_mul_impl_float(ROUND_UPWARD, hi1, lo2);+ float w = rounded_mul_impl_float(ROUND_UPWARD, hi1, hi2);+ if (isnan(x)) x = 0.0f; /* 0 * inf -> 0 */+ if (isnan(y)) y = 0.0f; /* 0 * inf -> 0 */+ if (isnan(z)) z = 0.0f; /* 0 * inf -> 0 */+ if (isnan(w)) w = 0.0f; /* 0 * inf -> 0 */+ float p = fast_fmax4_float(x, y, z, w);+ return rounded_add_impl_float(ROUND_UPWARD, p, hi3);+}++extern float rounded_hw_interval_mul_add_float_down(float lo1, float hi1, float lo2, float hi2, float lo3)+{+ float x = rounded_mul_impl_float(ROUND_DOWNWARD, lo1, lo2);+ float y = rounded_mul_impl_float(ROUND_DOWNWARD, lo1, hi2);+ float z = rounded_mul_impl_float(ROUND_DOWNWARD, hi1, lo2);+ float w = rounded_mul_impl_float(ROUND_DOWNWARD, hi1, hi2);+ if (isnan(x)) x = 0.0f; /* 0 * inf -> 0 */+ if (isnan(y)) y = 0.0f; /* 0 * inf -> 0 */+ if (isnan(z)) z = 0.0f; /* 0 * inf -> 0 */+ if (isnan(w)) w = 0.0f; /* 0 * inf -> 0 */+ float p = fast_fmin4_float(x, y, z, w);+ return rounded_add_impl_float(ROUND_DOWNWARD, p, lo3);+}++extern float rounded_hw_interval_div_float_up(float lo1, float hi1, float lo2, float hi2)+{+ float x = rounded_div_impl_float(ROUND_UPWARD, lo1, lo2);+ float y = rounded_div_impl_float(ROUND_UPWARD, lo1, hi2);+ float z = rounded_div_impl_float(ROUND_UPWARD, hi1, lo2);+ float w = rounded_div_impl_float(ROUND_UPWARD, hi1, hi2);+ if (isnan(x)) x = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(y)) y = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(z)) z = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(w)) w = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ return fast_fmax4_float(x, y, z, w);+}++extern float rounded_hw_interval_div_float_down(float lo1, float hi1, float lo2, float hi2)+{+ float x = rounded_div_impl_float(ROUND_DOWNWARD, lo1, lo2);+ float y = rounded_div_impl_float(ROUND_DOWNWARD, lo1, hi2);+ float z = rounded_div_impl_float(ROUND_DOWNWARD, hi1, lo2);+ float w = rounded_div_impl_float(ROUND_DOWNWARD, hi1, hi2);+ if (isnan(x)) x = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(y)) y = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(z)) z = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(w)) w = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ return fast_fmin4_float(x, y, z, w);+}++extern float rounded_hw_interval_div_add_float_up(float lo1, float hi1, float lo2, float hi2, float hi3)+{+ float x = rounded_div_impl_float(ROUND_UPWARD, lo1, lo2);+ float y = rounded_div_impl_float(ROUND_UPWARD, lo1, hi2);+ float z = rounded_div_impl_float(ROUND_UPWARD, hi1, lo2);+ float w = rounded_div_impl_float(ROUND_UPWARD, hi1, hi2);+ if (isnan(x)) x = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(y)) y = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(z)) z = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(w)) w = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ float p = fast_fmax4_float(x, y, z, w);+ return rounded_add_impl_float(ROUND_UPWARD, p, hi3);+}++extern float rounded_hw_interval_div_add_float_down(float lo1, float hi1, float lo2, float hi2, float lo3)+{+ float x = rounded_div_impl_float(ROUND_DOWNWARD, lo1, lo2);+ float y = rounded_div_impl_float(ROUND_DOWNWARD, lo1, hi2);+ float z = rounded_div_impl_float(ROUND_DOWNWARD, hi1, lo2);+ float w = rounded_div_impl_float(ROUND_DOWNWARD, hi1, hi2);+ if (isnan(x)) x = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(y)) y = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(z)) z = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ if (isnan(w)) w = 0.0f; /* 0 / 0, +-inf / +-inf -> 0 */+ float p = fast_fmin4_float(x, y, z, w);+ return rounded_add_impl_float(ROUND_DOWNWARD, p, lo3);+}++//+// Vector Operations+//++extern float rounded_hw_vector_sum_float(HsInt mode, HsInt length, HsInt offset, const float *a)+{+ switch (hs_rounding_mode_to_native(mode)) {+ case ROUND_TONEAREST:+ {+ float s = 0.0f;+ for (HsInt i = 0; i < length; ++i) {+ s = rounded_add_impl_float(ROUND_TONEAREST, s, a[offset + i]);+ }+ return s;+ }+ case ROUND_DOWNWARD:+ {+ float s = 0.0f;+ for (HsInt i = 0; i < length; ++i) {+ s = rounded_add_impl_float(ROUND_DOWNWARD, s, a[offset + i]);+ }+ return s;+ }+ case ROUND_UPWARD:+ {+ float s = 0.0f;+ for (HsInt i = 0; i < length; ++i) {+ s = rounded_add_impl_float(ROUND_UPWARD, s, a[offset + i]);+ }+ return s;+ }+ case ROUND_TOWARDZERO:+ {+ float s = 0.0f;+ for (HsInt i = 0; i < length; ++i) {+ s = rounded_add_impl_float(ROUND_TOWARDZERO, s, a[offset + i]);+ }+ return s;+ }+ default:+ UNREACHABLE();+ abort();+ }+}++extern void rounded_hw_vector_add_float(HsInt mode, HsInt length, HsInt offsetR, float * restrict result, HsInt offsetA, const float * restrict a, HsInt offsetB, const float * restrict b)+{+ // TODO: Use SIMD+ switch (hs_rounding_mode_to_native(mode)) {+ case ROUND_TONEAREST:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_add_impl_float(ROUND_TONEAREST, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_DOWNWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_add_impl_float(ROUND_DOWNWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_UPWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_add_impl_float(ROUND_UPWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_TOWARDZERO:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_add_impl_float(ROUND_TOWARDZERO, a[offsetA + i], b[offsetB + i]);+ }+ break;+ default:+ UNREACHABLE();+ abort();+ }+}++extern void rounded_hw_vector_sub_float(HsInt mode, HsInt length, HsInt offsetR, float * restrict result, HsInt offsetA, const float * restrict a, HsInt offsetB, const float * restrict b)+{+ // TODO: Use SIMD+ switch (hs_rounding_mode_to_native(mode)) {+ case ROUND_TONEAREST:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sub_impl_float(ROUND_TONEAREST, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_DOWNWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sub_impl_float(ROUND_DOWNWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_UPWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sub_impl_float(ROUND_UPWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_TOWARDZERO:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sub_impl_float(ROUND_TOWARDZERO, a[offsetA + i], b[offsetB + i]);+ }+ break;+ default:+ UNREACHABLE();+ abort();+ }+}++extern void rounded_hw_vector_mul_float(HsInt mode, HsInt length, HsInt offsetR, float * restrict result, HsInt offsetA, const float * restrict a, HsInt offsetB, const float * restrict b)+{+ // TODO: Use SIMD+ switch (hs_rounding_mode_to_native(mode)) {+ case ROUND_TONEAREST:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_mul_impl_float(ROUND_TONEAREST, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_DOWNWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_mul_impl_float(ROUND_DOWNWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_UPWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_mul_impl_float(ROUND_UPWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_TOWARDZERO:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_mul_impl_float(ROUND_TOWARDZERO, a[offsetA + i], b[offsetB + i]);+ }+ break;+ default:+ UNREACHABLE();+ abort();+ }+}++extern void rounded_hw_vector_fma_float(HsInt mode, HsInt length, HsInt offsetR, float * restrict result, HsInt offsetA, const float * restrict a, HsInt offsetB, const float * restrict b, HsInt offsetC, const float * restrict c)+{+ // TODO: Use SIMD+ switch (hs_rounding_mode_to_native(mode)) {+ case ROUND_TONEAREST:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_fma_impl_float(ROUND_TONEAREST, a[offsetA + i], b[offsetB + i], c[offsetC + i]);+ }+ break;+ case ROUND_DOWNWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_fma_impl_float(ROUND_DOWNWARD, a[offsetA + i], b[offsetB + i], c[offsetC + i]);+ }+ break;+ case ROUND_UPWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_fma_impl_float(ROUND_UPWARD, a[offsetA + i], b[offsetB + i], c[offsetC + i]);+ }+ break;+ case ROUND_TOWARDZERO:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_fma_impl_float(ROUND_TOWARDZERO, a[offsetA + i], b[offsetB + i], c[offsetC + i]);+ }+ break;+ default:+ UNREACHABLE();+ abort();+ }+}++extern void rounded_hw_vector_div_float(HsInt mode, HsInt length, HsInt offsetR, float * restrict result, HsInt offsetA, const float * restrict a, HsInt offsetB, const float * restrict b)+{+ // TODO: Use SIMD+ switch (hs_rounding_mode_to_native(mode)) {+ case ROUND_TONEAREST:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_div_impl_float(ROUND_TONEAREST, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_DOWNWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_div_impl_float(ROUND_DOWNWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_UPWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_div_impl_float(ROUND_UPWARD, a[offsetA + i], b[offsetB + i]);+ }+ break;+ case ROUND_TOWARDZERO:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_div_impl_float(ROUND_TOWARDZERO, a[offsetA + i], b[offsetB + i]);+ }+ break;+ default:+ UNREACHABLE();+ abort();+ }+}++extern void rounded_hw_vector_sqrt_float(HsInt mode, HsInt length, HsInt offsetR, float * restrict result, HsInt offsetA, const float * restrict a)+{+ // TODO: Use SIMD+ switch (hs_rounding_mode_to_native(mode)) {+ case ROUND_TONEAREST:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sqrt_impl_float(ROUND_TONEAREST, a[offsetA + i]);+ }+ break;+ case ROUND_DOWNWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sqrt_impl_float(ROUND_DOWNWARD, a[offsetA + i]);+ }+ break;+ case ROUND_UPWARD:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sqrt_impl_float(ROUND_UPWARD, a[offsetA + i]);+ }+ break;+ case ROUND_TOWARDZERO:+ for (HsInt i = 0; i < length; ++i) {+ result[offsetR + i] = rounded_sqrt_impl_float(ROUND_TOWARDZERO, a[offsetA + i]);+ }+ break;+ default:+ UNREACHABLE();+ abort();+ }+}
cbits/rounded.c view
@@ -30,6 +30,9 @@ #elif defined(__aarch64__) // If we are on AArch64, use the control register. #define USE_AARCH64_FPCR+#elif defined(__riscv)+// If we are on RISC-V, use the static rounding mode.+#define USE_RISCV #else // Otherwise, use C99's fesetround. #define USE_C99@@ -80,11 +83,26 @@ #include <x86intrin.h> typedef unsigned int fp_reg;+static inline ALWAYS_INLINE+fp_reg get_fp_reg(void)+{+ return _mm_getcsr();+}+static inline ALWAYS_INLINE+void restore_fp_reg(fp_reg reg)+{+ _mm_setcsr(reg);+}++#if 1+// Use raw MXCSR values.+// Unlike `fesetround`, we don't need to respect other FP flags (like FTZ),+// so we don't need to do something like `_mm_setcsr((reg & ~(3u << 13)) | (mode << 13))`. typedef unsigned int native_rounding_mode;-static const native_rounding_mode ROUND_TONEAREST = 0;-static const native_rounding_mode ROUND_DOWNWARD = 1;-static const native_rounding_mode ROUND_UPWARD = 2;-static const native_rounding_mode ROUND_TOWARDZERO = 3;+static const native_rounding_mode ROUND_TONEAREST = 0x1f80 | (0 << 13);+static const native_rounding_mode ROUND_DOWNWARD = 0x1f80 | (1 << 13);+static const native_rounding_mode ROUND_UPWARD = 0x1f80 | (2 << 13);+static const native_rounding_mode ROUND_TOWARDZERO = 0x1f80 | (3 << 13); static inline ALWAYS_INLINE native_rounding_mode hs_rounding_mode_to_native(HsInt mode)@@ -93,25 +111,55 @@ * The order of RoundingMode in Numeric.Rounded.Hardware.Internal.Rounding is * chosen so that the conversion here becomes trivial. */- return (native_rounding_mode)mode;+ return 0x1f80 | (mode << 13); } static inline ALWAYS_INLINE-fp_reg get_fp_reg(void)+void set_rounding(fp_reg reg, native_rounding_mode mode) {- return _mm_getcsr();+ (void)reg;+ _mm_setcsr(mode); }++#else+// Idea: Theoretically, `_mm_setcsr(*mode)` can be compiled to one instruction (LDMXCSR).+// Unfortunately, GCC/Clang do not do so (produce additional MOV instructions).+// Disabling for now.++typedef const volatile uint32_t *native_rounding_mode;+static const volatile uint32_t MXCSR_VALUES[4] = {+ 0x1f80 | (0 << 13),+ 0x1f80 | (1 << 13),+ 0x1f80 | (2 << 13),+ 0x1f80 | (3 << 13),+};++static const native_rounding_mode ROUND_TONEAREST = &MXCSR_VALUES[0];+static const native_rounding_mode ROUND_DOWNWARD = &MXCSR_VALUES[1];+static const native_rounding_mode ROUND_UPWARD = &MXCSR_VALUES[2];+static const native_rounding_mode ROUND_TOWARDZERO = &MXCSR_VALUES[3];+ static inline ALWAYS_INLINE-void set_rounding(fp_reg reg, native_rounding_mode mode)+native_rounding_mode hs_rounding_mode_to_native(HsInt mode) {- _mm_setcsr((reg & ~(3u << 13)) | (mode << 13));+ switch (mode) {+ case /* ToNearest */ 0: return ROUND_TONEAREST;+ case /* TowardNegInf */ 1: return ROUND_DOWNWARD;+ case /* TowardInf */ 2: return ROUND_UPWARD;+ case /* TowardZero */ 3: return ROUND_TOWARDZERO;+ default: UNREACHABLE(); return ROUND_TONEAREST;+ } }+ static inline ALWAYS_INLINE-void restore_fp_reg(fp_reg reg)+void set_rounding(fp_reg reg, native_rounding_mode mode) {- _mm_setcsr(reg);+ (void)reg;+ _mm_setcsr(*mode); } +#endif+ static const char backend_name[] = "SSE2"; #elif defined(USE_AARCH64_FPCR)@@ -185,6 +233,32 @@ static const char backend_name[] = "AArch64 FPCR"; +#elif defined(USE_RISCV)++/*+ * RISC-V dynamic rounding mode:+ * ROUND_TONEAREST = 0+ * ROUND_DOWNWARD = 2+ * ROUND_UPWARD = 3+ * ROUND_TOWARDZERO = 1+ * ROUND_TIESTOAWAY = 4+ * We don't use the dynamic rounding mode on RISC-V.+ */++typedef enum {+ /* The order is same as RoundingMode in Numeric.Rounded.Hardware.Internal.Rounding */+ ROUND_TONEAREST = 0,+ ROUND_DOWNWARD,+ ROUND_UPWARD,+ ROUND_TOWARDZERO+} native_rounding_mode;++static inline ALWAYS_INLINE+native_rounding_mode hs_rounding_mode_to_native(HsInt mode)+{ return (native_rounding_mode)mode; }++static const char backend_name[] = "RISC-V";+ #elif defined(USE_C99) #include <fenv.h>@@ -232,6 +306,8 @@ #if defined(USE_AVX512) #include "rounded-avx512.inl"+#elif defined(USE_RISCV)+#include "rounded-riscv.inl" #else #include "rounded-common.inl" #endif
rounded-hw.cabal view
@@ -2,7 +2,7 @@ -- asm-sources is Cabal 3.0 feature name: rounded-hw-version: 0.4.0.2+version: 0.4.0.3 synopsis: Directed rounding for built-in floating types description: Please see the README on GitHub at <https://github.com/minoki/haskell-floating-point/tree/master/rounded-hw#readme> category: Numeric, Math@@ -10,17 +10,18 @@ bug-reports: https://github.com/minoki/haskell-floating-point/issues author: ARATA Mizuki maintainer: minorinoki@gmail.com-copyright: 2020-2024 ARATA Mizuki+copyright: 2020-2026 ARATA Mizuki license: BSD-3-Clause license-file: LICENSE build-type: Custom tested-with:- GHC == 8.6.5, GHC == 8.8.4, GHC == 8.10.7, GHC == 9.0.2, GHC == 9.2.8, GHC == 9.4.8, GHC == 9.6.7, GHC == 9.8.4, GHC == 9.10.3, GHC == 9.12.2+ GHC == 8.6.5, GHC == 8.8.4, GHC == 8.10.7, GHC == 9.0.2, GHC == 9.2.8, GHC == 9.4.8, GHC == 9.6.7, GHC == 9.8.4, GHC == 9.10.3, GHC == 9.12.4, GHC == 9.14.1 extra-source-files: README.md ChangeLog.md cbits/rounded-common.inl cbits/rounded-avx512.inl+ cbits/rounded-riscv.inl cbits/interval-prim-x86_64-sse2.S cbits/interval-prim-x86_64-avx512.S @@ -32,8 +33,8 @@ -- Custom setup is required to allow assembly sources to #include "ghcconfig.h" on GHC < 9.12 custom-setup setup-depends:- Cabal >=3.0 && <3.17- , base >=4.12 && <4.23+ Cabal >=3.0 && <3.19+ , base >=4.12 && <4.24 flag pure-hs description: Disable FFI@@ -68,7 +69,7 @@ common deps build-depends: array >=0.5.2.0 && <0.6- , base >=4.12 && <4.23+ , base >=4.12 && <4.24 , deepseq >=1.4.4.0 && <1.6 , fp-ieee ==0.1.* , primitive >=0.6.1.1 && <0.10@@ -103,7 +104,7 @@ hs-source-dirs: src build-depends:- tagged >=0.8.6 && <0.9+ tagged >=0.8.6 && <0.10 ghc-options: -Wall -- Use FFI when flag(pure-hs) is off if !flag(pure-hs)@@ -180,6 +181,7 @@ KindSignatures MagicHash MultiParamTypeClasses+ MultiWayIf NumericUnderscores RankNTypes ScopedTypeVariables@@ -192,9 +194,9 @@ type: exitcode-stdio-1.0 main-is: doctests.hs build-depends:- doctest >=0.22.2 && <0.25+ doctest >=0.22.2 && <0.26 default-language: Haskell2010- if impl(ghc >= 9.14)+ if impl(ghc >= 10.0) buildable: False test-suite rounded-hw-test@@ -206,6 +208,7 @@ FromIntegerSpec FromRationalSpec IntervalArithmeticSpec+ IntervalArithmeticNESpec RoundedArithmeticSpec ShowFloatSpec Util@@ -214,10 +217,12 @@ test ghc-options: -threaded -rtsopts -with-rtsopts=-N build-depends:- QuickCheck >=2.14.3 && <2.18+ QuickCheck >=2.14.3 && <2.19 , hspec ^>=2.11.7 , random >=1.2.1.1 && <1.4 , rounded-hw+ if !flag(pure-hs)+ cpp-options: -DUSE_FFI if flag(x87-long-double) && (arch(i386) || (arch(x86_64) && !os(windows))) -- Support for 80-bit long double is not good on Win64, so don't test other-modules:@@ -259,3 +264,5 @@ FlexibleContexts QuantifiedConstraints default-language: Haskell2010+ if !flag(pure-hs)+ cpp-options: -DUSE_FFI
src/Numeric/Rounded/Hardware/Backend/Default.hs view
@@ -4,19 +4,22 @@ {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE StandaloneDeriving #-} {-# OPTIONS_GHC -Wno-orphans -Wno-unused-imports #-}++#include "ghcconfig.h"+ module Numeric.Rounded.Hardware.Backend.Default () where import qualified Numeric.Rounded.Hardware.Backend.ViaRational as VR import Numeric.Rounded.Hardware.Internal.Class-#ifdef USE_FFI+#if defined(USE_FFI) import qualified Numeric.Rounded.Hardware.Backend.C as C-#ifdef USE_GHC_PRIM+#if defined(USE_GHC_PRIM) && !defined(UnregisterisedCompiler) import qualified Numeric.Rounded.Hardware.Backend.FastFFI as FastFFI #endif-#ifdef USE_X87_LONG_DOUBLE+#if defined(USE_X87_LONG_DOUBLE) import Numeric.Rounded.Hardware.Backend.X87LongDouble () #endif-#ifdef USE_FLOAT128+#if defined(USE_FLOAT128) import Numeric.Rounded.Hardware.Backend.Float128 () #endif #endif@@ -26,8 +29,8 @@ import Numeric.Floating.IEEE import Unsafe.Coerce -#ifdef USE_FFI-#ifdef USE_GHC_PRIM+#if defined(USE_FFI)+#if defined(USE_GHC_PRIM) && !defined(UnregisterisedCompiler) type FloatImpl = C.CFloat -- TODO: Provide FastFFI.CFloat type DoubleImpl = FastFFI.CDouble #else
src/Numeric/Rounded/Hardware/Backend/ViaRational.hs view
@@ -96,7 +96,7 @@ instance (RealFloat a, RealFloatConstants a) => RoundedSqrt (ViaRational a) where roundedSqrt r (ViaRational x)- | r /= ToNearest && x >= 0 = ViaRational $+ | r /= ToNearest && x >= 0 && not (isInfinite x) = ViaRational $ case compare ((toRational y) ^ (2 :: Int)) (toRational x) of LT | r == TowardInf -> let z = nextUp y in assert (toRational x < (toRational z) ^ (2 :: Int)) z
src/Numeric/Rounded/Hardware/Internal/FloatUtil.hs view
@@ -14,8 +14,9 @@ distanceUlp x y | isInfinite x || isInfinite y || isNaN x || isNaN y = Nothing | otherwise = let m = min (abs x) (abs y)- m' = nextUp m- v = (toRational y - toRational x) / toRational (m' - m)+ ulp | m == 0 = minPositive+ | otherwise = encodeFloat 1 (max (fst (floatRange m)) (exponent m) - floatDigits m) `asTypeOf` m+ v = (toRational y - toRational x) / toRational ulp in if denominator v == 1 then Just (abs (numerator v)) else error "distanceUlp"
src/Numeric/Rounded/Hardware/Interval.hs view
@@ -2,6 +2,7 @@ {-# LANGUAGE DataKinds #-} {-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE MultiWayIf #-} {-# LANGUAGE RankNTypes #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeFamilies #-}@@ -87,10 +88,23 @@ {-# SPECIALIZE minI :: Interval Float -> Interval Float -> Interval Float #-} {-# SPECIALIZE minI :: Interval Double -> Interval Double -> Interval Double #-} +negP :: Num a => Rounded 'TowardInf a -> Rounded 'TowardNegInf a+negP (Rounded x) = Rounded (negate x)+{-# INLINE negP #-}++negN :: Num a => Rounded 'TowardNegInf a -> Rounded 'TowardInf a+negN (Rounded x) = Rounded (negate x)+{-# INLINE negN #-}+ powInt :: (Ord a, Num a, RoundedRing a) => Interval a -> Int -> Interval a-powInt (I a a') n | odd n || 0 <= a = I (a^n) (a'^n)- | a' <= 0 = I ((coerce (abs a'))^n) ((coerce (abs a))^n)- | otherwise = I 0 (max ((coerce (abs a))^n) (a'^n))+powInt (I a a') n+ | n == 0 = I 1 1+ | odd n = if | 0 <= a -> I (a^n) (a'^n)+ | a' <= 0 -> I (negP $ (negN a)^n) (negN $ (negP a')^n)+ | otherwise -> I (negP $ (negN a)^n) (a'^n) -- a < 0 < a'+ | otherwise = if | 0 <= a -> I (a^n) (a'^n)+ | a' <= 0 -> I ((coerce (abs a'))^n) ((coerce (abs a))^n)+ | otherwise -> I 0 (max ((coerce (abs a))^n) (a'^n)) -- a < 0 < a' powInt Empty _ = Empty {-# SPECIALIZE powInt :: Interval Float -> Int -> Interval Float #-} {-# SPECIALIZE powInt :: Interval Double -> Int -> Interval Double #-}
src/Numeric/Rounded/Hardware/Interval/ElementaryFunctions.hs view
@@ -165,7 +165,7 @@ else -- -pi <= x' <= pi, x' <= y' <= 3 * pi let include_minus_1 = minus_half_pi_iv `subset` t' || three_pi_2_iv `subset` t' include_plus_1 = pi_iv / 2 `subset` t' || five_pi_2_iv `subset` t'- u = hull (sinP $ singleton x') $ sinP (if y <= getRounded pi_down then singleton y' else singleton y' - 2 * pi_iv)+ u = hull (sinP $ singleton x') $ sinP (if y' <= getRounded pi_down then singleton y' else singleton y' - 2 * pi_iv) v | include_minus_1 = hull (-1) u | otherwise = u w | include_plus_1 = hull 1 v@@ -197,7 +197,7 @@ else -- -pi <= x' <= pi, x' <= y' <= 3 * pi let include_minus_1 = -pi_iv `subset` t' || pi_iv `subset` t' || three_pi_iv `subset` t' include_plus_1 = 0 `subset` t' || 2 * pi_iv `subset` t'- u = hull (cosP $ singleton x') $ cosP (if y <= getRounded pi_down then singleton y' else singleton y' - 2 * pi_iv)+ u = hull (cosP $ singleton x') $ cosP (if y' <= getRounded pi_down then singleton y' else singleton y' - 2 * pi_iv) v | include_minus_1 = hull (-1) u | otherwise = u w | include_plus_1 = hull 1 v
src/Numeric/Rounded/Hardware/Interval/NonEmpty.hs view
@@ -2,6 +2,7 @@ {-# LANGUAGE DataKinds #-} {-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE MultiWayIf #-} {-# LANGUAGE RankNTypes #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeFamilies #-}@@ -137,10 +138,23 @@ minI (I a a') (I b b') = I (min a b) (min a' b') {-# INLINE minI #-} +negP :: Num a => Rounded 'TowardInf a -> Rounded 'TowardNegInf a+negP (Rounded x) = Rounded (negate x)+{-# INLINE negP #-}++negN :: Num a => Rounded 'TowardNegInf a -> Rounded 'TowardInf a+negN (Rounded x) = Rounded (negate x)+{-# INLINE negN #-}+ powInt :: (Ord a, Num a, RoundedRing a) => Interval a -> Int -> Interval a-powInt (I a a') n | odd n || 0 <= a = I (a^n) (a'^n)- | a' <= 0 = I ((coerce (abs a'))^n) ((coerce (abs a))^n)- | otherwise = I 0 (max ((coerce (abs a))^n) (a'^n))+powInt (I a a') n+ | n == 0 = I 1 1+ | odd n = if | 0 <= a -> I (a^n) (a'^n)+ | a' <= 0 -> I (negP $ (negN a)^n) (negN $ (negP a')^n)+ | otherwise -> I (negP $ (negN a)^n) (a'^n) -- a < 0 < a'+ | otherwise = if | 0 <= a -> I (a^n) (a'^n)+ | a' <= 0 -> I ((coerce (abs a'))^n) ((coerce (abs a))^n)+ | otherwise -> I 0 (max ((coerce (abs a))^n) (a'^n)) -- a < 0 < a' {-# SPECIALIZE powInt :: Interval Float -> Int -> Interval Float #-} {-# SPECIALIZE powInt :: Interval Double -> Int -> Interval Double #-}
+ test/IntervalArithmeticNESpec.hs view
@@ -0,0 +1,44 @@+{-# LANGUAGE ScopedTypeVariables #-}+module IntervalArithmeticNESpec where+import Data.Proxy+import Numeric.Rounded.Hardware.Internal+import Numeric.Rounded.Hardware.Interval.NonEmpty+import Numeric.Rounded.Hardware.Interval.Class (makeInterval, equalAsSet, subset)+import Test.Hspec+import Test.Hspec.QuickCheck (prop)+import Test.QuickCheck++data OrdPair a = OrdPair a a deriving (Eq, Show)++instance (Arbitrary a, Ord a) => Arbitrary (OrdPair a) where+ arbitrary = do x <- arbitrary+ y <- arbitrary+ return $ OrdPair (min x y) (max x y)++verifyImplementation :: forall a. (Arbitrary a, Ord a, Show a, RoundedFractional a, RoundedSqrt a, RealFloatConstants a, RealFloat a) => Proxy a -> Spec+verifyImplementation _ = do+ prop "intervalAdd" $ \(OrdPair (x :: a) y) (OrdPair x' y') ->+ let iv1, iv2 :: Interval a+ iv1 = makeInterval (Rounded x) (Rounded y) + makeInterval (Rounded x') (Rounded y')+ iv2 = makeInterval (Rounded $ roundedAdd TowardNegInf x x') (Rounded $ roundedAdd TowardInf y y')+ in iv1 `equalAsSet` iv2+ prop "intervalSub" $ \(OrdPair (x :: a) y) (OrdPair x' y') ->+ let iv1, iv2 :: Interval a+ iv1 = makeInterval (Rounded x) (Rounded y) - makeInterval (Rounded x') (Rounded y')+ iv2 = makeInterval (Rounded $ roundedSub TowardNegInf x y') (Rounded $ roundedSub TowardInf y x')+ in iv1 `equalAsSet` iv2+ prop "intervalSqrt" $ \(OrdPair (NonNegative (x :: a)) (NonNegative y)) ->+ let iv1, iv2 :: Interval a+ iv1 = sqrt (makeInterval (Rounded x) (Rounded y))+ iv2 = makeInterval (Rounded $ roundedSqrt TowardNegInf x) (Rounded $ roundedSqrt TowardInf y)+ in iv1 `equalAsSet` iv2+ prop "powInt" $ \(OrdPair (x :: a) y) (NonNegative n) ->+ let iv :: Interval a+ iv = makeInterval (Rounded x) (Rounded y)+ iv_n = powInt iv n+ in fromRational ((toRational x)^n) `subset` iv_n .&&. fromRational ((toRational y)^n) `subset` iv_n++spec :: Spec+spec = do+ describe "Double" $ verifyImplementation (Proxy :: Proxy Double)+ describe "Float" $ verifyImplementation (Proxy :: Proxy Float)
test/IntervalArithmeticSpec.hs view
@@ -3,7 +3,7 @@ import Data.Proxy import Numeric.Rounded.Hardware.Internal import Numeric.Rounded.Hardware.Interval-import Numeric.Rounded.Hardware.Interval.Class (makeInterval, equalAsSet)+import Numeric.Rounded.Hardware.Interval.Class (makeInterval, equalAsSet, subset) import Test.Hspec import Test.Hspec.QuickCheck (prop) import Test.QuickCheck@@ -32,6 +32,11 @@ iv1 = sqrt (makeInterval (Rounded x) (Rounded y)) iv2 = makeInterval (Rounded $ roundedSqrt TowardNegInf x) (Rounded $ roundedSqrt TowardInf y) in iv1 `equalAsSet` iv2+ prop "powInt" $ \(OrdPair (x :: a) y) (NonNegative n) ->+ let iv :: Interval a+ iv = makeInterval (Rounded x) (Rounded y)+ iv_n = powInt iv n+ in fromRational ((toRational x)^n) `subset` iv_n .&&. fromRational ((toRational y)^n) `subset` iv_n spec :: Spec spec = do
test/RoundedArithmeticSpec.hs view
@@ -1,10 +1,13 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ScopedTypeVariables #-} module RoundedArithmeticSpec where import Data.Coerce import Data.Proxy+#if defined(USE_FFI) import qualified Numeric.Rounded.Hardware.Backend.C as Backend.C+#endif import Numeric.Rounded.Hardware.Backend.ViaRational import Numeric.Rounded.Hardware.Internal import Test.Hspec@@ -89,6 +92,7 @@ describe "Double default" $ verifyImplementation (Proxy :: Proxy Double) (Proxy :: Proxy Double) describe "Float default" $ verifyImplementation (Proxy :: Proxy Float) (Proxy :: Proxy Float) - -- TODO: Disable when `pure-hs` is on+#if defined(USE_FFI) describe "Double C" $ verifyImplementation (Proxy :: Proxy Double) (Proxy :: Proxy Backend.C.CDouble) describe "Float C" $ verifyImplementation (Proxy :: Proxy Float) (Proxy :: Proxy Backend.C.CFloat)+#endif
test/Spec.hs view
@@ -6,6 +6,7 @@ import qualified FromIntegerSpec import qualified FromRationalSpec import qualified IntervalArithmeticSpec+import qualified IntervalArithmeticNESpec import Numeric.Rounded.Hardware.Backend (backendName) import qualified RoundedArithmeticSpec import qualified ShowFloatSpec@@ -40,6 +41,7 @@ describe "showFloat" ShowFloatSpec.spec describe "rounded arithmetic" RoundedArithmeticSpec.spec describe "interval arithmetic" IntervalArithmeticSpec.spec+ describe "interval arithmetic (non-empty)" IntervalArithmeticNESpec.spec describe "Vector" VectorSpec.spec describe "Constants" ConstantsSpec.spec #ifdef TEST_X87_LONG_DOUBLE