Showing posts with label floating-point. Show all posts
Showing posts with label floating-point. Show all posts

Support of half precision floating-point numbers in RISC-V ISA: Zfh and Zfhmin

Floating-point support in RISC-V ISA

RISC-V does not mandate the support of any floating-point operations as part of the base ISA (RV32I and RV64I) but single precision and double-precision extensions (e.g. RV64F and RV64D) were among the first to be specified. The F extension added 32 floating-point registers, floating-point arithmetic operations, moves and conversions instructions. The D extension (which requires the F extension) extended this set of operation to also support double precision. There also exists a Q extension for quad-precision. 

In reality, there are multiple F and D extensions: RV32F (to extend the 32-bit base ISA) and RV64F (for the 64-bit base ISA), similarly there are RV32D and RV64D. RV32D has the particularity of adding 64-bit wide floating-point registers when the associated general purpose registers are only 32-bit wide. This flexibility of RISC-V ISA is reviewed in the post RISC-V Register Files.

Initially, smaller floaitng-point formats, such as half precision, were not supported. 

Half precision used to be lagging behind in term of ISA and hardware support. It was only specified has a storage format in the 2008 revision of the IEEE-754 standard (the IEEE standard specifying floating-point formats and operations). The momentum of deep learning and convolutional neural networks has kick-started a renewed interest for small number formats and in particular half precision (among many others). 

Note: In the IEEE-754 2019 revision, half precision is still not defined has a basic format but as an interchange format. Although the standard is more permissive in terms of which formats can admit arithmetic operations.

RISC-V International (the association behind the RISC-V specification(s)) recently ratified two extensions specifying sets of instructions for half-precision support: Zfh and Zfhmin. The later being defined as a subset of the former, we will review Zfh first.

Zfh: full half-precision support

Zfh extends floating-point registers (FRF) to support half precision. It requires the F extension, thus no new register file, nor any register width extension is required since the floating-point registers are already wide enough to contain single precision values.
As other single-precision values in RV32D/RV64D, half precision values are stored in larger FLEN registers using RISC-V NaN boxing scheme (more on that in a future post). 

This new extension defines:
  • Instruction to move data (with or without conversions):
    • flh: load from memory into an F register
    • fsh: store to memory from an F register
    • fmv.h.x, fmv.x.h: bit pattern move between X and F register files (unmodified)
    • fcvt.h.(w/l)[u], fcvt.(w/l)[u].h: conversions  between X and F register files (from/to integer)
  • Arithmetic operations: fadd.h, fsub.h, fmul.h, fsqrt.h, fdiv.h, fmin.h, fmax.h, f(n)madd/f(n)msub.h
  • Floating-point comparisons: fcmp.h
  • Conversions between half precision and other floating-point formats
  • Miscellaneous: fclass.hfsgn(/n/x).h

All arithmetic operations operate on uniform I/Os: all operands are half precision values, and the output is an half precision result. They share their opcode with the corresponding F/D instructions: the fmt fields (bits 25 and 26) value 2b'10 (2) encodes half precision (value 2b'00 for single precision, 2b'01 for double, and 2b'11 for quad but who need such a precision !).


The following diagram represents the data move instructions:



The availability of move to/from and conversion with 64-bit operands is conditioned on the availability of RV64 (general purpose move and integer conversions) and D extension (floating-point conversions).
Although certainly anecdotal, the availability of the Q extension (floating-point quad precision, a.k.a 128-bit format): instructions for conversions from/to quad precision are defined: fcvt.h.q and fcvt.q.h.

Zfhmin: reduced half-precision support

The extension Zfhmin can be seen as a subset of Zfh. It only mandates support for half-precision load and store, bit pattern moves from/to integer register file (no conversions with integer formats) and conversions with other floating-point formats. It represents a total of 6 instructions (extended to 8 with conversion from/to double precision format if the D extension is supported).
    Zfhmin constitutes a reduced set of instruction which can be used by platforms where computing with half-precision values directly is not required but which still require the capability to manipulate them, in particular in memory, before converting them to a larger format for an eventual computation.

Vector support for half precision: Zvfh and Zvfhmin

The RISC-V Vector extension (RVV) version 1.0 specified vector support for single and double precision (SEW=32 and 64 bits). A draft of the specification (link to source) introduces Zvfh which is an extension of all the floating-point instruction to half-precision , including conversions with other formats, and Zvfhmin which is a really reduced subset of operations. As is the case for other floating-point format, half precision is supported by a specific value of the vsew field of the vtype configuration register: vsew=1. This encoding corresponding to a Selected Element Width (SEW) of 16 bits and is identical for both integer and floating point formats.

