Index: sage/algebras/quaternion_algebra.py =================================================================== --- sage/algebras/quaternion_algebra.py (revision 4069) +++ sage/algebras/quaternion_algebra.py (revision 5522) @@ -36,4 +36,5 @@ from sage.modules.free_module import VectorSpace from sage.matrix.matrix_space import MatrixSpace +from sage.matrix.matrix0 import Matrix from sage.algebras.free_algebra import FreeAlgebra from sage.algebras.free_algebra_quotient import FreeAlgebraQuotient @@ -53,11 +54,19 @@ def fundamental_discriminant(D): + """ + Return the discriminant of the quadratic extension K=Q(\sqrt{D}), i.e. + an integer d congruent to either 0 or 1, mod 4, and such that, at most, + the only square dividing it is 4. + """ D = Integer(D) D = D.square_free_part() - if D%4 in (0,1): + if D%4 == 1: return D return 4*D def ramified_primes(a,b): + """ + Return a list of the finite primes ramifying in Q(a,b) + """ a = Integer(a); b = Integer(b) if a.is_square() or b.is_square() or (a+b).is_square(): @@ -212,5 +221,7 @@ D = T**2 - 4*N try: - S = D.square_root() + S = D.sqrt() + except AttributeError: + raise ValueError, "%s is not an integer"%D except ArithmeticError: raise ValueError, "Invalid inner product input." @@ -260,8 +271,7 @@ if not isinstance(gram,Matrix) or not gram.is_symmetric: raise AttributeError, "Argument gram (= %s) must be a symmetric matrix"%gram - (q0,t01,t02,t03,_,q1,t12,t13,_,_,q3,t23,_,_,_,q4) = gram.list() - n1 = q1/2 - n2 = q2/2 - n3 = q3/2 + (q0,t01,t02,t03,t10,q1,t12,t13,t20,t21,q2,t23,t30,t31,t32,q3) = gram.list() + norms = [q1/2, q2/2, q3/2] + traces = [t10, t20, t30, t12, t12, t23] return QuaternionAlgebraWithInnerProduct(K,norms,traces,names=names) Index: sage/all_cmdline.py =================================================================== --- sage/all_cmdline.py (revision 5292) +++ sage/all_cmdline.py (revision 5513) @@ -21,5 +21,5 @@ print t if 'type object' in str(msg): - msg = str(msg) + '\n\n** In SAGE, the easiest fix for this problem is to type "sage -ba"\n to rebuild all the SageX code (this takes several minutes).\n Alternatively, touch the last .pyx file in the traceback above. **\n' + msg = str(msg) + '\n\n** In SAGE, the easiest fix for this problem is to type "sage -ba"\n to rebuild all the Cython code (this takes several minutes).\n Alternatively, touch the last .pyx file in the traceback above. **\n' raise ValueError, msg Index: sage/calculus/calculus.py =================================================================== --- sage/calculus/calculus.py (revision 5482) +++ sage/calculus/calculus.py (revision 5507) @@ -734,8 +734,9 @@ return long(int(self)) - def numerical_approx(self, prec=53): + def numerical_approx(self, prec=None, digits=None): r""" Return a numerical approximation of self as either a real or - complex number. + complex number with at least the requested number of bits or + digits of precision. NOTE: You can use \code{foo.n()} as a shortcut for @@ -743,5 +744,6 @@ INPUT: - prec -- integer (default: 53): the number of bits of precision + prec -- an integer: the number of bits of precision + digits -- an integer: digits of precision OUTPUT: @@ -760,4 +762,6 @@ sage: cos(3).numerical_approx(200) -0.98999249660044545727157279473126130239367909661558832881409 + sage: numerical_approx(cos(3), digits=10) + -0.9899924966 sage: (i + 1).numerical_approx(32) 1.00000000 + 1.00000000*I @@ -765,6 +769,11 @@ 7.2740880444219335226246195788 """ + if prec is None: + if digits is None: + prec = 53 + else: + prec = int(digits * 3.4) + 2 + # make sure the field is of the right precision - prec = Integer(prec) field = RealField(prec) Index: sage/combinat/combinat.py =================================================================== --- sage/combinat/combinat.py (revision 5438) +++ sage/combinat/combinat.py (revision 5525) @@ -204,7 +204,10 @@ from sage.libs.all import pari +import expnums +import partitions as partitions_ext + ######### combinatorial sequences -def bell_number(n): +def bell_number(n, algorithm='sage'): r""" Returns the n-th Bell number (the number of ways to partition a @@ -213,11 +216,14 @@ INPUT: n -- an integer + algorithm -- 'sage': use N. Alexander's custom implementation in SAGE + 'gap': use Gap's Bell command (slow) If $n \leq 0$, returns $1$. - - Wraps GAP's Bell. + EXAMPLES: sage: bell_number(10) + 115975 + sage: bell_number(10, algorithm='gap') 115975 sage: bell_number(2) @@ -231,7 +237,20 @@ ... TypeError: no coercion of this rational to integer - """ - ans=gap.eval("Bell(%s)"%ZZ(n)) - return eval(ans) + + To compute all Bell numbers up to n, use expnums: + sage: expnums(10,1) + [1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147] + sage: [bell_number(n) for n in range(10)] + [1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147] + """ + n = ZZ(n) + if n <= 0: + return ZZ(1) + if algorithm == 'sage': + return expnums.expnums(n+1,1)[-1] + elif algorithm == 'gap': + return ZZ(gap.eval("Bell(%s)"%ZZ(n))) + else: + raise ValueError, "unknown algorithm '%s'"%algorithm ## def bernoulli_number(n): @@ -1168,5 +1187,7 @@ the positive integer n into sums with k summands. algorithm -- (default: 'gap') - 'gap' -- use GAP + 'gap' -- use GAP (VERY *slow*) + 'bober' -- use Jonathon Bober's implementation (*very* fast, + but new and not well tested yet). 'pari' -- use PARI. Speed seems the same as GAP until $n$ is in the thousands, in which case PARI is faster. *But* @@ -1193,4 +1214,6 @@ 7 sage: number_of_partitions(5, algorithm='pari') + 7 + sage: number_of_partitions(5, algorithm='bober') 7 @@ -1249,4 +1272,6 @@ ans=gap.eval("NrPartitions(%s,%s)"%(ZZ(n),ZZ(k))) return ZZ(ans) + elif algorithm == 'bober': + return partitions_ext.number_of_partitions(n) elif algorithm == 'pari': if not k is None: Index: sage/combinat/partitions.pyx =================================================================== --- sage/combinat/partitions.pyx (revision 5525) +++ sage/combinat/partitions.pyx (revision 5525) @@ -0,0 +1,38 @@ +""" +Number of partitions of integer + +AUTHOR: + -- William Stein (2007-07-28): initial version + -- Jonathon Bober (2007-07-28): wrote the program partitions_c.cc + that does all the actual heavy lifting. +""" + +import sys + +cdef extern from "gmp.h": + ctypedef void* mpz_t + +cdef extern from "partitions_c.h": + int part(mpz_t answer, unsigned int n) + +include "../ext/interrupt.pxi" + +from sage.rings.integer cimport Integer + +def number_of_partitions(n): + cdef unsigned int nn = n + if n != nn: + raise ValueError, "input must be a nonnegative integer less than %s"%(2*sys.maxint) + + if nn <= 0: + return Integer(0) + elif nn == 1: + return Integer(1) # part hangs on n=1 as input. + + cdef Integer ans = Integer(0) + + _sig_on + part(ans.value, nn) + _sig_off + + return ans Index: sage/combinat/partitions_c.cc =================================================================== --- sage/combinat/partitions_c.cc (revision 5526) +++ sage/combinat/partitions_c.cc (revision 5526) @@ -0,0 +1,694 @@ +/* + Program to compute the number of partitions of n. + + Author: Jonathan Bober + Version: .1 (7/28/2007) + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ +#include + +#include + +#include +#include +#include + +#include +#include + +const bool debug = false; +//const bool debug = true; + +const bool debugs = false; + +//const bool debugf = true; +const bool debugf = false; +//const bool debuga = true; +const bool debuga = false; +//const bool debugt = true; +const bool debugt = false; + +using std::cout; +using std::endl; + +const unsigned int min_prec = 53; +const unsigned int min_precision = 53; + +long GCD(long a, long b); + +unsigned int calc_precision(unsigned int n, unsigned int N); + +//mpfr versions + +mpz_t ztemp1, ztemp2, ztemp3; + +mpfr_t mp_one_over_12, mp_one_over_24, mp_sqrt2, mp_sqrt3, mp_pi; +mpfr_t mp_A, mp_B, mp_C, mp_D; +mp_rnd_t round_mode = GMP_RNDN; + +mpfr_t tempa1, tempa2, tempf1, tempf2, temps1, temps2, tempt1, tempt2; // temp variables for different functions, with precision set and cleared by mp_set_precision + +bool mp_vars_initialized = false; + +mp_prec_t mp_precision; + +void mp_init(unsigned int prec, unsigned int n); +//void mp_f_precompute(unsigned int n); +void mp_f(mpfr_t result, unsigned int k); +void mp_s(mpfr_t result, unsigned int j, unsigned int q); +void mp_a(mpfr_t result, unsigned int n, unsigned int k); +void mp_t(mpfr_t result, unsigned int n); + +unsigned int compute_initial_precision(unsigned int n); // computes the precision required to accurately compute p(n) +unsigned int compute_current_precision(unsigned int n, unsigned int N); // computed the precision required to + // accurately compute the tail of the rademacher series + // assuming that N terms have already been computed + +double compute_remainder(unsigned int n, unsigned int N); // Gives an upper bound on the error that occurs + // when only N terms of the Rademacher series have been + // computed. NOTE: should only be called when we already know + // that the error is small (ie, compute_current_precision returns + // a small number, eg, something < 32) + +//low precision (double) versions of the functions: + +double df_A, df_B, df_C, df_D; + +void d_f_precompute(unsigned int n); +double d_f(unsigned int k); +double d_s(unsigned int h,unsigned int q); +double d_A(unsigned int n, unsigned int k); +double d_t(unsigned int n, unsigned int N); + +unsigned int compute_initial_precision(unsigned int n) { + // We just want to know how many bits we will need to + // compute to get an accurate answer. + + // We know that + // + // p(n) ~ exp(pi * sqrt(2n/3))/(4n sqrt(3)), + // + // so for now we are assuming that p(n) < exp(pi * sqrt(2n/3))/n, + // so we need pi*sqrt(2n/3)/log(2) - log(n)/log(2) + EXTRA bits. + // + // EXTRA should depend on n, and should be something that ensures + // that the TOTAL ERROR in all computations is < (something small). + // This needs to be worked out carefully. EXTRA = log(n)/log(2) + 3 + // is probably good enough, and is convenient... + // + // but we really need: + // + // p(n) < something + // + // to be sure that we compute the correct answer + + /* I don't think that we need to be that careful about this. + * + mpfr_t bits; + mpfr_t t1; // we refer to t1 as temp1 + + mpfr_init2(bits, 32); // use 32 bits of precision to compute the approximation + + mpfr_init2(t1, 32); + + // Error = exp( pi* sqrt(2n/3)); + + mpfr_set_ui(t1, 2, round_mode); // temp1 = 2 + mpfr_mul_ui(t1, t1, n, round_mode); // temp1 = 2n + mpfr_div_ui(t1, t1, 3, round_mode); // temp1 = 2n/3 + mpfr_sqrt(t1, t1, round_mode); // temp1 = sqrt(2n/3) + + mpfr_const_pi(bits, round_mode); // bits = pi + + mpfr_mul(bits, bits, t1, round_mode); // bits = pi * sqrt(2n/3) + + mpfr_const_log2(t1, round_mode); // temp1 = log(2) + mpfr_div(bits, bits, t1, round_mode); // bits = pi * sqrt(2n/3)/log(2) + */ + unsigned int result = (unsigned int)(ceil(3.1415926535897931 * sqrt(2.0 * double(n)/ 3.0) / log(2))) + 3; + if(debug) + cout << "Using initial precision of " << result << " bits." << endl; + if(result > min_precision) { + return result; + } + else return min_precision; + +} + +unsigned int compute_current_precision(unsigned int n, unsigned int N) { + // we want to compute + // log(A/sqrt(N) + B*sqrt(N/(n-1))*sinh(C * sqrt(n) / N) / log(2) + // + // maybe there is a better way... + + // cout << N; + + mpfr_t A, B, C; + mpfr_init2(A, 32); + mpfr_init2(B, 32); + mpfr_init2(C, 32); + + mpfr_set_d(A,1.11431833485164,round_mode); + mpfr_set_d(B,0.059238439175445,round_mode); + mpfr_set_d(C,2.5650996603238,round_mode); + + mpfr_t error, t1, t2; + mpfr_init2(error, 32); // we shouldn't need much precision here since we just need the most significant bit + mpfr_init2(t1, 32); + mpfr_init2(t2, 32); + + mpfr_set(error, A, round_mode); // error = A + mpfr_sqrt_ui(t1, N, round_mode); // t1 = sqrt(N) + mpfr_div(error, error, t1, round_mode); // error = A/sqrt(N) + + + mpfr_sqrt_ui(t1, n, round_mode); // t1 = sqrt(n) + mpfr_mul(t1, t1, C, round_mode); // t1 = C * sqrt(n) + mpfr_div_ui(t1, t1, N, round_mode); // t1 = C * sqrt(n) / N + mpfr_sinh(t1, t1, round_mode); // t1 = sinh( ditto ) + mpfr_mul(t1, t1, B, round_mode); // t1 = B * sinh( ditto ) + + mpfr_set_ui(t2, N, round_mode); // t2 = N + mpfr_div_ui(t2, t2, n-1, round_mode); // t2 = N/(n-1) + mpfr_sqrt(t2, t2, round_mode); // t2 = sqrt( ditto ) + + mpfr_mul(t1, t1, t2, round_mode); // t1 = B * sqrt(N/(n-1)) * sinh(C * sqrt(n)/N) + + mpfr_add(error, error, t1, round_mode); // error = (ERROR ESTIMATE) + + unsigned int p = mpfr_get_exp(error) + 3; // I am almost certain that this does the right thing + // (The 3 is for good luck.) + + p = p + (unsigned int)ceil(log(n)/log(2)); + + if(debug) { + cout << "Error seems to be: "; + mpfr_out_str(stdout, 10, 0, error, round_mode); + cout << endl; + cout.flush(); + cout << "Switching to precision of " << p << " bits. " << endl; + } + + mpfr_clear(error); + mpfr_clear(t1); + mpfr_clear(t2); + mpfr_clear(A); + mpfr_clear(B); + mpfr_clear(C); + + + + + if(p > min_precision) { + return p; // don't want to return < 32 + // (for now, at least -- we can be more careful) + } + return min_precision; +} + +double compute_remainder(unsigned int n, unsigned int N) { + double A = 1.11431833485164; + double B = 0.059238439175445; + double C = 2.5650996603238; + double result; + result = A/sqrt(N) + B * sqrt(double(N)/double(n-1))*sinh(C * sqrt(double(n))/double(N)); + return result; +} + + +void mp_set_precision(unsigned int prec) { + static bool init = false; + if(init) { + mpfr_clear(tempa1); + mpfr_clear(tempa2); + mpfr_clear(tempf1); + mpfr_clear(tempf2); + mpfr_clear(tempt1); + mpfr_clear(tempt2); + mpfr_clear(temps1); + mpfr_clear(temps2); + } + mpfr_init2(tempa1, prec); + mpfr_init2(tempa2, prec); + mpfr_init2(tempf1, prec); + mpfr_init2(tempf2, prec); + mpfr_init2(tempt1, prec); + mpfr_init2(tempt2, prec); + mpfr_init2(temps1, prec); + mpfr_init2(temps2, prec); + + init = true; +} + +void mp_init(unsigned int prec, unsigned int n) { + static bool init = false; + mp_precision = prec; + mp_prec_t p = mp_precision; + + if(init) { + mpfr_clear(mp_one_over_12); mpfr_clear(mp_one_over_24); mpfr_clear(mp_sqrt2); mpfr_clear(mp_sqrt3); mpfr_clear(mp_pi); + mpfr_clear(mp_A); mpfr_clear(mp_B); mpfr_clear(mp_C); mpfr_clear(mp_D); + + mpz_clear(ztemp1); + mpz_clear(ztemp2); + mpz_clear(ztemp3); + } + mpfr_init2(mp_one_over_12,p); mpfr_init2(mp_one_over_24,p); mpfr_init2(mp_sqrt2,p); mpfr_init2(mp_sqrt3,p); mpfr_init2(mp_pi,p); + mpfr_init2(mp_A,p); mpfr_init2(mp_B,p); mpfr_init2(mp_C,p); mpfr_init2(mp_D,p); + mpz_init(ztemp1); + mpz_init(ztemp2); + mpz_init(ztemp3); + + init = true; + + mpfr_set_ui(mp_one_over_12, 1, round_mode); // mp_one_over_12 = 1/12 + mpfr_div_ui(mp_one_over_12, mp_one_over_12, 12, round_mode); // + + mpfr_set_ui(mp_one_over_24, 1, round_mode); // mp_one_over_24 = 1/24 + mpfr_div_ui(mp_one_over_24, mp_one_over_24, 24, round_mode); // + + mpfr_t n_minus; // + mpfr_init2(n_minus, p); // + mpfr_set_ui(n_minus, n, round_mode); // n_minus = n + mpfr_sub(n_minus, n_minus, mp_one_over_24,round_mode); // n_minus = n - 1/24 + + mpfr_t sqrt_n_minus; // + mpfr_init2(sqrt_n_minus, p); // + mpfr_sqrt(sqrt_n_minus, n_minus, round_mode); // n_minus = sqrt(n-1/24) + + + mpfr_sqrt_ui(mp_sqrt2, 2, round_mode); // mp_sqrt2 = sqrt(2) + mpfr_sqrt_ui(mp_sqrt3, 3, round_mode); // mp_sqrt3 = sqrt(3) + mpfr_const_pi(mp_pi, round_mode); // mp_pi = pi + + //mp_A = sqrt(2) * 3.1415926535897931 * sqrt(n - 1.0/24.0);----------- + // + mpfr_set(mp_A, mp_sqrt2, round_mode); // mp_A = sqrt(2) + mpfr_mul(mp_A, mp_A, mp_pi, round_mode); // mp_A = sqrt(2) * pi + mpfr_mul(mp_A, mp_A, sqrt_n_minus, round_mode); // mp_A = sqrt(2) * pi * sqrt(n - 1/24) + //-------------------------------------------------------------------- + + //cout << "n_minus_1/24: "; + //mpfr_out_str(stdout, 10, 20, n_minus, round_mode); + //cout << endl; + + //mp_B = 2.0 * sqrt(3) * (n - 1.0/24.0);------------------------------ + mpfr_set_ui(mp_B, 2, round_mode); // mp_A = 2 + mpfr_mul(mp_B, mp_B, mp_sqrt3, round_mode); // mp_A = 2*sqrt(3) + mpfr_mul(mp_B, mp_B, n_minus, round_mode); // mp_A = 2*sqrt(3)*(n-1/24) + //-------------------------------------------------------------------- + + //mp_C = sqrt(2) * pi * sqrt(n - 1.0/24.0) / sqrt(3);----------------- + mpfr_set(mp_C, mp_sqrt2, round_mode); // mp_C = sqrt(2) + mpfr_mul(mp_C, mp_C, mp_pi, round_mode); // mp_C = sqrt(2) * pi + mpfr_mul(mp_C, mp_C, sqrt_n_minus, round_mode); // mp_C = sqrt(2) * pi * sqrt(n - 1/24) + mpfr_div(mp_C, mp_C, mp_sqrt3, round_mode); // mp_C = sqrt(2) * pi * sqrt(n - 1/24) / sqrt3 + //-------------------------------------------------------------------- + + //mp_D = 2.0 * (n - 1.0/24.0) * sqrt(n - 1.0/24.0);------------------- + mpfr_set_ui(mp_D, 2, round_mode); // mp_D = 2 + mpfr_mul(mp_D, mp_D, n_minus, round_mode); // mp_D = 2 * (n - 1/24) + mpfr_mul(mp_D, mp_D, sqrt_n_minus, round_mode); // mp_D = 2 * (n - 1/24) * sqrt(n - 1/24) + //-------------------------------------------------------------------- + + mpfr_clear(n_minus); + mpfr_clear(sqrt_n_minus); + + + // some double precision versions of the above values + df_A = sqrt(2) * 3.1415926535897931 * sqrt(n - 1.0/24.0); + df_B = 2.0 * sqrt(3) * (n - 1.0/24.0); + df_C = sqrt(2) * 3.1415926535897931 * sqrt(n - 1.0/24.0) / sqrt(3); + df_D = 2.0 * (n - 1.0/24.0) * sqrt(n - 1.0/24.0); + +} + +void mp_f(mpfr_t result, unsigned int k) { + //result = pi * sqrt(2) * cosh(A/(sqrt(3)*k))/(B*k) - sinh(C/k)/D; + mpfr_set(result, mp_pi, round_mode); // result = pi + mpfr_mul(result, result, mp_sqrt2, round_mode); // result = sqrt(2) * pi + + mpfr_div(tempf1, mp_A, mp_sqrt3, round_mode); // temp1 = mp_A/sqrt(3) + mpfr_div_ui(tempf1, tempf1, k, round_mode); // temp1 = mp_A/(sqrt(3) * k) + mpfr_cosh(tempf1, tempf1, round_mode); // temp1 = cosh(mp_A/(sqrt(3) * k)) + mpfr_div(tempf1, tempf1, mp_B, round_mode); // temp1 = cosh(mp_A/(sqrt(3) * k))/mp_B + mpfr_div_ui(tempf1, tempf1, k, round_mode); // temp1 = cosh(mp_A/(sqrt(3) * k))/(mp_B*k) + + mpfr_mul(result, result, tempf1, round_mode); // result = sqrt(2) * pi * cosh(mp_A/(sqrt(3) * k))/(mp_B*k) + + mpfr_div_ui(tempf1, mp_C, k, round_mode); // temp1 = mp_C/k + mpfr_sinh(tempf1, tempf1, round_mode); // temp1 = sinh(mp_C/k) + mpfr_div(tempf1, tempf1, mp_D, round_mode); // temp1 = sinh(mp_C/k)/D + + mpfr_sub(result, result, tempf1, round_mode); // result = RESULT! + //return 3.1415926535897931 * sqrt(2) * cosh(df_A/(sqrt(3)*k))/(df_B*k) - sinh(df_C/k)/df_D; +} + +void mp_s(mpfr_t result, unsigned int h,unsigned int q) { + if(q < 3) { + mpfr_set_ui(result, 0, round_mode); + return; + } + + //double result, R1, R2, temp1, temp2; + unsigned int n, r1, r2, temp3 = 0; + + if(h == 1) { + if(q < (unsigned int)(sqrt(UINT_MAX))) { + //elternate, assuming no overflow + mpfr_set_ui(result, (q-1)*(q-2), round_mode); + mpfr_div_ui(result, result, 12*q, round_mode); + } + else { + mpz_set_ui(ztemp1, q-1); // temp = q-1 + mpz_mul_ui(ztemp1, ztemp1, q-2); // temp = (q-1)(q-2) + + mpfr_set_z(result, ztemp1, round_mode); // result = (q-1)(q-2) + mpz_set_ui(ztemp1, q); // temp = q + mpz_mul_ui(ztemp1, ztemp1, 12); // temp = 12q + mpfr_div_z(result, result, ztemp1, round_mode); // result = (q-1)(q-2)/12q + } +/* + //result = (q-1)*(q-2)/(12*q); + mpfr_set_ui(result, q-1, round_mode); + mpfr_mul_ui(result, result, q-2, round_mode); + mpfr_div_ui(result, result, q, round_mode); // in this step we don't want to assume that 12*q will not overflow + mpfr_div_ui(result, result, 12, round_mode); +*/ + + + return; + } + + mpfr_set_ui(result, 0, round_mode); // result = 0 + + if(debugs) + if(h > q) { + cout << "oops in mp_s()" << endl; + exit(1); + } + r1 = q; + r2 = h; + + n = 0; + while(r1 > 0 && r2 > 0) { + if(r1 < (unsigned int)(sqrt(UINT_MAX)/2.0)) { // if r1 is small enough we can use + // standard C integers + // NOTE: squareroot computation should be optimized by the compiler. + mpfr_set_ui(temps1, r1*r1 + r2*r2 + 1, round_mode); + mpfr_div_ui(temps1, temps1, r1 * r2, round_mode); + + } + else { + //temp1 = (R1*R1 + R2*R2 + 1.0)/(R1 * R2); // again, we are NOT going to assume no overflow + + mpfr_set_ui(temps1, r1, round_mode); // temp1 = r1 + mpfr_mul_ui(temps1, temps1, r1, round_mode); // temp1 = r1 * r1 + + mpfr_set_ui(temps2, r2, round_mode); // temp2 = r2 + mpfr_mul_ui(temps2, temps2, r2, round_mode); // temp2 = r2 * r2 + + mpfr_add(temps1, temps1, temps2, round_mode); // temp1 = r1*r1 + r2*r2 + mpfr_add_ui(temps1, temps1, 1, round_mode); // temp1 = r1*r1 + r2*r2 + 1 + + mpfr_div_ui(temps1, temps1, r1, round_mode); // temp1 = (r1*r1 + r2*r2 + 1)/r1 + mpfr_div_ui(temps1, temps1, r2, round_mode); // temp1 = (r1*r1 + r2*r2 + 1)/(r1 * r2) + } + if(n % 2 == 0){ + mpfr_add(result, result, temps1, round_mode); // result += temp1; + } + else { + mpfr_sub(result, result, temps1, round_mode); // result -= temp1; + } + temp3 = r1 % r2; + r1 = r2; + r2 = temp3; + n++; + } + + //cout << temps1 << endl; + //cout << result << endl; + + mpfr_mul(result, result, mp_one_over_12, round_mode); // result = result * 1.0/12.0; + + + if(n % 2 == 1) { + mpfr_set_d(temps1, .25, round_mode); + mpfr_sub(result, result, temps1, round_mode); // result = result - .25; + } + + + + //return result; +} + + +void mp_a(mpfr_t result, unsigned int n, unsigned int k) { + + if (k == 1) { + mpfr_set_ui(result, 1, round_mode); //result = 1 + return; + } + + mpfr_set_ui(result, 0, round_mode); + + unsigned int h = 0; + for(h = 1; h < k+1; h++) { + if(GCD(h,k) == 1) { + //result += cos(pi * ( s(h,k) - (2.0 * h * n)/k) ); + mp_s(tempa1, h, k); // temp1 = s(h,k) + + if(debugs) { + cout << "s(" << h << "," << k << ") = "; + mpfr_out_str(stdout, 10, 0, tempa1, round_mode); + cout << endl; + } + + mpfr_set_ui(tempa2, 2, round_mode); // temp2 = 2 + mpfr_mul_ui(tempa2, tempa2, h, round_mode); // temp2 = 2h + mpfr_mul_ui(tempa2, tempa2, n, round_mode); // temp2 = 2hn + mpfr_div_ui(tempa2, tempa2, k, round_mode); // temp2 = 2hn/k + + mpfr_sub(tempa1, tempa1, tempa2, round_mode); // temp1 = s(h,k) - 2hn/k + mpfr_mul(tempa1, tempa1, mp_pi, round_mode); // temp1 = pi * (s(h,k) - 2hn/k) + mpfr_cos(tempa1, tempa1, round_mode); // temp1 = cos( ditto ) + + mpfr_add(result, result, tempa1, round_mode); // result = result + temp1 + } + } + +} + + +void mp_t(mpfr_t result, unsigned int n) { + // NOTE: result should NOT have been initialized when this is called, + // as we initialize it to the proper precision in this function. + + unsigned int initial_precision = compute_initial_precision(n); + mpfr_t t1, t2; + mpfr_init2(t1, initial_precision); + mpfr_init2(t2, initial_precision); + mpfr_init2(result, initial_precision); + mpfr_set_ui(result, 0, round_mode); + + mp_init(initial_precision, n); + mp_set_precision(initial_precision); + + unsigned int current_precision = initial_precision; + unsigned int new_precision; + + double remainder = 100.0; +// d_f_precompute(n); + unsigned int k = 1; + for(k = 1; current_precision > min_precision; k++) { + mpfr_sqrt_ui(t1, k, round_mode); // t1 = sqrt(k) + + mp_a(t2, n, k); // t2 = A_k(n) + + if(debuga) { + cout << "a(" << k << ") = "; + mpfr_out_str(stdout, 10, 0, t2, round_mode); + cout << endl; + } + + mpfr_mul(t1, t1, t2, round_mode); // t1 = sqrt(k)*A_k(n) + + mp_f(t2, k); // t2 = f_k(n) + + if(debugf) { + cout << "f(" << k << ") = "; + mpfr_out_str(stdout, 10, 0, t2, round_mode); + cout << endl; + } + + mpfr_mul(t1, t1, t2, round_mode); // t1 = sqrt(k)*A_k(n)*f_k(n) + + mpfr_add(result, result, t1, round_mode); // result += summand + + if(debugt) { + cout << "Partial sum " << k << " = "; + mpfr_out_str(stdout, 10, 0, result, round_mode); + cout << endl; + } + //K++; + //temp = sqrt(K) * d_A(n,k) * d_f(k); + //result += temp; + + new_precision = compute_current_precision(n,k); + if(new_precision != current_precision) { + current_precision = new_precision; + mp_set_precision(current_precision); + mpfr_clear(t1); mpfr_clear(t2); + mpfr_init2(t1, current_precision); + mpfr_init2(t2, current_precision); + } + } + + double result2 = 0; + + for( ; remainder > .5; k++) { + result2 += sqrt(k) * d_A(n,k) * d_f(k); + remainder = compute_remainder(n,k); + } + + mpfr_set_d(t1, result2, round_mode); + mpfr_add(result, result, t1, round_mode); + +// result = result/(3.1415926535897931 * sqrt(2)); + mpfr_div(result, result, mp_pi, round_mode); + mpfr_div(result, result, mp_sqrt2, round_mode); +} + +// double versions of the functions + +void d_f_precompute(unsigned int n) { + df_A = sqrt(2) * 3.1415926535897931 * sqrt(n - 1.0/24.0); + df_B = 2.0 * sqrt(3) * (n - 1.0/24.0); + df_C = sqrt(2) * 3.1415926535897931 * sqrt(n - 1.0/24.0) / sqrt(3); + df_D = 2.0 * (n - 1.0/24.0) * sqrt(n - 1.0/24.0); +} + +double d_f(unsigned int k) { + return 3.1415926535897931 * sqrt(2) * cosh(df_A/(sqrt(3)*k))/(df_B*k) - sinh(df_C/k)/df_D; +} + +double d_s(unsigned int h,unsigned int q) { + if(q < 3) { + return 0.0; + } + + double result, R1, R2, temp1, temp2; + unsigned int n, r1, r2, temp3 = 0; + + if(h == 1) { + double Q; + Q = q; + result = (Q-1)*(Q-2)/(12*Q); + return result; + } + + result = 0; + R1 = q; + R2 = h; + + r1 = q; + r2 = h; + + n = 0; + while(r1 && r2) { + temp1 = (R1*R1 + R2*R2 + 1.0)/(R1 * R2); + if(n % 2 == 0){ + result += temp1; + } + else { + result -= temp1; + } + temp3 = r1 % r2; + r1 = r2; + r2 = temp3; + R1 = r1; + R2 = r2; + n++; + } + + result = result * 1.0/12.0; + + if(n % 2 == 1) { + result = result - .25; + } + + return result; +} + +double d_A(unsigned int n, unsigned int k) { + double result; + result = 0; + + if (k == 1) { + return 1.0; + } + + unsigned int h = 0; + for(h = 1; h < k+1; h++) { + if(GCD(h,k) == 1) { + // cout << "s(" << h << "," << k << ") = " << d_s(h,k) << endl; + // cout << "s(" << h << "," << k << ") = " << d_s(h,k) << endl; + result += cos( 3.1415926535897931 * ( d_s(h,k) - (2.0 * h * n)/double(k)) ); + } + } + return result; +} + +long GCD(long a, long b) +{ + long u, v, t, x; + + if (a < 0) { + a = -a; + } + + if (b < 0) { + b = -b; + } + + + if (b==0) + x = a; + else { + u = a; + v = b; + do { + t = u % v; + u = v; + v = t; + } while (v != 0); + + x = u; + } + + return x; +} + + +/* answer must have already been mpz_init'd. */ +int part(mpz_t answer, unsigned int n){ + mpfr_t result; + + mp_t(result, n); + + mpfr_get_z(answer, result, round_mode); + + return 0; +} Index: sage/combinat/partitions_c.h =================================================================== --- sage/combinat/partitions_c.h (revision 5520) +++ sage/combinat/partitions_c.h (revision 5520) @@ -0,0 +1,3 @@ +#include + +int part(mpz_t answer, unsigned int n); Index: sage/dsage/scripts/dsage_setup.py =================================================================== --- sage/dsage/scripts/dsage_setup.py (revision 5381) +++ sage/dsage/scripts/dsage_setup.py (revision 5524) @@ -79,10 +79,10 @@ def setup_server(template=None): check_dsage_dir() - template_dict = {'organization': 'SAGE (at %s)'%(socket.gethostname()), + template_dict = {'organization': 'SAGE (at %s)'%(socket.getfqdn()), 'unit': '389', 'locality': None, 'state': 'Washington', 'country': 'US', - 'cn': socket.gethostname(), + 'cn': socket.getfqdn(), 'uid': 'sage_user', 'dn_oid': None, Index: sage/graphs/graph.py =================================================================== --- sage/graphs/graph.py (revision 5445) +++ sage/graphs/graph.py (revision 5513) @@ -12,5 +12,5 @@ -- Emily Kirkmann (2007-02-11): added graph_border option to plot and show -- Robert L. Miller (2007-02-12): vertex color-maps, graph boundaries, - graph6 helper functions in SageX + graph6 helper functions in Cython SAGE Days 3 (2007-02-17--21): 3d plotting in Tachyon (2007-02-25): display a partition Index: sage/groups/perm_gps/permgroup_named.py =================================================================== --- sage/groups/perm_gps/permgroup_named.py (revision 5214) +++ sage/groups/perm_gps/permgroup_named.py (revision 5508) @@ -418,10 +418,11 @@ The projective general linear groups over GF(q). """ - def __init__(self, n, q): + def __init__(self, n, q, name='a'): """ INPUT: n -- positive integer; the degree q -- prime power; the size of the ground field - + name -- (default: 'a') variable name of indeterminate of finite field GF(q) + OUTPUT: PGL(n,q) @@ -437,4 +438,9 @@ 24 + + sage: G = PGL(2, 9, 'b'); G + Permutation Group with generators [(3,10,9,8,4,7,6,5), (1,2,4)(5,6,8)(7,9,10)] + sage: G.base_ring() + Finite Field in b of size 3^2 """ from sage.groups.perm_gps.permgroup import PermutationGroup, PermutationGroup_generic @@ -446,5 +452,5 @@ from_group=True, check=False) self._q = q - self._base_ring = GF(q) + self._base_ring = GF(q, name=name) self._n = n @@ -460,10 +466,11 @@ """ from sage.groups.perm_gps.permgroup import PermutationGroup, PermutationGroup_generic - def __init__(self, n, q): + def __init__(self, n, q, name='a'): """ INPUT: n -- positive integer; the degree q -- prime power; the size of the ground field - + name -- (default: 'a') variable name of indeterminate of finite field GF(q) + OUTPUT: PSL(n,q) @@ -479,4 +486,14 @@ The projective special linear group of degree 2 over Finite Field of size 3 + We create two groups over nontrivial finite fields: + sage: G = PSL(2, 4, 'b'); G + Permutation Group with generators [(3,4,5), (1,2,3)] + sage: G.base_ring() + Finite Field in b of size 2^2 + sage: G = PSL(2, 8); G + Permutation Group with generators [(3,8,6,4,9,7,5), (1,2,3)(4,7,5)(6,9,8)] + sage: G.base_ring() + Finite Field in a of size 2^3 + """ if n == 1: @@ -487,5 +504,5 @@ from_group=True, check=False) self._q = q - self._base_ring = GF(q) + self._base_ring = GF(q, name=name) self._n = n @@ -593,10 +610,11 @@ """ from sage.groups.perm_gps.permgroup import PermutationGroup, PermutationGroup_generic - def __init__(self, n, q): + def __init__(self, n, q, name='a'): """ INPUT: n -- positive integer; the degree q -- prime power; the size of the ground field - + name -- (default: 'a') variable name of indeterminate of finite field GF(q) + OUTPUT: PSp(n,q) @@ -618,4 +636,8 @@ Finite Field of size 3 + sage: G = PSp(2, 8, name='alpha'); G + Permutation Group with generators [(3,8,6,4,9,7,5), (1,2,3)(4,7,5)(6,9,8)] + sage: G.base_ring() + Finite Field in alpha of size 2^3 """ if n%2 == 1: @@ -626,5 +648,5 @@ from_group=True, check=False) self._q = q - self._base_ring = GF(q) + self._base_ring = GF(q, name=name) self._n = n @@ -644,5 +666,5 @@ n -- positive integer; the degree q -- prime power; the size of the ground field - + name -- (default: 'a') variable name of indeterminate of finite field GF(q) OUTPUT: PSU(n,q) @@ -651,13 +673,18 @@ sage: PSU(2,3) The projective special unitary group of degree 2 over Finite Field of size 3 - """ - from sage.groups.perm_gps.permgroup import PermutationGroup, PermutationGroup_generic - def __init__(self, n, q, var='a'): + + sage: G = PSU(2, 8, name='alpha'); G + The projective special unitary group of degree 2 over Finite Field in alpha of size 2^3 + sage: G.base_ring() + Finite Field in alpha of size 2^3 + """ + from sage.groups.perm_gps.permgroup import PermutationGroup, PermutationGroup_generic + def __init__(self, n, q, name='a'): id = 'PSU(%s,%s)'%(n,q) PermutationGroup_generic.__init__(self, id, from_group=True, check=False) self._q = q - self._base_ring = GF(q) - self._field_of_definition = GF(q**2, var) + self._base_ring = GF(q, name=name) + self._field_of_definition = GF(q**2, name) self._n = n @@ -678,5 +705,6 @@ n -- positive integer; the degree q -- prime power; the size of the ground field - + name -- (default: 'a') variable name of indeterminate of finite field GF(q) + OUTPUT: PGU(n,q) @@ -685,13 +713,18 @@ sage: PGU(2,3) The projective general unitary group of degree 2 over Finite Field of size 3 - """ - from sage.groups.perm_gps.permgroup import PermutationGroup, PermutationGroup_generic - def __init__(self, n, q, var='a'): + + sage: G = PGU(2, 8, name='alpha'); G + The projective general unitary group of degree 2 over Finite Field in alpha of size 2^3 + sage: G.base_ring() + Finite Field in alpha of size 2^3 + """ + from sage.groups.perm_gps.permgroup import PermutationGroup, PermutationGroup_generic + def __init__(self, n, q, name='a'): id = 'PGU(%s,%s)'%(n,q) PermutationGroup_generic.__init__(self, id, from_group=True, check=False) self._q = q - self._base_ring = GF(q) - self._field_of_definition = GF(q**2, var) + self._base_ring = GF(q, name=name) + self._field_of_definition = GF(q**2, name) self._n = n @@ -711,8 +744,12 @@ INPUT: - $q = 2^n$ -- an odd power of 2; the size of the ground field. (Strictly speaking, n should be - greater than 1, or else this group os not simple.) + q -- 2^n, an odd power of 2; the size of the ground + field. (Strictly speaking, n should be greater than 1, or + else this group os not simple.) + name -- (default: 'a') variable name of indeterminate of + finite field GF(q) + OUTPUT: - Sz(q) + A Suzuki group. EXAMPLE: @@ -726,4 +763,11 @@ The Suzuki group over Finite Field in a of size 2^3 + sage: G = SuzukiGroup(32, name='alpha') + sage: G.order() + 32537600 + sage: G.order().factor() + 2^10 * 5^2 * 31 * 41 + sage: G.base_ring() + Finite Field in alpha of size 2^5 REFERENCES: @@ -731,5 +775,5 @@ """ from sage.groups.perm_gps.permgroup import PermutationGroup, PermutationGroup_generic - def __init__(self, q, var='a'): + def __init__(self, q, name='a'): q = Integer(q) from sage.rings.arith import valuation @@ -741,5 +785,5 @@ from_group=True, check=False) self._q = q - self._base_ring = GF(q, var) + self._base_ring = GF(q, name=name) def base_ring(self): Index: sage/libs/ntl/ntl.pxi =================================================================== --- sage/libs/ntl/ntl.pxi (revision 3597) +++ sage/libs/ntl/ntl.pxi (revision 5511) @@ -1,5 +1,5 @@ +cdef extern object (make_GF2E(GF2E (*))) +cdef extern object (make_GF2EX(GF2EX (*))) +cdef extern object (make_GF2X(GF2X (*))) cdef extern object (make_ZZ(ntl_c_ZZ (*))) cdef extern object (make_ZZ_p(ZZ_p (*))) -cdef extern object (make_GF2X(GF2X (*))) -cdef extern object (make_GF2E(GF2E (*))) -cdef extern object (make_GF2EX(GF2EX (*))) Index: sage/matrix/docs.py =================================================================== --- sage/matrix/docs.py (revision 3521) +++ sage/matrix/docs.py (revision 5513) @@ -127,5 +127,5 @@ \begin{verbatim} -New matrices types can only be implemented in SageX. +New matrices types can only be implemented in Cython. *********** LEVEL 1 ********** Index: sage/matrix/matrix0.pyx =================================================================== --- sage/matrix/matrix0.pyx (revision 5493) +++ sage/matrix/matrix0.pyx (revision 5522) @@ -1306,4 +1306,19 @@ ################################################### + def is_symmetric(self): + """ + Returns True if this is a symmetric matrix. + """ + if self._ncols != self._nrows: return False + # could be bigger than an int on a 64-bit platform, this + # is the type used for indexing. + cdef Py_ssize_t i,j + + for i from 0 <= i < self._nrows: + for j from 0 <= j < i: + if self.get_unsafe(i,j) != self.get_unsafe(j,i): + return False + return True + def is_dense(self): """ Index: sage/misc/all.py =================================================================== --- sage/misc/all.py (revision 5263) +++ sage/misc/all.py (revision 5513) @@ -49,7 +49,8 @@ from sage_eval import sage_eval, sageobj -from sagex import sagex_lambda -from sagex_c import sagex -pyrex = sagex # synonym +from cython import cython_lambda +from cython_c import cython +pyrex = cython # synonym -- for now +sagex = cython # synonym -- for now from persist import save, load, dumps, loads, db, db_save @@ -103,4 +104,6 @@ norm, numerator, + numerical_approx, + N, objgens, objgen, Index: sage/misc/cython.py =================================================================== --- sage/misc/cython.py (revision 5513) +++ sage/misc/cython.py (revision 5513) @@ -0,0 +1,380 @@ +""" +Cython -- C-Extensions for Python + +AUTHORS: + -- William Stein (2006-01-18): initial version + -- William Stein (2007-07-28): update from sagex to cython +""" + +#***************************************************************************** +# Copyright (C) 2006 William Stein +# +# Distributed under the terms of the GNU General Public License (GPL) +# +# http://www.gnu.org/licenses/ +#***************************************************************************** + +import os, sys + +from misc import SPYX_TMP, SAGE_ROOT + +def cblas(): + if os.environ.has_key('SAGE_CBLAS'): + return os.environ['SAGE_CBLAS'] + elif os.path.exists('/usr/lib/libcblas.dylib') or \ + os.path.exists('/usr/lib/libcblas.so'): + return 'cblas' + elif os.path.exists('/usr/lib/libblas.dll.a'): # untested. + return 'blas' + else: + # This is very slow (?), but *guaranteed* to be available. + return 'gslcblas' + + +include_dirs = ['%s/local/include'%SAGE_ROOT, \ + '%s/local/include/python%s'%(SAGE_ROOT, sys.version[:3]), \ + '%s/devel/sage/sage/ext'%SAGE_ROOT, \ + '%s/devel/sage/'%SAGE_ROOT, \ + '%s/devel/sage/sage/gsl'%SAGE_ROOT] + + +standard_libs = ['mpfr', 'gmp', 'gmpxx', 'stdc++', 'pari', 'm', \ + 'mwrank', 'gsl', cblas(), 'ntl', 'csage'] + +offset = 0 + +def parse_keywords(kwd, s): + j = 0 + v = [] + while True: + i = s[j:].find(kwd) + if i == -1: break + j = i + j + s = s[:j] + '#' + s[j:] + j += len(kwd) + 1 + k = s[j:].find('\n') + if k == -1: + k = len(s) + for X in s[j:j+k].split(): + if X[0] == '#': # skip rest of line + break + v.append(X) + return v, s + +def environ_parse(s): + i = s.find('$') + if i == -1: + return s + j = s[i:].find('/') + if j == -1: + j = len(s) + else: + j = i + j + name = s[i+1:j] + if os.environ.has_key(name): + s = s[:i] + os.environ[name] + s[j:] + else: + return s + return environ_parse(s) + +def pyx_preparse(s): + lang = parse_keywords('clang', s) + if lang[0]: + lang = lang[0][0] + else: + lang = "c" + + v, s = parse_keywords('clib', s) + libs = v + standard_libs + + additional_source_files, s = parse_keywords('cfile', s) + + v, s = parse_keywords('cinclude', s) + inc = [environ_parse(x.replace('"','').replace("'","")) for x in v] + include_dirs + s = """ +include "cdefs.pxi" +""" + s + if lang != "c++": # has issues with init_csage() + s = """ +include "interrupt.pxi" # ctrl-c interrupt block support +include "stdsage.pxi" # ctrl-c interrupt block support +""" + s + return s, libs, inc, lang, additional_source_files + +################################################################ +# If the user attaches a .spyx file and changes it, we have +# to reload an .so. +# +# PROBLEM: Python does not allow one to reload an .so extension module. +# Solution, we create a different .so file and load that one, +# overwriting the definitions of everything in the original .so file. +# +# HOW: By using a sequence_number for each .spyx file; we keep +# these sequence numbers in a dict. +# +################################################################ + +sequence_number = {} + +def cython(filename, verbose=False, compile_message=False, + use_cache=False, create_local_c_file=False): + if filename[-5:] != '.spyx': + print "File (=%s) must have extension .spyx"%filename + + clean_filename = sanitize(filename) + base = os.path.split(os.path.splitext(clean_filename)[0])[1] + abs_base = os.path.abspath(os.path.split(os.path.splitext(clean_filename)[0])[0]) + + build_dir = '%s/%s'%(SPYX_TMP, base) + if os.path.exists(build_dir): + # There is already a module here. Maybe we do not have to rebuild? + # Find the name. + if use_cache: + prev_so = [F for F in os.listdir(build_dir) if F[-3:] == '.so'] + if len(prev_so) > 0: + prev_so = prev_so[0] # should have length 1 because of deletes below + if os.path.getmtime(filename) <= os.path.getmtime('%s/%s'%(build_dir, prev_so)): + # We do not have to rebuild. + return prev_so[:-3], build_dir + else: + os.makedirs(build_dir) + for F in os.listdir(build_dir): + G = '%s/%s'%(build_dir,F) + if not os.path.isdir(G): + os.unlink(G) + + cmd = 'cd "%s"; ln -s "%s"/* .'%(build_dir, abs_base) + os.system(cmd) + + if compile_message: + print "Compiling %s..."%filename + + F = open(filename).read() + + F, libs, includes, language, additional_source_files = pyx_preparse(F) + + # add the working directory to the includes so custom headers etc. work + includes.append(os.path.split(os.path.splitext(filename)[0])[0]) + + if language == 'c++': + extension = "cpp" + else: + extension = "c" + + global sequence_number + if not sequence_number.has_key(base): + sequence_number[base] = 0 + name = '%s_%s'%(base, sequence_number[base]) + + # increment the sequence number so will use a different one next time. + sequence_number[base] += 1 + + additional_source_files = ",".join(["'"+os.path.abspath(os.curdir)+"/"+fname+"'" \ + for fname in additional_source_files]) + + pyx = '%s/%s.pyx'%(build_dir, name) + open(pyx,'w').write(F) + setup=""" +# Build using 'python setup.py' +import distutils.sysconfig, os, sys +from distutils.core import setup, Extension + +if not os.environ.has_key('SAGE_ROOT'): + print " ERROR: The environment variable SAGE_ROOT must be defined." + sys.exit(1) +else: + SAGE_ROOT = os.environ['SAGE_ROOT'] + SAGE_LOCAL = SAGE_ROOT + '/local/' + +extra_link_args = ['-L' + SAGE_LOCAL + '/lib'] +extra_compile_args = ['-w'] + +ext_modules = [Extension('%s', sources=['%s.%s', %s], + libraries=%s, + library_dirs=[SAGE_LOCAL + '/lib/'], + extra_compile_args = extra_compile_args, + extra_link_args = extra_link_args, + language = '%s' )] + +setup(ext_modules = ext_modules, + include_dirs = %s) + """%(name, name, extension, additional_source_files, libs, language, includes) + open('%s/setup.py'%build_dir,'w').write(setup) + + cython_include = ' '.join(['-I %s'%x for x in includes if len(x.strip()) > 0 ]) + + cmd = 'cd %s && cython -p %s %s.pyx 1>log 2>err '%(build_dir, cython_include, name) + + if create_local_c_file: + target_c = '%s/_%s.c'%(os.path.abspath(os.curdir), base) + if language == 'c++': + target_c = target_c + "pp" + cmd += ' && cp %s.c %s'%(name, target_c) + + if verbose: + print cmd + if os.system(cmd): + log = open('%s/log'%build_dir).read() + err = subtract_from_line_numbers(open('%s/err'%build_dir).read(), offset) + raise RuntimeError, "Error converting %s to C:\n%s\n%s"%(filename, log, err) + + if language=='c++': + os.system("cd %s && mv %s.c %s.cpp"%(build_dir,name,name)) + +## if make_c_file_nice and os.path.exists(target_c): +## R = open(target_c).read() +## R = "/* THIS IS A PARSED TO MAKE READABLE VERSION OF THE C FILE. */" + R + +## # 1. Get rid of the annoying __pyx_'s before variable names. +## # R = R.replace('__pyx_v_','').replace('__pyx','') +## # 2. Replace the line number references by the actual code from the file, +## # since it is very painful to go back and forth, and the philosophy +## # of SAGE is that everything that can be very easy *is*. + +## pyx_file = os.path.abspath('%s/%s.pyx'%(build_dir,name)) +## S = '/* "%s":'%pyx_file +## n = len(S) +## last_i = -1 +## X = F.split('\n') +## stars = '*'*80 +## while True: +## i = R.find(S) +## if i == -1 or i == last_i: break +## last_i = i +## j = R[i:].find('\n') +## if j == -1: break +## line_number = int(R[i+n: i+j]) +## try: +## line = X[line_number-1] +## except IndexError: +## line = '(missing code)' +## R = R[:i+2] + '%s\n\n Line %s: %s\n\n%s'%(stars, line_number, line, stars) + R[i+j:] + +## open(target_c,'w').write(R) + + + cmd = 'cd %s && python setup.py build 1>log 2>err'%build_dir + if verbose: print cmd + if os.system(cmd): + log = open('%s/log'%build_dir).read() + err = open('%s/err'%build_dir).read() + raise RuntimeError, "Error compiling %s:\n%s\n%s"%(filename, log, err) + + # Move from lib directory. + cmd = 'mv %s/build/lib.*/* %s'%(build_dir, build_dir) + if verbose: print cmd + if os.system(cmd): + raise RuntimeError, "Error copying extension module for %s"%filename + + + return name, build_dir + + + +def subtract_from_line_numbers(s, n): + ans = [] + for X in s.split('\n'): + i = X.find(':') + j = i+1 + X[i+1:].find(':') + try: + ans.append('%s:%s:%s\n'%(X[:i], int(X[i+1:j]) - n, X[j+1:])) + except ValueError: + ans.append(X) + return '\n'.join(ans) + + +################################################################ +# COMPILE +################################################################ +def cython_lambda(vars, expr, + verbose=False, + compile_message=False, + use_cache=False): + """ + Create a compiled function which evaluates expr assuming machine + values for vars. + + WARNING: This implementation is not well tested. + + INPUT: + vars -- list of pairs (variable name, c-data type), where + the variable names and data types are strings. + OR -- a string such as + 'double x, int y, int z' + expr -- an expression involving the vars and constants; + You can access objects defined in the current + module scope globals() using sagobject_name. + See the examples below. + + EXAMPLES: + We create a Lambda function in pure Python (using the r to make sure the 3.2 + is viewed as a Python float): + sage: f = lambda x,y: x*x + y*y + x + y + 17r*x + 3.2r + + We make the same Lambda function, but in a compiled form. + sage: g = cython_lambda('double x, double y', 'x*x + y*y + x + y + 17*x + 3.2') + sage: g(2,3) + 55.200000000000003 + sage: g(0,0) + 3.2000000000000002 + + We access a global function and variable. + sage: a = 25 + sage: f = cython_lambda('double x', 'sage.math.sin(x) + sage.a') + sage: f(10) + 24.455978889110629 + sage: a = 50 + sage: f(10) + 49.455978889110632 + """ + if isinstance(vars, str): + v = vars + else: + v = ', '.join(['%s %s'%(typ,var) for typ, var in vars]) + + s = """ +class _s: + def __getattr__(self, x): + return globals()[x] + +sage = _s() + +def f(%s): + return %s + """%(v, expr) + if verbose: + print s + import sage.misc.misc + tmpfile = sage.misc.misc.tmp_filename() + ".spyx" + open(tmpfile,'w').write(s) + + import sage.server.support + d = {} + sage.server.support.cython_import_all(tmpfile, d, + verbose=verbose, compile_message=compile_message, + use_cache=use_cache, + create_local_c_file=False) + return d['f'] + + + +def sanitize(f): + """ + Given a filename f, replace it by a filename that is a valid Python + module name. + + This means that the characters are all alphanumeric or _'s and + doesn't begin with a numeral. + """ + s = '' + if f[0].isdigit(): + s += '_' + for a in f: + if a.isalnum(): + s += a + else: + s += '_' + return s + + + Index: sage/misc/cython_c.pyx =================================================================== --- sage/misc/cython_c.pyx (revision 5513) +++ sage/misc/cython_c.pyx (revision 5513) @@ -0,0 +1,37 @@ +import sage.misc.misc +import sage.server.support + +def cython(code, + verbose=False, compile_message=False, + make_c_file_nice=False, use_cache=False): + """ + Given a block of CYTHON (SAGE's variant of Pyrex) code (as a text + string), this function compiles it using your C compiler, and + includes it into the global scope of the module that called this + function. + + WARNING: Only use this from Python code, not from extension code, + since from extension code you would change the global scope (i.e., + of the SAGE interpreter). And it would be stupid, since you're + already writing Cython! + + Also, never use this in the standard SAGE library. Any code that + uses this can only run on a system that has a C compiler + installed, and we want to avoid making that assumption for casual + SAGE usage. Also, any code that uses this in the library would + greatly slow down startup time, since currently there is no + caching. + + AUTHOR: William Stein, 2006-10-31 + + TODO: Need to create a clever caching system so code only gets + compiled once. + """ + tmpfile = sage.misc.misc.tmp_filename() + ".spyx" + open(tmpfile,'w').write(code) + sage.server.support.cython_import_all(tmpfile, globals(), + verbose=verbose, compile_message=compile_message, + use_cache=use_cache, + create_local_c_file=False) + + Index: sage/misc/functional.py =================================================================== --- sage/misc/functional.py (revision 5369) +++ sage/misc/functional.py (revision 5507) @@ -30,7 +30,11 @@ import sage.interfaces.expect + from sage.rings.complex_double import CDF from sage.rings.real_double import RDF, RealDoubleElement + import sage.rings.real_mpfr +import sage.rings.complex_field +import sage.rings.integer import __builtin__ @@ -657,4 +661,47 @@ return x return x.numerator() + +def numerical_approx(x, prec=None, digits=None): + """ + Return a numerical approximation of x with at least prec bits of + precision. + + NOTE: N is an alias for numerical_approx. + + INPUT: + x -- an object that has a numerical_approx method, or can + be coerced into a real or complex field + prec -- an integer (bits of precision) + digits -- an integer (digits of precision) + + If neither the prec or digits are specified, the default + is 53 bits of precision. + + EXAMPLES: + sage: numerical_approx(pi, 10) + 3.1 + sage: numerical_approx(pi, digits=10) + 3.141592654 + sage: numerical_approx(pi^2 + e, digits=20) + 12.587886229548403854 + sage: N(pi^2 + e) + 12.5878862295484 + sage: N(pi^2 + e, digits=50) + 12.5878862295484038541947784712288136330709465009407 + """ + if prec is None: + if digits is None: + prec = 53 + else: + prec = int(digits * 3.4) + 2 + try: + return x.numerical_approx(prec) + except AttributeError: + try: + return sage.rings.real_mpfr.RealField(prec)(x) + except TypeError: + return sage.rings.complex_field.ComplexField(prec)(x) + +N = numerical_approx def objgens(x): Index: sage/misc/interpreter.py =================================================================== --- sage/misc/interpreter.py (revision 5137) +++ sage/misc/interpreter.py (revision 5513) @@ -109,5 +109,5 @@ from preparser import preparse_file -import sagex +import cython #import signal @@ -186,5 +186,5 @@ _ip.magic('run -i "%s"'%process_file(F)) elif F[-5:] == '.spyx': - X = load_sagex(F) + X = load_cython(F) __IPYTHON__.push(X) else: @@ -279,5 +279,5 @@ line = "" elif name[-5:] == '.spyx': - line = load_sagex(name) + line = load_cython(name) else: line = 'load("%s")'%name @@ -333,5 +333,5 @@ elif name[-5:] == '.spyx': try: - line = load_sagex(name) + line = load_cython(name) attached[name] = os.path.getmtime(name) except IOError, OSError: @@ -346,10 +346,10 @@ return line -def load_sagex(name): +def load_cython(name): cur = os.path.abspath(os.curdir) try: - mod, dir = sagex.sagex(name, compile_message=True, use_cache=True) + mod, dir = cython.cython(name, compile_message=True, use_cache=True) except (IOError, OSError, RuntimeError), msg: - print "Error compiling sagex file:\n%s"%msg + print "Error compiling cython file:\n%s"%msg return '' import sys Index: sage/misc/preparser.py =================================================================== --- sage/misc/preparser.py (revision 5137) +++ sage/misc/preparser.py (revision 5513) @@ -566,5 +566,5 @@ elif name_load[-5:] == '.spyx': import interpreter - L = interpreter.load_sagex(name_load) + L = interpreter.load_cython(name_load) else: #print "Loading of '%s' not implemented (load .py, .spyx, and .sage files)"%name_load Index: sage/misc/sageinspect.py =================================================================== --- sage/misc/sageinspect.py (revision 4399) +++ sage/misc/sageinspect.py (revision 5513) @@ -1,6 +1,6 @@ r"""nodoctest -Inspect Python, Sage, and Sagex objects. - -This module extends parts of Python's inspect module to Sagex objects. +Inspect Python, Sage, and Cython objects. + +This module extends parts of Python's inspect module to Cython objects. AUTHOR: @@ -13,7 +13,7 @@ sage: from sage.misc.sageinspect import * -Test introspection of modules defined in Python and Sagex files: - - Sagex modules: +Test introspection of modules defined in Python and Cython files: + + Cython modules: sage: sage_getfile(sage.rings.rational) '.../rational.pyx' @@ -30,12 +30,12 @@ sage: sage_getdoc(sage.misc.sageinspect).lstrip() - 'Inspect Python, Sage, and Sagex objects...' + 'Inspect Python, Sage, and Cython objects...' sage: sage_getsource(sage.misc.sageinspect).lstrip()[5:-1] - 'Inspect Python, Sage, and Sagex objects...' - -Test introspection of classes defined in Python and Sagex files: - - Sagex classes: + 'Inspect Python, Sage, and Cython objects...' + +Test introspection of classes defined in Python and Cython files: + + Cython classes: sage: sage_getfile(sage.rings.rational.Rational) '.../rational.pyx' @@ -57,7 +57,7 @@ 'class Attach:...' -Test introspection of functions defined in Python and Sagex files: - - Sagex functions: +Test introspection of functions defined in Python and Cython files: + + Cython functions: sage: sage_getdef(sage.rings.rational.make_rational, obj_name='mr') 'mr(s)' @@ -112,5 +112,5 @@ def _extract_embedded_position(docstring): r""" - If docstring has a Sagex embedded position, return a tuple (original_docstring, filename, line). If not, return None. + If docstring has a Cython embedded position, return a tuple (original_docstring, filename, line). If not, return None. AUTHOR: @@ -146,5 +146,5 @@ return inspect.getblock(lines[lineno:]) -def _sage_getargspec_sagex(source): +def _sage_getargspec_cython(source): r""" inspect.getargspec from source code. @@ -179,8 +179,8 @@ argname = s[0] - # Sagex often has type information; we split off the right most + # Cython often has type information; we split off the right most # identifier to discard this information argname = argname.split()[-1] - # Sagex often has C pointer symbols before variable names + # Cython often has C pointer symbols before variable names argname.lstrip('*') argnames.append(argname) @@ -198,5 +198,5 @@ return (argnames, None, None, argdefs) except: - raise ValueError, "Could not parse sagex argspec" + raise ValueError, "Could not parse cython argspec" def sage_getfile(obj): @@ -243,8 +243,8 @@ return sage_getargspec(obj.__call__) else: - # Perhaps it is binary and defined in a Sagex file + # Perhaps it is binary and defined in a Cython file source = sage_getsource(obj, is_binary=True) if source: - return _sage_getargspec_sagex(source) + return _sage_getargspec_cython(source) @@ -279,5 +279,5 @@ Return the docstring associated to obj as a string. - If obj is a Sagex object with an embedded position in its docstring, + If obj is a Cython object with an embedded position in its docstring, the embedded position is stripped. @@ -296,5 +296,5 @@ s = sage.misc.sagedoc.format(str(r)) - # If there is a Sagex embedded position, it needs to be stripped + # If there is a Cython embedded position, it needs to be stripped pos = _extract_embedded_position(s) if pos is not None: @@ -335,5 +335,5 @@ """ # If we can handle it, we do. This is because Python's inspect will - # happily dump binary for sagex extension source code. + # happily dump binary for cython extension source code. d = inspect.getdoc(obj) pos = _extract_embedded_position(d) @@ -370,5 +370,5 @@ sage: from sage.misc.sageinspect import * - sage: from sage.misc.sageinspect import _extract_source, _extract_embedded_position, _sage_getargspec_sagex, __internal_teststring + sage: from sage.misc.sageinspect import _extract_source, _extract_embedded_position, _sage_getargspec_cython, __internal_teststring If docstring is None, nothing bad happens: @@ -379,9 +379,9 @@ "...all..." - A sagex function with default arguments: + A cython function with default arguments: sage: sage_getdef(sage.rings.integer.Integer.factor, obj_name='factor') "factor(algorithm='pari')" - A sagex method without an embedded position can lead to surprising errors: + A cython method without an embedded position can lead to surprising errors: sage: sage_getsource(sage.rings.integer.Integer.__init__, is_binary=True) Traceback (most recent call last): @@ -413,10 +413,10 @@ pass # EOF # 14 - Test _sage_getargspec_sagex with multiple default arguments and a type: - sage: _sage_getargspec_sagex("def init(self, x=None, base=0):") + Test _sage_getargspec_cython with multiple default arguments and a type: + sage: _sage_getargspec_cython("def init(self, x=None, base=0):") (['self', 'x', 'base'], None, None, ('None', '0')) - sage: _sage_getargspec_sagex("def __init__(self, x=None, base=0):") + sage: _sage_getargspec_cython("def __init__(self, x=None, base=0):") (['self', 'x', 'base'], None, None, ('None', '0')) - sage: _sage_getargspec_sagex("def __init__(self, x=None, unsigned int base=0):") + sage: _sage_getargspec_cython("def __init__(self, x=None, unsigned int base=0):") (['self', 'x', 'base'], None, None, ('None', '0')) Index: age/misc/sagex.py =================================================================== --- sage/misc/sagex.py (revision 4802) +++ (revision ) @@ -1,376 +1,0 @@ -""" -SageX Compiled SAGE Code - -AUTHORS: - -- William Stein, 2006-01-18 - -""" - -#***************************************************************************** -# Copyright (C) 2006 William Stein -# -# Distributed under the terms of the GNU General Public License (GPL) -# -# http://www.gnu.org/licenses/ -#***************************************************************************** - -import os, sys - -from misc import SPYX_TMP, SAGE_ROOT - -def cblas(): - if os.environ.has_key('SAGE_CBLAS'): - return os.environ['SAGE_CBLAS'] - elif os.path.exists('/usr/lib/libcblas.dylib') or \ - os.path.exists('/usr/lib/libcblas.so'): - return 'cblas' - elif os.path.exists('/usr/lib/libblas.dll.a'): # untested. - return 'blas' - else: - # This is very slow (?), but *guaranteed* to be available. - return 'gslcblas' - - -include_dirs = ['%s/local/include'%SAGE_ROOT, \ - '%s/local/include/python%s'%(SAGE_ROOT, sys.version[:3]), \ - '%s/devel/sage/sage/ext'%SAGE_ROOT, \ - '%s/devel/sage/'%SAGE_ROOT, \ - '%s/devel/sage/sage/gsl'%SAGE_ROOT] - - -standard_libs = ['mpfr', 'gmp', 'gmpxx', 'stdc++', 'pari', 'm', \ - 'mwrank', 'gsl', cblas(), 'ntl', 'csage'] - -offset = 0 - -def parse_keywords(kwd, s): - j = 0 - v = [] - while True: - i = s[j:].find(kwd) - if i == -1: break - j = i + j - s = s[:j] + '#' + s[j:] - j += len(kwd) + 1 - k = s[j:].find('\n') - if k == -1: - k = len(s) - for X in s[j:j+k].split(): - if X[0] == '#': # skip rest of line - break - v.append(X) - return v, s - -def environ_parse(s): - i = s.find('$') - if i == -1: - return s - j = s[i:].find('/') - if j == -1: - j = len(s) - else: - j = i + j - name = s[i+1:j] - if os.environ.has_key(name): - s = s[:i] + os.environ[name] + s[j:] - else: - return s - return environ_parse(s) - -def pyx_preparse(s): - lang = parse_keywords('clang', s) - if lang[0]: - lang = lang[0][0] - else: - lang = "c" - - v, s = parse_keywords('clib', s) - libs = v + standard_libs - - additional_source_files, s = parse_keywords('cfile', s) - - v, s = parse_keywords('cinclude', s) - inc = [environ_parse(x.replace('"','').replace("'","")) for x in v] + include_dirs - s = """ -include "cdefs.pxi" -""" + s - if lang != "c++": # has issues with init_csage() - s = """ -include "interrupt.pxi" # ctrl-c interrupt block support -include "stdsage.pxi" # ctrl-c interrupt block support -""" + s - return s, libs, inc, lang, additional_source_files - -################################################################ -# If the user attaches a .spyx file and changes it, we have -# to reload an .so. -# -# PROBLEM: Python does not allow one to reload an .so extension module. -# Solution, we create a different .so file and load that one, -# overwriting the definitions of everything in the original .so file. -# -# HOW: By using a sequence_number for each .spyx file; we keep -# these sequence numbers in a dict. -# -################################################################ - -sequence_number = {} - -def sagex(filename, verbose=False, compile_message=False, - use_cache=False, create_local_c_file=False): - if filename[-5:] != '.spyx': - print "File (=%s) must have extension .spyx"%filename - - clean_filename = sanitize(filename) - base = os.path.split(os.path.splitext(clean_filename)[0])[1] - abs_base = os.path.abspath(os.path.split(os.path.splitext(clean_filename)[0])[0]) - - build_dir = '%s/%s'%(SPYX_TMP, base) - if os.path.exists(build_dir): - # There is already a module here. Maybe we do not have to rebuild? - # Find the name. - if use_cache: - prev_so = [F for F in os.listdir(build_dir) if F[-3:] == '.so'] - if len(prev_so) > 0: - prev_so = prev_so[0] # should have length 1 because of deletes below - if os.path.getmtime(filename) <= os.path.getmtime('%s/%s'%(build_dir, prev_so)): - # We do not have to rebuild. - return prev_so[:-3], build_dir - else: - os.makedirs(build_dir) - for F in os.listdir(build_dir): - G = '%s/%s'%(build_dir,F) - if not os.path.isdir(G): - os.unlink(G) - - cmd = 'cd "%s"; ln -s "%s"/* .'%(build_dir, abs_base) - os.system(cmd) - - if compile_message: - print "Compiling %s..."%filename - - F = open(filename).read() - - F, libs, includes, language, additional_source_files = pyx_preparse(F) - - # add the working directory to the includes so custom headers etc. work - includes.append(os.path.split(os.path.splitext(filename)[0])[0]) - - if language == 'c++': - extension = "cpp" - else: - extension = "c" - - global sequence_number - if not sequence_number.has_key(base): - sequence_number[base] = 0 - name = '%s_%s'%(base, sequence_number[base]) - - # increment the sequence number so will use a different one next time. - sequence_number[base] += 1 - - additional_source_files = ",".join(["'"+os.path.abspath(os.curdir)+"/"+fname+"'" \ - for fname in additional_source_files]) - - pyx = '%s/%s.pyx'%(build_dir, name) - open(pyx,'w').write(F) - setup=""" -# Build using 'python setup.py' -import distutils.sysconfig, os, sys -from distutils.core import setup, Extension - -if not os.environ.has_key('SAGE_ROOT'): - print " ERROR: The environment variable SAGE_ROOT must be defined." - sys.exit(1) -else: - SAGE_ROOT = os.environ['SAGE_ROOT'] - SAGE_LOCAL = SAGE_ROOT + '/local/' - -extra_link_args = ['-L' + SAGE_LOCAL + '/lib'] -extra_compile_args = ['-w'] - -ext_modules = [Extension('%s', sources=['%s.%s', %s], - libraries=%s, - library_dirs=[SAGE_LOCAL + '/lib/'], - extra_compile_args = extra_compile_args, - extra_link_args = extra_link_args, - language = '%s' )] - -setup(ext_modules = ext_modules, - include_dirs = %s) - """%(name, name, extension, additional_source_files, libs, language, includes) - open('%s/setup.py'%build_dir,'w').write(setup) - - sagex_include = ' '.join(['-I %s'%x for x in includes if len(x.strip()) > 0 ]) - - cmd = 'cd %s && sagexc -p %s %s.pyx 1>log 2>err '%(build_dir, sagex_include, name) - - if create_local_c_file: - target_c = '%s/_%s.c'%(os.path.abspath(os.curdir), base) - if language == 'c++': - target_c = target_c + "pp" - cmd += ' && cp %s.c %s'%(name, target_c) - - if verbose: - print cmd - if os.system(cmd): - log = open('%s/log'%build_dir).read() - err = subtract_from_line_numbers(open('%s/err'%build_dir).read(), offset) - raise RuntimeError, "Error converting %s to C:\n%s\n%s"%(filename, log, err) - - if language=='c++': - os.system("cd %s && mv %s.c %s.cpp"%(build_dir,name,name)) - -## if make_c_file_nice and os.path.exists(target_c): -## R = open(target_c).read() -## R = "/* THIS IS A PARSED TO MAKE READABLE VERSION OF THE C FILE. */" + R - -## # 1. Get rid of the annoying __pyx_'s before variable names. -## # R = R.replace('__pyx_v_','').replace('__pyx','') -## # 2. Replace the line number references by the actual code from the file, -## # since it is very painful to go back and forth, and the philosophy -## # of SAGE is that everything that can be very easy *is*. - -## pyx_file = os.path.abspath('%s/%s.pyx'%(build_dir,name)) -## S = '/* "%s":'%pyx_file -## n = len(S) -## last_i = -1 -## X = F.split('\n') -## stars = '*'*80 -## while True: -## i = R.find(S) -## if i == -1 or i == last_i: break -## last_i = i -## j = R[i:].find('\n') -## if j == -1: break -## line_number = int(R[i+n: i+j]) -## try: -## line = X[line_number-1] -## except IndexError: -## line = '(missing code)' -## R = R[:i+2] + '%s\n\n Line %s: %s\n\n%s'%(stars, line_number, line, stars) + R[i+j:] - -## open(target_c,'w').write(R) - - - cmd = 'cd %s && python setup.py build 1>log 2>err'%build_dir - if verbose: print cmd - if os.system(cmd): - log = open('%s/log'%build_dir).read() - err = open('%s/err'%build_dir).read() - raise RuntimeError, "Error compiling %s:\n%s\n%s"%(filename, log, err) - - # Move from lib directory. - cmd = 'mv %s/build/lib.*/* %s'%(build_dir, build_dir) - if verbose: print cmd - if os.system(cmd): - raise RuntimeError, "Error copying extension module for %s"%filename - - - return name, build_dir - - - -def subtract_from_line_numbers(s, n): - ans = [] - for X in s.split('\n'): - i = X.find(':') - j = i+1 + X[i+1:].find(':') - try: - ans.append('%s:%s:%s\n'%(X[:i], int(X[i+1:j]) - n, X[j+1:])) - except ValueError: - ans.append(X) - return '\n'.join(ans) - - -################################################################ -# COMPILE -################################################################ -def sagex_lambda(vars, expr, - verbose=False, - compile_message=False, - use_cache=False): - """ - Create a compiled function which evaluates expr assuming machine - values for vars. - - WARNING: This implementation is not well tested. - - INPUT: - vars -- list of pairs (variable name, c-data type), where - the variable names and data types are strings. - OR -- a string such as - 'double x, int y, int z' - expr -- an expression involving the vars and constants; - You can access objects defined in the current - module scope globals() using sagobject_name. - See the examples below. - - EXAMPLES: - We create a Lambda function in pure Python (using the r to make sure the 3.2 - is viewed as a Python float): - sage: f = lambda x,y: x*x + y*y + x + y + 17r*x + 3.2r - - We make the same Lambda function, but in a compiled form. - sage: g = sagex_lambda('double x, double y', 'x*x + y*y + x + y + 17*x + 3.2') - - We access a global function and variable. - sage: a = 25 - sage: f = sagex_lambda('double x', 'sage.math.sin(x) + sage.a') - sage: f(10) - 24.455978889110629 - sage: a = 50 - sage: f(10) - 49.455978889110632 - """ - if isinstance(vars, str): - v = vars - else: - v = ', '.join(['%s %s'%(typ,var) for typ, var in vars]) - - s = """ -class _s: - def __getattr__(self, x): - return globals()[x] - -sage = _s() - -def f(%s): - return %s - """%(v, expr) - if verbose: - print s - import sage.misc.misc - tmpfile = sage.misc.misc.tmp_filename() + ".spyx" - open(tmpfile,'w').write(s) - - import sage.server.support - d = {} - sage.server.support.sagex_import_all(tmpfile, d, - verbose=verbose, compile_message=compile_message, - use_cache=use_cache, - create_local_c_file=False) - return d['f'] - - - -def sanitize(f): - """ - Given a filename f, replace it by a filename that is a valid Python - module name. - - This means that the characters are all alphanumeric or _'s and - doesn't begin with a numeral. - """ - s = '' - if f[0].isdigit(): - s += '_' - for a in f: - if a.isalnum(): - s += a - else: - s += '_' - return s - - - Index: age/misc/sagex_c.pyx =================================================================== --- sage/misc/sagex_c.pyx (revision 1947) +++ (revision ) @@ -1,37 +1,0 @@ -import sage.misc.misc -import sage.server.support - -def sagex(code, - verbose=False, compile_message=False, - make_c_file_nice=False, use_cache=False): - """ - Given a block of SAGEX (SAGE's variant of Pyrex) code (as a text - string), this function compiles it using your C compiler, and - includes it into the global scope of the module that called this - function. - - WARNING: Only use this from Python code, not from extension code, - since from extension code you would change the global scope (i.e., - of the SAGE interpreter). And it would be stupid, since you're - already writing SageX! - - Also, never use this in the standard SAGE library. Any code that - uses this can only run on a system that has a C compiler - installed, and we want to avoid making that assumption for casual - SAGE usage. Also, any code that uses this in the library would - greatly slow down startup time, since currently there is no - caching. - - AUTHOR: William Stein, 2006-10-31 - - TODO: Need to create a clever caching system so code only gets - compiled once. - """ - tmpfile = sage.misc.misc.tmp_filename() + ".spyx" - open(tmpfile,'w').write(code) - sage.server.support.sagex_import_all(tmpfile, globals(), - verbose=verbose, compile_message=compile_message, - use_cache=use_cache, - create_local_c_file=False) - - Index: sage/modular/abvar/homology.py =================================================================== --- sage/modular/abvar/homology.py (revision 3544) +++ sage/modular/abvar/homology.py (revision 5510) @@ -130,8 +130,12 @@ return self.__abvar.dimension() * 2 - def submodule(self, U): + def submodule(self, U, check=True): r""" - Return the submodule of this homology group given by U, which should + Return the submodule of this homology group given by $U$, which should be a submodule of the free module associated to this homology group. + + INPUT: + U -- submodule of ambient free module + check -- currently ignored. NOTE: We do not check that U is invariant under all Hecke operators. Index: sage/modular/all.py =================================================================== --- sage/modular/all.py (revision 3512) +++ sage/modular/all.py (revision 5510) @@ -19,5 +19,6 @@ dimension_eis, dimension_modular_forms, - dimension_new_cusp_forms) + dimension_new_cusp_forms, + sturm_bound) from buzzard import buzzard_tpslopes Index: sage/modular/dims.py =================================================================== --- sage/modular/dims.py (revision 5145) +++ sage/modular/dims.py (revision 5510) @@ -2,12 +2,18 @@ Dimensions of spaces of modular forms -The dimension formulas and implementations in this module grew out of -a program that Bruce Caskel wrote (around 1996) in PARI, which Kevin -Buzzard subsequently extended. I (William Stein) then implemented it -in C++ for HECKE. I also implemented it in MAGMA. Also, the -functions for dimensions of spaces with nontrivial character are based -on a paper (that has no proofs) by Cohen and Oesterle (Springer -Lecture notes in math, volume 627, pages 69--78). The formulas for -GammaH(N) were found and implemented by Jordi Quer. +AUTHORS: + -- William Stein + -- Jordi Quer + +ACKNOWLEDGEMENT: + The dimension formulas and implementations in this module grew out + of a program that Bruce Caskel wrote (around 1996) in PARI, which + Kevin Buzzard subsequently extended. I (William Stein) then + implemented it in C++ for HECKE. I also implemented it in MAGMA. + Also, the functions for dimensions of spaces with nontrivial + character are based on a paper (that has no proofs) by Cohen and + Oesterle (Springer Lecture notes in math, volume 627, pages + 69--78). The formulas for GammaH(N) were found and implemented by + Jordi Quer. The formulas here are more complete than in HECKE or MAGMA. @@ -24,5 +30,5 @@ ########################################################################## - +import math from sage.rings.arith import (factor, euler_phi as phi, divisors, is_prime, @@ -847,2 +853,60 @@ return dimension_cusp_forms(X, k) + dimension_eis(X, k) +def sturm_bound(level, weight): + r""" + Returns the Sturm bound for modules forms with given level and + weight. + + EXAMPLES: + sage: sturm_bound(11,2) + 2 + sage: sturm_bound(389,2) + 65 + sage: sturm_bound(1,12) + 1 + sage: sturm_bound(100,2) + 30 + sage: sturm_bound(1,36) + 3 + + FURTHER DETAILS: Returns a positive integer~$n$ such that the + Hecke operators $T_1,\ldots, T_n$ acting on \emph{cusp forms} + generate the Hecke algebra as a $\Z$-module when the character is + trivial or quadratic. Otherwise, $T_1,\ldots,T_n$ generate the + Hecke algebra at least as a $\Z[\eps]$-module, where $\Z[\eps]$ is + the ring generated by the values of the Dirichlet character + $\eps$. Alternatively, this is a bound such that if two cusp + forms associated to this space of modular symbols are congruent + modulo $(\lambda, q^n)$, then they are congruent modulo $\lambda$. + + REFERENCES: + See the Agashe-Stein appendix to Lario and Schoof, \emph{Some + computations with Hecke rings and deformation rings}, + Experimental Math., 11 (2002), no. 2, 303-311. This result + originated in the paper Sturm, \emph{On the congruence of + modular forms}, Springer LNM 1240, 275--280, 1987. + + REMARK: + Kevin Buzzard pointed out to me (William Stein) in Fall 2002 + that the above bound is fine for $\Gamma_1(N)$ with character, + as one sees by taking a power of $f$. More precisely, if $f + \con 0 \pmod{p}$ for first $s$ coefficients, then $f^r \con 0 + \pmod{p}$ for first $sr$ coefficents. Since the weight of + $f^r$ is $r\cdot k(f)$, it follows that if $s \geq b$, where + $b$ is the Sturm bound for $\Gamma_0(N)$ at weight $k(f)$, then + $f^r$ has valuation large enough to be forced to be $0$ at + $r*k(f)$ by Sturm bound (which is valid if we choose $r$ + correctly). Thus $f \con 0 \pmod{p}$. Conclusion: For + $\Gamma_1(N)$ with fixed character, the Sturm bound is + \emph{exactly} the same as for $\Gamma_0(N)$. + + A key point is that we are finding $\Z[\eps]$ generators for the + Hecke algebra here, not $\Z$-generators. So if one wants + generators for the Hecke algebra over $\Z$, this bound must be + suitably modified (and I'm not sure what the modification is). + + AUTHOR: + -- William Stein + """ + return Integer(int(math.ceil((weight * idxG0(level) / ZZ(12))))) + Index: sage/modular/hecke/ambient_module.py =================================================================== --- sage/modular/hecke/ambient_module.py (revision 3569) +++ sage/modular/hecke/ambient_module.py (revision 5510) @@ -21,7 +21,5 @@ import degenmap - import module - import submodule @@ -30,7 +28,11 @@ import sage.rings.all +import sage.misc.misc as misc + import sage.rings.arith as arith import sage.matrix.matrix_space as matrix_space + +import sage.modular.dims as dims def is_AmbientHeckeModule(x): @@ -308,5 +310,6 @@ that we include the $n$ with $n$ not coprime to the level. """ - raise NotImplementedError + misc.verbose("WARNING: ambient.py -- hecke_bound; returning unproven guess.") + return dims.sturm_bound(self.level(), self.weight()) + 2*dims.dimension_eis(self.level(), self.weight()) + 5 def hecke_module_of_level(self, level): @@ -386,4 +389,6 @@ Returns True if and only if self is a submodule of V. """ + if not isinstance(V, module.HeckeModule_free_module): + raise TypeError, "V must be a Hecke module" if not V.is_ambient(): return False @@ -461,5 +466,5 @@ else: self.__is_new[p] = False - ns = self.submodule(d.kernel()) + ns = self.submodule(d.kernel(), check=False) ns.__is_new = {p:True} ns._is_full_hecke_module = True @@ -543,5 +548,5 @@ #end if #end for - os = self.submodule(d.image()) + os = self.submodule(d.image(), check=False) self.__is_old[p] = (os == self) @@ -558,10 +563,17 @@ if check: if not sage.modules.all.is_FreeModule(M): - raise TypeError, "M must be a free module" + V = self.free_module() + if isinstance(M, (list,tuple)): + M = V.span([V(x.element()) for x in M]) + else: + M = V.span(M) if not M.is_submodule(self.free_module()): raise TypeError, "M must be a submodule of the free module associated to this module." if M == self.free_module(): return self - return submodule.HeckeSubmodule(self, M, Mdual, check=False) + return self._submodule_class()(self, M, Mdual, check=check) + + def _submodule_class(self): + return submodule.HeckeSubmodule def submodule_from_nonembedded_module(self, V, Vdual=None, check=True): @@ -570,4 +582,7 @@ V -- submodule of ambient free module of the same rank as the rank of self. + Vdual -- used to pass in dual submodule + check -- whether to check that submodule is Hecke equivariant + OUTPUT: Hecke submodule of self Index: sage/modular/hecke/module.py =================================================================== --- sage/modular/hecke/module.py (revision 3550) +++ sage/modular/hecke/module.py (revision 5510) @@ -185,5 +185,15 @@ def is_zero(self): """ - Return True if this modular symbols space has dimension 0. + Return True if this Hecke module has dimension 0. + + EXAMPLES: + sage: ModularSymbols(11).is_zero() + False + sage: ModularSymbols(11).old_submodule().is_zero() + True + sage: CuspForms(10).is_zero() + True + sage: CuspForms(1,12).is_zero() + False """ return self.dimension() == 0 @@ -193,4 +203,11 @@ Return True if this space is invariant under all Hecke operators. + + Since self is guaranteed to be an anemic Hecke module, the + significance of this function is that it also ensures invariance + under Hecke operators of index that divide the level. + + EXAMPLES: + """ try: @@ -213,9 +230,28 @@ """ Return True if self is invariant under the Hecke operator $T_n$. + + Since self is guaranteed to be an anemic Hecke module it is + only interesting to call this function when $n$ is not coprime + to the level. + + EXAMPLES: + sage: M = ModularSymbols(22).cuspidal_subspace() + sage: M.is_hecke_invariant(2) + True + + We use check=False to create a nasty ``module'' that is not invariant under $T_2$: + sage: S = M.submodule(M.free_module().span([M.0.list()]), check=False); S + Modular Symbols subspace of dimension 1 of Modular Symbols space of dimension 7 for Gamma_0(22) of weight 2 with sign 0 over Rational Field + sage: S.is_hecke_invariant(2) + False + sage: [n for n in range(1,12) if S.is_hecke_invariant(n)] + [1, 3, 5, 7, 9, 11] """ if arith.gcd(n, self.level()) == 1: return True - try: - self.hecke_operator(n) + if self.is_ambient(): + return True + try: + self.hecke_operator(n).matrix() except ArithmeticError: return False @@ -318,11 +354,43 @@ return T.apply_sparse(self.gen(i)) - def _element_eigenvalue(self, x): + def _element_eigenvalue(self, x, name='alpha'): if not element.is_HeckeModuleElement(x): raise TypeError, "x must be a Hecke module element." if not x in self.ambient_hecke_module(): raise ArithmeticError, "x must be in the ambient Hecke module." - v = self.dual_eigenvector() + v = self.dual_eigenvector(name=name) return v.dot_product(x.element()) + + def _is_hecke_equivariant_free_module(self, submodule): + """ + Returns True if the given free submodule of the ambient free module + is invariant under all Hecke operators. + + EXAMPLES: + sage: M = ModularSymbols(11); V = M.free_module() + sage: M._is_hecke_equivariant_free_module(V.span([V.0])) + False + sage: M._is_hecke_equivariant_free_module(V) + True + sage: M._is_hecke_equivariant_free_module(M.cuspidal_submodule().free_module()) + True + + We do the same as above, but with a modular forms space: + sage: M = ModularForms(11); V = M.free_module() + sage: M._is_hecke_equivariant_free_module(V.span([V.0 + V.1])) + False + sage: M._is_hecke_equivariant_free_module(V) + True + sage: M._is_hecke_equivariant_free_module(M.cuspidal_submodule().free_module()) + True + """ + misc.verbose("Determining if free module is Hecke equivariant.") + bound = self.hecke_bound() + for p in arith.primes(bound+1): + try: + self.T(p).matrix().restrict(submodule, check=True) + except ArithmeticError: + return False + return True def _set_factor_number(self, i): @@ -500,5 +568,5 @@ W, is_irred = X[i] if is_irred: - A = self.submodule(W) + A = self.submodule(W, check=False) if compute_dual: if X[i][1] != is_irred: @@ -517,7 +585,7 @@ for i in range(len(U)): if compute_dual: - A = self.submodule(U[i], Udual[i]) + A = self.submodule(U[i], Udual[i], check=False) else: - A = self.submodule(U[i]) + A = self.submodule(U[i], check=False) D.append(A) for A in D: @@ -541,5 +609,5 @@ return self.free_module().degree() - def dual_eigenvector(self): + def dual_eigenvector(self, name='alpha'): """ Return an eigenvector for the Hecke operators acting on the @@ -550,4 +618,5 @@ INPUT: The input space must be simple. + name -- print name of generator for eigenvalue field. OUTPUT: @@ -555,4 +624,6 @@ ring. This vector is an eigenvector for all Hecke operators acting via their transpose. + + EXAMPLES: NOTES: @@ -567,10 +638,12 @@ modular symbols to a field. This functional phi is an eigenvector for the dual action of Hecke operators on - functionals. - """ - try: - return self.__dual_eigenvector - except AttributeError: + functionals. + """ + try: + return self.__dual_eigenvector[name] + except KeyError: pass + except AttributeError: + self.__dual_eigenvector = {} if not self.is_simple(): @@ -593,6 +666,6 @@ if n > 1: R = f.parent() - K = R.quotient(f, 'alpha') # Let K be the quotient R/(f), - # with generator printed "alpha". + K = R.quotient(f, name) # Let K be the quotient R/(f), + # with generator printed name alpha = K.gen() beta = ~alpha # multiplicative inverse of alpha @@ -632,6 +705,6 @@ w_lift = w_lift * (~alpha) - self.__dual_eigenvector = w_lift - return self.__dual_eigenvector + self.__dual_eigenvector[name] = w_lift + return w_lift def dual_hecke_matrix(self, n): @@ -650,8 +723,25 @@ return self._dual_hecke_matrices[n] - def eigenvalue(self, n): + def eigenvalue(self, n, name='alpha'): """ Assuming that self is a simple space, return the eigenvalue of the $n$th Hecke operator on self. + + INPUT: + n -- index of Hecke operator + name -- print representation of generator of eigenvalue field + + EXAMPLES: + sage: A = ModularSymbols(125,sign=1).new_subspace()[0] + sage: A.eigenvalue(7) + -3 + sage: A.eigenvalue(3) + -alpha - 2 + sage: A.eigenvalue(3,'w') + -w - 2 + sage: A.eigenvalue(3,'z').charpoly('x') + x^2 + 3*x + 1 + sage: A.hecke_polynomial(3) + x^2 + 3*x + 1 NOTES: @@ -672,5 +762,5 @@ n = int(n) try: - return self.__eigenvalues[n] + return self.__eigenvalues[n][name] except AttributeError: self.__eigenvalues = {} @@ -680,19 +770,10 @@ raise IndexError, "n must be a positive integer" - ## This was removed so that every element in - ## the return of system_of_eigenvalues had the same - ## parent. This was done below, so I just added the - ## n==1 case to the below if stmt. - ## -- Craig Citro, 07 Aug 06 - ## -## if n == 1: -## a1 = self.base_ring()(1) -## self.__eigenvalues[1] = a1 -## return a1 + ev = self.__eigenvalues if (arith.is_prime(n) or n==1): Tn_e = self._eigen_nonzero_element(n) - an = self._element_eigenvalue(Tn_e) - self.__eigenvalues[n] = an + an = self._element_eigenvalue(Tn_e, name=name) + dict_set(ev, n, name, an) return an @@ -707,5 +788,5 @@ (p, r) = (int(p), int(r)) pow = p**r - if not self.__eigenvalues.has_key(pow): + if not (ev.has_key(pow) and ev[pow].has_key(name)): # TODO: Optimization -- do something much more intelligent in case character is not defined. # For example, compute it using diamond operators @@ -713,18 +794,18 @@ if eps is None: Tn_e = self._eigen_nonzero_element(pow) - self.__eigenvalues[pow] = self._element_eigenvalue(Tn_e) + dict_set(ev, pow, name, self._element_eigenvalue(Tn_e)) else: # a_{p^r} := a_p * a_{p^{r-1}} - eps(p)p^{k-1} a_{p^{r-2}} - apr1 = self.eigenvalue(pow//p) - ap = self.eigenvalue(p) + apr1 = self.eigenvalue(pow//p, name=name) + ap = self.eigenvalue(p, name=name) k = self.weight() - apr2 = self.eigenvalue(pow//(p*p)) + apr2 = self.eigenvalue(pow//(p*p), name=name) apow = ap*apr1 - eps(p)*(p**(k-1)) * apr2 - self.__eigenvalues[pow] = apow - if prod == None: - prod = self.__eigenvalues[pow] + dict_set(ev, pow, name, apow) + if prod is None: + prod = ev[pow][name] else: - prod *= self.__eigenvalues[pow] - self.__eigenvalues[n] = prod + prod *= ev[pow][name] + dict_set(ev, n, name, prod) return prod @@ -888,9 +969,13 @@ return self.__projection - def system_of_eigenvalues(self, n): + def system_of_eigenvalues(self, n, name='alpha'): r""" Assuming that self is a simple space of modular symbols, return the eigenvalues $[a_1, \ldots, a_nmax]$ of the Hecke operators on self. See \code{self.eigenvalue(n)} for more details. + + INPUT: + n -- number of eigenvalues + alpha -- name of generate for eigenvalue field EXAMPLES: @@ -902,5 +987,4 @@ [[1, 1, 0], [1, -1, -alpha - 1]] - Next we define a function that does the above: sage: def b(N,k=2): @@ -921,14 +1005,21 @@ [1, alpha, -2*alpha - 1, -alpha - 1, 2*alpha, alpha - 2, 2*alpha + 2, -2*alpha - 1, 2, -2*alpha + 2] sage: v[0].parent() - Univariate Quotient Polynomial Ring in alpha over Rational Field with modulus x^2 + x - 1 - """ - return [self.eigenvalue(m) for m in range(1,n+1)] + Univariate Quotient Polynomial Ring in alpha over Rational Field with modulus x^2 + x - 1 + + This example illustrates setting the print name of the eigenvalue field. + sage: A = ModularSymbols(125,sign=1).new_subspace()[0] + sage: A.system_of_eigenvalues(10) + [1, alpha, -alpha - 2, -alpha - 1, 0, -alpha - 1, -3, -2*alpha - 1, 3*alpha + 2, 0] + sage: A.system_of_eigenvalues(10,'x') + [1, x, -x - 2, -x - 1, 0, -x - 1, -3, -2*x - 1, 3*x + 2, 0] + """ + return [self.eigenvalue(m, name=name) for m in range(1,n+1)] def weight(self): """ - Returns the weight of this modular symbols space. + Returns the weight of this Hecke module. INPUT: - ModularSymbols self -- an arbitrary space of modular symbols + self -- an arbitrary Hecke module OUTPUT: @@ -947,6 +1038,18 @@ """ Return the zero submodule of self. - """ - return self.submodule(self.free_module().zero_submodule()) + + EXAMPLES: + sage: ModularSymbols(11,4).zero_submodule() + Modular Symbols subspace of dimension 0 of Modular Symbols space of dimension 6 for Gamma_0(11) of weight 4 with sign 0 over Rational Field + sage: CuspForms(11,4).zero_submodule() + Modular Forms subspace of dimension 0 of Modular Forms space of dimension 4 for Congruence Subgroup Gamma0(11) of weight 4 over Rational Field + """ + return self.submodule(self.free_module().zero_submodule(), check=False) +def dict_set(v, n, key, val): + if v.has_key(n): + v[n][key] = val + else: + v[n] = {key:val} + Index: sage/modular/hecke/submodule.py =================================================================== --- sage/modular/hecke/submodule.py (revision 3515) +++ sage/modular/hecke/submodule.py (revision 5510) @@ -43,4 +43,9 @@ if not submodule.is_submodule(ambient.free_module()): raise ValueError, "submodule must be a submodule of the ambient free module" + + if check: + if not ambient._is_hecke_equivariant_free_module(submodule): + raise ValueError, "The submodule must be invariant under all Hecke operators." + self.__ambient = ambient self.__submodule = submodule @@ -68,5 +73,5 @@ # Neither is ambient M = self.free_module() + other.free_module() - return self.ambient().submodule(M) + return self.ambient().submodule(M, check=False) def __call__(self, x, check=True): @@ -190,5 +195,5 @@ if V.rank() + self.rank() == A.rank(): - C = A.submodule(V) + C = A.submodule(V, check=False) self.__complement = C return C @@ -323,5 +328,5 @@ V = self.free_module().intersection(other.free_module()) M = self.ambient_hecke_module() - return M.submodule(V) + return M.submodule(V, check=False) def is_ambient(self): @@ -458,9 +463,13 @@ Construct a submodule of self from the embedded free module M. """ - if not sage.modules.all.is_FreeModule(M): - raise TypeError, "M (=%s) must be a free module"%M - - if not M.is_submodule(self.free_module()): - raise TypeError, "M (=%s) must be a submodule of the free module (=%s) associated to this module."%(M, self.free_module()) + if check: + if not sage.modules.all.is_FreeModule(M): + V = self.ambient_module().free_module() + if isinstance(M, (list,tuple)): + M = V.span([V(x.element()) for x in M]) + else: + M = V.span(M) + if not M.is_submodule(self.free_module()): + raise TypeError, "M (=%s) must be a submodule of the free module (=%s) associated to this module."%(M, self.free_module()) return self.ambient().submodule(M, Mdual, check=check) @@ -471,4 +480,6 @@ V -- submodule of ambient free module of the same rank as the rank of self. + check -- whether to check that V is Hecke equivariant. + OUTPUT: Hecke submodule of self @@ -489,4 +500,4 @@ A = M_Vdual * M_E Vdual = A.row_space() - return self.ambient_hecke_module().submodule(V, Vdual) - + return self.ambient_hecke_module().submodule(V, Vdual, check=check) + Index: sage/modular/modform/ambient.py =================================================================== --- sage/modular/modform/ambient.py (revision 3616) +++ sage/modular/modform/ambient.py (revision 5510) @@ -91,4 +91,7 @@ return "Modular Forms space of dimension %s for %s of weight %s over %s"%( self.dimension(), self.group(), self.weight(), self.base_ring()) + + def _submodule_class(self): + return submodule.ModularFormsSubmodule def change_ring(self, base_ring): Index: sage/modular/modform/submodule.py =================================================================== --- sage/modular/modform/submodule.py (revision 4875) +++ sage/modular/modform/submodule.py (revision 5510) @@ -30,11 +30,14 @@ A submodule of an ambient space of modular forms. """ - def __init__(self, ambient_module, submodule): + def __init__(self, ambient_module, submodule, dual=None, check=False): """ ambient_module -- ModularFormsSpace - submodule -- a submodule of the ambient space. + submodule -- a submodule of the ambient space. + dual_module -- (default: None) ignored + check -- (default: False) whether to check that the + submodule is Hecke equivariant """ A = ambient_module - sage.modular.hecke.submodule.HeckeSubmodule.__init__(self, A, submodule) + sage.modular.hecke.submodule.HeckeSubmodule.__init__(self, A, submodule, check=check) space.ModularFormsSpace.__init__(self, A.group(), A.weight(), A.character(), A.base_ring()) @@ -54,4 +57,5 @@ +# TODO class ModularFormsSubmoduleWithBasis(ModularFormsSubmodule): pass Index: sage/modular/modsym/ambient.py =================================================================== --- sage/modular/modsym/ambient.py (revision 4877) +++ sage/modular/modsym/ambient.py (revision 5510) @@ -949,9 +949,4 @@ return self.factorization() - def hecke_bound(self): - # TODO - misc.verbose("WARNING: ambient.py -- hecke_bound; returning unproven guess.") - return 2*self.sturm_bound() + 10 - def is_cuspidal(self): try: @@ -1042,7 +1037,7 @@ if compute_dual: Vdual = S.transpose().kernel() - M = self.submodule(V, Vdual) + M = self.submodule(V, Vdual, check=False) else: - M = self.submodule(V) + M = self.submodule(V, check=False) M._set_sign(sign) return M @@ -1063,4 +1058,35 @@ def submodule(self, M, dual_free_module=None, check=True): + """ + Return the submdoule of M with given generators or free module. + + INPUT: + M -- a submodule of the ambient free module or generators for a submodule + dual_free_module (optional) -- this may be useful to speed up certain calculations; it + is the corresponding submodule of the ambient dual module + check -- bool (optional: default -- True); if True, actually check that M + is a module, which means it is invariant under all Hecke operators. + + EXAMPLES: + sage: M = ModularSymbols(11) + sage: M.submodule([M.0]) + Traceback (most recent call last): + ... + ValueError: The submodule must be invariant under all Hecke operators. + sage: M.eisenstein_submodule().basis() + ((1,0) - 1/5*(1,9),) + sage: M.basis() + ((1,0), (1,8), (1,9)) + sage: M.submodule([M.0 - 1/5*M.2]) + Modular Symbols subspace of dimension 1 of Modular Symbols space of dimension 3 for Gamma_0(11) of weight 2 with sign 0 over Rational Field + + NOTE: It would make more sense to only check that M is + invariant under the Hecke operators with index coprime to + the level. Unfortunately, I do not know a reasonable + algorithm for determining whether a module is invariant + under just the anemic Hecke algebra, since I do not know + an analogue of the Sturm bound for the anemic Hecke + algebra. -- William Stein, 2007-07-27 + """ if check: if not free_module.is_FreeModule(M): @@ -1070,8 +1096,6 @@ else: M = V.span(M) - elif not M.is_submodule(self.free_module()): - raise ArithmeticError, "M must be a submodule of the free module of self." return subspace.ModularSymbolsSubspace(self, M, dual_free_module=dual_free_module, check=check) - + ###################################################################### # Z-module of integral modular symbols. Index: sage/modular/modsym/space.py =================================================================== --- sage/modular/modsym/space.py (revision 4013) +++ sage/modular/modsym/space.py (revision 5510) @@ -1184,14 +1184,5 @@ Returns the Sturm bound for this space of modular symbols. - Returns a positive integer~$n$ such that the Hecke operators - $T_1,\ldots, T_n$ acting on \emph{cusp forms} generate the - Hecke algebra as a $\Z$-module when the character is trivial - or quadratic. Otherwise, $T_1,\ldots,T_n$ generate the Hecke - algebra at least as a $\Z[\eps]$-module, where $\Z[\eps]$ is - the ring generated by the values of the Dirichlet character - $\eps$. Alternatively, this is a bound such that if two cusp - forms associated to this space of modular symbols are - congruent modulo $(\lambda, q^n)$, then they are congruent - modulo $\lambda$. + Type \code{sturm\_bound?} for more details. EXAMPLES: @@ -1204,30 +1195,4 @@ sage: ModularSymbols(1,36).sturm_bound() 3 - - REFERENCES: - See the Agashe-Stein appendix to Lario and Schoof, \emph{Some - computations with Hecke rings and deformation rings}, - Experimental Math., 11 (2002), no. 2, 303-311. This result - originated in the paper Sturm, \emph{On the congruence of - modular forms}, Springer LNM 1240, 275--280, 1987. - - REMARK: - Kevin Buzzard pointed out to me (William Stein) in Fall 2002 - that the above bound is fine for $\Gamma_1(N)$ with character, - as one sees by taking a power of $f$. More precisely, if $f - \con 0 \pmod{p}$ for first $s$ coefficients, then $f^r \con 0 - \pmod{p}$ for first $sr$ coefficents. Since the weight of - $f^r$ is $r\cdot k(f)$, it follows that if $s \geq b$, where - $b$ is the Sturm bound for $\Gamma_0(N)$ at weight $k(f)$, then - $f^r$ has valuation large enough to be forced to be $0$ at - $r*k(f)$ by Sturm bound (which is valid if we choose $r$ - correctly). Thus $f \con 0 \pmod{p}$. Conclusion: For - $\Gamma_1(N)$ with fixed character, the Sturm bound is - \emph{exactly} the same as for $\Gamma_0(N)$. - - A key point is that we are finding $\Z[\eps]$ generators for - the Hecke algebra here, not $\Z$-generators. So if one wants - generators for the Hecke algebra over $\Z$, this bound must - be suitably modified. """ # For Gamma_0(N), n = \frac{k}{12}[\SL_2(\Z):\Gamma_0(N)] @@ -1235,7 +1200,7 @@ return self.__sturm_bound except: - self.__sturm_bound = \ - int(math.ceil((self.weight()*dims.idxG0(self.level()))/Integer(12))) + self.__sturm_bound = dims.sturm_bound(self.level(), self.weight()) return self.__sturm_bound + def plus_submodule(self, compute_dual=True): @@ -1272,7 +1237,7 @@ Y = self.dual_free_module() D = W.intersection(Y) - M = A.submodule(V, D) + M = A.submodule(V, D, check=False) else: - M = A.submodule(V) + M = A.submodule(V, check=False) M._set_sign(sign) return M @@ -1323,54 +1288,4 @@ """ raise NotImplementedError - - def hecke_bound(self): - r""" - Returns the Sturm bound for this space of modular symbols. - - Returns a positive integer~$n$ such that the Hecke operators - $T_1,\ldots, T_n$ acting on \emph{cusp forms} generate the - Hecke algebra as a $\Z$-module when the character is trivial - or quadratic. Otherwise, $T_1,\ldots,T_n$ generate the Hecke - algebra at least as a $\Z[\eps]$-module, where $\Z[\eps]$ is - the ring generated by the values of the Dirichlet character - $\eps$. Alternatively, this is a bound such that if two cusp - forms associated to this space of modular symbols are - congruent modulo $(\lambda, q^n)$, then they are congruent - modulo $\lambda$. - - REFERENCES: - See the Agashe-Stein appendix to Lario and Schoof's \emph{Some - computations with Hecke rings and deformation rings}, - Experimental Math., 11 (2002), no. 2, 303-311. This result - originated in the paper Sturm, \emph{On the congruence of - modular forms}, Springer LNM 1240, 275--280, 1987. - - REMARK: - - Kevin Buzzard pointed out to me (William Stein) in Fall 2002 - that the above bound is fine for $\Gamma_1(N)$ with character, - as one sees by taking a power of $f$. More precisely, if $f - \con 0 \pmod{p}$ for first $s$ coefficients, then $f^r \con 0 - \pmod{p}$ for first $sr$ coefficents. Since the weight of - $f^r$ is $r\cdot k(f)$, it follows that if $s \geq b$, where - $b$ is the Sturm bound for $\Gamma_0(N)$ at weight $k(f)$, then - $f^r$ has valuation large enough to be forced to be $0$ at - $r*k(f)$ by Sturm bound (which is valid if we choose $r$ - correctly). Thus $f \con 0 \pmod{p}$. Conclusion: For - $\Gamma_1(N)$ with fixed character, the Sturm bound is - \emph{exactly} the same as for $\Gamma_0(N)$. - - A key point is that we are finding $\Z[\eps]$ generators for - the Hecke algebra here, not $\Z$-generators. So if one wants - generators for the Hecke algebra over $\Z$, this bound must - be suitably modified. - """ - # For Gamma_0(N), n = \frac{k}{12}[\SL_2(\Z):\Gamma_0(N)] - try: - return self.__sturm_bound - except AttributeError: - self.__sturm_bound = \ - int(math.ceil((self.weight() + dims.idxG0(self.level()))/12.0)) - return self.__sturm_bound def modular_abelian_variety(self): Index: sage/modular/modsym/subspace.py =================================================================== --- sage/modular/modsym/subspace.py (revision 3564) +++ sage/modular/modsym/subspace.py (revision 5510) @@ -43,30 +43,24 @@ # Special Methods ################################ - def __init__(self, ambient_hecke_module, submodule, dual_free_module=None, check=True): - #self.__reduce = (ambient_hecke_module, submodule, dual_free_module, False) + def __init__(self, ambient_hecke_module, submodule, dual_free_module=None, check=False): + """ + INPUT: + ambient_hecke_module -- the ambient space of modular + symbols in which we're constructing a submodule + submodule -- the underlying free module of the submodule + dual_free_module -- underlying free module of the dual of + the submodule (optional) + check -- (default: False) whether to check that the + submodule is invariant under all Hecke operators T_p. + """ self.__ambient_hecke_module = ambient_hecke_module A = ambient_hecke_module sage.modular.modsym.space.ModularSymbolsSpace.__init__(self, A.group(), A.weight(), \ A.character(), A.sign(), A.base_ring()) - hecke.HeckeSubmodule.__init__(self, A, submodule, dual_free_module = dual_free_module) - - #def __reduce__(self): - # return self.__class__, self.__reduce, self.__dict__ + hecke.HeckeSubmodule.__init__(self, A, submodule, dual_free_module = dual_free_module, check=check) def _repr_(self): return "Modular Symbols subspace of dimension %s of %s"%( self.rank(), self.ambient_module()) - -## def __cmp__(self, other): -## if isinstance(other, ModularSymbolsSubspace): -## return cmp((self.__ambient_hecke_module, self.free_module()), -## (other.__ambient_hecke_module, other.free_module())) -## if not isinstance(other, sage.modular.modsym.space.ModularSymbolsSpace): -## return cmp(type(self), type(other)) -## c = cmp(self.ambient_hecke_module(), other.ambient_hecke_module()) -## if c: -## return c -## return cmp(self.free_module(), other.free_module()) - ################################ Index: sage/modules/matrix_morphism.py =================================================================== --- sage/modules/matrix_morphism.py (revision 2283) +++ sage/modules/matrix_morphism.py (revision 5510) @@ -181,8 +181,8 @@ E = self.__matrix.decomposition(is_diagonalizable=is_diagonalizable) if D.is_ambient(): - return Sequence([D.submodule(V) for V, _ in E], cr=True, check=False) + return Sequence([D.submodule(V, check=False) for V, _ in E], cr=True, check=False) else: B = D.basis_matrix() - return Sequence([D.submodule((V.basis_matrix() * B).row_space()) for V, _ in E], + return Sequence([D.submodule((V.basis_matrix() * B).row_space(), check=False) for V, _ in E], cr=True, check=False) @@ -210,5 +210,5 @@ B = V.basis_matrix() * D.basis_matrix() V = B.row_space() - return self.domain().submodule(V) + return self.domain().submodule(V, check=False) def image(self): @@ -221,5 +221,5 @@ B = V.basis_matrix() * D.basis_matrix() V = B.row_space() - return self.codomain().submodule(V) + return self.codomain().submodule(V, check=False) def matrix(self): Index: sage/rings/integer.pxi =================================================================== --- sage/rings/integer.pxi (revision 1319) +++ sage/rings/integer.pxi (revision 5511) @@ -1,2 +1,2 @@ +cdef extern mpz_t (*(get_value(Integer ))) cdef extern int (set_mpz(Integer ,mpz_t )) -cdef extern mpz_t (*(get_value(Integer ))) Index: sage/rings/integer_mod_ring.py =================================================================== --- sage/rings/integer_mod_ring.py (revision 4878) +++ sage/rings/integer_mod_ring.py (revision 5514) @@ -213,4 +213,7 @@ if cache: self._precompute_table() + + self._zero_element = integer_mod.IntegerMod(self, 0) + self._one_element = integer_mod.IntegerMod(self, 1) def krull_dimension(self): Index: sage/rings/number_field/number_field.py =================================================================== --- sage/rings/number_field/number_field.py (revision 5372) +++ sage/rings/number_field/number_field.py (revision 5506) @@ -474,4 +474,15 @@ OUTPUT: Integer if v is omitted, and Rational otherwise. + + EXAMPLES: + sage: K. = NumberField(x^3 + x^2 - 2*x + 8) + sage: K.disc() + -503 + sage: K.disc([1, t, t^2]) + -2012 + sage: K.disc([1/7, (1/5)*t, (1/3)*t^2]) + -2012/11025 + sage: (5*7*3)^2 + 11025 """ if v == None: @@ -479,8 +490,8 @@ return self.__disc except AttributeError: - self.__disc = ZZ(str(self.pari_nf()[2])) + self.__disc = QQ(str(self.pari_nf()[2])) return self.__disc else: - return Q(self.trace_pairing(v).det()) + return QQ(self.trace_pairing(v).det()) disc = discriminant Index: sage/rings/polynomial/multi_polynomial_libsingular.pyx =================================================================== --- sage/rings/polynomial/multi_polynomial_libsingular.pyx (revision 5421) +++ sage/rings/polynomial/multi_polynomial_libsingular.pyx (revision 5515) @@ -2396,6 +2396,7 @@ _p = pSubst(_p, mi, (parent._coerce_c(v))._poly) + gd = parent.gens_dict() for m,v in kw.iteritems(): - m = parent(m) + m = gd[m] for i from 0 < i <= _ring.N: if p_GetExp((m)._poly, i, _ring) != 0: Index: sage/rings/rational.pxi =================================================================== --- sage/rings/rational.pxi (revision 1319) +++ sage/rings/rational.pxi (revision 5511) @@ -1,3 +1,3 @@ +cdef extern void (set_from_Integer(Rational ,Integer )) +cdef extern void (set_from_Rational(Rational ,Rational )) cdef extern void (set_from_mpq(Rational ,mpq_t )) -cdef extern void (set_from_Rational(Rational ,Rational )) -cdef extern void (set_from_Integer(Rational ,Integer )) Index: sage/server/misc.py =================================================================== --- sage/server/misc.py (revision 5435) +++ sage/server/misc.py (revision 5523) @@ -46,10 +46,11 @@ -def open_page(address, port, secure, pause=6): +def open_page(address, port, secure): if secure: rsrc = 'https' else: rsrc = 'http' - cmd = 'echo `sleep %s; %s %s://%s:%s 1>&2 >/dev/null` &'%(pause, browser(), rsrc, address, port) + + cmd = 'echo `%s %s://%s:%s 1>&2 >/dev/null` &'%(browser(), rsrc, address, port) os.system(cmd) Index: sage/server/notebook/avatars.py =================================================================== --- sage/server/notebook/avatars.py (revision 5208) +++ sage/server/notebook/avatars.py (revision 5518) @@ -95,5 +95,5 @@ self._avatarId = avatarId if twist.OPEN_MODE: - rsrc = twist.AdminToplevel(self.cookie, 'root') + rsrc = twist.AdminToplevel(self.cookie, 'admin') return (iweb.IResource, rsrc, self.logout) Index: sage/server/notebook/cell.py =================================================================== --- sage/server/notebook/cell.py (revision 5287) +++ sage/server/notebook/cell.py (revision 5517) @@ -498,5 +498,5 @@ t = self.__in.rstrip() - if t.lstrip()[:5] == '%hide': + if t.lstrip().startswith('%hide'): cls = "cell_input_hide" else: @@ -518,14 +518,15 @@ r = max(1, number_of_rows(t.strip(), ncols)) - s += """ - - """%(cls, r, ncols, id, id, id, id,'readonly=1' if do_print else '', t) + if not do_print: + s += """ + + """%(cls, r, ncols, id, id, id, id,'readonly=1' if do_print else '', t) t = t.replace("<","<")+" " Index: sage/server/notebook/notebook.py =================================================================== --- sage/server/notebook/notebook.py (revision 5385) +++ sage/server/notebook/notebook.py (revision 5519) @@ -107,5 +107,5 @@ self.add_user('_sage_', '', '', account_type='user', force=True) self.add_user('guest', '', '', account_type='guest', force=True) - self.add_user('root', passwd, '', account_type='admin', force=True) + self.add_user('admin', passwd, '', account_type='admin', force=True) def user_exists(self, username): @@ -880,5 +880,5 @@ @@ -1509,4 +1509,5 @@ t = worksheet.plain_text(prompts=True, banner=False) + t = t.replace('<','<') body += """ 
%s
Index: sage/server/notebook/run_notebook.py
===================================================================
--- sage/server/notebook/run_notebook.py	(revision 5454)
+++ sage/server/notebook/run_notebook.py	(revision 5523)
@@ -21,10 +21,8 @@
 
 from sage.misc.misc import DOT_SAGE
