New test fails where long double is not long enough

[ginggs]

The following new test, introduced in 3c1205c, fails on architectures where char is unsigned by default (e.g. ARM, POWER).
test_mpileup($opts,in=>[qw(indel-AD.1)],out=>'mpileup/indel-AD.1.out',ref=>'indel-AD.1.fa',args=>q[-a AD]);

Diff of the output below:

--- test/mpileup/indel-AD.1.out
+++ test/mpileup/indel-AD.1.out.new
@@ -154,7 +154,7 @@
 000000F        524     .       T       G,<*>   0       .       DP=119;I16=64,44,1,0,4253,180647,27,729,6480,388800,60,3600,2460,59296,25,625;QS=0.993531,0.00646862,0;SGB=-0.379885;RPB=1;MQB=1;MQSB=1;BQB=1;MQ0F=0       PL:AD   0,255,255,255,255,255:108,1,0
 000000F        525     .       A       <*>     0       .       DP=119;I16=63,39,0,0,4126,179608,0,0,6120,367200,0,0,2330,56198,0,0;QS=1,0;MQSB=1;MQ0F=0        PL:AD   0,255,255:102,0
 000000F        526     .       T       <*>     0       .       DP=120;I16=62,45,0,0,4234,181460,0,0,6420,385200,0,0,2450,58910,0,0;QS=1,0;MQSB=1;MQ0F=0        PL:AD   0,255,255:107,0
-000000F        527     .       C       <*>     0       .       DP=120;I16=64,40,0,0,4275,186993,0,0,6240,374400,0,0,2382,57234,0,0;QS=1,0;MQSB=1;MQ0F=0        PL:AD   0,255,255:104,0
+000000F        527     .       C       <*>     0       .       DP=120;I16=64,40,0,0,4275,186993,0,0,6240,374400,0,0,2382,57234,0,0;QS=1,0;MQSB=1;MQ0F=0        PL:AD   0,255,0:104,0
 000000F        528     .       A       <*>     0       .       DP=121;I16=62,45,0,0,4346,189994,0,0,6420,385200,0,0,2426,58030,0,0;QS=1,0;MQSB=1;MQ0F=0        PL:AD   0,255,255:107,0
 000000F        529     .       T       <*>     0       .       DP=123;I16=64,43,0,0,4397,193617,0,0,6420,385200,0,0,2439,58361,0,0;QS=1,0;MQSB=1;MQ0F=0        PL:AD   0,255,255:107,0
 000000F        530     .       G       <*>     0       .       DP=123;I16=66,42,0,0,4392,192268,0,0,6480,388800,0,0,2486,59400,0,0;QS=1,0;MQSB=1;MQ0F=0        PL:AD   0,255,255:108,0

Any suggestions where I should be looking for the cause of the differing output, please?

[pd3]

The PLs are calculated here https://github.com/samtools/bcftools/blob/develop/bam2bcf.c#L249, errmod_cal is now part of htslib.

[ginggs]

@pd3 thanks for the pointer. I was able to dig a bit further and found the problem lies here:
https://github.com/samtools/htslib/blob/develop/errmod.c#L90
sum1 / sum tends towards 0/0 and beta[k] is set to NaN at least for the first 13 values of k when n=104 and q=37 (which triggers the differing test output above).

It seems to be purely coincidence that this fails on the same architectures where char is unsigned.

[ginggs]

Fixed in samtools/htslib#617