Unexpected loss of precision for powers of two

[stephane-chazelas]

IEEE 754 based float, double or long double should store powers of 2 exactly. Yet, here on Debian testing amd64:

$ ksh -c 'printf "%.30f\n" 0x1p{30..50}'
1073741824.000000000000000000000000000000
2147483648.000000000001519229000000000000
4294967296.000000000003038458000000000000
8589934592.000000000006076917000000000000
17179869184.000000000012153833000000000000
34359738368.000000000024307667000000000000
68719476736.000000000048615334000000000000
137438953472.000000000097230668000000000000
274877906944.000000000194461336000000000000
549755813888.000000000388922672000000000000
1099511627776.000000000777845344000000000000
2199023255552.000000001555690687000000000000
4398046511104.000000003111381375000000000000
8796093022208.000000006222762750000000000000
17592186044416.000000012445525499000000000000
35184372088832.000000024891050998000000000000
70368744177664.000000049782101996000000000000
140737488355328.000000099564203992000000000000
281474976710656.000000199128407985000000000000
562949953421312.000000398256815970000000000000
1125899906842624.000000796513631940000000000000

It's worse with typeset -F:

$ ksh -c 'typeset -F30 i; for i in 0x1p{30..50}; do print "$i"; done'
1073741824.000000000000000000000000000000
2147483648.000000044930857000000000000000
4294967296.000000089861715000000000000000
8589934592.000000179723429000000000000000
17179869184.000000359446858000000000000000
34359738368.000000718893716000000000000000
68719476736.000001437787432000000000000000
137438953472.000002875574864000000000000000
274877906944.000005751149729000000000000000
549755813888.000011502299458000000000000000
1099511627776.000023004598916000000000000000
2199023255552.000046009197831000000000000000
4398046511104.000092018395662000000000000000
8796093022208.000184036791325000000000000000
17592186044416.000368073582649000000000000000
35184372088832.000736147165298000000000000000
70368744177664.001472294330597000000000000000
140737488355328.002944588661194000000000000000
281474976710656.005889177322388000000000000000
562949953421312.011778354644775000000000000000
1125899906842624.023556709289551000000000000000

Converting that last one back to 0xHHHp+n notation:

$ ksh -c 'printf "%a\n" 0x1p50 1125899906842624.023556709289551000000000000000'
0x1.0000000000000000000000000000p+50
0x1.0000000000000182000000000000p+50

Other shells are fine for numbers up to 2¹⁰²³ as expected:

$ bash -c 'printf "%.500g\n" 0x1p1023'
89884656743115795386465259539451236680898848947115328636715040578866337902750481566354238661203768010560056939935696678829394884407208311246423715319737062188883946712432742638151109800623047059726541476042502884419075341171231440736956555270413618581675255342293149119973622969239858152417678164812112068608

Comparing with arbitrary precision:

$ perl -Mbigint -E 'say 1<<1023'
89884656743115795386465259539451236680898848947115328636715040578866337902750481566354238661203768010560056939935696678829394884407208311246423715319737062188883946712432742638151109800623047059726541476042502884419075341171231440736956555270413618581675255342293149119973622969239858152417678164812112068608

$ echo '2^1023' | BC_LINE_LENGTH=0 bc
89884656743115795386465259539451236680898848947115328636715040578866337902750481566354238661203768010560056939935696678829394884407208311246423715319737062188883946712432742638151109800623047059726541476042502884419075341171231440736956555270413618581675255342293149119973622969239858152417678164812112068608