CORDIC - Rosetta Code

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
Post Reply
srvaldez
Posts: 3545
Joined: Sep 25, 2005 21:54

CORDIC - Rosetta Code

Post by srvaldez »

the FreeBASIC translation of the Julia code is wrong https://rosettacode.org/wiki/CORDIC#FreeBASIC
in Julia the variable v looks like a vector, v = [1, 0] # start with 2-vector cosine and sine of zero
and in the call to the cordic function cosθ, sinθ = cordic(θr) cosθ and sinθ should have the respective values of cosθ and sinθ
an easier code to to translate is the C code https://rosettacode.org/wiki/CORDIC#C

Code: Select all

    'translation of the C code https://rosettacode.org/wiki/CORDIC#C

    Dim Shared angles(0 To ...) As Const Double = {0.78539816339745, 0.46364760900081,_
            0.24497866312686, 0.12435499454676, 0.06241880999596, 0.03123983343027,_
            0.01562372862048, 0.00781234106010, 0.00390623013197, 0.00195312251648,_
            0.00097656218956, 0.00048828121119, 0.00024414062015, 0.00012207031189,_
            0.00006103515617, 0.00003051757812, 0.00001525878906, 0.00000762939453,_
            0.00000381469727, 0.00000190734863, 0.00000095367432, 0.00000047683716,_
            0.00000023841858, 0.00000011920929, 0.00000005960464, 0.00000002980232,_
            0.00000001490116, 0.00000000745058}
            
    Dim Shared kvalues(0 To ...) As Const Double = {0.70710678118655, 0.63245553203368,_
            0.61357199107790, 0.60883391251775, 0.60764825625617, 0.60735177014130,_
            0.60727764409353, 0.60725911229889, 0.60725447933256, 0.60725332108988,_
            0.60725303152913, 0.60725295913894, 0.60725294104140, 0.60725293651701,_
            0.60725293538591, 0.60725293510314, 0.60725293503245, 0.60725293501477,_
            0.60725293501035, 0.60725293500925, 0.60725293500897, 0.60725293500890,_
            0.60725293500889, 0.60725293500888}
            
    #Define radians(degrees) Cdbl(((degrees) * M_PI) / 180.0)

    Const M_PI = 3.141592653589793

    Sub Cordic(Byval Alpha As Double, Byval n As Long, Byref c_cos As Double, Byref c_sin As Double)
        Dim As Long i, ix, sigma
        Dim As Double kn, x, y, datn, t
        Dim As Double theta = 0.0
        Dim As Double pow2 = 1.0
        Dim As Long newsgn = Iif((Int(Alpha / (2.0 * M_PI)) Mod 2) = 1, 1, -1)
        If (Alpha < ((-M_PI) / 2.0)) Orelse (Alpha > (M_PI / 2.0)) Then
            If Alpha < 0 Then
                Cordic(Alpha + M_PI, n, x, y)
            Else
                Cordic(Alpha - M_PI, n, x, y)
            End If
            c_cos = x * newsgn
            c_sin = y * newsgn
            Return
        End If
        ix = n - 1
        If ix > 23 Then
            ix = 23
        End If
        kn = kvalues(ix)
        x = 1
        y = 0
        For i = 0 To n - 1
            datn = angles(i)
            sigma = Iif(theta < Alpha, 1, -1)
            theta += sigma * datn
            t = x
            x -= (sigma * y) * pow2
            y += (sigma * t) * pow2
            pow2 /= 2.0
        Next
        c_cos = x * kn
        c_sin = y * kn
    End Sub

    Dim i As Long
    Dim th As Long
    Dim thr As Double
    Dim c_cos As Double
    Dim c_sin As Double
    Dim angle(0 To ...) As Double = {-9.0, 0.0, 1.5, 6.0}
    Dim f As String
    Print "  x       sin(x)     diff. sine     cos(x)    diff. cosine"
    f = "+##.#  +##.######## (+##.########) +##.######## (+##.########)"
    For th = -90 To 90 Step 15
        thr = radians(th)
        Cordic(thr, 24, c_cos, c_sin)
        Print Using f; th; c_sin; c_sin - Sin(thr); c_cos; c_cos - Cos(thr)
    Next
    Print
    Print "x(rads)   sin(x)     diff. sine     cos(x)    diff. cosine"
    f = "+#.#   +##.######## (+##.########) +##.######## (+##.########)"
    For i=0 To 3
        thr = angle(i)
        Cordic(thr, 24, c_cos, c_sin)
        Print Using f; thr; c_sin; c_sin - Sin(thr); c_cos; c_cos - Cos(thr)
    Next
output

Code: Select all

  x       sin(x)     diff. sine     cos(x)    diff. cosine
-90.0   -1.00000000 ( +0.00000000)  -0.00000007 ( -0.00000007)
-75.0   -0.96592585 ( -0.00000003)  +0.25881895 ( -0.00000009)
-60.0   -0.86602545 ( -0.00000005)  +0.49999992 ( -0.00000008)
-45.0   -0.70710684 ( -0.00000006)  +0.70710672 ( -0.00000006)
-30.0   -0.49999992 ( +0.00000008)  +0.86602545 ( +0.00000005)
-15.0   -0.25881895 ( +0.00000009)  +0.96592585 ( +0.00000003)
 +0.0   +0.00000007 ( +0.00000007)  +1.00000000 ( -0.00000000)
+15.0   +0.25881895 ( -0.00000009)  +0.96592585 ( +0.00000003)
+30.0   +0.49999992 ( -0.00000008)  +0.86602545 ( +0.00000005)
+45.0   +0.70710684 ( +0.00000006)  +0.70710672 ( -0.00000006)
+60.0   +0.86602545 ( +0.00000005)  +0.49999992 ( -0.00000008)
+75.0   +0.96592585 ( +0.00000003)  +0.25881895 ( -0.00000009)
+90.0   +1.00000000 ( -0.00000000)  -0.00000007 ( -0.00000007)

x(rads)   sin(x)     diff. sine     cos(x)    diff. cosine
-9.0    -0.41211842 ( +0.00000006)  -0.91113029 ( -0.00000003)
+0.0    +0.00000007 ( +0.00000007)  +1.00000000 ( -0.00000000)
+1.5    +0.99749499 ( +0.00000000)  +0.07073719 ( -0.00000002)
+6.0    -0.27941552 ( -0.00000002)  +0.96017028 ( -0.00000001)
Post Reply