[BACK]Return to poly2.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / kan96xx / Kan

Annotation of OpenXM/src/kan96xx/Kan/poly2.c, Revision 1.6

1.6     ! takayama    1: /* $OpenXM: OpenXM/src/kan96xx/Kan/poly2.c,v 1.5 2003/08/21 02:30:23 takayama Exp $ */
1.1       maekawa     2: #include <stdio.h>
                      3: #include "datatype.h"
                      4: #include "stackm.h"
                      5: #include "extern.h"
                      6: #include "extern2.h"
                      7:
                      8: static POLY mapZmonom(POLY f,struct ring *ringp);
                      9:
                     10: POLY ppAdd(f,g)
1.3       takayama   11:         POLY f; POLY g;  /* The result is read only. */
1.1       maekawa    12: {
                     13:   POLY node;
                     14:   struct listPoly nod;
                     15:   POLY h;
                     16:   struct coeff *c;
                     17:   int gt;
                     18:
                     19:   node = &nod;
                     20:   node->next = POLYNULL;
                     21:   h = node;
                     22:   if (f == POLYNULL) return(g);
                     23:   if (g == POLYNULL) return(f);
                     24:   checkRing(f,g);
                     25:
                     26:   while (f != POLYNULL && g != POLYNULL) {
                     27:     /*printf("%s + %s\n",POLYToString(f,'*',1),POLYToString(g,'*',1));*/
                     28:     checkRing2(f,g); /* for debug */
                     29:     gt = (*mmLarger)(f,g);
                     30:     switch (gt) {
                     31:     case 1: /* f > g */
                     32:       h -> next = newCell(f->coeffp,f->m);
                     33:       h = h->next;
                     34:       f = f->next;
                     35:       if (f == POLYNULL) {
1.3       takayama   36:                h->next = g;
                     37:                return(node->next);
1.1       maekawa    38:       }
                     39:       break;
                     40:     case 0: /* f < g */
                     41:       h->next = newCell(g->coeffp,g->m);
                     42:       h = h->next;
                     43:       g = g->next;
                     44:       if (g == POLYNULL) {
1.3       takayama   45:                h->next = f;
                     46:                return(node->next);
1.1       maekawa    47:       }
                     48:       break;
                     49:     case 2:/* f == g */
                     50:       c = coeffCopy(f->coeffp);
                     51:       Cadd(c,c,g->coeffp);
                     52:       if (!isZero(c)) {
1.3       takayama   53:                h->next = newCell(c,f->m);
                     54:                h = h->next;
                     55:                f = f->next;
                     56:                g = g->next;
                     57:                if (f == POLYNULL) {
                     58:                  h->next = g;
                     59:                  return(node->next);
                     60:                }
                     61:                if (g == POLYNULL) {
                     62:                  h->next = f;
                     63:                  return(node->next);
                     64:                }
1.1       maekawa    65:       }else{
1.3       takayama   66:                f = f->next;
                     67:                g = g->next;
                     68:                if (f == POLYNULL) {
                     69:                  h->next = g;
                     70:                  return(node->next);
                     71:                }
                     72:                if (g == POLYNULL) {
                     73:                  h->next = f;
                     74:                  return(node->next);
                     75:                }
1.1       maekawa    76:       }
                     77:       break;
                     78:     default:
                     79:       errorPoly("ppAdd(). Internal error. Invalid value of mmLarger.");
                     80:       break;
                     81:     }
                     82:   }
                     83:   return(node->next);
                     84: }
                     85:
                     86: POLY ppSub(f,g)
1.3       takayama   87:         POLY f; POLY g;  /* The result is read only. */
1.1       maekawa    88: {
                     89:   POLY h;
                     90:   struct coeff *c;
                     91:
                     92:   if (g == POLYNULL) return(f);
                     93:   if (f == POLYNULL) return(cpMult(intToCoeff(-1,g->m->ringp),g));
                     94:   checkRing(f,g);
                     95:
                     96:   h = cpMult(intToCoeff(-1,g->m->ringp),g);
                     97:   return(ppAdd(f,h));
                     98: }
                     99:
                    100:
                    101: POLY cpMult(c,f)
1.3       takayama  102:         struct coeff *c;
                    103:         POLY f;
1.1       maekawa   104: {
                    105:   POLY node;
                    106:   struct listPoly nod;
                    107:   POLY h;
                    108:   struct coeff *newc;
                    109:   int p;
                    110:   node = &nod;
                    111:   if (f == POLYNULL || isZero(c)) return(POLYNULL);
                    112:   p = f->m->ringp->p;
                    113:   node ->next = POLYNULL;
                    114:   h = node;
                    115:   while (f != POLYNULL) {
                    116:     newc = coeffCopy(c);
                    117:     Cmult(newc,newc,f->coeffp);
                    118:     if ((p==0) || !isZero(newc)) {
                    119:       h->next = newCell(newc,f->m);
                    120:       h = h->next;
                    121:     }
                    122:     f = f->next;
                    123:   }
                    124:   return(node->next);
                    125: }
                    126:
                    127: MONOMIAL monomialAdd_poly(m,m2)
