openpower/transcendentals.mdwn

   1 # DRAFT Scalar Transcendentals
   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 / Supercomputing
  89
  90 **The software requirements are**:
  91
  92 * Full public integration into GNU math libraries (libm)
  93 * Full public integration into well-known Numerical Computation systems (numpy)
  94 * Full public integration into upstream GNU and LLVM Compiler toolchains
  95 * Full public integration into Khronos OpenCL SPIR-V compatible Compilers
  96   seeking public Certification and Endorsement from the Khronos Group
  97   under their Trademarked Certification Programme.
  98
  99 # Proposed Opcodes vs Khronos OpenCL vs IEEE754-2019<a name="khronos_equiv"></a>
 100
 101 This list shows the (direct) equivalence between proposed opcodes,
 102 their Khronos OpenCL equivalents, and their IEEE754-2019 equivalents.
 103 98% of the opcodes in this proposal that are in the IEEE754-2019 standard
 104 are present in the Khronos Extended Instruction Set.
 105
 106 See
 107 <https://www.khronos.org/registry/spir-v/specs/unified1/OpenCL.ExtendedInstructionSet.100.html>
 108 and <https://ieeexplore.ieee.org/document/8766229>
 109
 110 * Special FP16 opcodes are *not* being proposed, except by indirect / inherent
 111   use of elwidth overrides that is already present in the SVP64 Specification.
 112 * "Native" opcodes are *not* being proposed: implementors will be expected
 113   to use the (equivalent) proposed opcode covering the same function.
 114 * "Fast" opcodes are *not* being proposed, because the Khronos Specification
 115   fast\_length, fast\_normalise and fast\_distance OpenCL opcodes require
 116   vectors (or can be done as scalar operations using other Power ISA
 117   instructions).
 118
 119 The OpenCL FP32 opcodes are **direct** equivalents to the proposed opcodes.
 120 Deviation from conformance with the Khronos Specification - including the
 121 Khronos Specification accuracy requirements - is not an option, as it
 122 results in non-compliance, and the vendor may not use the Trademarked words
 123 "Vulkan" etc. in conjunction with their product.
 124
 125 IEEE754-2019 Table 9.1 lists "additional mathematical operations".
 126 Interestingly the only functions missing when compared to OpenCL are
 127 compound, exp2m1, exp10m1, log2p1, log10p1, pown (integer power) and powr.
 128
 129 |opcode  |OpenCL FP32|OpenCL FP16|OpenCL native|OpenCL fast|IEEE754 |Power ISA |
 130 |------- |-----------|-----------|-------------|-----------|------- |--------- |
 131 |FSIN    |sin        |half\_sin  |native\_sin  |NONE       |sin     |NONE      |
 132 |FCOS    |cos        |half\_cos  |native\_cos  |NONE       |cos     |NONE      |
 133 |FTAN    |tan        |half\_tan  |native\_tan  |NONE       |tan     |NONE      |
 134 |NONE (1)|sincos     |NONE       |NONE         |NONE       |NONE    |NONE      |
 135 |FASIN   |asin       |NONE       |NONE         |NONE       |asin    |NONE      |
 136 |FACOS   |acos       |NONE       |NONE         |NONE       |acos    |NONE      |
 137 |FATAN   |atan       |NONE       |NONE         |NONE       |atan    |NONE      |
 138 |FSINPI  |sinpi      |NONE       |NONE         |NONE       |sinPi   |NONE      |
 139 |FCOSPI  |cospi      |NONE       |NONE         |NONE       |cosPi   |NONE      |
 140 |FTANPI  |tanpi      |NONE       |NONE         |NONE       |tanPi   |NONE      |
 141 |FASINPI |asinpi     |NONE       |NONE         |NONE       |asinPi  |NONE      |
 142 |FACOSPI |acospi     |NONE       |NONE         |NONE       |acosPi  |NONE      |
 143 |FATANPI |atanpi     |NONE       |NONE         |NONE       |atanPi  |NONE      |
 144 |FSINH   |sinh       |NONE       |NONE         |NONE       |sinh    |NONE      |
 145 |FCOSH   |cosh       |NONE       |NONE         |NONE       |cosh    |NONE      |
 146 |FTANH   |tanh       |NONE       |NONE         |NONE       |tanh    |NONE      |
 147 |FASINH  |asinh      |NONE       |NONE         |NONE       |asinh   |NONE      |
 148 |FACOSH  |acosh      |NONE       |NONE         |NONE       |acosh   |NONE      |
 149 |FATANH  |atanh      |NONE       |NONE         |NONE       |atanh   |NONE      |
 150 |FATAN2  |atan2      |NONE       |NONE         |NONE       |atan2   |NONE      |
 151 |FATAN2PI|atan2pi    |NONE       |NONE         |NONE       |atan2pi |NONE      |
 152 |FRSQRT  |rsqrt      |half\_rsqrt|native\_rsqrt|NONE       |rSqrt   |fsqrte, fsqrtes (4)   |
 153 |FCBRT   |cbrt       |NONE       |NONE         |NONE       |NONE (2)|NONE      |
 154 |FEXP2   |exp2       |half\_exp2 |native\_exp2 |NONE       |exp2    |NONE      |
 155 |FLOG2   |log2       |half\_log2 |native\_log2 |NONE       |log2    |NONE      |
 156 |FEXPM1  |expm1      |NONE       |NONE         |NONE       |expm1   |NONE      |
 157 |FLOG1P  |log1p      |NONE       |NONE         |NONE       |logp1   |NONE      |
 158 |FEXP    |exp        |half\_exp  |native\_exp  |NONE       |exp     |NONE      |
 159 |FLOG    |log        |half\_log  |native\_log  |NONE       |log     |NONE      |
 160 |FEXP10  |exp10      |half\_exp10|native\_exp10|NONE       |exp10   |NONE      |
 161 |FLOG10  |log10      |half\_log10|native\_log10|NONE       |log10   |NONE      |
 162 |FPOW    |pow        |NONE       |NONE         |NONE       |pow     |NONE      |
 163 |FPOWN   |pown       |NONE       |NONE         |NONE       |pown    |NONE      |
 164 |FPOWR   |powr       |half\_powr |native\_powr |NONE       |powr    |NONE      |
 165 |FROOTN  |rootn      |NONE       |NONE         |NONE       |rootn   |NONE      |
 166 |FHYPOT  |hypot      |NONE       |NONE         |NONE       |hypot   |NONE      |
 167 |FRECIP  |NONE       |half\_recip|native\_recip|NONE       |NONE (3)|fre, fres (4)        |
 168 |NONE    |NONE       |NONE       |NONE         |NONE       |compound|NONE      |
 169 |FEXP2M1   |NONE       |NONE       |NONE         |NONE       |exp2m1  |NONE      |
 170 |FEXP10M1   |NONE       |NONE       |NONE         |NONE       |exp10m1 |NONE      |
 171 |FLOG2P1  |NONE       |NONE       |NONE         |NONE       |log2p1  |NONE      |
 172 |FLOG10P1  |NONE       |NONE       |NONE         |NONE       |log10p1 |NONE      |
 173
 174 Note (1) FSINCOS is macro-op fused (see below).
 175
 176 Note (2) synthesised in IEEE754-2019 as "pown(x, 3)"
 177
 178 Note (3) synthesised in IEEE754-2019 using "1.0 / x"
 179
 180 Note (4) these are estimate opcodes that help accelerate
 181 software emulation
 182
 183 ## List of 2-arg opcodes
 184
 185 | opcode   | Description            | pseudocode                 | Extension   |
 186 | ------   | ----------------       | ----------------           | ----------- |
 187 | FATAN2   | atan2 arc tangent      | FRT = atan2(FRB, FRA)       | Zarctrignpi |
 188 | FATAN2PI | atan2 arc tangent / pi | FRT = atan2(FRB, FRA) / pi  | Zarctrigpi  |
 189 | FPOW     | x power of y           | FRT = pow(FRA, FRB)         | ZftransAdv  |
 190 | FPOWN    | x power of n (n int)   | FRT = pow(FRA, FRB)         | ZftransAdv  |
 191 | FPOWR    | x power of y (x +ve)   | FRT = exp(FRA log(FRB))     | ZftransAdv  |
 192 | FROOTN   | x power 1/n (n integer)| FRT = pow(FRA, 1/FRB)       | ZftransAdv  |
 193 | FHYPOT   | hypotenuse             | FRT = sqrt(FRA^2 + FRB^2)   | ZftransAdv  |
 194
 195 ## List of 1-arg transcendental opcodes
 196
 197 | opcode   | Description            | pseudocode                 | Extension   |
 198 | ------   | ----------------       | ----------------           | ----------- |
 199 | FRSQRT   | Reciprocal Square-root   | FRT = sqrt(FRA)          | Zfrsqrt    |
 200 | FCBRT    | Cube Root                | FRT = pow(FRA, 1.0 / 3)  | ZftransAdv |
 201 | FRECIP   | Reciprocal               | FRT = 1.0 / FRA          | Zftrans    |
 202 | FEXP2M1   | power-2 minus 1      | FRT = pow(2, FRA) - 1.0  | ZftransExt |
 203 | FLOG2P1   | log2 plus 1               | FRT = log(2, 1 + FRA)    | ZftransExt |
 204 | FEXP2    | power-of-2               | FRT = pow(2, FRA)        | Zftrans    |
 205 | FLOG2    | log2                     | FRT = log(2. FRA)        | Zftrans    |
 206 | FEXPM1   | exponential minus 1      | FRT = pow(e, FRA) - 1.0  | ZftransExt |
 207 | FLOG1P   | log plus 1               | FRT = log(e, 1 + FRA)    | ZftransExt |
 208 | FEXP     | exponential              | FRT = pow(e, FRA)        | ZftransExt |
 209 | FLOG     | natural log (base e)     | FRT = log(e, FRA)        | ZftransExt |
 210 | FEXP10M1   | power-10 minus 1      | FRT = pow(10, FRA) - 1.0  | ZftransExt |
 211 | FLOG10P1   | log10 plus 1               | FRT = log(10, 1 + FRA)    | ZftransExt |
 212 | FEXP10   | power-of-10              | FRT = pow(10, FRA)       | ZftransExt |
 213 | FLOG10   | log base 10              | FRT = log(10, FRA)       | ZftransExt |
 214
 215 ## List of 1-arg trigonometric opcodes
 216
 217 | opcode   | Description            | pseudocode                 | Extension   |
 218 | ------   | ----------------       | ----------------           | ----------- |
 219 | FSIN     | sin (radians)            | FRT = sin(FRA)           | Ztrignpi    |
 220 | FCOS     | cos (radians)            | FRT = cos(FRA)           | Ztrignpi    |
 221 | FTAN     | tan (radians)            | FRT = tan(FRA)           | Ztrignpi    |
 222 | FASIN    | arcsin (radians)         | FRT = asin(FRA)          | Zarctrignpi |
 223 | FACOS    | arccos (radians)         | FRT = acos(FRA)          | Zarctrignpi |
 224 | FATAN    | arctan (radians)         | FRT = atan(FRA)          | Zarctrignpi |
 225 | FSINPI   | sin times pi             | FRT = sin(pi * FRA)      | Ztrigpi |
 226 | FCOSPI   | cos times pi             | FRT = cos(pi * FRA)      | Ztrigpi |
 227 | FTANPI   | tan times pi             | FRT = tan(pi * FRA)      | Ztrigpi |
 228 | FASINPI  | arcsin / pi              | FRT = asin(FRA) / pi     | Zarctrigpi |
 229 | FACOSPI  | arccos / pi              | FRT = acos(FRA) / pi     | Zarctrigpi |
 230 | FATANPI  | arctan / pi              | FRT = atan(FRA) / pi     | Zarctrigpi |
 231 | FSINH    | hyperbolic sin (radians) | FRT = sinh(FRA)          | Zfhyp |
 232 | FCOSH    | hyperbolic cos (radians) | FRT = cosh(FRA)          | Zfhyp |
 233 | FTANH    | hyperbolic tan (radians) | FRT = tanh(FRA)          | Zfhyp |
 234 | FASINH   | inverse hyperbolic sin   | FRT = asinh(FRA)         | Zfhyp |
 235 | FACOSH   | inverse hyperbolic cos   | FRT = acosh(FRA)         | Zfhyp |
 236 | FATANH   | inverse hyperbolic tan   | FRT = atanh(FRA)         | Zfhyp |
 237
 238 [[!inline pages="openpower/power_trans_ops" raw=yes ]]
 239
 240 # Subsets
 241
 242 The full set is based on the Khronos OpenCL opcodes. If implemented
 243 entirely it would be too much for both Embedded and also 3D.
 244
 245 The subsets are organised by hardware complexity, need (3D, HPC), however
 246 due to synthesis producing inaccurate results at the range limits,
 247 the less common subsets are still required for IEEE754 HPC.
 248
 249 MALI Midgard, an embedded / mobile 3D GPU, for example only has the
 250 following opcodes:
 251
 252     E8 - fatan_pt2
 253     F0 - frcp (reciprocal)
 254     F2 - frsqrt (inverse square root, 1/sqrt(x))
 255     F3 - fsqrt (square root)
 256     F4 - fexp2 (2^x)
 257     F5 - flog2
 258     F6 - fsin1pi
 259     F7 - fcos1pi
 260     F9 - fatan_pt1
 261
 262 These in FP32 and FP16 only: no FP64 hardware, at all.
 263
 264 Vivante Embedded/Mobile 3D (etnaviv
 265 <https://github.com/laanwj/etna_viv/blob/master/rnndb/isa.xml>)
 266 only has the following:
 267
 268     sin, cos2pi
 269     cos, sin2pi
 270     log2, exp
 271     sqrt and rsqrt
 272     recip.
 273
 274 It also has fast variants of some of these, as a CSR Mode.
 275
 276 AMD's R600 GPU (R600\_Instruction\_Set\_Architecture.pdf) and the
 277 RDNA ISA (RDNA\_Shader\_ISA\_5August2019.pdf, Table 22, Section 6.3) have:
 278
 279     COS2PI (appx)
 280     EXP2
 281     LOG (IEEE754)
 282     RECIP
 283     RSQRT
 284     SQRT
 285     SIN2PI (appx)
 286
 287 AMD RDNA has F16 and F32 variants of all the above, and also has F64
 288 variants of SQRT, RSQRT and RECIP.  It is interesting that even the
 289 modern high-end AMD GPU does not have TAN or ATAN, where MALI Midgard
 290 does.
 291
 292 Also a general point, that customised optimised hardware targetting
 293 FP32 3D with less accuracy simply can neither be used for IEEE754 nor
 294 for FP64 (except as a starting point for hardware or software driven
 295 Newton Raphson or other iterative method).
 296
 297 Also in cost/area sensitive applications even the extra ROM lookup tables
 298 for certain algorithms may be too costly.
 299
 300 These wildly differing and incompatible driving factors lead to the
 301 subset subdivisions, below.
 302
 303 ## Transcendental Subsets
 304
 305 ### Zftrans
 306
 307 LOG2 EXP2 RECIP RSQRT
 308
 309 Zftrans contains the minimum standard transcendentals best suited to
 310 3D. They are also the minimum subset for synthesising log10, exp10,
 311 exp1m, log1p, the hyperbolic trigonometric functions sinh and so on.
 312
 313 They are therefore considered "base" (essential) transcendentals.
 314
 315 ### ZftransExt
 316
 317 LOG, EXP, EXP10, LOG10, LOGP1, EXP1M
 318
 319 These are extra transcendental functions that are useful, not generally
 320 needed for 3D, however for Numerical Computation they may be useful.
 321
 322 Although they can be synthesised using Ztrans (LOG2 multiplied
 323 by a constant), there is both a performance penalty as well as an
 324 accuracy penalty towards the limits, which for IEEE754 compliance is
 325 unacceptable. In particular, LOG(1+FRA) in hardware may give much better
 326 accuracy at the lower end (very small FRA) than LOG(FRA).
 327
 328 Their forced inclusion would be inappropriate as it would penalise
 329 embedded systems with tight power and area budgets.  However if they
 330 were completely excluded the HPC applications would be penalised on
 331 performance and accuracy.
 332
 333 Therefore they are their own subset extension.
 334
 335 ### Zfhyp
 336
 337 SINH, COSH, TANH, ASINH, ACOSH, ATANH
 338
 339 These are the hyperbolic/inverse-hyperbolic functions. Their use in 3D
 340 is limited.
 341
 342 They can all be synthesised using LOG, SQRT and so on, so depend
 343 on Zftrans.  However, once again, at the limits of the range, IEEE754
 344 compliance becomes impossible, and thus a hardware implementation may
 345 be required.
 346
 347 HPC and high-end GPUs are likely markets for these.
 348
 349 ### ZftransAdv
 350
 351 CBRT, POW, POWN, POWR, ROOTN
 352
 353 These are simply much more complex to implement in hardware, and typically
 354 will only be put into HPC applications.
 355
 356 * **Zfrsqrt**: Reciprocal square-root.
 357
 358 ## Trigonometric subsets
 359
 360 ### Ztrigpi vs Ztrignpi
 361
 362 * **Ztrigpi**: SINPI COSPI TANPI
 363 * **Ztrignpi**: SIN COS TAN
 364
 365 Ztrignpi are the basic trigonometric functions through which all others
 366 could be synthesised, and they are typically the base trigonometrics
 367 provided by GPUs for 3D, warranting their own subset.
 368
 369 In the case of the Ztrigpi subset, these are commonly used in for loops
 370 with a power of two number of subdivisions, and the cost of multiplying
 371 by PI inside each loop (or cumulative addition, resulting in cumulative
 372 errors) is not acceptable.
 373
 374 In for example CORDIC the multiplication by PI may be moved outside of
 375 the hardware algorithm as a loop invariant, with no power or area penalty.
 376
 377 Again, therefore, if SINPI (etc.) were excluded, programmers would be
 378 penalised by being forced to divide by PI in some circumstances. Likewise
 379 if SIN were excluded, programmers would be penaslised by being forced
 380 to *multiply* by PI in some circumstances.
 381
 382 Thus again, a slightly different application of the same general argument
 383 applies to give Ztrignpi and Ztrigpi as subsets.  3D GPUs will almost
 384 certainly provide both.
 385
 386 ### Zarctrigpi and Zarctrignpi
 387
 388 * **Zarctrigpi**: ATAN2PI ASINPI ACOSPI
 389 * **Zarctrignpi**: ATAN2 ACOS ASIN
 390
 391 These are extra trigonometric functions that are useful in some
 392 applications, but even for 3D GPUs, particularly embedded and mobile class
 393 GPUs, they are not so common and so are typically synthesised, there.
 394
 395 Although they can be synthesised using Ztrigpi and Ztrignpi, there is,
 396 once again, both a performance penalty as well as an accuracy penalty
 397 towards the limits, which for IEEE754 compliance is unacceptable, yet
 398 is acceptable for 3D.
 399
 400 Therefore they are their own subset extensions.
 401
 402 # Synthesis, Pseudo-code ops and macro-ops
 403
 404 The pseudo-ops are best left up to the compiler rather than being actual
 405 pseudo-ops, by allocating one scalar FP register for use as a constant
 406 (loop invariant) set to "1.0" at the beginning of a function or other
 407 suitable code block.
 408
 409 * FSINCOS - fused macro-op between FSIN and FCOS (issued in that order).
 410 * FSINCOSPI - fused macro-op between FSINPI and FCOSPI (issued in that order).
 411
 412 FATANPI example pseudo-code:
 413
 414     fmvis ft0, 0x3F800 // upper bits of f32 1.0 (BF16)
 415     fatan2pis FRT, FRA, ft0
 416
 417 Hyperbolic function example (obviates need for Zfhyp except for
 418 high-performance or correctly-rounding):
 419
 420     ASINH( x ) = ln( x + SQRT(x**2+1))
 421
 422 # Evaluation and commentary
 423
 424 Moved to [[discussion]]
 425