Comment puis-je écrire une fonction de puissance moi-même?

Les pouvoirs négatifs ne sont pas un problème, ils ne sont que l’inverse ( 1/x ) du pouvoir positif.

Les puissances en virgule flottante sont un peu plus compliquées. Comme vous le savez, une puissance fractionnaire est équivalente à une racine (par exemple, x^(1/2) == sqrt(x) ) et vous savez également que la multiplication des puissances avec la même base équivaut à append leurs exposants.

Avec tout ce qui précède, vous pouvez:

• Décomposez l’exposant dans une partie entière et une partie rationnelle .
• Calculez la puissance entière avec une boucle (vous pouvez l’optimiser en décomposant les facteurs et en réutilisant les calculs partiels).
• Calculez la racine avec n’importe quel algorithme (toute approximation itérative comme la bissection ou la méthode de Newton pourrait fonctionner).
• Multipliez le résultat.
• Si l’exposant était négatif, appliquer l’inverse.

Exemple:

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

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

Edit: oui, cette définition fournit vraiment quelque chose d’utile. Par exemple, sur un x86, cela se traduit presque directement en FYL2X (Y * Log ₂ (X)) et 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) 

Le code se termine un peu plus longtemps que prévu, principalement parce que F2XM1 ne fonctionne qu’avec des nombres compris entre -1.0 et 1.1.0. Le fld st(0)/frndint/fsubr st(1),st piece soustrait la partie entière, nous ne disposons donc que de la fraction. Nous y appliquons F2XM1 , ajoutons le 1, puis utilisons FSCALE pour gérer la partie entière de l’exponentiation.

Généralement, l’implémentation de la fonction pow(double, double) dans les bibliothèques mathématiques est basée sur l’identité:

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

En utilisant cette identité, il vous suffit de savoir comment élever un seul nombre a à un exposant arbitraire et comment utiliser une base logarithmique a . Vous avez effectivement transformé une fonction multi-variable complexe en deux fonctions d’une seule variable, et une multiplication, ce qui est assez simple à mettre en œuvre. Les valeurs les plus communément choisies de a sont e ou 2 – e car e^x et log_e(1+x) ont de très bonnes propriétés mathématiques, et 2 parce qu’elles ont de bonnes propriétés pour l’implémentation en arithmétique à virgule flottante.

Le fait de le faire de cette façon est que (si vous voulez obtenir une précision totale), vous devez calculer le terme log_a(x) (et son produit avec y ) avec une précision supérieure à la représentation à virgule flottante de x et y . Par exemple, si x et y sont des doublons et que vous souhaitez obtenir un résultat de haute précision, vous devrez trouver un moyen de stocker les résultats intermédiaires (et les calculs arithmétiques) dans un format plus précis. Le format Intel x87 est un choix courant, de même que les entiers 64 bits (même si vous voulez vraiment une implémentation de qualité supérieure, vous devrez faire quelques calculs sur des entiers 96 bits, ce qui est un peu pénible dans certains cas). langues). Il est beaucoup plus facile de gérer cela si vous implémentez powf(float,float) , car vous pouvez simplement utiliser le double pour les calculs intermédiaires. Je recommande de commencer par cela si vous souhaitez utiliser cette approche.

L’algorithme que j’ai décrit n’est pas le seul moyen possible de calculer le pow . Il est simplement le plus approprié pour délivrer un résultat à grande vitesse qui satisfait à une limite de précision a priori fixe. Il est moins adapté dans certains autres contextes et est certainement beaucoup plus difficile à mettre en œuvre que l’algorithme à carrés répétés [root] que d’autres ont suggéré.

Si vous voulez essayer l’algorithme carré répété [root], commencez par écrire une fonction de puissance entière non signée qui n’utilise que des carrés répétés. Une fois que vous avez bien compris l’algorithme pour ce cas réduit, vous trouverez qu’il est assez simple de l’étendre pour gérer les exposants fractionnaires.

Il existe deux cas distincts à traiter: les exposants entiers et les exposants fractionnaires.

Pour les exposants entiers, vous pouvez utiliser une exponentiation en faisant un carré.

 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) 

