/opengles/src/fixed.cpp

  1// ==========================================================================
  2//
  3// fixed.cpp	Implementation of Fixed Point Arithmetic
  4//
  5//		If Intel Graphics Performance Primitives (GPP) are available,
  6//		the functions defined in this header should be mapped to
  7//		the assembler primitives provided in this library.
  8//
  9//		Fixed point numbers are represented in 16.16 format, where
 10//		the high 16 bits represent the signed integer part, while the
 11//		lower 16 bits represent the fractional part of a number.
 12//
 13// --------------------------------------------------------------------------
 14//
 15// Copyright (C) 2003  David Blythe   All Rights Reserved.
 16//
 17// Permission is hereby granted, free of charge, to any person obtaining a
 18// copy of this software and associated documentation files (the "Software"),
 19// to deal in the Software without restriction, including without limitation
 20// the rights to use, copy, modify, merge, publish, distribute, sublicense,
 21// and/or sell copies of the Software, and to permit persons to whom the
 22// Software is furnished to do so, subject to the following conditions:
 23//
 24// The above copyright notice and this permission notice shall be included
 25// in all copies or substantial portions of the Software.
 26//
 27// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 28// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 29// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
 30// DAVID BLYTHE BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN
 31// AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
 32// CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 33//
 34// ==========================================================================
 35
 36
 37#include "stdafx.h"
 38#include "fixed.h"
 39
 40
 41#ifndef EGL_USE_GPP
 42
 43
 44// --------------------------------------------------------------------------
 45// Calculate inverse of fixed point number
 46//
 47// Parameters:
 48//	a		-	the number whose inverse should be calculated
 49// --------------------------------------------------------------------------
 50EGL_Fixed EGL_Inverse(EGL_Fixed a) {
 51    I32 exp;
 52    EGL_Fixed x;
 53    /* 1/(4x) */
 54    static const I32 __gl_rcp_tab[] = { /* domain 0.5 .. 1.0-1/16 */
 55		0x8000, 0x71c7, 0x6666, 0x5d17, 0x5555, 0x4ec4, 0x4924, 0x4444
 56    };
 57    if (a == EGL_ZERO) return 0x7fffffff;
 58	bool sign = false;
 59
 60	if (a < 0) {
 61		sign = true;
 62		a = -a;
 63	}
 64
 65#ifdef EGL_USE_CLZ
 66	exp = _CountLeadingZeros(a);
 67#else
 68    x = a;
 69    exp = 31;
 70    if (x & 0xffff0000) { exp -= 16; x >>= 16; }
 71    if (x & 0xff00) { exp -= 8; x >>= 8; }
 72    if (x & 0xf0) { exp -= 4; x >>= 4; }
 73    if (x & 0xc) { exp -= 2; x >>= 2; }
 74    if (x & 0x2) { exp -= 1; }
 75#endif
 76    x = __gl_rcp_tab[(a>>(28-exp))&0x7]<<2;
 77    exp -= 16;
 78    if (exp <= 0)
 79		x >>= -exp;
 80    else
 81		x <<= exp;
 82//printf("est %f\n", __GL_F_2_FLOAT(x));
 83    /* two iterations of newton-raphson  x = x(2-ax) */
 84	x = EGL_Mul(x,(EGL_ONE*2 - EGL_Mul(a,x)));
 85	x = EGL_Mul(x,(EGL_ONE*2 - EGL_Mul(a,x)));
 86//printf("recip %f %f\n", __GL_F_2_FLOAT(a), __GL_F_2_FLOAT(x));
 87
 88	if (sign)
 89		return -x;
 90	else
 91		return x;
 92}
 93
 94
 95// --------------------------------------------------------------------------
 96// Calculate the inverse of the square root of fixed point number
 97//
 98// Parameters:
 99//	a		-	the numbers whose inverse of square root should be 
