Correctly exponentiating a rational number by a real

Help
1

I am working on a project to implement a rational number type which overloads as many intrinsic operators and procedures as I can stand to write tests for. I want to solicit opinions on how best to implement operator(**) for raising a rational number to a real power. I can think of two natural ways to do it

  1. Convert the rational number to a real (divide numerator by denominator) and then exponentiate the result.
  2. Exponentiate the numerator and denominator of the rational number separately, and then divide the two results. In the case of negative powers, swap the numerator and denominator first and raise them to a positive power.

Due to the finite precision of floating point, these two are not equivalent. I can come up with test cases which make either one look superior.

For example, consider the expression rational(2, 3)**(-3.0) representing raising 2/3 to the -3 power, which “ought to” result in 27/8. The true-mathematics result is exactly representable in floating point, but method #1 will not give it due to the inexact intermediate step in computing a floating-point approximation to 2/3. Method #2 will give the exact result, since the intermediate steps will all be exactly representable.

As a second example, consider the expression rational(1, 2)**3.14. Here, the rational number has an exact floating-point representation, but the exponent does not. The answer to 17 decimal places, compared with the above two methods is

  • exact: 0.11343989441464511
  • #1 (divide first): 0.11343989441464510
  • #2 (powers first): 0.11343989441464508

which points toward method #1 as being the more appropriate method. It would seem that when the exponent is not exactly representable, one is better off doing one inexact exponentiation rather than two.

I checked what Julia does, since it has native rational numbers, and determined that it performs the division first (in the implementation, the rational is “promoted” to the type of the exponent). I’m curious what other languages with native rational numbers (or widely used libraries) elect to do. I’m also curious if any readers here have an intuitive sense for which behavior seems better.

2 Likes
2

Another question would be: is it possible to select automatically the best algorithm depending on the values?
knowing that numbers that can be represented exactly in IEEE 754 are:

  • integers (not too big),
  • rationals whose denominator is a power of 2 (with not too many decimals).
3

For your second example with a non-rational exponent, another approach is to perform the exponentiation in two steps: first for the integer part of the exponent using integer exponentiation and second for the remaining floating point part using floating point exponentiation.

A quick test shows that this gives: 0.113439894414645109366190922628… which is the closest double-precision floating-point number to the result when using more precision.

program test
  use iso_fortran_env, only: dp=>real64, int32
  implicit none

  real(dp) :: a, b, c

  a = 1.0_dp
  b = 2.0_dp

  c = a/b
  c = (c**int(3,int32))*(c**(0.14_dp))
  write(*,'(F70.50)') c

end program test
2 Likes
5

With floating-point computations involving irrational numbers especially, usually one can only achieve 2 of the following 3:

  • accuracy,
  • portability,
  • performance

If portability and performance are not as much of a concern (original post didn’t seem to indicate they were) as accuracy, then the simpleton approach can help with cases like the couple listed above:

   r = exp( a * log( p/q ) )

while it’s possible with this approach the portability won’t suffer as much either if one is using IEEE floating-point arithmetic only.

2 Likes
7

This is an old topic, but I wonder if one is splitting hairs here?

The program

program ratexp
implicit none
double precision x,y,z,p
integer j(2)
x = 1.0
y = 2.0
p = 3.14
z = (x/y)**p             ! 1st variant
j = transfer(z,j)
print '(1x2Z8)',j(2),j(1)
z = x**p/y**p            ! 2nd variant
j = transfer(z,j)
print '(1x2Z8)',j(2),j(1)
end program

seems to give the same results, whether the compiler generates x87 code or SSE2 code:

 3FBD0A65792512C5
 3FBD0A65792512C5

Even if one or two of the least significant bits differed, we should not regard that as bad.