-from   sage.server.misc import print_open_msg, find_next_available_port, open_page
+from   sage.server.misc import print_open_msg, find_next_available_port
 import os, shutil, socket
 
 import notebook
-
-PAUSE = 7
 
 conf_path       = os.path.join(DOT_SAGE, 'notebook')
@@ -63,7 +61,7 @@
 
     nb = notebook.load_notebook(directory)
-    if reset or not nb.user_exists('root'):
+    if reset or not nb.user_exists('admin'):
         while True:
-            print "Setting password for the root user."
+            print "Setting password for the admin user."
             passwd = getpass.getpass("Enter new password: ")
             passwd2 = getpass.getpass("Retype new password: ")
@@ -73,16 +71,15 @@
                 break
         if reset:
-            nb.user('root').set_password(passwd)
-            print "Password changed for user 'root'."
+            nb.user('admin').set_password(passwd)
+            print "Password changed for user 'admin'."
         else:
             nb.create_default_users(passwd)
-            print "User root created with the password you specified."
+            print "User admin created with the password you specified."
             print "\n\n"
             print "*"*70
             print "\n"
             if secure:
-                print "Login to the SAGE notebook as root with the password you specified above."
+                print "Login to the SAGE notebook as admin with the password you specified above."
             
-
     if not server_pool is None:
         nb.set_server_pool(server_pool)
