Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagelib
 Path: blob/master/sage/rings/bernmm.pyx
4072 views
1
r"""

2
Cython wrapper for bernmm library

3


4
AUTHOR: 

5
    - David Harvey (2008-06): initial version

6
"""

7


8
#*****************************************************************************

9
#       Copyright (C) 2008 William Stein <[email protected]>

10
#                     2008 David Harvey <[email protected]>

11
#

12
#  Distributed under the terms of the GNU General Public License (GPL)

13
#                  http://www.gnu.org/licenses/

14
#*****************************************************************************

15


16
include "../ext/cdefs.pxi"

17
include "../ext/interrupt.pxi"

18


19


20
cdef extern from "bernmm/bern_rat.h":

21
    void bern_rat "bernmm::bern_rat" (mpq_t res, long k, int num_threads)

22


23
cdef extern from "bernmm/bern_modp.h":

24
    long bern_modp "bernmm::bern_modp" (long p, long k)

25


26


27


28
from sage.rings.rational cimport Rational

29


30


31
def bernmm_bern_rat(long k, int num_threads = 1):

32
    r"""

33
    Computes k-th Bernoulli number using a multimodular algorithm.

34
    (Wrapper for bernmm library.)

35


36
    INPUT:

37
        k -- non-negative integer

38
        num_threads -- integer >= 1, number of threads to use

39


40
    COMPLEXITY:

41
        Pretty much quadratic in $k$. See the paper ``A multimodular algorithm

42
        for computing Bernoulli numbers'', David Harvey, 2008, for more details.

43


44
    EXAMPLES:

45
        sage: from sage.rings.bernmm import bernmm_bern_rat

46


47
        sage: bernmm_bern_rat(0)

48
        1

49
        sage: bernmm_bern_rat(1)

50
        -1/2

51
        sage: bernmm_bern_rat(2)

52
        1/6

53
        sage: bernmm_bern_rat(3)

54
        0

55
        sage: bernmm_bern_rat(100)

56
        -94598037819122125295227433069493721872702841533066936133385696204311395415197247711/33330

57
        sage: bernmm_bern_rat(100, 3)

58
        -94598037819122125295227433069493721872702841533066936133385696204311395415197247711/33330

59
        

60
    TESTS:

61
        sage: lst1 = [ bernoulli(2*k, algorithm='bernmm', num_threads=2) for k in [2932, 2957, 3443, 3962, 3973] ]

62
        sage: lst2 = [ bernoulli(2*k, algorithm='pari') for k in [2932, 2957, 3443, 3962, 3973] ]

63
        sage: lst1 == lst2

64
        True

65
        sage: [ Zmod(101)(t) for t in lst1 ]

66
        [77, 72, 89, 98, 86]

67
        sage: [ Zmod(101)(t) for t in lst2 ]

68
        [77, 72, 89, 98, 86]

69
    """

70
    cdef Rational x

71


72
    if k < 0:

73
        raise ValueError, "k must be non-negative"

74


75
    x = Rational()

76
    sig_on()

77
    bern_rat(x.value, k, num_threads)

78
    sig_off()

79


80
    return x

81


82


83
def bernmm_bern_modp(long p, long k):

84
    r"""

85
    Computes $B_k \mod p$, where $B_k$ is the k-th Bernoulli number.

86


87
    If $B_k$ is not $p$-integral, returns -1.

88


89
    INPUT:

90
        p -- a prime

91
        k -- non-negative integer

92


93
    COMPLEXITY:

94
        Pretty much linear in $p$.

95


96
    EXAMPLES:

97
        sage: from sage.rings.bernmm import bernmm_bern_modp

98


99
        sage: bernoulli(0) % 5, bernmm_bern_modp(5, 0)

100
        (1, 1)

101
        sage: bernoulli(1) % 5, bernmm_bern_modp(5, 1)

102
        (2, 2)

103
        sage: bernoulli(2) % 5, bernmm_bern_modp(5, 2)

104
        (1, 1)

105
        sage: bernoulli(3) % 5, bernmm_bern_modp(5, 3)

106
        (0, 0)

107
        sage: bernoulli(4), bernmm_bern_modp(5, 4)

108
        (-1/30, -1)

109
        sage: bernoulli(18) % 5, bernmm_bern_modp(5, 18)

110
        (4, 4)

111
        sage: bernoulli(19) % 5, bernmm_bern_modp(5, 19)

112
        (0, 0)

113


114
        sage: p = 10000019; k = 1000

115
        sage: bernoulli(k) % p

116
        1972762

117
        sage: bernmm_bern_modp(p, k)

118
        1972762

119


120
    """

121
    cdef long x

122


123
    if k < 0:

124
        raise ValueError, "k must be non-negative"

125


126
    sig_on()

127
    x = bern_modp(p, k)

128
    sig_off()

129


130
    return x

131
    

132


133
# ============ end of file

134


135