`Core.Intrinsics.fpext` cannot extend Float16 to Float16

[eschnett]

With Julia nightly:

julia> Core.Intrinsics.fpext(Float16, Float16(3.3))
ERROR: fpext: output bitsize must be >= input bitsize
Stacktrace:
 [1] top-level scope
   @ REPL[3]:1
 [2] top-level scope
   @ ~/.julia/packages/Metal/N2ABH/src/initialization.jl:79

This works for other types (Float32, Float64), and should be working according to the error message.

[oscardssmith]

relatedly:

julia> Base.fptrunc(Float16, Float16(2.1))
Float16(0.0)

julia> Base.fptrunc(Float16, Float32(2.1))
Float16(2.1)

[eschnett]

@oscardssmith That's with Julia 1.11 I assume? Julia nightly doesn't have Base.fptrunc.

The same problem is there for BFloat16.

[oscardssmith]

@eschnett are you sure?

julia/src/intrinsics.h

Line 68 in 090e291

ADD_I(fptrunc, 2) \

It really looks like it's still there...

[eschnett]

Silly me, it's just a warning, it's still there.

[vchuravy]

I looked at this, but it isn't easy to fix.

Both intrinsics are defined here:

julia/src/runtime_intrinsics.c

Lines 1632 to 1657 in 88c71dd

	#define fptrunc(tr, pr, a) \
	if (!(osize < 8 * sizeof(a))) \
	jl_error("fptrunc: output bitsize must be < input bitsize"); \
	else if (osize == 16) { \
	if ((jl_datatype_t*)tr == jl_float16_type) \
	*(uint16_t*)pr = float_to_half(a); \
	else /*if ((jl_datatype_t*)tr == jl_bfloat16_type)*/ \
	*(uint16_t*)pr = float_to_bfloat(a); \
	} \
	else if (osize == 32) \
	*(float*)pr = a; \
	else if (osize == 64) \
	*(double*)pr = a; \
	else \
	jl_error("fptrunc: runtime floating point intrinsics are not implemented for bit sizes other than 16, 32 and 64");
	#define fpext(tr, pr, a) \
	if (!(osize >= 8 * sizeof(a))) \
	jl_error("fpext: output bitsize must be >= input bitsize"); \
	if (osize == 32) \
	*(float*)pr = a; \
	else if (osize == 64) \
	*(double*)pr = a; \
	else \
	jl_error("fpext: runtime floating point intrinsics are not implemented for bit sizes other than 32 and 64");
	un_fintrinsic_withtype(fptrunc,fptrunc)
	un_fintrinsic_withtype(fpext,fpext)

Thanks to C macros it is a bit hard to know what pr and a is in this context.

So we have to look at un_fintrinisc_withtype:

julia/src/runtime_intrinsics.c

Lines 1059 to 1067 in 88c71dd

	#define un_fintrinsic_withtype(OP, name) \
	un_fintrinsic_bfloat(OP, jl_##name##bf16) \
	un_fintrinsic_half(OP, jl_##name##16) \
	un_fintrinsic_ctype(OP, jl_##name##32, float) \
	un_fintrinsic_ctype(OP, jl_##name##64, double) \
	JL_DLLEXPORT jl_value_t *jl_##name(jl_value_t *ty, jl_value_t *a) \
	{ \
	return jl_fintrinsic_1(ty, a, #name, jl_##name##bf16, jl_##name##16, jl_##name##32, jl_##name##64); \
	}

This leads to:

julia/src/runtime_intrinsics.c

Lines 1091 to 1096 in 88c71dd

	case 2:
	if (jl_typeof(a) == (jl_value_t*)jl_float16_type)
	halfop(sz2 * host_char_bit, ty, pa, pr);
	else /*if (jl_typeof(a) == (jl_value_t*)jl_bfloat16_type)*/
	bfloatop(sz2 * host_char_bit, ty, pa, pr);
	break;

and

julia/src/runtime_intrinsics.c

Lines 764 to 776 in 88c71dd

	#define un_fintrinsic_half(OP, name) \
	static inline void name(unsigned osize, jl_value_t *ty, void *pa, void *pr) JL_NOTSAFEPOINT \
	{ \
	uint16_t a = *(uint16_t*)pa; \
	float A = half_to_float(a); \
	if (osize == 16) { \
	float R; \
	OP(ty, &R, A); \
	*(uint16_t*)pr = float_to_half(R); \
	} else { \
	OP(ty, (uint16_t*)pr, A); \
	} \
	}

So for half intrinsics we always convert a to float?

Which is why this is false for fpext:

if (!(osize >= 8 * sizeof(a))) \ 
             jl_error("fpext: output bitsize must be >= input bitsize"); \ 

But confusingly, when osize == 16 then we are also converting pr to float?

So, I don't know how the code for fptrunc works?

         else if (osize == 16) { \ 
             if ((jl_datatype_t*)tr == jl_float16_type) \ 
                 *(uint16_t*)pr = float_to_half(a); \ 
             else /*if ((jl_datatype_t*)tr == jl_bfloat16_type)*/ \ 
                 *(uint16_t*)pr = float_to_bfloat(a); \ 
         } \ 

Maybe @maleadt remembers enough to fix this?

[xal-0]

The control flow is somewhat confusing. For both fpext and fptrunc it works like this:

• jl_fintrinsic_1 looks at the type provided as the first argument, allocates space for a result of that size, and then calls the corresponding function jl_fptrunc{b16,16,32,64} generated by un_fintrinsic_{bfloat,half,ctype,ctype}. osize gets the size of the result type in bits, ty gets the type of the float argument.
• un_fintrinsic_half converts from float16 to float32, and the macro is expanded with tr getting the type of the float input, pr getting a pointer to the result, and a getting the float32 version of the input. osize is in scope in the macro body.
• fpext and fptrunc write to pr assuming it expects a float16 when osize == 16, but this is not true, because un_fintrinsic_half assume a float32 result and convert it back to a float16.
• Additionally, fptrunc confuses tr for the type of the result, but it is in fact the type of the float argument.

So the result is fptrunc only hits the lines with float_to_half/float_to_bfloat when both the input and result types are some from of float16, and it writes the result into the lowest 16 bits of a float32, which get truncated away by un_fintrinsic_half, leaving whatever was originally in the memory allocated for the result to be interpreted as the upper half of a float32.

When truncating from float32 or float64 to float16, it works properly because un_fintrinsic_ctype doesn't touch the result.

I'm working on reorganizing some of runtime_intrinsics.c to make this less confusing.