CrystalSpace

Public API Reference

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__

Generated for Crystal Space 1.2.1 by doxygen 1.5.3