Zvfhmin only mandates the support of conversion between half and single precision: it extends the support of vfwcvt.f.f.v (widening half-to-float) and vfncvt.f.f.w (narrowing float-to-half) to SEW=16-bit.

Zvfh extends to half-precision all the floating-point vector operations, floating-point reductions, floating-point moves., a brief description of those operations can be found on Part 2 of our series RISC-V Vector Extension in a Nutshell.

Both Zvfhmin and Zvfh mandates the support of single precision element in vectors. On top of this support Zvfh mandates at least Zfhmin on the scalar floating-point side.

Half precision in RVA22 profile

The RISC-V consortium defines profiles. These profiles aim at defining a common set of mandatory extensions and a reduced set of optional extensions which can be used by hardware and software providers to build a compatible ecosystem without having to deal with more specialized ISA extensions. Profile descriptions can be found on RISC-V github.

RVA22 is the most recent profile, it is dedicated for 64-bit application processors

Zfhmin is part of the mandatory extensions of the RVA22 profile, while Zfh is an optional extension (which supersede Zfhmin when selected). This means that all application processors targeting compatibility with the RISC-V ecosystem must have a minimal support for half precision, and than extended support is part of the extended profile. Neither of the vector extensions Zvfh nor Zvfhmin are required in the RVA22 profile.

Conclusion

RISC-V Half-precision support in scalar operations has already been ratified (as extensions Zfh and Zfhmin) and his part of the latest application profile (RVA22). There exist draft specifications for the support of half-precision in vector operations: Zvfh and Zvfhmin. Other formats, such as BFloat16 or more esoteric number formats, should follow in the coming years.
RISC-V is an open community so do not hesitate to sign-in to stay up-to-date and participate to the effort: http://riscv.org.

Reference:

Determining numerical bound for floating-point reciprocal overflow


This article is both an example of pythonsollya use and a small study of floating-point property of the reciprocal function.

Context and objectives

Our goal is to find a bound $b$ such that when computing $\frac{1}{x}$ for any $x \lt b $  an IEEE overflow exception is raised (and the value corresponding to an overflow in the given rounding mode is returned). $b^*$ designs the exact value and $b$ the floating-point number of precision $p$ which verify the above property. Let us note that $b$ is not necessarily equal to $b^*$ rounded in a pre-defined rounding mode.

We will use the following notations: $\circ$ is the rounding operation which can be specialized by a rounding mode, $RNE$ is rounding to nearest (tie-to-even), $RU$ is rounding up (to $+\infty$) and $RD$ is rounding downward (to $-\infty)$, $next(x)$ is the minimal floating-point number strictly greater than $x$ and $prev(x)$ is the minimal floating-point number strictly lower than $x$.

Let us first define an overflow bound:

Defintion 1: $b$, (respectively $b^*$) is defined as the floating-point (resp. exact) overflow bound of the reciprocal function if it is the minimal floating-point (resp. real) number which verifies: $ \circ(\frac{1}{b})$ does not overflow.

Section 7.4.0 of the IEEE-754 standard defines overflow as: "The overflow exception shall be signaled if and only if the destination format’s largest finite number is exceeded in magnitude by what would have been the rounded floating-point result (...) were the exponent range unbounded. " This means that $b$ shoud depend on the rounding mode.

Let us first consider $b_{RNE}$ the bound associated to the round to nearest (tie-to-even) mode and defined as the maximal $x$ such that: $ \frac{1}{x} \ge B_{RNE} = \omega + \frac{1}{2} ulp(\omega) $.
$B_{RNE}$ is the mid-point between the largest floating-point number $\omega$ and the first value defined with an unbounded exponent range which will be rounded to $+\infty$ as per Section 7.4.
Note that the inequality does not need to be strict since $\omega$ has an odd mantissa ($1 ... 1$) if $\frac{1}{x} = B_{RNE}$ it should be rounded to a value greater than $\omega$ because of the tie-to-even rule, thus it will overflow.

$p$ is the format precision (e.g. $p=24$ for binary32 format).
$ b^*_{RNE} = \frac{1}{B_{RNE}} = \frac{1}{\omega + 2^{emax - p}} $
$ b^*_{RU} = \frac{1}{B_{RU}} = \frac{1}{\omega} $

Evaluating bound values

