Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
############################################################# # # Sparse Vector over mpz_t (the GMP integers) # #############################################################
from cysignals.memory cimport sig_malloc, sig_free
from sage.libs.gmp.mpz cimport * from sage.data_structures.binary_search cimport * from sage.rings.integer cimport Integer
cdef int allocate_mpz_vector(mpz_vector* v, Py_ssize_t num_nonzero) except -1: """ Allocate memory for a mpz_vector -- most user code won't call this. num_nonzero is the number of nonzero entries to allocate memory for.
This function *does* call mpz_init on each mpz_t that is allocated. It does *not* clear the entries of v, if there are any. """ cdef Py_ssize_t i raise MemoryError("Error allocating memory") for i from 0 <= i < num_nonzero: mpz_clear(v.entries[i]) sig_free(v.entries) v.entries = NULL raise MemoryError("Error allocating memory")
cdef int mpz_vector_init(mpz_vector* v, Py_ssize_t degree, Py_ssize_t num_nonzero) except -1: """ Initialize a mpz_vector -- most user code *will* call this. """
cdef void mpz_vector_clear(mpz_vector* v): cdef Py_ssize_t i # Free all mpz objects allocated in creating v # Free entries and positions of those entries. # These were allocated from the Python heap. # If mpz_vector_init was not called, then this # will (of course!) cause a core dump.
cdef Py_ssize_t mpz_binary_search0(mpz_t* v, Py_ssize_t n, mpz_t x): """ Find the position of the integers x in the array v, which has length n. Returns -1 if x is not in the array v. """ cdef Py_ssize_t i, j, k, c if n == 0: return -1 i = 0 j = n-1 while i<=j: if i == j: if mpz_cmp(v[i],x) == 0: return i return -1 k = (i+j)/2 c = mpz_cmp(v[k],x) if c > 0: # v[k] > x j = k-1 elif c < 0: # v[k] < x i = k+1 else: # only possibility is that v[k] == x return k return -1
cdef Py_ssize_t mpz_binary_search(mpz_t* v, Py_ssize_t n, mpz_t x, Py_ssize_t* ins): """ Find the position of the integer x in the array v, which has length n. Returns -1 if x is not in the array v, and in this case ins is set equal to the position where x should be inserted in order to obtain an ordered array.
INPUT: v -- array of mpz_t (integer) n -- integer (length of array v) x -- mpz_t (integer) OUTPUT: position of x (as an Py_ssize_t) ins -- (call be pointer), the insertion point if x is not found. """ cdef Py_ssize_t i, j, k, c if n == 0: ins[0] = 0 return -1 i = 0 j = n-1 while i<=j: if i == j: c = mpz_cmp(v[i],x) if c == 0: # v[i] == x ins[0] = i return i if c < 0: # v[i] < x ins[0] = i + 1 else: ins[0] = i return -1 k = (i+j)/2 c = mpz_cmp(v[k], x) if c > 0: # v[k] > x j = k-1 elif c < 0: # v[k] < x i = k+1 else: # only possibility is that v[k] == x ins[0] = k return k # end while ins[0] = j+1 return -1
cdef int mpz_vector_get_entry(mpz_t ans, mpz_vector* v, Py_ssize_t n) except -1: """ Returns the n-th entry of the sparse vector v. This would be v[n] in Python syntax.
The return is done using the pointer ans, which is to an mpz_t that *must* have been initialized using mpz_init. """ raise IndexError("Index (=%s) must be between 0 and %s." % (n, v.degree - 1)) cdef Py_ssize_t m
cdef object mpz_vector_to_list(mpz_vector* v): """ Returns a Python list of 2-tuples (i,x), where x=v[i] runs through the nonzero elements of x, in order. """ cdef object X cdef Integer a cdef Py_ssize_t i X = [] for i from 0 <= i < v.num_nonzero: a = Integer() a.set_from_mpz(v.entries[i]) X.append( (v.positions[i], a) ) return X
cdef int mpz_vector_set_entry(mpz_vector* v, Py_ssize_t n, mpz_t x) except -1: """ Set the n-th component of the sparse vector v equal to x. This would be v[n] = x in Python syntax. """ raise IndexError("Index (=%s) must be between 0 and %s." % (n, v.degree - 1)) cdef Py_ssize_t i, m, ins cdef Py_ssize_t m2, ins2 cdef Py_ssize_t *pos cdef mpz_t *e
# The position n was found in the array of positions. # Now there are two cases: # 1. x =/= 0, which is easy, and # 2. x = 0, in which case we have to recopy # positions and entries, without the m-th # element, and change num_nonzero. # v.entries[m] = x else: # case 2 # v.entries[i] = e[i] # v.entries[i-1] = e[i] else: # Allocate new memory and copy over elements from the # old array. This is similar to case 2 above, # except we are inserting a new entry rather than # deleting an old one. The new entry should be inserted # at position ins, which was computed using binary search. # # There is one exception -- if the new entry is 0, we # do nothing and return. # v.entries[i] = e[i] # v.entries[ins] = x
cdef mpz_t mpz_set_tmp mpz_init(mpz_set_tmp) cdef int mpz_vector_set_entry_str(mpz_vector* v, Py_ssize_t n, char *x_str) except -1: """ Set the n-th component of the sparse vector v equal to x. This would be v[n] = x in Python syntax. """ mpz_set_str(mpz_set_tmp, x_str, 0) mpz_vector_set_entry(v, n, mpz_set_tmp)
cdef int add_mpz_vector_init(mpz_vector* sum, mpz_vector* v, mpz_vector* w, mpz_t multiple) except -1: """ Initialize sum and set sum = v + multiple*w. """ print("Can't add vectors of degree %s and %s" % (v.degree, w.degree)) raise ArithmeticError("The vectors must have the same degree.")
cdef Py_ssize_t nz, i, j, k, do_multiply cdef mpz_vector* z cdef mpz_t tmp mpz_vector_init(sum, v.degree, 0) return 0
# Do not do the multiply if the multiple is 1.
# ALGORITHM: # 1. Allocate enough memory to hold the union of the two # lists of positions. We allocate the sum of the number # of positions of both (up to the degree), to avoid # having to make two passes. This might be slightly wasteful of # memory, but is faster. # 2. Move along the entries of v and w, copying them into the # new position / entry array. When position are the same, # add modulo p. # 3. Set num_nonzero and return success code.
# 1. Allocate memory: # 2. Merge entries
# This means: z.entries[k] = (multiple*w.entries[j]) else: # This means: z.entries[k] = v.entries[i] # This means: z.entries[k] = v.entries[i] # This means: tmp = multiple*w.entries[j] # This means: z.entries[k] = tmp else: else: # equal, so add and copy # This means: tmp = v.entries[i] + multiple*w.entries[j] else: # This means: z.entries[k] = tmp #end if # end while
cdef int mpz_vector_scale(mpz_vector* v, mpz_t scalar) except -1: if mpz_sgn(scalar) == 0: # scalar = 0 mpz_vector_clear(v) mpz_vector_init(v, v.degree, 0) return 0 cdef Py_ssize_t i for i from 0 <= i < v.num_nonzero: # v.entries[i] = scalar * v.entries[i] mpz_mul(v.entries[i], v.entries[i], scalar) return 0
cdef int mpz_vector_scalar_multiply(mpz_vector* v, mpz_vector* w, mpz_t scalar) except -1: """ v = w * scalar """ cdef Py_ssize_t i # rescale self return mpz_vector_scale(v, scalar) else: v.positions = NULL raise MemoryError("error allocating rational sparse vector mpz's") sig_free(v.entries) v.entries = NULL raise MemoryError("error allocating rational sparse vector positions")
cdef int mpz_vector_cmp(mpz_vector* v, mpz_vector* w): if v.degree < w.degree: return -1 elif v.degree > w.degree: return 1 cdef Py_ssize_t i cdef int c for i from 0 <= i < v.num_nonzero: c = mpz_cmp(v.entries[i], w.entries[i]) if c < 0: return -1 elif c > 0: return 1 return 0
|