$p1 suldgb b32 $r3 cv zero u8 g[$r4d] $r2 $p0
long mov b32 $r3 0x3f800000
long nop
+sched 0x00 0x00 0x00 0x00 0x00 0x00 0x00
+long nop
long ret
// SIZE: 9 * 8 bytes
//
gk104_rcp_f64:
- long nop
+ // Step 1: classify input according to exponent and value, and calculate
+ // result for 0/inf/nan. $r2 holds the exponent value, which starts at
+ // bit 52 (bit 20 of the upper half) and is 11 bits in length
+ ext u32 $r2 $r1 0xb14
+ add b32 $r3 $r2 0xffffffff
+ joinat #rcp_rejoin
+ // We want to check whether the exponent is 0 or 0x7ff (i.e. NaN, inf,
+ // denorm, or 0). Do this by substracting 1 from the exponent, which will
+ // mean that it's > 0x7fd in those cases when doing unsigned comparison
+ set $p0 0x1 gt u32 $r3 0x7fd
+ // $r3: 0 for norms, 0x36 for denorms, -1 for others
+ long mov b32 $r3 0x0
+ sched 0x2f 0x04 0x2d 0x2b 0x2f 0x28 0x28
+ join (not $p0) nop
+ // Process all special values: NaN, inf, denorm, 0
+ mov b32 $r3 0xffffffff
+ // A number is NaN if its abs value is greater than or unordered with inf
+ set $p0 0x1 gtu f64 abs $r0d 0x7ff0000000000000
+ (not $p0) bra #rcp_inf_or_denorm_or_zero
+ // NaN -> NaN, the next line sets the "quiet" bit of the result. This
+ // behavior is both seen on the CPU and the blob
+ join or b32 $r1 $r1 0x80000
+rcp_inf_or_denorm_or_zero:
+ and b32 $r4 $r1 0x7ff00000
+ // Other values with nonzero in exponent field should be inf
+ set $p0 0x1 eq s32 $r4 0x0
+ sched 0x2b 0x04 0x2f 0x2d 0x2b 0x2f 0x20
+ $p0 bra #rcp_denorm_or_zero
+ // +/-Inf -> +/-0
+ xor b32 $r1 $r1 0x7ff00000
+ join mov b32 $r0 0x0
+rcp_denorm_or_zero:
+ set $p0 0x1 gtu f64 abs $r0d 0x0
+ $p0 bra #rcp_denorm
+ // +/-0 -> +/-Inf
+ join or b32 $r1 $r1 0x7ff00000
+rcp_denorm:
+ // non-0 denorms: multiply with 2^54 (the 0x36 in $r3), join with norms
+ mul rn f64 $r0d $r0d 0x4350000000000000
+ sched 0x2f 0x28 0x2b 0x28 0x28 0x04 0x28
+ join mov b32 $r3 0x36
+rcp_rejoin:
+ // All numbers with -1 in $r3 have their result ready in $r0d, return them
+ // others need further calculation
+ set $p0 0x1 lt s32 $r3 0x0
+ $p0 bra #rcp_end
+ // Step 2: Before the real calculation goes on, renormalize the values to
+ // range [1, 2) by setting exponent field to 0x3ff (the exponent of 1)
+ // result in $r6d. The exponent will be recovered later.
+ ext u32 $r2 $r1 0xb14
+ and b32 $r7 $r1 0x800fffff
+ add b32 $r7 $r7 0x3ff00000
+ long mov b32 $r6 $r0
+ sched 0x2b 0x04 0x28 0x28 0x2a 0x2b 0x2e
+ // Step 3: Convert new value to float (no overflow will occur due to step
+ // 2), calculate rcp and do newton-raphson step once
+ cvt rz f32 $r5 f64 $r6d
+ long rcp f32 $r4 $r5
+ mov b32 $r0 0xbf800000
+ fma rn f32 $r5 $r4 $r5 $r0
+ fma rn f32 $r0 neg $r4 $r5 $r4
+ // Step 4: convert result $r0 back to double, do newton-raphson steps
+ cvt f64 $r0d f32 $r0
+ cvt f64 $r6d neg f64 $r6d
+ sched 0x2e 0x29 0x29 0x29 0x29 0x29 0x29
+ cvt f64 $r8d f32 0x3f800000
+ // 4 Newton-Raphson Steps, tmp in $r4d, result in $r0d
+ // The formula used here (and above) is:
+ // RCP_{n + 1} = 2 * RCP_{n} - x * RCP_{n} * RCP_{n}
+ // The following code uses 2 FMAs for each step, and it will basically
+ // looks like:
+ // tmp = -src * RCP_{n} + 1
+ // RCP_{n + 1} = RCP_{n} * tmp + RCP_{n}
+ fma rn f64 $r4d $r6d $r0d $r8d
+ fma rn f64 $r0d $r0d $r4d $r0d
+ fma rn f64 $r4d $r6d $r0d $r8d
+ fma rn f64 $r0d $r0d $r4d $r0d
+ fma rn f64 $r4d $r6d $r0d $r8d
+ fma rn f64 $r0d $r0d $r4d $r0d
+ sched 0x29 0x20 0x28 0x28 0x28 0x28 0x28
+ fma rn f64 $r4d $r6d $r0d $r8d
+ fma rn f64 $r0d $r0d $r4d $r0d
+ // Step 5: Exponent recovery and final processing
+ // The exponent is recovered by adding what we added to the exponent.
+ // Suppose we want to calculate rcp(x), but we have rcp(cx), then
+ // rcp(x) = c * rcp(cx)
+ // The delta in exponent comes from two sources:
+ // 1) The renormalization in step 2. The delta is:
+ // 0x3ff - $r2
+ // 2) (For the denorm input) The 2^54 we multiplied at rcp_denorm, stored
+ // in $r3
+ // These 2 sources are calculated in the first two lines below, and then
+ // added to the exponent extracted from the result above.
+ // Note that after processing, the new exponent may >= 0x7ff (inf)
+ // or <= 0 (denorm). Those cases will be handled respectively below
+ subr b32 $r2 $r2 0x3ff
+ long add b32 $r4 $r2 $r3
+ ext u32 $r3 $r1 0xb14
+ // New exponent in $r3
+ long add b32 $r3 $r3 $r4
+ add b32 $r2 $r3 0xffffffff
+ sched 0x28 0x2b 0x28 0x2b 0x28 0x28 0x2b
+ // (exponent-1) < 0x7fe (unsigned) means the result is in norm range
+ // (same logic as in step 1)
+ set $p0 0x1 lt u32 $r2 0x7fe
+ (not $p0) bra #rcp_result_inf_or_denorm
+ // Norms: convert exponents back and return
+ shl b32 $r4 $r4 clamp 0x14
+ long add b32 $r1 $r4 $r1
+ bra #rcp_end
+rcp_result_inf_or_denorm:
+ // New exponent >= 0x7ff means that result is inf
+ set $p0 0x1 ge s32 $r3 0x7ff
+ (not $p0) bra #rcp_result_denorm
+ sched 0x20 0x25 0x28 0x2b 0x23 0x25 0x2f
+ // Infinity
+ and b32 $r1 $r1 0x80000000
+ long mov b32 $r0 0x0
+ add b32 $r1 $r1 0x7ff00000
+ bra #rcp_end
+rcp_result_denorm:
+ // Denorm result comes from huge input. The greatest possible fp64, i.e.
+ // 0x7fefffffffffffff's rcp is 0x0004000000000000, 1/4 of the smallest
+ // normal value. Other rcp result should be greater than that. If we
+ // set the exponent field to 1, we can recover the result by multiplying
+ // it with 1/2 or 1/4. 1/2 is used if the "exponent" $r3 is 0, otherwise
+ // 1/4 ($r3 should be -1 then). This is quite tricky but greatly simplifies
+ // the logic here.
+ set $p0 0x1 ne u32 $r3 0x0
+ and b32 $r1 $r1 0x800fffff
+ // 0x3e800000: 1/4
+ $p0 cvt f64 $r6d f32 0x3e800000
+ sched 0x2f 0x28 0x2c 0x2e 0x2a 0x20 0x27
+ // 0x3f000000: 1/2
+ (not $p0) cvt f64 $r6d f32 0x3f000000
+ add b32 $r1 $r1 0x00100000
+ mul rn f64 $r0d $r0d $r6d
+rcp_end:
long ret
// RSQ F64: Newton Raphson rsqrt(x): r_{i+1} = r_i * (1.5 - 0.5 * x * r_i * r_i)
// SIZE: 14 * 8 bytes
//
gk104_rsq_f64:
- long nop
+ // Before getting initial result rsqrt64h, two special cases should be
+ // handled first.
+ // 1. NaN: set the highest bit in mantissa so it'll be surely recognized
+ // as NaN in rsqrt64h
+ set $p0 0x1 gtu f64 abs $r0d 0x7ff0000000000000
+ $p0 or b32 $r1 $r1 0x00080000
+ and b32 $r2 $r1 0x7fffffff
+ sched 0x27 0x20 0x28 0x2c 0x25 0x28 0x28
+ // 2. denorms and small normal values: using their original value will
+ // lose precision either at rsqrt64h or the first step in newton-raphson
+ // steps below. Take 2 as a threshold in exponent field, and multiply
+ // with 2^54 if the exponent is smaller or equal. (will multiply 2^27
+ // to recover in the end)
+ ext u32 $r3 $r1 0xb14
+ set $p1 0x1 le u32 $r3 0x2
+ long or b32 $r2 $r0 $r2
+ $p1 mul rn f64 $r0d $r0d 0x4350000000000000
+ rsqrt64h $r5 $r1
+ // rsqrt64h will give correct result for 0/inf/nan, the following logic
+ // checks whether the input is one of those (exponent is 0x7ff or all 0
+ // except for the sign bit)
+ set b32 $r6 ne u32 $r3 0x7ff
+ long and b32 $r2 $r2 $r6
+ sched 0x28 0x2b 0x20 0x27 0x28 0x2e 0x28
+ set $p0 0x1 ne u32 $r2 0x0
+ $p0 bra #rsq_norm
+ // For 0/inf/nan, make sure the sign bit agrees with input and return
+ and b32 $r1 $r1 0x80000000
+ long mov b32 $r0 0x0
+ long or b32 $r1 $r1 $r5
+ long ret
+rsq_norm:
+ // For others, do 4 Newton-Raphson steps with the formula:
+ // RSQ_{n + 1} = RSQ_{n} * (1.5 - 0.5 * x * RSQ_{n} * RSQ_{n})
+ // In the code below, each step is written as:
+ // tmp1 = 0.5 * x * RSQ_{n}
+ // tmp2 = -RSQ_{n} * tmp1 + 0.5
+ // RSQ_{n + 1} = RSQ_{n} * tmp2 + RSQ_{n}
+ long mov b32 $r4 0x0
+ sched 0x2f 0x29 0x29 0x29 0x29 0x29 0x29
+ // 0x3f000000: 1/2
+ cvt f64 $r8d f32 0x3f000000
+ mul rn f64 $r2d $r0d $r8d
+ mul rn f64 $r0d $r2d $r4d
+ fma rn f64 $r6d neg $r4d $r0d $r8d
+ fma rn f64 $r4d $r4d $r6d $r4d
+ mul rn f64 $r0d $r2d $r4d
+ fma rn f64 $r6d neg $r4d $r0d $r8d
+ sched 0x29 0x29 0x29 0x29 0x29 0x29 0x29
+ fma rn f64 $r4d $r4d $r6d $r4d
+ mul rn f64 $r0d $r2d $r4d
+ fma rn f64 $r6d neg $r4d $r0d $r8d
+ fma rn f64 $r4d $r4d $r6d $r4d
+ mul rn f64 $r0d $r2d $r4d
+ fma rn f64 $r6d neg $r4d $r0d $r8d
+ fma rn f64 $r4d $r4d $r6d $r4d
+ sched 0x29 0x20 0x28 0x2e 0x00 0x00 0x00
+ // Multiply 2^27 to result for small inputs to recover
+ $p1 mul rn f64 $r4d $r4d 0x41a0000000000000
+ long mov b32 $r1 $r5
+ long mov b32 $r0 $r4
long ret
//
projects / mesa.git / commitdiff