rogr: Your problem is that [tt]SINGLE[/tt] and [tt]DOUBLE[/tt] do not use base 10 to store the number, but rather base 2. There are some decimal values -- fairly simple ones, even -- that cannot be stored in base 2. An example of how this could be impossible:
The number 1/3 could be stored as 0.1 is base 3 (2/3 would be 0.2 and 3/3 would be 1.0). However, in base 10, you would need an infinite number of digits to store the 1/3 (0.33333333....).
An example of a decimal number that can't be stored by a [tt]SINGLE[/tt] or even a [tt]DOUBLE[/tt] is 0.99. It seems like a fairly simple number, but the closest a [tt]SINGLE[/tt] can get is roughly 0.9899999.
As for your code, your method for finding cubic roots most likely will not work, because if the roots are not exact multiples of 1/100, the expression will never be exactly 0. What you instead want to check for is the expression crossing the boundary. You should also check its derivative, especially if it has a zero that immediately bounces back (the way that y=x^2 has a zero at x=0 but doesn't cross the X axis). Basically, if the derivative changes signs, then you should do a binary search within the range to find the local minimum or maximum. Then, since even doing a search this tightly will likely not yield exactly 0, you need to check the result against an epsilon -- an error threshold -- and report it as a probable zero if it is less than that threshold.
The following code implements this:
[tt]
'f(x) = (x - a)(x - b)(x - c)
'f(x) = (x + 5.7293)(x - 3)(x - 3)
a# = -5.7293 'these are the roots
b# = 3
c# = 3
'f(x) = x^3 - (a + b + c)x^2 + (ab + ac + bc)x - abc
'f(x) = Ax^3 + Bx^2 + Cx + D
cubicA# = 1#
cubicB# = -(a# + b# + c#)
cubicC# = a# * b# + a# * c# + b# * c#
cubicD# = -a# * b# * c#
DIM root#(3), rootDouble%(3)
numRoots% = 0
'Derivative f'(x) = 3Ax^2 + 2Bx + C
quadraticA# = 3# * cubicA#
quadraticB# = 2# * cubicB#
quadraticC# = cubicC#
epsilon# = .00000001# 'error threshold
x# = -100#
thisCubic# = cubicA# * x# ^ 3# + cubicB# * x# ^ 2# + cubicC# * x# + cubicD#
thisQuadratic# = quadraticA# * x# ^ 2# + quadraticB# * x# + quadraticC#
FOR i% = -9999 TO 10000
lastX# = x#
x# = i% / 100# 'this way, we preserve accuracy
lastCubic# = thisCubic#
lastQuadratic# = thisQuadratic#
thisCubic# = cubicA# * x# ^ 3# + cubicB# * x# ^ 2# + cubicC# * x# + cubicD#
thisQuadratic# = quadraticA# * x# ^ 2# + quadraticB# * x# + quadraticC#
IF thisCubic# = 0# THEN 'very unlikely unless the roots are multiples
numRoots% = numRoots% + 1 'of 1/2 or 1/5 , but could happen
root#(numRoots%) = x#
IF numRoots% = 3 THEN EXIT FOR
ELSEIF SGN(thisCubic#) * SGN(lastCubic#) < 0 THEN
lowBoundary# = lastX#
highBoundary# = x#
FOR j% = 1 TO 25 'higher is more precise, but unnecessary
testX# = lowBoundary# / 2# + highBoundary# / 2# '<-- slower, more accurate
testCubic# = cubicA# * testX# ^ 3# + cubicB# * testX# ^ 2# + cubicC# * testX# + cubicD#
IF SGN(testCubic#) = SGN(lastCubic#) THEN
lowBoundary# = testX#
ELSE
highBoundary# = testX#
END IF
NEXT j%
numRoots% = numRoots% + 1
root#(numRoots%) = lowBoundary# / 2# + highBoundary# / 2#
IF numRoots% = 3 THEN EXIT FOR
ELSEIF SGN(thisQuadratic#) * SGN(lastQuadratic#) < 0 THEN
IF lastCubic# = 0 THEN
rootDouble%(numRoots%) = -1
ELSE
testX# = lastX# / 2# + x# / 2
testCubic# = cubicA# * testX# ^ 3# + cubicB# * testX# ^ 2# + cubicC# * testX# + cubicD#
IF ABS(testCubic#) < ABS(thisCubic#) THEN 'it is closer to zero in the middle
lowBoundary# = lastX#
highBoundary# = x#
FOR j% = 1 TO 25 'higher is more precise, but unnecessary
testX# = lowBoundary# / 2# + highBoundary# / 2# '<-- slower, more accurate
testQuadratic# = quadraticA# * testX# ^ 2# + quadraticB# * testX# + quadraticC#
IF SGN(testQuadratic#) = SGN(lastQuadratic#) THEN
lowBoundary# = testX#
ELSE
highBoundary# = testX#
END IF
NEXT j%
testX# = lowBoundary# / 2# + highBoundary# / 2#
testCubic# = cubicA# * testX# ^ 3# + cubicB# * testX# ^ 2# + cubicC# * testX# + cubicD#
IF ABS(testCubic#) < epsilon# THEN
numRoots% = numRoots% + 1
root#(numRoots%) = testX#
rootDouble%(numRoots%) = -1
IF numRoots% = 3 THEN EXIT FOR
END IF
END IF
END IF
END IF
NEXT i%
PRINT "Roots found:"; numRoots%
FOR i% = 1 TO numRoots%
IF rootDouble%(i%) THEN
PRINT "Double root at: ";
ELSE
PRINT "Root at: ";
END IF
PRINT root#(i%)
NEXT i%
[/tt]