Newbie to Forums, SIMD question
Peter Dove
I am working with Delphi and converting a lot of my Floating Point
code to assembler due to the lack of optimization done my Delphi in
that area. I am new to assembler but picking it up fast. I have a
peice of code that I want to convert to SIMD to get even more speed. I
was wondering if anyone would do the conversion to SIMD for me so I
can understand what is going on. I have included it as the current
code along with the commented out Delphi code. Any comments as to
speeding things up or how to convert to SIMD would be wonderful.
Thanks in anticipation
Peter Dove
Current Code
asm
//gray := ( val.Red + val.blue + val.Green ) * preThird;
fld dword ptr[val.red];
fadd dword ptr[val.blue];
fadd dword ptr[val.green];
fmul dword ptr[prethird];
fstp dword ptr[gray]
//Load up Saturation and Gray into regsters
fld dword ptr[a];
fld dword ptr[b];
fld dword ptr[gray];
fld dword ptr[saturation];
//Red := Gray+(((val.Red-Gray)*(saturation)));
fld dword ptr[val.red];
fsub st(0),st(2);
fmul st(0),st(1);
fadd st(0),st(2);
//If val.Red > 1 then val.Red := 1;
fcom st(3);
fstsw ax;
sahf;
FCMovNB st(0),st(3);
//If val.Red < 0 then val.Red := 0;
fcom st(4);
fstsw ax;
sahf;
FCMovB st(0),st(4);
fstp dword ptr[val.red];
//val.Green := val.Gray+(((val.Green-Gray)*(saturation)));
fld dword ptr[val.Green];
fsub st(0),st(2);
fmul st(0),st(1);
fadd st(0),st(2);
//If val.Green > 1 then val.Green := 1;
fcom st(3);
fstsw ax;
sahf;
FCMovNB st(0),st(3);
//If val.Green < 0 then val.Green := 0;
fcom st(4);
fstsw ax;
sahf;
FCMovB st(0),st(4);
fstp dword ptr[val.Green];
//val.Blue := Gray+(((val.Blue-Gray)*(saturation)));
fld dword ptr[val.Blue];
fsub st(0),st(2);
fmul st(0),st(1);
fadd st(0),st(2);
//If val.Blue > 1 then val.Blue := 1;
fcom st(3);
fstsw ax;
sahf;
FCMovNB st(0),st(3);
//If val.Blue < 0 then val.Blue := 0;
fcom st(4);
fstsw ax;
sahf;
FCMovB st(0),st(4);
fstp dword ptr[val.Blue];
ffree st(4);
ffree st(3);
ffree st(2);
ffree st(1);
ffree st(0);
end;
Robert Wessel
This is pretty easy with SSE.
Assume that you load the following values into XMM registers:
xmm1 = {0,0,0,0}
xmm2 = {1,1,1,1)
xmm3 = {gray,gray,gray,0} ;assumes you've computed gray
xmm4 = {red,green,blue,0}
xmm5 = {sat,sat,sat,0}
(And I'll leave you the task of shuffling stuff around to fit in the right places).
Then execute something like:
subps xmm4,xmm3
mulps xmm4,xmm5
addps xmm4,xmm3
minps xmm4,xmm1
maxps xmm4,xmm2
And then xmm4 ought to have the scaled red, green and blue values.
Peter Dove
>
> Assume that you load the following values into XMM registers:
>
> xmm1 = {0,0,0,0}
> xmm2 = {1,1,1,1)
> xmm3 = {gray,gray,gray,0} ;assumes you've computed gray
> xmm4 = {red,green,blue,0}
> xmm5 = {sat,sat,sat,0}
>
> (And I'll leave you the task of shuffling stuff around to fit in the right places).
>
> Then execute something like:
> subps xmm4,xmm3
> mulps xmm4,xmm5
> addps xmm4,xmm3
> minps xmm4,xmm1
> maxps xmm4,xmm2
>
> And then xmm4 ought to have the scaled red, green and blue values.
Hi,
Thanks very much for your reply. I have now been able to shave off
some 150ms for the speed which was taking about 1 sec. I have one
further question.
One of the things which is taking the time is that I must do a Power
for every pixel in order to apply a gamma value. This is of course a
problem because its hard to efficiently mix code with SSE. So my
question is, is there any way we can combine the code in some way... I
have pasted an example of the code which is embedded into the main
loop. because the array of data is so large it makes no sense to make
2 passes, one for the power ( R by g by B ) and one for the rest of
the options. Any ideas?.... hope this makes sense.
Peter
// First I pasre the Delphi Code and then the conversion to ASM
{ If red > dblBlackPoint then
begin
red := Exp(preGamma * Ln(red)); //Power ( pout, preGamma );
red := (red - c) * PreLevel;
end
else
red := 0;}
//Now the asm for the above
fld dword ptr[val.red];
fcom dword ptr[dblBlackPoint];
fstsw ax;
sahf;
jbe @SPIsZero;
fldln2;
fxch st(1);
fyl2x;
fmul dword ptr[preGamma];
FLDL2E
FMUL
FLD ST(0)
FRNDINT
FSUB ST(1), ST
FXCH ST(1)
F2XM1
FLD1
FADD
FSCALE
FSTP ST(1)
fsub dword ptr[c];
fmul dword ptr[preLevel]
fstp dword ptr[val.red];
jmp @RedEnd;
@SPIsZero:
xor eax, eax;
mov [val.red], eax;
ffree st(0);
@RedEnd:
Hope you can help
Peter
Robert Wessel
You don't ask for much, do ya? ;-)
OK, my first suggestion is to figure out how to get the pow() function
out of your inner loop by rewriting your algorithm. That will likely
have a much bigger effect that a low-level optimization of its
evaluation.
If you can't...
Pow() is a dreadfully complex function to evaluate well. You've
already simplified it to exp(a*ln(b)) so I assume precision isn't much
of an issue here.
Both exp() and ln() can be evaluated with a polynomial approximation
(there are several choices, including the simple Taylor or Maclaurin
series). If you can restrict your input ranges, and can live with
some error in the low bits, it should be pretty straight forward to
come up with some fixed length polynomials that you can parallelize
for each color plane and execute with the vector instructions.
For example, you can evaluate exp(x) with the Taylor series:
exp(x)=1+x+(x**2/2!)+(x**3/3!)+(x**4/4!)...
You compute that by repeated multiplying the term by x, and then
dividing (invariably a multiply by the reciprocal) by the additional
constant and adding the new term to the sum (so each additional term
evaluated is just one multiplication, one division and one addition).
The critical issue is convergence. For example, if you can guarantee
that x is less than .5, evaluating the five terms above guarantees a
maximum error of a bit more than 1/3840 (which is the maximum value of
the next term). There are series for ln() as well (see any basic
Calculus text book, plus the example below*). So the idea is to come
up with a fixed polynomial for the overall expression, that achieves
your precision requirements, that you can then parallelize.
To do that it's important to understand the exact input ranges, and
the required precision. And to determine how to restrict the
requirements as much as possible, and matching those limits to various
series. For example, if you can guarantee that x is close to one (and
especially close to, but greater than one), there's a series
particularly well suited for computing ln(x). Or if you can instead
compute ln(1+x) with -1<x<1 (especially with x close to zero), there's
a good series for that.
And then you can often do various transformations to keep the range of
the inputs inside the range you require for convergence.
If you can set some of those limits, I can probably make some more
specific suggestions.
*(One) series for ln(x) (x>0):
t=(x-1)/(x+1)
ln(x)/2=t+(t**3/3)+(t**5/5)+(t**7/7)...
Note that this converges most quickly when x is close to one.
arjan de lumens
Color calculations? Then I will make the following assumptions:
- 'red' is in the range 0.0 to 1.0, never outside
- 'preGamma' has the same value for a very large number of pixels
- you can live with 8-12 bits of precision.
If these assumptions are correct, you can use a lookup table, which
should be much faster than either pow, exp/log or series expansions.
To generate the lookup table (C code for a fixed 1001-entry table,
I am not familiar with Delphi):
float GammaLUT[1000+1];
void initGammaLUT( float preGamma )
{
for(i=0;i<=1000;i++)
GammaLUT[i] = pow(i/1000.0, preGamma);
}
To actually use the table, use the following assembly code:
fld dword ptr [val.red]
fmul dword ptr [const_1000] ;; const_1000 = 1000
fistp dword ptr [temp_val] ;; temporary variable
mov eax,[temp_val]
mov eax,[GammaLUT + 4*eax]
mov [val.red], eax
This should leave you with a value of ~8-10 bits of precision.
If you need more precision, you can make the lookup table larger or
apply interpolation between the entries of the lookup table
(beware: if the table is too large, you will get cache misses when
accessing it, eliminating the performance advantage); if you need
more performance, interleave instructions for several table lookups.
As long as SSE lacks approximate exp and log instructions, this
is probably the fastest available method. However, if you can accept
approximations to the gamma value, you can run sequences such as
(preGamma = 2.0)
movaps xmm0, [your_color_value]
mulps xmm0, xmm0
movaps [your_color_value], xmm0
(preGamma = 2.25)
movaps xmm0, [your_color_value]
rsqrtps xmm1, xmm0 ;; xmm1 <= (color ^ -0.5)
mulps xmm0, xmm0 ;; xmm0 <= (color ^ 2)
rsqrtps xmm1, xmm1 ;; xmm1 <= (color ^ 0.25)
mulps xmm0, xmm1 ;; xmm0 <= (color ^ 2.25)
movaps [your_color_value], xmm0
Peter Dove
Thanks again, your help is very much appreciated, I wish I had the
time to spend some study of maths!
OK so the dblGamma is put in as 2.2 or 1.5 etc and then the following
calculation is made
preGamma := 1/dblGamma;
The value of any pixel ( red, green, or blue ) is between 0 and 1
And thats the crux of it.
Thanks
Peter
Peter Dove
Hi,
Thanks, thanks and thanks again. I was thinking already about LUT's
for this function but its nice to see the actual code.
Peter
Robert Wessel
OK, I spent some time thinking about this. You've got a somewhat
difficult problem since you don't have much constraint on the input
ranges. Anyway here goes.
We'll use two series to compute exp(y*ln(x)):
(1) exp(x)=1+x+(x**2/2!)+(x**3/3!)+(x**4/4!)…
(2) ln(x+1)=x-(x**2/2)+(x**3/3)-(x**4/4)…
Both of these converse most rapidly when x is near zero. So the main
complexity is in scaling the arguments to friendly values.
The following section will ignore zeros, denormals, infinities and
NaNs, all of which require special handling.
First we need to scale the input x to a value very close to one.
Ideally, if we can scale x so that |x-1|<1/128 (perhaps the bounds
need to be a bit tighter than that), we'll need to evaluate only three
terms of (2) to exceed to precision of a single (IOW, the fourth term
would be a factor of at least 2**23 smaller than the first). I'm
going to be a little less ambitious (I've not been thorough enough in
the error analysis) and go for four terms. If you don't mind loosing
an extra bit of two of precision, or you can improve the error
analysis, you may be able to drop the fourth term.
We'll scale the input x by using the relation:
ln(a*b*x)=ln(a)+ln(b)+ln(x)
We'll use three lookup tables to perform the initial scaling of x.
We'll compute ln(a) as ln(2**n), where n is the exponent of the input,
which we can just extract from the exponent field of the input x
floating point number:
mov ecx,vector_input_x[#]
mov ebx,ecx
shr ebx,21
and ebx,0x000003fc
mov eax,scale_ln_2_x[ebx]
mov vector_lna[#],eax
and ecx,0x807fffff
or ecx,0x40000000
mov vector_input_x[#],ecx
The table scale_ln_2_x contains 256 singles, each of which containing
the precomputed values ln(2*-128) through ln(2**127). The above code
then reduces the exponent field to zero (remember the +128 bias), thus
factoring out a 2**n term for ln(a). The precomputed ln(2**n) is
stored in an vector for SSE use later.
An example: let's say that the input exponent is 130 (+2 after the
bias). The value we'll pull out of the table will be ln(2**2), and
we'll change the exponent in the number to 128 (+0 with the bias). So
we've divided the original value by 2**2, and used the precomputed
lookup table to get the ln(a) value.
This leaves input x in the range 1<=x<2. We'll now reduce than by
using the six leading bits of the mantissa as an index into a second
lookup table to compute ln(b).
mov ebx,ecx
shr ebx,15
and ebx,0x000000fc
mov eax,scale_ln_1_plus_64th[ebx]
mov vector_lnb[#],eax
mov eax,scale_1_over_1_plus_64th[ebx]
mov vector_scale_lnb[#],eax
The 64 entry lookup table (which is split into two sections for
convenience) consists of pairs of precomputed values for the points:
(1+0/64+1/128), (1+1/64+1/128), (1+2/64+1/128)… (1+63/64+1/128). The
two precomputed values are ln(b) (for scale_ln_1_plus_64th), and 1/x
(for scale_1_over_1_plus_64th).
We store both values in vectors for later use. The ln(b) value is
used as described above, and the reciprocal value is used to finally
scale x to our desired range (we use the reciprocal to avoid a slow
division).
Again an example: Assume that the six initial bits of mantissa are
110000. So the number, remembering the implied one bit, is 1.110000…
(binary), or someplace between 1+48/64 and 1+49/64. We'll split the
difference, and use 1+48/64+1/128 to get a scaling factor the will
reduce that number to the range 127/128..129/128, or within about
1/128th of one. The lookup tables will result in the values
ln(1+48/64+1/128) (aka ln(b)) and 1/(1+48/64+1/128).
This gets repeated for each of the inputs to the vector portion (up to
four), with the updated (partially scaled) x, and the three values
from the lookup tables stored in four vectors.
Next we load up the partially scaled x's, and finish the scaling
process:
movaps xmm1,vector_input_x
mulps xmm1,vector_scale_lnb
Now we compute four terms of (2):
subps xmm1,{1,1,1,1}
movaps xmm2,xmm1 ;sum=x
movaps xmm3,xmm1
mulps xmm3,xmm1 ;x**2
movaps xmm4,{1/2,1/2,1/2,1/2)
mulps xmm4,xmm3 ;x**2/2
subps xmm2,xmm4
mulps xmm3,xmm1 ;x**3
movaps xmm4,{1/3,1/3,1/3,1/3)
mulps xmm4,xmm3 ;x**3/3
subps xmm2,xmm4
mulps xmm3,xmm1 ;x**4
movaps xmm4,{1/4,1/4,1/4,1/4)
mulps xmm4,xmm3 ;x**4/4
subps xmm2,xmm4
Xmm2 now contains ln(scaled_x), to which we then add ln(a) and ln(b),
and then multiply by y:
addps xmm2,vector_lna
addps xmm2,vector_lnb
mulps xmm2,vector_input_y
Xmm2 now contains y*ln(x) for four inputs.
Now we're going to scale the partial results (which we'll call z from
here on) to a magnitude less than 1/64 with the relation e**(a+b+x) =
(e**a)*(e**b)*(e**x). We'll take advantage of the fact that the
partial result cannot be larger than 89 without driving the final
result to infinity (or zero) – IOW, e**89 cannot be represented in a
single.
We'll start by limiting the range:
mov xmm1,{-90,-90,-90,-90}
maxps xmm2,xmm1
mov xmm1,{90,90,90,90}
minps xmm2,xmm1
We then compute the integer n, where n*(1/64) is largest value still
smaller (in magnitude) than z, and then reduce z by that amount.
movaps xmm3,xmm2
mulps xmm3,{64,64,64,64} ;z*64
cvttss2si eax,xmm3 ;convert low single to integer (w/truncation)
cvtsi2ss xmm3,eax ;and back to single
mov vector_expscale+0,eax ;and save scaling factor
shufps xmm3,xmm3,0x39 ;rotate xmm3 right one slot
cvttss2si eax,xmm3 ;repeat for 2nd single (now in low slot)
cvtsi2ss xmm3,eax
mov vector_expscale+4,eax
shufps xmm3,xmm3,0x39
cvttss2si eax,xmm3 ;repeat for 3rd single (now in low slot)
cvtsi2ss xmm3,eax
mov vector_expscale+8,eax
shufps xmm3,xmm3,0x39
cvttss2si eax,xmm3 ;repeat for 4th single (now in low slot)
cvtsi2ss xmm3,eax
mov vector_expscale+12,eax
shufps xmm3,xmm3,0x39 ;now back in "normal" order
mulps xmm3,(1/64,1/64,1/64,1/64} ;and scale back to 1/64ths
subps xmm2,xmm3 ;z scaled to < 1/64
Now we compute e**expscale. Expscale is an integer in the range of
(–90*64) to (+90*64). We split the computation into the two pieces,
one with the integer part and one with the fractional part, and then
use a lookup table to compute e**a and e**b. For example, if we had
the value –3 5/64, we'd compute e**(-3) and e**(-5/64).
mov eax,expscale[#]
mov ebx,eax ;-90*64..90*64
sar ebx,6 ;-90..90 (a)
mov edx,ebx
add ebx,90 ;0..180
mov ecx,exp_int[ebx*4] ;exp(a)
mov vector_expa[#],ecx
sal edx,6 ;-90*64..90*64 rounded to exclude 64ths
sub eax,edx ;-63..63 (b)
add eax,63
mov ecx,exp_64ths[eax*4] ;exp(b)
mov vector_expb[#],ecx
The two lookup tables contain the values e**(-90)..e**90 (for exp_int)
and e**(-63/64)..e**(63/64) (for exp_64ths). The first contains 181
entries and the second 127. Note that the first and last couple of
entries in the exp_int table should be zero and +Inf, respectively.
This gets repeated for each of the four singles in the vector.
Finally we compute exp(scaled_z) using (1) Again, the convergence
allows us to limit evaluation to five terms total (counting the "1" as
the first term).
movaps xmm1,{1,1,1,1}
movaps xmm3,xmm2 ;x
addps xmm1,xmm3
mulps xmm3,xmm2 ;x**2
mulps xmm3,{1/2,1/2,1/2,1/2} ;x**2/2!
addps xmm1,xmm3
mulps xmm3,xmm2 ;x**3/2!
mulps xmm3,{1/3,1/3,1/3,1/3} ;x**3/3!
addps xmm1,xmm3
mulps xmm3,xmm2 ;x**4/3!
mulps xmm3,{1/4,1/4,1/4,1/4} ;x**4/4!
addps xmm1,xmm3
And finally we factor in exp(a) and exp(b):
mulps xmm1,vector_expa
mulps xmm1,vector_expb
And, we're done, with four parallel evaluations of pow(x,y) in xmm1.
As always, debugging and testing is left to the reader.
There are a total of five tables containing 692 precomputed singles
(2768 bytes), plus a few work fields, in about 169 instructions
without any branches, so this should all cache and predict nicely.
Compared to four of the original calculations, this should be rather
faster. A fair bit of optimization should be possible to expose ILP
with the usual rearrangements and register unrolling. There are also
a couple of spots where some SSE2 instructions would be useful. I'm
estimating that this should get you a solid 17 bits of precision in
the results, perhaps as many as 19.
Are you sorry you asked? ;-)
So far we've ignored zeros, denormals, infinities and NaNs. In short,
the above will not work if either operand is a zero, denormal,
infinity or NaN, nor will it work for a negative x. At the stage
where the initial inputs are prepared you'll need to insert code to
handle those cases. You can detect denormals and zero because they
have an exponent field of zero, and infinities and NaNs have an
exponent of 255. Negative x's are detected by inspecting the sign
bit.
Some values can be scaled to ones within range, while others just need
to be computed as special cases. Also note that I'm going to ignore
the different types of NaN (SNaN, QnaN, Indefinites, etc.), and assume
they're all the same. If you need to distinguish between them, you'll
have to muck through it yourself.
The special cases should be handled by creating a mask vector along
with a vector containing all the special-case precomputed results (in
their correct slots). The "normal" calculation must still proceed for
those slots, so forcing in (1,1) is an easy way to avoid problems.
The mask vector should contain 0xffffffff in each slot corresponding
to the presence of a precomputed result. At the end of the above
process, the precomputed results can be merged with a "real" results
with:
movaps xmm2,vector_mask
andnps xmm1,xmm2
andps xmm2,vector_precomputed
orps xmm1,xmm2
I'm not going to go through all the combinations of special cases you
need to handle. I suggest you make a table with all the special cases
for x on one axis, and the special cases for y on the other, and make
sure you've covered all the combinations. Anyway, some rules:
X may not be negative, and should return a NaN.
If either x or y is a NaN, the result should be a NaN.
Denormal x's should be scaled to a non-denormal range by taking the
square root, and then halving y, and running the usual calculation.
Denormal y's should be treated as zero, and return a 1, unless x is
+Infinity or 1, in which case 0, 1 or infinity should be returned, as
appropriate, based on the sign of y.
An x of zero should always return 0, a y of zero should always return
1, pow(0,0) should return 0, 1 or NaN, as appropriate for your
application.
An infinite x should return zero or infinity depending on the sign of
y, or one (probably) if y is zero.
A +Infinity y should return 1 or +Infinity depending on whether x is
above or below one.
A –Infinity y should return 0 or 1 depending on whether x is above or
below one.
I'm sure I've missed one or two cases. Also be aware that these can
combine in "interesting" ways (consider, for example, all the cases
where x and y are various combinations of zeros and infinities), which
is why you should construct a table to make sure you catch all the
cases. If you choose to distinguish between types of NaNs, you'll
have to expand things rather further.
Fun, fun, fun!