100//				calculated
101// --------------------------------------------------------------------------
102EGL_Fixed EGL_InvSqrt(EGL_Fixed a) {
103    EGL_Fixed x /*= (a+EGL_ONE)>>1*//*a>>1*/;
104    /* 1/(2sqrt(x)) - extra divide by 2 to scale to 16 bits */
105    static const U16 __gl_rsq_tab[] = { /* domain 0.5 .. 1.0-1/16 */
106		0xb504, 0xaaaa, 0xa1e8, 0x9a5f, 0x93cd, 0x8e00, 0x88d6, 0x8432,
107    };
108    I32 i, exp;
109    if (a == EGL_ZERO) return 0x7fffffff;
110    if (a == EGL_ONE) return a;
111
112#ifdef EGL_USE_CLZ
113	exp = _CountLeadingZeros(a);
114#else
115    x = a;
116    exp = 31;
117    if (x & 0xffff0000) { exp -= 16; x >>= 16; }
118    if (x & 0xff00) { exp -= 8; x >>= 8; }
119    if (x & 0xf0) { exp -= 4; x >>= 4; }
120    if (x & 0xc) { exp -= 2; x >>= 2; }
121    if (x & 0x2) { exp -= 1; }
122#endif
123    x = __gl_rsq_tab[(a>>(28-exp))&0x7]<<1;
124//printf("t %f %x %f %d %d\n", __GL_F_2_FLOAT(a), a, __GL_F_2_FLOAT(x), exp, 28-exp);
125    exp -= 16;
126    if (exp <= 0)
127		x >>= -exp>>1;
128    else
129		x <<= (exp>>1)+(exp&1);
130    if (exp&1) x = EGL_Mul(x, __gl_rsq_tab[0]);
131//printf("nx %f\n", __GL_F_2_FLOAT(x));
132
133    /* newton-raphson */
134    /* x = x/2*(3-(a*x)*x) */
135    i = 0;
136    do {
137//printf("%f\n", __GL_F_2_FLOAT(x));
138		x = EGL_Mul((x>>1),(EGL_ONE*3 - EGL_Mul(EGL_Mul(a,x),x)));
139    } while(++i < 3);
140//printf("rsqrt %f %f\n", __GL_F_2_FLOAT(a), __GL_F_2_FLOAT(x));
141    return x;
142}
143
144static const U16 __gl_sin_tab[] = {
145#include "gl_sin.h"
146};
147
148// --------------------------------------------------------------------------
149// Calculate cosine of fixed point number
150//
151// Parameters:
152//	a		-	the numbers whose cosine should be calculated
153// --------------------------------------------------------------------------
154EGL_Fixed EGL_Cos(EGL_Fixed a) {
155    EGL_Fixed v;
156    /* reduce to [0,1) */
157    while (a < EGL_ZERO) a += EGL_2PI;
158    a *= EGL_R2PI;
159    a >>= EGL_PRECISION;
160    a += 0x4000;
161
162    /* now in the range [0, 0xffff], reduce to [0, 0xfff] */
163    a >>= 4;
164
165    v = (a & 0x400) ? __gl_sin_tab[0x3ff - (a & 0x3ff)] : __gl_sin_tab[a & 0x3ff];
166    v = EGL_Mul(v,EGL_ONE);
167    return (a & 0x800) ? -v : v;
168}
169
170// --------------------------------------------------------------------------
171// Calculate sine of fixed point number
172//
173// Parameters:
174//	value		-	the numbers whose sine should be calculated
175// --------------------------------------------------------------------------
176EGL_Fixed EGL_Sin(EGL_Fixed a) {
177    EGL_Fixed v;
178
179    /* reduce to [0,1) */
180    while (a < EGL_ZERO) a += EGL_2PI;
181    a *= EGL_R2PI;
182    a >>= EGL_PRECISION;
183
184    /* now in the range [0, 0xffff], reduce to [0, 0xfff] */
185    a >>= 4;
186
187    v = (a & 0x400) ? __gl_sin_tab[0x3ff - (a & 0x3ff)] : __gl_sin_tab[a & 0x3ff];
188    v = EGL_Mul(v,EGL_ONE);
189    return (a & 0x800) ? -v : v;
190}
191
192// --------------------------------------------------------------------------
193// Calculate square root of fixed point number
194//
195// Parameters:
196//	value		-	the number whose square root should be calculated
197// --------------------------------------------------------------------------
198EGL_Fixed EGL_Sqrt(EGL_Fixed a) {
199    EGL_Fixed s;
200    I32 i;
201    s = (a + EGL_ONE) >> 1;
202    /* 6 iterations to converge */
203    for (i = 0; i < 6; i++)
204		s = (s + EGL_Div(a, s)) >> 1;
205    return s;
206}
207
208#endif //ndef EGL_USE_GPP
209
210/* assume 0 <= x <= 1 and y >= 0 */
211static __inline EGL_Fixed
212xpow(EGL_Fixed x, EGL_Fixed y) {
213    I32 exp;
214    EGL_Fixed p, f = x;
215//#define __POW_LOW_PREC
216    /* absolute error < 0.009, more than good enough for 5-bit color */
217    static const U16 __gl_log_tab[] = { /* domain .5 - 1.0 */
218		0xffff, 0xd47f, 0xad96, 0x8a62, 0x6a3f, 0x4caf, 0x3151, 0x17d6, 0x0000
219    };
220    static const U16 __gl_alog_tab[] = { /* domain 0 - 1.0 */
221		0xffff, 0xeac0, 0xd744, 0xc567, 0xb504, 0xa5fe, 0x9837, 0x8b95, 0x8000
222    };
223
224    if (y == EGL_ZERO || x == EGL_ONE)
225		return EGL_ONE;
226    if (x == EGL_ZERO)
227		return EGL_ZERO;
228
229    exp = 15;
230    if (f & 0xff00) { exp -= 8; f >>= 8; }
231    if (f & 0xf0) { exp -= 4; f >>= 4; }
232    if (f & 0xc) { exp -= 2; f >>= 2; }
233    if (f & 0x2) { exp -= 1; }
234    f = x << exp;
235//printf("f/x/exp %f %f %d %x\n", __GL_F_2_FLOAT(f), __GL_F_2_FLOAT(x), exp, x);
236    /* compute log base 2 */
237    x = (f&0x0fff)<<4;
238    f = (f >> 12)&0x7;
239    p = EGL_Mul(__gl_log_tab[f+1]-__gl_log_tab[f],x) + __gl_log_tab[f];
240//printf("p %f x %f\n", __GL_F_2_FLOAT(p), __GL_F_2_FLOAT(x));
241    p = EGL_Mul(p, y) + y*exp;
242//printf("np %f \n", __GL_F_2_FLOAT(p));
243
244    /* compute antilog (2^p) */
245    exp = EGL_IntFromFixed(p);
246    p = EGL_FractionFromFixed(p);
247//printf("exp = %d frac %f\n", exp, __GL_F_2_FLOAT(p));
248    x = (p&0x1fff)<<3;
249    p = p >> 13;
250    return (EGL_Mul(__gl_alog_tab[p+1]-__gl_alog_tab[p],x)+
251		__gl_alog_tab[p])>>exp;
252}
253
254// --------------------------------------------------------------------------
255// Exponentiation function
256//
257// Parameters:
258//	value		-	the basis; assume 0 <= value <= 1
259//	exponent	-	the exponent, exponent >= 0
260// --------------------------------------------------------------------------
261OGLES_API EGL_Fixed EGL_Power(EGL_Fixed a, EGL_Fixed n) {
262#if (defined(ARM) || defined(_ARM_))
263    return xpow(a, n);
264#else
265	return EGL_FixedFromFloat(pow(EGL_FloatFromFixed(a), EGL_FloatFromFixed(n)));
266#endif
267}
268