[BACK]Return to H.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

Diff for /OpenXM_contrib2/asir2000/engine/H.c between version 1.3 and 1.4

version 1.3, 2000/08/22 05:04:04 version 1.4, 2001/07/03 01:41:25
Line 45 
Line 45 
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,   * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.   * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *   *
  * $OpenXM: OpenXM_contrib2/asir2000/engine/H.c,v 1.2 2000/08/21 08:31:25 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/H.c,v 1.3 2000/08/22 05:04:04 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "inline.h"  #include "inline.h"
Line 79  void minipoly_mod(int,UM,UM,UM);
Line 79  void minipoly_mod(int,UM,UM,UM);
 int find_root(int,UM,int *);  int find_root(int,UM,int *);
 void showum(UM);  void showum(UM);
 void showumat(int **,int);  void showumat(int **,int);
   
   void henmain2(LUM,UM,UM,UM,UM,int,int,LUM *);
   void addtolum(int,int,LUM,LUM);
   void clearlum(int,int,LUM);
   void henprep2(int,int,int,UM,UM,UM,UM,UM,UM,UM);
   void henprep(P,ML,ML,ML *,ML *);
   
 #if 1  #if 1
 #define Mulum mulum  #define Mulum mulum
 #define Divum divum  #define Divum divum
Line 277  ML *listp;
Line 284  ML *listp;
         }          }
 }  }
   
   void hensel2(index,count,f,listp)
   int index,count;
   P f;
   ML *listp;
   {
           register int i,j;
           int mod,q,n,bound,dx;
           int *p;
           ML blist,clist,bqlist,cqlist,rlist;
           UM fm,qfm,gm,qgm,hm,qhm,qam,qbm,w;
           UM *b;
           LUM fl,tl;
           LUM *l;
           int dr,k;
   
           dx = UDEG(f);
           if ( dx == 1 ) {
                   *listp = blist = MLALLOC(1); blist->n = 1; blist->c[0] = 0;
                   return;
           }
           berle(index,count,f,&blist);
           n = blist->n;
           mod = blist->mod;
   
           if ( n == 1 ) {
                   *listp = blist = MLALLOC(1); blist->n = 1; blist->c[0] = 0;
                   return;
           }
   
           /* find k s.t. mod^k <= 2^27 < mod^(k+1); set q = mod^k */
           for ( k = 1, q = mod; q <= ((1<<27)/mod); q *= mod, k++ );
   
           /* mignotte bound */
           bound = mignotte(q,f);
   
           *listp = rlist = MLALLOC(n);
           rlist->n = n;
           rlist->mod = q;
           rlist->bound = bound;
   
           if ( bound == 1 ) {
                   gcdgen(f,blist,&clist);
                   henprep(f,blist,clist,&bqlist,&cqlist);
   
                   for ( i = 0, b = (UM *)bqlist->c; i < n; i++ ) {
                           COEF(rlist)[i] = tl = LUMALLOC(DEG(b[i]),1);
                           for ( j = 0; j <= DEG(tl); j++ )
                                           COEF(tl)[j][0] = COEF(b[i])[j];
                           COEF(rlist)[i] = tl;
                   }
           } else {
                   /* fl = f mod q */
                   fl = LUMALLOC(dx,bound);
                   ptolum(q,bound,f,fl);
                   /* fm = f mod mod */
                   fm = W_UMALLOC(dx);
                   ptoum(mod,f,fm);
                   /* fm = f mod q */
                   qfm = W_UMALLOC(dx);
                   ptoum(q,f,qfm);
   
                   gm = W_UMALLOC(dx); qgm = W_UMALLOC(dx);
                   hm = W_UMALLOC(dx); qhm = W_UMALLOC(dx);
                   qam = W_UMALLOC(dx); qbm = W_UMALLOC(dx);
                   w = W_UMALLOC(dx);
                   for ( i = 0; i < n-1; i++ ) {
                           cpyum(COEF(blist)[i],gm);
                           cpyum(fm,w);
                           divum(mod,w,gm,hm);
   
                           /* find am,bm s.t. qam*qgm+qbm*qhm=1 mod q, qgm=gm mod mod, qhm=hm mod mod */
                           henprep2(mod,q,k,qfm,gm,hm,qgm,qhm,qam,qbm);
   
                           henmain2(fl,qgm,qhm,qam,qbm,q,bound,&tl);
                           rlist->c[i] = (pointer)tl;
                           cpyum(hm,fm);
                           cpyum(qhm,qfm);
                   }
                   rlist->c[i] = fl;
           }
   }
   
   /*
           f = g0*h0 mod m -> f = gk*hk mod m^(bound), f is replaced by hk
   */
   
   void henmain2(f,g0,h0,a0,b0,m,bound,gp)
   LUM f;
   UM g0,h0,a0,b0;
   int m,bound;
   LUM *gp;
   {
           int n,dg,dh,i,k,j,dg1,dh1;
           UM wu,wr,ws,wt,q,wh1,wg1,wc,wd,we,wz,w1,w2;
           LUM wb0,wb1,wb2,fk,gk,hk;
   
           n = DEG(f); dg = DEG(g0); dh = DEG(h0);
   
           W_LUMALLOC(n,bound,wb0);
           W_LUMALLOC(n,bound,wb1);
           W_LUMALLOC(n,bound,wb2);
   
           wt = W_UMALLOC(2*n); ws = W_UMALLOC(2*n);
           wr = W_UMALLOC(2*n); wu = W_UMALLOC(2*n);
           q = W_UMALLOC(2*n);
           wh1 = W_UMALLOC(2*n); wg1 = W_UMALLOC(2*n);
   
           /* gk = g0 */
           gk = LUMALLOC(n,bound);
           DEG(gk) = dg;
           for ( i = 0; i <= dg; i++ )
                   COEF(gk)[i][0] = COEF(g0)[i];
   
           /* hk = h0 */
           W_LUMALLOC(n,bound,hk);
           DEG(hk) = dh;
           for ( i = 0; i <= dh; i++ )
                   COEF(hk)[i][0] = COEF(h0)[i];
   
           /* fk = gk*hk */
           W_LUMALLOC(n,bound,fk);
           mullum(m,bound,gk,hk,fk);
   
           wc = W_UMALLOC(2*n); wd = W_UMALLOC(2*n);
           we = W_UMALLOC(2*n); wz = W_UMALLOC(2*n);
   
   #if 0
           mulum(m,a0,g0,wc);
           mulum(m,b0,h0,wd);
           addum(m,wc,wd,wz);
           if ( DEG(wz) != 0 || COEF(wz)[0] != 1 )
                   error("henmain2 : cannot happen(extgcd)");
   #endif
   
           fprintf(stderr,"bound=%d\n",bound);
           for ( k = 1; k < bound; k++ ) {
                   fprintf(stderr,".");
                   /* at this point, f = gk*hk mod y^k */
   
   #if 0
                   for ( j = 0; j < k; j++ )
                           for ( i = 0; i <= n; i++ )
                                   if ( COEF(f)[i][j] != COEF(f)[i][j] )
                                           error("henmain2 : cannot happen(f=fk)");
   #endif
   
                   /* wt = (f-gk*hk)/y^k */
                   for ( i = 0; i <= n; i++ )
                           COEF(ws)[i] = COEF(f)[i][k];
                   degum(ws,n);
                   for ( i = 0; i <= n; i++ )
                           COEF(wu)[i] = COEF(fk)[i][k];
                   degum(wu,n);
                   subum(m,ws,wu,wt);
   
                   /* compute wf1,wg1 s.t. wh1*g0+wg1*h0 = wt */
                   mulum(m,a0,wt,wh1); DEG(wh1) = divum(m,wh1,h0,q);
                   mulum(m,wh1,g0,wc); subum(m,wt,wc,wd); DEG(wd) = divum(m,wd,h0,wg1);
   
                   /* check */
   #if 0
                   if ( DEG(wd) >= 0 || DEG(wg1) > dg )
                           error("henmain2 : cannot happen(adj)");
   
                   mulum(m,wg1,h0,wc); mulum(m,wh1,g0,wd); addum(m,wc,wd,we);
                   subum(m,we,wt,wz);
                   if ( DEG(wz) >= 0 )
                           error("henmain2 : cannot happen(coef)");
   #endif
   
                   /* fk += ((wg1*hk+wh1*gk)*y^k+wg1*wh1*y^(2*k) mod m^bound */
   
                   /* wb0 = wh1*y^k */
                   clearlum(n,bound,wb0);
                   DEG(wb0) = dh1 = DEG(wh1);
                   for ( i = 0; i <= dh1; i++ )
                           COEF(wb0)[i][k] = COEF(wh1)[i];
   
                   /* wb2 = gk*wb0 mod y^bound */
                   clearlum(n,bound,wb2);
                   mullum(m,bound,gk,wb0,wb2);
   
                   /* fk += wb2 */
                   addtolum(m,bound,wb2,fk);
   
                   /* wb1 = wg1*y^k */
                   clearlum(n,bound,wb1);
                   DEG(wb1) = dg1 = DEG(wg1);
                   for ( i = 0; i <= n; i++ )
                           COEF(wb1)[i][k] = COEF(wg1)[i];
   
                   /* wb2 = hk*wb1 mod y^bound */
                   clearlum(n,bound,wb2);
                   mullum(m,bound,hk,wb1,wb2);
   
                   /* fk += wb2 */
                   addtolum(m,bound,wb2,fk);
   
                   /* fk += wg1*wh1*y^(2*k) mod y^bound) */
                   if ( 2*k < bound ) {
                           clearlum(n,bound,wb2);
                           mullum(m,bound,wb0,wb1,wb2);
                           addtolum(m,bound,wb2,fk);
                   }
   
                   /* gk += wg1*y^k, hk += wh1*y^k */
                   for ( i = 0; i <= DEG(wg1); i++ )
                           COEF(gk)[i][k] = COEF(wg1)[i];
                   for ( i = 0; i <= DEG(wh1); i++ )
                           COEF(hk)[i][k] = COEF(wh1)[i];
           }
           fprintf(stderr,"\n");
           *gp = gk;
           clearlum(n,bound,f);
           DEG(f) = dh;
           for ( i = 0; i <= dh; i++ )
                   for ( j = 0; j < bound; j++ )
                           COEF(f)[i][j] = COEF(hk)[i][j];
   }
   
   void clearlum(n,bound,f)
   int n,bound;
   LUM f;
   {
           int i;
   
           for ( i = 0; i <= n; i++ )
                   bzero(COEF(f)[i],bound*sizeof(int));
   }
   
   /* g += f */
   
   void addtolum(m,bound,f,g)
   int m,bound;
   LUM f;
   LUM g;
   {
           int n,i;
   
           n = DEG(f);
           for ( i = 0; i <= n; i++ )
                   addpadic(m,bound,COEF(f)[i],COEF(g)[i]);
   }
   
 void hsq(index,count,f,nindex,dcp)  void hsq(index,count,f,nindex,dcp)
 int index,count,*nindex;  int index,count,*nindex;
 P f;  P f;
