Hide keyboard shortcuts

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

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

############################################################# 

# 

# Sparse Vector over mpq_t (the GMP rationals) 

# 

############################################################# 

  

from cysignals.memory cimport sig_malloc, sig_free 

  

from sage.libs.gmp.mpq cimport * 

from sage.data_structures.binary_search cimport * 

  

from sage.rings.rational cimport Rational 

  

cdef int reallocate_mpq_vector(mpq_vector* v, Py_ssize_t num_nonzero) except -1: 

mpq_vector_clear(v) 

allocate_mpq_vector(v, num_nonzero) 

return 0 

  

cdef int allocate_mpq_vector(mpq_vector* v, Py_ssize_t num_nonzero) except -1: 

""" 

Allocate memory for a mpq_vector -- most user code won't call this. 

num_nonzero is the number of nonzero entries to allocate memory for. 

  

This function *does* call mpq_init on each mpq_t that is allocated. 

It does *not* clear the entries of v, if there are any. 

""" 

cdef Py_ssize_t i 

v.entries = <mpq_t *> sig_malloc(num_nonzero*sizeof(mpq_t)) 

if v.entries == NULL: 

raise MemoryError("Error allocating memory") 

for i from 0 <= i < num_nonzero: 

mpq_init(v.entries[i]) 

v.positions = <Py_ssize_t*>sig_malloc(num_nonzero*sizeof(Py_ssize_t)) 

if v.positions == NULL: 

for i from 0 <= i < num_nonzero: 

mpq_clear(v.entries[i]) 

sig_free(v.entries) 

v.entries = NULL 

raise MemoryError("Error allocating memory") 

return 0 

  

cdef int mpq_vector_init(mpq_vector* v, Py_ssize_t degree, Py_ssize_t num_nonzero) except -1: 

""" 

Initialize a mpq_vector -- most user code *will* call this. 

""" 

allocate_mpq_vector(v, num_nonzero) 

v.num_nonzero = num_nonzero 

v.degree = degree 

  

cdef void mpq_vector_clear(mpq_vector* v): 

cdef Py_ssize_t i 

if v.entries == NULL: 

return 

# Free all mpq objects allocated in creating v 

for i from 0 <= i < v.num_nonzero: 

mpq_clear(v.entries[i]) 

# Free entries and positions of those entries. 

# These were allocated from the Python heap. 

# If mpq_vector_init was not called, then this 

# will (of course!) cause a core dump. 

sig_free(v.entries) 

sig_free(v.positions) 

  

cdef Py_ssize_t mpq_binary_search0(mpq_t* v, Py_ssize_t n, mpq_t x): 

""" 

Find the position of the rational 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 mpq_equal(v[i],x): 

return i 

return -1 

k = (i+j)/2 

c = mpq_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 mpq_binary_search(mpq_t* v, Py_ssize_t n, mpq_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 mpq_t (rational) 

n -- integer (length of array v) 

x -- mpq_t (rational) 

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 = mpq_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 = mpq_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 mpq_vector_get_entry(mpq_t ans, mpq_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 mpq_t 

that *must* have been initialized using mpq_init. 

""" 

if n >= v.degree: 

raise IndexError("Index must be between 0 and %s." % (v.degree - 1)) 

cdef Py_ssize_t m 

m = binary_search0(v.positions, v.num_nonzero, n) 

if m == -1: 

mpq_set_si(ans, 0,1) 

return 0 

mpq_set(ans, v.entries[m]) 

return 0 

  

