diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -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. diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright ARATA Mizuki (c) 2020-2024 +Copyright ARATA Mizuki (c) 2020-2026 All rights reserved. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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. diff --git a/benchmark/Conversion.hs b/benchmark/Conversion.hs --- a/benchmark/Conversion.hs +++ b/benchmark/Conversion.hs @@ -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) diff --git a/cbits/interval-prim-x86_64-avx512.S b/cbits/interval-prim-x86_64-avx512.S --- a/cbits/interval-prim-x86_64-avx512.S +++ b/cbits/interval-prim-x86_64-avx512.S @@ -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 diff --git a/cbits/interval-prim-x86_64-sse2.S b/cbits/interval-prim-x86_64-sse2.S --- a/cbits/interval-prim-x86_64-sse2.S +++ b/cbits/interval-prim-x86_64-sse2.S @@ -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 diff --git a/cbits/rounded-riscv.inl b/cbits/rounded-riscv.inl new file mode 100644 --- /dev/null +++ b/cbits/rounded-riscv.inl @@ -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(); + } +} diff --git a/cbits/rounded.c b/cbits/rounded.c --- a/cbits/rounded.c +++ b/cbits/rounded.c @@ -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 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 @@ -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 diff --git a/rounded-hw.cabal b/rounded-hw.cabal --- a/rounded-hw.cabal +++ b/rounded-hw.cabal @@ -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 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 diff --git a/src/Numeric/Rounded/Hardware/Backend/Default.hs b/src/Numeric/Rounded/Hardware/Backend/Default.hs --- a/src/Numeric/Rounded/Hardware/Backend/Default.hs +++ b/src/Numeric/Rounded/Hardware/Backend/Default.hs @@ -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 diff --git a/src/Numeric/Rounded/Hardware/Backend/ViaRational.hs b/src/Numeric/Rounded/Hardware/Backend/ViaRational.hs --- a/src/Numeric/Rounded/Hardware/Backend/ViaRational.hs +++ b/src/Numeric/Rounded/Hardware/Backend/ViaRational.hs @@ -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 diff --git a/src/Numeric/Rounded/Hardware/Internal/FloatUtil.hs b/src/Numeric/Rounded/Hardware/Internal/FloatUtil.hs --- a/src/Numeric/Rounded/Hardware/Internal/FloatUtil.hs +++ b/src/Numeric/Rounded/Hardware/Internal/FloatUtil.hs @@ -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" diff --git a/src/Numeric/Rounded/Hardware/Interval.hs b/src/Numeric/Rounded/Hardware/Interval.hs --- a/src/Numeric/Rounded/Hardware/Interval.hs +++ b/src/Numeric/Rounded/Hardware/Interval.hs @@ -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 #-} diff --git a/src/Numeric/Rounded/Hardware/Interval/ElementaryFunctions.hs b/src/Numeric/Rounded/Hardware/Interval/ElementaryFunctions.hs --- a/src/Numeric/Rounded/Hardware/Interval/ElementaryFunctions.hs +++ b/src/Numeric/Rounded/Hardware/Interval/ElementaryFunctions.hs @@ -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 diff --git a/src/Numeric/Rounded/Hardware/Interval/NonEmpty.hs b/src/Numeric/Rounded/Hardware/Interval/NonEmpty.hs --- a/src/Numeric/Rounded/Hardware/Interval/NonEmpty.hs +++ b/src/Numeric/Rounded/Hardware/Interval/NonEmpty.hs @@ -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 #-} diff --git a/test/IntervalArithmeticNESpec.hs b/test/IntervalArithmeticNESpec.hs new file mode 100644 --- /dev/null +++ b/test/IntervalArithmeticNESpec.hs @@ -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) diff --git a/test/IntervalArithmeticSpec.hs b/test/IntervalArithmeticSpec.hs --- a/test/IntervalArithmeticSpec.hs +++ b/test/IntervalArithmeticSpec.hs @@ -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 diff --git a/test/RoundedArithmeticSpec.hs b/test/RoundedArithmeticSpec.hs --- a/test/RoundedArithmeticSpec.hs +++ b/test/RoundedArithmeticSpec.hs @@ -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 diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -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