Ask Your Question
0

Convert from int to double and back in Cython

asked 2011-04-20 11:06:34 -0500

kcrisman gravatar image

This is a really basic question, I'm sure, for cython types.

In the notebook, I have the following cell to compute a famous but totally inefficient $\pi(x)$ (prime counting function). It already is about 100 times faster than the equivalent in native Python.

Question: The result += (fact - j*int(floor(t))) line still is bright yellow. I don't know if there is any easy way to mix these types. The problem is I need for t to take two ints and make them something I can apply floor to, which appears to be a double. Then I need to make the floor an int again to multiply it by j. Any ideas? Again, I'm sure this is simple, but I haven't earned my Cython stripes yet.

%cython
cdef extern from "math.h":
    double floor(double x)

cdef int primeish(int n):
    cdef int fact
    cdef int result
    cdef int j
    cdef double t
    if n==1:
        return 0
    elif n==2:
        return 1
    elif n==3:
        return 2
    else:
        result = -1
        fact = 1
        for j in range(3,n+1):
            fact = fact*(j-2)
            t = fact/j
            result += (fact - j*int(floor(t)))
        return result
edit retag flag offensive close merge delete

1 answer

Sort by ยป oldest newest most voted
2

answered 2011-04-20 12:40:20 -0500

benjaminfjones gravatar image

updated 2011-04-20 19:16:15 -0500

The statement

result += (fact - j*(<int>t))

will cast t to an int (dropping the fractional part like floor does) without calling floor. In the annotated cython code this change reduces the "yellowness" of that line. I'm not sure why it's yellow at all though, it seems like this statement should translate easily to pure C code.

If you want to use the C math library version of floor I think you need an appropriate include statement.


Edit:

If you use the @cython.cdivision(True) decorator, cython will not add exception checking for division by zero, this allows the division statement to be converted to one line of pure C code. In your code j is never zero, so this is okay. Also, you can skip the double t entirely if you just write fact/j. Since both are integers, the division occurs as integers and the result is the floor of the rational number fact/j.

Here is new code with the changes:

%cython
import cython

@cython.cdivision(True)    # turn division by zero checking off
def primeish(int n):
    cdef int fact
    cdef int result
    cdef int j
    cdef int t
    if n==1:
        return 0
    elif n==2:
        return 1
    elif n==3:
        return 2
    else:
        result = -1
        fact = 1
        for j in range(3,n+1):
            fact = fact*(j-2)
            t = fact / j    # integer division
            result += fact - t*j
        return result
edit flag offensive delete link more

Comments

Makes sense; I remember something about 'casting' from reading C tutorials...

kcrisman gravatar imagekcrisman ( 2011-04-20 15:06:38 -0500 )edit

Followup to make you earn your check mark :) - any way to make the t = fact/j turn white too? Obviously this can be done in C, but I don't know how to make Cython do it right. Thank you!

kcrisman gravatar imagekcrisman ( 2011-04-20 15:07:56 -0500 )edit

Actually, it doesn't matter, because I need too long of numbers for the code to work right. No choice but to use Sage integers, I guess :) This is not for high performance!

kcrisman gravatar imagekcrisman ( 2011-04-20 15:24:56 -0500 )edit

To answer your second comment: there is a cython compiler directive you can use `cdivision=True`. This turns off exception checking for division by zero. It's the exception checking that adds more C code and why the annotated code there is yellow. Look at http://docs.cython.org/src/reference/compilation.html

benjaminfjones gravatar imagebenjaminfjones ( 2011-04-20 18:35:59 -0500 )edit

Thanks, that makes sense. I'm just never sure what will behave like Python and what will behave like what I think (not know) C will do. As I said, even with longs this code becomes wrong after about n=22, and as it's just for demonstration I have no interest in importing mpz etc :)

kcrisman gravatar imagekcrisman ( 2011-04-21 03:17:24 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2011-04-20 11:06:34 -0500

Seen: 4,071 times

Last updated: Apr 20 '11