@@ -116,4 +113,10 @@
         notebook_opts = '"%s",address="%s",port=%s,secure=%s' % (os.path.abspath(directory),
                 address, port, secure)
+
+        if open_viewer:
+            open_page = "from sage.server.misc import open_page; open_page('%s', %s, %s)"%(address, port, secure)
+        else:
+            open_page = ''
+        
         config = open(conf, 'w')
         config.write("""
@@ -163,6 +166,7 @@
 application = service.Application("SAGE Notebook")
 s = strports.service('%s', factory)
+%s
 s.setServiceParent(application)
-"""%(notebook_opts, not secure, os.path.abspath(directory), strport))
+"""%(notebook_opts, not secure, os.path.abspath(directory), strport, open_page))
 
 
@@ -172,5 +176,5 @@
         print_open_msg(address, port, secure=secure)
         if secure:
-            print "There is a root account.  If you do not remember the password,"
+            print "There is an admin account.  If you do not remember the password,"
             print "quit the notebook and type notebook(reset=True)."
         e = os.system('sage -twistd --pidfile="%s"/twistd.pd -ny "%s"/twistedconf.tac'%(directory, directory))
Index: sage/server/notebook/worksheet.py
===================================================================
--- sage/server/notebook/worksheet.py	(revision 5289)
+++ sage/server/notebook/worksheet.py	(revision 5518)
@@ -885,5 +885,5 @@
             s += cell.html(ncols, do_print=do_print) + '\n'
 
-        if not published:
+        if not do_print and not published:
             s += '\n\n'
             s += '\n
 
\n'
@@ -1870,8 +1870,8 @@
 
     ##########################################################
-    # Parsing the %sagex, %jsmath, %python, etc., extension.
-    ##########################################################
-
-    def sagex_import(self, cmd, C):
+    # Parsing the %cython, %jsmath, %python, etc., extension.
+    ##########################################################
+
+    def cython_import(self, cmd, C):
         # Choice: Can use either C.relative_id() or self.next_block_id().
         # C.relative_id() has the advantage that block evals are cached, i.e.,
@@ -1884,5 +1884,5 @@
         if not (os.path.exists(spyx) and open(spyx).read() == cmd):
             open(spyx,'w').write(cmd)
-        s  = '_support_.sagex_import_all("%s", globals())'%spyx
+        s  = '_support_.cython_import_all("%s", globals())'%spyx
         return s
 
@@ -1915,6 +1915,6 @@
                 return False, t
             s = t
-        if s.startswith("%pyrex") or s.startswith("%sagex"):  # a block of Sagex code.
-            return True, self.sagex_import(after_first_word(s).lstrip(), C)
+        if s.startswith("%cython") or s.startswith("%pyrex") or s.startswith("%sagex"):  # a block of Cython code.
+            return True, self.cython_import(after_first_word(s).lstrip(), C)
             
         i = s.find('\n')
Index: sage/server/support.py
===================================================================
--- sage/server/support.py	(revision 5135)
+++ sage/server/support.py	(revision 5513)
@@ -257,22 +257,22 @@
 
 ######################################################################
-# Sagex
-######################################################################
-import sage.misc.sagex
+# Cython
+######################################################################
+import sage.misc.cython
 import sys
 import __builtin__
 
-def sagex_import(filename, verbose=False, compile_message=False,
+def cython_import(filename, verbose=False, compile_message=False,
                  use_cache=False, create_local_c_file=True):
     """
     INPUT:
-        filename -- name of a file that contains sagex code
+        filename -- name of a file that contains cython code
         
     OUTPUT:
-        module -- the module that contains the compiled sagex code.
+        module -- the module that contains the compiled cython code.
 
     Raises an ImportError exception if anything goes wrong.
     """
-    name, build_dir = sage.misc.sagex.sagex(filename, verbose=verbose,
+    name, build_dir = sage.misc.cython.cython(filename, verbose=verbose,
                                             compile_message=compile_message,
                                             use_cache=use_cache,
@@ -282,17 +282,17 @@
 
 
-def sagex_import_all(filename, globals, verbose=False, compile_message=False,
+def cython_import_all(filename, globals, verbose=False, compile_message=False,
                      use_cache=False, create_local_c_file=True):
     """
     INPUT:
-        filename -- name of a file that contains sagex code
+        filename -- name of a file that contains cython code
         
     OUTPUT:
-        changes globals using the attributes of the Sagex module
+        changes globals using the attributes of the Cython module
         that do not begin with an underscore. 
 
     Raises an ImportError exception if anything goes wrong.
     """
-    m = sagex_import(filename, verbose=verbose, compile_message=compile_message,
+    m = cython_import(filename, verbose=verbose, compile_message=compile_message,
                      use_cache=use_cache,
                      create_local_c_file=create_local_c_file)
Index: sage/version.py
===================================================================
--- sage/version.py	(revision 5504)
+++ sage/version.py	(revision 5528)
@@ -1,2 +1,2 @@
 """nodoctests"""
-version='2.7.1'; date='2007-07-24'
+version='2.7.2'; date='2007-07-28'
Index: setup.py
===================================================================
--- setup.py	(revision 5486)
+++ setup.py	(revision 5520)
@@ -6,5 +6,5 @@
 
 
-## Choose cblas library -- note -- make sure to update sage/misc/sagex.py
+## Choose cblas library -- note -- make sure to update sage/misc/cython.py
 ## if you change this!! 
 if os.environ.has_key('SAGE_BLAS'):
@@ -446,6 +446,6 @@
               sources = ['sage/ext/interactive_constructors_c.pyx']), \
 
-    Extension('sage.misc.sagex_c',
-              sources = ['sage/misc/sagex_c.pyx']), \
+    Extension('sage.misc.cython_c',
+              sources = ['sage/misc/cython_c.pyx']), \
 
     Extension('sage.rings.real_mpfr',
@@ -598,4 +598,11 @@
               ), \
                             
+    Extension('sage.combinat.partitions',
+              ['sage/combinat/partitions.pyx',
+               'sage/combinat/partitions_c.cc'],
+              libraries = ['gmp', 'mpfr'],
+              language='c++'              
+              ), \
+
     Extension('sage.graphs.graph_fast',
               ['sage/graphs/graph_fast.pyx'],
@@ -631,6 +638,6 @@
 
 ######################################################################
-# CODE for generating C/C++ code from Pyrex and doing dependency
-# checking, etc.  In theory distutils would run Pyrex, but I don't
+# CODE for generating C/C++ code from Cython and doing dependency
+# checking, etc.  In theory distutils would run Cython, but I don't
 # trust it at all, and it won't have the more sophisticated dependency
 # checking that we need.
@@ -652,8 +659,8 @@
 
 
-def need_to_pyrex(filename, outfile):
+def need_to_cython(filename, outfile):
     """
     INPUT:
-        filename -- The name of a pyrex file in the SAGE source tree.
+        filename -- The name of a cython file in the SAGE source tree.
         outfile -- The name of the corresponding c or cpp file in the build directory.
 
@@ -688,5 +695,5 @@
         # The multiple imports with parens, e.g.,
         #        import (a.b.c.d, e.f.g)
-        # of Python are not allowed in Pyrex.
+        # of Python are not allowed in Cython.
         # In both cases, the module name is the second word if we split on whitespace.
         try:
@@ -740,39 +747,5 @@
     return False
 
-def process_pyrexembed_file(f, m):
-    # This is a pyrexembed file, so process accordingly.
-    dir, base = os.path.split(f[:-5])
-    tmp = '%s/.tmp_pyrexembed'%dir
-    if os.path.exists(tmp) and not os.path.isdir(tmp):
-        print "Please delete file '%s' in %s"%(tmp, dir)
-        sys.exit(1)
-    if not os.path.exists(tmp):
-        os.makedirs(tmp)
-    pyxe_file = "%s/%s.pyxe"%(tmp, base)
-
-    # The following files will be produced by pyrexembed.
-    cpp_file = "%s/%s_embed.cpp"%(dir, base)
-    c_file = "%s/%s.c"%(dir, base)
-    pyx_file = "%s/%s.pyx"%(dir,base)
-    pyx_embed_file = "%s/%s.pyx"%(tmp, base)
-    h_file = "%s/%s_embed.h"%(tmp, base)
-    if is_older(f, pyx_file) or is_older(f, cpp_file) or is_older(f, h_file):
-        os.system('cp -p %s %s'%(f, pyxe_file))
-        os.system('cp -p %s/*.pxi %s'%(dir, tmp))
-        os.system('cp -p %s/*.pxd %s'%(dir, tmp))
-        os.system('cp -p %s/*.h %s'%(dir, tmp))
-        cmd = "pyrexembed %s"%pyxe_file
-        print cmd
-        ret = os.system(cmd)
-        if ret != 0:
-            print "sage: Error running pyrexembed."
-            sys.exit(1)
-        process_pyrex_file(pyx_embed_file, m)
-        cmd = 'cp -p %s/*.pyx %s/; cp -p %s/*.c %s/; cp -p %s/*.h %s/; cp -p %s/*.cpp %s/'%(tmp, dir, tmp, dir, tmp, dir, tmp, dir)
-        print cmd
-        os.system(cmd)
-    return [cpp_file, c_file]
-
-def process_pyrex_file(f, m):
+def process_cython_file(f, m):
     """
     INPUT:
@@ -780,5 +753,5 @@
         m -- Extension module description (i.e., object of type Extension).
     """