Line 435  ML blist,*clistp;
Line 686  ML blist,*clistp;
         }          }
 }  }
   
   /* find a,b s.t. qa*qg+qb*qh=1 mod q, qg=g mod mod, qh=h mod mod */
   /* q = mod^k */
   
   void henprep2(mod,q,k,f,g,h,qg,qh,qa,qb)
   int mod,q,k;
   UM f,g,h,qg,qh,qa,qb;
   {
           int n;
           UM wg,wh,wa,wb;
           UM wt,ws,wu;
           ML bl,cl,bql,cql;
           P ff;
   
           n = DEG(f);
           wg = W_UMALLOC(2*n); wh = W_UMALLOC(2*n);
           wa = W_UMALLOC(2*n); wb = W_UMALLOC(2*n);
           cpyum(g,wg); cpyum(h,wh);
   
           /* wa*g+wb*h = 1 mod mod */
           eucum(mod,wg,wh,wa,wb);
   
   #if 0
           /* check */
           wt = W_UMALLOC(2*n); ws = W_UMALLOC(2*n); wu = W_UMALLOC(2*n);
           mulum(mod,wa,g,wt);
           mulum(mod,wb,h,ws);
           addum(mod,wt,ws,wu);
           if ( DEG(wu) != 0 || COEF(wu)[0] != 1 )
                   error("henprep 1");
   #endif
   
           bl = MLALLOC(2); bl->n = 2; bl->mod = mod; bl->c[0] = g; bl->c[1] = h;
           cl = MLALLOC(2); cl->n = 2; cl->mod = mod; cl->c[0] = wb; cl->c[1] = wa;
           umtop(CO->v,f,&ff); /* XXX */
           henprep(ff,bl,cl,&bql,&cql); /* XXX */
   
           cpyum(bql->c[0],qg); cpyum(bql->c[1],qh);
           cpyum(cql->c[0],qb); cpyum(cql->c[1],qa);
   
   #if 0
           /* check */
           mulum(q,qa,qg,wt);
           mulum(q,qb,qh,ws);
           addum(q,wt,ws,wu);
           if ( DEG(wu) != 0 || COEF(wu)[0] != 1 )
                   error("henprep 2");
   #endif
   }
   
 /*  /*
         henprep(f,&blist,&clist,&bqlist,&cqlist);          henprep(f,blist,clist,&bqlist,&cqlist);
  */   */
   
 void henprep(f,blist,clist,bqlistp,cqlistp)  void henprep(f,blist,clist,bqlistp,cqlistp)
