1 # Transcendental operations
3 To be updated to OpenPOWER.
7 *This proposal extends OpenPOWER scalar floating point operations to
8 add IEEE754 transcendental functions (pow, log etc) and trigonometric
9 functions (sin, cos etc). These functions are also 98% shared with the
10 Khronos Group OpenCL Extended Instruction Set.*
19 * Luis Vitorio Cargnini
25 * <http://bugs.libre-riscv.org/show_bug.cgi?id=127>
26 * <https://www.khronos.org/registry/spir-v/specs/unified1/OpenCL.ExtendedInstructionSet.100.html>
27 * Discussion: <http://lists.libre-riscv.org/pipermail/libre-riscv-dev/2019-August/002342.html>
28 * [[rv_major_opcode_1010011]] for opcode listing.
29 * [[zfpacc_proposal]] for accuracy settings proposal
33 * **Zftrans**: standard transcendentals (best suited to 3D)
34 * **ZftransExt**: extra functions (useful, not generally needed for 3D,
35 can be synthesised using Ztrans)
36 * **Ztrigpi**: trig. xxx-pi sinpi cospi tanpi
37 * **Ztrignpi**: trig non-xxx-pi sin cos tan
38 * **Zarctrigpi**: arc-trig. a-xxx-pi: atan2pi asinpi acospi
39 * **Zarctrignpi**: arc-trig. non-a-xxx-pi: atan2, asin, acos
40 * **Zfhyp**: hyperbolic/inverse-hyperbolic. sinh, cosh, tanh, asinh,
41 acosh, atanh (can be synthesised - see below)
42 * **ZftransAdv**: much more complex to implement in hardware
43 * **Zfrsqrt**: Reciprocal square-root.
45 Minimum recommended requirements for 3D: Zftrans, Ztrignpi,
46 Zarctrignpi, with Ztrigpi and Zarctrigpi as augmentations.
48 Minimum recommended requirements for Mobile-Embedded 3D: Ztrignpi, Zftrans, with Ztrigpi as an augmentation.
52 * Decision on accuracy, moved to [[zfpacc_proposal]]
53 <http://lists.libre-riscv.org/pipermail/libre-riscv-dev/2019-August/002355.html>
54 * Errors **MUST** be repeatable.
55 * How about four Platform Specifications? 3DUNIX, UNIX, 3DEmbedded and Embedded?
56 <http://lists.libre-riscv.org/pipermail/libre-riscv-dev/2019-August/002361.html>
57 Accuracy requirements for dual (triple) purpose implementations must
58 meet the higher standard.
59 * Reciprocal Square-root is in its own separate extension (Zfrsqrt) as
60 it is desirable on its own by other implementors. This to be evaluated.
62 # Requirements <a name="requirements"></a>
64 This proposal is designed to meet a wide range of extremely diverse needs,
65 allowing implementors from all of them to benefit from the tools and hardware
66 cost reductions associated with common standards adoption in RISC-V (primarily IEEE754 and Vulkan).
68 **There are *four* different, disparate platform's needs (two new)**:
70 * 3D Embedded Platform (new)
72 * 3D UNIX Platform (new)
75 **The use-cases are**:
78 * Numerical Computation
79 * (Potentially) A.I. / Machine-learning (1)
81 (1) although approximations suffice in this field, making it more likely
82 to use a custom extension. High-end ML would inherently definitely
85 **The power and die-area requirements vary from**:
87 * Ultra-low-power (smartwatches where GPU power budgets are in milliwatts)
88 * Mobile-Embedded (good performance with high efficiency for battery life)
92 (2) Supercomputing is left out of the requirements as it is traditionally
93 covered by Supercomputer Vectorisation Standards (such as RVV).
95 **The software requirements are**:
97 * Full public integration into GNU math libraries (libm)
98 * Full public integration into well-known Numerical Computation systems (numpy)
99 * Full public integration into upstream GNU and LLVM Compiler toolchains
100 * Full public integration into Khronos OpenCL SPIR-V compatible Compilers
101 seeking public Certification and Endorsement from the Khronos Group
102 under their Trademarked Certification Programme.
104 **The "contra"-requirements are**:
106 Ultra Low Power Embedded platforms (smart watches) are sufficiently
107 resource constrained that Vectorisation (of any kind) is likely to be
108 unnecessary and inappropriate.
109 * The requirements are **not** for the purposes of developing a full custom
110 proprietary GPU with proprietary firmware driven by *hardware* centric
111 optimised design decisions as a priority over collaboration.
112 * A full custom proprietary GPU ASIC Manufacturer *may* benefit from
113 this proposal however the fact that they typically develop proprietary
114 software that is not shared with the rest of the community likely to
115 use this proposal means that they have completely different needs.
116 * This proposal is for *sharing* of effort in reducing development costs
118 # Proposed Opcodes vs Khronos OpenCL vs IEEE754-2019<a name="khronos_equiv"></a>
120 This list shows the (direct) equivalence between proposed opcodes,
121 their Khronos OpenCL equivalents, and their IEEE754-2019 equivalents.
122 98% of the opcodes in this proposal that are in the IEEE754-2019 standard
123 are present in the Khronos Extended Instruction Set.
125 For RISCV opcode encodings see [[rv_major_opcode_1010011]]
126 **TODO** replace with OpenPOWER
129 <https://www.khronos.org/registry/spir-v/specs/unified1/OpenCL.ExtendedInstructionSet.100.html>
130 and <https://ieeexplore.ieee.org/document/8766229>
132 * Special FP16 opcodes are *not* being proposed, except by indirect / inherent
133 use of the "fmt" field that is already present in the RISC-V Specification.
134 * "Native" opcodes are *not* being proposed: implementors will be expected
135 to use the (equivalent) proposed opcode covering the same function.
136 * "Fast" opcodes are *not* being proposed, because the Khronos Specification
137 fast\_length, fast\_normalise and fast\_distance OpenCL opcodes require
138 vectors (or can be done as scalar operations using other RISC-V instructions).
140 The OpenCL FP32 opcodes are **direct** equivalents to the proposed opcodes.
141 Deviation from conformance with the Khronos Specification - including the
142 Khronos Specification accuracy requirements - is not an option, as it
143 results in non-compliance, and the vendor may not use the Trademarked words
144 "Vulkan" etc. in conjunction with their product.
146 IEEE754-2019 Table 9.1 lists "additional mathematical operations".
147 Interestingly the only functions missing when compared to OpenCL are
148 compound, exp2m1, exp10m1, log2p1, log10p1, pown (integer power) and powr.
151 opcode | OpenCL FP32 | OpenCL FP16 | OpenCL native | OpenCL fast | IEEE754 |
152 FSIN | sin | half\_sin | native\_sin | NONE | sin |
153 FCOS | cos | half\_cos | native\_cos | NONE | cos |
154 FTAN | tan | half\_tan | native\_tan | NONE | tan |
155 NONE (1) | sincos | NONE | NONE | NONE | NONE |
156 FASIN | asin | NONE | NONE | NONE | asin |
157 FACOS | acos | NONE | NONE | NONE | acos |
158 FATAN | atan | NONE | NONE | NONE | atan |
159 FSINPI | sinpi | NONE | NONE | NONE | sinPi |
160 FCOSPI | cospi | NONE | NONE | NONE | cosPi |
161 FTANPI | tanpi | NONE | NONE | NONE | tanPi |
162 FASINPI | asinpi | NONE | NONE | NONE | asinPi |
163 FACOSPI | acospi | NONE | NONE | NONE | acosPi |
164 FATANPI | atanpi | NONE | NONE | NONE | atanPi |
165 FSINH | sinh | NONE | NONE | NONE | sinh |
166 FCOSH | cosh | NONE | NONE | NONE | cosh |
167 FTANH | tanh | NONE | NONE | NONE | tanh |
168 FASINH | asinh | NONE | NONE | NONE | asinh |
169 FACOSH | acosh | NONE | NONE | NONE | acosh |
170 FATANH | atanh | NONE | NONE | NONE | atanh |
171 FATAN2 | atan2 | NONE | NONE | NONE | atan2 |
172 FATAN2PI | atan2pi | NONE | NONE | NONE | atan2pi |
173 FRSQRT | rsqrt | half\_rsqrt | native\_rsqrt | NONE | rSqrt |
174 FCBRT | cbrt | NONE | NONE | NONE | NONE (2) |
175 FEXP2 | exp2 | half\_exp2 | native\_exp2 | NONE | exp2 |
176 FLOG2 | log2 | half\_log2 | native\_log2 | NONE | log2 |
177 FEXPM1 | expm1 | NONE | NONE | NONE | expm1 |
178 FLOG1P | log1p | NONE | NONE | NONE | logp1 |
179 FEXP | exp | half\_exp | native\_exp | NONE | exp |
180 FLOG | log | half\_log | native\_log | NONE | log |
181 FEXP10 | exp10 | half\_exp10 | native\_exp10 | NONE | exp10 |
182 FLOG10 | log10 | half\_log10 | native\_log10 | NONE | log10 |
183 FPOW | pow | NONE | NONE | NONE | pow |
184 FPOWN | pown | NONE | NONE | NONE | pown |
185 FPOWR | powr | half\_powr | native\_powr | NONE | powr |
186 FROOTN | rootn | NONE | NONE | NONE | rootn |
187 FHYPOT | hypot | NONE | NONE | NONE | hypot |
188 FRECIP | NONE | half\_recip | native\_recip | NONE | NONE (3) |
189 NONE | NONE | NONE | NONE | NONE | compound |
190 NONE | NONE | NONE | NONE | NONE | exp2m1 |
191 NONE | NONE | NONE | NONE | NONE | exp10m1 |
192 NONE | NONE | NONE | NONE | NONE | log2p1 |
193 NONE | NONE | NONE | NONE | NONE | log10p1 |
196 Note (1) FSINCOS is macro-op fused (see below).
198 Note (2) synthesised in IEEE754-2019 as "pown(x, 3)"
200 Note (3) synthesised in IEEE754-2019 using "1.0 / x"
202 ## List of 2-arg opcodes
205 opcode | Description | pseudocode | Extension |
206 FATAN2 | atan2 arc tangent | rd = atan2(rs2, rs1) | Zarctrignpi |
207 FATAN2PI | atan2 arc tangent / pi | rd = atan2(rs2, rs1) / pi | Zarctrigpi |
208 FPOW | x power of y | rd = pow(rs1, rs2) | ZftransAdv |
209 FPOWN | x power of n (n int) | rd = pow(rs1, rs2) | ZftransAdv |
210 FPOWR | x power of y (x +ve) | rd = exp(rs1 log(rs2)) | ZftransAdv |
211 FROOTN | x power 1/n (n integer)| rd = pow(rs1, 1/rs2) | ZftransAdv |
212 FHYPOT | hypotenuse | rd = sqrt(rs1^2 + rs2^2) | ZftransAdv |
215 ## List of 1-arg transcendental opcodes
218 opcode | Description | pseudocode | Extension |
219 FRSQRT | Reciprocal Square-root | rd = sqrt(rs1) | Zfrsqrt |
220 FCBRT | Cube Root | rd = pow(rs1, 1.0 / 3) | ZftransAdv |
221 FRECIP | Reciprocal | rd = 1.0 / rs1 | Zftrans |
222 FEXP2 | power-of-2 | rd = pow(2, rs1) | Zftrans |
223 FLOG2 | log2 | rd = log(2. rs1) | Zftrans |
224 FEXPM1 | exponential minus 1 | rd = pow(e, rs1) - 1.0 | ZftransExt |
225 FLOG1P | log plus 1 | rd = log(e, 1 + rs1) | ZftransExt |
226 FEXP | exponential | rd = pow(e, rs1) | ZftransExt |
227 FLOG | natural log (base e) | rd = log(e, rs1) | ZftransExt |
228 FEXP10 | power-of-10 | rd = pow(10, rs1) | ZftransExt |
229 FLOG10 | log base 10 | rd = log(10, rs1) | ZftransExt |
232 ## List of 1-arg trigonometric opcodes
235 opcode | Description | pseudo-code | Extension |
236 FSIN | sin (radians) | rd = sin(rs1) | Ztrignpi |
237 FCOS | cos (radians) | rd = cos(rs1) | Ztrignpi |
238 FTAN | tan (radians) | rd = tan(rs1) | Ztrignpi |
239 FASIN | arcsin (radians) | rd = asin(rs1) | Zarctrignpi |
240 FACOS | arccos (radians) | rd = acos(rs1) | Zarctrignpi |
241 FATAN | arctan (radians) | rd = atan(rs1) | Zarctrignpi |
242 FSINPI | sin times pi | rd = sin(pi * rs1) | Ztrigpi |
243 FCOSPI | cos times pi | rd = cos(pi * rs1) | Ztrigpi |
244 FTANPI | tan times pi | rd = tan(pi * rs1) | Ztrigpi |
245 FASINPI | arcsin / pi | rd = asin(rs1) / pi | Zarctrigpi |
246 FACOSPI | arccos / pi | rd = acos(rs1) / pi | Zarctrigpi |
247 FATANPI | arctan / pi | rd = atan(rs1) / pi | Zarctrigpi |
248 FSINH | hyperbolic sin (radians) | rd = sinh(rs1) | Zfhyp |
249 FCOSH | hyperbolic cos (radians) | rd = cosh(rs1) | Zfhyp |
250 FTANH | hyperbolic tan (radians) | rd = tanh(rs1) | Zfhyp |
251 FASINH | inverse hyperbolic sin | rd = asinh(rs1) | Zfhyp |
252 FACOSH | inverse hyperbolic cos | rd = acosh(rs1) | Zfhyp |
253 FATANH | inverse hyperbolic tan | rd = atanh(rs1) | Zfhyp |
258 The full set is based on the Khronos OpenCL opcodes. If implemented
259 entirely it would be too much for both Embedded and also 3D.
261 The subsets are organised by hardware complexity, need (3D, HPC), however
262 due to synthesis producing inaccurate results at the range limits,
263 the less common subsets are still required for IEEE754 HPC.
265 MALI Midgard, an embedded / mobile 3D GPU, for example only has the
269 F0 - frcp (reciprocal)
270 F2 - frsqrt (inverse square root, 1/sqrt(x))
271 F3 - fsqrt (square root)
278 These in FP32 and FP16 only: no FP32 hardware, at all.
280 Vivante Embedded/Mobile 3D (etnaviv
281 <https://github.com/laanwj/etna_viv/blob/master/rnndb/isa.xml>)
282 only has the following:
290 It also has fast variants of some of these, as a CSR Mode.
292 AMD's R600 GPU (R600\_Instruction\_Set\_Architecture.pdf) and the
293 RDNA ISA (RDNA\_Shader\_ISA\_5August2019.pdf, Table 22, Section 6.3) have:
303 AMD RDNA has F16 and F32 variants of all the above, and also has F64
304 variants of SQRT, RSQRT and RECIP. It is interesting that even the
305 modern high-end AMD GPU does not have TAN or ATAN, where MALI Midgard
308 Also a general point, that customised optimised hardware targetting
309 FP32 3D with less accuracy simply can neither be used for IEEE754 nor
310 for FP64 (except as a starting point for hardware or software driven
311 Newton Raphson or other iterative method).
313 Also in cost/area sensitive applications even the extra ROM lookup tables
314 for certain algorithms may be too costly.
316 These wildly differing and incompatible driving factors lead to the
317 subset subdivisions, below.
319 ## Transcendental Subsets
323 LOG2 EXP2 RECIP RSQRT
325 Zftrans contains the minimum standard transcendentals best suited to
326 3D. They are also the minimum subset for synthesising log10, exp10,
327 exp1m, log1p, the hyperbolic trigonometric functions sinh and so on.
329 They are therefore considered "base" (essential) transcendentals.
333 LOG, EXP, EXP10, LOG10, LOGP1, EXP1M
335 These are extra transcendental functions that are useful, not generally
336 needed for 3D, however for Numerical Computation they may be useful.
338 Although they can be synthesised using Ztrans (LOG2 multiplied
339 by a constant), there is both a performance penalty as well as an
340 accuracy penalty towards the limits, which for IEEE754 compliance is
341 unacceptable. In particular, LOG(1+rs1) in hardware may give much better
342 accuracy at the lower end (very small rs1) than LOG(rs1).
344 Their forced inclusion would be inappropriate as it would penalise
345 embedded systems with tight power and area budgets. However if they
346 were completely excluded the HPC applications would be penalised on
347 performance and accuracy.
349 Therefore they are their own subset extension.
353 SINH, COSH, TANH, ASINH, ACOSH, ATANH
355 These are the hyperbolic/inverse-hyperbolic functions. Their use in 3D
358 They can all be synthesised using LOG, SQRT and so on, so depend
359 on Zftrans. However, once again, at the limits of the range, IEEE754
360 compliance becomes impossible, and thus a hardware implementation may
363 HPC and high-end GPUs are likely markets for these.
367 CBRT, POW, POWN, POWR, ROOTN
369 These are simply much more complex to implement in hardware, and typically
370 will only be put into HPC applications.
372 * **Zfrsqrt**: Reciprocal square-root.
374 ## Trigonometric subsets
376 ### Ztrigpi vs Ztrignpi
378 * **Ztrigpi**: SINPI COSPI TANPI
379 * **Ztrignpi**: SIN COS TAN
381 Ztrignpi are the basic trigonometric functions through which all others
382 could be synthesised, and they are typically the base trigonometrics
383 provided by GPUs for 3D, warranting their own subset.
385 In the case of the Ztrigpi subset, these are commonly used in for loops
386 with a power of two number of subdivisions, and the cost of multiplying
387 by PI inside each loop (or cumulative addition, resulting in cumulative
388 errors) is not acceptable.
390 In for example CORDIC the multiplication by PI may be moved outside of
391 the hardware algorithm as a loop invariant, with no power or area penalty.
393 Again, therefore, if SINPI (etc.) were excluded, programmers would be
394 penalised by being forced to divide by PI in some circumstances. Likewise
395 if SIN were excluded, programmers would be penaslised by being forced
396 to *multiply* by PI in some circumstances.
398 Thus again, a slightly different application of the same general argument
399 applies to give Ztrignpi and Ztrigpi as subsets. 3D GPUs will almost
400 certainly provide both.
402 ### Zarctrigpi and Zarctrignpi
404 * **Zarctrigpi**: ATAN2PI ASINPI ACOSPI
405 * **Zarctrignpi**: ATAN2 ACOS ASIN
407 These are extra trigonometric functions that are useful in some
408 applications, but even for 3D GPUs, particularly embedded and mobile class
409 GPUs, they are not so common and so are typically synthesised, there.
411 Although they can be synthesised using Ztrigpi and Ztrignpi, there is,
412 once again, both a performance penalty as well as an accuracy penalty
413 towards the limits, which for IEEE754 compliance is unacceptable, yet
414 is acceptable for 3D.
416 Therefore they are their own subset extensions.
418 # Synthesis, Pseudo-code ops and macro-ops
420 The pseudo-ops are best left up to the compiler rather than being actual
421 pseudo-ops, by allocating one scalar FP register for use as a constant
422 (loop invariant) set to "1.0" at the beginning of a function or other
425 * FSINCOS - fused macro-op between FSIN and FCOS (issued in that order).
426 * FSINCOSPI - fused macro-op between FSINPI and FCOSPI (issued in that order).
428 FATANPI example pseudo-code:
430 lui t0, 0x3F800 // upper bits of f32 1.0
432 fatan2pi.s rd, rs1, ft0
434 Hyperbolic function example (obviates need for Zfhyp except for
435 high-performance or correctly-rounding):
437 ASINH( x ) = ln( x + SQRT(x**2+1))
439 # Evaluation and commentary
441 This section will move later to discussion.
445 Used to be an alias. Some implementors may wish to implement divide as
448 Others may have shared hardware for recip and divide, others may not.
450 To avoid penalising one implementor over another, recip stays.
452 ## To evaluate: should LOG be replaced with LOG1P (and EXP with EXPM1)?
454 RISC principle says "exclude LOG because it's covered by LOGP1 plus an ADD".
455 Research needed to ensure that implementors are not compromised by such
457 <http://lists.libre-riscv.org/pipermail/libre-riscv-dev/2019-August/002358.html>
459 > > correctly-rounded LOG will return different results than LOGP1 and ADD.
460 > > Likewise for EXP and EXPM1
462 > ok, they stay in as real opcodes, then.
464 ## ATAN / ATAN2 commentary
466 Discussion starts here:
467 <http://lists.libre-riscv.org/pipermail/libre-riscv-dev/2019-August/002470.html>
471 would like to point out that the general implementations of ATAN2 do a
472 bunch of special case checks and then simply call ATAN.
474 double ATAN2( double y, double x )
475 { // IEEE 754-2008 quality ATAN2
478 if( ISNAN( x ) ) return x;
479 if( ISNAN( y ) ) return y;
481 // deal with infinities
482 if( x == +∞ && |y|== +∞ ) return copysign( π/4, y );
483 if( x == +∞ ) return copysign( 0.0, y );
484 if( x == -∞ && |y|== +∞ ) return copysign( 3π/4, y );
485 if( x == -∞ ) return copysign( π, y );
486 if( |y|== +∞ ) return copysign( π/2, y );
488 // deal with signed zeros
489 if( x == 0.0 && y != 0.0 ) return copysign( π/2, y );
490 if( x >=+0.0 && y == 0.0 ) return copysign( 0.0, y );
491 if( x <=-0.0 && y == 0.0 ) return copysign( π, y );
493 // calculate ATAN2 textbook style
494 if( x > 0.0 ) return ATAN( |y / x| );
495 if( x < 0.0 ) return π - ATAN( |y / x| );
499 Yet the proposed encoding makes ATAN2 the primitive and has ATAN invent
500 a constant and then call/use ATAN2.
502 When one considers an implementation of ATAN, one must consider several
503 ranges of evaluation::
505 x [ -∞, -1.0]:: ATAN( x ) = -π/2 + ATAN( 1/x );
506 x (-1.0, +1.0]:: ATAN( x ) = + ATAN( x );
507 x [ 1.0, +∞]:: ATAN( x ) = +π/2 - ATAN( 1/x );
509 I should point out that the add/sub of π/2 can not lose significance
510 since the result of ATAN(1/x) is bounded 0..π/2
512 The bottom line is that I think you are choosing to make too many of
513 these into OpCodes, making the hardware function/calculation unit (and
514 sequencer) more complicated that necessary.
516 --------------------------------------------------------
518 We therefore I think have a case for bringing back ATAN and including ATAN2.
520 The reason is that whilst a microcode-like GPU-centric platform would
521 do ATAN2 in terms of ATAN, a UNIX-centric platform would do it the other
524 (that is the hypothesis, to be evaluated for correctness. feedback requested).
526 This because we cannot compromise or prioritise one platfrom's
527 speed/accuracy over another. That is not reasonable or desirable, to
528 penalise one implementor over another.
530 Thus, all implementors, to keep interoperability, must both have both
531 opcodes and may choose, at the architectural and routing level, which
532 one to implement in terms of the other.
534 Allowing implementors to choose to add either opcode and let traps sort it
535 out leaves an uncertainty in the software developer's mind: they cannot
536 trust the hardware, available from many vendors, to be performant right
543 I might suggest that if there were a way for a calculation to be performed
544 and the result of that calculation chained to a subsequent calculation
545 such that the precision of the result-becomes-operand is wider than
546 what will fit in a register, then you can dramatically reduce the count
547 of instructions in this category while retaining
553 can be calculated as::
557 Where 1/y has about 26-to-32 bits of fraction. No, it's not IEEE 754-2008
558 accurate, but GPUs want speed and
560 1/y is fully pipelined (F32) while x/y cannot be (at reasonable area). It
561 is also not "that inaccurate" displaying 0.625-to-0.52 ULP.
563 Given that one has the ability to carry (and process) more fraction bits,
564 one can then do high precision multiplies of π or other transcendental
567 And GPUs have been doing this almost since the dawn of 3D.
569 // calculate ATAN2 high performance style
570 // Note: at this point x != y
574 if( y < 0.0 && |y| < |x| ) return - π/2 - ATAN( x / y );
575 if( y < 0.0 && |y| > |x| ) return + ATAN( y / x );
576 if( y > 0.0 && |y| < |x| ) return + ATAN( y / x );
577 if( y > 0.0 && |y| > |x| ) return + π/2 - ATAN( x / y );
581 if( y < 0.0 && |y| < |x| ) return + π/2 + ATAN( x / y );
582 if( y < 0.0 && |y| > |x| ) return + π - ATAN( y / x );
583 if( y > 0.0 && |y| < |x| ) return + π - ATAN( y / x );
584 if( y > 0.0 && |y| > |x| ) return +3π/2 + ATAN( x / y );
587 This way the adds and subtracts from the constant are not in a precision