Wer ist online?
Buchempfehlung
FreeBASIC-Chat
FreeBASIC-Nachrichten jetzt auch über Twitter erhalten. Follow us!fb:porticula NoPaste
t_cubert.bas
| Uploader: |  Volta |
| Datum/Zeit: | 30.09.2010 20:42:05 |
'PROGRAM test_cube_root
#Include "quad_1x.bi"
Using quad_precision
#define EPSILON_HALF 2.22044604925031308E-016
Declare Function cube_root(y As quad) As quad
Dim As quad x, lhs, rhs, diff
Dim As Double r
Dim As Integer i
Screen 0
Width 90,25
Color 15,1
Cls
Randomize 4711
For i = 1 To 20
x.hi=Rnd
r=Rnd
x.hi = (x.hi - 0.5) / r
x.lo = x.hi * EPSILON_HALF
lhs = cube_root(x)
If (x.hi >= 0.0) Then
rhs = Exp( Log(x) / 3.0)
Else
rhs = - Exp( Log(- x) / 3.0)
End If
diff = lhs - rhs
Print "lhs ="; lhs, " Diff. ="; diff
Next
Sleep
Function cube_root(y As quad) As quad 'RESULT(x)
' Calculates the cube root of y in quadruple-precision.
' Programmer : Alan Miller, CSIRO Mathematical Information Sciences
' Latest revision - 6 December 1996
Dim As quad f, xsq, x
Dim As Double one_third = 0.33333333333333330
' First approximation is the double-precision cube root.
' Only one iteration is necessary.
x.hi = Abs(y.hi)^ one_third
If Sgn(y.hi) < 0.0 Then x.hi= -x.hi
' Put f(x) = x**3 - y
' Then Newton's iterative method uses:
' new x = x - f(x) / f'(x)
' Halley's method uses:
' new x = x - f / (f' - 0.5 * f * f'' / f')
' In this case
' new x = x - x * f / (3.x^3 - f)
x.lo = 0.0
xsq = exactmul2(x.hi, x.hi)
f = x * xsq - y
x = x - x * f / (3.0 * x * xsq - f)
Return x
End Function' cube_root