ip(x) and int(x) bug

TB 6.007
Win 7 pro
AMD FX8350
ASUS 990FX

I have a program in which I need to take the cube root of a series of positive integers and then find the integer portion of the result. So, the cube root of 3 is 1.4422.... and the value I want is the integer 1. However, using both ip(x) and int(x) I get anomalous outcomes for perfect cubes 8,27 and 64. In each case the value returned using both ip(x) and int(x) is 1 less than the correct value, e.g. the cube root of 8 is the integer 2 and yet when the value returned by the cube root equation is fed into these I get the result 1.

Below is sample code that invariably reproduces this anomaly on my computer. You will note that at 125 it produces the correct answer, but produces the wrong answer again at 343. I have not tried it at higher values.

FOR i=1 to 343
LET x=i^(1/3)
PRINT "i=";i, "Cube root of i=";x,"ip(cr of i)=";ip(x),"int(cr of i)=";int(x)
NEXT i
END

Comments

ip(x) and int(x) "bug"

This is not a bug in True Basic. Pretty much any other programming language would give the same result, including assembler.

The phenomenon here is the result of finite-precision computation. Suppose you could store only 4 decimal digits. Then the stored value of "1/3" would be .3333, which is obviously not equal to 1/3. Computer storage of numbers has the same limitation, just a lot more digits. So 64^(1/3) yields 7.9999... and the integer part of that is 7.

To allow for the way the FPU works, use a "very close to" function such as ROUND, or an exponent just a bit larger than 1/3. Both of the following codes give the result you want. (Incidentally, because True Basic is an interpreter rather than an optimizing compiler, storing 1/3 in inv3 saves a divide every iteration.)

LET inv3 = 1/3
FOR k = 1 to 30
LET k3 = k*k*k
LET cb = ROUND(k3^inv3,10) ! round to billionths
PRINT " cube =";k3, "Cube root =";cb, "int =";k-int(cb)
NEXT k

or

LET inv3 = 1/3
LET inv3 = inv3 + EPS(inv3) ! just a bit bigger
FOR k = 1 to 30
LET k3 = k*k*k
LET cb = k3^inv3
PRINT " cube =";k3, "Cube root =";cb, "int =";k-int(cb)
NEXT k

ought to be a better way

I'm not sure what kind of algorithm they are using to do this... I certainly understand things like the accumulation of rounding errors and so on. Fractions, even tiny ones, do add up. But in the case of something like ip(x), the way a human does it, we just drop off the decimals. No build up of unwanted decimals there!

Which makes me think it should be fairly easy to program a version that uses some of TrueBasic's string handling functions to make a number a string, find the decimal point, eliminate it and everything after it and then turn it back into a number.

Seems like it oughta work....

RE: ought to be a better way

It's not a bug (as other people have pointed out)

this example will show better what's happening:


print 8^(1/3) ! This rounds the value for display
print using "##.##############":8^(1/3) ! Now you can see the actual value computed
print ip(8^(1/3)) ! ip() is doing exactly what it's supposed to do. It works on the actual value
! not the displayed value.
end

Tom L

ip(x) and int(x) bug

Clearly, for LET x=i^(1/3), an exponentiation algorithm is used. The main problem is that you're supplying an exponent, "1/3", which cannot be represented exactly, so is a little smaller than the fraction 1/3. Hence the result of the exponentiation is a little smaller than the mathematical cube root and INT(x) is operating as defined by truncating this to an integer one smaller than what you think it should be. If you want the closest integer when you're already close, use ROUND(x,n). That's what STR$(x) will do, with n=8 I think.

Your basic numerical methods class covers many of the ways that finite-precision computation differs from the mathematical operations it approximates; this is only a simple example. What the computer does can differ from what you think it does.

There are ways to get a cube root algorithm whose results are always correct to the last bit, but exponentiation by a power not equal to 1/3 is not one of them. Solving x^3 - i = 0 numerically can do it, for instance, as "^3" implemented as x*x*x is exact to the last bit. More convenient is a hardware (OK, firmware) cube root, unfortunately much rarer than a hardware square root. The BASIC standard specifies a SQR function but not a CBR, a pity. Fortran doesn't have one either.

ip(x) and int(x) bug and digital crumbs

Hi,