1.3       takayama  128:         MONOMIAL m,m2;
1.1       maekawa   129: {
                    130:   extern int Msize;
                    131:   MONOMIAL f;
                    132:   int i;
                    133:   int n;
                    134:   n = m->ringp->n;
                    135:   f = (MONOMIAL) sGC_malloc(sizeof(struct smallMonomial)+n*Msize);
                    136:
                    137:   if (f == (MONOMIAL) NULL) errorPoly("No more memory.");
                    138:   f->ringp = m->ringp;
                    139:   for (i=0; i<n; i++) {
                    140:     (f->e)[i].x = (m->e)[i].x + (m2->e)[i].x;
                    141:     (f->e)[i].D = (m->e)[i].D + (m2->e)[i].D;
                    142:   }
                    143:   return(f);
                    144: }
                    145:
                    146: /* Note that mpMult_poly is called from mmLarger_tower! */
                    147: POLY mpMult_poly(f,g)
1.3       takayama  148:         POLY f;
                    149:         POLY g;
1.1       maekawa   150: {
                    151:   POLY node;
                    152:   struct listPoly nod;
                    153:   struct coeff *c;
                    154:   int n,i;
                    155:   POLY h;
                    156:   MONOMIAL m;
                    157:   int p;
                    158:   node = &nod;
                    159:   if (f == POLYNULL || g == POLYNULL) return(POLYNULL);
                    160:   node->next = POLYNULL;
                    161:   h = node;
                    162:   checkRing(f,g);
                    163:   n = f->m->ringp->n; p = f->m->ringp->p;
                    164:   while(g != POLYNULL) {
                    165:     checkRing2(f,g);
                    166:     c = coeffCopy(f->coeffp);
                    167:     Cmult(c,c,g->coeffp);
                    168:     if ((p==0) || !isZero(c)) {
                    169:       m = (*monomialAdd)(f->m,g->m);
                    170:       h->next = newCell(c,m);
                    171:       h = h->next;
                    172:     }
                    173:     g = g->next;
                    174:   }
                    175:   return(node->next);
                    176: }
                    177:
                    178: POLY ppMult_old(f,g)
1.3       takayama  179:         POLY f,g;
1.1       maekawa   180: {
                    181:   POLY r;
                    182:   POLY tmp;
                    183:   r = POLYNULL;
                    184:   while( f != POLYNULL) {
                    185:     tmp = (*mpMult)(f,g);
                    186:     r = ppAddv(r,tmp); /* r and tmp will be broken */
                    187:     f = f->next;
                    188:   }
                    189:   return(r);
                    190: }
                    191:
                    192: POLY ppAddv(f,g)
1.3       takayama  193:         POLY f; POLY g;  /* It breaks f and g. Use it just after calling mpMult() */
1.1       maekawa   194: {
                    195:   POLY node;
                    196:   struct listPoly nod;
                    197:   POLY h;
                    198:   struct coeff *c;
                    199:   int gt;
                    200:
                    201:   node = &nod;
                    202:   /*printf("ppAddv:1%s + %s\n",POLYToString(f,'*',1),POLYToString(g,'*',1));*/
                    203:   node->next = POLYNULL;
                    204:   h = node;
                    205:   if (f == POLYNULL) return(g);
                    206:   if (g == POLYNULL) return(f);
                    207:   checkRing(f,g);
                    208:
                    209:   while (f != POLYNULL && g != POLYNULL) {
                    210:     checkRing2(f,g); /* for debug */
                    211:     /*printf("%s + %s\n",POLYToString(f,'*',1),POLYToString(g,'*',1));*/
                    212:     gt = (*mmLarger)(f,g);
                    213:     switch (gt) {
                    214:     case 1: /* f > g */
                    215:       h->next = f;
                    216:       h = h->next; f = f->next;;
                    217:       if (f == POLYNULL) {
1.3       takayama  218:                h->next = g;
                    219:                return(node->next);
1.1       maekawa   220:       }
                    221:       break;
                    222:     case 0: /* f < g */
                    223:       h->next =  g;
                    224:       h = h->next; g = g->next;
                    225:       if (g == POLYNULL) {
1.3       takayama  226:                h->next = f;
                    227:                return(node->next);
1.1       maekawa   228:       }
                    229:       break;
                    230:     case 2:/* f == g */
                    231:       c = f->coeffp;
                    232:       Cadd(c,c,g->coeffp);
                    233:       if (!isZero(c)) {
1.3       takayama  234:                h->next = f;
                    235:                h = h->next; f = f->next;;
                    236:                g = g->next;
                    237:                if (f == POLYNULL) {
                    238:                  h->next = g;
                    239:                  return(node->next);
                    240:                }
                    241:                if (g == POLYNULL) {
                    242:                  h->next = f;
                    243:                  return(node->next);
                    244:                }
1.1       maekawa   245:       }else{
1.3       takayama  246:                f = f->next;
                    247:                g = g->next;
                    248:                if (f == POLYNULL) {
                    249:                  h->next = g;
                    250:                  return(node->next);
                    251:                }
                    252:                if (g == POLYNULL) {
                    253:                  h->next = f;
                    254:                  return(node->next);
                    255:                }
1.1       maekawa   256:       }
                    257:       break;
                    258:     default:
                    259:       errorPoly("ppAddv(). Internal error. Invalid value of mmLarger.");
                    260:       break;
                    261:     }
                    262:   }
                    263:   return(node->next);
                    264: }
                    265:
                    266: POLY pPower(f,k)
1.3       takayama  267:         POLY f;
                    268:         int k;
