Sharedwww / talks / bernoulli / factorial.spyxOpen in CoCalc
Author: William A. Stein
1
r"""
2
This example illustrates how to use external includes and libraries
3
from spyx files. In particular, we show how to use GMP.
4
5
NOTE -- in the compilation anything from the following libraries
6
is linked in:
7
'mpfr', 'gmp', 'gmpxx', 'stdc++', 'pari', 'm', 'mwrank'
8
Modify sage/misc/pyrex.py to add more.
9
10
Note that *every* function from an external library that you use must
11
be explicitly declared. The declaration has to be good enough for the
12
Pyrex program to correctly convert this code to C. However, your C
13
compiler and the actual C code will only see the actual headers from
14
the original include files.
15
16
# TODO: The factorial below is actually three times as fast as
17
what's built into SAGE... which means I should move
18
this into SAGE.
19
20
EXAMPLES:
21
sage: time n = factorial(100000)
22
CPU times: user 6.93 s, sys: 0.00 s, total: 6.93 s
23
Wall time: 6.95
24
25
sage: time n=factorial_ZZ(100000)
26
CPU times: user 8.56 s, sys: 2.22 s, total: 10.79 s
27
Wall time: 10.83
28
29
sage: time n = magma('Factorial(100000)')
30
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
31
Wall time: 9.02
32
33
sage: time n=factorial_pure_python(100000)
34
CPU times: user 78.72 s, sys: 2.64 s, total: 81.37 s
35
Wall time: 94.89
36
37
"""
38
39
# If __no_preparse__ is in the file, the SAGE pre-processor
40
# is not run on the file, so the ^'s are xor and int literals
41
# are ints, etc.
42
43
#__no_preparse__
44
45
46
cdef extern from "gmp.h":
47
ctypedef void* mpz_t
48
void mpz_init(mpz_t integer)
49
void mpz_set_si(mpz_t integer, signed long int n)
50
void mpz_mul_si (mpz_t rop, mpz_t op1, long int op2)
51
char *mpz_get_str(char *str, int base, mpz_t op)
52
53
cdef extern from "stdlib.h":
54
void free(void *ptr)
55
56
57
def factorial_gmp(int n):
58
cdef mpz_t f
59
cdef int i
60
cdef char* s
61
62
mpz_init(f)
63
mpz_set_si(f, n)
64
65
# The _sig_on / _sig_off commands make it so
66
# this calculation is Ctrl-C interruptable.
67
# These should *only* be used to wrap code
68
# that compiles to pure C, i.e., doesn't access
69
# any of the Python library.
70
_sig_on
71
for i from 2 <= i <= n-1:
72
mpz_mul_si(f, f, i)
73
s = mpz_get_str(NULL, 32, f)
74
_sig_off
75
76
r = ZZ(0)
77
r.set_str(s, 32)
78
free(s)
79
return r
80
81
82
83
84