Fast sine approxymation

[angros47]

I found this algorithm here (http://www.musicdsp.org/archive.php?classid=0#261) and I ported to FreeBasic:

Code: Select all

function fastsin(phase as unsigned integer) as single
    const frf3 = -1.0 / 6.0
    const frf5 = 1.0 / 120.0
    const frf7 = -1.0 / 5040.0
    const frf9 = 1.0 / 362880.0
    const f0pi5 = 1.570796327
    dim as single x, x2, absin
    dim as unsigned integer tmp = &H3f800000 or (phase shr 7)
    if (phase and &H40000000) then tmp = tmp xor &H007fffff
    x = (*(cast(single ptr,@tmp)) - 1.0) * f0pi5
    x2 = x * x
    absin = ((((frf9 * x2 + frf7) * x2 + frf5) * x2 + frf3) * x2 + 1.0) * x
    return iif((phase and &H80000000),-absin , absin)
end function


screenres 800,600,32


for i as integer=0 to 800
	line -(i,fastsin(i*10000000)*100 + 300),rgb(255,255,255)

next

sleep

Phase is in 0 -> (2^32)-1 range and maps to 0 -> ~2PI


Try this:

Code: Select all

dim y as single=0
var x=timer

for i as integer=0 to 80000000
	y=y+sin(i)
next

?timer-x

vs this:

Code: Select all

dim y as single=0
var x=timer

for i as integer=0 to 80000000
	y=y+fastsin(i)
next

?timer-x

[srvaldez]

something is wrong, try the following test

Code: Select all

for i as integer=0 to 10
   ?sin(i)
next

?"----------------"

for i as integer=0 to 10
   ?fastsin(i)
next

sleep

fastsin just returns zeros

[fxm]

.....
Phase is in 0 -> (2^32)-1 range and maps to 0 -> ~2PI

[D.J.Peters]

Code: Select all

#include "crt/math.bi"
const as single  PI        = ATN(1.0)*4.0
const as single  DEG2RAD   = PI/180.0
const as single  PI2       = PI*2.0
const as longint PHASE_MAX = (2^32)-1
const as single  RAD2PHASE = PHASE_MAX / PI2

function fastsin(rad as single) as single
    const frf3 = -1.0 / 6.0
    const frf5 =  1.0 / 120.0
    const frf7 = -1.0 / 5040.0
    const frf9 =  1.0 / 362880.0
    const f0pi5 = 1.570796327
    dim as unsigned integer phase = RAD2PHASE*rad
    dim as single x, x2, absin
    dim as uinteger tmp = &H3f800000 or (phase shr 7)
    if (phase and &H40000000) then tmp = tmp xor &H007fffff
    x = (*(cast(single ptr,@tmp)) - 1.0) * f0pi5
    x2 = x * x
    absin = ((((frf9 * x2 + frf7) * x2 + frf5) * x2 + frf3) * x2 + 1.0) * x
    return iif((phase and &H80000000),-absin , absin)
end function

function nakedsin naked(byref rad as single) as single
  asm mov eax,[esp+4]
  asm fld dword ptr [eax]
  asm fsin
  asm ret 4
end function

dim as single rad 
dim as double t4,t3,t2,t1=timer

const n=10000000

print "please wait ..."
t1=timer
for i as integer=1 to n
  rad = (i mod 360)*DEG2RAD
  dim as single r = sin(rad)
next
t1=timer-t1
print "sin()      = " & t1

t2=timer
for i as integer=1 to n
  rad = (i mod 360)*DEG2RAD
  dim as single r = sinf(rad)
next
t2=timer-t2
print "sinf()     = " & t2

t3=timer
for i as integer=1 to n
  rad = (i mod 360)*DEG2RAD
  dim as single r = nakedsin(rad)
next
t3=timer-t3
print "nakedsin() = " & t3

t4=timer
for i as integer=1 to n
  rad = (i mod 360)*DEG2RAD
  dim as single r = fastsin(rad)
next
t4=timer-t4
print "fastsin()  = " & t4

sleep

[fxm]

Very small detail, I think rather that:
2^31 -> PI
2^30 -> PI / 2
2^29 -> PI / 4
...................
...................
2^0 -> PI / 2^31

therefore:
RAD2PHASE = (2^31) / PI
or
RAD2PHASE = (2^32) / PI2

[sean_vn]

Maybe this one:
Universal SIMD-Mathlibrary
http://webuser.hs-furtwangen.de/~dersch/