| ## Public API Reference

- Main Page
- Modules
- Namespaces
- Classes
- Files
- Related Pages

- File List
- File Members

# csqsqrt.h

Go to the documentation of this file.

    00001 /*
    00002     Fast computation of sqrt(x) and 1/sqrt(x)
    00003     Copyright (C) 2002 by Matthew Reda <reda@mac.com> (PowerPC version)
    00004
    00005     This library is free software; you can redistribute it and/or
    00006     modify it under the terms of the GNU Library General Public
    00007     License as published by the Free Software Foundation; either
    00008     version 2 of the License, or (at your option) any later version.
    00009
    00010     This library is distributed in the hope that it will be useful,
    00011     but WITHOUT ANY WARRANTY; without even the implied warranty of
    00012     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    00013     Library General Public License for more details.
    00014
    00015     You should have received a copy of the GNU Library General Public
    00016     License along with this library; if not, write to the Free
    00017     Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
    00018 */
    00019
    00026 #ifndef __CS_CSQSQRT_H__
    00027 #define __CS_CSQSQRT_H__
    00028
    00029 #include "cssysdef.h"
    00030 #include <math.h>
    00031
    00039 static CS_FORCEINLINE float csQsqrt (float x);
    00040
    00048 static CS_FORCEINLINE float csQisqrt (float x);
    00049
    00052 #if !defined(CS_NO_QSQRT) && \
    00053     defined(CS_PROCESSOR_POWERPC) && defined(CS_COMPILER_GCC)
    00054
    00055 /*
    00056  * Use the PowerPC fsqrte to get an estimate of 1/sqrt(x) Then apply two
    00057  * Newton-Rhaphson refinement steps to get a more accurate response Finally
    00058  * multiply by x to get x/sqrt(x) = sqrt(x).  Add additional refinement steps
    00059  * to get a more accurate result.  Zero is treated as a special case, otherwise
    00060  * we end up returning NaN (Not a Number).
    00061  */
    00062 static CS_FORCEINLINE float csQsqrt(float x)
    00063 {
    00064   float y0 = 0.0;
    00065
    00066   if (x != 0.0)
    00067   {
    00068     float x0 = x * 0.5f;
    00069
    00070     __asm__ __volatile__ ("frsqrte %0,%1" : "=f" (y0) : "f" (x));
    00071
    00072     y0 = y0 * (1.5f - x0 * y0 * y0);
    00073     y0 = (y0 * (1.5f - x0 * y0 * y0)) * x;
    00074   }
    00075
    00076   return y0;
    00077 }
    00078
    00079 /*
    00080  * Similar to csQsqrt() above, except we do not multiply by x at the end, and
    00081  * return 1/sqrt(x).
    00082  */
    00083 static inline float csQisqrt(float x)
    00084 {
    00085   float x0 = x * 0.5f;
    00086   float y0;
    00087   __asm__ __volatile__ ("frsqrte %0,%1" : "=f" (y0) : "f" (x));
    00088
    00089   y0 = y0 * (1.5f - x0 * y0 * y0);
    00090   y0 = y0 * (1.5f - x0 * y0 * y0);
    00091
    00092   return y0;
    00093 }
    00094
    00095 #else
    00096
    00097 static CS_FORCEINLINE float csQsqrt (float x) { return sqrtf(x); }
    00098 static CS_FORCEINLINE float csQisqrt(float x) { return 1.0f / sqrtf(x); }
    00099
    00100 #endif
    00101
    00102 #endif // __CS_CSQSQRT_H__