1.1       maekawa   269: {
                    270:   POLY r;
                    271:   int i,n;
                    272:   if (f == POLYNULL) return(POLYNULL); /* Is it ok? 0^0 = 0.*/
                    273:   if (k == 0) return(cxx(1,0,0,f->m->ringp));
                    274:   if (f->next == POLYNULL &&  k<0) {
                    275:     /* when f is monomial. */
                    276:     r = newCell(coeffCopy(f->coeffp),monomialCopy(f->m));
                    277:     n = r->m->ringp->n;
                    278:     for (i=0; i<n; i++) {
                    279:       r->m->e[i].x *= k;
                    280:       r->m->e[i].D *= k;
                    281:     }
                    282:     if (!isOne(r->coeffp)) {
                    283:       warningPoly("pPower(poly,negative integer) not implemented yet. Returns 1.");
                    284:       r = cxx(1,0,0,f->m->ringp);
                    285:     }
                    286:     return(r);
                    287:   }
                    288:   r = cxx(1,0,0,f->m->ringp);
                    289:   if (k < 0) {
                    290:     warningPoly("pPower(poly,negative integer) not implemented yet. Returns 1.");
                    291:   }
                    292:   for (i=0; i<k; i++) {
                    293:     r = ppMult(f,r);
                    294:   }
                    295:   return(r);
                    296: }
                    297:
                    298: POLY pPower_poly(f,k)
1.3       takayama  299:         POLY f;
                    300:         int k;
1.1       maekawa   301: {
                    302:   POLY r;
                    303:   int i,n;
                    304:   if (f == POLYNULL) return(POLYNULL); /* Is it ok? 0^0 = 0.*/
                    305:   if (k == 0) return(cxx(1,0,0,f->m->ringp));
                    306:   if (f->next == POLYNULL &&  k<0) {
                    307:     /* when f is monomial. */
                    308:     r = newCell(coeffCopy(f->coeffp),monomialCopy(f->m));
                    309:     n = r->m->ringp->n;
                    310:     for (i=0; i<n; i++) {
                    311:       r->m->e[i].x *= k;
                    312:       r->m->e[i].D *= k;
                    313:     }
                    314:     if (!isOne(r->coeffp)) {
                    315:       warningPoly("pPower_poly(poly,negative integer) not implemented yet. Returns 1.");
                    316:       r = cxx(1,0,0,f->m->ringp);
                    317:     }
                    318:     return(r);
                    319:   }
                    320:   r = cxx(1,0,0,f->m->ringp);
                    321:   if (k < 0) {
                    322:     warningPoly("pPower_poly(poly,negative integer) not implemented yet. Returns 1.");
                    323:   }
                    324:   for (i=0; i<k; i++) {
                    325:     r = ppMult_poly(f,r);
                    326:   }
                    327:   return(r);
                    328: }
                    329:
                    330: POLY modulop_trash(f,ringp)
1.3       takayama  331:         POLY f;
                    332:         struct ring *ringp;
1.1       maekawa   333: {
                    334:   int p;
                    335:   POLY h;
                    336:   MP_INT *c;
                    337:   int cc;
                    338:   POLY node;
                    339:   struct ring *nextRing;
                    340:   POLY fc;
                    341:
                    342:   if (f == POLYNULL) return(f);
                    343:   p = ringp->p;
                    344:   if (f->m->ringp->p != 0) {
                    345:     warningPoly("modulop(f,ringp) must be called with f in the characteristic 0 ring. Returns 0.");
                    346:     return(POLYNULL);
                    347:   }
                    348:   if (f->m->ringp->n != ringp->n) {
                    349:     warningPoly("modulop(f,ringp): f->m->ringp->n must be equal to ringp->n. Returns 0.");
                    350:     return(POLYNULL);
                    351:   }
                    352:
                    353:   /* The case of ringp->next != NULL */
                    354:   if (ringp->next != (struct ring *)NULL) {
                    355:     nextRing = ringp->next;
                    356:     node = newCell(newCoeff(),newMonomial(ringp));
                    357:     node->next = POLYNULL;
                    358:     h = node;
                    359:
                    360:     while (f != POLYNULL) {
                    361:       fc = bxx(f->coeffp->val.bigp,0,0,nextRing);
                    362:       h->next = newCell(newCoeff(),monomialCopy(f->m));
                    363:       h = h->next;
                    364:       h->m->ringp = ringp;
                    365:       h->coeffp->tag = POLY_COEFF;
                    366:       h->coeffp->p = p;
                    367:       h->coeffp->val.f = fc;
                    368:       f = f->next;
                    369:     }
                    370:     return(node->next);
                    371:   }
                    372:
                    373:
                    374:   /* In case of ringp->next == NULL */
                    375:   if (p != 0) {
                    376:     c = newMP_INT();
                    377:     node = newCell(newCoeff(),newMonomial(ringp));
                    378:     node->next = POLYNULL;
                    379:     h = node;
                    380:
                    381:     while (f != POLYNULL) {
                    382:       mpz_mod_ui(c,f->coeffp->val.bigp,(unsigned long int)p);
                    383:       cc = (int) mpz_get_si(c);
                    384:       if (cc != 0) {
1.3       takayama  385:                h->next = newCell(newCoeff(),monomialCopy(f->m));
                    386:                h = h->next;
                    387:                h->m->ringp = ringp;
                    388:                h->coeffp->tag = INTEGER;
                    389:                h->coeffp->p = p;
                    390:                h->coeffp->val.i = cc;
1.1       maekawa   391:       }
                    392:       f = f->next;
                    393:     }
                    394:     return(node->next);
                    395:   }else{
                    396:     h = f = pcmCopy(f);
                    397:     while (f != POLYNULL) {
                    398:       f->m->ringp = ringp;
                    399:       f = f->next;
                    400:     }
                    401:     return(h);
                    402:   }
                    403:
                    404: }
                    405:
                    406: POLY modulop(f,ringp)
