rhpvorderman
commented
9 months ago Never mind what I said about the speed. I was looking at loop unrolling. Specifically at how I can write a loop such that the compiler unrolls it. Then I found one of the best answers on stackoverflow ever: https://stackoverflow.com/a/36591776/16437839
This answer notes that when there is a dependency in a loop, GCC does not unroll. As example an addition was given, which is exactly the calculation here. In order to leverage out of order execution it is best to assign the calculation to different variables to calculate intermediate results and sum them afterwards.
So I did, and the result is faster than my previous float implementation. So it is both faster and more accurate. How awesome is that? I am elated!
marcelm
Hi Marcel,
It seems I made an error of judgement in #726. In that thread I state that single precision numbers is enough. That statement is incorrect, or a half-truth.
Since phred scores are between 0 and 93 they have only two significant decimal digits. Two decimal digits means 100 units, so that can be approached by a floating point number that has a significand of at least 7 bits significance. Floating point has 23 bits, so that is ample.
However this does not properly look at the calculation required to add two phred units. Single precision numbers have a 23-bit significand meaning that they can represent numbers that are 2^23 ~= 8 million units apart. Adding a phred of 0 to a phred of 70 (10 million units) is therefore barely possible. Adding a phred of 93 to a phred of 0 simply does not register. That means that floating point numbers can never be sufficient as an intermediate when summing large numbers of phred probabilities as some phreds simply get lost.
This is still true when looking at more practical use cases. In practice the highest illumina phred is 37. So 10**-3.7 verses 1.0 that is roughly 13 bits. Combined with a read length of 151 (9 bits) that makes for 22bits of precision needed to ensure all probabilities. Float has that. Going beyond illumina to nanopore however it immediately becomes clear that there aren't enough bits as nanopore can reach sizes up to 2 million (21 bits) and has a phred range between 0 and 50ish (17 bits).
A double precision number with a 53-bit significand is ample however for most practical use cases. I have sharpened the test case again so that it is correct even when using the extremes of the phred range.
Speed does not mean anything when the result is not correct. Luckily the vectorization can be kept and expected errors is still vectorized and still faster than 4.4.