Le deuxième "elif" est ce qui le distingue de la fonction naïve. Il permet à la fonction de faire des appels récursifs O (log n) au lieu de O (n).

Pour les exposants fractionnaires, vous pouvez utiliser l'identité a ^ b = C ^ (b * log_C (a)). Il est pratique de prendre C = 2, donc a ^ b = 2 ^ (b * log2 (a)). Cela réduit le problème d'écriture des fonctions pour 2 ^ x et log2 (x).

La raison pour laquelle il est pratique de prendre C = 2 est que les nombres à virgule flottante sont stockés en virgule flottante en base 2. log2 (a * 2 ^ b) = log2 (a) + b. Cela facilite l'écriture de votre fonction log2: vous n'avez pas besoin de l'exactitude pour chaque nombre positif, juste sur l'intervalle [1, 2). De même, pour calculer 2 ^ x, vous pouvez multiplier 2 ^ (partie entière de x) * 2 ^ (partie fractionnaire de x). La première partie est sortingviale à stocker dans un nombre à virgule flottante, pour la deuxième partie, vous avez juste besoin d'une fonction 2 ^ x sur l'intervalle [0, 1).

La partie difficile consiste à trouver une bonne approximation de 2 ^ x et log2 (x). Une approche simple consiste à utiliser les séries de Taylor .

Par définition:

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

où 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! + ...

où n! = 1 * 2 * ... * n n! = 1 * 2 * ... * n .

En pratique, vous pouvez stocker un tableau des 10 premières valeurs de 1/n! , puis approximativement

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

parce que 10! est un nombre énorme, donc 1/10! est très petit (2.7557319224⋅10 ^ -7).

Les fonctions Wolfram offrent une grande variété de formules pour calculer les puissances. Certains d’entre eux seraient très simples à mettre en œuvre.

Pour les puissances entières positives, examinez l’ exponentiation par la quadrature et l’ exponentiation par chaîne d’addition .

En utilisant trois fonctions auto-implémentées iPow(x, n) , Ln(x) et Exp(x) , je suis capable de calculer fPow(x, a) , x et un être double . Aucune des fonctions ci-dessous n’utilise des fonctions de bibliothèque, mais seulement une itération.

Quelques explications sur les fonctions implémentées:

(1) iPow(x, n) : x est double , n est int . Ceci est une itération simple, car n est un entier.

(2) Ln(x) : cette fonction utilise l’itération Taylor Series. La série utilisée dans l’itération est Σ (from int i = 0 to n) {(1 / (2 * i + 1)) * ((x - 1) / (x + 1)) ^ (2 * n + 1)} Le symbole ^ désigne la fonction de puissance Pow(x, n) implémentée dans la 1ère fonction, qui utilise une itération simple.

(3) Exp(x) : cette fonction utilise à nouveau l’itération Taylor Series. La série utilisée dans l’itération est Σ (from int i = 0 to n) {x^i / i!} . Ici, le ^ désigne la fonction de puissance, mais il n’est pas calculé en appelant la 1ère fonction Pow(x, n) ; à la place, il est implémenté dans la 3ème fonction, en même temps que la factorielle, en utilisant d *= x / i . J’ai senti que je devais utiliser cette astuce , car dans cette fonction, l’itération prend plus d’étapes par rapport aux autres fonctions et le factoriel (i!) Déborde la plupart du temps. Afin de s’assurer que l’itération ne déborde pas, la fonction de puissance dans cette partie est itérée en même temps que la factorielle. De cette façon, j’ai surmonté le débordement.

(4) fPow(x, a) : x et a sont les deux doubles . Cette fonction ne fait rien mais appelle simplement les trois autres fonctions implémentées ci-dessus. L’idée principale de cette fonction dépend de certains calculs: fPow(x, a) = Exp(a * Ln(x)) . Et maintenant, j’ai déjà toutes les fonctions iPow , Ln et Exp avec itération.