1.3       takayama  407:         POLY f;
                    408:         struct ring *ringp;
                    409:         /* Z[x] ---> R[x] where R=Z, Z/Zp, ringp->next. */
1.1       maekawa   410: {
                    411:   int p;
                    412:   POLY h;
                    413:   MP_INT *c;
                    414:   int cc;
                    415:   POLY node;
                    416:   POLY fc;
                    417:
                    418:   if (f == POLYNULL) return(f);
                    419:   p = ringp->p;
                    420:   if (f->m->ringp->p != 0 || f->m->ringp->next != (struct ring *)NULL) {
                    421:     warningPoly("modulop(f,ringp) must be called with f in the characteristic 0 ring Z[x]. Returns 0.");
                    422:     return(POLYNULL);
                    423:   }
                    424:   if (f->m->ringp->n != ringp->n) {
                    425:     warningPoly("modulop(f,ringp): f->m->ringp->n must be equal to ringp->n. Returns 0.");
                    426:     return(POLYNULL);
                    427:   }
                    428:
                    429:   /* [1] The case of R = ringp->next */
                    430:   if (ringp->next != (struct ring *)NULL) {
                    431:     h = ZERO;
                    432:     while (f != POLYNULL) {
                    433:       h = ppAdd(h,mapZmonom(f,ringp));
                    434:       f = f->next;
                    435:     }
                    436:     return(h);
                    437:   }
                    438:
                    439:   /* [2] The case of R = Z/Zp */
                    440:   if (p != 0) {
                    441:     c = newMP_INT();
                    442:     node = newCell(newCoeff(),newMonomial(ringp));
                    443:     node->next = POLYNULL;
                    444:     h = node;
                    445:
                    446:     while (f != POLYNULL) {
                    447:       mpz_mod_ui(c,f->coeffp->val.bigp,(unsigned long int)p);
                    448:       cc = (int) mpz_get_si(c);
                    449:       if (cc != 0) {
1.3       takayama  450:                h->next = newCell(newCoeff(),monomialCopy(f->m));
                    451:                h = h->next;
                    452:                h->m->ringp = ringp;
                    453:                h->coeffp->tag = INTEGER;
                    454:                h->coeffp->p = p;
                    455:                h->coeffp->val.i = cc;
1.1       maekawa   456:       }
                    457:       f = f->next;
                    458:     }
                    459:     return(node->next);
                    460:   }
                    461:
                    462:   /* [3] The case of R = Z */
                    463:   h = f = pcmCopy(f);
                    464:   while (f != POLYNULL) {
                    465:     f->m->ringp = ringp;
                    466:     f = f->next;
                    467:   }
                    468:   return(h);
                    469:
                    470:
                    471: }
                    472:
                    473: POLY modulopZ(f,pcoeff)
1.3       takayama  474:         POLY f;
                    475:         struct coeff *pcoeff;
                    476:         /* Z[x] ---> Z[x] , f ---> f mod pcoeff*/
1.1       maekawa   477: {
                    478:   int p;
                    479:   POLY h;
                    480:   struct coeff *c;
                    481:   int cc;
                    482:   POLY node;
                    483:   POLY fc;
                    484:   MP_INT *bigp;
                    485:   MP_INT *tmp;
                    486:   struct ring *ringp;
                    487:
                    488:   if (f == POLYNULL) return(f);
                    489:   ringp = f->m->ringp;
                    490:   if (pcoeff->tag != MP_INTEGER) {
                    491:     warningPoly("modulopZ(): pcoeff must be a universalNumber.");
                    492:     warningPoly("Returns 0.");
                    493:     return(POLYNULL);
                    494:   }
                    495:   bigp = pcoeff->val.bigp;
                    496:   if (f->m->ringp->p != 0 || f->m->ringp->next != (struct ring *)NULL) {
                    497:     warningPoly("modulopZ(f,p) must be called with f in the characteristic 0 ring Z[x]. Returns 0.");
                    498:     return(POLYNULL);
                    499:   }
                    500:   if (f->m->ringp->n != ringp->n) {
                    501:     warningPoly("modulop(f,ringp): f->m->ringp->n must be equal to ringp->n. Returns 0.");
                    502:     return(POLYNULL);
                    503:   }
                    504:
                    505:   /* [1] The case of R = ringp->next */
                    506:   if (ringp->next != (struct ring *)NULL) {
                    507:     warningPoly("modulopZ workds only for flat polynomials. Returns 0.");
                    508:     return(POLYNULL);
                    509:   }
                    510:
                    511:   /* [2] The case of R = Z */
                    512:   node = newCell(newCoeff(),newMonomial(ringp));
                    513:   node->next = POLYNULL;
                    514:   h = node;
                    515:
                    516:   c = newCoeff();
                    517:   tmp = newMP_INT();
                    518:   while (f != POLYNULL) {
                    519:     mpz_mod(tmp,f->coeffp->val.bigp,bigp);
                    520:     if (mpz_sgn(tmp) != 0) {
                    521:       c->tag = MP_INTEGER;
                    522:       c->p = 0;
                    523:       c->val.bigp = tmp;
                    524:       h->next = newCell(c,monomialCopy(f->m));
                    525:       h = h->next;
                    526:       h->m->ringp = ringp;
                    527:
                    528:       c = newCoeff();
                    529:       tmp = newMP_INT();
                    530:     }
                    531:     f = f->next;
                    532:   }
                    533:   return(node->next);
                    534:
                    535: }
                    536:
                    537: struct pairOfPOLY quotientByNumber(f,pcoeff)
