(x, prec, rnd=round_fast, which=2)
return ci, si
Calculation of Ci(x), Si(x) for real x.

    which = 0 -- returns (Ci(x), -)
    which = 1 -- returns (Si(x), -)
    which = 2 -- returns (Ci(x), Si(x))

    Note: if x < 0, Ci(x) needs an additional imaginary term, pi*i.
Calculation of Ci(x), Si(x) for real x.
Calculation of Ci ( x ) Si ( x ) for real x .
def mpf_ci_si(x, prec, rnd=round_fast, which=2):
    """
    Calculation of Ci(x), Si(x) for real x.

    which = 0 -- returns (Ci(x), -)
    which = 1 -- returns (Si(x), -)
    which = 2 -- returns (Ci(x), Si(x))

    Note: if x < 0, Ci(x) needs an additional imaginary term, pi*i.
    """
    wp = prec + 20
    sign, man, exp, bc = x
    ci, si = None, None
    if not man:
        if x == fzero:
            return (fninf, fzero)
        if x == fnan:
            return (x, x)
        ci = fzero
        if which != 0:
            if x == finf:
                si = mpf_shift(mpf_pi(prec, rnd), -1)
            if x == fninf:
                si = mpf_neg(mpf_shift(mpf_pi(prec, negative_rnd[rnd]), -1))
        return (ci, si)
    # For small x: Ci(x) ~ euler + log(x), Si(x) ~ x
    mag = exp+bc
    if mag < -wp:
        if which != 0:
            si = mpf_perturb(x, 1-sign, prec, rnd)
        if which != 1:
            y = mpf_euler(wp)
            xabs = mpf_abs(x)
            ci = mpf_add(y, mpf_log(xabs, wp), prec, rnd)
        return ci, si
    # For huge x: Ci(x) ~ sin(x)/x, Si(x) ~ pi/2
    elif mag > wp:
        if which != 0:
            if sign:
                si = mpf_neg(mpf_pi(prec, negative_rnd[rnd]))
            else:
                si = mpf_pi(prec, rnd)
            si = mpf_shift(si, -1)
        if which != 1:
            ci = mpf_div(mpf_sin(x, wp), x, prec, rnd)
        return ci, si
    else:
        wp += abs(mag)
    # Use an asymptotic series? The smallest value of n!/x^n
    # occurs for n ~ x, where the magnitude is ~ exp(-x).
    asymptotic = mag-1 > math.log(wp, 2)
    # Case 1: convergent series near 0
    if not asymptotic:
        if which != 0:
            si = mpf_pos(mpf_ci_si_taylor(x, wp, 1), prec, rnd)
        if which != 1:
            ci = mpf_ci_si_taylor(x, wp, 0)
            ci = mpf_add(ci, mpf_euler(wp), wp)
            ci = mpf_add(ci, mpf_log(mpf_abs(x), wp), prec, rnd)
        return ci, si
    x = mpf_abs(x)
    # Case 2: asymptotic series for x >> 1
    xf = to_fixed(x, wp)
    xr = (MPZ_ONE<<(2*wp)) // xf   # 1/x
    s1 = (MPZ_ONE << wp)
    s2 = xr
    t = xr
    k = 2
    while t:
        t = -t
        t = (t*xr*k)>>wp
        k += 1
        s1 += t
        t = (t*xr*k)>>wp
        k += 1
        s2 += t
    s1 = from_man_exp(s1, -wp)
    s2 = from_man_exp(s2, -wp)
    s1 = mpf_div(s1, x, wp)
    s2 = mpf_div(s2, x, wp)
    cos, sin = mpf_cos_sin(x, wp)
    # Ci(x) = sin(x)*s1-cos(x)*s2
    # Si(x) = pi/2-cos(x)*s1-sin(x)*s2
    if which != 0:
        si = mpf_add(mpf_mul(cos, s1), mpf_mul(sin, s2), wp)
        si = mpf_sub(mpf_shift(mpf_pi(wp), -1), si, wp)
        if sign:
            si = mpf_neg(si)
        si = mpf_pos(si, prec, rnd)
    if which != 1:
        ci = mpf_sub(mpf_mul(sin, s1), mpf_mul(cos, s2), prec, rnd)
    return ci, si
def mpf_ci_si ( x , prec , rnd = round_fast , which = 2 ) : wp = prec + 20 sign , man , exp , bc = x ci , si = None , None if not man : if x == fzero : return ( fninf , fzero ) if x == fnan : return ( x , x ) ci = fzero if which != 0 : if x == finf : si = mpf_shift ( mpf_pi ( prec , rnd ) , - 1 ) if x == fninf : si = mpf_neg ( mpf_shift ( mpf_pi ( prec , negative_rnd [ rnd ] ) , - 1 ) ) return ( ci , si ) # For small x: Ci(x) ~ euler + log(x), Si(x) ~ x mag = exp + bc if mag < - wp : if which != 0 : si = mpf_perturb ( x , 1 - sign , prec , rnd ) if which != 1 : y = mpf_euler ( wp ) xabs = mpf_abs ( x ) ci = mpf_add ( y , mpf_log ( xabs , wp ) , prec , rnd ) return ci , si # For huge x: Ci(x) ~ sin(x)/x, Si(x) ~ pi/2 elif mag > wp : if which != 0 : if sign : si = mpf_neg ( mpf_pi ( prec , negative_rnd [ rnd ] ) ) else : si = mpf_pi ( prec , rnd ) si = mpf_shift ( si , - 1 ) if which != 1 : ci = mpf_div ( mpf_sin ( x , wp ) , x , prec , rnd ) return ci , si else : wp += abs ( mag ) # Use an asymptotic series? The smallest value of n!/x^n # occurs for n ~ x, where the magnitude is ~ exp(-x). asymptotic = mag - 1 > math . log ( wp , 2 ) # Case 1: convergent series near 0 if not asymptotic : if which != 0 : si = mpf_pos ( mpf_ci_si_taylor ( x , wp , 1 ) , prec , rnd ) if which != 1 : ci = mpf_ci_si_taylor ( x , wp , 0 ) ci = mpf_add ( ci , mpf_euler ( wp ) , wp ) ci = mpf_add ( ci , mpf_log ( mpf_abs ( x ) , wp ) , prec , rnd ) return ci , si x = mpf_abs ( x ) # Case 2: asymptotic series for x >> 1 xf = to_fixed ( x , wp ) xr = ( MPZ_ONE << ( 2 * wp ) ) // xf # 1/x s1 = ( MPZ_ONE << wp ) s2 = xr t = xr k = 2 while t : t = - t t = ( t * xr * k ) >> wp k += 1 s1 += t t = ( t * xr * k ) >> wp k += 1 s2 += t s1 = from_man_exp ( s1 , - wp ) s2 = from_man_exp ( s2 , - wp ) s1 = mpf_div ( s1 , x , wp ) s2 = mpf_div ( s2 , x , wp ) cos , sin = mpf_cos_sin ( x , wp ) # Ci(x) = sin(x)*s1-cos(x)*s2 # Si(x) = pi/2-cos(x)*s1-sin(x)*s2 if which != 0 : si = mpf_add ( mpf_mul ( cos , s1 ) , mpf_mul ( sin , s2 ) , wp ) si = mpf_sub ( mpf_shift ( mpf_pi ( wp ) , - 1 ) , si , wp ) if sign : si = mpf_neg ( si ) si = mpf_pos ( si , prec , rnd ) if which != 1 : ci = mpf_sub ( mpf_mul ( sin , s1 ) , mpf_mul ( cos , s2 ) , prec , rnd ) return ci , si