# Ticket #54: freegb.lib

File freegb.lib, 51.4 KB (added by , 14 years ago) |
---|

Line | |
---|---|

1 | ////////////////////////////////////////////////////////////////////////////// |

2 | version="$Id: freegb.lib,v 1.12 2008/10/06 17:04:27 Singular Exp $"; |

3 | category="Noncommutative"; |

4 | info=" |

5 | LIBRARY: freegb.lib Twosided Noncommutative Groebner bases in Free Algebras |

6 | AUTHOR: Viktor Levandovskyy, levandov@math.rwth-aachen.de |

7 | |

8 | BACKGROUND: |

9 | TODO: WRITE DOWN SOME THEORY BEHIND (DIVISIBILITY, ORDERING ON FREE ASSOC. ALG., GB) |

10 | TODO: ADD REFERENCES |

11 | TODO: EXPLAIN YOUR ASSUMPTIONS/IMPLICITE AGREEMENTS/INPUT, OUTPUT FORMATS ETC. |

12 | TODO: BAD NAMES, SEE 3.9.1 Procedures in a library, 5^th rule. PLEASE RENAME! |

13 | |

14 | |

15 | MAIN PROCEDURES: |

16 | |

17 | freegbRing(int d); creates a ring with d blocks of shifted original variables |

18 | freegbasis(list L, int n); computes two-sided GB of an ideal up to degree n |

19 | |

20 | AUXILIARY PROCEDURES: |

21 | |

22 | lp2lstr(ideal K, def save); converts letter-place ideal to a list of modules |

23 | lst2str(list L[,int n]); converts a list of modules into polys in free alg. |

24 | mod2str(module M[,int n]); converts a module into a polynomial in free algebra |

25 | vct2str(module M[,int n]); converts a vector into a word in free algebra |

26 | " |

27 | |

28 | // TODO: LOTS OF GLOBAL UNLISTED PROCS WITHOUT HELP/EXAMPLES... -> |

29 | // -> MAKE THEM STATIC OR DOCUMENT THEM! |

30 | // TODO: AT LEAST SOME (E.G. Serre MUST BE LISTED IN THE HEADER) |

31 | |

32 | // TODO: MOVE THE FOLLOWING TO THE ABOVE BACKGROUND... |

33 | |

34 | // this library computes two-sided GB of an ideal |

35 | // in a free associative algebra |

36 | |

37 | // a monomial is encoded via a vector V |

38 | // where V[1] = coefficient |

39 | // V[1+i] = the corresponding symbol |

40 | |

41 | |

42 | |

43 | |

44 | |

45 | LIB "discretize.lib"; // for replace |

46 | LIB "qhmoduli.lib"; // for Max |

47 | |

48 | |

49 | // obsolete? |

50 | |

51 | proc lshift(module M, int s, string varing, def lpring) |

52 | { |

53 | // FINALLY IMPLEMENTED AS A PART OT THE CODE |

54 | // shifts a poly from the ring @R to s positions |

55 | // M lives in varing, the result in lpring |

56 | // to be run from varing |

57 | int i, j, k, sm, sv; |

58 | vector v; |

59 | // execute("setring "+lpring); |

60 | setring lpring; |

61 | poly @@p; |

62 | ideal I; |

63 | execute("setring "+varing); |

64 | sm = ncols(M); |

65 | for (i=1; i<=s; i++) |

66 | { |

67 | // modules, e.g. free polynomials |

68 | for (j=1; j<=sm; j++) |

69 | { |

70 | //vectors, e.g. free monomials |

71 | v = M[j]; |

72 | sv = size(v); |

73 | sp = "@@p = @@p + "; |

74 | for (k=2; k<=sv; k++) |

75 | { |

76 | sp = sp + string(v[k])+"("+string(k-1+s)+")*"; |

77 | } |

78 | sp = sp + string(v[1])+";"; // coef; |

79 | setring lpring; |

80 | // execute("setring "+lpring); |

81 | execute(sp); |

82 | execute("setring "+varing); |

83 | } |

84 | setring lpring; |

85 | // execute("setring "+lpring); |

86 | I = I,@@p; |

87 | @@p = 0; |

88 | } |

89 | setring lpring; |

90 | //execute("setring "+lpring); |

91 | export(I); |

92 | // setring varing; |

93 | execute("setring "+varing); |

94 | } |

95 | |

96 | proc skip0(vector v) |

97 | { |

98 | // skips zeros in a vector, producing another vector |

99 | int sv = nrows(v); |

100 | int sw = size(v); |

101 | if (sv == sw) |

102 | { |

103 | return(v); |

104 | } |

105 | int i; |

106 | int j=1; |

107 | vector w; |

108 | for (i=1; i<=sv; i++) |

109 | { |

110 | if (v[i] != 0) |

111 | { |

112 | w = w + v[i]*gen(j); |

113 | j++; |

114 | } |

115 | } |

116 | return(w); |

117 | } |

118 | |

119 | |

120 | // TODO: ANY ASSUMPTIONS ON THE BASERING? |

121 | // TODO: EXPLAIN WHICH FREE ALGEBRA IS MEANT? |

122 | |