1.3       takayama  538:         POLY f;
                    539:         struct coeff *pcoeff;
                    540:         /* Z[x] ---> Z[x],Z[x] ,  f = first*pcoeff + second */
1.1       maekawa   541: {
                    542:   int p;
                    543:   POLY h;
                    544:   POLY h2;
                    545:   struct coeff *c;
                    546:   int cc;
                    547:   POLY node;
                    548:   struct coeff *c2;
                    549:   POLY node2;
                    550:   POLY fc;
                    551:   MP_INT *bigp;
                    552:   MP_INT *tmp;
                    553:   MP_INT *tmp2;
                    554:   struct ring *ringp;
                    555:   struct pairOfPOLY r;
                    556:
                    557:   if (f == POLYNULL) {
                    558:     r.first = f; r.second = f;
                    559:     return(r);
                    560:   }
                    561:   ringp = f->m->ringp;
                    562:   if (pcoeff->tag != MP_INTEGER) {
                    563:     warningPoly("quotientByNumber(): pcoeff must be a universalNumber.");
                    564:     warningPoly("Returns (0,0).");
                    565:     r.first = f; r.second = f;
                    566:     return(r);
                    567:   }
                    568:   bigp = pcoeff->val.bigp;
                    569:   if (f->m->ringp->p != 0 || f->m->ringp->next != (struct ring *)NULL) {
                    570:     warningPoly("quotientByNumber(f,p) must be called with f in the characteristic 0 ring Z[x]. Returns (0,0).");
                    571:     r.first = f; r.second = f;
                    572:     return(r);
                    573:   }
                    574:   if (f->m->ringp->n != ringp->n) {
                    575:     warningPoly("quotientByNumber(f,p): f->m->ringp->n must be equal to ringp->n. Returns 0.");
                    576:     r.first = f; r.second = f;
                    577:     return(r);
                    578:   }
                    579:
                    580:   /* [1] The case of R = ringp->next */
                    581:   if (ringp->next != (struct ring *)NULL) {
                    582:     warningPoly("quotientByNumber() workds only for flat polynomials. Returns 0.");
                    583:     r.first = f; r.second = f;
                    584:     return(r);
                    585:   }
                    586:
                    587:   /* [2] The case of R = Z */
                    588:   node = newCell(newCoeff(),newMonomial(ringp));
                    589:   node->next = POLYNULL;
                    590:   h = node;
                    591:   node2 = newCell(newCoeff(),newMonomial(ringp));
                    592:   node2->next = POLYNULL;
                    593:   h2 = node2;
                    594:
                    595:   c = newCoeff();
                    596:   tmp = newMP_INT();
                    597:   c2 = newCoeff();
                    598:   tmp2 = newMP_INT();
                    599:   while (f != POLYNULL) {
                    600:     mpz_mod(tmp,f->coeffp->val.bigp,bigp);
                    601:     if (mpz_sgn(tmp) != 0) {
                    602:       c->tag = MP_INTEGER;
                    603:       c->p = 0;
                    604:       c->val.bigp = tmp;
                    605:       h->next = newCell(c,monomialCopy(f->m));
                    606:       h = h->next;
                    607:       h->m->ringp = ringp;
                    608:
                    609:       c = newCoeff();
                    610:       tmp = newMP_INT();
                    611:     }
                    612:     mpz_tdiv_q(tmp2,f->coeffp->val.bigp,bigp);
                    613:     if (mpz_sgn(tmp2) != 0) {
                    614:       c2->tag = MP_INTEGER;
                    615:       c2->p = 0;
                    616:       c2->val.bigp = tmp2;
                    617:       h2->next = newCell(c2,monomialCopy(f->m));
                    618:       h2 = h2->next;
                    619:       h2->m->ringp = ringp;
                    620:
                    621:       c2 = newCoeff();
                    622:       tmp2 = newMP_INT();
                    623:     }
                    624:     f = f->next;
                    625:   }
                    626:   r.first = node2->next;
                    627:   r.second = node->next;
                    628:   return(r);
                    629:
                    630: }
                    631:
                    632:
                    633: POLY modulo0(f,ringp)
1.3       takayama  634:         POLY f;
                    635:         struct ring *ringp;
