openpower/transcendentals.mdwn

   1 # Transcendental operations
   2
   3 Summary:
   4
   5 *This proposal extends Power ISA scalar floating point operations to
   6 add IEEE754 transcendental functions (pow, log etc) and trigonometric
   7 functions (sin, cos etc). These functions are also 98% shared with the
   8 Khronos Group OpenCL Extended Instruction Set.*
   9
  10 With thanks to:
  11
  12 * Jacob Lifshay
  13 * Dan Petroski
  14 * Mitch Alsup
  15 * Allen Baum
  16 * Andrew Waterman
  17 * Luis Vitorio Cargnini
  18
  19 [[!toc levels=2]]
  20
  21 See:
  22
  23 * <http://bugs.libre-soc.org/show_bug.cgi?id=127>
  24 * <https://www.khronos.org/registry/spir-v/specs/unified1/OpenCL.ExtendedInstructionSet.100.html>
  25 * [[power_trans_ops]] for opcode listing.
  26
  27 Extension subsets:
  28
  29 * **Zftrans**: standard transcendentals (best suited to 3D)
  30 * **ZftransExt**: extra functions (useful, not generally needed for 3D,
  31   can be synthesised using Ztrans)
  32 * **Ztrigpi**: trig. xxx-pi sinpi cospi tanpi
  33 * **Ztrignpi**: trig non-xxx-pi sin cos tan
  34 * **Zarctrigpi**: arc-trig. a-xxx-pi: atan2pi asinpi acospi
  35 * **Zarctrignpi**: arc-trig. non-a-xxx-pi: atan2, asin, acos
  36 * **Zfhyp**: hyperbolic/inverse-hyperbolic.  sinh, cosh, tanh, asinh,
  37   acosh, atanh (can be synthesised - see below)
  38 * **ZftransAdv**: much more complex to implement in hardware
  39 * **Zfrsqrt**: Reciprocal square-root.
  40
  41 Minimum recommended requirements for 3D: Zftrans, Ztrignpi,
  42 Zarctrignpi, with Ztrigpi and Zarctrigpi as augmentations.
  43
  44 Minimum recommended requirements for Mobile-Embedded 3D:
  45 Ztrignpi, Zftrans, with Ztrigpi as an augmentation.
  46
  47 # TODO:
  48
  49 * Decision on accuracy, moved to [[zfpacc_proposal]]
  50 <http://lists.libre-riscv.org/pipermail/libre-riscv-dev/2019-August/002355.html>
  51 * Errors **MUST** be repeatable.
  52 * How about four Platform Specifications? 3DUNIX, UNIX, 3DEmbedded and Embedded?
  53 <http://lists.libre-riscv.org/pipermail/libre-riscv-dev/2019-August/002361.html>
  54   Accuracy requirements for dual (triple) purpose implementations must
  55   meet the higher standard.
  56 * Reciprocal Square-root is in its own separate extension (Zfrsqrt) as
  57   it is desirable on its own by other implementors.  This to be evaluated.
  58
  59 # Requirements <a name="requirements"></a>
  60
  61 This proposal is designed to meet a wide range of extremely diverse
  62 needs, allowing implementors from all of them to benefit from the tools
  63 and hardware cost reductions associated with common standards adoption
  64 in Power ISA (primarily IEEE754 and Vulkan).
  65
  66 **There are *four* different, disparate platform's needs (two new)**:
  67
  68 * 3D Embedded Platform (new)
  69 * Embedded Platform
  70 * 3D UNIX Platform (new)
  71 * UNIX Platform
  72
  73 **The use-cases are**:
  74
  75 * 3D GPUs
  76 * Numerical Computation
  77 * (Potentially) A.I. / Machine-learning (1)
  78
  79 (1) although approximations suffice in this field, making it more likely
  80 to use a custom extension.  High-end ML would inherently definitely
  81 be excluded.
  82
  83 **The power and die-area requirements vary from**:
  84
  85 * Ultra-low-power (smartwatches where GPU power budgets are in milliwatts)
  86 * Mobile-Embedded (good performance with high efficiency for battery life)
  87 * Desktop Computing
  88 * Server / HPC (2)
  89
  90 (2) Supercomputing is left out of the requirements as it is traditionally
  91 covered by Supercomputer Vectorisation Standards (such as RVV).
  92
  93 **The software requirements are**:
  94
  95 * Full public integration into GNU math libraries (libm)
  96 * Full public integration into well-known Numerical Computation systems (numpy)
  97 * Full public integration into upstream GNU and LLVM Compiler toolchains
  98 * Full public integration into Khronos OpenCL SPIR-V compatible Compilers
  99   seeking public Certification and Endorsement from the Khronos Group
 100   under their Trademarked Certification Programme.
 101
 102 **The "contra"-requirements are**:
 103
 104   Ultra Low Power Embedded platforms (smart watches) are sufficiently
 105   resource constrained that Vectorisation (of any kind) is likely to be
 106   unnecessary and inappropriate.
 107 * The requirements are **not** for the purposes of developing a full custom
 108   proprietary GPU with proprietary firmware driven by *hardware* centric
 109   optimised design decisions as a priority over collaboration.
 110 * A full custom proprietary GPU ASIC Manufacturer *may* benefit from
 111   this proposal however the fact that they typically develop proprietary
 112   software that is not shared with the rest of the community likely to
 113   use this proposal means that they have completely different needs.
 114 * This proposal is for *sharing* of effort in reducing development costs
 115
 116 # Proposed Opcodes vs Khronos OpenCL vs IEEE754-2019<a name="khronos_equiv"></a>
 117
 118 This list shows the (direct) equivalence between proposed opcodes,
 119 their Khronos OpenCL equivalents, and their IEEE754-2019 equivalents.
 120 98% of the opcodes in this proposal that are in the IEEE754-2019 standard
 121 are present in the Khronos Extended Instruction Set.
 122
 123 See
 124 <https://www.khronos.org/registry/spir-v/specs/unified1/OpenCL.ExtendedInstructionSet.100.html>
 125 and <https://ieeexplore.ieee.org/document/8766229>
 126
 127 * Special FP16 opcodes are *not* being proposed, except by indirect / inherent
 128   use of elwidth overrides that is already present in the SVP64 Specification.
 129 * "Native" opcodes are *not* being proposed: implementors will be expected
 130   to use the (equivalent) proposed opcode covering the same function.
 131 * "Fast" opcodes are *not* being proposed, because the Khronos Specification
 132   fast\_length, fast\_normalise and fast\_distance OpenCL opcodes require
 133   vectors (or can be done as scalar operations using other Power ISA
 134   instructions).
 135
 136 The OpenCL FP32 opcodes are **direct** equivalents to the proposed opcodes.
 137 Deviation from conformance with the Khronos Specification - including the
 138 Khronos Specification accuracy requirements - is not an option, as it
 139 results in non-compliance, and the vendor may not use the Trademarked words
 140 "Vulkan" etc. in conjunction with their product.
 141
 142 IEEE754-2019 Table 9.1 lists "additional mathematical operations".
 143 Interestingly the only functions missing when compared to OpenCL are
 144 compound, exp2m1, exp10m1, log2p1, log10p1, pown (integer power) and powr.
 145
 146 |opcode  |OpenCL FP32|OpenCL FP16|OpenCL native|OpenCL fast|IEEE754 |Power ISA |
 147 |------- |-----------|-----------|-------------|-----------|------- |--------- |
 148 |FSIN    |sin        |half\_sin  |native\_sin  |NONE       |sin     |NONE      |
 149 |FCOS    |cos        |half\_cos  |native\_cos  |NONE       |cos     |NONE      |
 150 |FTAN    |tan        |half\_tan  |native\_tan  |NONE       |tan     |NONE      |
 151 |NONE (1)|sincos     |NONE       |NONE         |NONE       |NONE    |NONE      |
 152 |FASIN   |asin       |NONE       |NONE         |NONE       |asin    |NONE      |
 153 |FACOS   |acos       |NONE       |NONE         |NONE       |acos    |NONE      |
 154 |FATAN   |atan       |NONE       |NONE         |NONE       |atan    |NONE      |
 155 |FSINPI  |sinpi      |NONE       |NONE         |NONE       |sinPi   |NONE      |
 156 |FCOSPI  |cospi      |NONE       |NONE         |NONE       |cosPi   |NONE      |
 157 |FTANPI  |tanpi      |NONE       |NONE         |NONE       |tanPi   |NONE      |
 158 |FASINPI |asinpi     |NONE       |NONE         |NONE       |asinPi  |NONE      |
 159 |FACOSPI |acospi     |NONE       |NONE         |NONE       |acosPi  |NONE      |
 160 |FATANPI |atanpi     |NONE       |NONE         |NONE       |atanPi  |NONE      |
 161 |FSINH   |sinh       |NONE       |NONE         |NONE       |sinh    |NONE      |
 162 |FCOSH   |cosh       |NONE       |NONE         |NONE       |cosh    |NONE      |
 163 |FTANH   |tanh       |NONE       |NONE         |NONE       |tanh    |NONE      |
 164 |FASINH  |asinh      |NONE       |NONE         |NONE       |asinh   |NONE      |
 165 |FACOSH  |acosh      |NONE       |NONE         |NONE       |acosh   |NONE      |
 166 |FATANH  |atanh      |NONE       |NONE         |NONE       |atanh   |NONE      |
 167 |FATAN2  |atan2      |NONE       |NONE         |NONE       |atan2   |NONE      |
 168 |FATAN2PI|atan2pi    |NONE       |NONE         |NONE       |atan2pi |NONE      |
 169 |FRSQRT  |rsqrt      |half\_rsqrt|native\_rsqrt|NONE       |rSqrt   |fsqrte, fsqrtes (4)   |
 170 |FCBRT   |cbrt       |NONE       |NONE         |NONE       |NONE (2)|NONE      |
 171 |FEXP2   |exp2       |half\_exp2 |native\_exp2 |NONE       |exp2    |NONE      |
 172 |FLOG2   |log2       |half\_log2 |native\_log2 |NONE       |log2    |NONE      |
 173 |FEXPM1  |expm1      |NONE       |NONE         |NONE       |expm1   |NONE      |
 174 |FLOG1P  |log1p      |NONE       |NONE         |NONE       |logp1   |NONE      |
 175 |FEXP    |exp        |half\_exp  |native\_exp  |NONE       |exp     |NONE      |
 176 |FLOG    |log        |half\_log  |native\_log  |NONE       |log     |NONE      |
 177 |FEXP10  |exp10      |half\_exp10|native\_exp10|NONE       |exp10   |NONE      |
 178 |FLOG10  |log10      |half\_log10|native\_log10|NONE       |log10   |NONE      |
 179 |FPOW    |pow        |NONE       |NONE         |NONE       |pow     |NONE      |
 180 |FPOWN   |pown       |NONE       |NONE         |NONE       |pown    |NONE      |
 181 |FPOWR   |powr       |half\_powr |native\_powr |NONE       |powr    |NONE      |
 182 |FROOTN  |rootn      |NONE       |NONE         |NONE       |rootn   |NONE      |
 183 |FHYPOT  |hypot      |NONE       |NONE         |NONE       |hypot   |NONE      |
 184 |FRECIP  |NONE       |half\_recip|native\_recip|NONE       |NONE (3)|fre, fres (4)        |
 185 |NONE    |NONE       |NONE       |NONE         |NONE       |compound|NONE      |
 186 |NONE    |NONE       |NONE       |NONE         |NONE       |exp2m1  |NONE      |
 187 |NONE    |NONE       |NONE       |NONE         |NONE       |exp10m1 |NONE      |
 188 |NONE    |NONE       |NONE       |NONE         |NONE       |log2p1  |NONE      |
 189 |NONE    |NONE       |NONE       |NONE         |NONE       |log10p1 |NONE      |
 190
 191 Note (1) FSINCOS is macro-op fused (see below).
 192
 193 Note (2) synthesised in IEEE754-2019 as "pown(x, 3)"
 194
 195 Note (3) synthesised in IEEE754-2019 using "1.0 / x"
 196
 197 Note (4) these are estimate opcodes that help accelerate
 198 software emulation
 199
 200 ## List of 2-arg opcodes
 201
 202 | opcode   | Description            | pseudocode                 | Extension   |
 203 | ------   | ----------------       | ----------------           | ----------- |
 204 | FATAN2   | atan2 arc tangent      | rd = atan2(rs2, rs1)       | Zarctrignpi |
 205 | FATAN2PI | atan2 arc tangent / pi | rd = atan2(rs2, rs1) / pi  | Zarctrigpi  |
 206 | FPOW     | x power of y           | rd = pow(rs1, rs2)         | ZftransAdv  |
 207 | FPOWN    | x power of n (n int)   | rd = pow(rs1, rs2)         | ZftransAdv  |
 208 | FPOWR    | x power of y (x +ve)   | rd = exp(rs1 log(rs2))     | ZftransAdv  |
 209 | FROOTN   | x power 1/n (n integer)| rd = pow(rs1, 1/rs2)       | ZftransAdv  |
 210 | FHYPOT   | hypotenuse             | rd = sqrt(rs1^2 + rs2^2)   | ZftransAdv  |
 211
 212 ## List of 1-arg transcendental opcodes
 213
 214 | opcode   | Description            | pseudocode                 | Extension   |
 215 | ------   | ----------------       | ----------------           | ----------- |
 216 | FRSQRT   | Reciprocal Square-root   | rd = sqrt(rs1)          | Zfrsqrt    |
 217 | FCBRT    | Cube Root                | rd = pow(rs1, 1.0 / 3)  | ZftransAdv |
 218 | FRECIP   | Reciprocal               | rd = 1.0 / rs1          | Zftrans    |
 219 | FEXP2    | power-of-2               | rd = pow(2, rs1)        | Zftrans    |
 220 | FLOG2    | log2                     | rd = log(2. rs1)        | Zftrans    |
 221 | FEXPM1   | exponential minus 1      | rd = pow(e, rs1) - 1.0  | ZftransExt |
 222 | FLOG1P   | log plus 1               | rd = log(e, 1 + rs1)    | ZftransExt |
 223 | FEXP     | exponential              | rd = pow(e, rs1)        | ZftransExt |
 224 | FLOG     | natural log (base e)     | rd = log(e, rs1)        | ZftransExt |
 225 | FEXP10   | power-of-10              | rd = pow(10, rs1)       | ZftransExt |
 226 | FLOG10   | log base 10              | rd = log(10, rs1)       | ZftransExt |
 227
 228 ## List of 1-arg trigonometric opcodes
 229
 230 | opcode   | Description            | pseudocode                 | Extension   |
 231 | ------   | ----------------       | ----------------           | ----------- |
 232 | FSIN     | sin (radians)            | rd = sin(rs1)           | Ztrignpi    |
 233 | FCOS     | cos (radians)            | rd = cos(rs1)           | Ztrignpi    |
 234 | FTAN     | tan (radians)            | rd = tan(rs1)           | Ztrignpi    |
 235 | FASIN    | arcsin (radians)         | rd = asin(rs1)          | Zarctrignpi |
 236 | FACOS    | arccos (radians)         | rd = acos(rs1)          | Zarctrignpi |
 237 | FATAN    | arctan (radians)         | rd = atan(rs1)          | Zarctrignpi |
 238 | FSINPI   | sin times pi             | rd = sin(pi * rs1)      | Ztrigpi |
 239 | FCOSPI   | cos times pi             | rd = cos(pi * rs1)      | Ztrigpi |
 240 | FTANPI   | tan times pi             | rd = tan(pi * rs1)      | Ztrigpi |
 241 | FASINPI  | arcsin / pi              | rd = asin(rs1) / pi     | Zarctrigpi |
 242 | FACOSPI  | arccos / pi              | rd = acos(rs1) / pi     | Zarctrigpi |
 243 | FATANPI  | arctan / pi              | rd = atan(rs1) / pi     | Zarctrigpi |
 244 | FSINH    | hyperbolic sin (radians) | rd = sinh(rs1)          | Zfhyp |
 245 | FCOSH    | hyperbolic cos (radians) | rd = cosh(rs1)          | Zfhyp |
 246 | FTANH    | hyperbolic tan (radians) | rd = tanh(rs1)          | Zfhyp |
 247 | FASINH   | inverse hyperbolic sin   | rd = asinh(rs1)         | Zfhyp |
 248 | FACOSH   | inverse hyperbolic cos   | rd = acosh(rs1)         | Zfhyp |
 249 | FATANH   | inverse hyperbolic tan   | rd = atanh(rs1)         | Zfhyp |
 250
 251 [[!inline pages="openpower/power_trans_ops" raw=yes ]]
 252
 253 # Subsets
 254
 255 The full set is based on the Khronos OpenCL opcodes. If implemented
 256 entirely it would be too much for both Embedded and also 3D.
 257
 258 The subsets are organised by hardware complexity, need (3D, HPC), however
 259 due to synthesis producing inaccurate results at the range limits,
 260 the less common subsets are still required for IEEE754 HPC.
 261
 262 MALI Midgard, an embedded / mobile 3D GPU, for example only has the
 263 following opcodes:
 264
 265     E8 - fatan_pt2
 266     F0 - frcp (reciprocal)
 267     F2 - frsqrt (inverse square root, 1/sqrt(x))
 268     F3 - fsqrt (square root)
 269     F4 - fexp2 (2^x)
 270     F5 - flog2
 271     F6 - fsin1pi
 272     F7 - fcos1pi
 273     F9 - fatan_pt1
 274
 275 These in FP32 and FP16 only: no FP64 hardware, at all.
 276
 277 Vivante Embedded/Mobile 3D (etnaviv
 278 <https://github.com/laanwj/etna_viv/blob/master/rnndb/isa.xml>)
 279 only has the following:
 280
 281     sin, cos2pi
 282     cos, sin2pi
 283     log2, exp
 284     sqrt and rsqrt
 285     recip.
 286
 287 It also has fast variants of some of these, as a CSR Mode.
 288
 289 AMD's R600 GPU (R600\_Instruction\_Set\_Architecture.pdf) and the
 290 RDNA ISA (RDNA\_Shader\_ISA\_5August2019.pdf, Table 22, Section 6.3) have:
 291
 292     COS2PI (appx)
 293     EXP2
 294     LOG (IEEE754)
 295     RECIP
 296     RSQRT
 297     SQRT
 298     SIN2PI (appx)
 299
 300 AMD RDNA has F16 and F32 variants of all the above, and also has F64
 301 variants of SQRT, RSQRT and RECIP.  It is interesting that even the
 302 modern high-end AMD GPU does not have TAN or ATAN, where MALI Midgard
 303 does.
 304
 305 Also a general point, that customised optimised hardware targetting
 306 FP32 3D with less accuracy simply can neither be used for IEEE754 nor
 307 for FP64 (except as a starting point for hardware or software driven
 308 Newton Raphson or other iterative method).
 309
 310 Also in cost/area sensitive applications even the extra ROM lookup tables
 311 for certain algorithms may be too costly.
 312
 313 These wildly differing and incompatible driving factors lead to the
 314 subset subdivisions, below.
 315
 316 ## Transcendental Subsets
 317
 318 ### Zftrans
 319
 320 LOG2 EXP2 RECIP RSQRT
 321
 322 Zftrans contains the minimum standard transcendentals best suited to
 323 3D. They are also the minimum subset for synthesising log10, exp10,
 324 exp1m, log1p, the hyperbolic trigonometric functions sinh and so on.
 325
 326 They are therefore considered "base" (essential) transcendentals.
 327
 328 ### ZftransExt
 329
 330 LOG, EXP, EXP10, LOG10, LOGP1, EXP1M
 331
 332 These are extra transcendental functions that are useful, not generally
 333 needed for 3D, however for Numerical Computation they may be useful.
 334
 335 Although they can be synthesised using Ztrans (LOG2 multiplied
 336 by a constant), there is both a performance penalty as well as an
 337 accuracy penalty towards the limits, which for IEEE754 compliance is
 338 unacceptable. In particular, LOG(1+rs1) in hardware may give much better
 339 accuracy at the lower end (very small rs1) than LOG(rs1).
 340
 341 Their forced inclusion would be inappropriate as it would penalise
 342 embedded systems with tight power and area budgets.  However if they
 343 were completely excluded the HPC applications would be penalised on
 344 performance and accuracy.
 345
 346 Therefore they are their own subset extension.
 347
 348 ### Zfhyp
 349
 350 SINH, COSH, TANH, ASINH, ACOSH, ATANH
 351
 352 These are the hyperbolic/inverse-hyperbolic functions. Their use in 3D
 353 is limited.
 354
 355 They can all be synthesised using LOG, SQRT and so on, so depend
 356 on Zftrans.  However, once again, at the limits of the range, IEEE754
 357 compliance becomes impossible, and thus a hardware implementation may
 358 be required.
 359
 360 HPC and high-end GPUs are likely markets for these.
 361
 362 ### ZftransAdv
 363
 364 CBRT, POW, POWN, POWR, ROOTN
 365
 366 These are simply much more complex to implement in hardware, and typically
 367 will only be put into HPC applications.
 368
 369 * **Zfrsqrt**: Reciprocal square-root.
 370
 371 ## Trigonometric subsets
 372
 373 ### Ztrigpi vs Ztrignpi
 374
 375 * **Ztrigpi**: SINPI COSPI TANPI
 376 * **Ztrignpi**: SIN COS TAN
 377
 378 Ztrignpi are the basic trigonometric functions through which all others
 379 could be synthesised, and they are typically the base trigonometrics
 380 provided by GPUs for 3D, warranting their own subset.
 381
 382 In the case of the Ztrigpi subset, these are commonly used in for loops
 383 with a power of two number of subdivisions, and the cost of multiplying
 384 by PI inside each loop (or cumulative addition, resulting in cumulative
 385 errors) is not acceptable.
 386
 387 In for example CORDIC the multiplication by PI may be moved outside of
 388 the hardware algorithm as a loop invariant, with no power or area penalty.
 389
 390 Again, therefore, if SINPI (etc.) were excluded, programmers would be
 391 penalised by being forced to divide by PI in some circumstances. Likewise
 392 if SIN were excluded, programmers would be penaslised by being forced
 393 to *multiply* by PI in some circumstances.
 394
 395 Thus again, a slightly different application of the same general argument
 396 applies to give Ztrignpi and Ztrigpi as subsets.  3D GPUs will almost
 397 certainly provide both.
 398
 399 ### Zarctrigpi and Zarctrignpi
 400
 401 * **Zarctrigpi**: ATAN2PI ASINPI ACOSPI
 402 * **Zarctrignpi**: ATAN2 ACOS ASIN
 403
 404 These are extra trigonometric functions that are useful in some
 405 applications, but even for 3D GPUs, particularly embedded and mobile class
 406 GPUs, they are not so common and so are typically synthesised, there.
 407
 408 Although they can be synthesised using Ztrigpi and Ztrignpi, there is,
 409 once again, both a performance penalty as well as an accuracy penalty
 410 towards the limits, which for IEEE754 compliance is unacceptable, yet
 411 is acceptable for 3D.
 412
 413 Therefore they are their own subset extensions.
 414
 415 # Synthesis, Pseudo-code ops and macro-ops
 416
 417 The pseudo-ops are best left up to the compiler rather than being actual
 418 pseudo-ops, by allocating one scalar FP register for use as a constant
 419 (loop invariant) set to "1.0" at the beginning of a function or other
 420 suitable code block.
 421
 422 * FSINCOS - fused macro-op between FSIN and FCOS (issued in that order).
 423 * FSINCOSPI - fused macro-op between FSINPI and FCOSPI (issued in that order).
 424
 425 FATANPI example pseudo-code:
 426
 427     fmvis ft0, 0x3F800 // upper bits of f32 1.0 (BF16)
 428     fatan2pis FRT, FRA, ft0
 429
 430 Hyperbolic function example (obviates need for Zfhyp except for
 431 high-performance or correctly-rounding):
 432
 433     ASINH( x ) = ln( x + SQRT(x**2+1))
 434
 435 # Evaluation and commentary
 436
 437 Moved to [[discussion]]
 438