I wrote a implementation in MMIX:

          % Constants
  FISRCON GREG #5FE6EB50C7B537A9
  THREHAF GREG #3FF8000000000000
          % Save half of the original number
          OR $2,$0,0
          INCH $2,#FFF0
          % Bit level hacking
          SRU $1,$0,1
          SUBU $0,FISRCON,$1
          % First iteration
          FMUL $1,$2,$0
          FMUL $1,$1,$0
          FSUB $1,THREHAF,$1
          FMUL $0,$0,$1
          % Second iteration
          FMUL $1,$2,$0
          FMUL $1,$1,$0
          FSUB $1,THREHAF,$1
          FMUL $0,$0,$1

This implementation makes an assumption that the original number is greater than 2^-1021.