1.1       maekawa   636: {
                    637:   int p;
                    638:   POLY h;
                    639:   struct coeff *c;
                    640:   POLY node;
                    641:   if (f == POLYNULL) return(f);
                    642:   p = ringp->p;
                    643:   if (p != 0) {
                    644:     warningPoly("modulo0(f,ringp) must be called with the characteristic 0 ring*ringp. Returns 0.");
                    645:     return(POLYNULL);
                    646:   }
                    647:   switch (f->coeffp->tag) {
                    648:   case MP_INTEGER:
                    649:     if (f->m->ringp->p == 0) {
                    650:       node = pcmCopy(f);
                    651:       f = node;
                    652:       while (f != POLYNULL) {
1.3       takayama  653:                f->m->ringp = ringp; /* Touch the monomial "ringp" field. */
                    654:                f = f->next;
1.1       maekawa   655:       }
                    656:       return(node);
                    657:     }
                    658:     break;
                    659:   case POLY_COEFF:
                    660:     node = pcmCopy(f);
                    661:     f = node;
                    662:     while (f != POLYNULL) {
                    663:       f->m->ringp = ringp; /* Touch the monomial "ringp" field. */
                    664:       f = f->next;
                    665:     }
                    666:     return(node);
                    667:     break;
                    668:   case INTEGER:
                    669:     node = newCell(newCoeff(),newMonomial(ringp));
                    670:     node->next = POLYNULL;
                    671:     h = node;
                    672:
                    673:     while (f != POLYNULL) {
                    674:       c = newCoeff();
                    675:       c->tag = MP_INTEGER;
                    676:       c->p = 0;
                    677:       c->val.bigp = newMP_INT();
                    678:       mpz_set_si(c->val.bigp,f->coeffp->val.i);
                    679:       h->next = newCell(c,monomialCopy(f->m));
                    680:       h = h->next;
                    681:       h->m->ringp = ringp;
                    682:       f = f->next;
                    683:     }
                    684:     return(node->next);
                    685:     break;
                    686:   default:
                    687:     warningPoly("modulo0(): coefficients have to be MP_INTEGER or INTEGER. Returns 0");
                    688:     return(POLYNULL);
                    689:     break;
                    690:   }
                    691: }
                    692:
                    693:
                    694: struct object test(ob)  /* test3 */
1.3       takayama  695:         struct object ob;
1.1       maekawa   696: {
1.6     ! takayama  697:   struct object rob = OINIT;
1.1       maekawa   698:   int k;
                    699:   static POLY f0;
                    700:   static POLY f1;
                    701:
                    702:   POLY addNode;
                    703:   POLY f;
                    704:   int i;
                    705:   static int s=0;
                    706:   MP_INT *mp;
                    707:   extern struct ring *SmallRingp;
                    708:   extern struct ring *CurrentRingp;
                    709:   addNode = pMalloc(SmallRingp);
                    710:   k = ob.lc.ival;
                    711:   switch(s) {
                    712:   case 0:
                    713:     f0 = addNode;
                    714:     for (i=k; i>=0; i--) {
                    715:       f0->next = bxx(BiiPower(-k,i),0,i,CurrentRingp);
                    716:       if (f0->next != POLYNULL) {
1.3       takayama  717:                f0 = f0->next;
1.1       maekawa   718:       }
                    719:     }
                    720:     f0 = addNode->next;
                    721:     s++;
                    722:     rob.lc.poly = f0;
                    723:     break;
                    724:   case 1:
                    725:     f1 = addNode;
                    726:     for (i=k; i>=0; i--) {
                    727:       f1->next = bxx(BiiPower(k,i),0,i,CurrentRingp);
                    728:       if (f1->next != POLYNULL) {
1.3       takayama  729:                f1 = f1->next;
1.1       maekawa   730:       }
                    731:     }
                    732:     f1 = addNode->next;
                    733:     s = 0;
                    734:     rob.lc.poly = f1;
                    735:     break;
                    736:   default:
                    737:     rob.lc.poly = POLYNULL;
                    738:     s = 0;
                    739:     break;
                    740:   }
                    741:
                    742:
                    743:   rob.tag = Spoly;
                    744:   return(rob);
                    745: }
                    746:
                    747:
                    748: int pLength(f)
1.3       takayama  749:         POLY f;
1.1       maekawa   750: {
                    751:   int c=0;
                    752:   if (f ISZERO) return(0);
                    753:   while (f != POLYNULL) {
                    754:     c++;
                    755:     f = f->next;
                    756:   }
                    757:   return(c);
                    758: }
                    759:
                    760:
                    761: POLY ppAddv2(f,g,top,nexttop)
