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