example - math.h en dev c++

¿Cómo puedo escribir una función de potencia yo mismo? (13)

A ^B = Log ^-1 (Log (A) * B)

Editar: sí, esta definición realmente proporciona algo útil. Por ejemplo, en un x86, se traduce casi directamente a FYL2X (Y * Log ₂ (X)) y F2XM1 (2 ^x -1):

fyl2x fld st(0) frndint fsubr st(1),st fxch st(1) fchs f2xmi fld1 faddp st(1),st fscale fstp st(1)

El código termina un poco más de lo esperado, principalmente porque F2XM1 solo funciona con números en el rango -1.0..1.0. El fld st(0)/frndint/fsubr st(1),st pieza resta la parte entera, por lo que nos queda solo la fracción. Aplicamos F2XM1 a eso, agregamos el 1 de nuevo, luego usamos FSCALE para manejar la parte entera de la exponenciación.

Siempre me pregunté cómo puedo hacer una función que calcula la potencia (por ejemplo, 2 ³ ) yo mismo. En la mayoría de los idiomas, estos están incluidos en la biblioteca estándar, principalmente como pow(double x, double y) , pero ¿cómo puedo escribirlo yo mismo?

Estaba pensando for loops , pero creo que mi cerebro se metió en un bucle (cuando quería hacer una potencia con un exponente no entero, como 5 ^4.5 o negativos 2 ^-21 ) y me volví loco;)

Entonces, ¿cómo puedo escribir una función que calcule el poder de un número real? Gracias

Ah, tal vez sea importante tener en cuenta: no puedo usar funciones que usen poderes (por ejemplo, exp ), lo que lo haría inútil en última instancia.

Es un ejercicio interesante. Aquí hay algunas sugerencias, que debes probar en este orden:

Use un bucle
Use recursión (no mejor, pero interesante de todos modos)
Optimice enormemente su recursividad utilizando técnicas de divide y vencerás
Usa logaritmos

Estoy usando aritmética de punto fijo largo y mi pow está basado en log2 / exp2. Los números consisten en:

int sig = { -1; +1 } int sig = { -1; +1 } signum
DWORD a[A+B]
A es el número de DWORD s para la parte entera del número
B es el número de DWORD para la parte fraccionaria

Mi solución simplificada es esta:

//--------------------------------------------------------------------------- longnum exp2 (const longnum &x) { int i,j; longnum c,d; c.one(); if (x.iszero()) return c; i=x.bits()-1; for(d=2,j=_longnum_bits_b;j<=i;j++,d*=d) if (x.bitget(j)) c*=d; for(i=0,j=_longnum_bits_b-1;i<_longnum_bits_b;j--,i++) if (x.bitget(j)) c*=_longnum_log2[i]; if (x.sig<0) {d.one(); c=d/c;} return c; } //--------------------------------------------------------------------------- longnum log2 (const longnum &x) { int i,j; longnum c,d,dd,e,xx; c.zero(); d.one(); e.zero(); xx=x; if (xx.iszero()) return c; //**** error: log2(0) = infinity if (xx.sig<0) return c; //**** error: log2(negative x) ... no result possible if (d.geq(x,d)==0) {xx=d/xx; xx.sig=-1;} i=xx.bits()-1; e.bitset(i); i-=_longnum_bits_b; for (;i>0;i--,e>>=1) // integer part { dd=d*e; j=dd.geq(dd,xx); if (j==1) continue; // dd> xx c+=i; d=dd; if (j==2) break; // dd==xx } for (i=0;i<_longnum_bits_b;i++) // fractional part { dd=d*_longnum_log2[i]; j=dd.geq(dd,xx); if (j==1) continue; // dd> xx c.bitset(_longnum_bits_b-i-1); d=dd; if (j==2) break; // dd==xx } c.sig=xx.sig; c.iszero(); return c; } //--------------------------------------------------------------------------- longnum pow (const longnum &x,const longnum &y) { //x^y = exp2(y*log2(x)) int ssig=+1; longnum c; c=x; if (y.iszero()) {c.one(); return c;} // ?^0=1 if (c.iszero()) return c; // 0^?=0 if (c.sig<0) { c.overflow(); c.sig=+1; if (y.isreal()) {c.zero(); return c;} //**** error: negative x ^ noninteger y if (y.bitget(_longnum_bits_b)) ssig=-1; } c=exp2(log2(c)*y); c.sig=ssig; c.iszero(); return c; } //---------------------------------------------------------------------------

dónde:

_longnum_bits_a = A*32 _longnum_bits_b = B*32 _longnum_log2[i] = 2 ^ (1/(2^i)) ... precomputed sqrt table _longnum_log2[0]=sqrt(2) _longnum_log2[1]=sqrt[tab[0]) _longnum_log2[i]=sqrt(tab[i-1]) longnum::zero() sets *this=0 longnum::one() sets *this=+1 bool longnum::iszero() returns (*this==0) bool longnum::isnonzero() returns (*this!=0) bool longnum::isreal() returns (true if fractional part !=0) bool longnum::isinteger() returns (true if fractional part ==0) int longnum::bits() return num of used bits in number counted from LSB longnum::bitget()/bitset()/bitres()/bitxor() are bit access longnum.overflow() rounds number if there was a overflow X.FFFFFFFFFF...FFFFFFFFF??h -> (X+1).0000000000000...000000000h int longnum::geq(x,y) is comparition |x|,|y| returns 0,1,2 for (<,>,==)

Todo lo que necesitas para entender este código es que los números en forma binaria consisten en la suma de potencias de 2, cuando necesitas calcular 2 ^ num entonces se puede reescribir como esto

2^(b(-n)*2^(-n) + ... + b(+m)*2^(+m))

donde n son bits fraccionarios m son bits enteros. la multiplicación / división por 2 en forma binaria es un simple cambio de bit así que si lo pones todo junto obtienes un código para exp2 similar a mi. log2 se basa en la búsqueda de binaru ... cambiando los bits de resultado de MSB a LSB hasta que coincida con el valor buscado (algoritmo muy similar al del cálculo rápido de sqrt). Espero que esto ayude a aclarar las cosas ...

Hay dos casos distintos para tratar: exponentes de enteros y exponentes fraccionarios.

Para exponentes enteros, puedes usar exponenciación cuadrando.

def pow(base, exponent): if exponent == 0: return 1 elif exponent < 0: return 1 / pow(base, -exponent) elif exponent % 2 == 0: half_pow = pow(base, exponent // 2) return half_pow * half_pow else: return base * pow(base, exponent - 1)

El segundo "elif" es lo que lo distingue de la función pow ingenua. Permite que la función realice O (log n) llamadas recursivas en lugar de O (n).

Para exponentes fraccionarios, puede usar la identidad a ^ b = C ^ (b * log_C (a)). Es conveniente tomar C = 2, entonces a ^ b = 2 ^ (b * log2 (a)). Esto reduce el problema a funciones de escritura para 2 ^ x y log2 (x).

La razón por la que es conveniente tomar C = 2 es que los números de punto flotante se almacenan en punto flotante base-2. log2 (a * 2 ^ b) = log2 (a) + b. Esto hace que sea más fácil escribir su función log2: no es necesario que sea precisa para cada número positivo, solo en el intervalo [1, 2). De manera similar, para calcular 2 ^ x, puedes multiplicar 2 ^ (parte entera de x) * 2 ^ (parte fraccional de x). La primera parte es trivial para almacenar en un número de punto flotante, para la segunda parte, solo necesita una función 2 ^ x durante el intervalo [0, 1).

La parte difícil es encontrar una buena aproximación de 2 ^ x y log2 (x). Un enfoque simple es usar la serie de Taylor .

Los poderes negativos no son un problema, solo son el inverso ( 1/x ) de la potencia positiva.

Los poderes de punto flotante son un poco más complicados; como sabes, una potencia fraccionaria es equivalente a una raíz (por ejemplo, x^(1/2) == sqrt(x) ) y también sabes que multiplicar los poderes con la misma base es equivalente a agregar sus exponentes.

Con todo lo anterior, puedes:

Descompón el exponente en una parte entera y una parte racional .
Calcule la potencia del entero con un bucle (puede optimizar su descomposición en factores y reutilizar cálculos parciales).
Calcule la raíz con cualquier algoritmo que le guste (cualquier aproximación iterativa como la bisección o el método de Newton podría funcionar).
Multiplique el resultado.
Si el exponente fue negativo, aplica el inverso.

Ejemplo:

2^(-3.5) = (2^3 * 2^(1/2)))^-1 = 1 / (2*2*2 * sqrt(2))

Muchos enfoques se dan en otras respuestas. Aquí hay algo que pensé que podría ser útil en caso de poderes integrales.

En el caso del poder entero x de n ^x , el enfoque directo tomaría x-1 multiplicaciones. Para optimizar esto, podemos usar programación dinámica y reutilizar un resultado de multiplicación anterior para evitar todas las multiplicaciones x. Por ejemplo, en 5 ⁹ , podemos, por ejemplo, hacer lotes de 3 , es decir, calcular 5 ³ una vez, obtener 125 y luego el cubo 125 usando la misma lógica, tomando solo 4 multiplicaciones en el proceso, en lugar de 8 multiplicaciones con la manera directa .