Hmmm. That's the first bug I've ever encountered in True Basic. It's not just "version 6", since the underlying core code for any and all versions is exactly the same. I have Silver 5.33 and Gold 5.5b19, and both give this same anomaly. I'm sure the old DOS version would also behave the same, but I can't run that one now.

Since there is no one who knows how to change/modify the core code, it's not ever going to get fixed. This bug is disturbing.

Anyway, to get an idea of the problem, add the following lines and note that some digital "crumbs" are left over and somehow this may be causing the problem.

print
print,"1-cube root of 1=";1-(1^(1/3))
print,"2-cube root of 8=";2-(8^(1/3))
print,"3-cube root of 27=";3-(27^(1/3))
print,"4-cube root of 64=";4-(64^(1/3))
print,"5-cube root of 125=";5-(125^(1/3))
print,"6-cube root of 216=";6-(216^(1/3))
print,"7-cube root of 343=";7-(343^(1/3))
print,"8-cube root of 512=";8-(512^(1/3))
print,"9-cube root of 729=";9-(729^(1/3))
print,"10-cube root of 1000=";10-(1000^(1/3))
print,"11-cube root of 1331=";11-(1331^(1/3))
print,"12-cube root of 1728=";12-(1728^(1/3))
print,"13-cube root of 2197=";13-(2197^(1/3)) ! etc

print,"3-exp(log(27)/3)=";3-exp(log(27)/3) ! still has crumbs...

! BUT: use this instead:

LET x=2^(log2(i)/3) ! or: LET x=10^(log10(i)/3), or: LET x=exp(log(i)/3) ; all of which which seem to work correctly.

fascinating. (Never had this problem with a slide rule back in the day!)

Regards,
Mike C.

Other methods seem buggy too...

Hi Mike,

I only go back to just after the end of the slide rule era, so you may have a few years on me. I started college about the time the TI-58C came out. I started using True Basic when they came out with the Atari ST version. I bought the "Gold" version for the Atari back then. At the time, I believe it came with "lifetime" upgrades. A few years ago when I decided to do some True Basic programming on a Windows PC. I tried to get them to give me a break on the Windows version, since the Atari ST is "legacy"... no dice. The sad thing is that I haven't been able to access my old Atari disks and get all of the add-on libraries off of it. I could've gotten the Bronze version for the PC and just recompiled the old libs on the PC...

I tried some of your alternates for finding cube roots and then found an anomaly on one, so I created the program below to test them all. I also realized that I could try just subtracting the floating point part to create my own "int" function, so I tried that with each of the methods. Nothing seems to work perfectly. After the fact I realize I could have tried your log method with a variety of different logs by just using a loop rather than hardcoding the logs. Oh, well. Just goes to show I'm not a real programmer.

Cheers,
Lee

FOR x=1 to 10
let i=x^3
If x<>ip(i^(1/3)) then print "anomaly for ip(i^(1/3)) when i=";i
If x<>int(i^(1/3)) then print "anomaly for int(i^(1/3)) when i=";i
If x<>(i^(1/3))-fp(i^(1/3)) then print "anomaly for (i^(1/3))-fp(i^(1/3)) when i=";i
If x<>ip(exp(log(i)/3)) then print "anomaly for ip(exp(log(i)/3)) when i=";i
If x<>int(exp(log(i)/3)) then print "anomaly for int(exp(log(i)/3)) when i=";i
If x<>(exp(log(i)/3))-fp(exp(log(i)/3)) then print "anomaly for (exp(log(i)/3))-fp(exp(log(i)/3)) when i=";i
If x<>ip(10^(log10(i)/3)) then print "anomaly for ip(10^(log10(i)/3)) when i=";i
If x<>int(10^(log10(i)/3)) then print "anomaly for int(10^(log10(i)/3)) when i=";i
If x<>(10^(log10(i)/3))-fp(10^(log10(i)/3)) then print "anomaly for (10^(log10(i)/3))-fp(10^(log10(i)/3)) when=";i
If x<>ip(2^(log10(i)/3)) then print "anomaly for ip(2^(log2(i)/3)) when i=";i
If x<>int(2^(log10(i)/3)) then print "anomaly for int(2^(log2(i)/3)) when i=";i
If x<>(2^(log10(i)/3))-fp(2^(log10(i)/3)) then print "anomaly for (2^(log10(i)/3))-fp(2^(log10(i)/3)) when i=";i
next x
end