Power function with fractional exponents

[ntran]

Hi everyone,

I'm trying to make a program to compute a^n, where a and n are rational numbers. Ex. 5^(0.345)
Without using predefined functions like pow() or sqrt().

I have found this method:

- I read both, the base and the exponent (user input).
- Then convert the exponent from decimal to fraction:
- base^(exponent) = base^(num/den) = x
- Elevating both sides to (den) gives: x^(den) = base^(num)
- So I have to solve: x^(den) - base^(num) = 0; where x is unknown.
- I'm using Newton's method to find the root...

This is my problem: Most of the times, any of the two factors of the equation: x^(den) or base^(num), become huge numbers and go out of range, and I get inf and nan errors.

Can someone suggest a solution to this problem? or
Can someone suggest another method to achieve a power function?

Thanks.

#include <iostream>
using namespace std;

//  Computes a^n
//  Where a and n are rational numbers

main(){
	
	float a;   				// Base
	float n;     				// Exponent
	int num;	 			// Numerator of the exponent
	int den = 1; 				// Denominator of the exponent
	
	cout << "Base: ";
	cin >> a;
	cout << "Exponent: ";
	cin >> n;
	
	//----------------------------------------------------------------
	
	while (n != (int) n) {			// Convert n from decimal to fraction:
		n = n*10;			// Multiplying numerator and denominator x10 
		den = den*10;			// until both become whole integers
	}
	num = n;				// Now: a^n --> a^(num/den)

	//----------------------------------------------------------------

	int t_x = num;
	int t_y = den;
	int temp;
	
	while (t_y != 0){
		temp = t_x%t_y;
		t_x = t_y;
		t_y = temp;
	}
	
	num = num/t_x;
	den = den/t_x;				// Simplifying the (num/den) expression to lowest terms

	//----------------------------------------------------------------

						//  Solve x = a^(num/den)
						//  Rising both sides to the (den) power: x^(den) = a^(num)
						//  Passing all the terms to one side of the equation:
						//  x^(den) - a^(num) = 0
						//  Finding the root with Newton's method
	
	float x;				// Next x - Newton's method
	float x0 = 1;				// Current x - Newton's method, initial value set to 1
	float tol = 1;				// Tolerance - Newton's method
	float atonum = a;			// Variable for computing a^(num)
	float xtoden;				// Variable for computing x0^(den)
	
	for (int i=1; i<num; i++) atonum = atonum*a;
	cout << endl <<"  Solving:  x^" << den << " - " << atonum << " = 0" << endl << endl;
	
	while (tol > 0.001){							//
		xtoden = x0;							//
		for (int i=1; i<den; i++) xtoden = xtoden*x0;			//
		x = x0 - (xtoden-atonum) / ((den)*(xtoden/x0));			//
		tol = ((x-x0)/x0)*100;						//  Newton's Method Iterations
		if (tol < 0) tol = tol*(-1);					//
		cout << "x0 = " << x0 << endl;					//
		cout << "x = " << x << endl;					//
		cout << "tol = " << tol << endl << endl;			//
		x0 = x;								//
	}
	
	cout << "Result = " << x;   		//  Displaying the result
}

[funprogrammer]

I think the problem is related to atonum xtoden variables going out of range. Did you try using long double instead of float for these variables? I am not saying this will fix the issue, but its definitely worth the try

[ntran]

Yes funprogrammer, the problem is exactly what you point out.
Those variables are going out of range. I have already try using long double, but it won't help much.

Ex. If I try 5.16^(0.3333), the program will break, even with long double.

Thanks for your answer.

[funprogrammer]

@ntran: I stepped through your code, even though I dont understand the Newton's Method Iterations. The atonum variable in line 56 is what is causing the issue. I don't know if this method of a^num with the for loop is the correct way of doing it. I would check your algorithm again to see if you are following the correct one and also have you ported the correct algorithm to C++ code.