123 | proc lst2str(list L, list #) |

124 | "USAGE: lst2str(L[,n]); L a list of modules, n an optional integer |

125 | RETURN: list of strings |

126 | PURPOSE: converts a list of modules into a list of strings representing |

127 | corresponding module as an element in the free algebra |

128 | EXAMPLE: example lst2str; shows examples |

129 | NOTE: if an optional integer is not 0, stars signs are used in multiplication |

130 | " |

131 | { |

132 | // returns a list of strings |

133 | // being sentences in words built from L |

134 | // if #[1] = 1, use * between generators |

135 | int useStar = 0; |

136 | if ( size(#)>0 ) |

137 | { |

138 | if (typeof(#[1]) == "int") |

139 | { |

140 | useStar = 1; |

141 | } |

142 | } |

143 | int i; |

144 | int s = size(L); |

145 | list N; |

146 | for(i=1; i<=s; i++) |

147 | { |

148 | if ((typeof(L[i]) == "module") || (typeof(L[i]) == "matrix") ) |

149 | { |

150 | N[i] = mod2str(L[i],useStar); |

151 | } |

152 | else |

153 | { |

154 | ERROR("module or matrix expected in the list"); |

155 | brake; |

156 | // return(N); |

157 | } |

158 | } |

159 | return(N); |

160 | } |

161 | example |

162 | { |

163 | "EXAMPLE:"; echo = 2; |

164 | ring r = 0,(x,y,z),(dp(1),dp(2)); |

165 | module M = [-1,x,y],[-7,y,y],[3,x,x]; |

166 | module N = [1,x,y,x,y],[-2,y,x,y,x],[6,x,y,y,x,y]; |

167 | list L; L[1] = M; L[2] = N; |

168 | lst2str(L); |

169 | lst2str(L[1],1); |

170 | } |

171 | |

172 | |

173 | // TODO: ANY ASSUMPTIONS ON THE BASERING/INPUT!? |

174 | // TODO: EXPLAIN WHICH FREE ALGEBRA IS MEANT? |

175 | |

176 | proc mod2str(module M, list #) |

177 | "USAGE: mod2str(M[,n]); M a module, n an optional integer |

178 | RETURN: string |

179 | PURPOSE: Converts the module M into a string representing M |

180 | as an element in the free algebra |

181 | EXAMPLE: example mod2str; shows examples |

182 | NOTE: if an optional integer is not 0, stars signs are used in multiplication |

183 | " |

184 | { |

185 | // returns a string |

186 | // a sentence in words built from M |

187 | // if #[1] = 1, use * between generators |

188 | int useStar = 0; |

189 | if ( size(#)>0 ) |

190 | { |

191 | if (typeof(#[1]) == "int") |

192 | { |

193 | useStar = 1; |

194 | } |

195 | } |

196 | int i; |

197 | int s = ncols(M); |

198 | string t; |

199 | string mp; |

200 | for(i=1; i<=s; i++) |

201 | { |

202 | mp = vct2str(M[i],useStar); |

203 | if (mp[1] == "-") |

204 | { |

205 | t = t + mp; |

206 | } |

207 | else |

208 | { |

209 | t = t + "+" + mp; |

210 | } |

211 | } |

212 | if (t[1]=="+") |

213 | { |

214 | t = t[2..size(t)]; // remove first "+" |

215 | } |

216 | return(t); |

217 | } |

218 | example |

219 | { |

220 | "EXAMPLE:"; echo = 2; |

221 | ring r = 0,(x,y,z),(dp); |

222 | module M = [1,x,y,x,y],[-2,y,x,y,x],[6,x,y,y,x,y]; |

223 | mod2str(M); |

224 | mod2str(M,1); |

225 | // TODO: WHAT ABOUT ZERO? |

226 | M = 0; |

227 | mod2str(M); // TODO: FIX BUG: MUST BE JUST A ZERO! |

228 | M = [0], [0], [1], [0, x]; |

229 | mod2str(M); // TODO: FIX THIS BUGS! |

230 | M = [1, x, y, z], [0, x], [2, y], [0]; |

231 | mod2str(M); |

232 | } |

233 | |

234 | // TODO: ANY ASSUMPTIONS ON THE BASERING/INPUT!? |

235 | // TODO: EXPLAIN WHICH FREE ALGEBRA IS MEANT? |

236 | |

237 | proc vct2str(vector v, list #) |

238 | "USAGE: vct2str(v[,n]); v a vector, n an optional integer |

239 | RETURN: string |

240 | PURPOSE: Converts a vector into a string representing a word in the free alg. |

241 | EXAMPLE: example vct2str; shows examples |

242 | NOTE: if an optional integer is not 0, stars signs are used in multiplication |

243 | " |

244 | { |

245 | // if #[1] = 1, use * between generators |

246 | int useStar = 0; |

247 | if ( size(#)>0 ) |

248 | { |

249 | if (typeof(#[1]) == "int") |

250 | { |

251 | useStar = 1; |

252 | } |

253 | } |

254 | int ppl = printlevel-voice+2; |

255 | // for a word, encoded by v |

256 | // produces a string for it |

257 | v = skip0(v); |

258 | number cf = leadcoef(v[1]); |

259 | int s = size(v); |

260 | string vs,vv,vp,err; |

261 | int i,j,p,q; |

262 | for (i=1; i<=s-1; i++) |

263 | { |

264 | p = IsVar(v[i+1]); |

265 | if (p==0) |

266 | { |

267 | err = "Error: monomial expected at" + string(i+1); |

268 | dbprint(ppl,err); |

269 | return("_"); |

270 | } |

271 | if (p==1) |

272 | { |

273 | if (useStar && (size(vs) >0)) { vs = vs + "*"; } |

274 | vs = vs + string(v[i+1]); |

275 | } |

276 | else //power |

277 | { |

278 | vv = string(v[i+1]); |

279 | q = find(vv,"^"); |

280 | if (q==0) |

281 | { |

282 | q = find(vv,string(p)); |

283 | if (q==0) |

284 | { |

285 | err = "error in find for string "+vv; |

286 | dbprint(ppl,err); |

287 | return("_"); |

288 | } |

289 | } |

290 | // q>0 |

291 | vp = vv[1..q-1]; |

292 | for(j=1;j<=p;j++) |

293 | { |

294 | if (useStar && (size(vs) >0)) { vs = vs + "*"; } |

295 | vs = vs + vp; |

296 | } |

297 | } |

298 | } |

299 | string scf; |

300 | if (cf == -1) |

301 | { |

302 | scf = "-"; |

303 | } |

304 | else |

305 | { |

306 | scf = string(cf); |

307 | if (cf == 1) |

308 | { |

309 | scf = ""; |

310 | } |

311 | } |

312 | if (useStar && (size(scf) >0) && (scf!="-") ) { scf = scf + "*"; } |

313 | vs = scf + vs; |

314 | return(vs); |

315 | } |

316 | example |

317 | { |

318 | "EXAMPLE:"; echo = 2; |

319 | ring r = (0,a),(x,y3,z(1)),dp; |

320 | vector v = [-7,x,y3^4,x2,z(1)^3]; |

321 | vct2str(v); |

322 | vct2str(v,1); |

323 | vector w = [-7a^5+6a,x,y3,y3,x,z(1),z(1)]; |

324 | vct2str(w); |

325 | vct2str(w,1); |

326 | // TODO: WHAT ABOUT ZERO? |

327 | vector M = 0; |

328 | vct2str(M); |

329 | M = [0]; |

330 | vct2str(M); |

331 | M = [0, x, y3, z(1)]; |

332 | vct2str(M,1); // TODO: FIX BUG: MUST BE ZERO! |

333 | M = [1, x, 0, z(1)]; |

334 | vct2str(M,1); |

335 | M = [1, 0, z(1), 0, 0, x*y3, 0, 666]; |

336 | vct2str(M,1); // TODO: FIX BUG: WHAT IS THIS??? |

337 | |

338 | } |

339 | |

340 | proc IsVar(poly p) |

341 | { |

342 | // checks whether p is a variable indeed |

343 | // if it's a power of a variable, returns the power |

344 | if (p==0) { return(0); } //"p=0"; |

345 | poly q = leadmonom(p); |

346 | if ( (p-lead(p)) !=0 ) { return(0); } // "p-lm(p)>0"; |

347 | intvec v = leadexp(p); |

348 | int s = size(v); |

349 | int i=1; |

350 | int cnt = 0; |

351 | int pwr = 0; |

352 | for (i=1; i<=s; i++) |

353 | { |

354 | if (v[i] != 0) |

355 | { |

356 | cnt++; |

357 | pwr = v[i]; |

358 | } |

359 | } |

360 | // "cnt:"; cnt; |

361 | if (cnt==1) { return(pwr); } |

362 | else { return(0); } |

363 | } |

364 | example |

365 | { |

366 | "EXAMPLE:"; echo = 2; |

367 | ring r = 0,(x,y),dp; |

368 | poly f = xy+1; |

369 | IsVar(f); |

370 | poly g = xy; |

371 | IsVar(g); |

372 | poly h = y^3; |

373 | IsVar(h); |

374 | poly i = 1; |

375 | IsVar(i); |

376 | } |

377 | |

378 | // new conversion routines |

379 | |

380 | proc id2words(ideal I, int d) |

381 | { |

382 | // NOT FINISHED |

383 | // input: ideal I of polys in letter-place notation |

384 | // in the ring with d real vars |

385 | // output: the list of strings: associative words |

386 | // extract names of vars |

387 | int i,m,n; string s; string place = "(1)"; |

388 | list lv; |

389 | for(i=1; i<=d; i++) |

390 | { |

391 | s = string(var(i)); |

392 | // get rid of place |

393 | n = find(s, place); |

394 | if (n>0) |

395 | { |

396 | s = s[1..n-1]; |

397 | } |

398 | lv[i] = s; |

399 | } |

400 | poly p,q; |

401 | for (i=1; i<=ncols(I); i++) |

402 | { |

403 | if (I[i] != 0) |

404 | { |

405 | p = I[i]; |

406 | while (p!=0) |

407 | { |

408 | q = leadmonom(p); |

409 | } |

410 | } |

411 | } |

412 | |

413 | return(lv); |

414 | } |

415 | example |

416 | { |

417 | "EXAMPLE:"; echo = 2; |

418 | ring r = 0,(x(1),y(1),z(1),x(2),y(2),z(2)),dp; |

419 | ideal I = x(1)*y(2) -z(1)*x(2); |

420 | id2words(I,3); |

421 | } |

422 | |

423 | proc mono2word(poly p, int d) |

424 | { |

425 | } |

426 | |

427 | // given the element -7xy^2x, it is represented as [-7,x,y^2,x] or as [-7,x,y,y,x] |

428 | // use the orig ord on (x,y,z) and expand it blockwise to (x(i),y(i),z(i)) |

429 | |

430 | // the correspondences: |

431 | // monomial in K<x,y,z> <<--->> vector in R |

432 | // polynomial in K<x,y,z> <<--->> list of vectors (matrix/module) in R |

433 | // ideal in K<x,y,z> <<--->> list of matrices/modules in R |

434 | |

435 | |

436 | // 1. form a new ring |

437 | // 2. NOP |

438 | // 3. compute GB -> with the kernel stuff |

439 | // 4. skip shifted elts (check that no such exist?) |

440 | // 5. go back to orig vars, produce strings/modules |

441 | // 6. return the result |

442 | |

443 | |

444 | // TODO: BAD NAME -> RENAME! |

445 | // TODO: NO ASSUMPTIONS? WHAT ABOUT NON-COMM. INPUT RING? NON-HOMOG. INPUT? |

446 | // TODO: ADD STEP-BY-STEP COMMENTS. |

447 | |

448 | proc freegbasis(list LM, int d) |

449 | "USAGE: freegbasis(L, d); L a list of modules, d an integer |

450 | RETURN: ring, |

451 | TODO: EXPLAIN OUTPUT FORMAT |

452 | PURPOSE: compute the two-sided Groebner basis of an ideal, encoded by L |

453 | (TODO: EXPLAIN FORMAT OF L) |

454 | in the free associative algebra, |

455 | (TODO: WHICH free associative algebra? SYMBOLS? GROUND FIELD?) |

456 | up to degree d |

457 | BACKGROUND: |

458 | TODO: EXPLAIN ALGORITHM, GIVE REFERENCES! |

459 | TODO: WHAT ABOUT TERMINATION? |

460 | DISPLAY: |

461 | TODO: EXPLAIN! |

462 | EXAMPLE: example freegbasis; shows examples |

463 | " |

464 | { |

465 | // d = up to degree, will be shifted to d+1 |

466 | if (d<1) {"bad d"; return(0);} |

467 | |

468 | int ppl = printlevel-voice+2; |

469 | string err = ""; |

470 | |

471 | int i,j,s; |

472 | def save = basering; |

473 | // determine max no of places in the input |

474 | int slm = size(LM); // numbers of polys in the ideal |

475 | int sm; |

476 | intvec iv; |

477 | module M; |

478 | for (i=1; i<=slm; i++) |

479 | { |

480 | // modules, e.g. free polynomials |

481 | M = LM[i]; |

482 | sm = ncols(M); |

483 | for (j=1; j<=sm; j++) |

484 | { |

485 | //vectors, e.g. free monomials |

486 | iv = iv, size(M[j])-1; // 1 place is reserved by the coeff |

487 | } |

488 | } |

489 | |

490 | // TODO: THE FOLLOWING SEEMS TO BE THE CONTENT OF 'freegbRing'! |

491 | // AND THUS IS OVERSIMPLIFIED: WHAT IF VARIABLES/PARAMETERS HAVE BRACKETS? |

492 | // SEE EXAMPLE FOR 'freegbRing'! |

493 | |

494 | int D = Max(iv); // max size of input words |

495 | if (d<D) {"bad d"; return(LM);} |

496 | D = D + d-1; |

497 | // D = d; |

498 | list LR = ringlist(save); |

499 | list L, tmp; |

500 | L[1] = LR[1]; // ground field |

501 | L[4] = LR[4]; // quotient ideal |

502 | tmp = LR[2]; // varnames |

503 | s = size(LR[2]); |

504 | for (i=1; i<=D; i++) |

505 | { |

506 | for (j=1; j<=s; j++) |

507 | { |

508 | tmp[i*s+j] = string(tmp[j])+"("+string(i+1)+")"; |

509 | } |

510 | } |

511 | for (i=1; i<=s; i++) |

512 | { |

513 | tmp[i] = string(tmp[i])+"("+string(1)+")"; |

514 | } |

515 | L[2] = tmp; |

516 | list OrigNames = LR[2]; |

517 | // ordering: d blocks of the ord on r |

518 | // try to get whether the ord on r is blockord itself |

519 | s = size(LR[3]); |

520 | if (s==2) |

521 | { |

522 | // not a blockord, 1 block + module ord |

523 | tmp = LR[3][s]; // module ord |

524 | for (i=1; i<=D; i++) |

525 | { |

526 | LR[3][s-1+i] = LR[3][1]; |

527 | } |

528 | LR[3][s+D] = tmp; |

529 | } |

530 | if (s>2) |

531 | { |

532 | // there are s-1 blocks |

533 | int nb = s-1; |

534 | tmp = LR[3][s]; // module ord |

535 | for (i=1; i<=D; i++) |

536 | { |

537 | for (j=1; j<=nb; j++) |

538 | { |

539 | LR[3][i*nb+j] = LR[3][j]; |

540 | } |

541 | } |

542 | // size(LR[3]); |

543 | LR[3][nb*(D+1)+1] = tmp; |

544 | } |

545 | L[3] = LR[3]; |

546 | def @R = ring(L); |

547 | |

548 | // HERE END 'freegbRing'! ;-) |

549 | |

550 | |

551 | // TODO: THE FOLLOWING CONVERSION DESERVES TO BE AN AUXILIARY PROCEDURE! |

552 | |

553 | setring @R; // ! |

554 | ideal I; |

555 | poly @p; |

556 | s = size(OrigNames); |

557 | // "s:";s; |

558 | // convert LM to canonical vectors (no powers) |

559 | setring save; // !! |

560 | kill M; // M was defined earlier |

561 | module M; |

562 | slm = size(LM); // numbers of polys in the ideal |

563 | int sv,k,l; |

564 | vector v; |

565 | // poly p; |

566 | string sp; |

567 | setring @R; // !!! |

568 | poly @@p=0; |

569 | setring save; // !!!! |

570 | for (l=1; l<=slm; l++) |

571 | { |

572 | // modules, e.g. free polynomials |

573 | M = LM[l]; |

574 | sm = ncols(M); // in intvec iv the sizes are stored |

575 | // modules, e.g. free polynomials |

576 | for (j=1; j<=sm; j++) |

577 | { |

578 | //vectors, e.g. free monomials |

579 | v = M[j]; |

580 | sv = size(v); |

581 | // "sv:";sv; |

582 | sp = "@@p = @@p + "; |

583 | for (k=2; k<=sv; k++) |

584 | { |

585 | sp = sp + string(v[k])+"("+string(k-1)+")*"; |

586 | } |

587 | sp = sp + string(v[1])+";"; // coef; |

588 | setring @R; // !!!!! |

589 | execute(sp); |

590 | setring save; // !!!!! ! |

591 | } |

592 | setring @R; // !!!!! !! |

593 | // "@@p:"; @@p; |

594 | I = I,@@p; |

595 | @@p = 0; |

596 | setring save; // !!!!! !!! |

597 | } |

598 | kill sp; |

599 | // 3. compute GB |

600 | setring @R; // !!!!! !!!! |

601 | |

602 | // TODO: UNBELIEVABLE PLENITUDE OF SETRING ABOVE! |

603 | // YOU CAN DO BETTER AS FOLLOWS: |

604 | // 1. GENERATE ONLY 1 (ONE!!!!) STRING sp FOR THE WHOLE(!!!) L |

605 | // BEFORE(!!!!) SWITCHING TO @R |

606 | // 2. EXECUTE IT IN @R |

607 | |

608 | dbprint(ppl,"computing GB"); |

609 | ideal J = system("freegb",I,d,nvars(save)); |

610 | // ideal J = slimgb(I); |

611 | dbprint(ppl,J); |

612 | // 4. skip shifted elts |

613 | ideal K = select1(J,1..s); // s = size(OrigNames) |

614 | dbprint(ppl,K); |

615 | dbprint(ppl, "done with GB"); |

616 | |

617 | // TODO: THE FOLLOWING DESERVES TO BE AN AUXILIARY PROCEDURE AS WELL! IS IT lp2lstr? |

618 | |

619 | // K contains vars x(1),...z(1) = images of originals |

620 | // 5. go back to orig vars, produce strings/modules |

621 | if (K[1] == 0) |

622 | { |

623 | "no reasonable output, GB gives 0"; |

624 | return(0); |

625 | } |

626 | |

627 | // YOU HAVE A LOT OF SETRINGS BELOW, ARE YOU SURE THAT YOU CANNOT ELIMINATE SOME...? |

628 | |

629 | int sk = size(K); |

630 | int sp, sx, a, b; |

631 | intvec x; |

632 | poly p,q; |

633 | poly pn; |

634 | // vars in 'save' |

635 | setring save; // ! |

636 | module N; |

637 | list LN; |

638 | vector V; |

639 | poly pn; |

640 | // test and skip exponents >=2 |

641 | setring @R; // !! |

642 | for(i=1; i<=sk; i++) |

643 | { |

644 | p = K[i]; |

645 | while (p!=0) |

646 | { |

647 | q = lead(p); |

648 | // "processing q:";q; |

649 | x = leadexp(q); |

650 | sx = size(x); |

651 | for(k=1; k<=sx; k++) |

652 | { |

653 | if ( x[k] >= 2 ) |

654 | { |

655 | err = "skip: the value x[k] is " + string(x[k]); |

656 | dbprint(ppl,err); |

657 | // return(0); |

658 | K[i] = 0; |

659 | p = 0; |

660 | q = 0; |

661 | break; |

662 | } |

663 | } |

664 | p = p - q; |

665 | } |

666 | } |

667 | K = simplify(K,2); |

668 | sk = size(K); |

669 | |

670 | for(i=1; i<=sk; i++) |

671 | { |

672 | // setring save; |

673 | // V = 0; |

674 | setring @R; // !!! |

675 | p = K[i]; |

676 | while (p!=0) |

677 | { |

678 | q = lead(p); |

679 | err = "processing q:" + string(q); |

680 | dbprint(ppl,err); |

681 | x = leadexp(q); |

682 | sx = size(x); |

683 | pn = leadcoef(q); |

684 | setring save; // !!!! |

685 | pn = imap(@R,pn); |

686 | V = V + leadcoef(pn)*gen(1); |

687 | for(k=1; k<=sx; k++) |

688 | { |

689 | if (x[k] ==1) |

690 | { |

691 | a = k / s; // block number=a+1, a!=0 |

692 | b = k % s; // remainder |

693 | // printf("a: %s, b: %s",a,b); |

694 | if (b == 0) |

695 | { |

696 | // that is it's the last var in the block |

697 | b = s; |

698 | a = a-1; |

699 | } |

700 | V = V + var(b)*gen(a+2); |

701 | } |

702 | // else |

703 | // { |

704 | // printf("error: the value x[k] is %s", x[k]); |

705 | // return(0); |

706 | // } |

707 | } |

708 | err = "V: " + string(V); |

709 | dbprint(ppl,err); |

710 | // printf("V: %s", string(V)); |

711 | N = N,V; |

712 | V = 0; |

713 | setring @R; // !!!!! |

714 | p = p - q; |

715 | pn = 0; |

716 | } |

717 | setring save; // !!!!! ! |

718 | LN[i] = simplify(N,2); |

719 | N = 0; |

720 | } |

721 | setring save; // !!!!! !! |

722 | return(LN); |

723 | } |

724 | example |

725 | { |

726 | "EXAMPLE:"; echo = 2; |

727 | ring r = 0,(x,y,z),(dp(1),dp(2)); |

728 | module M = [-1,x,y],[-7,y,y],[3,x,x]; |

729 | module N = [1,x,y,x],[-1,y,x,y]; |

730 | list L; L[1] = M; L[2] = N; |

731 | lst2str(L); |

732 | def U = freegbasis(L,5); |

733 | lst2str(U); |

734 | kill U; |

735 | setring r; // non-homog. input: |

736 | M = [-1,x,y],[-7,y,y,z],[3,x,x,x,z,z]; |

737 | N = [1,x,y,z,x],[-1,y,x,y]; |

738 | L[1] = M; L[2] = N; |

739 | lst2str(L); |

740 | def U = freegbasis(L,5); |

741 | lst2str(U); // OK |

742 | kill U,r; |

743 | ring r = 0,(x,y,z),(dp(1),dp(2)); |

744 | def R = nc_algebra(1,0); setring R; // should be the same as 1st! |

745 | module M = [-1,x,y],[-7,y,y],[3,x,x]; |

746 | module N = [1,x,y,x],[-1,y,x,y]; |

747 | list L; L[1] = M; L[2] = N; |

748 | lst2str(L); |

749 | def U = freegbasis(L,5); |

750 | lst2str(U); // OK |

751 | kill U, R; |

752 | setring r; |

753 | def R = nc_algebra(-1,1); setring R; // some non-commutativity |

754 | R; |

755 | module M = [-1,x,y],[-7,y,y],[3,x,x]; |

756 | module N = [1,x,y,x],[-1,y,x,y]; |

757 | list L; L[1] = M; L[2] = N; |

758 | lst2str(L); |

759 | def U = freegbasis(L,5); |

760 | lst2str(U); |

761 | |

762 | |

763 | } |

764 | |

765 | proc crs(list LM, int d) |

766 | "USAGE: crs(L, d); L a list of modules, d an integer |

767 | RETURN: ring |

768 | PURPOSE: create a ring and shift the ideal |

769 | EXAMPLE: example crs; shows examples |

770 | " |

771 | { |

772 | // d = up to degree, will be shifted to d+1 |

773 | if (d<1) {"bad d"; return(0);} |

774 | |

775 | int ppl = printlevel-voice+2; |

776 | string err = ""; |

777 | |

778 | int i,j,s; |

779 | def save = basering; |

780 | // determine max no of places in the input |

781 | int slm = size(LM); // numbers of polys in the ideal |

782 | int sm; |

783 | intvec iv; |

784 | module M; |

785 | for (i=1; i<=slm; i++) |

786 | { |

787 | // modules, e.g. free polynomials |

788 | M = LM[i]; |

789 | sm = ncols(M); |

790 | for (j=1; j<=sm; j++) |

791 | { |

792 | //vectors, e.g. free monomials |

793 | iv = iv, size(M[j])-1; // 1 place is reserved by the coeff |

794 | } |

795 | } |

796 | int D = Max(iv); // max size of input words |

797 | if (d<D) {"bad d"; return(LM);} |

798 | D = D + d-1; |

799 | // D = d; |

800 | list LR = ringlist(save); |

801 | list L, tmp; |

802 | L[1] = LR[1]; // ground field |

803 | L[4] = LR[4]; // quotient ideal |

804 | tmp = LR[2]; // varnames |

805 | s = size(LR[2]); |

806 | for (i=1; i<=D; i++) |

807 | { |

808 | for (j=1; j<=s; j++) |

809 | { |

810 | tmp[i*s+j] = string(tmp[j])+"("+string(i)+")"; |

811 | } |

812 | } |

813 | for (i=1; i<=s; i++) |

814 | { |

815 | tmp[i] = string(tmp[i])+"("+string(0)+")"; |

816 | } |

817 | L[2] = tmp; |

818 | list OrigNames = LR[2]; |

819 | // ordering: d blocks of the ord on r |

820 | // try to get whether the ord on r is blockord itself |

821 | s = size(LR[3]); |

822 | if (s==2) |

823 | { |

824 | // not a blockord, 1 block + module ord |

825 | tmp = LR[3][s]; // module ord |

826 | for (i=1; i<=D; i++) |

827 | { |

828 | LR[3][s-1+i] = LR[3][1]; |

829 | } |

830 | LR[3][s+D] = tmp; |

831 | } |

832 | if (s>2) |

833 | { |

834 | // there are s-1 blocks |

835 | int nb = s-1; |

836 | tmp = LR[3][s]; // module ord |

837 | for (i=1; i<=D; i++) |

838 | { |

839 | for (j=1; j<=nb; j++) |

840 | { |

841 | LR[3][i*nb+j] = LR[3][j]; |

842 | } |

843 | } |

844 | // size(LR[3]); |

845 | LR[3][nb*(D+1)+1] = tmp; |

846 | } |

847 | L[3] = LR[3]; |

848 | def @R = ring(L); |

849 | setring @R; |

850 | ideal I; |

851 | poly @p; |

852 | s = size(OrigNames); |

853 | // "s:";s; |

854 | // convert LM to canonical vectors (no powers) |

855 | setring save; |

856 | kill M; // M was defined earlier |

857 | module M; |

858 | slm = size(LM); // numbers of polys in the ideal |

859 | int sv,k,l; |

860 | vector v; |

861 | // poly p; |

862 | string sp; |

863 | setring @R; |

864 | poly @@p=0; |

865 | setring save; |

866 | for (l=1; l<=slm; l++) |

867 | { |

868 | // modules, e.g. free polynomials |

869 | M = LM[l]; |

870 | sm = ncols(M); // in intvec iv the sizes are stored |

871 | for (i=0; i<=d-iv[l]; i++) |

872 | { |

873 | // modules, e.g. free polynomials |

874 | for (j=1; j<=sm; j++) |

875 | { |

876 | //vectors, e.g. free monomials |

877 | v = M[j]; |

878 | sv = size(v); |

879 | // "sv:";sv; |

880 | sp = "@@p = @@p + "; |

881 | for (k=2; k<=sv; k++) |

882 | { |

883 | sp = sp + string(v[k])+"("+string(k-2+i)+")*"; |

884 | } |

885 | sp = sp + string(v[1])+";"; // coef; |

886 | setring @R; |

887 | execute(sp); |

888 | setring save; |

889 | } |

890 | setring @R; |

891 | // "@@p:"; @@p; |

892 | I = I,@@p; |

893 | @@p = 0; |

894 | setring save; |

895 | } |

896 | } |

897 | setring @R; |

898 | export I; |

899 | return(@R); |

900 | } |

901 | example |

902 | { |

903 | "EXAMPLE:"; echo = 2; |

904 | ring r = 0,(x,y,z),(dp(1),dp(2)); |

905 | module M = [-1,x,y],[-7,y,y],[3,x,x]; |

906 | module N = [1,x,y,x],[-1,y,x,y]; |

907 | list L; L[1] = M; L[2] = N; |

908 | lst2str(L); |

909 | def U = crs(L,5); |

910 | setring U; U; |

911 | I; |

912 | } |

913 | |

914 | proc polylen(ideal I) |

915 | { |

916 | // returns the ideal of length of polys |

917 | int i; |

918 | intvec J; |

919 | number s = 0; |

920 | for(i=1;i<=ncols(I);i++) |

921 | { |

922 | J[i] = size(I[i]); |

923 | s = s + J[i]; |

924 | } |

925 | printf("the sum of length %s",s); |

926 | // print(s); |

927 | return(J); |

928 | } |

929 | |

930 | |

931 | // TODO: ASSUMPTIONS? |

932 | // TODO: PROC FORGETS ABOUT NON-COMM. RELATIONS IN THE |

933 | //// BASERING! DOCUMENT IT OR FIX... |

934 | // TODO: WRITE DOWN WHAT IS THE ORDERING IN THE OUTPUT RING? |

935 | // TODO: OVERSIMPLIFIED: WHAT IF PARAMETERS HAVE BRACKETS |

936 | //// AND COINCIDE WITH GENERATED VARIABLES? SEE EXAMPLE |

937 | |

938 | |

939 | proc freegbRing(int d) |

940 | "USAGE: freegbRing(d); d an integer |

941 | RETURN: ring |

942 | PURPOSE: creates a ring with d blocks of shifted original variables |

943 | EXAMPLE: example freegbRing; shows examples |

944 | " |

945 | { |

946 | // d = up to degree, will be shifted to d+1 |

947 | if (d<1) {"bad d"; return(0);} |

948 | |

949 | int ppl = printlevel-voice+2; |

950 | string err = ""; |

951 | |

952 | int i,j,s; |

953 | def save = basering; |

954 | int D = d-1; |

955 | list LR = ringlist(save); |

956 | list L, tmp; |

957 | L[1] = LR[1]; // ground field |

958 | L[4] = LR[4]; // quotient ideal |

959 | tmp = LR[2]; // varnames |

960 | s = size(LR[2]); |

961 | for (i=1; i<=D; i++) |

962 | { |

963 | for (j=1; j<=s; j++) |

964 | { |

965 | tmp[i*s+j] = string(tmp[j])+"("+string(i+1)+")"; |

966 | } |

967 | } |

968 | for (i=1; i<=s; i++) |

969 | { |

970 | tmp[i] = string(tmp[i])+"("+string(1)+")"; |

971 | } |

972 | L[2] = tmp; |

973 | list OrigNames = LR[2]; |

974 | // ordering: d blocks of the ord on r |

975 | // try to get whether the ord on r is blockord itself |

976 | // TODO: make L(2) ordering! exponent is maximally 2 |

977 | s = size(LR[3]); |

978 | if (s==2) |

979 | { |

980 | // not a blockord, 1 block + module ord |

981 | tmp = LR[3][s]; // module ord |

982 | for (i=1; i<=D; i++) |

983 | { |

984 | LR[3][s-1+i] = LR[3][1]; |

985 | } |

986 | LR[3][s+D] = tmp; |

987 | } |

988 | if (s>2) |

989 | { |

990 | // there are s-1 blocks |

991 | int nb = s-1; |

992 | tmp = LR[3][s]; // module ord |

993 | for (i=1; i<=D; i++) |

994 | { |

995 | for (j=1; j<=nb; j++) |

996 | { |

997 | LR[3][i*nb+j] = LR[3][j]; |

998 | } |

999 | } |

1000 | // size(LR[3]); |

1001 | LR[3][nb*(D+1)+1] = tmp; |

1002 | } |

1003 | L[3] = LR[3]; |

1004 | def @R = ring(L); |

1005 | // setring @R; |

1006 | return (@R); |

1007 | } |

1008 | example |

1009 | { |

1010 | "EXAMPLE:"; echo = 2; |

1011 | ring r = 0,(x,y,z),(dp(1),dp(2)); |

1012 | r; |

1013 | def A = freegbRing(2); setring A; |

1014 | A; |

1015 | kill A, r; |

1016 | ring r = 0,(x(1..3)),(dp(1),lp(2)); |

1017 | r; |

1018 | def A = freegbRing(2); setring A; |

1019 | A; // OK |

1020 | kill A, r; |

1021 | ring r = (0,a(1),b(1)),(a, b),(lp(1),dp(1)); |

1022 | r; |

1023 | def A = freegbRing(2); setring A; |

1024 | A; // BUG: parameter should not be named as a variable and vice verse! |

1025 | a(1); typeof(a(1)); |

1026 | kill A, r; |

1027 | ring r = 0,(x,y,z),dp; |

1028 | def R = nc_algebra(-1, 1); setring R; |

1029 | R; |

1030 | def A = freegbRing(2); setring A; |

1031 | A; // NOTE: the putput is a purely commutative ring! |

1032 | kill A, R, r; |

1033 | } |

1034 | |

1035 | proc ex_shift() |

1036 | { |

1037 | LIB "freegb.lib"; |

1038 | ring r = 0,(x,y,z),(dp(1),dp(2)); |

1039 | module M = [-1,x,y],[-7,y,y],[3,x,x]; |

1040 | module N = [1,x,y,x],[-1,y,x,y]; |

1041 | list L; L[1] = M; L[2] = N; |

1042 | lst2str(L); |

1043 | def U = crs(L,5); |

1044 | setring U; U; |

1045 | I; |

1046 | poly p = I[2]; // I[8]; |

1047 | p; |

1048 | system("stest",p,7,7,3); // error -> the world is ok |

1049 | poly q1 = system("stest",p,1,7,3); //ok |

1050 | poly q6 = system("stest",p,6,7,3); //ok |

1051 | system("btest",p,3); //ok |

1052 | system("btest",q1,3); //ok |

1053 | system("btest",q6,3); //ok |

1054 | } |

1055 | |

1056 | proc test_shrink() |

1057 | { |

1058 | LIB "freegb.lib"; |

1059 | ring r =0,(x,y,z),dp; |

1060 | int d = 5; |

1061 | def R = freegbRing(d); |

1062 | setring R; |

1063 | poly p1 = x(1)*y(2)*z(3); |

1064 | poly p2 = x(1)*y(4)*z(5); |

1065 | poly p3 = x(1)*y(1)*z(3); |

1066 | poly p4 = x(1)*y(2)*z(2); |

1067 | poly p5 = x(3)*z(5); |

1068 | poly p6 = x(1)*y(1)*x(3)*z(5); |

1069 | poly p7 = x(1)*y(2)*x(3)*y(4)*z(5); |

1070 | poly p8 = p1+p2+p3+p4+p5 + p6 + p7; |

1071 | p1; system("shrinktest",p1,3); |

1072 | p2; system("shrinktest",p2,3); |

1073 | p3; system("shrinktest",p3,3); |

1074 | p4; system("shrinktest",p4,3); |

1075 | p5; system("shrinktest",p5,3); |

1076 | p6; system("shrinktest",p6,3); |

1077 | p7; system("shrinktest",p7,3); |

1078 | p8; system("shrinktest",p8,3); |

1079 | poly p9 = p1 + 2*p2 + 5*p5 + 7*p7; |

1080 | p9; system("shrinktest",p9,3); |

1081 | } |

1082 | |

1083 | proc ex2() |

1084 | { |

1085 | option(prot); |

1086 | LIB "freegb.lib"; |

1087 | ring r = 0,(x,y),dp; |

1088 | module M = [-1,x,y],[3,x,x]; // 3x^2 - xy |

1089 | def U = freegb(M,7); |

1090 | lst2str(U); |

1091 | } |

1092 | |

1093 | proc ex_nonhomog() |

1094 | { |

1095 | option(prot); |

1096 | LIB "freegb.lib"; |

1097 | ring r = 0,(x,y,h),dp; |

1098 | list L; |

1099 | module M; |

1100 | M = [-1,y,y],[1,x,x,x]; // x3-y2 |

1101 | L[1] = M; |

1102 | M = [1,x,h],[-1,h,x]; // xh-hx |

1103 | L[2] = M; |

1104 | M = [1,y,h],[-1,h,y]; // yh-hy |

1105 | L[3] = M; |

1106 | def U = freegb(L,4); |

1107 | lst2str(U); |

1108 | // strange elements in the basis |

1109 | } |

1110 | |

1111 | proc ex_nonhomog_comm() |

1112 | { |

1113 | option(prot); |

1114 | LIB "freegb.lib"; |

1115 | ring r = 0,(x,y),dp; |

1116 | module M = [-1,y,y],[1,x,x,x]; |

1117 | def U = freegb(M,5); |

1118 | lst2str(U); |

1119 | } |

1120 | |

1121 | proc ex_nonhomog_h() |

1122 | { |

1123 | option(prot); |

1124 | LIB "freegb.lib"; |

1125 | ring r = 0,(x,y,h),(a(1,1),dp); |

1126 | module M = [-1,y,y,h],[1,x,x,x]; // x3 - y2h |

1127 | def U = freegb(M,6); |

1128 | lst2str(U); |

1129 | } |

1130 | |

1131 | proc ex_nonhomog_h2() |

1132 | { |

1133 | option(prot); |

1134 | LIB "freegb.lib"; |

1135 | ring r = 0,(x,y,h),(dp); |

1136 | list L; |

1137 | module M; |

1138 | M = [-1,y,y,h],[1,x,x,x]; // x3 - y2h |

1139 | L[1] = M; |

1140 | M = [1,x,h],[-1,h,x]; // xh - hx |

1141 | L[2] = M; |

1142 | M = [1,y,h],[-1,h,y]; // yh - hy |

1143 | L[3] = M; |

1144 | def U = freegbasis(L,3); |

1145 | lst2str(U); |

1146 | // strange answer CHECK |

1147 | } |

1148 | |

1149 | |

1150 | proc ex_nonhomog_3() |

1151 | { |

1152 | option(prot); |

1153 | LIB "./freegb.lib"; |

1154 | ring r = 0,(x,y,z),(dp); |

1155 | list L; |

1156 | module M; |

1157 | M = [1,z,y],[-1,x]; // zy - x |

1158 | L[1] = M; |

1159 | M = [1,z,x],[-1,y]; // zx - y |

1160 | L[2] = M; |

1161 | M = [1,y,x],[-1,z]; // yx - z |

1162 | L[3] = M; |

1163 | lst2str(L); |

1164 | list U = freegb(L,4); |

1165 | lst2str(U); |

1166 | // strange answer CHECK |

1167 | } |

1168 | |

1169 | proc ex_densep_2() |

1170 | { |

1171 | option(prot); |

1172 | LIB "freegb.lib"; |

1173 | ring r = (0,a,b,c),(x,y),(Dp); // deglex |

1174 | module M = [1,x,x], [a,x,y], [b,y,x], [c,y,y]; |

1175 | lst2str(M); |

1176 | list U = freegb(M,5); |

1177 | lst2str(U); |

1178 | // a=b is important -> finite basis!!! |

1179 | module M = [1,x,x], [a,x,y], [a,y,x], [c,y,y]; |

1180 | lst2str(M); |

1181 | list U = freegb(M,5); |

1182 | lst2str(U); |

1183 | } |

1184 | |

1185 | |

1186 | // 1. form a new ring |

1187 | // 2. produce shifted generators |

1188 | // 3. compute GB |

1189 | // 4. skip shifted elts |

1190 | // 5. go back to orig vars, produce strings/modules |

1191 | // 6. return the result |

1192 | |

1193 | proc freegbold(list LM, int d) |

1194 | "USAGE: freegbold(L, d); L a list of modules, d an integer |

1195 | RETURN: ring |

1196 | PURPOSE: compute the two-sided Groebner basis of an ideal, encoded by L in |

1197 | the free associative algebra, up to degree d |

1198 | EXAMPLE: example freegbold; shows examples |

1199 | " |

1200 | { |

1201 | // d = up to degree, will be shifted to d+1 |

1202 | if (d<1) {"bad d"; return(0);} |

1203 | |

1204 | int ppl = printlevel-voice+2; |

1205 | string err = ""; |

1206 | |

1207 | int i,j,s; |

1208 | def save = basering; |

1209 | // determine max no of places in the input |

1210 | int slm = size(LM); // numbers of polys in the ideal |

1211 | int sm; |

1212 | intvec iv; |

1213 | module M; |

1214 | for (i=1; i<=slm; i++) |

1215 | { |

1216 | // modules, e.g. free polynomials |

1217 | M = LM[i]; |

1218 | sm = ncols(M); |

1219 | for (j=1; j<=sm; j++) |

1220 | { |

1221 | //vectors, e.g. free monomials |

1222 | iv = iv, size(M[j])-1; // 1 place is reserved by the coeff |

1223 | } |

1224 | } |

1225 | int D = Max(iv); // max size of input words |

1226 | if (d<D) {"bad d"; return(LM);} |

1227 | D = D + d-1; |

1228 | // D = d; |

1229 | list LR = ringlist(save); |

1230 | list L, tmp; |

1231 | L[1] = LR[1]; // ground field |

1232 | L[4] = LR[4]; // quotient ideal |

1233 | tmp = LR[2]; // varnames |

1234 | s = size(LR[2]); |

1235 | for (i=1; i<=D; i++) |

1236 | { |

1237 | for (j=1; j<=s; j++) |

1238 | { |

1239 | tmp[i*s+j] = string(tmp[j])+"("+string(i+1)+")"; |

1240 | } |

1241 | } |

1242 | for (i=1; i<=s; i++) |

1243 | { |

1244 | tmp[i] = string(tmp[i])+"("+string(1)+")"; |

1245 | } |

1246 | L[2] = tmp; |

1247 | list OrigNames = LR[2]; |

1248 | // ordering: d blocks of the ord on r |

1249 | // try to get whether the ord on r is blockord itself |

1250 | // TODO: make L(2) ordering! exponent is maximally 2 |

1251 | s = size(LR[3]); |

1252 | if (s==2) |

1253 | { |

1254 | // not a blockord, 1 block + module ord |

1255 | tmp = LR[3][s]; // module ord |

1256 | for (i=1; i<=D; i++) |

1257 | { |

1258 | LR[3][s-1+i] = LR[3][1]; |

1259 | } |

1260 | LR[3][s+D] = tmp; |

1261 | } |

1262 | if (s>2) |

1263 | { |

1264 | // there are s-1 blocks |

1265 | int nb = s-1; |

1266 | tmp = LR[3][s]; // module ord |

1267 | for (i=1; i<=D; i++) |

1268 | { |

1269 | for (j=1; j<=nb; j++) |

1270 | { |

1271 | LR[3][i*nb+j] = LR[3][j]; |

1272 | } |

1273 | } |

1274 | // size(LR[3]); |

1275 | LR[3][nb*(D+1)+1] = tmp; |

1276 | } |

1277 | L[3] = LR[3]; |

1278 | def @R = ring(L); |

1279 | setring @R; |

1280 | ideal I; |

1281 | poly @p; |

1282 | s = size(OrigNames); |

1283 | // "s:";s; |

1284 | // convert LM to canonical vectors (no powers) |

1285 | setring save; |

1286 | kill M; // M was defined earlier |

1287 | module M; |

1288 | slm = size(LM); // numbers of polys in the ideal |

1289 | int sv,k,l; |

1290 | vector v; |

1291 | // poly p; |

1292 | string sp; |

1293 | setring @R; |

1294 | poly @@p=0; |

1295 | setring save; |

1296 | for (l=1; l<=slm; l++) |

1297 | { |

1298 | // modules, e.g. free polynomials |

1299 | M = LM[l]; |

1300 | sm = ncols(M); // in intvec iv the sizes are stored |

1301 | for (i=0; i<=d-iv[l]; i++) |

1302 | { |

1303 | // modules, e.g. free polynomials |

1304 | for (j=1; j<=sm; j++) |

1305 | { |

1306 | //vectors, e.g. free monomials |

1307 | v = M[j]; |

1308 | sv = size(v); |

1309 | // "sv:";sv; |

1310 | sp = "@@p = @@p + "; |

1311 | for (k=2; k<=sv; k++) |

1312 | { |

1313 | sp = sp + string(v[k])+"("+string(k-1+i)+")*"; |

1314 | } |

1315 | sp = sp + string(v[1])+";"; // coef; |

1316 | setring @R; |

1317 | execute(sp); |

1318 | setring save; |

1319 | } |

1320 | setring @R; |

1321 | // "@@p:"; @@p; |

1322 | I = I,@@p; |

1323 | @@p = 0; |

1324 | setring save; |

1325 | } |

1326 | } |

1327 | kill sp; |

1328 | // 3. compute GB |

1329 | setring @R; |

1330 | dbprint(ppl,"computing GB"); |

1331 | // ideal J = groebner(I); |

1332 | ideal J = slimgb(I); |

1333 | dbprint(ppl,J); |

1334 | // 4. skip shifted elts |

1335 | ideal K = select1(J,1..s); // s = size(OrigNames) |

1336 | dbprint(ppl,K); |

1337 | dbprint(ppl, "done with GB"); |

1338 | // K contains vars x(1),...z(1) = images of originals |

1339 | // 5. go back to orig vars, produce strings/modules |

1340 | if (K[1] == 0) |

1341 | { |

1342 | "no reasonable output, GB gives 0"; |

1343 | return(0); |

1344 | } |

1345 | int sk = size(K); |

1346 | int sp, sx, a, b; |

1347 | intvec x; |

1348 | poly p,q; |

1349 | poly pn; |

1350 | // vars in 'save' |

1351 | setring save; |

1352 | module N; |

1353 | list LN; |

1354 | vector V; |

1355 | poly pn; |

1356 | // test and skip exponents >=2 |

1357 | setring @R; |

1358 | for(i=1; i<=sk; i++) |

1359 | { |

1360 | p = K[i]; |

1361 | while (p!=0) |

1362 | { |

1363 | q = lead(p); |

1364 | // "processing q:";q; |

1365 | x = leadexp(q); |

1366 | sx = size(x); |

1367 | for(k=1; k<=sx; k++) |

1368 | { |

1369 | if ( x[k] >= 2 ) |

1370 | { |

1371 | err = "skip: the value x[k] is " + string(x[k]); |

1372 | dbprint(ppl,err); |

1373 | // return(0); |

1374 | K[i] = 0; |

1375 | p = 0; |

1376 | q = 0; |

1377 | break; |

1378 | } |

1379 | } |

1380 | p = p - q; |

1381 | } |

1382 | } |

1383 | K = simplify(K,2); |

1384 | sk = size(K); |

1385 | for(i=1; i<=sk; i++) |

1386 | { |

1387 | // setring save; |

1388 | // V = 0; |

1389 | setring @R; |

1390 | p = K[i]; |

1391 | while (p!=0) |

1392 | { |

1393 | q = lead(p); |

1394 | err = "processing q:" + string(q); |

1395 | dbprint(ppl,err); |

1396 | x = leadexp(q); |

1397 | sx = size(x); |

1398 | pn = leadcoef(q); |

1399 | setring save; |

1400 | pn = imap(@R,pn); |

1401 | V = V + leadcoef(pn)*gen(1); |

1402 | for(k=1; k<=sx; k++) |

1403 | { |

1404 | if (x[k] ==1) |

1405 | { |

1406 | a = k / s; // block number=a+1, a!=0 |

1407 | b = k % s; // remainder |

1408 | // printf("a: %s, b: %s",a,b); |

1409 | if (b == 0) |

1410 | { |

1411 | // that is it's the last var in the block |

1412 | b = s; |

1413 | a = a-1; |

1414 | } |

1415 | V = V + var(b)*gen(a+2); |

1416 | } |

1417 | // else |

1418 | // { |

1419 | // printf("error: the value x[k] is %s", x[k]); |

1420 | // return(0); |

1421 | // } |

1422 | } |

1423 | err = "V: " + string(V); |

1424 | dbprint(ppl,err); |

1425 | // printf("V: %s", string(V)); |

1426 | N = N,V; |

1427 | V = 0; |

1428 | setring @R; |

1429 | p = p - q; |

1430 | pn = 0; |

1431 | } |

1432 | setring save; |

1433 | LN[i] = simplify(N,2); |

1434 | N = 0; |

1435 | } |

1436 | setring save; |

1437 | return(LN); |

1438 | } |

1439 | example |

1440 | { |

1441 | "EXAMPLE:"; echo = 2; |

1442 | ring r = 0,(x,y,z),(dp(1),dp(2)); |

1443 | module M = [-1,x,y],[-7,y,y],[3,x,x]; |

1444 | module N = [1,x,y,x],[-1,y,x,y]; |

1445 | list L; L[1] = M; L[2] = N; |

1446 | lst2str(L); |

1447 | def U = freegbold(L,5); |

1448 | lst2str(U); |

1449 | } |

1450 | |

1451 | proc sgb(ideal I, int d) |

1452 | { |

1453 | // new code |

1454 | // map x_i to x_i(1) via map() |

1455 | //LIB "freegb.lib"; |

1456 | def save = basering; |

1457 | //int d =7;// degree |

1458 | int nv = nvars(save); |

1459 | def R = freegbRing(d); |

1460 | setring R; |

1461 | int i; |

1462 | ideal Imap; |

1463 | for (i=1; i<=nv; i++) |

1464 | { |

1465 | Imap[i] = var(i); |

1466 | } |

1467 | //ideal I = x(1)*y(2), y(1)*x(2)+z(1)*z(2); |

1468 | ideal I = x(1)*x(2),x(1)*y(2) + z(1)*x(2); |

1469 | option(prot); |

1470 | //option(teach); |

1471 | ideal J = system("freegb",I,d,nv); |

1472 | } |

1473 | |

1474 | |

1475 | |

1476 | static proc checkCeq() |

1477 | { |

1478 | ring r = 0,(x,y),Dp; |

1479 | def A = freegbRing(4); |

1480 | setring A; |

1481 | A; |

1482 | // I = x2-xy |

1483 | ideal I = x(1)*x(2) - x(1)*y(2), x(2)*x(3) - x(2)*y(3), x(3)*x(4) - x(3)*y(4); |

1484 | ideal C = x(2)-x(1),x(3)-x(2),x(4)-x(3),y(2)-y(1),y(3)-y(2),y(4)-y(3); |

1485 | ideal K = I,C; |

1486 | groebner(K); |

1487 | } |

1488 | |

1489 | |

1490 | proc exHom1() |

1491 | { |

1492 | // we start with |

1493 | // z*y - x, z*x - y, y*x - z |

1494 | LIB "freegb.lib"; |

1495 | LIB "elim.lib"; |

1496 | ring r = 0,(x,y,z,h),dp; |

1497 | list L; |

1498 | module M; |

1499 | M = [1,z,y],[-1,x,h]; // zy - xh |

1500 | L[1] = M; |

1501 | M = [1,z,x],[-1,y,h]; // zx - yh |

1502 | L[2] = M; |

1503 | M = [1,y,x],[-1,z,h]; // yx - zh |

1504 | L[3] = M; |

1505 | lst2str(L); |

1506 | def U = crs(L,4); |

1507 | setring U; |

1508 | I = I, |

1509 | y(2)*h(3)+z(2)*x(3), y(3)*h(4)+z(3)*x(4), |

1510 | y(2)*x(3)-z(2)*h(3), y(3)*x(4)-z(3)*h(4); |

1511 | I = simplify(I,2); |

1512 | ring r2 = 0,(x(0..4),y(0..4),z(0..4),h(0..4)),dp; |

1513 | ideal J = imap(U,I); |

1514 | // ideal K = homog(J,h); |

1515 | option(redSB); |

1516 | option(redTail); |

1517 | ideal L = groebner(J); //(K); |

1518 | ideal LL = sat(L,ideal(h))[1]; |

1519 | ideal M = subst(LL,h,1); |

1520 | M = simplify(M,2); |

1521 | setring U; |

1522 | ideal M = imap(r2,M); |

1523 | lst2str(U); |

1524 | } |

1525 | |

1526 | static proc test1() |

1527 | { |

1528 | LIB "freegb.lib"; |

1529 | ring r = 0,(x,y),Dp; |

1530 | int d = 10; // degree |

1531 | def R = freegbRing(d); |

1532 | setring R; |

1533 | ideal I = x(1)*x(2) - y(1)*y(2); |

1534 | option(prot); |

1535 | option(teach); |

1536 | ideal J = system("freegb",I,d,2); |

1537 | J; |

1538 | } |

1539 | |

1540 | static proc test2() |

1541 | { |

1542 | LIB "freegb.lib"; |

1543 | ring r = 0,(x,y),Dp; |

1544 | int d = 10; // degree |

1545 | def R = freegbRing(d); |

1546 | setring R; |

1547 | ideal I = x(1)*x(2) - x(1)*y(2); |

1548 | option(prot); |

1549 | option(teach); |

1550 | ideal J = system("freegb",I,d,2); |

1551 | J; |

1552 | } |

1553 | |

1554 | static proc test3() |

1555 | { |

1556 | LIB "freegb.lib"; |

1557 | ring r = 0,(x,y,z),dp; |

1558 | int d =5; // degree |

1559 | def R = freegbRing(d); |

1560 | setring R; |

1561 | ideal I = x(1)*y(2), y(1)*x(2)+z(1)*z(2); |

1562 | option(prot); |

1563 | option(teach); |

1564 | ideal J = system("freegb",I,d,3); |

1565 | } |

1566 | |

1567 | proc schur2-3() |

1568 | { |

1569 | // nonhomog: |

1570 | // h^4-10*h^2+9,f*e-e*f+h, h*2-e*h-2*e,h*f-f*h+2*f |

1571 | // homogenized with t |

1572 | // h^4-10*h^2*t^2+9*t^4,f*e-e*f+h*t, h*2-e*h-2*e*t,h*f-f*h+2*f*t, |

1573 | // t*h - h*t, t*f - f*t, t*e - e*t |

1574 | } |

1575 | |

1576 | proc adem(int i, int j) |

1577 | { |

1578 | // produces Adem relations for i<2j in char 0 |

1579 | // assume: 0<i<2j |

1580 | // requires presence of vars up to i+j |

1581 | if ( (i<0) || (i >= 2*j) ) |

1582 | { |

1583 | ERROR("arguments out of range"); return(0); |

1584 | } |

1585 | ring @r = 0,(s(i+j..0)),lp; |

1586 | poly p,q; |

1587 | number n; |

1588 | int ii = i div 2; int k; |

1589 | // k=0 => s(0)=1 |

1590 | n = binomial(j-1,i); |

1591 | q = n*s(i+j)*s(0); |

1592 | printf("k=0, term=%s",q); |

1593 | p = p + q; |

1594 | for (k=1; k<= ii; k++) |

1595 | { |

1596 | n = binomial(j-k-1,i-2*k); |

1597 | q = n*s(i+j-k)*s(k);; |

1598 | printf("k=%s, term=%s",k,q); |

1599 | p = p + q; |

1600 | } |

1601 | poly AdemRel = p; |

1602 | export AdemRel; |

1603 | return(@r); |

1604 | } |

1605 | example |

1606 | { |

1607 | "EXAMPLE:"; echo = 2; |

1608 | def A = adem(2,5); |

1609 | setring A; |

1610 | AdemRel; |

1611 | } |

1612 | |

1613 | /* |

1614 | 1,1: 0 |

1615 | 1,2: s(3)*s(0) == s(3) -> def for s(3):=s(1)s(2) |

1616 | 2,1: adm |

1617 | 2,2: s(3)*s(1) == s(1)s(2)s(1) |

1618 | 1,3: 0 ( since 2*s(4)*s(0) = 0 mod 2) |

1619 | 3,1: adm |

1620 | 2,3: s(5)*s(0)+s(4)*s(1) == s(5)+s(4)*s(1) |

1621 | 3,2: 0 |

1622 | 3,3: s(5)*s(1) |

1623 | 1,4: 3*s(5)*s(0) == s(5) -> def for s(5):=s(1)*s(4) |

1624 | 4,1: adm |

1625 | 2,4: 3*s(6)*s(0)+s(5)*s(1) == s(6) + s(5)*s(1) == s(6) + s(1)*s(4)*s(1) |

1626 | 4,2: adm |

1627 | 4,3: s(5)*s(2) |

1628 | 3,4: s(7)*s(0)+2*s(6)*s(1) == s(7) -> def for s(7):=s(3)*s(4) |

1629 | 4,4: s(7)*s(1)+s(6)*s(2) |

1630 | */ |

1631 | |

1632 | /* s1,s2: |

1633 | s1*s1 =0, s2*s2 = s1*s2*s1 |

1634 | */ |

1635 | |

1636 | /* |

1637 | try char 0: |

1638 | s1,s2: |

1639 | s1*s1 =0, s2*s2 = s1*s2*s1, s(1)*s(3)== s(1)*s(1)*s(3) == 0 = 2*s(4) ->def for s(4) |

1640 | hence 2==0! only in char 2 |

1641 | */ |

1642 | |

1643 | proc adem2mod(int n) |

1644 | { |

1645 | // Adem rels modulo 2 |

1646 | } |

1647 | |

1648 | proc stringpoly2lplace(string s) |

1649 | { |

1650 | // decomposes sentence into terms |

1651 | s = replace(s,newline,""); // get rid of newlines |

1652 | s = replace(s," ",""); // get rid of empties |

1653 | //arith symbols: +,- |

1654 | // decompose into words with coeffs |

1655 | list LS; |

1656 | int i,j,ie,je,k,cnt; |

1657 | // s[1]="-" situation |

1658 | if (s[1]=="-") |

1659 | { |

1660 | LS = stringpoly2lplace(string(s[2..size(s)])); |

1661 | LS[1] = string("-"+string(LS[1])); |

1662 | return(LS); |

1663 | } |

1664 | i = find(s,"-",2); |

1665 | // i==1 might happen if the 1st symbol coeff is negative |

1666 | j = find(s,"+"); |

1667 | list LL; |

1668 | if (i==j) |

1669 | { |

1670 | "return a monomial"; |

1671 | // that is both are 0 -> s is a monomial |

1672 | LS[1] = s; |

1673 | return(LS); |

1674 | } |

1675 | if (i==0) |

1676 | { |

1677 | "i==0 situation"; |

1678 | // no minuses at all => pluses only |

1679 | cnt++; |

1680 | LS[cnt] = string(s[1..j-1]); |

1681 | s = s[j+1..size(s)]; |

1682 | while (s!= "") |

1683 | { |

1684 | j = find(s,"+"); |

1685 | cnt++; |

1686 | if (j==0) |

1687 | { |

1688 | LS[cnt] = string(s); |

1689 | s = ""; |

1690 | } |

1691 | else |

1692 | { |

1693 | LS[cnt] = string(s[1..j-1]); |

1694 | s = s[j+1..size(s)]; |

1695 | } |

1696 | } |

1697 | return(LS); |

1698 | } |

1699 | if (j==0) |

1700 | { |

1701 | "j==0 situation"; |

1702 | // no pluses at all except the lead coef => the rest are minuses only |

1703 | cnt++; |

1704 | LS[cnt] = string(s[1..i-1]); |

1705 | s = s[i..size(s)]; |

1706 | while (s!= "") |

1707 | { |

1708 | i = find(s,"-",2); |

1709 | cnt++; |

1710 | if (i==0) |

1711 | { |

1712 | LS[cnt] = string(s); |

1713 | s = ""; |

1714 | } |

1715 | else |

1716 | { |

1717 | LS[cnt] = string(s[1..i-1]); |

1718 | s = s[i..size(s)]; |

1719 | } |

1720 | } |

1721 | return(LS); |

1722 | } |

1723 | // now i, j are nonzero |

1724 | if (i>j) |

1725 | { |

1726 | "i>j situation"; |

1727 | // + comes first, at place j |

1728 | cnt++; |

1729 | // "cnt:"; cnt; "j:"; j; |

1730 | LS[cnt] = string(s[1..j-1]); |

1731 | s = s[j+1..size(s)]; |

1732 | LL = stringpoly2lplace(s); |

1733 | LS = LS + LL; |

1734 | kill LL; |

1735 | return(LS); |

1736 | } |

1737 | else |

1738 | { |

1739 | "j>i situation"; |

1740 | // - might come first, at place i |

1741 | if (i>1) |

1742 | { |

1743 | cnt++; |

1744 | LS[cnt] = string(s[1..i-1]); |

1745 | s = s[i..size(s)]; |

1746 | } |

1747 | else |

1748 | { |

1749 | // i==1-> minus at leadcoef |

1750 | ie = find(s,"-",i+1); |

1751 | je = find(s,"+",i+1); |

1752 | if (je == ie) |

1753 | { |

1754 | "ie=je situation"; |

1755 | //monomial |

1756 | cnt++; |

1757 | LS[cnt] = s; |

1758 | return(LS); |

1759 | } |

1760 | if (je < ie) |

1761 | { |

1762 | "je<ie situation"; |

1763 | // + comes first |

1764 | cnt++; |

1765 | LS[cnt] = s[1..je-1]; |

1766 | s = s[je+1..size(s)]; |

1767 | } |

1768 | else |

1769 | { |

1770 | // ie < je |

1771 | "ie<je situation"; |

1772 | cnt++; |

1773 | LS[cnt] = s[1..ie-1]; |

1774 | s = s[ie..size(s)]; |

1775 | } |

1776 | } |

1777 | "going into recursion with "+s; |

1778 | LL = stringpoly2lplace(s); |

1779 | LS = LS + LL; |

1780 | return(LS); |

1781 | } |

1782 | } |

1783 | example |

1784 | { |

1785 | "EXAMPLE:"; echo = 2; |

1786 | string s = "x*y+y*z+z*t"; // + only |

1787 | stringpoly2lplace(s); |

1788 | string s2 = "x*y - y*z-z*t*w*w"; // +1, - only |

1789 | stringpoly2lplace(s2); |

1790 | string s3 = "-x*y + y - 2*x +7*w*w*w"; |

1791 | stringpoly2lplace(s3); |

1792 | } |

1793 | |

1794 | proc addplaces(list L) |

1795 | { |

1796 | // adds places to the list of strings |

1797 | // according to their order in the list |

1798 | int s = size(L); |

1799 | int i; |

1800 | for (i=1; i<=s; i++) |

1801 | { |

1802 | if (typeof(L[i]) == "string") |

1803 | { |

1804 | L[i] = L[i] + "(" + string(i) + ")"; |

1805 | } |

1806 | else |

1807 | { |

1808 | ERROR("entry of type string expected"); |

1809 | return(0); |

1810 | } |

1811 | } |

1812 | return(L); |

1813 | } |

1814 | example |

1815 | { |

1816 | "EXAMPLE:"; echo = 2; |

1817 | string a = "f1"; string b = "f2"; |

1818 | list L = a,b,a; |

1819 | addplaces(L); |

1820 | } |

1821 | |

1822 | proc sent2lplace(string s) |

1823 | { |

1824 | // SENTence of words TO LetterPLACE |

1825 | list L = stringpoly2lplace(s); |

1826 | int i; int ss = size(L); |

1827 | for(i=1; i<=ss; i++) |

1828 | { |

1829 | L[i] = str2lplace(L[i]); |

1830 | } |

1831 | return(L); |

1832 | } |

1833 | example |

1834 | { |

1835 | "EXAMPLE:"; echo = 2; |

1836 | ring r = 0,(f2,f1),dp; |

1837 | string s = "f2*f1*f1 - 2*f1*f2*f1+ f1*f1*f2"; |

1838 | sent2lplace(s); |

1839 | } |

1840 | |

1841 | proc testnumber(string s) |

1842 | { |

1843 | string t; |

1844 | if (s[1]=="-") |

1845 | { |

1846 | // two situations: either there's a negative number |

1847 | t = s[2..size(s)]; |

1848 | if (testnumber(t)) |

1849 | { |

1850 | //a negative number |

1851 | } |

1852 | else |

1853 | { |

1854 | // a variable times -1 |

1855 | } |

1856 | // or just a "-" for -1 |

1857 | } |

1858 | t = "ring @r=("; |

1859 | t = t + charstr(basering)+"),"; |

1860 | t = t + string(var(1))+",dp;"; |

1861 | // write(":w tstnum.tst",t); |

1862 | t = t+ "number @@Nn = " + s + ";"+"$"; |

1863 | write(":w tstnum.tst",t); |

1864 | string runsing = system("Singular"); |

1865 | int k; |

1866 | t = runsing+ " -teq <tstnum.tst >tstnum.out"; |

1867 | k = system("sh",t); |

1868 | if (k!=0) |

1869 | { |

1870 | ERROR("Problems running Singular"); |

1871 | } |

1872 | int i = system("sh", "grep error tstnum.out > /dev/NULL"); |

1873 | if (i!=0) |

1874 | { |

1875 | // no error: s is a number |

1876 | i = 1; |

1877 | } |

1878 | k = system("sh","rm tstnum.tst tstnum.out > /dev/NULL"); |

1879 | return(i); |

1880 | } |

1881 | example |

1882 | { |

1883 | "EXAMPLE:"; echo = 2; |

1884 | ring r = (0,a),x,dp; |

1885 | string s = "a^2+7*a-2"; |

1886 | testnumber(s); |

1887 | s = "b+a"; |

1888 | testnumber(s); |

1889 | } |

1890 | |

1891 | proc str2lplace(string s) |

1892 | { |

1893 | // converts a word (monomial) with coeff into letter-place |

1894 | // string: coef*var1^exp1*var2^exp2*...varN^expN |

1895 | s = strpower2rep(s); // expand powers |

1896 | if (size(s)==0) { return(0); } |

1897 | int i,j,k,insC; |

1898 | string a,b,c,d,t; |

1899 | // 1. get coeff |

1900 | i = find(s,"*"); |

1901 | if (i==0) { return(s); } |

1902 | list VN; |

1903 | c = s[1..i-1]; // incl. the case like (-a^2+1) |

1904 | int tn = testnumber(c); |

1905 | if (tn == 0) |

1906 | { |

1907 | // failed test |

1908 | if (c[1]=="-") |

1909 | { |

1910 | // two situations: either there's a negative number |

1911 | t = c[2..size(c)]; |

1912 | if (testnumber(t)) |

1913 | { |

1914 | //a negative number |

1915 | // nop here |

1916 | } |

1917 | else |

1918 | { |

1919 | // a variable times -1 |

1920 | c = "-1"; |

1921 | j++; VN[j] = t; //string(c[2..size(c)]); |

1922 | insC = 1; |

1923 | } |

1924 | } |

1925 | else |

1926 | { |

1927 | // just a variable with coeff 1 |

1928 | j++; VN[j] = string(c); |

1929 | c = "1"; |

1930 | insC = 1; |

1931 | } |

1932 | } |

1933 | // get vars |

1934 | t = s; |

1935 | // t = s[i+1..size(s)]; |

1936 | k = size(t); //j = 0; |

1937 | while (k>0) |

1938 | { |

1939 | t = t[i+1..size(t)]; //next part |

1940 | i = find(t,"*"); // next * |

1941 | if (i==0) |

1942 | { |

1943 | // last monomial |

1944 | j++; |

1945 | VN[j] = t; |

1946 | k = size(t); |

1947 | break; |

1948 | } |

1949 | b = t[1..i-1]; |

1950 | // print(b); |

1951 | j++; |

1952 | VN[j] = b; |

1953 | k = size(t); |

1954 | } |

1955 | VN = addplaces(VN); |

1956 | VN[size(VN)+1] = string(c); |

1957 | return(VN); |

1958 | } |

1959 | example |

1960 | { |

1961 | "EXAMPLE:"; echo = 2; |

1962 | ring r = (0,a),(f2,f1),dp; |

1963 | str2lplace("-2*f2^2*f1^2*f2"); |

1964 | str2lplace("-f1*f2"); |

1965 | str2lplace("(-a^2+7a)*f1*f2"); |

1966 | } |

1967 | |

1968 | proc strpower2rep(string s) |

1969 | { |

1970 | // makes x*x*x*x out of x^4 ., rep statys for repetitions |

1971 | // looks for "-" problem |

1972 | // exception: "-" as coeff |

1973 | string ex,t; |

1974 | int i,j,k; |

1975 | |

1976 | i = find(s,"^"); // first ^ |

1977 | if (i==0) { return(s); } // no ^ signs |

1978 | |

1979 | if (s[1] == "-") |

1980 | { |

1981 | // either -coef or -1 |

1982 | // got the coeff: |

1983 | i = find(s,"*"); |

1984 | if (i==0) |

1985 | { |

1986 | // no *'s => coef == -1 or s == -23 |

1987 | i = size(s)+1; |

1988 | } |

1989 | t = string(s[2..i-1]); // without "-" |

1990 | if ( testnumber(t) ) |

1991 | { |

1992 | // a good number |

1993 | t = strpower2rep(string(s[2..size(s)])); |

1994 | t = "-"+t; |

1995 | return(t); |

1996 | } |

1997 | else |

1998 | { |

1999 | // a variable |

2000 | t = strpower2rep(string(s[2..size(s)])); |

2001 | t = "-1*"+ t; |

2002 | return(t); |

2003 | } |

2004 | } |

2005 | // the case when leadcoef is a number in () |

2006 | if (s[1] == "(") |

2007 | { |

2008 | i = find(s,")",2); // must be nonzero |

2009 | t = s[2..i-1]; |

2010 | if ( testnumber(t) ) |

2011 | { |

2012 | // a good number |

2013 | } |

2014 | else {"strpower2rep: bad number as coef";} |

2015 | ex = string(s[i+2..size(s)]); // 2 because of * |

2016 | ex = strpower2rep(ex); |

2017 | t = "("+t+")*"+ex; |

2018 | return(t); |

2019 | } |

2020 | |

2021 | i = find(s,"^"); // first ^ |

2022 | j = find(s,"*",i+1); // next * == end of ^ |

2023 | if (j==0) |

2024 | { |

2025 | ex = s[i+1..size(s)]; |

2026 | } |

2027 | else |

2028 | { |

2029 | ex = s[i+1..j-1]; |

2030 | } |

2031 | execute("int @exp = " + ex + ";"); //@exp = exponent |

2032 | // got varname |

2033 | for (k=i-1; k>0; k--) |

2034 | { |

2035 | if (s[k] == "*") break; |

2036 | } |

2037 | string varn = s[k+1..i-1]; |

2038 | // "varn:"; varn; |

2039 | string pref; |

2040 | if (k>0) |

2041 | { |

2042 | pref = s[1..k]; // with * on the k-th place |

2043 | } |

2044 | // "pref:"; pref; |

2045 | string suf; |

2046 | if ( (j>0) && (j+1 <= size(s)) ) |

2047 | { |

2048 | suf = s[j+1..size(s)]; // without * on the 1st place |

2049 | } |

2050 | // "suf:"; suf; |

2051 | string toins; |

2052 | for (k=1; k<=@exp; k++) |

2053 | { |

2054 | toins = toins + varn+"*"; |

2055 | } |

2056 | // "toins: "; toins; |

2057 | if (size(suf) == 0) |

2058 | { |

2059 | toins = toins[1..size(toins)-1]; // get rid of trailing * |

2060 | } |

2061 | else |

2062 | { |

2063 | suf = strpower2rep(suf); |

2064 | } |

2065 | ex = pref + toins + suf; |

2066 | return(ex); |

2067 | // return(strpower2rep(ex)); |

2068 | } |

2069 | example |

2070 | { |

2071 | "EXAMPLE:"; echo = 2; |

2072 | ring r = (0,a),(x,y,z,t),dp; |

2073 | strpower2rep("-x^4"); |

2074 | strpower2rep("-2*x^4*y^3*z*t^2"); |

2075 | strpower2rep("-a^2*x^4"); |

2076 | } |

2077 | |

2078 | proc Liebr(poly a, poly b, list #) |

2079 | { |

2080 | // alias ppLiebr; |

2081 | //if int N is given compute [a,[...[a,b]]]] left normed bracket |

2082 | poly q; |

2083 | int N=1; |

2084 | if (size(#)>0) |

2085 | { |

2086 | if (typeof(#[1])=="int") |

2087 | { |

2088 | N = int(#[1]); |

2089 | } |

2090 | } |

2091 | if (N<=0) { return(q); } |

2092 | while (b!=0) |

2093 | { |

2094 | q = q + pmLiebr(a,lead(b)); |

2095 | b = b - lead(b); |

2096 | } |

2097 | int i; |

2098 | if (N >1) |

2099 | { |

2100 | for(i=1; i<=N; i++) |

2101 | { |

2102 | q = Liebr(a,q); |

2103 | } |

2104 | } |

2105 | return(q); |

2106 | } |

2107 | example |

2108 | { |

2109 | "EXAMPLE:"; echo = 2; |

2110 | ring r = 0,(x(1),y(1),x(2),y(2),x(3),y(3),x(4),y(4)),dp; |

2111 | poly a = x(1)*y(2); poly b = y(1); |

2112 | int uptodeg=4; int lV=2; |

2113 | export uptodeg; export lV; |

2114 | Liebr(a,b); |

2115 | Liebr(x(1),y(1),2); |

2116 | } |

2117 | |

2118 | proc pmLiebr(poly a, poly b) |

2119 | { |

2120 | // a poly, b mono |

2121 | poly s; |

2122 | while (a!=0) |

2123 | { |

2124 | s = s + mmLiebr(lead(a),lead(b)); |

2125 | a = a - lead(a); |

2126 | } |

2127 | return(s); |

2128 | } |

2129 | |

2130 | //proc pshift(poly a, int i, int uptodeg, int lV) |

2131 | proc pshift(poly a, int i) |

2132 | { |

2133 | // shifts a monomial a by i |

2134 | // calls pLPshift(p,sh,uptodeg,lVblock); |

2135 | return(system("stest",a,i,uptodeg,lV)); |

2136 | } |

2137 | |

2138 | proc mmLiebr(poly a, poly b) |

2139 | { |

2140 | // a,b, monomials |

2141 | a = lead(a); |

2142 | b = lead(b); |

2143 | int sa = deg(a); |

2144 | int sb = deg(b); |

2145 | poly v = a*pshift(b,sa) - b*pshift(a,sb); |

2146 | return(v); |

2147 | } |

2148 | |

2149 | static proc test_shift() |

2150 | { |

2151 | LIB "freegb.lib"; |

2152 | ring r = 0,(a,b),dp; |

2153 | int d =5; |

2154 | def R = freegbRing(d); |

2155 | setring R; |

2156 | int uptodeg = d; export uptodeg; |

2157 | int lV = 2; export lV; |

2158 | poly p = mmLiebr(a(1),b(1)); |

2159 | poly p = Liebr(a(1),b(1)); |

2160 | } |

2161 | |

2162 | proc Serre(intmat A, int zu) |

2163 | { |

2164 | // zu = 1 -> with commutators [f_i,f_j]; zu == 0 without them |

2165 | // suppose that A is cartan matrix |

2166 | // then Serre's relations are |

2167 | // (ad f_j)^{1-A_{ij}} ( f_i) |

2168 | int ppl = printlevel-voice+2; |

2169 | int n = ncols(A); // hence n variables |

2170 | int i,j,k,l; |

2171 | poly p,q; |

2172 | ideal I; |

2173 | for (i=1; i<=n; i++) |

2174 | { |

2175 | for (j=1; j<=n; j++) |

2176 | { |

2177 | l = 1 - A[i,j]; |

2178 | // printf("i:%s, j: %s, l: %s",i,j,l); |

2179 | dbprint(ppl,"i, j, l: ",i,j,l); |

2180 | // if ((i!=j) && (l >0)) |

2181 | // if ( (i!=j) && ( ((zu ==0) && (l >=2)) || ((zu ==1) && (l >=1)) ) ) |

2182 | if ((i!=j) && (l >0)) |

2183 | { |

2184 | q = Liebr(var(j),var(i)); |

2185 | // printf("first bracket: %s",q); |

2186 | dbprint(ppl,"first bracket: ",q); |

2187 | // if (l >=2) |

2188 | // { |

2189 | for (k=1; k<=l-1; k++) |

2190 | { |

2191 | q = Liebr(var(j),q); |

2192 | // printf("further bracket: %s",q); |

2193 | dbprint(ppl,"further bracket:",q); |

2194 | } |

2195 | // } |

2196 | } |

2197 | if (q!=0) { I = I,q; q=0;} |

2198 | } |

2199 | } |

2200 | I = simplify(I,2); |

2201 | return(I); |

2202 | } |

2203 | example |

2204 | { |

2205 | "EXAMPLE:"; echo = 2; |

2206 | intmat A[2][2] = 2, -1, -1, 2; // sl_3 == A_2 |

2207 | ring r = 0,(f1,f2),dp; |

2208 | int uptodeg = 3; |

2209 | int lV = 2; |

2210 | export uptodeg; export lV; |

2211 | def R = freegbRing(uptodeg); |

2212 | setring R; |

2213 | ideal I = Serre(A,1); |

2214 | I; |

2215 | Serre(A,0); |

2216 | } |

2217 | |

2218 | // TODO: EXPLAIN CONVERTION RULE, ENCODING OF 'LN'! |

2219 | // THIS SEEMS TO BE THE LAST PART OF FREEGBASIS? |

2220 | // TODO: ASSUMPTIONS! WHY NO SQUARES IN ELEMENTS FROM K? |

2221 | |

2222 | proc lp2lstr(ideal K, def save) |

2223 | "USAGE: lp2lstr(K,save); K an ideal, save a ring |

2224 | RETURN: nothing (exports object LN into save) |

2225 | PURPOSE: converts letter-place ideal to list of modules |

2226 | EXAMPLE: example lp2lstr; shows examples |

2227 | " |

2228 | { |

2229 | def @R = basering; |

2230 | string err; |

2231 | int s = nvars(save); |

2232 | int i,j,k; |

2233 | // K contains vars x(1),...z(1) = images of originals |

2234 | // 5. go back to orig vars, produce strings/modules |

2235 | int sk = size(K); |

2236 | int sp, sx, a, b; |

2237 | intvec x; |

2238 | poly p,q; |

2239 | poly pn; |

2240 | // vars in 'save' |

2241 | setring save; |

2242 | module N; |

2243 | list LN; |

2244 | vector V; |

2245 | poly pn; |

2246 | // test and skip exponents >=2 |

2247 | setring @R; |

2248 | |

2249 | // TODO: YOU CAN FOLD TWO FOLLOWING LOOPS TOGETHER! |

2250 | // THE FOLLOWING CHECK/PREPROCESSING MIGHT BE EXPENSIVE! |

2251 | for(i=1; i<=sk; i++) |

2252 | { |

2253 | p = K[i]; |

2254 | while (p!=0) |

2255 | { |

2256 | q = lead(p); |

2257 | // "processing q:";q; |

2258 | x = leadexp(q); |

2259 | sx = size(x); |

2260 | for(k=1; k<=sx; k++) |

2261 | { |

2262 | if ( x[k] >= 2 ) |

2263 | { |

2264 | err = "skip: the value x[k] is " + string(x[k]); |

2265 | dbprint(ppl,err); |

2266 | // return(0); |

2267 | K[i] = 0; |

2268 | p = 0; |

2269 | q = 0; |

2270 | break; |

2271 | } |

2272 | } |

2273 | p = p - q; |

2274 | } |

2275 | } |

2276 | K = simplify(K,2); |

2277 | sk = size(K); |

2278 | for(i=1; i<=sk; i++) |

2279 | { |

2280 | // setring save; |

2281 | // V = 0; |

2282 | setring @R; |

2283 | p = K[i]; |

2284 | while (p!=0) |

2285 | { |

2286 | q = lead(p); |

2287 | err = "processing q:" + string(q); |

2288 | dbprint(ppl,err); |

2289 | x = leadexp(q); |

2290 | sx = size(x); |

2291 | pn = leadcoef(q); |

2292 | setring save; |

2293 | pn = imap(@R,pn); |

2294 | V = V + leadcoef(pn)*gen(1); |

2295 | for(k=1; k<=sx; k++) |

2296 | { |

2297 | if (x[k] ==1) |

2298 | { |

2299 | a = k / s; // block number=a+1, a!=0 |

2300 | b = k % s; // remainder |

2301 | // printf("a: %s, b: %s",a,b); |

2302 | if (b == 0) |

2303 | { |

2304 | // that is it's the last var in the block |

2305 | b = s; |

2306 | a = a-1; |

2307 | } |

2308 | V = V + var(b)*gen(a+2); |

2309 | } |

2310 | } |

2311 | err = "V: " + string(V); |

2312 | dbprint(ppl,err); |

2313 | // printf("V: %s", string(V)); |

2314 | N = N,V; |

2315 | V = 0; |

2316 | setring @R; |

2317 | p = p - q; |

2318 | pn = 0; |

2319 | } |

2320 | setring save; |

2321 | LN[i] = simplify(N,2); |

2322 | N = 0; |

2323 | } |

2324 | setring save; |

2325 | export LN; |

2326 | // return(LN); |

2327 | } |

2328 | example |

2329 | { |

2330 | "EXAMPLE:"; echo = 2; |

2331 | intmat A[2][2] = 2, -1, -1, 2; // sl_3 == A_2 |

2332 | ring r = 0,(f1,f2),dp; |

2333 | int uptodeg = 3; |

2334 | def R = freegbRing(uptodeg); |

2335 | setring R; |

2336 | ideal I = Serre(A,1); // TODO: BUG: WHY IS THIS FAILED??? HIDING ARGUMENTS IS A BAD STYLE!!! |

2337 | |

2338 | int lV = 2; // NOT CLEAR WHERE DO YOU NEED THIS? |

2339 | export uptodeg; export lV; // TODO: BUG: THEY WHERE NOR DEFINED BEFORE IN CODE!!! |

2340 | I = Serre(A,1); |

2341 | |

2342 | lp2lstr(I,r); |

2343 | setring r; |

2344 | lst2str(LN,1); |

2345 | |

2346 | |

2347 | kill uptodeg; |

2348 | kill lV; |

2349 | } |

2350 | |

2351 | proc strList2poly(list L) |

2352 | { |

2353 | // list L comes from sent2lplace (which takes a poly on the input) |

2354 | // each entry of L is a sublist with the coef on the last place |

2355 | int s = size(L); int t; |

2356 | int i,j; |

2357 | list M; |

2358 | poly p,q; |

2359 | string Q; |

2360 | for(i=1; i<=s; i++) |

2361 | { |

2362 | M = L[i]; |

2363 | t = size(M); |

2364 | // q = M[t]; // a constant |

2365 | Q = string(M[t]); |

2366 | for(j=1; j<t; j++) |

2367 | { |

2368 | // q = q*M[j]; |

2369 | Q = Q+"*"+string(M[j]); |

2370 | } |

2371 | execute("q="+Q+";"); |

2372 | // q; |

2373 | p = p + q; |

2374 | } |

2375 | kill Q; |

2376 | return(p); |

2377 | } |

2378 | example |

2379 | { |

2380 | "EXAMPLE:"; echo = 2; |

2381 | ring r =0,(x,y,z,t),Dp; |

2382 | def A = freegbRing(4); |

2383 | setring A; |

2384 | string t = "-2*y*z*y*z + y*t*z*z - z*x*x*y + 2*z*y*z*y"; |

2385 | list L = sent2lplace(t); |

2386 | L; |

2387 | poly p = strList2poly(L); |

2388 | p; |

2389 | } |

2390 | |

2391 | proc file2lplace(string fname) |

2392 | { |

2393 | // format: from the usual string to letterplace |

2394 | string s = read(fname); |

2395 | // assume: file is a comma-sep list of polys |

2396 | // the vars are declared before |

2397 | // the file ends with ";" |

2398 | string t; int i; |

2399 | ideal I; |

2400 | list tst; |

2401 | while (s!="") |

2402 | { |

2403 | i = find(s,","); |

2404 | "i"; i; |

2405 | if (i==0) |

2406 | { |

2407 | i = find(s,";"); |

2408 | if (i==0) |

2409 | { |

2410 | // no ; ?? |

2411 | "no colon or semicolon found anymore"; |

2412 | return(I); |

2413 | } |

2414 | // no "," but ";" on the i-th place |

2415 | t = s[1..i-1]; |

2416 | s = ""; |

2417 | "processing: "; t; |

2418 | tst = sent2lplace(t); |

2419 | tst; |

2420 | I = I, strList2poly(tst); |

2421 | return(I); |

2422 | } |

2423 | // here i !=0 |

2424 | t = s[1..i-1]; |

2425 | s = s[i+1..size(s)]; |

2426 | "processing: "; t; |

2427 | tst = sent2lplace(t); |

2428 | tst; |

2429 | I = I, strList2poly(tst); |

2430 | } |

2431 | return(I); |

2432 | } |

2433 | example |

2434 | { |

2435 | "EXAMPLE:"; echo = 2; |

2436 | ring r =0,(x,y,z,t),dp; |

2437 | def A = freegbRing(4); |

2438 | setring A; |

2439 | string fn = "myfile"; |

2440 | string s1 = "z*y*y*y - 3*y*z*x*y + 3*y*y*z*y - y*x*y*z,"; |

2441 | string s2 = "-2*y*x*y*z + y*y*z*z - z*z*y*y + 2*z*y*z*y,"; |

2442 | string s3 = "z*y*x*t - 2*y*z*y*t + y*y*z*t - t*z*y*y + 2*t*y*z*y - t*x*y*z;"; |

2443 | write(":w "+fn,s1); write(":a "+fn,s2); write(":a "+fn,s3); |

2444 | read(fn); |

2445 | ideal I = file2lplace(fn); |

2446 | I; |

2447 | } |

2448 | |

2449 | static proc get_ls3nilp() |

2450 | { |

2451 | //first app of file2lplace |

2452 | ring r =0,(x,y,z,t),dp; |

2453 | int d = 10; |

2454 | def A = freegbRing(d); |

2455 | setring A; |

2456 | ideal I = file2lplace("./ls3nilp.bg"); |

2457 | // and now test the correctness: go back from lplace to strings |

2458 | lp2lstr(I,r); |

2459 | setring r; |

2460 | lst2str(LN,1); // agree! |

2461 | } |

2462 | |

2463 | static proc doc_example |

2464 | { |

2465 | LIB "freegb.lib"; |

2466 | ring r = 0,(x,y,z),dp; |

2467 | int d =4; // degree bound |

2468 | def R = freegbRing(d); |

2469 | setring R; |

2470 | ideal I = x(1)*y(2) + y(1)*z(2), x(1)*x(2) + x(1)*y(2) - y(1)*x(2) - y(1)*y(2); |

2471 | option(redSB);option(redTail); |

2472 | ideal J = system("freegb",I,d,nvars(r)); |

2473 | J; |

2474 | // visualization: |

2475 | lp2lstr(J,r); // export an object called LN to the ring r |

2476 | setring r; // change to the ring r |

2477 | lst2str(LN,1); // output the strings |

2478 | } |

2479 | |

2480 | |

2481 | |

2482 | // TODO: |

2483 | // multiply two letterplace polynomials, lpMult |

2484 | // reduction/ Normalform? needs kernel stuff |

2485 |