Fast float-to-int formula with rounding toward -Inf: how does it work?

Question

I've been playing around with x87 FPU programming,and I just came across the following magic spell for converting from float (32-bit) to int (32-bit):

    flds    from             # number to convert from
    fadds   magic1           # add magic number 1 (???)
    fstps   to               # store in destination as single-precision binary32
    movl    to,%eax          # load float result as an int
    subl    $magic2,%eax     # subtract magic number 2 (???)
    sarl    $1,%eax          # divide by 2
    movl    %eax,to          # and look, here's the result!

 .section .rodata
magic1:  .float 6291456.0    # 3 * 2^21
.equ magic2, 0x4ac00000      # ???

   .data
from:    .float 31415.9265   # pick a number, any number...
to:      .long  0            # result ends up here

(AT&T syntax with GAS directives)

I have tried this out, and it seems to work (rounding towards -infinity) but I have absolutely no idea why! Can anyone explain how it works?

Answer 1

Quick answer:

The first magical number forces rescaling of the argument so that the LSbit of the integer part becomes rightmost-but-one in the fraction.

Then adding the second magical will erase the exponent bits.

Regarding the necessity of the last division by 2, IDK, there must be an extra technicality (probably related to the addition of 3.2^21 rather than 1.2^21).