-    # This is a pyrex file, so process accordingly.
+    # This is a cython file, so process accordingly.
     pyx_inst_file = '%s/%s'%(SITE_PACKAGES, f)
     if is_older(f, pyx_inst_file):
@@ -789,14 +762,14 @@
         outfile += 'pp'
 
-    if need_to_pyrex(f, outfile):
-        cmd = "pyrexc --embed-positions -I%s %s"%(os.getcwd(), f)
+    if need_to_cython(f, outfile):
+        cmd = "cython --embed-positions -I%s %s"%(os.getcwd(), f)
         print cmd
         ret = os.system(cmd)
         if ret != 0:
-            print "sage: Error running pyrexc."
+            print "sage: Error running cython."
             sys.exit(1)
         # If the language for the extension is C++,
         # then move the resulting output file to have the correct extension.
-        # (I don't know how to tell Pyrex to do this automatically.)
+        # (I don't know how to tell Cython to do this automatically.)
         if m.language == 'c++':
             os.system('mv %s.c %s'%(f[:-4], outfile))
@@ -804,5 +777,5 @@
     
 
-def pyrex(ext_modules):
+def cython(ext_modules):
     for m in ext_modules:
         m.extra_compile_args += extra_compile_args
@@ -812,17 +785,14 @@
             f = m.sources[i]
             s = open(f).read()
-            if f[-5:] == '.pyxe':# and s.find("#embed") != -1 and s.find('#}embed') != -1:
-                new_sources = process_pyrexembed_file(f, m)
-            elif f[-4:] == ".pyx":
-                new_sources += process_pyrex_file(f, m)
+            if f[-4:] == ".pyx":
+                new_sources += process_cython_file(f, m)
             else:
                 new_sources.append(f)
         m.sources = new_sources
 
-
         
 
 if not sdist:
-    pyrex(ext_modules)
+    cython(ext_modules)
 
 setup(name        = 'sage',
@@ -838,5 +808,5 @@
       author_email= 'wstein@gmail.com',
       
-      url         = 'http://modular.math.washington.edu/sage',
+      url         = 'http://www.sagemath.org',
       
       packages    = ['sage',