La pregunta es ¿cuál es el tamaño ideal del lote b para que el número de multiplicaciones sea mínimo? Así que vamos a escribir la ecuación para esto. Si f (x, b) es la función que representa el número de multiplicaciones que conlleva el cálculo de n ^x utilizando el método anterior, entonces

Explicación: Un producto de lote de números p tomará multiplicaciones p-1. Si dividimos x multiplicaciones en b lotes, habría (x / b) -1 multiplicaciones requeridas dentro de cada lote, y b-1 multiplicaciones requeridas para todos los lotes b.

Ahora podemos calcular la primera derivada de esta función con respecto a b y equipararla a 0 para obtener la b para el menor número de multiplicaciones.

Ahora vuelva a colocar este valor de b en la función f (x, b) para obtener el menor número de multiplicaciones:

Para todas las x positivas, este valor es menor que las multiplicaciones por el camino directo.

Normalmente, la implementación de la función pow(double, double) en las bibliotecas matemáticas se basa en la identidad:

pow(x,y) = pow(a, y * log_a(x))

Usando esta identidad, solo necesita saber cómo elevar un único número a a un exponente arbitrario, y cómo tomar una base de logaritmo a . Ha convertido efectivamente una función multivariable compleja en dos funciones de una sola variable y una multiplicación, que es bastante fácil de implementar. Los valores más comúnmente elegidos de a son e o 2 - e porque e^x y log_e(1+x) tienen algunas propiedades matemáticas muy agradables, y 2 porque tiene algunas propiedades agradables para la implementación en aritmética de coma flotante.

El truco de hacerlo de esta manera es que (si desea obtener una precisión total) necesita calcular el término log_a(x) (y su producto con y ) con una precisión mayor que la representación de coma flotante de y . Por ejemplo, si y son dobles, y desea obtener un resultado de alta precisión, deberá encontrar alguna forma de almacenar los resultados intermedios (y hacer aritmética) en un formato de mayor precisión. El formato Intel x87 es una opción común, al igual que los enteros de 64 bits (aunque si realmente desea una implementación de alta calidad, tendrá que hacer un par de cálculos enteros de 96 bits, que son un poco dolorosos en algunos idiomas). Es mucho más fácil lidiar con esto si implementa powf(float,float) , porque entonces puede usar double para cálculos intermedios. Recomendaría comenzar con eso si quiere usar este enfoque.

El algoritmo que describí no es la única forma posible de calcular pow . Es simplemente el más adecuado para entregar un resultado de alta velocidad que satisfaga un límite fijo de precisión a priori . Es menos adecuado en algunos otros contextos, y es mucho más difícil de implementar que el algoritmo de repetición de cuadrados [raíz] que otros han sugerido.

Si desea probar el algoritmo cuadrado repetido [raíz], comience por escribir una función de potencia entera sin signo que solo use la cuadratura repetida. Una vez que tenga una buena comprensión del algoritmo para ese caso reducido, le resultará bastante sencillo extenderlo para manejar exponentes fraccionarios.

Para las potencias enteras positivas, observe la exponenciación mediante cuadratura y la exponenciación de la cadena de adición .

Por definición:

a ^ b = exp (b ln (a))

donde exp(x) = 1 + x + x^2/2 + x^3/3! + x^4/4! + x^5/5! + ... exp(x) = 1 + x + x^2/2 + x^3/3! + x^4/4! + x^5/5! + ...

donde n! = 1 * 2 * ... * n n! = 1 * 2 * ... * n .

En la práctica, ¡podría almacenar una matriz de los primeros 10 valores de 1/n! , y luego aproximar

exp(x) = 1 + x + x^2/2 + x^3/3! + ... + x^10/10!

porque 10! es un número enorme, ¡así que 1/10! es muy pequeño (2.7557319224⋅10 ^ -7).

Puedes encontrar la función pow como esta:

static double pows (double p_nombre, double p_puissance) { double nombre = p_nombre; double i=0; for(i=0; i < (p_puissance-1);i++){ nombre = nombre * p_nombre; } return (nombre); }

Puede encontrar la función de piso de esta manera:

static double floors(double p_nomber) { double x = p_nomber; long partent = (long) x; if (x<0) { return (partent-1); } else { return (partent); } }

Atentamente

Un algoritmo mejor para calcular de manera eficiente las potencias enteras positivas repetidamente cuadra la base, mientras se hace un seguimiento del resto multiplicando. Aquí hay una solución de muestra en Python que debería ser relativamente fácil de entender y traducir a su idioma preferido:

def power(base, exponent): remaining_multiplicand = 1 result = base while exponent > 1: remainder = exponent % 2 if remainder > 0: remaining_multiplicand = remaining_multiplicand * result exponent = (exponent - remainder) / 2 result = result * result return result * remaining_multiplicand