cdef object mpq_vector_to_list(mpq_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 Rational a 

cdef Py_ssize_t i 

X = [] 

for i from 0 <= i < v.num_nonzero: 

a = Rational() 

a.set_from_mpq(v.entries[i]) 

X.append( (v.positions[i], a) ) 

return X 

  

  

cdef int mpq_vector_set_entry(mpq_vector* v, Py_ssize_t n, mpq_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. 

""" 

if n >= v.degree or n < 0: 

raise IndexError("Index must be between 0 and the degree minus 1.") 

cdef Py_ssize_t i, m, ins 

cdef Py_ssize_t m2, ins2 

cdef Py_ssize_t *pos 

cdef mpq_t *e 

  

m = binary_search(v.positions, v.num_nonzero, n, &ins) 

if m != -1: 

# 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. 

if mpq_sgn(x) != 0: # x != 0, case 1 

# v.entries[m] = x 

mpq_set(v.entries[m], x) 

else: # case 2 

e = v.entries 

pos = v.positions 

allocate_mpq_vector(v, v.num_nonzero - 1) 

for i from 0 <= i < m: 

# v.entries[i] = e[i] 

mpq_set(v.entries[i], e[i]) 

v.positions[i] = pos[i] 

mpq_clear(e[i]) 

mpq_clear(e[m]) 

for i from m < i < v.num_nonzero: 

# v.entries[i-1] = e[i] 

mpq_set(v.entries[i-1], e[i]) 

mpq_clear(e[i]) 

v.positions[i-1] = pos[i] 

sig_free(e) 

sig_free(pos) 

v.num_nonzero = v.num_nonzero - 1 

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. 

if mpq_sgn(x) == 0: 

return 0 

v.num_nonzero = v.num_nonzero + 1 

e = v.entries 

pos = v.positions 

allocate_mpq_vector(v, v.num_nonzero) 

for i from 0 <= i < ins: 

# v.entries[i] = e[i] 

mpq_set(v.entries[i], e[i]) 

mpq_clear(e[i]) 

v.positions[i] = pos[i] 

# v.entries[ins] = x 

mpq_set(v.entries[ins], x) 

v.positions[ins] = n 

for i from ins < i < v.num_nonzero: 

mpq_set(v.entries[i], e[i-1]) 

mpq_clear(e[i-1]) 

v.positions[i] = pos[i-1] 

sig_free(e) 

sig_free(pos) 

  

  

  

cdef mpq_t mpq_set_tmp 

mpq_init(mpq_set_tmp) 

cdef int mpq_vector_set_entry_str(mpq_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. 

""" 

mpq_set_str(mpq_set_tmp, x_str, 0) 

mpq_vector_set_entry(v, n, mpq_set_tmp) 

  

  

cdef int add_mpq_vector_init(mpq_vector* sum, 

mpq_vector* v, 

mpq_vector* w, 

mpq_t multiple) except -1: 

""" 

Initialize sum and set sum = v + multiple*w. 

""" 

if v.degree != w.degree: 

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 mpq_vector* z 

cdef mpq_t tmp 

if mpq_cmp_si(multiple, 0, 1) == 0: 

mpq_vector_init(sum, v.degree, 0) 

return 0 

  

mpq_init(tmp) 

# Do not do the multiply if the multiple is 1. 

do_multiply = mpq_cmp_si(multiple, 1,1) 

  

z = sum 

# 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: 

nz = v.num_nonzero + w.num_nonzero 

if nz > v.degree: nz = v.degree 

mpq_vector_init(z, v.degree, nz) 

# 2. Merge entries 

i = 0 # index into entries of v 

j = 0 # index into entries of w 

k = 0 # index into z (the vector we are creating) 

while i < v.num_nonzero or j < w.num_nonzero: 

if i >= v.num_nonzero: # just copy w in 

z.positions[k] = w.positions[j] 

  

if do_multiply: 

# This means: z.entries[k] = (multiple*w.entries[j]) 

mpq_mul(z.entries[k], multiple, w.entries[j]) 

else: 

mpq_set(z.entries[k], w.entries[j]) 

j = j + 1 

k = k + 1 

elif j >= w.num_nonzero: # just copy v in 

z.positions[k] = v.positions[i] 

# This means: z.entries[k] = v.entries[i] 

mpq_set(z.entries[k], v.entries[i]) 

i = i + 1 

k = k + 1 

elif v.positions[i] < w.positions[j]: # copy entry from v in 

z.positions[k] = v.positions[i] 

# This means: z.entries[k] = v.entries[i] 

mpq_set(z.entries[k], v.entries[i]) 

i = i + 1 

k = k + 1 

elif v.positions[i] > w.positions[j]: # copy entry from w in 

if do_multiply: 

# This means: tmp = multiple*w.entries[j] 

mpq_mul(tmp, multiple, w.entries[j]) 

# This means: z.entries[k] = tmp 

mpq_set(z.entries[k], tmp) 

else: 

mpq_set(z.entries[k], w.entries[j]) 

z.positions[k] = w.positions[j] 

k = k + 1 

j = j + 1 

else: # equal, so add and copy 

if do_multiply: 

# This means: tmp = v.entries[i] + multiple*w.entries[j] 

mpq_mul(tmp, multiple, w.entries[j]) 

mpq_add(tmp, tmp, v.entries[i]) 

else: 

mpq_add(tmp, v.entries[i], w.entries[j]) 

if mpq_sgn(tmp) != 0: 

z.positions[k] = v.positions[i] 

# This means: z.entries[k] = tmp 

mpq_set(z.entries[k], tmp) 

k = k + 1 # only increment if sum is nonzero! 

i = i + 1 

j = j + 1 

#end if 

# end while 

z.num_nonzero = k 

for i from k <= i < z.num_nonzero: 

mpq_clear(z.entries[i]) 

mpq_clear(tmp) 

return 0 

  

cdef int mpq_vector_scale(mpq_vector* v, mpq_t scalar) except -1: 

if mpq_sgn(scalar) == 0: # scalar = 0 

mpq_vector_clear(v) 

mpq_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] 

mpq_mul(v.entries[i], v.entries[i], scalar) 

return 0 

  

cdef int mpq_vector_scalar_multiply(mpq_vector* v, mpq_vector* w, mpq_t scalar) except -1: 

""" 

v = w * scalar 

""" 

cdef Py_ssize_t i 

if v == w: 

# rescale self 

return mpq_vector_scale(v, scalar) 

else: 

mpq_vector_clear(v) 

v.entries = <mpq_t*> sig_malloc(w.num_nonzero * sizeof(mpq_t)) 

if v.entries == NULL: 

v.positions = NULL 

raise MemoryError("error allocating rational sparse vector mpq's") 

v.positions = <Py_ssize_t*> sig_malloc(w.num_nonzero * sizeof(Py_ssize_t)) 

if v.positions == NULL: 

sig_free(v.entries) 

v.entries = NULL 

raise MemoryError("error allocating rational sparse vector positions") 

v.num_nonzero = w.num_nonzero 

v.degree = w.degree 

for i from 0 <= i < v.num_nonzero: 

mpq_init(v.entries[i]) 

mpq_mul(v.entries[i], w.entries[i], scalar) 

v.positions[i] = w.positions[i] 

return 0 

  

cdef int mpq_vector_cmp(mpq_vector* v, mpq_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 = mpq_cmp(v.entries[i], w.entries[i]) 

if c < 0: 

return -1 

elif c > 0: 

return 1 

return 0