2

I have a problem with what i guess is a rounding error with floating-points in OpenEdge ABL / Progress 4GL

display  truncate(log(4) / log(2) , 0) .

This returns 1.0 but should give me a 2.0

if i do this pseudo solution it gives me the right answer in most cases which hints to floating-points.

display  truncate(log(4) / log(2)  + 0.00000001, 0) .

What I am after is this 

find the largest x where 

p^x < n, p is prime, n and x is natural numbers.

=>

x = log(n) / log(p)

Any takes on this one?

Tom Bascom
  • 13,405
  • 2
  • 27
  • 33
Zesar
  • 566
  • 2
  • 14

4 Answers4

2

No numerical arithmetic system is exact. The natural logarithms of 4 and 2 cannot be represented exactly. Since the log function can only return a representable value, it returns an approximation of the exact mathematical result.

Sometimes this approximation will be slightly higher than the mathematical result. Sometimes it will be slightly lower. Therefore, you cannot generally expect that log(x*x) will be exactly twice log(x).

Ideally, a high-quality log implementation would return the representable value that is closest to the exact mathematical value. (This is called a “correctly rounded” result.) In that case, and if you are using binary floating-point (which is common), then log(4) would always be exactly twice log(2). Since this does not happen for you, it seems the log implementation you are using does not provide correctly rounded results.

However, for this problem, you also need log(8) to be exactly three times log(2), and so on for additional powers. Even if the log implementation did return correctly rounded results, this would not necessarily be true for all the values you need. For some y = x5, log(y) might not be exactly five times log(x), because rounding log(y) to the closest representable value might round down while rounding log(x) rounds up, just because of where the exact values happen to lie relative to the nearest representable values.

Therefore, you cannot rely on even a best-possible log implementation to tell you exactly how many powers of x divide some number y. You can get close, and then you can test the result by confirming or denying it with integer arithmetic. There are likely other approaches depending upon the needs specific to your situation.

Eric Postpischil
  • 195,579
  • 13
  • 168
  • 312
  • I'm a little confused by this. With a correctly rounded `log` and an `x` where `x*x` is exactly representable and no overflow or underflow shenanigans happening, surely `log(x*x) == 2*log(x)`? – tmyklebu Oct 02 '13 at 15:23
  • @tmyklebu: Yes, that would be true with binary floating-point. However, for this problem, we also need `log(x*x*x) == 3*log(x)`, and that does not hold in binary floating-point. I will update the answer. – Eric Postpischil Oct 02 '13 at 15:31
0

I think you want:

/* find the largest x where p^x < n, p is prime, n and x is natural numbers.
 */

define variable p as integer no-undo format ">,>>>,>>>,>>9".
define variable x as integer no-undo format ">>9".
define variable n as integer no-undo format ">,>>>,>>>,>>9".

define variable i as integer no-undo format "->>9".

define variable z as decimal no-undo format ">>9.9999999999".

update p n with side-labels.

/* approximate x
 */

z = log( n ) / log( p ).
display z.

x = integer( truncate( z, 0 )).  /* estimate x                */

/* is p^x < n ?
 */

if exp( p, x ) >= n then
  do while exp( p, x ) >= n:    /* was the estimate too high? */
    assign
      i = i - 1
      x = x - 1
    .
  end.
 else
  do while exp( p, x + 1 ) < n:  /* was the estimate too low? */
    assign
      i = i + 1
      x = x + 1
    .
  end.

display
  x skip
  exp( p, x ) label "p^x" format ">,>>>,>>>,>>9" skip
  i skip
  log( n ) skip
  log( p ) skip
  z skip
 with
  side-labels
.
Tom Bascom
  • 13,405
  • 2
  • 27
  • 33
  • this creates an error when the result is close to a integer solution, which will give you a false positive. – Zesar Oct 07 '13 at 08:40
  • This does not ever throw an error. It is not any more sensitive to false positives than any other solution being offered. – Tom Bascom Oct 07 '13 at 11:18
  • 1
    if you round the value at position 0 you will get the wrong value with eg. log(9)/log(8) => 1 and not 0 as is requested by the definition. when I said error I meant in the math, sorry for the unspecificness ;-) – Zesar Oct 07 '13 at 11:36
  • display log(8) format "9.9999999999" log(9) format "9.9999999999" truncate( round( log(9) / log(8), 0 ), 0 ) . 2.0794415417 2.1972245773 1.00 – Tom Bascom Oct 07 '13 at 14:25
  • @Mark: Perhaps you should re-read your own answer? You say "Absent any other preferable options, like a different rounding function...", well there is another rounding function. It is called ROUND(). – Tom Bascom Oct 07 '13 at 18:09
  • @Mark -- I do understand the comments about representing real numbers. Obviously simply adding some constant to a result is not going to correctly address those problems. Likewise, failing to understand that the Progress 4GL TRUNCATE() function is not the same as the ROUND() function isn't going to solve anything either. – Tom Bascom Oct 07 '13 at 20:11
  • Given that LOG() is inexact and that the error in the result could fall on either side of the nearest integer I've put together some code that tests p^x < n and adjusts x up or down as needed (as Eric suggests). I tossed in a few debugging results too so that you can more readily follow the calculations and see what is going on. – Tom Bascom Oct 07 '13 at 20:46
0

The root of the problem is that the log function, susceptible to floating point truncation error, is being used to address a question in the realm of natural numbers. First, I should point out that actually, in the example given, 1 really is the correct answer. We are looking for the largest x such that p^x < n; not p^x <= n. 2^1 < 4, but 2^2 is not. That said, we still have a problem, because when p^x = n for some x, log(n) divided by log(p) could probably just as well land slightly above the whole number rather than below, unless there is some systemic bias in the implementation of the log function. So in this case where there is some x for which p^x=n, we actually want to be sure to round down to the next lower whole value for x.

So even a solution like this will not correct this problem:

display  truncate(round(log(4) / log(2), 10) , 0) .

I see two ways to deal with this. One is similar to what you already tried, except that because we actually want to round down to the next lower natural number, we would subtract rather than add:

display  truncate(log(4) / log(2)  - 0.00000001, 0) .

This will work as long as n is less than 10^16, but a more tidy solution would be to settle the boundary conditions with actual integer math. Of course, this will fail too if you get to numbers that are higher than the maximum integer value. But if this is not a concern, you can just use your first solution get the approximate solution:

display  truncate(log(4) / log(2) , 0) .

And then test whether the result works in the equation p^x < n. If it isn't less than n, subtract one and try again.

On a side note, by the way, the definition of natural numbers does not include zero, so if the lowest possible value for x is 1, then the lowest possible value for p^x is p, so if n is less than or equal to p, there is no natural number solution.

Mark Bailey
  • 1,091
  • 10
  • 25
0

Most calculators can not calculate sqrt{2}*sqrt{2} either. The problem is that we usually do not have that many decimals.

Work around: Avoid TRUNCATE use ROUND like

 ROUND(log(4) / log(2), 0).

Round(a,b) rounds up the decimal a to closest number having b decimals.

AD - Stop Putin -
  • 295
  • 3
  • 13