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 |NONE    |NONE       |NONE       |NONE         |NONE       |exp2m1  |NONE      |
 170 |NONE    |NONE       |NONE       |NONE         |NONE       |exp10m1 |NONE      |
 171 |NONE    |NONE       |NONE       |NONE         |NONE       |log2p1  |NONE      |
 172 |NONE    |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      | rd = atan2(rs2, rs1)       | Zarctrignpi |
 188 | FATAN2PI | atan2 arc tangent / pi | rd = atan2(rs2, rs1) / pi  | Zarctrigpi  |
 189 | FPOW     | x power of y           | rd = pow(rs1, rs2)         | ZftransAdv  |
 190 | FPOWN    | x power of n (n int)   | rd = pow(rs1, rs2)         | ZftransAdv  |
 191 | FPOWR    | x power of y (x +ve)   | rd = exp(rs1 log(rs2))     | ZftransAdv  |
 192 | FROOTN   | x power 1/n (n integer)| rd = pow(rs1, 1/rs2)       | ZftransAdv  |
 193 | FHYPOT   | hypotenuse             | rd = sqrt(rs1^2 + rs2^2)   | ZftransAdv  |
 194
 195 ## List of 1-arg transcendental opcodes
 196
 197 | opcode   | Description            | pseudocode                 | Extension   |
 198 | ------   | ----------------       | ----------------           | ----------- |
 199 | FRSQRT   | Reciprocal Square-root   | rd = sqrt(rs1)          | Zfrsqrt    |
 200 | FCBRT    | Cube Root                | rd = pow(rs1, 1.0 / 3)  | ZftransAdv |
 201 | FRECIP   | Reciprocal               | rd = 1.0 / rs1          | Zftrans    |
 202 | FEXP2    | power-of-2               | rd = pow(2, rs1)        | Zftrans    |
 203 | FLOG2    | log2                     | rd = log(2. rs1)        | Zftrans    |
 204 | FEXPM1   | exponential minus 1      | rd = pow(e, rs1) - 1.0  | ZftransExt |
 205 | FLOG1P   | log plus 1               | rd = log(e, 1 + rs1)    | ZftransExt |
 206 | FEXP     | exponential              | rd = pow(e, rs1)        | ZftransExt |
 207 | FLOG     | natural log (base e)     | rd = log(e, rs1)        | ZftransExt |
 208 | FEXP10   | power-of-10              | rd = pow(10, rs1)       | ZftransExt |
 209 | FLOG10   | log base 10              | rd = log(10, rs1)       | ZftransExt |
 210
 211 ## List of 1-arg trigonometric opcodes
 212
 213 | opcode   | Description            | pseudocode                 | Extension   |
 214 | ------   | ----------------       | ----------------           | ----------- |
 215 | FSIN     | sin (radians)            | rd = sin(rs1)           | Ztrignpi    |
 216 | FCOS     | cos (radians)            | rd = cos(rs1)           | Ztrignpi    |
 217 | FTAN     | tan (radians)            | rd = tan(rs1)           | Ztrignpi    |
 218 | FASIN    | arcsin (radians)         | rd = asin(rs1)          | Zarctrignpi |
 219 | FACOS    | arccos (radians)         | rd = acos(rs1)          | Zarctrignpi |
 220 | FATAN    | arctan (radians)         | rd = atan(rs1)          | Zarctrignpi |
 221 | FSINPI   | sin times pi             | rd = sin(pi * rs1)      | Ztrigpi |
 222 | FCOSPI   | cos times pi             | rd = cos(pi * rs1)      | Ztrigpi |
 223 | FTANPI   | tan times pi             | rd = tan(pi * rs1)      | Ztrigpi |
 224 | FASINPI  | arcsin / pi              | rd = asin(rs1) / pi     | Zarctrigpi |
 225 | FACOSPI  | arccos / pi              | rd = acos(rs1) / pi     | Zarctrigpi |
 226 | FATANPI  | arctan / pi              | rd = atan(rs1) / pi     | Zarctrigpi |
 227 | FSINH    | hyperbolic sin (radians) | rd = sinh(rs1)          | Zfhyp |
 228 | FCOSH    | hyperbolic cos (radians) | rd = cosh(rs1)          | Zfhyp |
 229 | FTANH    | hyperbolic tan (radians) | rd = tanh(rs1)          | Zfhyp |
 230 | FASINH   | inverse hyperbolic sin   | rd = asinh(rs1)         | Zfhyp |
 231 | FACOSH   | inverse hyperbolic cos   | rd = acosh(rs1)         | Zfhyp |
 232 | FATANH   | inverse hyperbolic tan   | rd = atanh(rs1)         | Zfhyp |
 233
 234 [[!inline pages="openpower/power_trans_ops" raw=yes ]]
 235
 236 # Subsets
 237
 238 The full set is based on the Khronos OpenCL opcodes. If implemented
 239 entirely it would be too much for both Embedded and also 3D.
 240
 241 The subsets are organised by hardware complexity, need (3D, HPC), however
 242 due to synthesis producing inaccurate results at the range limits,
 243 the less common subsets are still required for IEEE754 HPC.
 244
 245 MALI Midgard, an embedded / mobile 3D GPU, for example only has the
 246 following opcodes:
 247
 248     E8 - fatan_pt2
 249     F0 - frcp (reciprocal)
 250     F2 - frsqrt (inverse square root, 1/sqrt(x))
 251     F3 - fsqrt (square root)
 252     F4 - fexp2 (2^x)
 253     F5 - flog2
 254     F6 - fsin1pi
 255     F7 - fcos1pi
 256     F9 - fatan_pt1
 257
 258 These in FP32 and FP16 only: no FP64 hardware, at all.
 259
 260 Vivante Embedded/Mobile 3D (etnaviv
 261 <https://github.com/laanwj/etna_viv/blob/master/rnndb/isa.xml>)
 262 only has the following:
 263
 264     sin, cos2pi
 265     cos, sin2pi
 266     log2, exp
 267     sqrt and rsqrt
 268     recip.
 269
 270 It also has fast variants of some of these, as a CSR Mode.
 271
 272 AMD's R600 GPU (R600\_Instruction\_Set\_Architecture.pdf) and the
 273 RDNA ISA (RDNA\_Shader\_ISA\_5August2019.pdf, Table 22, Section 6.3) have:
 274
 275     COS2PI (appx)
 276     EXP2
 277     LOG (IEEE754)
 278     RECIP
 279     RSQRT
 280     SQRT
 281     SIN2PI (appx)
 282
 283 AMD RDNA has F16 and F32 variants of all the above, and also has F64
 284 variants of SQRT, RSQRT and RECIP.  It is interesting that even the
 285 modern high-end AMD GPU does not have TAN or ATAN, where MALI Midgard
 286 does.
 287
 288 Also a general point, that customised optimised hardware targetting
 289 FP32 3D with less accuracy simply can neither be used for IEEE754 nor
 290 for FP64 (except as a starting point for hardware or software driven
 291 Newton Raphson or other iterative method).
 292
 293 Also in cost/area sensitive applications even the extra ROM lookup tables
 294 for certain algorithms may be too costly.
 295
 296 These wildly differing and incompatible driving factors lead to the
 297 subset subdivisions, below.
 298
 299 ## Transcendental Subsets
 300
 301 ### Zftrans
 302
 303 LOG2 EXP2 RECIP RSQRT
 304
 305 Zftrans contains the minimum standard transcendentals best suited to
 306 3D. They are also the minimum subset for synthesising log10, exp10,
 307 exp1m, log1p, the hyperbolic trigonometric functions sinh and so on.
 308
 309 They are therefore considered "base" (essential) transcendentals.
 310
 311 ### ZftransExt
 312
 313 LOG, EXP, EXP10, LOG10, LOGP1, EXP1M
 314
 315 These are extra transcendental functions that are useful, not generally
 316 needed for 3D, however for Numerical Computation they may be useful.
 317
 318 Although they can be synthesised using Ztrans (LOG2 multiplied
 319 by a constant), there is both a performance penalty as well as an
 320 accuracy penalty towards the limits, which for IEEE754 compliance is
 321 unacceptable. In particular, LOG(1+rs1) in hardware may give much better
 322 accuracy at the lower end (very small rs1) than LOG(rs1).
 323
 324 Their forced inclusion would be inappropriate as it would penalise
 325 embedded systems with tight power and area budgets.  However if they
 326 were completely excluded the HPC applications would be penalised on
 327 performance and accuracy.
 328
 329 Therefore they are their own subset extension.
 330
 331 ### Zfhyp
 332
 333 SINH, COSH, TANH, ASINH, ACOSH, ATANH
 334
 335 These are the hyperbolic/inverse-hyperbolic functions. Their use in 3D
 336 is limited.
 337
 338 They can all be synthesised using LOG, SQRT and so on, so depend
 339 on Zftrans.  However, once again, at the limits of the range, IEEE754
 340 compliance becomes impossible, and thus a hardware implementation may
 341 be required.
 342
 343 HPC and high-end GPUs are likely markets for these.
 344
 345 ### ZftransAdv
 346
 347 CBRT, POW, POWN, POWR, ROOTN
 348
 349 These are simply much more complex to implement in hardware, and typically
 350 will only be put into HPC applications.
 351
 352 * **Zfrsqrt**: Reciprocal square-root.
 353
 354 ## Trigonometric subsets
 355
 356 ### Ztrigpi vs Ztrignpi
 357
 358 * **Ztrigpi**: SINPI COSPI TANPI
 359 * **Ztrignpi**: SIN COS TAN
 360
 361 Ztrignpi are the basic trigonometric functions through which all others
 362 could be synthesised, and they are typically the base trigonometrics
 363 provided by GPUs for 3D, warranting their own subset.
 364
 365 In the case of the Ztrigpi subset, these are commonly used in for loops
 366 with a power of two number of subdivisions, and the cost of multiplying
 367 by PI inside each loop (or cumulative addition, resulting in cumulative
 368 errors) is not acceptable.
 369
 370 In for example CORDIC the multiplication by PI may be moved outside of
 371 the hardware algorithm as a loop invariant, with no power or area penalty.
 372
 373 Again, therefore, if SINPI (etc.) were excluded, programmers would be
 374 penalised by being forced to divide by PI in some circumstances. Likewise
 375 if SIN were excluded, programmers would be penaslised by being forced
 376 to *multiply* by PI in some circumstances.
 377
 378 Thus again, a slightly different application of the same general argument
 379 applies to give Ztrignpi and Ztrigpi as subsets.  3D GPUs will almost
 380 certainly provide both.
 381
 382 ### Zarctrigpi and Zarctrignpi
 383
 384 * **Zarctrigpi**: ATAN2PI ASINPI ACOSPI
 385 * **Zarctrignpi**: ATAN2 ACOS ASIN
 386
 387 These are extra trigonometric functions that are useful in some
 388 applications, but even for 3D GPUs, particularly embedded and mobile class
 389 GPUs, they are not so common and so are typically synthesised, there.
 390
 391 Although they can be synthesised using Ztrigpi and Ztrignpi, there is,
 392 once again, both a performance penalty as well as an accuracy penalty
 393 towards the limits, which for IEEE754 compliance is unacceptable, yet
 394 is acceptable for 3D.
 395
 396 Therefore they are their own subset extensions.
 397
 398 # Synthesis, Pseudo-code ops and macro-ops
 399
 400 The pseudo-ops are best left up to the compiler rather than being actual
 401 pseudo-ops, by allocating one scalar FP register for use as a constant
 402 (loop invariant) set to "1.0" at the beginning of a function or other
 403 suitable code block.
 404
 405 * FSINCOS - fused macro-op between FSIN and FCOS (issued in that order).
 406 * FSINCOSPI - fused macro-op between FSINPI and FCOSPI (issued in that order).
 407
 408 FATANPI example pseudo-code:
 409
 410     fmvis ft0, 0x3F800 // upper bits of f32 1.0 (BF16)
 411     fatan2pis FRT, FRA, ft0
 412
 413 Hyperbolic function example (obviates need for Zfhyp except for
 414 high-performance or correctly-rounding):
 415
 416     ASINH( x ) = ln( x + SQRT(x**2+1))
 417
 418 # Evaluation and commentary
 419
 420 Moved to [[discussion]]
 421