Line 549  ML bqlist,cqlist,*listp;
Line 849  ML bqlist,cqlist,*listp;
                 for ( j = DEG(b[i]), pp = COEF(l[i]), px = COEF(b[i]); j >= 0; j-- )                  for ( j = DEG(b[i]), pp = COEF(l[i]), px = COEF(b[i]); j >= 0; j-- )
                         pp[j][0] = px[j];                          pp[j][0] = px[j];
         }          }
           fprintf(stderr,"bound=%d\n",bound);
         for ( i = 1; i < bound; i++ ) {          for ( i = 1; i < bound; i++ ) {
                   fprintf(stderr,".");
                 mullum(mod,i+1,l[0],l[1],wb0);                  mullum(mod,i+1,l[0],l[1],wb0);
                 for ( j = 2; j < np; j++ ) {                  for ( j = 2; j < np; j++ ) {
                         mullum(mod,i+1,l[j],wb0,wb1);                          mullum(mod,i+1,l[j],wb0,wb1);
Line 578  ML bqlist,cqlist,*listp;
Line 880  ML bqlist,cqlist,*listp;
                 for ( j = n, px = COEF(wq0); j >= 0; j-- )                  for ( j = n, px = COEF(wq0); j >= 0; j-- )
                         px[j] = 0;                          px[j] = 0;
         }          }
           fprintf(stderr,"\n");
 }  }
   
 static double M;  static double M;
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>