There is a good way to compute rotation number of a circle homeomorphism (this was the way Poincaré thinked of it): you calculate the rotation number buy its continued fraction in a direct way.

You start from a point $x$ and $f(x)$: this gives you a decomposition of the circle into points that are on the right side of $x$ (in $]x,f(x)[$)  and points which are on its left side (in $]f(x),x[$). You look at $f^2(x)$ and you write $R$ if it is on the right side of $x$, $L$ otherwise. Iteranting $f$ you find a sequence of $R$'s and $L$'s. If you get $LLLLR$, for example, you record 4 (this is the number of $L$'s) and you approximate the rotation number of $f$ by $1/4$.

Renormalizing $f$, you iterate this process finding $\rho=[0,a_1,a_2,\ldots,a_k]$.

I won't be more precise here.

Every detail is very well explained in de Melo & van Strien's [*One-Dimensional Dynamics*][1], section I.1.

You can find a paper by Bruin ([*Numerical determination of the continued fraction expansion of the rotation number*][2]) in which he compares different methods on Arnold tongues.

EDIT[update]:
Recently, I wrote for myself some [sage][3] lines implementing the algorithm I described you. 
This is my second version, now working for rational numbers too : I was originally interested only in irrational rotation numbers (comments are welcome to improve it!).

    L=8 #length for cf-expansion, depending on your computer, 8 or 9 suggested for a try run
    A=100000 #maximum size of single element of the sequence
    
    def partfrac(x):
        return x-floor(x)

    ##### computing rational approximations given continued fraction expansion
    # input b=a continued fraction expansion
    # input l=L length of computed expansion
    def rational_approximation(b,l):
        p=[0,1]
        q=[1,b[1]]
        for i in range(1,l+1):
            p.append(b[i+1]*p[i]+p[i-1])
            q.append(b[i+1]*q[i]+q[i-1])
        return simplify(p[l+1]/q[l+1])


    #computing rotation number of a given circle map f
    def rotation(f):
        a=[0]
        orbit=[]
        orbit.append(partfrac(f(0)))
        if orbit[0]==0 :
            print 'map with a fixed point'
            return 0
        
        def shift(x):  #set f(0) as the origin + 1
            if partfrac(x)>orbit[0]:
                return partfrac(x)-1
            return partfrac(x)
            
        def first_return(p,pre_p,y):
            x=shift(f(y))
            while x<pre_p or x>p:
                x=shift(f(x))
            return x
    
        a.append(1)
        x=orbit[0]
        
        if shift(f(orbit[0]))==0:
            print 'map with periodic point of order 2'
            return 1/2
        
        if shift(f(orbit[0]))<0:
            while shift(f(x))<0:
                a[1]=a[1]+1
                x = shift(f(x))
                if a[1]>A:
                    print 'approximatively 0'
                    return 0
                if shift(f(x))==0:
                    print 'periodic point'
                    a[1]=a[1]+1
                    return 1/a[1]
            orbit.append(shift(x))
            z = shift(f(x))
            a.append(0)
            while z>0:
                y = z
                z = first_return(shift(orbit[0]),shift(orbit[1]),z)
                a[2]=a[2]+1
                if a[2]>A:
                    print 'approximatively rational'
                    return 1/a[1]
                if z==0:
                    print 'periodic point'
                    a[2]=a[2]+1
                    return rational_approximation(a,1)
            orbit.append(y)
        
        if shift(f(orbit[0]))>0:
            def shift(y):  #set f(0) as the origin
                if partfrac(y)>=orbit[0]:
                    return partfrac(y)-1
                    return partfrac(y)
            orbit.append(orbit[0]-1)
            a.append(0)
            while shift(f(x))>0:
                a[2] = a[2] + 1
                x = shift(f(x))
                if a[2]>A:
                    print 'approximatively rational'
                    return 1/a[1]
                if shift(f(x))==0:
                    print 'periodic point'
                    a[2]=a[2]+1
                    return rational_approximation(a,1)
            orbit.append(shift(x))    
            z = shift(f(x))
        
        for i in range(1,L):
            a.append(0)
            if shift(orbit[i+1])<shift(orbit[i]):
                while z>0:
                    y = z
                    z = first_return(shift(orbit[i]),shift(orbit[i+1]),z)
                    a[i+2]=a[i+2]+1
                    if a[i+2]>A:
                        print 'approximatively rational'
                        return rational_approximation(a,i)
                    if z==0:
                        print 'periodic point'
                        a[i+2]=a[i+2]+1
                        return rational_approximation(a,i+1)
            if shift(orbit[i+1])>shift(orbit[i]):
                 while z<0:
                     y = z
                     z = first_return(shift(orbit[i+1]),shift(orbit[i]),z)
                     a[i+2]=a[i+2]+1
                     if a[i+2]>A:
                         print 'approximatively rational'
                         return rational_approximation(a,i)
                    if z==0:
                         print 'periodic point'
                         a[i+2]=a[i+2]+1
                         return rational_approximation(a,i+1)
             orbit.append(y)
        
        print a
        return rational_approximation(a,L)
    

  [1]: http://www2.warwick.ac.uk/fac/sci/maths/people/staff/sebastian_van_strien/demelo-strien.pdf
  [2]: http://www.ams.org/mathscinet/search/publdoc.html?arg3=&co4=AND&co5=AND&co6=AND&co7=AND&dr=all&pg4=AUCN&pg5=TI&pg6=PC&pg7=ALLF&pg8=ET&r=1&review_format=html&s4=bruin&s5=numerical%20determination&s6=&s7=&s8=All&vfpref=html&yearRangeFirst=&yearRangeSecond=&yrop=eq
  [3]: http://www.sagemath.org/