español - sqrt c++ sintaxis
Punto fijo rápido pow, log, exp y sqrt (4)
A continuación se muestra un ejemplo de implementación en C del algoritmo de base 2 de registro de punto fijo de Clay S. Turner [ 1 ]. El algoritmo no requiere ningún tipo de tabla de consulta. Esto puede ser útil en sistemas donde las restricciones de memoria son limitadas y el procesador carece de una FPU, como es el caso de muchos microcontroladores. La base de registro e y la base de registro 10 también son compatibles mediante el uso de la propiedad de logaritmos que, para cualquier base n :
log (x)
y
log (x) = _______
n log (n)
y
donde, para este algoritmo, y es igual a 2.
Una buena característica de esta implementación es que admite precisión variable: la precisión se puede determinar en tiempo de ejecución, a expensas del rango. De la forma en que lo he implementado, el procesador (o compilador) debe ser capaz de hacer cálculos de 64 bits para mantener algunos resultados intermedios. Puede adaptarse fácilmente para no requerir soporte de 64 bits, pero el rango se reducirá.
Cuando se usan estas funciones, se espera que x
sea un valor de punto fijo escalado de acuerdo con la precision
especificada. Por ejemplo, si la precision
es 16, entonces se debe escalar x
en 2 ^ 16 (65536). El resultado es un valor de punto fijo con el mismo factor de escala que la entrada. Un valor de retorno de INT32_MIN
representa infinito negativo. Un valor de retorno de INT32_MAX
indica un error y errno
se establecerá en EINVAL
, lo que indica que la precisión de entrada no fue válida.
#include <errno.h>
#include <stddef.h>
#include "log2fix.h"
#define INV_LOG2_E_Q1DOT31 UINT64_C(0x58b90bfc) // Inverse log base 2 of e
#define INV_LOG2_10_Q1DOT31 UINT64_C(0x268826a1) // Inverse log base 2 of 10
int32_t log2fix (uint32_t x, size_t precision)
{
int32_t b = 1U << (precision - 1);
int32_t y = 0;
if (precision < 1 || precision > 31) {
errno = EINVAL;
return INT32_MAX; // indicates an error
}
if (x == 0) {
return INT32_MIN; // represents negative infinity
}
while (x < 1U << precision) {
x <<= 1;
y -= 1U << precision;
}
while (x >= 2U << precision) {
x >>= 1;
y += 1U << precision;
}
uint64_t z = x;
for (size_t i = 0; i < precision; i++) {
z = z * z >> precision;
if (z >= 2U << precision) {
z >>= 1;
y += b;
}
b >>= 1;
}
return y;
}
int32_t logfix (uint32_t x, size_t precision)
{
uint64_t t;
t = log2fix(x, precision) * INV_LOG2_E_Q1DOT31;
return t >> 31;
}
int32_t log10fix (uint32_t x, size_t precision)
{
uint64_t t;
t = log2fix(x, precision) * INV_LOG2_10_Q1DOT31;
return t >> 31;
}
El código para esta implementación también se encuentra en Github , junto con un programa de muestra / prueba que ilustra cómo usar esta función para calcular y mostrar logaritmos de números leídos de la entrada estándar.
[ 1 ] CS Turner, 1 , IEEE Signal Processing Mag. , pp. 124,140, sep. 2010.
Tengo una clase de punto fijo (10.22) y necesito un pow, un sqrt, una exp y una función de registro.
Lamentablemente, no tengo ni idea de por dónde empezar. ¿Puede alguien proporcionarme algunos enlaces a artículos útiles o, mejor aún, proporcionarme un código?
Supongo que una vez que tengo una función exp, se vuelve relativamente fácil implementar pow y sqrt como se han convertido.
pow (x, y) => exp (y * log (x)) sqrt (x) => pow (x, 0.5)
Esas son las funciones exp y log que me resultan difíciles (como si recordara algunas de mis reglas de registro, no puedo recordar mucho más sobre ellas).
Es de suponer que, por cierto, también habría un método más rápido para sqrt y pow, por lo que cualquier puntero en ese frente sería apreciado incluso si es solo para usar los métodos que describo anteriormente :)
Tenga en cuenta: esto TIENE que ser multiplataforma y en código C / C ++ puro, por lo que no puedo usar ninguna optimización de ensamblador.
Un buen punto de partida es el libro de Jack Crenshaw, "Math Toolkit for Real-Time Programming" . Tiene una buena discusión de algoritmos e implementaciones para varias funciones trascendentales.
Una solución muy simple es usar una aproximación decente basada en tablas. En realidad, no necesita una gran cantidad de datos si reduce sus entradas correctamente. exp(a)==exp(a/2)*exp(a/2)
, lo que significa que realmente solo necesita calcular exp(x)
para 1 < x < 2
. En ese rango, una aproximación de runga-kutta daría resultados razonables con ~ 16 entradas IIRC.
De manera similar, sqrt(a) == 2 * sqrt(a/4)
que significa que solo necesita entradas de tabla para 1 < a < 4
. El registro (a) es un poco más difícil: log(a) == 1 + log(a/e)
. Esta es una iteración bastante lenta, pero el registro (1024) es solo de 6.9, por lo que no tendrá muchas iteraciones.
Usaría un algoritmo "entero primero" similar para pow: pow(x,y)==pow(x, floor(y)) * pow(x, frac(y))
. Esto funciona porque pow(double, int)
es trivial (divide y vencerás).
[editar] Para el componente integral del log(a)
, puede ser útil almacenar una tabla 1, e, e^2, e^3, e^4, e^5, e^6, e^7
para que puede reducir el log(a) == n + log(a/e^n)
mediante una simple búsqueda binaria codificada en un archivo en esa tabla. La mejora de 7 a 3 pasos no es tan grande, pero significa que solo tiene que dividir una vez por e^n
lugar de n
veces por e
.
[edit 2] Y para el último término de log(a/e^n)
, puedes usar log(a/e^n) = log((a/e^n)^8)/8
- cada iteración produce 3 más bits por tabla de búsqueda . Eso mantiene su código y tamaño de tabla pequeños. Normalmente, este es el código para sistemas integrados y no tienen cachés grandes.
[editar 3] Eso no es muy inteligente de mi lado. log(a) = log(2) + log(a/2)
. Solo puede almacenar el valor de punto fijo log2=0.30102999566
, contar el número de ceros iniciales, cambiar a
en el rango utilizado para su tabla de búsqueda y multiplicar ese cambio (entero) por la constante de punto fijo log2
. Puede ser tan bajo como 3 instrucciones.
El uso de e
para el paso de reducción solo le da un "buen" log(e)=1.0
constante, pero eso es una falsa optimización. 0.30102999566 es una constante tan buena como 1.0; ambas son constantes de 32 bits en 10.22 puntos fijos. Usar 2 como la constante para la reducción de rango le permite usar un cambio de bit para una división.
Aún obtienes el truco de la edición 2, log(a/2^n) = log((a/2^n)^8)/8
. Básicamente, esto le da un resultado (a + b/8 + c/64 + d/512) * 0.30102999566
- con b, c, d en el rango [0,7]. a.bcd
realmente es un número octal. No es una sorpresa ya que usamos 8 como potencia. (El truco funciona igualmente bien con el poder 2, 4 o 16).
[editar 4] Todavía tenía un final abierto. pow(x, frac(y)
es solo pow(sqrt(x), 2 * frac(y))
y tenemos un 1/sqrt(x)
decente. Eso nos da un enfoque mucho más eficiente. Diga frac(y)=0.101
binario, es decir, 1/2 más 1/8. Entonces eso significa que x^0.101
es (x^1/2 * x^1/8)
1/2 (x^1/2 * x^1/8)
1/8 (x^1/2 * x^1/8)
. Pero x^1/2
1/2 es solo sqrt(x)
x^1/8
es (sqrt(sqrt(sqrt(x)))
. Al guardar una operación más, Newton-Raphson NR(x)
nos da 1/sqrt(x)
por lo que calculamos 1.0/(NR(x)*NR((NR(NR(x)))
. Solo invertimos el resultado final, no usamos la función sqrt directamente.
Verifique la implementación de sqrt de mi punto fijo utilizando solo operaciones de enteros. Fue divertido inventar. Muy viejo ahora.
De lo contrario, compruebe el conjunto de algoritmos CORDIC . Esa es la forma de implementar todas las funciones que enumeró y las funciones trigonométricas.
EDITAR: He publicado la fuente revisada en GitHub here