src/nouveau/codegen/lib/gk110.asm - third_party/mesa - Git at Google

 .section #gk110_builtin_code
 // DIV U32
 //
 // UNR recurrence (q = a / b):
 // look for z such that 2^32 - b <= b * z < 2^32
 // then q - 1 <= (a * z) / 2^32 <= q
 //
 // INPUT:   $r0: dividend, $r1: divisor
 // OUTPUT:  $r0: result, $r1: modulus
 // CLOBBER: $r2 - $r3, $p0 - $p1
 // SIZE:    22 / 14 * 8 bytes
 //
 gk110_div_u32:
    sched 0x28 0x04 0x28 0x04 0x28 0x28 0x28
    bfind u32 $r2 $r1
    xor b32 $r2 $r2 0x1f
    mov b32 $r3 0x1
    shl b32 $r2 $r3 clamp $r2
    cvt u32 $r1 neg u32 $r1
    mul $r3 u32 $r1 u32 $r2
    add $r2 (mul high u32 $r2 u32 $r3) $r2
    sched 0x28 0x28 0x28 0x28 0x28 0x28 0x28
    mul $r3 u32 $r1 u32 $r2
    add $r2 (mul high u32 $r2 u32 $r3) $r2
    mul $r3 u32 $r1 u32 $r2
    add $r2 (mul high u32 $r2 u32 $r3) $r2
    mul $r3 u32 $r1 u32 $r2
    add $r2 (mul high u32 $r2 u32 $r3) $r2
    mul $r3 u32 $r1 u32 $r2
    sched 0x04 0x28 0x04 0x28 0x28 0x2c 0x04
    add $r2 (mul high u32 $r2 u32 $r3) $r2
    mov b32 $r3 $r0
    mul high $r0 u32 $r0 u32 $r2
    cvt u32 $r2 neg u32 $r1
    add $r1 (mul u32 $r1 u32 $r0) $r3
    set $p0 0x1 ge u32 $r1 $r2
    $p0 sub b32 $r1 $r1 $r2
    sched 0x28 0x2c 0x04 0x20 0x2e 0x28 0x20
    $p0 add b32 $r0 $r0 0x1
    $p0 set $p0 0x1 ge u32 $r1 $r2
    $p0 sub b32 $r1 $r1 $r2
    $p0 add b32 $r0 $r0 0x1
    ret

 // DIV S32, like DIV U32 after taking ABS(inputs)
 //
 // INPUT:   $r0: dividend, $r1: divisor
 // OUTPUT:  $r0: result, $r1: modulus
 // CLOBBER: $r2 - $r3, $p0 - $p3
 //
 gk110_div_s32:
    set $p2 0x1 lt s32 $r0 0x0
    set $p3 0x1 lt s32 $r1 0x0 xor $p2
    sched 0x20 0x28 0x28 0x04 0x28 0x04 0x28
    cvt s32 $r0 abs s32 $r0
    cvt s32 $r1 abs s32 $r1
    bfind u32 $r2 $r1
    xor b32 $r2 $r2 0x1f
    mov b32 $r3 0x1
    shl b32 $r2 $r3 clamp $r2
    cvt u32 $r1 neg u32 $r1
    sched 0x28 0x28 0x28 0x28 0x28 0x28 0x28
    mul $r3 u32 $r1 u32 $r2
    add $r2 (mul high u32 $r2 u32 $r3) $r2
    mul $r3 u32 $r1 u32 $r2
    add $r2 (mul high u32 $r2 u32 $r3) $r2
    mul $r3 u32 $r1 u32 $r2
    add $r2 (mul high u32 $r2 u32 $r3) $r2
    mul $r3 u32 $r1 u32 $r2
    sched 0x28 0x28 0x04 0x28 0x04 0x28 0x28
    add $r2 (mul high u32 $r2 u32 $r3) $r2
    mul $r3 u32 $r1 u32 $r2
    add $r2 (mul high u32 $r2 u32 $r3) $r2
    mov b32 $r3 $r0
    mul high $r0 u32 $r0 u32 $r2
    cvt u32 $r2 neg u32 $r1
    add $r1 (mul u32 $r1 u32 $r0) $r3
    sched 0x2c 0x04 0x28 0x2c 0x04 0x28 0x20
    set $p0 0x1 ge u32 $r1 $r2
    $p0 sub b32 $r1 $r1 $r2
    $p0 add b32 $r0 $r0 0x1
    $p0 set $p0 0x1 ge u32 $r1 $r2
    $p0 sub b32 $r1 $r1 $r2
    $p0 add b32 $r0 $r0 0x1
    $p3 cvt s32 $r0 neg s32 $r0
    sched 0x04 0x2e 0x28 0x04 0x28 0x28 0x28
    $p2 cvt s32 $r1 neg s32 $r1
    ret

 // RCP F64
 //
 // INPUT:   $r0d
 // OUTPUT:  $r0d
 // CLOBBER: $r2 - $r9, $p0
 //
 // The core of RCP and RSQ implementation is Newton-Raphson step, which is
 // used to find successively better approximation from an imprecise initial
 // value (single precision rcp in RCP and rsqrt64h in RSQ).
 //
 gk110_rcp_f64:
    // 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 subtracting 1 from the exponent, which will
    // mean that it's > 0x7fd in those cases when doing unsigned comparison
    set b32 $p0 0x1 gt u32 $r3 0x7fd
    // $r3: 0 for norms, 0x36 for denorms, -1 for others
    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 b32 $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 b32 $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
    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
    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 f64 neg $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
    add b32 $r4 $r2 $r3
    ext u32 $r3 $r1 0xb14
    // New exponent in $r3
    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 b32 $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
    add b32 $r1 $r4 $r1
    bra #rcp_end
 rcp_result_inf_or_denorm:
    // New exponent >= 0x7ff means that result is inf
    set b32 $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
    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 b32 $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:
    ret

 // RSQ F64
 //
 // INPUT:   $r0d
 // OUTPUT:  $r0d
 // CLOBBER: $r2 - $r9, $p0 - $p1
 //
 gk110_rsq_f64:
    // 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 b32 $p1 0x1 le u32 $r3 0x2
    or b32 $r2 $r0 $r2
    $p1 mul rn f64 $r0d $r0d 0x4350000000000000
    rsqrt64h f32 $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
    and b32 $r2 $r2 $r6
    sched 0x28 0x2b 0x20 0x27 0x28 0x2e 0x28
    set b32 $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
    mov b32 $r0 0x0
    or b32 $r1 $r1 $r5
    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}
    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
    mov b32 $r1 $r5
    mov b32 $r0 $r4
    ret

 .section #gk110_builtin_offsets
 .b64 #gk110_div_u32
 .b64 #gk110_div_s32
 .b64 #gk110_rcp_f64
 .b64 #gk110_rsq_f64
	.section #gk110_builtin_code
	// DIV U32
	//
	// UNR recurrence (q = a / b):
	// look for z such that 2^32 - b <= b * z < 2^32
	// then q - 1 <= (a * z) / 2^32 <= q
	//
	// INPUT: $r0: dividend, $r1: divisor
	// OUTPUT: $r0: result, $r1: modulus
	// CLOBBER: $r2 - $r3, $p0 - $p1
	// SIZE: 22 / 14 * 8 bytes
	//
	gk110_div_u32:
	sched 0x28 0x04 0x28 0x04 0x28 0x28 0x28
	bfind u32 $r2 $r1
	xor b32 $r2 $r2 0x1f
	mov b32 $r3 0x1
	shl b32 $r2 $r3 clamp $r2
	cvt u32 $r1 neg u32 $r1
	mul $r3 u32 $r1 u32 $r2
	add $r2 (mul high u32 $r2 u32 $r3) $r2
	sched 0x28 0x28 0x28 0x28 0x28 0x28 0x28
	mul $r3 u32 $r1 u32 $r2
	add $r2 (mul high u32 $r2 u32 $r3) $r2
	mul $r3 u32 $r1 u32 $r2
	add $r2 (mul high u32 $r2 u32 $r3) $r2
	mul $r3 u32 $r1 u32 $r2
	add $r2 (mul high u32 $r2 u32 $r3) $r2
	mul $r3 u32 $r1 u32 $r2
	sched 0x04 0x28 0x04 0x28 0x28 0x2c 0x04
	add $r2 (mul high u32 $r2 u32 $r3) $r2
	mov b32 $r3 $r0
	mul high $r0 u32 $r0 u32 $r2
	cvt u32 $r2 neg u32 $r1
	add $r1 (mul u32 $r1 u32 $r0) $r3
	set $p0 0x1 ge u32 $r1 $r2
	$p0 sub b32 $r1 $r1 $r2
	sched 0x28 0x2c 0x04 0x20 0x2e 0x28 0x20
	$p0 add b32 $r0 $r0 0x1
	$p0 set $p0 0x1 ge u32 $r1 $r2
	$p0 sub b32 $r1 $r1 $r2
	$p0 add b32 $r0 $r0 0x1
	ret

	// DIV S32, like DIV U32 after taking ABS(inputs)
	//
	// INPUT: $r0: dividend, $r1: divisor
	// OUTPUT: $r0: result, $r1: modulus
	// CLOBBER: $r2 - $r3, $p0 - $p3
	//
	gk110_div_s32:
	set $p2 0x1 lt s32 $r0 0x0
	set $p3 0x1 lt s32 $r1 0x0 xor $p2
	sched 0x20 0x28 0x28 0x04 0x28 0x04 0x28
	cvt s32 $r0 abs s32 $r0
	cvt s32 $r1 abs s32 $r1
	bfind u32 $r2 $r1
	xor b32 $r2 $r2 0x1f
	mov b32 $r3 0x1
	shl b32 $r2 $r3 clamp $r2
	cvt u32 $r1 neg u32 $r1
	sched 0x28 0x28 0x28 0x28 0x28 0x28 0x28
	mul $r3 u32 $r1 u32 $r2
	add $r2 (mul high u32 $r2 u32 $r3) $r2
	mul $r3 u32 $r1 u32 $r2
	add $r2 (mul high u32 $r2 u32 $r3) $r2
	mul $r3 u32 $r1 u32 $r2
	add $r2 (mul high u32 $r2 u32 $r3) $r2
	mul $r3 u32 $r1 u32 $r2
	sched 0x28 0x28 0x04 0x28 0x04 0x28 0x28
	add $r2 (mul high u32 $r2 u32 $r3) $r2
	mul $r3 u32 $r1 u32 $r2
	add $r2 (mul high u32 $r2 u32 $r3) $r2
	mov b32 $r3 $r0
	mul high $r0 u32 $r0 u32 $r2
	cvt u32 $r2 neg u32 $r1
	add $r1 (mul u32 $r1 u32 $r0) $r3
	sched 0x2c 0x04 0x28 0x2c 0x04 0x28 0x20
	set $p0 0x1 ge u32 $r1 $r2
	$p0 sub b32 $r1 $r1 $r2
	$p0 add b32 $r0 $r0 0x1
	$p0 set $p0 0x1 ge u32 $r1 $r2
	$p0 sub b32 $r1 $r1 $r2
	$p0 add b32 $r0 $r0 0x1
	$p3 cvt s32 $r0 neg s32 $r0
	sched 0x04 0x2e 0x28 0x04 0x28 0x28 0x28
	$p2 cvt s32 $r1 neg s32 $r1
	ret

	// RCP F64
	//
	// INPUT: $r0d
	// OUTPUT: $r0d
	// CLOBBER: $r2 - $r9, $p0
	//
	// The core of RCP and RSQ implementation is Newton-Raphson step, which is
	// used to find successively better approximation from an imprecise initial
	// value (single precision rcp in RCP and rsqrt64h in RSQ).
	//
	gk110_rcp_f64:
	// 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 subtracting 1 from the exponent, which will
	// mean that it's > 0x7fd in those cases when doing unsigned comparison
	set b32 $p0 0x1 gt u32 $r3 0x7fd
	// $r3: 0 for norms, 0x36 for denorms, -1 for others
	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 b32 $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 b32 $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
	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
	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 f64 neg $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
	add b32 $r4 $r2 $r3
	ext u32 $r3 $r1 0xb14
	// New exponent in $r3
	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 b32 $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
	add b32 $r1 $r4 $r1
	bra #rcp_end
	rcp_result_inf_or_denorm:
	// New exponent >= 0x7ff means that result is inf
	set b32 $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
	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 b32 $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:
	ret

	// RSQ F64
	//
	// INPUT: $r0d
	// OUTPUT: $r0d
	// CLOBBER: $r2 - $r9, $p0 - $p1
	//
	gk110_rsq_f64:
	// 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 b32 $p1 0x1 le u32 $r3 0x2
	or b32 $r2 $r0 $r2
	$p1 mul rn f64 $r0d $r0d 0x4350000000000000
	rsqrt64h f32 $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
	and b32 $r2 $r2 $r6
	sched 0x28 0x2b 0x20 0x27 0x28 0x2e 0x28
	set b32 $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
	mov b32 $r0 0x0
	or b32 $r1 $r1 $r5
	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}
	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
	mov b32 $r1 $r5
	mov b32 $r0 $r4
	ret

	.section #gk110_builtin_offsets
	.b64 #gk110_div_u32
	.b64 #gk110_div_s32
	.b64 #gk110_rcp_f64
	.b64 #gk110_rsq_f64