Changeset 2c7f28 in git
- Timestamp:
- May 27, 2011, 4:25:28 PM (13 years ago)
- Branches:
- (u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
- Children:
- 600ea94eb3125257cff2947d86f88bde3f272148
- Parents:
- 9bb545791780a99344a26d6177decb3870f6da87
- git-author:
- Frank Seelisch <seelisch@mathematik.uni-kl.de>2011-05-27 16:25:28+02:00
- git-committer:
- Mohamed Barakat <mohamed.barakat@rwth-aachen.de>2011-11-09 12:39:06+01:00
- Location:
- libpolys/polys
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
libpolys/polys/Makefile.am
r9bb5457 r2c7f28 49 49 kbuckets.cc sbuckets.cc weight.cc weight0.c simpleideals.cc matpol.cc \ 50 50 ${USE_P_PROCS_STATIC_CC} ${USE_P_PROCS_DYNAMIC_CC} templates/mod_raw.cc \ 51 ext_fields/algext.cc clapsing.cc clapconv.cc51 ext_fields/algext.cc ext_fields/transext.cc clapsing.cc clapconv.cc 52 52 53 53 BUILT_SOURCES = templates/p_Procs.inc -
libpolys/polys/ext_fields/algext.cc
r9bb5457 r2c7f28 88 88 if (a == NULL) return TRUE; 89 89 p_Test((poly)a, naRing); 90 assume(p_ Deg((poly)a, naRing) <= p_Deg(naMinpoly, naRing));90 assume(p_Totaldegree((poly)a, naRing) <= p_Totaldegree(naMinpoly, naRing)); 91 91 return TRUE; 92 92 } … … 128 128 } 129 129 130 BOOLEAN naEqual 130 BOOLEAN naEqual(number a, number b, const coeffs cf) 131 131 { 132 132 naTest(a); naTest(b); … … 139 139 /// deg test 140 140 int aDeg = 0; 141 if (a != NULL) aDeg = p_ Deg((poly)a, naRing);141 if (a != NULL) aDeg = p_Totaldegree((poly)a, naRing); 142 142 int bDeg = 0; 143 if (b != NULL) bDeg = p_ Deg((poly)b, naRing);143 if (b != NULL) bDeg = p_Totaldegree((poly)b, naRing); 144 144 if (aDeg != bDeg) return FALSE; 145 145 … … 172 172 { 173 173 naTest(a); 174 if (p_GetExp((poly)a, 1, naRing) != 0) return FALSE; 175 return n_IsOne(p_GetCoeff((poly)a, naRing), naCoeffs); 174 poly aAsPoly = (poly)a; 175 if (!p_IsConstant(aAsPoly, naRing)) return FALSE; 176 return n_IsOne(p_GetCoeff(aAsPoly, naRing), naCoeffs); 176 177 } 177 178 … … 179 180 { 180 181 naTest(a); 181 if (p_GetExp((poly)a, 1, naRing) != 0) return FALSE; 182 return n_IsMOne(p_GetCoeff((poly)a, naRing), naCoeffs); 182 poly aAsPoly = (poly)a; 183 if (!p_IsConstant(aAsPoly, naRing)) return FALSE; 184 return n_IsMOne(p_GetCoeff(aAsPoly, naRing), naCoeffs); 183 185 } 184 186 … … 206 208 { 207 209 naTest(a); 208 if (p_GetExp((poly)a, 1, naRing) != 0) return 0; 209 return n_Int(p_GetCoeff((poly)a, naRing), naCoeffs); 210 poly aAsPoly = (poly)a; 211 if (!p_IsConstant(aAsPoly, naRing)) return 0; 212 return n_Int(p_GetCoeff(aAsPoly, naRing), naCoeffs); 210 213 } 211 214 … … 213 216 BOOLEAN naGreater(number a, number b, const coeffs cf) 214 217 { 218 naTest(a); naTest(b); 215 219 if (naIsZero(a, cf)) return FALSE; 216 220 if (naIsZero(b, cf)) return TRUE; 217 221 int aDeg = 0; 218 if (a != NULL) aDeg = p_ Deg((poly)a, naRing);222 if (a != NULL) aDeg = p_Totaldegree((poly)a, naRing); 219 223 int bDeg = 0; 220 if (b != NULL) bDeg = p_ Deg((poly)b, naRing);224 if (b != NULL) bDeg = p_Totaldegree((poly)b, naRing); 221 225 return (aDeg > bDeg); 222 226 } … … 228 232 if (a == NULL) return FALSE; 229 233 if (n_GreaterZero(p_GetCoeff((poly)a, naRing), naCoeffs)) return TRUE; 230 if (p_ Deg((poly)a, naRing) > 0)return TRUE;234 if (p_Totaldegree((poly)a, naRing) > 0) return TRUE; 231 235 return FALSE; 232 236 } … … 298 302 for |exp| >= 8 compute power along binary presentation of |exp|, e.g. 299 303 p^13 = p^1 * p^4 * p^8, where we utilise that 300 p^(2^(k+1)) = p^(2^k) * p^(2^k) 304 p^(2^(k+1)) = p^(2^k) * p^(2^k); 301 305 intermediate reduction modulo the minimal polynomial is controlled by 302 306 the in-place method heuristicReduce(poly, poly, coeffs); see there. … … 318 322 int expAbs = exp; if (expAbs < 0) expAbs = -expAbs; 319 323 320 /* now compute 'a' to the 'expAbs'-th power*/324 /* now compute a^expAbs */ 321 325 poly pow; poly aAsPoly = (poly)a; 322 326 if (expAbs <= 7) … … 363 367 } 364 368 365 /* may reduce p modul ethe reducer by calling definiteReduce;369 /* may reduce p modulo the reducer by calling definiteReduce; 366 370 the decision is made based on the following heuristic 367 371 (which should also only be changed here in this method): … … 374 378 p_Test((poly)reducer, naRing); 375 379 #endif 376 if (p_ Deg(p, naRing) > 10 * p_Deg(reducer, naRing))380 if (p_Totaldegree(p, naRing) > 10 * p_Totaldegree(reducer, naRing)) 377 381 definiteReduce(p, reducer, cf); 378 382 } … … 424 428 AlgExtInfo *e = (AlgExtInfo *)param; 425 429 /* for extension coefficient fields we expect the underlying 426 polynomial srings to be IDENTICAL, i.e. the SAME OBJECT;430 polynomial rings to be IDENTICAL, i.e. the SAME OBJECT; 427 431 this expectation is based on the assumption that we have properly 428 432 registered cf and perform reference counting rather than creating … … 445 449 { 446 450 noOfTerms++; 447 int d = 0; 448 for (int i = 1; i <= rVar(naRing); i++) 449 d += p_GetExp(aAsPoly, i, naRing); 451 int d = p_GetExp(aAsPoly, 1, naRing); 450 452 if (d > theDegree) theDegree = d; 451 453 pIter(aAsPoly); … … 496 498 number naMap00(number a, const coeffs src, const coeffs dst) 497 499 { 500 if (n_IsZero(a, src)) return NULL; 498 501 assume(src == dst->extRing->cf); 499 502 poly result = p_One(dst->extRing); … … 505 508 number naMapP0(number a, const coeffs src, const coeffs dst) 506 509 { 510 if (n_IsZero(a, src)) return NULL; 507 511 /* mapping via intermediate int: */ 508 512 int n = n_Int(a, src); … … 514 518 515 519 /* assumes that either src = Q(a), dst = Q(a), or 516 src = Z p(a), dst = Zp(a) */520 src = Z/p(a), dst = Z/p(a) */ 517 521 number naCopyMap(number a, const coeffs src, const coeffs dst) 518 522 { … … 523 527 number naMap0P(number a, const coeffs src, const coeffs dst) 524 528 { 529 if (n_IsZero(a, src)) return NULL; 525 530 int p = rChar(dst->extRing); 526 531 int n = nlModP(a, p, src); … … 534 539 number naMapPP(number a, const coeffs src, const coeffs dst) 535 540 { 541 if (n_IsZero(a, src)) return NULL; 536 542 assume(src == dst->extRing->cf); 537 543 poly result = p_One(dst->extRing); … … 543 549 number naMapUP(number a, const coeffs src, const coeffs dst) 544 550 { 551 if (n_IsZero(a, src)) return NULL; 545 552 /* mapping via intermediate int: */ 546 553 int n = n_Int(a, src); … … 554 561 { 555 562 /* dst is expected to be an algebraic field extension */ 556 assume(getCoeffType(dst) == n _algExt);563 assume(getCoeffType(dst) == naID); 557 564 558 565 int h = 0; /* the height of the extension tower given by dst */ … … 585 592 if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst)) 586 593 { 587 if (strcmp(r Parameter(src->extRing)[0],588 r Parameter(dst->extRing)[0]) == 0)594 if (strcmp(rRingVar(0, src->extRing), 595 rRingVar(0, dst->extRing)) == 0) 589 596 return naCopyMap; /// Q(a) --> Q(a) 590 597 else … … 605 612 606 613 BOOLEAN naInitChar(coeffs cf, void * infoStruct) 607 { 608 assume( getCoeffType(cf) == naID ); 609 614 { 610 615 AlgExtInfo *e = (AlgExtInfo *)infoStruct; 611 616 /// first check whether cf->extRing != NULL and delete old ring??? -
libpolys/polys/ext_fields/transext.cc
r9bb5457 r2c7f28 4 4 /* $Id$ */ 5 5 /* 6 * ABSTRACT: numbers in a transcendental extension field K(t_1, .., t_s) 7 * Assuming that we have a coeffs object cf, then these numbers 8 * are quotients of polynomials in the polynomial ring 9 * K[t_1, .., t_s] represented by cf->algring. 6 * ABSTRACT: numbers in a rational function field K(t_1, .., t_s) with 7 * transcendental variables t_1, ..., t_s, where s >= 1. 8 * Denoting the implemented coeffs object by cf, then these numbers 9 * are represented as quotients of polynomials in the polynomial 10 * ring K[t_1, .., t_s] represented by cf->algring. 10 11 */ 11 12 … … 55 56 number ntGcd(number a, number b, const coeffs cf); 56 57 number ntLcm(number a, number b, const coeffs cf); 57 numberntSize(number a, const coeffs cf);58 int ntSize(number a, const coeffs cf); 58 59 void ntDelete(number * a, const coeffs cf); 59 60 void ntCoeffWrite(const coeffs cf); … … 62 63 static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param); 63 64 65 void heuristicGcdCancellation(number a, const coeffs cf); 66 void definiteGcdCancellation(number a, const coeffs cf); 67 64 68 #ifdef LDEBUG 65 69 BOOLEAN ntDBTest(number a, const char *f, const int l, const coeffs cf) 66 70 { 67 71 assume(getCoeffType(cf) == ntID); 68 fraction f= (fraction)a;69 if (is0( f)) return TRUE;70 assume(num( f) != NULL); /**< f != 0 ==> numerator(f) != 0 */71 p_Test(num( f), ntRing);72 if (!denIs1( f)) p_Test(den(f), ntRing);72 fraction t = (fraction)a; 73 if (is0(t)) return TRUE; 74 assume(num(t) != NULL); /**< t != 0 ==> numerator(t) != 0 */ 75 p_Test(num(t), ntRing); 76 if (!denIs1(t)) p_Test(den(t), ntRing); 73 77 return TRUE; 74 78 } 75 79 #endif 76 77 void heuristicReduce(poly &p, poly reducer, const coeffs cf);78 void definiteReduce(poly &p, poly reducer, const coeffs cf);79 80 80 81 /* returns the bottom field in this field extension tower; if the tower … … 96 97 } 97 98 98 BOOLEAN naIsZero(number a, const coeffs cf) 99 { 100 naTest(a); 101 return (a == NULL); 102 } 103 104 void naDelete(number * a, const coeffs cf) 105 { 106 if (*a == NULL) return; 107 poly aAsPoly = (poly)(*a); 108 p_Delete(&aAsPoly, naRing); 99 BOOLEAN ntIsZero(number a, const coeffs cf) 100 { 101 ntTest(a); 102 return (is0(a)); 103 } 104 105 void ntDelete(number * a, const coeffs cf) 106 { 107 fraction f = (fraction)(*a); 108 if (is0(f)) return; 109 p_Delete(&num(f), ntRing); 110 if (!denIs1(f)) p_Delete(&den(f), ntRing); 111 omFreeBin((ADDRESS)f, fractionObjectBin); 109 112 *a = NULL; 110 113 } 111 114 112 BOOLEAN n aEqual(number a, number b, const coeffs cf)113 { 114 n aTest(a); naTest(b);115 BOOLEAN ntEqual(number a, number b, const coeffs cf) 116 { 117 ntTest(a); ntTest(b); 115 118 116 119 /// simple tests 117 120 if (a == b) return TRUE; 118 if ((a == NULL) && (b != NULL)) return FALSE; 119 if ((b == NULL) && (a != NULL)) return FALSE; 120 121 /// deg test 122 int aDeg = 0; 123 if (a != NULL) aDeg = p_Deg((poly)a, naRing); 124 int bDeg = 0; 125 if (b != NULL) bDeg = p_Deg((poly)b, naRing); 126 if (aDeg != bDeg) return FALSE; 127 128 /// subtraction test 129 number c = naSub(a, b, cf); 130 BOOLEAN result = naIsZero(c, cf); 131 naDelete(&c, naCoeffs); 132 return result; 133 } 134 135 number naCopy(number a, const coeffs cf) 136 { 137 naTest(a); 138 if (a == NULL) return NULL; 139 return (number)p_Copy((poly)a, naRing); 140 } 141 142 number naGetNumerator(number &a, const coeffs cf) 143 { 144 return naCopy(a, cf); 145 } 146 147 number naGetDenom(number &a, const coeffs cf) 148 { 149 naTest(a); 150 return naInit(1, cf); 151 } 152 153 BOOLEAN naIsOne(number a, const coeffs cf) 154 { 155 naTest(a); 156 if (p_GetExp((poly)a, 1, naRing) != 0) return FALSE; 157 return n_IsOne(p_GetCoeff((poly)a, naRing), naCoeffs); 158 } 159 160 BOOLEAN naIsMOne(number a, const coeffs cf) 161 { 162 naTest(a); 163 if (p_GetExp((poly)a, 1, naRing) != 0) return FALSE; 164 return n_IsMOne(p_GetCoeff((poly)a, naRing), naCoeffs); 121 if ((is0(a)) && (!is0(b))) return FALSE; 122 if ((is0(b)) && (!is0(a))) return FALSE; 123 124 /// cheap test if gcd's have been cancelled in both numbers 125 fraction fa = (fraction)a; 126 fraction fb = (fraction)b; 127 if ((c(fa) == 1) && (c(fb) == 1)) 128 { 129 poly f = p_Add_q(p_Copy(num(fa), ntRing), 130 p_Neg(p_Copy(num(fb), ntRing), ntRing), 131 ntRing); 132 if (f != NULL) { p_Delete(&f, ntRing); return FALSE; } 133 if (denIs1(fa) && denIs1(fb)) return TRUE; 134 if (denIs1(fa) && !denIs1(fb)) return FALSE; 135 if (!denIs1(fa) && denIs1(fb)) return FALSE; 136 f = p_Add_q(p_Copy(den(fa), ntRing), 137 p_Neg(p_Copy(den(fb), ntRing), ntRing), 138 ntRing); 139 if (f != NULL) { p_Delete(&f, ntRing); return FALSE; } 140 return TRUE; 141 } 142 143 /* default: the more expensive multiplication test 144 a/b = c/d <==> a*d = b*c */ 145 poly f = p_Copy(num(fa), ntRing); 146 if (!denIs1(fb)) f = p_Mult_q(f, p_Copy(den(fb), ntRing), ntRing); 147 poly g = p_Copy(num(fb), ntRing); 148 if (!denIs1(fa)) g = p_Mult_q(g, p_Copy(den(fa), ntRing), ntRing); 149 poly h = p_Add_q(f, p_Neg(g, ntRing), ntRing); 150 if (h == NULL) return TRUE; 151 else 152 { 153 p_Delete(&h, ntRing); 154 return FALSE; 155 } 156 } 157 158 number ntCopy(number a, const coeffs cf) 159 { 160 ntTest(a); 161 if (is0(a)) return NULL; 162 fraction f = (fraction)a; 163 poly g = p_Copy(num(f), ntRing); 164 poly h = NULL; if (!denIs1(f)) h = p_Copy(den(f), ntRing); 165 fraction result = (fraction)omAlloc0Bin(fractionObjectBin); 166 num(result) = g; 167 den(result) = h; 168 c(result) = c(f); 169 return (number)result; 170 } 171 172 number ntGetNumerator(number &a, const coeffs cf) 173 { 174 ntTest(a); 175 definiteGcdCancellation(a, cf); 176 if (is0(a)) return NULL; 177 fraction f = (fraction)a; 178 poly g = p_Copy(num(f), ntRing); 179 fraction result = (fraction)omAlloc0Bin(fractionObjectBin); 180 num(result) = g; 181 den(result) = NULL; 182 c(result) = 0; 183 return (number)result; 184 } 185 186 number ntGetDenom(number &a, const coeffs cf) 187 { 188 ntTest(a); 189 definiteGcdCancellation(a, cf); 190 fraction f = (fraction)a; 191 poly g; 192 if (is0(f) || denIs1(f)) g = p_One(ntRing); 193 else g = p_Copy(den(f), ntRing); 194 fraction result = (fraction)omAlloc0Bin(fractionObjectBin); 195 num(result) = g; 196 den(result) = NULL; 197 c(result) = 0; 198 return (number)result; 199 } 200 201 BOOLEAN ntIsOne(number a, const coeffs cf) 202 { 203 ntTest(a); 204 definiteGcdCancellation(a, cf); 205 fraction f = (fraction)a; 206 if (!denIs1(f)) return FALSE; 207 poly g = num(f); 208 if (!p_IsConstant(g, ntRing)) return FALSE; 209 return n_IsOne(p_GetCoeff(g, ntRing), ntCoeffs); 210 } 211 212 BOOLEAN ntIsMOne(number a, const coeffs cf) 213 { 214 ntTest(a); 215 definiteGcdCancellation(a, cf); 216 fraction f = (fraction)a; 217 if (!denIs1(f)) return FALSE; 218 poly g = num(f); 219 if (!p_IsConstant(g, ntRing)) return FALSE; 220 return n_IsMOne(p_GetCoeff(g, ntRing), ntCoeffs); 165 221 } 166 222 167 223 /// this is in-place, modifies a 168 number naNeg(number a, const coeffs cf) 169 { 170 naTest(a); 171 if (a != NULL) a = (number)p_Neg((poly)a, naRing); 224 number ntNeg(number a, const coeffs cf) 225 { 226 ntTest(a); 227 if (!is0(a)) 228 { 229 fraction f = (fraction)a; 230 num(f) = p_Neg(num(f), ntRing); 231 } 172 232 return a; 173 233 } 174 234 175 number n aImPart(number a, const coeffs cf)176 { 177 n aTest(a);235 number ntImPart(number a, const coeffs cf) 236 { 237 ntTest(a); 178 238 return NULL; 179 239 } 180 240 181 number n aInit(int i, const coeffs cf)241 number ntInit(int i, const coeffs cf) 182 242 { 183 243 if (i == 0) return NULL; 184 else return (number)p_ISet(i, naRing); 185 } 186 187 int naInt(number &a, const coeffs cf) 188 { 189 naTest(a); 190 if (p_GetExp((poly)a, 1, naRing) != 0) return 0; 191 return n_Int(p_GetCoeff((poly)a, naRing), naCoeffs); 192 } 193 194 /* TRUE iff (a != 0 and (b == 0 or deg(a) > deg(b))) */ 195 BOOLEAN naGreater(number a, number b, const coeffs cf) 196 { 197 if (naIsZero(a, cf)) return FALSE; 198 if (naIsZero(b, cf)) return TRUE; 199 int aDeg = 0; 200 if (a != NULL) aDeg = p_Deg((poly)a, naRing); 201 int bDeg = 0; 202 if (b != NULL) bDeg = p_Deg((poly)b, naRing); 203 return (aDeg > bDeg); 204 } 205 206 /* TRUE iff a != 0 and (LC(a) > 0 or deg(a) > 0) */ 207 BOOLEAN naGreaterZero(number a, const coeffs cf) 208 { 209 naTest(a); 210 if (a == NULL) return FALSE; 211 if (n_GreaterZero(p_GetCoeff((poly)a, naRing), naCoeffs)) return TRUE; 212 if (p_Deg((poly)a, naRing) > 0) return TRUE; 213 return FALSE; 214 } 215 216 void naCoeffWrite(const coeffs cf) 217 { 218 char *x = rRingVar(0, naRing); 219 Print("// Coefficients live in the extension field K[%s]/<f(%s)>\n", x, x); 220 Print("// with the minimal polynomial f(%s) = %s\n", x, 221 p_String(naMinpoly, naRing)); 222 PrintS("// and K: "); n_CoeffWrite(cf->extRing->cf); 223 } 224 225 number naPar(int i, const coeffs cf) 226 { 227 assume(i == 1); // there is only one parameter in this extension field 228 poly p = p_ISet(1, naRing); // p = 1 229 p_SetExp(p, 1, 1, naRing); // p = the sole extension variable 230 p_Setm(p, naRing); 231 return (number)p; 232 } 233 234 number naAdd(number a, number b, const coeffs cf) 235 { 236 naTest(a); naTest(b); 237 if (a == NULL) return naCopy(b, cf); 238 if (b == NULL) return naCopy(a, cf); 239 poly aPlusB = p_Add_q(p_Copy((poly)a, naRing), 240 p_Copy((poly)b, naRing), naRing); 241 definiteReduce(aPlusB, naMinpoly, cf); 242 return (number)aPlusB; 243 } 244 245 number naSub(number a, number b, const coeffs cf) 246 { 247 naTest(a); naTest(b); 248 if (b == NULL) return naCopy(a, cf); 249 poly minusB = p_Neg(p_Copy((poly)b, naRing), naRing); 250 if (a == NULL) return (number)minusB; 251 poly aMinusB = p_Add_q(p_Copy((poly)a, naRing), minusB, naRing); 252 definiteReduce(aMinusB, naMinpoly, cf); 253 return (number)aMinusB; 254 } 255 256 number naMult(number a, number b, const coeffs cf) 257 { 258 naTest(a); naTest(b); 259 if (a == NULL) return NULL; 260 if (b == NULL) return NULL; 261 poly aTimesB = p_Mult_q(p_Copy((poly)a, naRing), 262 p_Copy((poly)b, naRing), naRing); 263 definiteReduce(aTimesB, naMinpoly, cf); 264 return (number)aTimesB; 265 } 266 267 number naDiv(number a, number b, const coeffs cf) 268 { 269 naTest(a); naTest(b); 270 if (b == NULL) WerrorS(nDivBy0); 271 if (a == NULL) return NULL; 272 poly bInverse = (poly)naInvers(b, cf); 273 poly aDivB = p_Mult_q(p_Copy((poly)a, naRing), bInverse, naRing); 274 definiteReduce(aDivB, naMinpoly, cf); 275 return (number)aDivB; 244 else 245 { 246 fraction result = (fraction)omAlloc0Bin(fractionObjectBin); 247 num(result) = p_ISet(i, ntRing); 248 den(result) = NULL; 249 c(result) = 0; 250 return (number)result; 251 } 252 } 253 254 int ntInt(number &a, const coeffs cf) 255 { 256 ntTest(a); 257 if (is0(a)) return 0; 258 definiteGcdCancellation(a, cf); 259 fraction f = (fraction)a; 260 if (!denIs1(f)) return 0; 261 if (!p_IsConstant(num(f), ntRing)) return 0; 262 return n_Int(p_GetCoeff(num(f), ntRing), ntCoeffs); 263 } 264 265 /* This method will only consider the numerators of a and b, without 266 cancelling gcd's before. 267 Moreover it may return TRUE only if one or both numerators 268 are zero or if their degrees are equal. Then TRUE is returned iff 269 coeff(numerator(a)) > coeff(numerator(b)); 270 In all other cases, FALSE will be returned. */ 271 BOOLEAN ntGreater(number a, number b, const coeffs cf) 272 { 273 ntTest(a); ntTest(b); 274 number aNumCoeff = NULL; int aNumDeg = 0; 275 number bNumCoeff = NULL; int bNumDeg = 0; 276 if (!is0(a)) 277 { 278 fraction fa = (fraction)a; 279 aNumDeg = p_Totaldegree(num(fa), ntRing); 280 aNumCoeff = p_GetCoeff(num(fa), ntRing); 281 } 282 if (!is0(b)) 283 { 284 fraction fb = (fraction)b; 285 bNumDeg = p_Totaldegree(num(fb), ntRing); 286 bNumCoeff = p_GetCoeff(num(fb), ntRing); 287 } 288 if (aNumDeg != bNumDeg) return FALSE; 289 else return n_Greater(aNumCoeff, bNumCoeff, ntCoeffs); 290 } 291 292 /* this method will only consider the numerator of a, without cancelling 293 the gcd before; 294 returns TRUE iff the leading coefficient of the numerator of a is > 0 295 or the leading term of the numerator of a is not a 296 constant */ 297 BOOLEAN ntGreaterZero(number a, const coeffs cf) 298 { 299 ntTest(a); 300 if (is0(a)) return FALSE; 301 fraction f = (fraction)a; 302 poly g = num(f); 303 return (n_GreaterZero(p_GetCoeff(g, ntRing), ntCoeffs) || 304 (!p_LmIsConstant(g, ntRing))); 305 } 306 307 void ntCoeffWrite(const coeffs cf) 308 { 309 PrintS("// Coefficients live in the rational function field\n"); 310 Print("// K("); 311 for (int i = 0; i < rVar(ntRing); i++) 312 { 313 if (i > 0) PrintS(", "); 314 Print("%s", rRingVar(i, ntRing)); 315 } 316 PrintS(") with\n"); 317 PrintS("// K: "); n_CoeffWrite(cf->extRing->cf); 318 } 319 320 /* the i-th parameter */ 321 number ntPar(int i, const coeffs cf) 322 { 323 assume((1 <= i) && (i <= rVar(ntRing))); 324 poly p = p_ISet(1, ntRing); 325 p_SetExp(p, i, 1, ntRing); 326 p_Setm(p, ntRing); 327 fraction result = (fraction)omAlloc0Bin(fractionObjectBin); 328 num(result) = p; 329 den(result) = NULL; 330 c(result) = 0; 331 return (number)result; 332 } 333 334 number ntAdd(number a, number b, const coeffs cf) 335 { 336 ntTest(a); ntTest(b); 337 if (is0(a)) return ntCopy(b, cf); 338 if (is0(b)) return ntCopy(a, cf); 339 340 fraction fa = (fraction)a; 341 fraction fb = (fraction)b; 342 poly f; 343 if (denIs1(fa) && denIs1(fb)) f = NULL; 344 else if (!denIs1(fa) && denIs1(fb)) f = p_Copy(den(fa), ntRing); 345 else if (denIs1(fa) && !denIs1(fb)) f = p_Copy(den(fb), ntRing); 346 else /* both denom's are != 1 */ f = p_Mult_q(p_Copy(den(fa), ntRing), 347 p_Copy(den(fb), ntRing), 348 ntRing); 349 poly g = p_Copy(num(fa), ntRing); 350 if (!denIs1(fb)) g = p_Mult_q(g, p_Copy(den(fb), ntRing), ntRing); 351 poly h = p_Copy(num(fb), ntRing); 352 if (!denIs1(fa)) h = p_Mult_q(h, p_Copy(den(fa), ntRing), ntRing); 353 g = p_Add_q(g, h, ntRing); 354 fraction result = (fraction)omAlloc0Bin(fractionObjectBin); 355 num(result) = g; 356 den(result) = f; 357 c(result) = c(fa) + c(fb) + ADD_COMPLEXITY; 358 heuristicGcdCancellation((number)result, cf); 359 return (number)result; 360 } 361 362 number ntSub(number a, number b, const coeffs cf) 363 { 364 ntTest(a); ntTest(b); 365 if (is0(a)) return ntNeg(ntCopy(b, cf), cf); 366 if (is0(b)) return ntCopy(a, cf); 367 368 fraction fa = (fraction)a; 369 fraction fb = (fraction)b; 370 poly f; 371 if (denIs1(fa) && denIs1(fb)) f = NULL; 372 else if (!denIs1(fa) && denIs1(fb)) f = p_Copy(den(fa), ntRing); 373 else if (denIs1(fa) && !denIs1(fb)) f = p_Copy(den(fb), ntRing); 374 else /* both den's are != 1 */ f = p_Mult_q(p_Copy(den(fa), ntRing), 375 p_Copy(den(fb), ntRing), 376 ntRing); 377 poly g = p_Copy(num(fa), ntRing); 378 if (!denIs1(fb)) g = p_Mult_q(g, p_Copy(den(fb), ntRing), ntRing); 379 poly h = p_Copy(num(fb), ntRing); 380 if (!denIs1(fa)) h = p_Mult_q(h, p_Copy(den(fa), ntRing), ntRing); 381 g = p_Add_q(g, p_Neg(h, ntRing), ntRing); 382 fraction result = (fraction)omAlloc0Bin(fractionObjectBin); 383 num(result) = g; 384 den(result) = f; 385 c(result) = c(fa) + c(fb) + ADD_COMPLEXITY; 386 heuristicGcdCancellation((number)result, cf); 387 return (number)result; 388 } 389 390 number ntMult(number a, number b, const coeffs cf) 391 { 392 ntTest(a); ntTest(b); 393 if (is0(a) || is0(b)) return NULL; 394 395 fraction fa = (fraction)a; 396 fraction fb = (fraction)b; 397 poly f; 398 if (denIs1(fa) && denIs1(fb)) f = NULL; 399 else if (!denIs1(fa) && denIs1(fb)) f = p_Copy(den(fa), ntRing); 400 else if (denIs1(fa) && !denIs1(fb)) f = p_Copy(den(fb), ntRing); 401 else /* both den's are != 1 */ f = p_Mult_q(p_Copy(den(fa), ntRing), 402 p_Copy(den(fb), ntRing), 403 ntRing); 404 poly g = p_Copy(num(fa), ntRing); 405 poly h = p_Copy(num(fb), ntRing); 406 g = p_Mult_q(g, h, ntRing); 407 fraction result = (fraction)omAlloc0Bin(fractionObjectBin); 408 num(result) = g; 409 den(result) = f; 410 c(result) = c(fa) + c(fb) + MULT_COMPLEXITY; 411 heuristicGcdCancellation((number)result, cf); 412 return (number)result; 413 } 414 415 number ntDiv(number a, number b, const coeffs cf) 416 { 417 ntTest(a); ntTest(b); 418 if (is0(a)) return NULL; 419 if (is0(b)) WerrorS(nDivBy0); 420 421 fraction fa = (fraction)a; 422 fraction fb = (fraction)b; 423 poly f = p_Copy(num(fb), ntRing); 424 if (!denIs1(fa)) f = p_Mult_q(f, p_Copy(den(fa), ntRing), ntRing); 425 poly g = p_Copy(num(fa), ntRing); 426 if (!denIs1(fb)) g = p_Mult_q(g, p_Copy(den(fb), ntRing), ntRing); 427 fraction result = (fraction)omAlloc0Bin(fractionObjectBin); 428 num(result) = g; 429 den(result) = f; 430 c(result) = c(fa) + c(fb) + MULT_COMPLEXITY; 431 heuristicGcdCancellation((number)result, cf); 432 return (number)result; 276 433 } 277 434 … … 280 437 for |exp| >= 8 compute power along binary presentation of |exp|, e.g. 281 438 p^13 = p^1 * p^4 * p^8, where we utilise that 282 p^(2^(k+1)) = p^(2^k) * p^(2^k) 283 intermediate reduction modulo the minimal polynomial is controlled by284 the in-place method heuristicReduce(poly, poly, coeffs); see there.439 p^(2^(k+1)) = p^(2^k) * p^(2^k); 440 intermediate cancellation is controlled by the in-place method 441 heuristicGcdCancellation; see there. 285 442 */ 286 void n aPower(number a, int exp, number *b, const coeffs cf)287 { 288 n aTest(a);443 void ntPower(number a, int exp, number *b, const coeffs cf) 444 { 445 ntTest(a); 289 446 290 447 /* special cases first */ 291 if ( a == NULL)448 if (is0(a)) 292 449 { 293 450 if (exp >= 0) *b = NULL; 294 451 else WerrorS(nDivBy0); 295 452 } 296 else if (exp == 0) *b = n aInit(1, cf);297 else if (exp == 1) *b = n aCopy(a, cf);298 else if (exp == -1) *b = n aInvers(a, cf);453 else if (exp == 0) *b = ntInit(1, cf); 454 else if (exp == 1) *b = ntCopy(a, cf); 455 else if (exp == -1) *b = ntInvers(a, cf); 299 456 300 457 int expAbs = exp; if (expAbs < 0) expAbs = -expAbs; 301 458 302 /* now compute 'a' to the 'expAbs'-th power*/303 poly pow; poly aAsPoly = (poly)a;459 /* now compute a^expAbs */ 460 number pow; number t; 304 461 if (expAbs <= 7) 305 462 { 306 pow = p_Copy(aAsPoly, naRing);463 pow = ntCopy(a, cf); 307 464 for (int i = 2; i <= expAbs; i++) 308 465 { 309 pow = p_Mult_q(pow, p_Copy(aAsPoly, naRing), naRing); 310 heuristicReduce(pow, naMinpoly, cf); 466 t = ntMult(pow, a, cf); 467 ntDelete(&pow, cf); 468 pow = t; 469 heuristicGcdCancellation(pow, cf); 311 470 } 312 definiteReduce(pow, naMinpoly, cf);313 471 } 314 472 else 315 473 { 316 pow = p_ISet(1, naRing);317 poly factor = p_Copy(aAsPoly, naRing);474 pow = ntInit(1, cf); 475 number factor = ntCopy(a, cf); 318 476 while (expAbs != 0) 319 477 { 320 478 if (expAbs & 1) 321 479 { 322 pow = p_Mult_q(pow, p_Copy(factor, naRing), naRing); 323 heuristicReduce(pow, naMinpoly, cf); 480 t = ntMult(pow, factor, cf); 481 ntDelete(&pow, cf); 482 pow = t; 483 heuristicGcdCancellation(pow, cf); 324 484 } 325 485 expAbs = expAbs / 2; 326 486 if (expAbs != 0) 327 487 { 328 factor = p_Mult_q(factor, factor, naRing); 329 heuristicReduce(factor, naMinpoly, cf); 488 t = ntMult(factor, factor, cf); 489 ntDelete(&factor, cf); 490 factor = t; 491 heuristicGcdCancellation(factor, cf); 330 492 } 331 493 } 332 p_Delete(&factor, naRing); 333 definiteReduce(pow, naMinpoly, cf); 494 ntDelete(&factor, cf); 334 495 } 335 496 336 497 /* invert if original exponent was negative */ 337 number n = (number)pow;338 498 if (exp < 0) 339 499 { 340 number m = naInvers(n, cf); 341 naDelete(&n, cf); 342 n = m; 343 } 344 *b = n; 345 } 346 347 /* may reduce p module the reducer by calling definiteReduce; 348 the decision is made based on the following heuristic 349 (which should also only be changed here in this method): 350 if (deg(p) > 10*deg(reducer) then perform reduction; 351 modifies p */ 352 void heuristicReduce(poly &p, poly reducer, const coeffs cf) 353 { 354 #ifdef LDEBUG 355 p_Test((poly)p, naRing); 356 p_Test((poly)reducer, naRing); 357 #endif 358 if (p_Deg(p, naRing) > 10 * p_Deg(reducer, naRing)) 359 definiteReduce(p, reducer, cf); 360 } 361 362 void naWrite(number &a, const coeffs cf) 363 { 364 naTest(a); 365 if (a == NULL) 500 t = ntInvers(pow, cf); 501 ntDelete(&pow, cf); 502 pow = t; 503 } 504 *b = pow; 505 } 506 507 /* modifies a */ 508 void heuristicGcdCancellation(number a, const coeffs cf) 509 { 510 ntTest(a); 511 if (is0(a)) return; 512 fraction f = (fraction)a; 513 if (denIs1(f) || n_IsOne(p_GetCoeff(num(f), ntRing), ntCoeffs)) 514 { c(f) = 0; return; } 515 if (c(f) <= BOUND_COMPLEXITY) return; 516 else definiteGcdCancellation(a, cf); 517 } 518 519 /* modifies a */ 520 void definiteGcdCancellation(number a, const coeffs cf) 521 { 522 ntTest(a); 523 if (is0(a)) return; 524 fraction f = (fraction)a; 525 if (denIs1(f) || n_IsOne(p_GetCoeff(num(f), ntRing), ntCoeffs)) 526 { c(f) = 0; return; } 527 if (c(f) > BOUND_COMPLEXITY) 528 { 529 /* TO BE IMPLEMENTED! 530 for the time, cancellation of gcd's does not take place */ 531 Print("// TO BE IMPLEMENTED: transext.cc:definiteGcdCancellation\n"); 532 Print("// (complexity of number = %d exceeds the bound = %d\n", 533 c(f), BOUND_COMPLEXITY); 534 } 535 } 536 537 void ntWrite(number &a, const coeffs cf) 538 { 539 ntTest(a); 540 definiteGcdCancellation(a, cf); 541 if (is0(a)) 366 542 StringAppendS("0"); 367 543 else 368 544 { 369 poly aAsPoly = (poly)a; 370 /* basically, just write aAsPoly using p_Write, 371 but use brackets around the output, if a is not 372 a constant living in naCoeffs = cf->extRing->cf */ 373 BOOLEAN useBrackets = !(p_IsConstant(aAsPoly, naRing)); 545 fraction f = (fraction)a; 546 BOOLEAN useBrackets = !(p_IsConstant(num(f), ntRing)) && (!denIs1(f)); 374 547 if (useBrackets) StringAppendS("("); 375 p_String0( aAsPoly, naRing, naRing);548 p_String0(num(f), ntRing, ntRing); 376 549 if (useBrackets) StringAppendS(")"); 377 } 378 } 379 380 const char * naRead(const char *s, number *a, const coeffs cf) 381 { 382 poly aAsPoly; 383 const char * result = p_Read(s, aAsPoly, naRing); 384 *a = (number)aAsPoly; 385 return result; 386 } 387 388 /* implemented by the rule lcm(a, b) = a * b / gcd(a, b) */ 389 number naLcm(number a, number b, const coeffs cf) 390 { 391 naTest(a); naTest(b); 392 if (a == NULL) return NULL; 393 if (b == NULL) return NULL; 394 number theProduct = (number)p_Mult_q(p_Copy((poly)a, naRing), 395 p_Copy((poly)b, naRing), naRing); 396 /* note that theProduct needs not be reduced w.r.t. naMinpoly; 397 but the final division will take care of the necessary reduction */ 398 number theGcd = naGcd(a, b, cf); 399 return naDiv(theProduct, theGcd, cf); 400 } 401 402 /* expects *param to be castable to AlgExtInfo */ 403 static BOOLEAN naCoeffIsEqual(const coeffs cf, n_coeffType n, void * param) 404 { 405 if (naID != n) return FALSE; 406 AlgExtInfo *e = (AlgExtInfo *)param; 407 /* for extension coefficient fields we expect the underlying 408 polynomials rings to be IDENTICAL, i.e. the SAME OBJECT; 550 if (!denIs1(f)) 551 { 552 StringAppendS("/"); 553 useBrackets = !p_IsConstant(den(f), ntRing); 554 if (useBrackets) StringAppendS("("); 555 p_String0(den(f), ntRing, ntRing); 556 if (useBrackets) StringAppendS(")"); 557 } 558 } 559 } 560 561 /* this just reads polynomials; 562 the rest will be taken care of by the interpreter: When it encounters 563 a "/", then it will impose a division that will finally lead to a 564 rational function */ 565 const char * ntRead(const char *s, number *a, const coeffs cf) 566 { 567 poly p; 568 const char * result = p_Read(s, p, ntRing); 569 if (p == NULL) { *a = NULL; return result; } 570 else 571 { 572 fraction f = (fraction)omAlloc0Bin(fractionObjectBin); 573 num(f) = p; 574 den(f) = NULL; 575 c(f) = 0; 576 *a = (number)f; 577 return result; 578 } 579 } 580 581 /* expects *param to be castable to TransExtInfo */ 582 static BOOLEAN ntCoeffIsEqual(const coeffs cf, n_coeffType n, void * param) 583 { 584 if (ntID != n) return FALSE; 585 TransExtInfo *e = (TransExtInfo *)param; 586 /* for rational function fields we expect the underlying 587 polynomial rings to be IDENTICAL, i.e. the SAME OBJECT; 409 588 this expectation is based on the assumption that we have properly 410 589 registered cf and perform reference counting rather than creating 411 590 multiple copies of the same coefficient field/domain/ring */ 412 return (naRing == e->r); 413 /* (Note that then also the minimal ideals will necessarily be 414 the same, as they are attached to the ring.) */ 415 } 416 417 int naSize(number a, const coeffs cf) 418 { 419 if (a == NULL) return -1; 591 return (ntRing == e->r); 592 } 593 594 number ntLcm(number a, number b, const coeffs cf) 595 { 596 ntTest(a); ntTest(b); 597 /* TO BE IMPLEMENTED! 598 for the time, we simply return NULL, representing the number zero */ 599 Print("// TO BE IMPLEMENTED: transext.cc:ntLcm\n"); 600 return NULL; 601 } 602 603 number ntGcd(number a, number b, const coeffs cf) 604 { 605 ntTest(a); ntTest(b); 606 /* TO BE IMPLEMENTED! 607 for the time, we simply return NULL, representing the number zero */ 608 Print("// TO BE IMPLEMENTED: transext.cc:ntGcd\n"); 609 return NULL; 610 } 611 612 int ntSize(number a, const coeffs cf) 613 { 614 ntTest(a); 615 if (is0(a)) return -1; 420 616 /* this has been taken from the old implementation of field extensions, 421 where we computed the sum of the degree and the number of terms in 422 (poly)a; so we leave it at that, for the time being; 423 maybe, the number of terms alone is a better measure? */ 424 poly aAsPoly = (poly)a; 425 int theDegree = 0; int noOfTerms = 0; 426 while (aAsPoly != NULL) 617 where we computed the sum of the degrees and the numbers of terms in 618 the numerator and denominator of a; so we leave it at that, for the 619 time being */ 620 fraction f = (fraction)a; 621 poly p = num(f); 622 int noOfTerms = 0; 623 int numDegree = 0; 624 while (p != NULL) 427 625 { 428 626 noOfTerms++; 429 627 int d = 0; 430 for (int i = 1; i <= rVar(naRing); i++) 431 d += p_GetExp(aAsPoly, i, naRing); 432 if (d > theDegree) theDegree = d; 433 pIter(aAsPoly); 434 } 435 return theDegree + noOfTerms; 436 } 437 438 /* performs polynomial division and overrides p by the remainder 439 of division of p by the reducer; 440 modifies p */ 441 void definiteReduce(poly &p, poly reducer, const coeffs cf) 442 { 443 #ifdef LDEBUG 444 p_Test((poly)p, naRing); 445 p_Test((poly)reducer, naRing); 446 #endif 447 p_PolyDiv(p, reducer, FALSE, naRing); 448 } 449 450 /* IMPORTANT NOTE: Since an algebraic field extension is again a field, 451 the gcd of two elements is not very interesting. (It 452 is actually any unit in the field, i.e., any non- 453 zero element.) Note that the below method does not operate 454 in this strong sense but rather computes the gcd of 455 two given elements in the underlying polynomial ring. */ 456 number naGcd(number a, number b, const coeffs cf) 457 { 458 naTest(a); naTest(b); 459 if ((a == NULL) && (b == NULL)) WerrorS(nDivBy0); 460 return (number)p_Gcd((poly)a, (poly)b, naRing); 461 } 462 463 number naInvers(number a, const coeffs cf) 464 { 465 naTest(a); 466 if (a == NULL) WerrorS(nDivBy0); 467 poly aFactor = NULL; poly mFactor = NULL; 468 poly theGcd = p_ExtGcd((poly)a, aFactor, naMinpoly, mFactor, naRing); 469 naTest((number)theGcd); naTest((number)aFactor); naTest((number)mFactor); 470 /* the gcd must be 1 since naMinpoly is irreducible and a != NULL: */ 471 assume(naIsOne((number)theGcd, cf)); 472 p_Delete(&theGcd, naRing); 473 p_Delete(&mFactor, naRing); 474 return (number)(aFactor); 475 } 476 477 /* assumes that src = Q, dst = Q(a) */ 478 number naMap00(number a, const coeffs src, const coeffs dst) 479 { 628 for (int i = 1; i <= rVar(ntRing); i++) 629 d += p_GetExp(p, i, ntRing); 630 if (d > numDegree) numDegree = d; 631 pIter(p); 632 } 633 int denDegree = 0; 634 if (!denIs1(f)) 635 { 636 p = den(f); 637 while (p != NULL) 638 { 639 noOfTerms++; 640 int d = 0; 641 for (int i = 1; i <= rVar(ntRing); i++) 642 d += p_GetExp(p, i, ntRing); 643 if (d > denDegree) denDegree = d; 644 pIter(p); 645 } 646 } 647 return numDegree + denDegree + noOfTerms; 648 } 649 650 number ntInvers(number a, const coeffs cf) 651 { 652 ntTest(a); 653 if (is0(a)) WerrorS(nDivBy0); 654 fraction f = (fraction)a; 655 fraction result = (fraction)omAlloc0Bin(fractionObjectBin); 656 poly g; 657 if (denIs1(f)) g = p_One(ntRing); 658 else g = p_Copy(den(f), ntRing); 659 num(result) = g; 660 den(result) = p_Copy(num(f), ntRing); 661 c(result) = c(f); 662 return (number)result; 663 } 664 665 /* assumes that src = Q, dst = Q(t_1, ..., t_s) */ 666 number ntMap00(number a, const coeffs src, const coeffs dst) 667 { 668 if (n_IsZero(a, src)) return NULL; 480 669 assume(src == dst->extRing->cf); 481 poly result = p_One(dst->extRing); 482 p_SetCoeff(result, naCopy(a, src), dst->extRing); 483 return (number)result; 484 } 485 486 /* assumes that src = Z/p, dst = Q(a) */ 487 number naMapP0(number a, const coeffs src, const coeffs dst) 488 { 670 poly p = p_One(dst->extRing); 671 p_SetCoeff(p, ntCopy(a, src), dst->extRing); 672 fraction f = (fraction)omAlloc0Bin(fractionObjectBin); 673 num(f) = p; den(f) = NULL; c(f) = 0; 674 return (number)f; 675 } 676 677 /* assumes that src = Z/p, dst = Q(t_1, ..., t_s) */ 678 number ntMapP0(number a, const coeffs src, const coeffs dst) 679 { 680 if (n_IsZero(a, src)) return NULL; 489 681 /* mapping via intermediate int: */ 490 682 int n = n_Int(a, src); 491 683 number q = n_Init(n, dst->extRing->cf); 492 poly result = p_One(dst->extRing); 493 p_SetCoeff(result, q, dst->extRing); 494 return (number)result; 495 } 496 497 /* assumes that either src = Q(a), dst = Q(a), or 498 src = Zp(a), dst = Zp(a) */ 499 number naCopyMap(number a, const coeffs src, const coeffs dst) 500 { 501 return naCopy(a, dst); 502 } 503 504 /* assumes that src = Q, dst = Z/p(a) */ 505 number naMap0P(number a, const coeffs src, const coeffs dst) 506 { 684 poly p; 685 if (n_IsZero(q, dst->extRing->cf)) 686 { 687 n_Delete(&q, dst->extRing->cf); 688 return NULL; 689 } 690 p = p_One(dst->extRing); 691 p_SetCoeff(p, q, dst->extRing); 692 fraction f = (fraction)omAlloc0Bin(fractionObjectBin); 693 num(f) = p; den(f) = NULL; c(f) = 0; 694 return (number)f; 695 } 696 697 /* assumes that either src = Q(t_1, ..., t_s), dst = Q(t_1, ..., t_s), or 698 src = Z/p(t_1, ..., t_s), dst = Z/p(t_1, ..., t_s) */ 699 number ntCopyMap(number a, const coeffs src, const coeffs dst) 700 { 701 return ntCopy(a, dst); 702 } 703 704 /* assumes that src = Q, dst = Z/p(t_1, ..., t_s) */ 705 number ntMap0P(number a, const coeffs src, const coeffs dst) 706 { 707 if (n_IsZero(a, src)) return NULL; 507 708 int p = rChar(dst->extRing); 508 709 int n = nlModP(a, p, src); 509 710 number q = n_Init(n, dst->extRing->cf); 510 poly result = p_One(dst->extRing); 511 p_SetCoeff(result, q, dst->extRing); 512 return (number)result; 513 } 514 515 /* assumes that src = Z/p, dst = Z/p(a) */ 516 number naMapPP(number a, const coeffs src, const coeffs dst) 517 { 711 poly g; 712 if (n_IsZero(q, dst->extRing->cf)) 713 { 714 n_Delete(&q, dst->extRing->cf); 715 return NULL; 716 } 717 g = p_One(dst->extRing); 718 p_SetCoeff(g, q, dst->extRing); 719 fraction f = (fraction)omAlloc0Bin(fractionObjectBin); 720 num(f) = g; den(f) = NULL; c(f) = 0; 721 return (number)f; 722 } 723 724 /* assumes that src = Z/p, dst = Z/p(t_1, ..., t_s) */ 725 number ntMapPP(number a, const coeffs src, const coeffs dst) 726 { 727 if (n_IsZero(a, src)) return NULL; 518 728 assume(src == dst->extRing->cf); 519 poly result = p_One(dst->extRing); 520 p_SetCoeff(result, naCopy(a, src), dst->extRing); 521 return (number)result; 522 } 523 524 /* assumes that src = Z/u, dst = Z/p(a), where u != p */ 525 number naMapUP(number a, const coeffs src, const coeffs dst) 526 { 729 poly p = p_One(dst->extRing); 730 p_SetCoeff(p, ntCopy(a, src), dst->extRing); 731 fraction f = (fraction)omAlloc0Bin(fractionObjectBin); 732 num(f) = p; den(f) = NULL; c(f) = 0; 733 return (number)f; 734 } 735 736 /* assumes that src = Z/u, dst = Z/p(t_1, ..., t_s), where u != p */ 737 number ntMapUP(number a, const coeffs src, const coeffs dst) 738 { 739 if (n_IsZero(a, src)) return NULL; 527 740 /* mapping via intermediate int: */ 528 741 int n = n_Int(a, src); 529 742 number q = n_Init(n, dst->extRing->cf); 530 poly result = p_One(dst->extRing); 531 p_SetCoeff(result, q, dst->extRing); 532 return (number)result; 533 } 534 535 nMapFunc naSetMap(const coeffs src, const coeffs dst) 536 { 537 /* dst is expected to be an algebraic field extension */ 538 assume(getCoeffType(dst) == n_algExt); 743 poly p; 744 if (n_IsZero(q, dst->extRing->cf)) 745 { 746 n_Delete(&q, dst->extRing->cf); 747 return NULL; 748 } 749 p = p_One(dst->extRing); 750 p_SetCoeff(p, q, dst->extRing); 751 fraction f = (fraction)omAlloc0Bin(fractionObjectBin); 752 num(f) = p; den(f) = NULL; c(f) = 0; 753 return (number)f; 754 } 755 756 nMapFunc ntSetMap(const coeffs src, const coeffs dst) 757 { 758 /* dst is expected to be a rational function field */ 759 assume(getCoeffType(dst) == ntID); 539 760 540 761 int h = 0; /* the height of the extension tower given by dst */ … … 546 767 if ((!nCoeff_is_Zp(bDst)) && (!nCoeff_is_Q(bDst))) return NULL; 547 768 769 /* Let T denote the sequence of transcendental extension variables, i.e., 770 K[t_1, ..., t_s] =: K[T]; 771 Let moreover, for any such sequence T, T' denote any subsequence of T 772 of the form t_1, ..., t_w with w <= s. */ 773 548 774 if (nCoeff_is_Q(src) && nCoeff_is_Q(bDst)) 549 return n aMap00; /// Q --> Q(a)775 return ntMap00; /// Q --> Q(T) 550 776 551 777 if (nCoeff_is_Zp(src) && nCoeff_is_Q(bDst)) 552 return n aMapP0; /// Z/p --> Q(a)778 return ntMapP0; /// Z/p --> Q(T) 553 779 554 780 if (nCoeff_is_Q(src) && nCoeff_is_Zp(bDst)) 555 return n aMap0P; /// Q --> Z/p(a)781 return ntMap0P; /// Q --> Z/p(T) 556 782 557 783 if (nCoeff_is_Zp(src) && nCoeff_is_Zp(bDst)) 558 784 { 559 if (src->ch == dst->ch) return n aMapPP; /// Z/p --> Z/p(a)560 else return n aMapUP; /// Z/u --> Z/p(a)785 if (src->ch == dst->ch) return ntMapPP; /// Z/p --> Z/p(T) 786 else return ntMapUP; /// Z/u --> Z/p(T) 561 787 } 562 788 … … 567 793 if (nCoeff_is_Q(bSrc) && nCoeff_is_Q(bDst)) 568 794 { 569 if ( strcmp(rParameter(src->extRing)[0],570 rParameter(dst->extRing)[0]) == 0)571 return naCopyMap; /// Q(a) --> Q(a)572 else573 return NULL; /// Q(b) --> Q(a)795 if (rVar(src->extRing) > rVar(dst->extRing)) return NULL; 796 for (int i = 0; i < rVar(src->extRing); i++) 797 if (strcmp(rRingVar(i, src->extRing), 798 rRingVar(i, dst->extRing)) != 0) return NULL; 799 return ntCopyMap; /// Q(T') --> Q(T) 574 800 } 575 801 576 802 if (nCoeff_is_Zp(bSrc) && nCoeff_is_Zp(bDst)) 577 803 { 578 if ( strcmp(rParameter(src->extRing)[0],579 rParameter(dst->extRing)[0]) == 0)580 return naCopyMap; /// Z/p(a) --> Z/p(a)581 else582 return NULL; /// Z/p(b) --> Z/p(a)804 if (rVar(src->extRing) > rVar(dst->extRing)) return NULL; 805 for (int i = 0; i < rVar(src->extRing); i++) 806 if (strcmp(rRingVar(i, src->extRing), 807 rRingVar(i, dst->extRing)) != 0) return NULL; 808 return ntCopyMap; /// Z/p(T') --> Z/p(T) 583 809 } 584 810 … … 586 812 } 587 813 588 BOOLEAN naInitChar(coeffs cf, void * infoStruct) 589 { 590 assume( getCoeffType(cf) == naID ); 591 592 AlgExtInfo *e = (AlgExtInfo *)infoStruct; 814 BOOLEAN ntInitChar(coeffs cf, void * infoStruct) 815 { 816 TransExtInfo *e = (TransExtInfo *)infoStruct; 593 817 /// first check whether cf->extRing != NULL and delete old ring??? 594 818 cf->extRing = e->r; 595 cf->extRing->minideal = e->i; 596 597 assume(cf->extRing != NULL); // extRing; 598 assume((cf->extRing->minideal != NULL) && // minideal has one 599 (IDELEMS(cf->extRing->minideal) != 0) && // non-zero generator 600 (cf->extRing->minideal->m[0] != NULL) ); // at m[0]; 601 assume(cf->extRing->cf != NULL); // extRing->cf; 602 assume(getCoeffType(cf) == naID); // coeff type; 819 cf->extRing->minideal = NULL; 820 821 assume(cf->extRing != NULL); // extRing; 822 assume(cf->extRing->cf != NULL); // extRing->cf; 823 assume(getCoeffType(cf) == ntID); // coeff type; 603 824 604 825 /* propagate characteristic up so that it becomes … … 606 827 cf->ch = cf->extRing->cf->ch; 607 828 608 #ifdef LDEBUG 609 p_Test((poly)naMinpoly, naRing); 610 #endif 611 612 cf->cfGreaterZero = naGreaterZero; 613 cf->cfGreater = naGreater; 614 cf->cfEqual = naEqual; 615 cf->cfIsZero = naIsZero; 616 cf->cfIsOne = naIsOne; 617 cf->cfIsMOne = naIsMOne; 618 cf->cfInit = naInit; 619 cf->cfInt = naInt; 620 cf->cfNeg = naNeg; 621 cf->cfPar = naPar; 622 cf->cfAdd = naAdd; 623 cf->cfSub = naSub; 624 cf->cfMult = naMult; 625 cf->cfDiv = naDiv; 626 cf->cfExactDiv = naDiv; 627 cf->cfPower = naPower; 628 cf->cfCopy = naCopy; 629 cf->cfWrite = naWrite; 630 cf->cfRead = naRead; 631 cf->cfDelete = naDelete; 632 cf->cfSetMap = naSetMap; 633 cf->cfGetDenom = naGetDenom; 634 cf->cfGetNumerator = naGetNumerator; 635 cf->cfRePart = naCopy; 636 cf->cfImPart = naImPart; 637 cf->cfCoeffWrite = naCoeffWrite; 638 cf->cfDBTest = naDBTest; 639 cf->cfGcd = naGcd; 640 cf->cfLcm = naLcm; 641 cf->cfSize = naSize; 642 cf->nCoeffIsEqual = naCoeffIsEqual; 643 cf->cfInvers = naInvers; 644 cf->cfIntDiv = naDiv; 829 cf->cfGreaterZero = ntGreaterZero; 830 cf->cfGreater = ntGreater; 831 cf->cfEqual = ntEqual; 832 cf->cfIsZero = ntIsZero; 833 cf->cfIsOne = ntIsOne; 834 cf->cfIsMOne = ntIsMOne; 835 cf->cfInit = ntInit; 836 cf->cfInt = ntInt; 837 cf->cfNeg = ntNeg; 838 cf->cfPar = ntPar; 839 cf->cfAdd = ntAdd; 840 cf->cfSub = ntSub; 841 cf->cfMult = ntMult; 842 cf->cfDiv = ntDiv; 843 cf->cfExactDiv = ntDiv; 844 cf->cfPower = ntPower; 845 cf->cfCopy = ntCopy; 846 cf->cfWrite = ntWrite; 847 cf->cfRead = ntRead; 848 cf->cfDelete = ntDelete; 849 cf->cfSetMap = ntSetMap; 850 cf->cfGetDenom = ntGetDenom; 851 cf->cfGetNumerator = ntGetNumerator; 852 cf->cfRePart = ntCopy; 853 cf->cfImPart = ntImPart; 854 cf->cfCoeffWrite = ntCoeffWrite; 855 cf->cfDBTest = ntDBTest; 856 cf->cfGcd = ntGcd; 857 cf->cfLcm = ntLcm; 858 cf->cfSize = ntSize; 859 cf->nCoeffIsEqual = ntCoeffIsEqual; 860 cf->cfInvers = ntInvers; 861 cf->cfIntDiv = ntDiv; 645 862 646 863 return FALSE; -
libpolys/polys/ext_fields/transext.h
r9bb5457 r2c7f28 6 6 /* $Id$ */ 7 7 /* 8 * ABSTRACT: numbers in a transcendental extension field K(t_1, .., t_s) 9 * Assuming that we have a coeffs object cf, then these numbers 10 * are quotients of polynomials in the polynomial ring 11 * K[t_1, .., t_s] represented by cf->algring. 8 * ABSTRACT: numbers in a rational function field K(t_1, .., t_s) with 9 * transcendental variables t_1, ..., t_s, where s >= 1. 10 * Denoting the implemented coeffs object by cf, then these numbers 11 * are represented as quotients of polynomials in the polynomial 12 * ring K[t_1, .., t_s] represented by cf->algring. 12 13 */ 13 14 … … 20 21 typedef struct sip_sideal * ideal; 21 22 23 struct spolyrec; 24 typedef struct spolyrec polyrec; 25 typedef polyrec * poly; 26 22 27 /// struct for passing initialization parameters to naInitChar 23 typedef struct { ring r; } transExtInfo;28 typedef struct { ring r; } TransExtInfo; 24 29 25 30 /* a number in K(t_1, .., t_s) is represented by either NULL … … 28 33 if the denominator is 1, the member 'denominator' is NULL; 29 34 as a consequence of the above we get: if some number n is not NULL, then 30 n->numerator cannot be NULL */ 31 struct { poly numerator; poly denominator; } ip_fraction; 32 typedef struct ip_fraction * fraction; 33 /* some useful accessors: */ 34 #define is0(n) (n == NULL) /**< TRUE iff n represents 0 in K(t_1, .., t_s) */ 35 #define num(n) n->numerator 36 #define den(n) n->denominator 37 #define denIs1(n) (n->denominator == NULL) /**< TRUE iff the denom. is 1 */ 35 n->numerator cannot be NULL; 36 The member 'complexity' attempts to capture the complexity of any given 37 number n, i.e., starting with a bunch of numbers n_i that have their gcd's 38 cancelled out, n may be constructed from the n_i's by using field 39 arithmetics (+, -, *, /). If we never cancel out gcd's during this process, 40 n will become rather complex. The larger the attribute 'complexity' of n 41 is, the more likely it is that n contains some non-trivial gcd. Thus, this 42 attribute will be used by a heuristic method to cancel out gcd's from time 43 to time. (This heuristic may be set up such that cancellation can be 44 enforced after each arithmetic operation, or such that it will never take 45 place.) Moreover, the 'complexity' of n is zero iff the gcd in n (that is, 46 the gcd of its numerator and denominator) is trivial. */ 47 struct fractionObject 48 { 49 poly numerator; 50 poly denominator; 51 int complexity; 52 }; 53 typedef struct fractionObject * fraction; 54 55 omBin fractionObjectBin = omGetSpecBin(sizeof(fractionObject)); 56 57 /* constants for controlling the complexity of numbers */ 58 #define ADD_COMPLEXITY 1 /**< complexity increase due to + and - */ 59 #define MULT_COMPLEXITY 2 /**< complexity increase due to * and / */ 60 #define BOUND_COMPLEXITY 10 /**< maximum complexity of a number */ 61 62 /* some useful accessors for fractions: */ 63 #define is0(f) (f == NULL) /**< TRUE iff n represents 0 in K(t_1, .., t_s) */ 64 #define num(f) f->numerator 65 #define den(f) f->denominator 66 #define denIs1(f) (f->denominator == NULL) /**< TRUE iff den. represents 1 */ 67 #define c(f) f->complexity 38 68 39 69 /// Get a mapping function from src into the domain of this type (n_transExt) … … 68 98 number ntGcd(number a, number b, const coeffs cf); 69 99 number ntLcm(number a, number b, const coeffs cf); 70 numberntSize(number a, const coeffs cf);100 int ntSize(number a, const coeffs cf); 71 101 void ntDelete(number * a, const coeffs cf); 72 102 void ntCoeffWrite(const coeffs cf);
Note: See TracChangeset
for help on using the changeset viewer.