[BACK]Return to base5.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari / src / basemath

Annotation of OpenXM_contrib/pari/src/basemath/base5.c, Revision 1.1.1.1

1.1       maekawa     1: /*******************************************************************/
                      2: /*                                                                 */
                      3: /*                       BASIC NF OPERATIONS                       */
                      4: /*                          (continued 2)                          */
                      5: /*                                                                 */
                      6: /*******************************************************************/
                      7: /* $Id: base5.c,v 1.1.1.1 1999/09/16 13:47:22 karim Exp $ */
                      8: #include "pari.h"
                      9: GEN mat_to_vecpol(GEN x, long v);
                     10:
                     11: GEN
                     12: matbasistoalg(GEN nf,GEN x)
                     13: {
                     14:   long i,j,lx,li;
                     15:   GEN p1,z;
                     16:
                     17:   if (typ(x)!=t_MAT)
                     18:     err(talker,"argument must be a matrix in matbasistoalg");
                     19:   lx=lg(x); z=cgetg(lx,t_MAT); if (lx==1) return z;
                     20:
                     21:   li=lg(x[1]);
                     22:   for (j=1; j<lx; j++)
                     23:   {
                     24:     p1=cgetg(li,t_COL); z[j]=(long)p1;
                     25:     for (i=1; i<li; i++) p1[i]=(long)basistoalg(nf,gcoeff(x,i,j));
                     26:   }
                     27:   return z;
                     28: }
                     29:
                     30: GEN
                     31: matalgtobasis(GEN nf,GEN x)
                     32: {
                     33:   long i,j,lx,li;
                     34:   GEN p1,z;
                     35:
                     36:   if (typ(x)!=t_MAT)
                     37:     err(talker,"argument must be a matrix in matalgtobasis");
                     38:   lx=lg(x); z=cgetg(lx,t_MAT); if (lx==1) return z;
                     39:
                     40:   li=lg(x[1]);
                     41:   for (j=1; j<lx; j++)
                     42:   {
                     43:     p1=cgetg(li,t_COL); z[j]=(long)p1;
                     44:     for (i=1; i<li; i++) p1[i]=(long)algtobasis(nf,gcoeff(x,i,j));
                     45:   }
                     46:   return z;
                     47: }
                     48:
                     49: static GEN
                     50: rnfmakematrices(GEN rnf)
                     51: {
                     52:   long i,j,k,n,r1,r2,ru,ruk,r1rel,r2rel;
                     53:   GEN nf,pol,rac,base,base1,r1r2,racnf,sig,vecmat,vecM,vecMC,vecT2,rack;
                     54:   GEN M,p2,p3,MC,sigk,T2,T,p1,MD,TI,MDI;
                     55:
                     56:   nf=(GEN)rnf[10]; racnf=(GEN)nf[6]; pol=(GEN)rnf[1];
                     57:   n=lgef(pol)-3;
                     58:   base=(GEN)rnf[7]; base1=(GEN)base[1]; rac=(GEN)rnf[6]; sig=(GEN)rnf[2];
                     59:   r1r2=(GEN)nf[2]; r1=itos((GEN)r1r2[1]); r2=itos((GEN)r1r2[2]); ru=r1+r2;
                     60:   vecmat=cgetg(8,t_VEC);
                     61:   vecM=cgetg(ru+1,t_VEC); vecmat[1]=(long)vecM;
                     62:   vecMC=cgetg(ru+1,t_VEC); vecmat[2]=(long)vecMC;
                     63:   vecT2=cgetg(ru+1,t_VEC); vecmat[3]=(long)vecT2;
                     64:   for (k=1; k<=ru; k++)
                     65:   {
                     66:     rack=(GEN)rac[k]; ruk=lg(rack)-1;
                     67:     M=cgetg(n+1,t_MAT); vecM[k]=(long)M;
                     68:     for (j=1; j<=n; j++)
                     69:     {
                     70:       p2=cgetg(ruk+1,t_COL); M[j]=(long)p2; p3=lift((GEN)base1[j]);
                     71:       p3=gsubst(p3,varn(nf[1]),(GEN)racnf[k]);
                     72:       for (i=1; i<=ruk; i++) p2[i]=lsubst(p3,varn(rnf[1]),(GEN)rack[i]);
                     73:     }
                     74:     MC=gconj(gtrans(M)); vecMC[k]=(long)MC;
                     75:     if (k<=r1)
                     76:     {
                     77:       sigk=(GEN)sig[k]; r1rel=itos((GEN)sigk[1]); r2rel=itos((GEN)sigk[2]);
                     78:       if (r1rel+r2rel != lg(MC)-1) err(talker,"bug in rnfmakematrices");
                     79:       for (j=r1rel+1; j<=r1rel+r2rel; j++) MC[j]=lmul2n((GEN)MC[j],1);
                     80:     }
                     81:     T2=gmul(MC,M); vecT2[k]=(long)T2;
                     82:   }
                     83:   T=cgetg(n+1,t_MAT); vecmat[4]=(long)T;
                     84:   for (j=1; j<=n; j++)
                     85:   {
                     86:     p1=cgetg(n+1,t_COL); T[j]=(long)p1;
                     87:     for (i=1; i<=n; i++)
                     88:       p1[i]=ltrace(gmodulcp(gmul((GEN)base1[i],(GEN)base1[j]),pol));
                     89:   }
                     90:   MD=cgetg(1,t_MAT); vecmat[5]=(long)MD; /* matrice de la differente */
                     91:   TI=cgetg(1,t_MAT); vecmat[6]=(long)TI; /* matrice .... ? */
                     92:   MDI=cgetg(1,t_MAT); vecmat[7]=(long)MDI; /* matrice .... ? */
                     93:   return vecmat;
                     94: }
                     95:
                     96: GEN
                     97: rnfinitalg(GEN nf,GEN pol,long prec)
                     98: {
                     99:   long av=avma,tetpil,m,n,r1,r2,vnf,i,j,k,vpol,v1,r1j,r2j,lfac,degabs;
                    100:   GEN RES,sig,r1r2,rac,p1,p2,liftpol,delta,RAC,ro,p3,bas;
                    101:   GEN f,f2,fac,fac1,fac2,id,p4,p5;
                    102:
                    103:   if (typ(pol)!=t_POL) err(notpoler,"rnfinitalg");
                    104:   nf=checknf(nf); n=lgef(pol)-3; vpol=varn(pol);
                    105:   vnf=0;
                    106:   for (i=0; i<=n; i++)
                    107:   {
                    108:     long tp1;
                    109:
                    110:     p1=(GEN)pol[i+2];
                    111:     tp1=typ(p1);
                    112:     if (! is_const_t(tp1))
                    113:     {
                    114:       if (tp1!=t_POLMOD) err(typeer,"rnfinitalg");
                    115:       if (!gegal((GEN)p1[1],(GEN)nf[1]))
                    116:        err(talker,"incompatible number fields in rnfinitalg");
                    117:       p1=(GEN)p1[2];
                    118:       if (! is_const_t(typ(p1)))
                    119:       {
                    120:        v1=varn(p1);
                    121:        if (vnf && vnf!=v1) err(talker,"different variables in rnfinitalg");
                    122:        if (!vnf) vnf=v1;
                    123:       }
                    124:     }
                    125:   }
                    126:   if (!vnf) vnf=varn(nf[1]);
                    127:   if (vpol>=vnf)
                    128:     err(talker,"main variable must be of higher priority in rnfinitalg");
                    129:   RES=cgetg(12,t_VEC);
                    130:   RES[1]=(long)pol;
                    131:   m=lgef(nf[1])-3; degabs=n*m;
                    132:   r1r2=(GEN)nf[2]; r1=itos((GEN)r1r2[1]); r2=itos((GEN)r1r2[2]);
                    133:   sig=cgetg(r1+r2+1,t_VEC); RES[2]=(long)sig;
                    134:   rac=(GEN)nf[6]; liftpol=lift(pol);
                    135:   RAC=cgetg(r1+r2+1,t_VEC); RES[6]=(long)RAC;
                    136:   for (j=1; j<=r1; j++)
                    137:   {
                    138:     p1=gsubst(liftpol,vnf,(GEN)rac[j]);
                    139:     ro=roots(p1,prec);
                    140:     r1j=0;
                    141:     while (r1j<n && gcmp0(gimag((GEN)ro[r1j+1]))) r1j++;
                    142:     p2=cgetg(3,t_VEC); p2[1]=lstoi(r1j); p2[2]=lstoi(r2j=((n-r1j)>>1));
                    143:     sig[j]=(long)p2;
                    144:     p3=cgetg(r1j+r2j+1,t_VEC);
                    145:     for (i=1; i<=r1j; i++) p3[i]=lreal((GEN)ro[i]);
                    146:     for (; i<=r1j+r2j; i++) p3[i]=(long)ro[(i<<1)-r1j];
                    147:     RAC[j]=(long)p3;
                    148:   }
                    149:   for (; j<=r1+r2; j++)
                    150:   {
                    151:     p2=cgetg(3,t_VEC); p2[1]=zero; p2[2]=lstoi(n); sig[j]=(long)p2;
                    152:     p1=gsubst(liftpol,vnf,(GEN)rac[j]);
                    153:     RAC[j]=(long)roots(p1,prec);
                    154:   }
                    155:   p1 = rnfpseudobasis(nf,pol);
                    156:
                    157:   delta = cgetg(3,t_VEC);
                    158:     delta[1]=p1[3];
                    159:     delta[2]=p1[4];
                    160:   RES[3]=(long)delta;
                    161:   p2 = matbasistoalg(nf,(GEN)p1[1]);
                    162:   bas = cgetg(3,t_VEC);
                    163:     bas[1]=(long)mat_to_vecpol(p2,vpol);
                    164:     bas[2]=(long)p1[2];
                    165:   RES[7]=(long)bas;
                    166:   RES[8]=linvmat(p2);
                    167:
                    168:   f2=idealdiv(nf,discsr(pol),(GEN)p1[3]);
                    169:   fac=idealfactor(nf,f2);
                    170:   fac1=(GEN)fac[1]; fac2=(GEN)fac[2]; lfac=lg(fac1)-1;
                    171:   f=idmat(m);
                    172:   for (i=1; i<=lfac; i++)
                    173:   {
                    174:     if (mpodd((GEN)fac2[i])) err(bugparier,"rnfinitalg (odd exponent)");
                    175:     f=idealmul(nf,f,idealpow(nf,(GEN)fac1[i],gmul2n((GEN)fac2[i],-1)));
                    176:   }
                    177:   RES[4]=(long)f;
                    178:   RES[10]=(long)nf;
                    179:   RES[5]=(long)rnfmakematrices(RES);
                    180:   if (DEBUGLEVEL>1) msgtimer("matrices");
                    181:   RES[9]=lgetg(1,t_VEC); /* table de multiplication */
                    182:   p2=cgetg(6,t_VEC); RES[11]=(long)p2;
                    183:   p1=rnfequation2(nf,pol); for (i=1; i<=3; i++) p2[i]=p1[i];
                    184:   p4=cgetg(degabs+1,t_MAT);
                    185:   for (i=1; i<=n; i++)
                    186:   { /* removing denominators speeds up multiplication */
                    187:     GEN cop3,com, om = rnfelementreltoabs(RES,gmael(bas,1,i));
                    188:
                    189:     if (DEBUGLEVEL>1) msgtimer("i = %ld",i);
                    190:     com = content(om); om = gdiv(om,com);
                    191:     id=gmael(bas,2,i);
                    192:     for (j=1; j<=m; j++)
                    193:     {
                    194:       p5=cgetg(degabs+1,t_COL); p4[(i-1)*m+j]=(long)p5;
                    195:       p1=gmul((GEN)nf[7],(GEN)id[j]);
                    196:       p3 = gsubst(p1,varn(nf[1]), (GEN)p2[2]);
                    197:       cop3 = content(p3); p3 = gdiv(p3,cop3);
                    198:       p3 = gmul(gmul(com,cop3), lift_intern(gmul(om,p3)));
                    199:
                    200:       for (k=1; k<lgef(p3)-1; k++) p5[k]=p3[k+1];
                    201:       for (   ; k<=degabs;    k++) p5[k]=zero;
                    202:     }
                    203:   }
                    204:   if (DEBUGLEVEL>1) msgtimer("p4");
                    205:   p3=denom(p4); if (gcmp1(p3)) p3=NULL; else p4=gmul(p3,p4);
                    206:   p4=hnfmod(p4,detint(p4));
                    207:   if (DEBUGLEVEL>1) msgtimer("hnfmod");
                    208:   for (j=degabs-1; j>0; j--)
                    209:     if (cmpis(gcoeff(p4,j,j),2) > 0)
                    210:     {
                    211:       p1=shifti(gcoeff(p4,j,j),-1);
                    212:       for (k=j+1; k<=degabs; k++)
                    213:         if (cmpii(gcoeff(p4,j,k),p1) > 0)
                    214:           for (i=1; i<=j; i++)
                    215:             coeff(p4,i,k)=lsubii(gcoeff(p4,i,k),gcoeff(p4,i,j));
                    216:     }
                    217:   if (p3) p4=gdiv(p4,p3);
                    218:   p2[4]=(long)mat_to_vecpol(p4,vpol);
                    219:   p2[5]=linvmat(p4);
                    220:   tetpil=avma; return gerepile(av,tetpil,gcopy(RES));
                    221: }
                    222:
                    223: GEN
                    224: rnfbasistoalg(GEN rnf,GEN x)
                    225: {
                    226:   long tx=typ(x),lx=lg(x),av=avma,tetpil,i,n;
                    227:   GEN p1,z,nf;
                    228:
                    229:   checkrnf(rnf); nf=(GEN)rnf[10];
                    230:   switch(tx)
                    231:   {
                    232:     case t_VEC:
                    233:       x=gtrans(x); /* fall through */
                    234:     case t_COL:
                    235:       n=lg(x)-1; p1=cgetg(n+1,t_COL);
                    236:       for (i=1; i<=n; i++)
                    237:       {
                    238:        if (typ(x[i])==t_COL) p1[i]=(long)basistoalg(nf,(GEN)x[i]);
                    239:        else p1[i]=x[i];
                    240:       }
                    241:       p1=gmul(gmael(rnf,7,1),p1); tetpil=avma;
                    242:       return gerepile(av,tetpil,gmodulcp(p1,(GEN)rnf[1]));
                    243:
                    244:     case t_MAT:
                    245:       z=cgetg(lx,tx);
                    246:       for (i=1; i<lx; i++) z[i]=(long)rnfbasistoalg(rnf,(GEN)x[i]);
                    247:       return z;
                    248:
                    249:     case t_POLMOD:
                    250:       return gcopy(x);
                    251:
                    252:     default:
                    253:       z=cgetg(3,t_POLMOD); z[1]=lcopy((GEN)rnf[1]);
                    254:       z[2]=lmul(x,polun[varn(rnf[1])]); return z;
                    255:   }
                    256: }
                    257:
                    258: long polegal_spec(GEN x, GEN y);
                    259:
                    260: /* assuem x is a t_POLMOD */
                    261: GEN
                    262: lift_to_pol(GEN x)
                    263: {
                    264:   GEN y = (GEN)x[2];
                    265:   return (typ(y) != t_POL)? gtopoly(y,varn(x[1])): y;
                    266: }
                    267:
                    268: GEN
                    269: rnfalgtobasis(GEN rnf,GEN x)
                    270: {
                    271:   long av=avma,tetpil,tx=typ(x),lx=lg(x),i,N;
                    272:   GEN z;
                    273:
                    274:   checkrnf(rnf);
                    275:   switch(tx)
                    276:   {
                    277:     case t_VEC: case t_COL: case t_MAT:
                    278:       z=cgetg(lx,tx);
                    279:       for (i=1; i<lx; i++) z[i]=(long)rnfalgtobasis(rnf,(GEN)x[i]);
                    280:       return z;
                    281:
                    282:     case t_POLMOD:
                    283:       if (!polegal_spec((GEN)rnf[1],(GEN)x[1]))
                    284:        err(talker,"not the same number field in rnfalgtobasis");
                    285:       x=lift_to_pol(x); /* fall through */
                    286:     case t_POL:
                    287:       N=lgef(rnf[1])-3;
                    288:       if (tx==t_POL && lgef(x)-3 >= N) x=gmod(x,(GEN)rnf[1]);
                    289:       z=cgetg(N+1,t_COL); for (i=1; i<=N; i++) z[i]=(long)truecoeff(x,i-1);
                    290:       tetpil=avma; return gerepile(av,tetpil,gmul((GEN)rnf[8],z));
                    291:   }
                    292:   return gscalcol(x, lgef(rnf[1])-3);
                    293: }
                    294:
                    295: /* x doit etre un polymod ou un polynome ou un vecteur de tels objets... */
                    296: GEN
                    297: rnfelementreltoabs(GEN rnf,GEN x)
                    298: {
                    299:   long av=avma,tx,i,lx,va,tp3;
                    300:   GEN z,p1,p2,p3,polabs,teta,alpha,s,k;
                    301:
                    302:   checkrnf(rnf); tx=typ(x); lx=lg(x); va=varn((GEN)rnf[1]);
                    303:   switch(tx)
                    304:   {
                    305:     case t_VEC: case t_COL: case t_MAT:
                    306:       z=cgetg(lx,tx);
                    307:       for (i=1; i<lx; i++) z[i]=(long)rnfelementreltoabs(rnf,(GEN)x[i]);
                    308:       return z;
                    309:
                    310:     case t_POLMOD:
                    311:       x=lift_to_pol(x); /* fall through */
                    312:     case t_POL:
                    313:       if (gvar(x) > va)
                    314:       {
                    315:        if (gcmp0(x)) {x=cgetg(2,t_POL); x[1]=evalvarn(va) | evallgef(2);}
                    316:        else
                    317:        {
                    318:          p1=cgetg(3,t_POL); p1[1]=evalvarn(va) | evallgef(3) | evalsigne(1);
                    319:          p1[2]=(long)x; x=p1;
                    320:        }
                    321:       }
                    322:       p1=(GEN)rnf[11]; polabs=(GEN)p1[1]; alpha=(GEN)p1[2]; k=(GEN)p1[3];
                    323:       teta=gmodulcp(gsub(polx[va],gmul(k,(GEN)alpha[2])),polabs);
                    324:       s=gzero;
                    325:       for (i=lgef(x)-1; i>1; i--)
                    326:       {
                    327:        p3=(GEN)x[i]; tp3=typ(p3);
                    328:        if (is_const_t(tp3)) p2 = p3;
                    329:        else
                    330:          switch(tp3)
                    331:          {
                    332:            case t_POLMOD:
                    333:              p3 = (GEN)p3[2]; /* fall through */
                    334:            case t_POL:
                    335:              p2 = poleval(p3,alpha);
                    336:          }
                    337:        s=gadd(p2,gmul(teta,s));
                    338:       }
                    339:       return gerepileupto(av,s);
                    340:
                    341:     default: return gcopy(x);
                    342:   }
                    343: }
                    344:
                    345: GEN
                    346: rnfelementabstorel(GEN rnf,GEN x)
                    347: {
                    348:   long av=avma,tx,i,lx;
                    349:   GEN z,p1,s,tetap,k,nf;
                    350:
                    351:   checkrnf(rnf); tx=typ(x); lx=lg(x);
                    352:   switch(tx)
                    353:   {
                    354:     case t_VEC: case t_COL: case t_MAT:
                    355:       z=cgetg(lx,tx);
                    356:       for (i=1; i<lx; i++) z[i]=(long)rnfelementabstorel(rnf,(GEN)x[i]);
                    357:       return z;
                    358:
                    359:     case t_POLMOD:
                    360:       x=lift_to_pol(x); /* fall through */
                    361:     case t_POL:
                    362:       p1=(GEN)rnf[11]; k=(GEN)p1[3]; nf=(GEN)rnf[10];
                    363:       tetap=gmodulcp(gadd(polx[varn(rnf[1])],
                    364:            gmul(k,gmodulcp(polx[varn(nf[1])],(GEN)nf[1]))),(GEN)rnf[1]);
                    365:       s=gzero;
                    366:       for (i=lgef(x)-1; i>1; i--) s=gadd((GEN)x[i],gmul(tetap,s));
                    367:       return gerepileupto(av,s);
                    368:
                    369:     default: return gcopy(x);
                    370:   }
                    371: }
                    372:
                    373: /* x doit etre un polymod ou un polynome ou un vecteur de tels objets... */
                    374: GEN
                    375: rnfelementup(GEN rnf,GEN x)
                    376: {
                    377:   long i,lx,tx;
                    378:   GEN z;
                    379:
                    380:   checkrnf(rnf); tx=typ(x); lx=lg(x);
                    381:   switch(tx)
                    382:   {
                    383:     case t_VEC: case t_COL: case t_MAT:
                    384:       z=cgetg(lx,tx);
                    385:       for (i=1; i<lx; i++) z[i]=(long)rnfelementup(rnf,(GEN)x[i]);
                    386:       return z;
                    387:
                    388:     case t_POLMOD:
                    389:       x=(GEN)x[2]; /* fall through */
                    390:     case t_POL:
                    391:       return poleval(x,gmael(rnf,11,2));
                    392:
                    393:     default: return gcopy(x);
                    394:   }
                    395: }
                    396:
                    397: /* x doit etre un polymod ou un polynome ou un vecteur de tels objets..*/
                    398: GEN
                    399: rnfelementdown(GEN rnf,GEN x)
                    400: {
                    401:   long av=avma,tetpil,i,lx,tx;
                    402:   GEN p1,z;
                    403:
                    404:   checkrnf(rnf); tx=typ(x); lx=lg(x);
                    405:   switch(tx)
                    406:   {
                    407:     case t_VEC: case t_COL: case t_MAT:
                    408:       z=cgetg(lx,tx);
                    409:       for (i=1; i<lx; i++) z[i]=(long)rnfelementdown(rnf,(GEN)x[i]);
                    410:       return z;
                    411:
                    412:     case t_POLMOD:
                    413:       x=(GEN)x[2]; /* fall through */
                    414:     case t_POL:
                    415:       if (gcmp0(x)) return gzero;
                    416:
                    417:       p1=rnfelementabstorel(rnf,x);
                    418:       if (typ(p1)==t_POLMOD && varn(p1[1])==varn(rnf[1])) p1=(GEN)p1[2];
                    419:       if (gvar(p1)>varn(rnf[1]))
                    420:       {
                    421:        tetpil=avma;
                    422:        return gerepile(av,tetpil,gcopy(p1));
                    423:       }
                    424:       if (lgef(p1)==3)
                    425:       {
                    426:        tetpil=avma;
                    427:         return gerepile(av,tetpil,gcopy((GEN)p1[2]));
                    428:       }
                    429:       err(talker,"element is not in the base field in rnfelementdown");
                    430:
                    431:     default: return gcopy(x);
                    432:   }
                    433: }
                    434:
                    435: /* x est exprime sur la base relative */
                    436: static GEN
                    437: rnfprincipaltohermite(GEN rnf,GEN x)
                    438: {
                    439:   long av=avma,tetpil;
                    440:   GEN nf,bas,bas1,p1,z;
                    441:
                    442:   x=rnfbasistoalg(rnf,x); nf=(GEN)rnf[10];
                    443:   bas=(GEN)rnf[7]; bas1=(GEN)bas[1];
                    444:   p1=rnfalgtobasis(rnf,gmul(x,gmodulcp(bas1,(GEN)rnf[1])));
                    445:   z=cgetg(3,t_VEC); z[2]=bas[2];
                    446:   settyp(p1,t_MAT); z[1]=(long)matalgtobasis(nf,p1);
                    447:
                    448:   tetpil=avma;
                    449:   return gerepile(av,tetpil,nfhermite(nf,z));
                    450: }
                    451:
                    452: GEN
                    453: rnfidealhermite(GEN rnf,GEN x)
                    454: {
                    455:   long tx=typ(x),lx=lg(x),av=avma,tetpil,i,j,n,m;
                    456:   GEN z,p1,p2,x1,x2,x1j,nf,bas,unnf,zeronf;
                    457:
                    458:   checkrnf(rnf);
                    459:   n=lgef(rnf[1])-3; nf=(GEN)rnf[10]; bas=(GEN)rnf[7];
                    460:
                    461:   switch(tx)
                    462:   {
                    463:     case t_INT: case t_FRAC: case t_FRACN: z=cgetg(3,t_VEC);
                    464:       m=lgef(nf[1])-3; zeronf=gscalcol_i(gzero,m); unnf=gscalcol_i(gun,m);
                    465:       p1=cgetg(n+1,t_MAT); z[1]=(long)p1;
                    466:       for (j=1; j<=n; j++)
                    467:       {
                    468:        p2=cgetg(n+1,t_COL); p1[j]=(long)p2;
                    469:        for (i=1; i<=n; i++) p2[i]=(i==j)?(long)unnf:(long)zeronf;
                    470:       }
                    471:       z[2]=lmul(x,(GEN)bas[2]); return z;
                    472:
                    473:     case t_POLMOD: case t_POL:
                    474:       p1=rnfalgtobasis(rnf,x); tetpil=avma;
                    475:       return gerepile(av,tetpil,rnfprincipaltohermite(rnf,p1));
                    476:
                    477:     case t_VEC:
                    478:       switch(lx)
                    479:       {
                    480:        case 3:
                    481:          x1=(GEN)x[1];
                    482:          if (typ(x1)!=t_MAT || lg(x1) < n+1 || lg(x1[1]) != n+1)
                    483:            err(talker,"incorrect type in rnfidealhermite");
                    484:          p1=cgetg(n+1,t_MAT);
                    485:          for (j=1; j<=n; j++)
                    486:          {
                    487:            p2=cgetg(n+1,t_COL); p1[j]=(long)p2; x1j=(GEN)x1[j];
                    488:            for (i=1; i<=n; i++)
                    489:            {
                    490:               tx = typ(x1j[i]);
                    491:               if (is_const_t(tx)) p2[i] = x1j[i];
                    492:               else
                    493:                 switch(tx)
                    494:                 {
                    495:                   case t_POLMOD: case t_POL:
                    496:                     p2[i]=(long)algtobasis(nf,(GEN)x1j[i]); break;
                    497:                   case t_COL:
                    498:                     p2[i]=x1j[i]; break;
                    499:                   default: err(talker,"incorrect type in rnfidealhermite");
                    500:                 }
                    501:            }
                    502:          }
                    503:          x2=(GEN)x[2];
                    504:          if (typ(x2)!=t_VEC || lg(x2)!=lg(x1))
                    505:            err(talker,"incorrect type in rnfidealhermite");
                    506:          tetpil=avma; z=cgetg(3,t_VEC); z[1]=lcopy(p1); z[2]=lcopy(x2);
                    507:          z=gerepile(av,tetpil,nfhermite(nf,z));
                    508:          if (lg(z[1]) != n+1)
                    509:            err(talker,"not an ideal in rnfidealhermite");
                    510:          return z;
                    511:
                    512:        case 6:
                    513:          err(impl,"rnfidealhermite for prime ideals");
                    514:        default:
                    515:          err(typeer,"rnfidealhermite");
                    516:       }
                    517:
                    518:     case t_COL:
                    519:       if (lx!=(n+1)) err(typeer,"rnfidealhermite");
                    520:       return rnfprincipaltohermite(rnf,x);
                    521:
                    522:     case t_MAT:
                    523:       return rnfidealabstorel(rnf,x);
                    524:   }
                    525:   err(typeer,"rnfidealhermite");
                    526:   return NULL; /* not reached */
                    527: }
                    528:
                    529: GEN
                    530: rnfidealnormrel(GEN rnf,GEN id)
                    531: {
                    532:   long av=avma,i,n;
                    533:   GEN z,id2,nf;
                    534:
                    535:   checkrnf(rnf);
                    536:   id=rnfidealhermite(rnf,id); id2=(GEN)id[2];
                    537:   n=lgef(rnf[1])-3; nf=(GEN)rnf[10];
                    538:   if (n==1) { avma=av; return idmat(lgef(nf[1]-3)); }
                    539:   z=(GEN)id2[1]; for (i=2; i<=n; i++) z=idealmul(nf,z,(GEN)id2[i]);
                    540:   return gerepileupto(av,z);
                    541: }
                    542:
                    543: GEN
                    544: rnfidealnormabs(GEN rnf,GEN id)
                    545: {
                    546:   long av=avma,i,n;
                    547:   GEN z,id2;
                    548:
                    549:   checkrnf(rnf);
                    550:   id=rnfidealhermite(rnf,id); id2=(GEN)id[2];
                    551:   n=lgef(rnf[1])-3;
                    552:   z=gun; for (i=1; i<=n; i++) z=gmul(z,dethnf((GEN)id2[i]));
                    553:   return gerepileupto(av,z);
                    554: }
                    555:
                    556: GEN
                    557: rnfidealreltoabs(GEN rnf,GEN x)
                    558: {
                    559:   long av=avma,tetpil,i,j,k,n,m;
                    560:   GEN nf,basinv,om,id,p1,p2,p3,p4,p5;
                    561:
                    562:   checkrnf(rnf);
                    563:   x=rnfidealhermite(rnf,x);
                    564:   n=lgef(rnf[1])-3; nf=(GEN)rnf[10]; m=lgef(nf[1])-3;
                    565:   basinv=(GEN)((GEN)rnf[11])[5];
                    566:   p3=cgetg(n*m+1,t_MAT); p2=gmael(rnf,11,2);
                    567:   for (i=1; i<=n; i++)
                    568:   {
                    569:     om=rnfbasistoalg(rnf,gmael(x,1,i)); id=gmael(x,2,i);
                    570:     om=rnfelementreltoabs(rnf,om);
                    571:     for (j=1; j<=m; j++)
                    572:     {
                    573:       p1=gmul((GEN)nf[7],(GEN)id[j]);
                    574:       p4=lift_intern(gmul(om,gsubst(p1,varn(nf[1]),p2)));
                    575:       p5=cgetg(n*m+1,t_COL);
                    576:       for (k=0; k<n*m; k++) p5[k+1]=(long)truecoeff(p4,k);
                    577:       p3[(i-1)*m+j]=(long)p5;
                    578:     }
                    579:   }
                    580:   p1=gmul(basinv,p3); p2=detint(p1);
                    581:   tetpil=avma; return gerepile(av,tetpil,hnfmod(p1,p2));
                    582: }
                    583:
                    584: GEN
                    585: rnfidealabstorel(GEN rnf,GEN x)
                    586: {
                    587:   long av=avma,tetpil,n,m,j,k;
                    588:   GEN nf,basabs,ma,xj,p1,p2,id;
                    589:
                    590:   checkrnf(rnf); n=lgef(rnf[1])-3; nf=(GEN)rnf[10]; m=lgef(nf[1])-3;
                    591:   if (typ(x)!=t_MAT || lg(x)!=(n*m+1) || lg(x[1])!=(n*m+1))
                    592:     err(impl,"rnfidealabstorel for an ideal not in HNF");
                    593:   basabs=gmael(rnf,11,4); ma=cgetg(n*m+1,t_MAT);
                    594:   for (j=1; j<=n*m; j++)
                    595:   {
                    596:     p2=cgetg(n+1,t_COL); ma[j]=(long)p2;
                    597:     xj=gmul(basabs,(GEN)x[j]);
                    598:     xj=lift_intern(rnfelementabstorel(rnf,xj));
                    599:     for (k=0; k<n; k++)
                    600:       p2[k+1]=(long)truecoeff(xj,k);
                    601:   }
                    602:   ma=gmul((GEN)rnf[8],ma);
                    603:   ma=matalgtobasis(nf,ma);
                    604:   p1=cgetg(n*m+1,t_VEC); id=idmat(m);
                    605:   for (j=1; j<=n*m; j++) p1[j]=(long)id;
                    606:   p2=cgetg(3,t_VEC); p2[1]=(long)ma; p2[2]=(long)p1;
                    607:   tetpil=avma; return gerepile(av,tetpil,nfhermite(nf,p2));
                    608: }
                    609:
                    610: GEN
                    611: rnfidealdown(GEN rnf,GEN x)
                    612: {
                    613:   long av=avma,tetpil;
                    614:
                    615:   checkrnf(rnf); x=rnfidealhermite(rnf,x);
                    616:   tetpil=avma; return gerepile(av,tetpil,gcopy(gmael(x,2,1)));
                    617: }
                    618:
                    619: /* lift ideal x to the relative extension, returns a Z-basis */
                    620: GEN
                    621: rnfidealup(GEN rnf,GEN x)
                    622: {
                    623:   long av=avma,tetpil,i,n,m;
                    624:   GEN nf,bas,bas2,p1,p2,zeronf,unnf;
                    625:
                    626:   checkrnf(rnf);
                    627:   bas=(GEN)rnf[7]; bas2=(GEN)bas[2];
                    628:   n=lg(bas2)-1; nf=(GEN)rnf[10]; m=lgef((GEN)nf[1])-3;
                    629:   zeronf=zerocol(m); unnf=gscalcol_i(gun,m);
                    630:   p2=cgetg(3,t_VEC); p1=cgetg(n+1,t_VEC);
                    631:   p2[1]=(long)idmat_intern(n,unnf,zeronf);
                    632:   p2[2]=(long)p1;
                    633:   for (i=1; i<=n; i++) p1[i]=(long)idealmul(nf,x,(GEN)bas2[i]);
                    634:   tetpil=avma; return gerepile(av,tetpil,rnfidealreltoabs(rnf,p2));
                    635: }
                    636:
                    637: /* x a relative HNF ---> vector of 2 generators (relative polymods) */
                    638: GEN
                    639: rnfidealtwoelement(GEN rnf,GEN x)
                    640: {
                    641:   long av=avma,tetpil,j;
                    642:   GEN p1,p2,p3,res,polabs,nfabs,z;
                    643:
                    644:   res=(GEN)rnf[11]; polabs=(GEN)res[1];
                    645:   nfabs=cgetg(10,t_VEC); nfabs[1]=(long)polabs;
                    646:   for (j=2; j<=9; j++) nfabs[j]=zero;
                    647:   nfabs[7]=(long)lift((GEN)res[4]); nfabs[8]=res[5];
                    648:   p1=rnfidealreltoabs(rnf,x);
                    649:   p2=ideal_two_elt(nfabs,p1);
                    650:   p3=rnfelementabstorel(rnf,gmul((GEN)res[4],(GEN)p2[2]));
                    651:   tetpil=avma; z=cgetg(3,t_VEC); z[1]=lcopy((GEN)p2[1]);
                    652:   z[2]=(long)rnfalgtobasis(rnf,p3);
                    653:   return gerepile(av,tetpil,z);
                    654: }
                    655:
                    656: GEN
                    657: rnfidealmul(GEN rnf,GEN x,GEN y) /* x et y sous HNF relative uniquement */
                    658: {
                    659:   long av=avma,tetpil,i,j,n;
                    660:   GEN z,nf,x1,x2,p1,p2,p3,p4,p5,res;
                    661:
                    662:   z=rnfidealtwoelement(rnf,y);
                    663:   nf=(GEN)rnf[10]; n=lgef(rnf[1])-3;
                    664:   x=rnfidealhermite(rnf,x);
                    665:   x1=gmodulcp(gmul(gmael(rnf,7,1),matbasistoalg(nf,(GEN)x[1])),(GEN) rnf[1]);
                    666:   x2=(GEN)x[2]; p1=gmul((GEN)z[1],(GEN)x[1]);
                    667:   p2=lift_intern(gmul(rnfbasistoalg(rnf,(GEN)z[2]),x1));
                    668:   p3=cgetg(n+1,t_MAT);
                    669:   for (j=1; j<=n; j++)
                    670:   {
                    671:     p4=cgetg(n+1,t_COL); p3[j]=(long)p4; p5=(GEN)p2[j];
                    672:     for (i=1; i<=n; i++)
                    673:       p4[i]=(long)algtobasis(nf,truecoeff((GEN)p5,i-1));
                    674:   }
                    675:   res=cgetg(3,t_VEC);
                    676:   res[1]=(long)concatsp(p1,p3);
                    677:   res[2]=(long)concatsp(x2,x2);
                    678:   tetpil=avma; return gerepile(av,tetpil,nfhermite(nf,res));
                    679: }
                    680:
                    681: /*********************************************************************/
                    682: /**                                                                 **/
                    683: /**         LIBRARY FOR POLYNOMIALS WITH COEFFS. IN Z_K/P          **/
                    684: /**  An element in Z_K/P is a t_COL with degree(nf[1]) components.  **/
                    685: /**  These are integers modulo the prime p under prime ideal P      **/
                    686: /**  (only f(P/p) elements are non zero). These components are      **/
                    687: /**  given on the integer basis of K.                               **/
                    688: /**                                                                 **/
                    689: /*********************************************************************/
                    690:
                    691: /* K.B: What follows is not meant to work (yet?) */
                    692:
                    693: GEN
                    694: polnfmulscal(GEN nf,GEN s,GEN x)
                    695: {
                    696:   long i,lx=lgef(x);
                    697:   GEN z;
                    698:
                    699:   if (lx<3) return gcopy(x);
                    700:   if (gcmp0(s))
                    701:   {
                    702:     z=cgetg(2,t_POL); z[1]=evallgef(2) | evalvarn(varn(x));
                    703:     return z;
                    704:   }
                    705:   z=cgetg(lx,t_POL); z[1]=x[1];
                    706:   for (i=2; i<lx; i++) z[i]=(long)element_mul(nf,s,(GEN)x[i]);
                    707:   return z;
                    708: }
                    709:
                    710: GEN
                    711: polnfmul(GEN nf, GEN x, GEN y)
                    712: {
                    713:   long av,tetpil,m,i,d,imin,imax,lx,ly,lz;
                    714:   GEN p1,z,zeronf;
                    715:
                    716:   if (gcmp0(x)||gcmp0(y))
                    717:   {
                    718:     z=cgetg(2,t_POL); z[1]=evallgef(2) | evalvarn(varn(x));
                    719:     return z;
                    720:   }
                    721:   m=lgef(nf[1])-3; av=avma;
                    722:   lx=lgef(x)-3; ly=lgef(y)-3; lz=lx+ly;
                    723:   zeronf=gscalcol_i(gzero,m);
                    724:   z=cgetg(lz+3,t_POL);
                    725:   z[1] = evallgef(lz+3) | evalvarn(x) | evalsigne(1);
                    726:   for (d=0; d<=lz; d++)
                    727:   {
                    728:     p1=zeronf; imin=max(0,d-ly); imax=min(d,lx);
                    729:     for (i=imin; i<=imax; i++)
                    730:       p1=gadd(p1,element_mul(nf,(GEN)x[i+2],(GEN)y[d-i+2]));
                    731:     z[d+2]=(long)p1;
                    732:   }
                    733:   tetpil=avma; return gerepile(av,tetpil,gcopy(z));
                    734: }
                    735:
                    736: /* division euclidienne */
                    737: GEN
                    738: polnfdeuc(GEN nf, GEN x, GEN y, GEN *ptr)
                    739: {
                    740:   long av=avma,m,i,d,tx,lx,ly,lz,fl;
                    741:   GEN z,unnf,lcy,r;
                    742:   GEN *gptr[2];
                    743:
                    744:   if (gcmp0(y)) err(talker,"division by zero in polnfdiv");
                    745:   lx=lgef(x); ly=lgef(y); lz=lx-ly;
                    746:   if (gcmp0(x) || lz<0) { *ptr=gcopy(x); return zeropol(varn(x)); }
                    747:
                    748:   m=lgef(nf[1])-3; unnf=gscalcol_i(gun,m);
                    749:   x=dummycopy(x); y=dummycopy(y);
                    750:   for (i=2; i<lx; i++)
                    751:   {
                    752:     tx=typ(x[i]);
                    753:     if (is_intreal_t(tx) || tx == t_INTMOD || is_frac_t(tx))
                    754:       x[i]=lmul((GEN)x[i],unnf);
                    755:   }
                    756:   for (i=2; i<ly; i++)
                    757:   {
                    758:     tx=typ(y[i]);
                    759:     if (is_intreal_t(tx) || tx == t_INTMOD || is_frac_t(tx))
                    760:       y[i]=lmul((GEN)y[i],unnf);
                    761:   }
                    762:
                    763:   lz += 3;
                    764:   z=cgetg(lz,t_POL); z[1]=evallgef(lz) | evalvarn(x) | evalsigne(1);
                    765:   lcy=(GEN)y[ly-1];
                    766:   if (gegal(lift(lcy),unnf)) fl=0;
                    767:   else
                    768:   {
                    769:     fl=1; y=polnfmulscal(nf,element_inv(nf,lcy),y);
                    770:   }
                    771:   for (d=lz-1; d>=2; d--)
                    772:   {
                    773:     z[d]=x[d+ly-3];
                    774:     for (i=d; i<d+ly-3; i++)
                    775:       x[i]=lsub((GEN)x[i],element_mul(nf,(GEN)z[d],(GEN)y[i-d-2]));
                    776:   }
                    777:   if (fl) z=polnfmulscal(nf,lcy,z);
                    778:
                    779:   for(;;)
                    780:   {
                    781:     if (!gcmp0((GEN)x[d]))
                    782:     {
                    783:       r=cgetg(d,t_POL);
                    784:       r[1] = evallgef(d) | evalvarn(varn(x)) | evalsigne(1);
                    785:       for (i=2; i<d; i++) r[i]=x[i];
                    786:       break;
                    787:     }
                    788:     if (d==2) { r = zeropol(varn(x)); break; }
                    789:     d--;
                    790:   }
                    791:   *ptr=r; gptr[0]=ptr; gptr[1]=&z;
                    792:   gerepilemany(av,gptr,2); return z;
                    793: }
                    794:
                    795: GEN
                    796: polnfpow(GEN nf,GEN x,GEN k)
                    797: {
                    798:   long s,av=avma,m;
                    799:   GEN y,z;
                    800:
                    801:   m=lgef(nf[1])-3;
                    802:   if (typ(k)!=t_INT) err(talker,"not an integer exponent in nfpow");
                    803:   s=signe(k); if (s<0) err(impl,"polnfpow for negative exponents");
                    804:
                    805:   z=x; y=cgetg(3,t_POL);
                    806:   y[1] = evallgef(3) | evalvarn(varn(x)) | evalsigne(1);
                    807:   y[2] = (long)gscalcol_i(gun,m);
                    808:   for(;;)
                    809:   {
                    810:     if (mpodd(k)) y=polnfmul(nf,z,y);
                    811:     k=shifti(k,-1);
                    812:     if (!signe(k)) { cgiv(k); return gerepileupto(av,y); }
                    813:     z=polnfmul(nf,z,z);
                    814:   }
                    815: }

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