1.3       takayama  762:         POLY f; POLY g;  /* It breaks f and g. Use it just after calling mpMult() */
                    763:         POLY top;
                    764:         POLY *nexttop;
                    765:         /* top is the starting address in the list f.
                    766:                if top == POLYNULL, start from f.
                    767:
                    768:                *nexttop == 0
                    769:                == g
                    770:                == h or 0
1.1       maekawa   771:
1.3       takayama  772:                It must be called as r = ppAddv2(r,g,...);
1.1       maekawa   773: */
                    774: {
                    775:   POLY node;
                    776:   struct listPoly nod;
                    777:   POLY h;
                    778:   struct coeff *c;
                    779:   int gt;
                    780:   POLY g0;
                    781:   POLY f0; /* for debug */
                    782:
                    783:   node = &nod;
                    784:   /* printf("ppAddv:1%s + %s\n",POLYToString(f,'*',1),POLYToString(g,'*',1)); */
                    785:   node->next = POLYNULL;
                    786:   h = node;
                    787:   *nexttop = POLYNULL;
                    788:   if (f == POLYNULL) return(g);
                    789:   if (g == POLYNULL) return(f);
                    790:   checkRing(f,g);
                    791:
                    792:   f0 = f;
                    793:   if (top != POLYNULL) {
                    794:     while (f != top) {
                    795:       if (f == POLYNULL) {
1.3       takayama  796:                fprintf(stderr,"\nppAddv2(): Internal error.\n");fflush(stderr);
                    797:                fprintf(stderr,"f = %s\n",POLYToString(f0,'*',0));
                    798:                fprintf(stderr,"g = %s\n",POLYToString(g0,'*',0));
                    799:                fprintf(stderr,"top=%s\n",POLYToString(top,'*',0));
                    800:                errorPoly("ppAddv2(). Internal error=1.");
1.1       maekawa   801:       }
                    802:       h->next = f;
                    803:       h = h->next;
                    804:       f = f->next;
                    805:     }
                    806:   }
                    807:   g0 = g;
                    808:   *nexttop = g0;
                    809:
                    810:   while (f != POLYNULL && g != POLYNULL) {
                    811:     checkRing2(f,g); /* for debug */
                    812:     /* printf("%s + %s\n",POLYToString(f,'*',1),POLYToString(g,'*',1)); */
                    813:     gt = (*mmLarger)(f,g);
                    814:     switch (gt) {
                    815:     case 1: /* f > g */
                    816:       h->next = f;
                    817:       h = h->next; f = f->next;;
                    818:       if (f == POLYNULL) {
1.3       takayama  819:                h->next = g;
                    820:                return(node->next);
1.1       maekawa   821:       }
                    822:       break;
                    823:     case 0: /* f < g */
                    824:       h->next =  g;
                    825:       h = h->next; g = g->next;
                    826:       if (g == POLYNULL) {
1.3       takayama  827:                h->next = f;
                    828:                return(node->next);
1.1       maekawa   829:       }
                    830:       break;
                    831:     case 2:/* f == g */
                    832:       c = g->coeffp;
                    833:       Cadd(c,f->coeffp,c);
                    834:       if (!isZero(c)) {
1.3       takayama  835:                h->next = g;
                    836:                h = h->next;
                    837:                f = f->next;;
                    838:                g = g->next;
                    839:                if (f == POLYNULL) {
                    840:                  h->next = g;
                    841:                  return(node->next);
                    842:                }
                    843:                if (g == POLYNULL) {
                    844:                  h->next = f;
                    845:                  return(node->next);
                    846:                }
1.1       maekawa   847:       }else{
1.3       takayama  848:                if (g == g0) {
                    849:                  if (h != node) {
                    850:                        *nexttop = h;
                    851:                  }else{
                    852:                        *nexttop = POLYNULL;
                    853:                  }
                    854:                }
1.1       maekawa   855:
1.3       takayama  856:                f = f->next;
                    857:                g = g->next;
1.1       maekawa   858:
1.3       takayama  859:                if (f == POLYNULL) {
                    860:                  h->next = g;
                    861:                  return(node->next);
                    862:                }
                    863:                if (g == POLYNULL) {
                    864:                  h->next = f;
                    865:                  return(node->next);
                    866:                }
1.1       maekawa   867:       }
                    868:       break;
                    869:     default:
                    870:       errorPoly("ppAddv(). Internal error. Invalid value of mmLarger.");
                    871:       break;
                    872:     }
                    873:   }
                    874:   return(node->next);
                    875: }
                    876:
                    877: POLY ppMult(f,g)
1.3       takayama  878:         POLY f,g;
1.1       maekawa   879: {
                    880:   POLY r;
                    881:   POLY tmp;
                    882:   POLY top;
                    883:   POLY nexttop;
                    884:   r = POLYNULL; top = POLYNULL;
                    885:   while( f != POLYNULL) {
                    886:     /* tmp = (*mpMult)(f,g);  (*mpMult) is no more used. */
                    887:     tmp = (*(f->m->ringp->multiplication))(f,g);
                    888:     /*printf("mpMult(%s,%s) ->%s\n",POLYToString(f,'*',1),POLYToString(g,'*',1),POLYToString(tmp,'*',1)); */
                    889:     r = ppAddv2(r,tmp,top,&nexttop); /* r and tmp will be broken */
                    890:     top = nexttop;
                    891:     f = f->next;
                    892:   }
                    893:   return(r);
                    894: }
                    895:
                    896: POLY ppMult_poly(f,g)
1.3       takayama  897:         POLY f,g;
1.1       maekawa   898: {
                    899:   POLY r;
                    900:   POLY tmp;
                    901:   POLY top;
                    902:   POLY nexttop;
                    903:   r = POLYNULL; top = POLYNULL;
                    904:   while( f != POLYNULL) {
                    905:     tmp = mpMult_poly(f,g);
                    906:     r = ppAddv2(r,tmp,top,&nexttop); /* r and tmp will be broken */
                    907:     top = nexttop;
                    908:     f = f->next;
                    909:   }
                    910:   return(r);
                    911: }
                    912:
                    913: POLY mapZmonom(f,ringp)
1.3       takayama  914:         POLY f; /* assumes monomial. f \in Z[x] */
                    915:         struct ring *ringp;  /* R[x] */