Let us compute (approximation to) $b^*_{RNE}$ and $b^*_{RU}$ using pythonsollya:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# importing pythonsollya
>>> import sollya
# forcing display setting to hexadecimal value
# (easier to read)
>>> sollya.settings.display = hexadecimal
# building 2.0 as a Sollya Object
>>> S2 = sollya.SollyaObject(2)
# building binary32's omega value (greatest floating-point number)
>>> omega = S2**127 * (S2**24 - 1) * S2**-23
>>> omega
0x1.fffffep127
>>> B_RU = omega
# omega + 1/2 ulp(omega)
>>> B_RNE = omega  + S2**103
# rounding 1 / B_RU with 100-bit precision
>>> sollya.round(1 / B_RU, 100, sollya.RD)
0x1.000001000001000001000001p-128
# rounding B_RNE to 100-bit precision
>>> sollya.round(1 / B_RNE, 100, sollya.RD)
0x1.0000008000004000002p-128
# rounding 1 / B_RNE to binary32 precision
>>> sollya.round(1 / B_RNE, sollya.binary32, sollya.RN)
0x1p-128
# rounding 1 / B_RU to binary32 precision
>>> sollya.round(1 / B_RU, sollya.binary32, sollya.RN)
0x1p-128

The first interestings values are given lines 16 and 19: correctly rounded approximation on 100 bits of $b_{RU} = \frac{1}{B_{RU}}$ and $b_{RNE} = \frac{1}{B_{RNE}}$. Let us first notice that those two values are greater than they counterpart rounded to binary32 (visible on lines 22 and 25).
Let us call $c_- = 2^{-128}$ and $c_+ = next(c_-) = 2^{-128} + 2^{-149} $ , we have $c_- \lt b_{RU} \lt c_+ $ and $c_- \lt b_{RNE}\lt c_+ $ which means $ \frac{1}{c_-} \gt \frac{1}{b_{RNE}} \gt \frac{1}{c_+} $ and $ \frac{1}{c_-} \gt \frac{1}{b_{RU}} \gt \frac{1}{c_+} $ and $\forall x , x \le c, \frac{1}{x} \gt \frac{1}{b_{RNE}}$

Claim 1: $ \forall mode, \circ_{mode}(\frac{1}{c_+}) $ is a normal value.
Claim 2: $ \forall mode, \circ_{mode}(\frac{1}{c_-}) $ overflows.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
>>> c_minus = S2**-128
>>> c_plus = S2**-128 + S2**-149
>>> c_minus
0x1p-128
>>> c_plus
0x1.000008p-128
>>> for rnd_mode in [sollya.RN, sollya.RU, sollya.RD, sollya.RZ]:
...     print sollya.round(1 / c_minus, sollya.binary32, rnd_mode), sollya.round(1 / c_plus, sollya.binary32, rnd_mode)
infty 0x1.fffffp127
infty 0x1.fffff2p127
0x1.fffffep127 0x1.fffffp127
0x1.fffffep127 0x1.fffffp127

$c_+ $ is greater than both $b^*_{RU}$ and $b^*_{RNE}$, and $c_-$ is lower (strictly) than both of them. Thus, as there is no floating point number between $c_-$ and $c_+$ we can chose $c_-$ as the overflow bound (included) or $c_+$(excluded). There is no need to distinguish several values of $c$ based on the rounding mode. $\forall \ x \ge c_+$, $\frac{1}{x}$ will not overflow and $\forall \ x \le c_-$, $\frac{1}{x}$ will overflow in any rounding mode.

Let us now consider the remaining rounding modes (rounding down and towards zero). Without loss of generality we can reduce the study to positive overflow and thus to a single rounding mode, let us say rounding down.
$B_{RD} = next(\omega) = 2^{128} $ (assuming unbounded exponent range)


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
>>> B_RD = S2**128
>>> b_RD = S2**-128
>>> sollya.round(1 / B_RD, sollya.binary32, sollya.RD)
0x1p-128
>>> sollya.round(1 / B_RD, sollya.binary32, sollya.RU)
0x1p-128
>>> sollya.round(1 / c_minus, sollya.binary32, sollya.RD)
0x1.fffffep127
>>> sollya.round(1 / c_plus, sollya.binary32, sollya.RD)
0x1.fffffp127

$\frac{1}{B_{RD}}$ is exactly equal to the floating point number $c_-$ and $\frac{1}{c_-}$ will indeed lead to an overflow even in rounded down mode. This overflow will return the $\omega$ value as output.
$\circ_{RD}(\frac{1}{c_+}) = \omega - 2^{104} - 2^{105} - 2^{106} $ which is a normal number strictly lower than omega. Thus we get the same behavior as for $RNE$ and $RU$: $c_-$ is the overflow bound (included) and $c_+$ is the excluded overflow bound $\square$.

References:

  • Sollya: a tool for manipulating arbitrary precision numbers, computing approximation polynomial ... (Sollya's website)
  • pythonsollya: a python wrapper for sollya (pythonsollya's gitlab)
  • Quick introduction to IEEE-754 floating-point arithmetic standard (web)