SQRT using FPU

[AlfaOmega08]

Sometimes ago I found a small asm inline code for gcc.
I've unfortunately lost it.
I need it to "double sqrt(double x);"...

Does anyone know how to do it?
Else, I've to implement the Newton or Babilonian method :'(

Thanks in advance...

[JamesM]

I thought there was a sqrt instruction... did you check the manuals?

[AlfaOmega08]

I remember there are lots of fpu instructions like this:
fsqrt
fsin
fcos
ftan
[...]
But you have to store and read values in fpu registers, and there is specific inline asm to do that

[lukem95]

the linux kernel has a whole section in arch/i386/math-emu

but it looks complex as hell:

Code: Select all

static void fsqrt_(FPU_REG *st0_ptr)
{
  char st0_tag = st0_ptr->tag;

  clear_C1();
  if ( !(st0_tag ^ TW_Valid) )
    {
      int expon;
      
      if (st0_ptr->sign == SIGN_NEG)
	{
	  arith_invalid(st0_ptr);  /* sqrt(negative) is invalid */
	  return;
	}

#ifdef DENORM_OPERAND
      if ( (st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
	return;
#endif DENORM_OPERAND

      expon = st0_ptr->exp - EXP_BIAS;
      st0_ptr->exp = EXP_BIAS + (expon & 1);  /* make st(0) in  [1.0 .. 4.0) */
      
      wm_sqrt(st0_ptr, control_word);	/* Do the computation */
      
      st0_ptr->exp += expon >> 1;
      st0_ptr->sign = SIGN_POS;
    }
  else if ( st0_tag == TW_Zero )
    return;
  else if ( st0_tag == TW_Infinity )
    {
      if ( st0_ptr->sign == SIGN_NEG )
	arith_invalid(st0_ptr);  /* sqrt(-Infinity) is invalid */
      return;
    }
  else
    { single_arg_error(st0_ptr); return; }

}

[nick8325]

glibc uses the following (public domain) code:

Code: Select all

double sqrt (double x)
{
  double res;
  asm ("fsqrt" : "=t" (res) : "0" (x));
  return res;
}

[AlfaOmega08]

Thanks nick
It's exactly wath was I looking for

[bluecode]

the linux kernel has a whole section in arch/i386/math-emu

Well, from the pathname it should be obvious, that the code you presented tries to emulate the FPU behaviour (if no FPU is present or if it is deactivated through the bit in Cr0).

[lukem95]

i did wonder now you mention it, i thought it had something to do with... i don't know actually.

it may have been useful anyway

[Solar]

Note that sqrt() is defined to either set errno = EDOM or to raise an invalid floating point exception in case that the parameter is < 0. This can be quite tricky to get right and efficient.