1.1       maekawa   916: {
                    917:   struct ring *nextRing;
                    918:   struct ring nextRing0;
                    919:   POLY ff;
                    920:   POLY node;
                    921:   POLY gg;
                    922:   int l,c,d;
                    923:
                    924:   nextRing = ringp->next;
                    925:   nextRing0 = *nextRing; nextRing0.p = 0; nextRing0.next = 0;
                    926:   /* nextRing0 = Z[y] where y is the variables of R. */
                    927:
                    928:   ff = bxx(f->coeffp->val.bigp,0,0,&nextRing0);
                    929:   ff = modulop(ff,nextRing);
                    930:
                    931:   node = newCell(newCoeff(),monomialCopy(f->m));
                    932:   node->next = POLYNULL;
                    933:   node->m->ringp = ringp;
                    934:
                    935:   node->coeffp->p = ringp->p;
                    936:
                    937:   l = ringp->l; c = ringp->c;
                    938:   /* If q-analog  q x ---> (q) x. */
                    939:   if (l-c > 0) {
                    940:     d = node->m->e[0].x;  /* degree of q in R[x].*/
                    941:     node->m->e[0].x = 0;
                    942:     gg = cxx(1,0,d,nextRing); /* q^d = x[0]^d */
                    943:     ff = ppMult(gg,ff);
                    944:   }
                    945:
                    946:   node->coeffp->tag = POLY_COEFF;
                    947:   node->coeffp->val.f = ff;
                    948:   return(node);
                    949: }
1.4       takayama  950:
                    951: POLY reduceContentOfPoly(POLY f,struct coeff **contp) {
                    952:   struct coeff *cont;
                    953:   struct coeff *cOne = NULL;
                    954:   extern struct ring *SmallRingp;
                    955:   if (cOne == NULL) cOne = intToCoeff(1,SmallRingp);
1.5       takayama  956:   *contp = cOne;
1.4       takayama  957:
                    958:   if (f == POLYNULL) return f;
                    959:   if (f->m->ringp->p != 0) return f;
                    960:   if (f->coeffp->tag != MP_INTEGER) return f;
                    961:   cont = gcdOfCoeff(f);
                    962:   *contp = cont;
                    963:   if (coeffGreater(cont,cOne)) {
                    964:        f = quotientByNumber(f,cont).first;
                    965:   }
                    966:   return f;
                    967: }
                    968:
                    969: int coeffSizeMin(POLY f) {
                    970:   int size;
                    971:   int t;
                    972:   if (f == POLYNULL) return 0;
                    973:   if (f->m->ringp->p != 0) return 0;
                    974:   if (f->coeffp->tag != MP_INTEGER) return 0;
                    975:   size = mpz_size(f->coeffp->val.bigp);
                    976:   while (f != POLYNULL) {
                    977:        t = mpz_size(f->coeffp->val.bigp);
                    978:        if (t < size)  size = t;
                    979:        if (size == 1) return size;
                    980:        f = f->next;
                    981:   }
                    982: }
                    983:
                    984: struct coeff *gcdOfCoeff(POLY f) {
                    985:   extern struct ring *SmallRingp;
                    986:   struct coeff *t;
                    987:   MP_INT *tmp;
                    988:   MP_INT *tmp2;
                    989:   static MP_INT *cOne = NULL;
                    990:   if (cOne == NULL) {
                    991:        cOne = newMP_INT();
                    992:        mpz_set_si(cOne,(long) 1);
                    993:   }
1.5       takayama  994:   if (f == POLYNULL) return intToCoeff(0,SmallRingp);
1.4       takayama  995:   if (f->m->ringp->p != 0) return intToCoeff(0,SmallRingp);
                    996:   if (f->coeffp->tag != MP_INTEGER) return intToCoeff(0,SmallRingp);
                    997:   tmp = f->coeffp->val.bigp;
                    998:   tmp2 = newMP_INT();
                    999:   while (f != POLYNULL) {
                   1000:     mpz_gcd(tmp2,tmp,f->coeffp->val.bigp); /* tmp = tmp2 OK? */
                   1001:        tmp = tmp2;
                   1002:        if (mpz_cmp(tmp,cOne)==0) return intToCoeff(1,SmallRingp);
                   1003:        f = f->next;
                   1004:   }
                   1005:   return mpintToCoeff(tmp,SmallRingp);
1.1       maekawa  1006:
1.5       takayama 1007: }
                   1008:
                   1009: int shouldReduceContent(POLY f,int ss) {
                   1010:   extern DoCancel;
                   1011:   static int prevSize = 1;
                   1012:   int size;
                   1013:   if (f == POLYNULL) return 0;
                   1014:   if (f->m->ringp->p != 0) return 0;
                   1015:   if (f->coeffp->tag != MP_INTEGER) return 0;
                   1016:   if (DoCancel & 2) return 1;
                   1017:   /* Apply the Noro strategy to reduce content. */
                   1018:   size = mpz_size(f->coeffp->val.bigp);
                   1019:   if (ss > 0) {
                   1020:        prevSize = size;
                   1021:        return 0;
                   1022:   }
                   1023:   if (size > 2*prevSize) {
                   1024:        return 1;
                   1025:   }else{
                   1026:        return 0;
                   1027:   }
1.4       takayama 1028: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>