nb J’ai utilisé une constant MAX_DELTA_DOUBLE pour décider à quelle étape arrêter l’itération. Je l’ai mis à 1.0E-15 , ce qui semble raisonnable pour les doubles. Donc, l’itération s’arrête si (delta < MAX_DELTA_DOUBLE) Si vous avez besoin de plus de précision, vous pouvez utiliser long double et diminuer la valeur constante pour MAX_DELTA_DOUBLE , à 1.0E-18 par exemple (1.0E-18 serait le minimum).

Voici le code qui fonctionne pour moi.

 #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); } 

Un meilleur algorithme pour calculer efficacement les puissances entières positives est la base carrée de façon répétée, tout en gardant une trace des multiplicands de rest supplémentaires. Voici un exemple de solution en Python qui devrait être relativement facile à comprendre et à traduire dans votre langue préférée:

 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 

Pour le faire traiter des exposants négatifs, il suffit de calculer la version positive et de diviser 1 par le résultat, cela devrait donc être une simple modification du code ci-dessus. Les exposants fractionnaires sont considérablement plus difficiles, car cela signifie essentiellement calculer une n-racine de la base, où n = 1/abs(exponent % 1) et multiplier le résultat par le résultat du calcul de la puissance de la partie entière:

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

Vous pouvez calculer les racines à un niveau de précision souhaité en utilisant la méthode de Newton. Consultez l’ article de Wikipedia sur l’algorithme .

C’est un exercice intéressant. Voici quelques suggestions que vous devriez essayer dans cet ordre:

1. Utilisez une boucle.
2. Utilisez la récursivité (pas mieux, mais néanmoins intéressante)
3. Optimisez considérablement votre récursion en utilisant des techniques de division et de conquête
4. Utiliser des logarithmes

Vous pouvez trouver la fonction de pow comme ceci:

 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); } 

Vous pouvez trouver la fonction de plancher comme ceci:

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

Meilleures salutations

J’utilise l’arithmétique longue en virgule fixe et mon pow est basé sur log2 / exp2. Les nombres consistent en:

• int sig = { -1; +1 } int sig = { -1; +1 } signum
• DWORD a[A+B]
• A est le nombre de DWORD pour la partie entière du nombre
• B est le nombre de DWORD pour la partie fractionnaire

Ma solution simplifiée est la suivante:

 //--------------------------------------------------------------------------- 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; } //--------------------------------------------------------------------------- 

où:

 _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 (<,>,==) 

Tout ce dont vous avez besoin pour comprendre ce code est que les nombres sous forme binary consistent en une sum de puissances de 2, quand vous avez besoin de calculer 2 ^ num alors il peut être réécrit comme ceci

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

où n sont des bits fractionnaires et m sont des bits entiers. la multiplication / division par 2 en forme binary est simple, donc si vous mettez tout cela ensemble, vous obtenez un code pour exp2 similaire à mon. log2 est basé sur la recherche binaru ... en changeant les bits de résultat de MSB à LSB jusqu'à ce qu'il corresponde à la valeur recherchée (algorithme très similaire à celui du calcul rapide de sqrt). J'espère que cela aide à clarifier les choses ...

Un grand nombre d’approches sont données dans d’autres réponses. Voici quelque chose que je pensais être utile dans le cas des pouvoirs intégraux.

Dans le cas de la puissance entière x de n ^x , l’approche simple prendrait des multiplications de x-1. Afin d’optimiser cela, nous pouvons utiliser la programmation dynamic et réutiliser un résultat de multiplication antérieur pour éviter toutes les multiplications x. Par exemple, dans 5 ⁹ , on peut, par exemple, faire des lots de 3 , c’est-à-dire calculer 5 ³ , obtenir 125, puis le cube 125 en utilisant la même logique, en ne prenant que 4 multiplications au lieu de 8 .

La question est de savoir quelle est la taille idéale du lot b pour que le nombre de multiplications soit minimum. Alors, écrivons l’équation pour cela. Si f (x, b) est la fonction représentant le nombre de multiplications impliquées dans le calcul de n ^{x en} utilisant la méthode ci-dessus, alors

Remets maintenant cette valeur de b dans la fonction f (x, b) pour obtenir le plus petit nombre de multiplications:

Pour tous les x positifs, cette valeur est inférieure aux multiplications de la manière la plus simple.