Para hacer que maneje exponentes negativos, todo lo que tiene que hacer es calcular la versión positiva y dividir 1 por el resultado, por lo que debería ser una modificación simple del código anterior. Los exponentes fraccionarios son considerablemente más difíciles, ya que significa esencialmente calcular una enésima raíz de la base, donde n = 1/abs(exponent % 1) y multiplicar el resultado por el resultado del cálculo de la potencia de la parte entera:

power(base, exponent - (exponent % 1))

Puede calcular las raíces hasta un nivel deseado de precisión utilizando el método de Newton. Echa un vistazo al artículo de wikipedia sobre el algoritmo .

Usando tres funciones auto implementadas iPow(x, n) , Ln(x) y Exp(x) , puedo calcular fPow(x, a) , x y a ser dobles . Ninguna de las siguientes funciones usa funciones de biblioteca, solo iteración.

Alguna explicación sobre las funciones implementadas:

(1) iPow(x, n) : x es double , n es int . Esta es una iteración simple, ya que n es un número entero.

(2) Ln(x) : esta función usa la iteración de la serie Taylor. La serie utilizada en la iteración es Σ (from int i = 0 to n) {(1 / (2 * i + 1)) * ((x - 1) / (x + 1)) ^ (2 * n + 1)} . El símbolo ^ denota la función de potencia Pow(x, n) implementada en la primera función, que usa iteración simple.

(3) Exp(x) : Esta función, nuevamente, usa la iteración de la serie Taylor. La serie utilizada en la iteración es Σ (from int i = 0 to n) {x^i / i!} . Aquí, el ^ denota la función de potencia, pero no se calcula llamando a la función 1st Pow(x, n) ; en su lugar, se implementa dentro de la 3ª función, al mismo tiempo que el factorial, utilizando d *= x / i . Sentí que tenía que usar este truco , porque en esta función, la iteración toma más pasos relativos a las otras funciones y el factorial (i!) desborda la mayor parte del tiempo. Para asegurarse de que la iteración no se desborde, la función de potencia en esta parte se itera al mismo tiempo que el factorial. De esta manera, superé el desbordamiento.

(4) fPow(x, a) : x y a son ambos dobles . Esta función no hace más que llamar a las otras tres funciones implementadas anteriormente. La idea principal de esta función depende de algunos cálculos: fPow(x, a) = Exp(a * Ln(x)) . Y ahora, tengo todas las funciones iPow , Ln y Exp con iteración ya.

nb Utilicé un constant MAX_DELTA_DOUBLE para decidir en qué paso detener la iteración. Lo configuré en 1.0E-15 , lo que parece razonable para los dobles. Por lo tanto, la iteración se detiene si (delta < MAX_DELTA_DOUBLE) Si necesita más precisión, puede usar el long double y disminuir el valor constante para MAX_DELTA_DOUBLE , por MAX_DELTA_DOUBLE , para 1.0E-18 (1.0E-18 sería el mínimo).

Aquí está el código, que funciona para mí.

#define MAX_DELTA_DOUBLE 1.0E-15 #define EULERS_NUMBER 2.718281828459045 double MathAbs_Double (double x) { return ((x >= 0) ? x : -x); } int MathAbs_Int (int x) { return ((x >= 0) ? x : -x); } double MathPow_Double_Int(double x, int n) { double ret; if ((x == 1.0) || (n == 1)) { ret = x; } else if (n < 0) { ret = 1.0 / MathPow_Double_Int(x, -n); } else { ret = 1.0; while (n--) { ret *= x; } } return (ret); } double MathLn_Double(double x) { double ret = 0.0, d; if (x > 0) { int n = 0; do { int a = 2 * n + 1; d = (1.0 / a) * MathPow_Double_Int((x - 1) / (x + 1), a); ret += d; n++; } while (MathAbs_Double(d) > MAX_DELTA_DOUBLE); } else { printf("/nerror: x < 0 in ln(x)/n"); exit(-1); } return (ret * 2); } double MathExp_Double(double x) { double ret; if (x == 1.0) { ret = EULERS_NUMBER; } else if (x < 0) { ret = 1.0 / MathExp_Double(-x); } else { int n = 2; double d; ret = 1.0 + x; do { d = x; for (int i = 2; i <= n; i++) { d *= x / i; } ret += d; n++; } while (d > MAX_DELTA_DOUBLE); } return (ret); } double MathPow_Double_Double(double x, double a) { double ret; if ((x == 1.0) || (a == 1.0)) { ret = x; } else if (a < 0) { ret = 1.0 / MathPow_Double_Double(x, -a); } else { ret = MathExp_Double(a * MathLn_Double(x)); } return (ret); }

Las funciones de Wolfram ofrecen una amplia variedad de fórmulas para calcular los poderes. Algunos de ellos serían muy sencillos de implementar.