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

Diff for /OpenXM_contrib2/asir2000/engine/A.c between version 1.5 and 1.8

version 1.5, 2001/10/09 01:36:09 version 1.8, 2018/03/29 01:32:51
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/A.c,v 1.4 2001/06/07 04:54:39 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/A.c,v 1.7 2013/06/13 18:40:31 ohara Exp $
 */  */
 #include "ca.h"  #include "ca.h"
   #include "parse.h"
   
 int get_lprime();  int get_lprime();
   
Line 56  VL vl;
Line 57  VL vl;
 P p1,p2,p0;  P p1,p2,p0;
 P *pr;  P *pr;
 {  {
         Q m,d;    Q m,d;
         P t,q,r;    P t,q,r;
   
         if ( VR(p2) == VR(p0) ) {    if ( VR(p2) == VR(p0) ) {
                 *pr = p1;      *pr = p1;
                 return;      return;
         }    }
         if ( VR(p1) == VR(p0) ) {    if ( VR(p1) == VR(p0) ) {
                 error("pdiva : ???");      error("pdiva : ???");
                 return;      return;
         }    }
         subq(DEG(DC(p1)),DEG(DC(p2)),&m); addq(m,ONE,&d);    subq(DEG(DC(p1)),DEG(DC(p2)),&m); addq(m,ONE,&d);
         pwrq((Q)COEF(DC(p2)),d,&m); mulp(vl,p1,(P)m,&t);    pwrq((Q)COEF(DC(p2)),d,&m); mulp(vl,p1,(P)m,&t);
         divsrp(vl,t,p2,&q,&r); remsdcp(vl,q,p0,pr);    divsrp(vl,t,p2,&q,&r); remsdcp(vl,q,p0,pr);
 }  }
   
 void rema(vl,p0,p1,p2,pr)  void rema(vl,p0,p1,p2,pr)
 VL vl;  VL vl;
 P p1,p2,p0,*pr;  P p1,p2,p0,*pr;
 {  {
         P m,m1,m2,g,lc,x,t;    P m,m1,m2,g,lc,x,t;
         Q q;    Q q;
         V v;    V v;
         int d1,d2;    int d1,d2;
   
         if ( !p1 || !p2 || NUM(p2) || ( VR(p2) == VR(p0) ) )    if ( !p1 || !p2 || NUM(p2) || ( VR(p2) == VR(p0) ) )
                 *pr = 0;      *pr = 0;
         else if ( NUM(p1) )    else if ( NUM(p1) )
                 *pr = p1;      *pr = p1;
         else if ( ( v = VR(p1) ) == VR(p0) )    else if ( ( v = VR(p1) ) == VR(p0) )
                 remsdcp(vl,p1,p0,pr);      remsdcp(vl,p1,p0,pr);
         else if ( ( d1 = deg(v,p1) ) < ( d2 = deg(v,p2) ) )    else if ( ( d1 = deg(v,p1) ) < ( d2 = deg(v,p2) ) )
                 *pr = p1;      *pr = p1;
         else if ( d1 == d2 ) {    else if ( d1 == d2 ) {
                 mulp(vl,p1,LC(p2),&m); mulp(vl,p2,LC(p1),&m1);      mulp(vl,p1,LC(p2),&m); mulp(vl,p2,LC(p1),&m1);
                 subp(vl,m,m1,&m2);      subp(vl,m,m1,&m2);
                 remsdcp(vl,m2,p0,pr);      remsdcp(vl,m2,p0,pr);
         } else {    } else {
                 g = p1; lc = LC(p2); MKV(v,x);      g = p1; lc = LC(p2); MKV(v,x);
                 while ( ( d1 = deg(v,g) ) >= d2 ) {      while ( ( d1 = deg(v,g) ) >= d2 ) {
                         mulp(vl,g,lc,&m); mulp(vl,p2,LC(g),&m1);        mulp(vl,g,lc,&m); mulp(vl,p2,LC(g),&m1);
                         STOQ(d1-d2,q); pwrp(vl,x,q,&t);        STOQ(d1-d2,q); pwrp(vl,x,q,&t);
                         mulp(vl,m1,t,&m2); subp(vl,m,m2,&m1);        mulp(vl,m1,t,&m2); subp(vl,m,m2,&m1);
                         remsdcp(vl,m1,p0,&g);        remsdcp(vl,m1,p0,&g);
                 }      }
                 *pr = g;      *pr = g;
         }    }
 }  }
   
 void gcda(vl,p0,p1,p2,pr)  void gcda(vl,p0,p1,p2,pr)
 VL vl;  VL vl;
 P p1,p2,p0,*pr;  P p1,p2,p0,*pr;
 {  {
         P g1,g2,r,t,c;    P g1,g2,r,t,c;
   
         remsdcp(vl,p1,p0,&t); pcp(vl,t,&g1,&c); remsdcp(vl,p2,p0,&g2);    remsdcp(vl,p1,p0,&t); pcp(vl,t,&g1,&c); remsdcp(vl,p2,p0,&g2);
         while ( 1 ) {    while ( 1 ) {
                 if ( NUM(g2) || (VR(g2) == VR(p0)) ) {      if ( NUM(g2) || (VR(g2) == VR(p0)) ) {
                         *pr = (P)ONE;        *pr = (P)ONE;
                         return;        return;
                 }      }
                 pcp(vl,g2,&t,&c); pmonic(vl,p0,t,&g2); rema(vl,p0,g1,g2,&r);      pcp(vl,g2,&t,&c); pmonic(vl,p0,t,&g2); rema(vl,p0,g1,g2,&r);
                 if ( !r ) {      if ( !r ) {
                         pmonic(vl,p0,g2,pr);        pmonic(vl,p0,g2,pr);
                         return;        return;
                 }      }
                 g1 = g2; g2 = r;      g1 = g2; g2 = r;
         }    }
 }  }
   
 void sqa(vl,p0,p,dcp)  void sqa(vl,p0,p,dcp)
Line 131  VL vl;
Line 132  VL vl;
 P p,p0;  P p,p0;
 DCP *dcp;  DCP *dcp;
 {  {
         V v;    V v;
         int i,j,n;    int i,j,n;
         P f,f1,f2,gcd,flat,q,r,t;    P f,f1,f2,gcd,flat,q,r,t;
         P *a;    P *a;
         int *b;    int *b;
         DCP dc,dc0;    DCP dc,dc0;
   
         if ( ( v = VR(p) ) == VR(p0) ) {    if ( ( v = VR(p) ) == VR(p0) ) {
                 NEWDC(dc); *dcp = dc; COEF(dc) = (P)ONE; DEG(dc) = ONE;      NEWDC(dc); *dcp = dc; COEF(dc) = (P)ONE; DEG(dc) = ONE;
                 NEXT(dc) = 0;      NEXT(dc) = 0;
         } else {    } else {
                 n = deg(v,p); W_CALLOC(n,P,a); W_CALLOC(n,int,b);      n = deg(v,p); W_CALLOC(n,P,a); W_CALLOC(n,int,b);
                 for ( i = 0, f = p; ;i++ ) {      for ( i = 0, f = p; ;i++ ) {
                         if ( VR(f) != v )        if ( VR(f) != v )
                                 break;          break;
                         remsdcp(vl,f,p0,&f1); diffp(vl,f1,VR(f1),&f2);        remsdcp(vl,f,p0,&f1); diffp(vl,f1,VR(f1),&f2);
                         while ( 1 ) {        while ( 1 ) {
                                 pmonic(vl,p0,f2,&t); f2 = t; rema(vl,p0,f1,f2,&r);          pmonic(vl,p0,f2,&t); f2 = t; rema(vl,p0,f1,f2,&r);
                                 if ( !r ) {          if ( !r ) {
                                         pmonic(vl,p0,f2,&gcd); pqra(vl,p0,f,gcd,&flat,&r);            pmonic(vl,p0,f2,&gcd); pqra(vl,p0,f,gcd,&flat,&r);
                                         break;            break;
                                 }          }
                                 f1 = f2; f2 = r;          f1 = f2; f2 = r;
                         }        }
                         if ( VR(gcd) != v ) {        if ( VR(gcd) != v ) {
                                 a[i] = f; b[i] = 1;          a[i] = f; b[i] = 1;
                                 break;          break;
                         }        }
                         for ( j = 1, f = gcd; ; j++ ) {        for ( j = 1, f = gcd; ; j++ ) {
                                 pqra(vl,p0,f,flat,&q,&r);          pqra(vl,p0,f,flat,&q,&r);
                                 if ( r )          if ( r )
                                         break;            break;
                                 else          else
                                         f = q;            f = q;
                         }        }
                         a[i] = flat; b[i] = j;        a[i] = flat; b[i] = j;
                 }      }
                 for ( i = 0, j = 0, dc0 = 0; a[i]; i++ ) {      for ( i = 0, j = 0, dc0 = 0; a[i]; i++ ) {
                         NEXTDC(dc0,dc); j += b[i]; STOQ(j,DEG(dc));        NEXTDC(dc0,dc); j += b[i]; STOQ(j,DEG(dc));
                         if ( a[i+1] )        if ( a[i+1] )
                                 pqra(vl,p0,a[i],a[i+1],&COEF(dc),&t);          pqra(vl,p0,a[i],a[i+1],&COEF(dc),&t);
                         else        else
                                 COEF(dc) = a[i];          COEF(dc) = a[i];
                 }      }
                 NEXT(dc) = 0; *dcp = dc0;      NEXT(dc) = 0; *dcp = dc0;
         }    }
 }  }
   
 void pqra(vl, p0, p1, p2, p, pr)  void pqra(vl, p0, p1, p2, p, pr)
Line 184  VL vl;
Line 185  VL vl;
 P p0,p1,p2;  P p0,p1,p2;
 P *p,*pr;  P *p,*pr;
 {  {
         V v,v1;    V v,v1;
         P  m,m1,q,r;    P  m,m1,q,r;
         Q tq;    Q tq;
         int n;    int n;
   
         if ( ( v = VR(p1) ) != ( v1 = VR(p2) ) ) {    if ( ( v = VR(p1) ) != ( v1 = VR(p2) ) ) {
                 while ( ( VR(vl) != v ) && ( VR(vl) != v1 ) )      while ( ( VR(vl) != v ) && ( VR(vl) != v1 ) )
                         vl = NEXT(vl);        vl = NEXT(vl);
                 if ( VR(vl) == v ) {      if ( VR(vl) == v ) {
                         *p = p1; *pr = 0;        *p = p1; *pr = 0;
                 } else {      } else {
                         *p = 0; *pr = p1;        *p = 0; *pr = p1;
                 }      }
         } else if ( (n = deg(v,p1) - deg(v,p2)) < 0 ) {    } else if ( (n = deg(v,p1) - deg(v,p2)) < 0 ) {
                 *pr = p1; *p = 0;      *pr = p1; *p = 0;
         } else {    } else {
                 n++; STOQ(n,tq);      n++; STOQ(n,tq);
                 pwrp(vl,LC(p2),tq,&m); mulp(vl,m,p1,&m1); divsrp(vl,m1,p2,&q,&r);      pwrp(vl,LC(p2),tq,&m); mulp(vl,m,p1,&m1); divsrp(vl,m1,p2,&q,&r);
                 remsdcp(vl,q,p0,p); remsdcp(vl,r,p0,pr);      remsdcp(vl,q,p0,p); remsdcp(vl,r,p0,pr);
         }    }
 }  }
   
 void pmonic(vl, p0, p, pr)  void pmonic(vl, p0, p, pr)
 VL vl;  VL vl;
 P p,p0,*pr;  P p,p0,*pr;
 {  {
         P d,m,pp,c,t;    P d,m,pp,c,t;
   
         if ( NUM(p) || ( VR(p) == VR(p0) ) )    if ( NUM(p) || ( VR(p) == VR(p0) ) )
                 *pr = (P)ONE;      *pr = (P)ONE;
         else if ( NUM(COEF(DC(p))) )    else if ( NUM(COEF(DC(p))) )
                 *pr = p;      *pr = p;
         else {    else {
                 ptozp(COEF(DC(p)),1,(Q *)&c,&pp); pinva(p0,pp,&d);      ptozp(COEF(DC(p)),1,(Q *)&c,&pp); pinva(p0,pp,&d);
                 mulp(vl,p,d,&m); remsdcp(vl,m,p0,&t); ptozp(t,1,(Q *)&c,pr);      mulp(vl,p,d,&m); remsdcp(vl,m,p0,&t); ptozp(t,1,(Q *)&c,pr);
         }    }
 }  }
   
 void remsdcp(vl, p, p0, pr)  void remsdcp(vl, p, p0, pr)
 VL vl;  VL vl;
 P p,p0,*pr;  P p,p0,*pr;
 {  {
         P q,r,c;    P q,r,c;
         DCP dc,dcr,dcr0;    DCP dc,dcr,dcr0;
   
         if ( !p )    if ( !p )
                 *pr = 0;      *pr = 0;
         else if ( NUM(p) )    else if ( NUM(p) )
                 *pr = (P)ONE;      *pr = (P)ONE;
         else if ( VR(p) == VR(p0) ) {    else if ( VR(p) == VR(p0) ) {
                 divsrp(vl,p,p0,&q,&r); ptozp(r,1,(Q *)&c,pr);      divsrp(vl,p,p0,&q,&r); ptozp(r,1,(Q *)&c,pr);
         } else {    } else {
                 for ( dcr = dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {      for ( dcr = dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
                         divsrp(vl,COEF(dc),p0,&q,&r);        divsrp(vl,COEF(dc),p0,&q,&r);
                         if ( r ) {        if ( r ) {
                                 NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = r;          NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = r;
                         }        }
                 }      }
                 if ( !dcr0 )      if ( !dcr0 )
                         *pr = 0;        *pr = 0;
                 else {      else {
                         NEXT(dcr) = 0; MKP(VR(p),dcr0,r); ptozp(r,1,(Q *)&c,pr);        NEXT(dcr) = 0; MKP(VR(p),dcr0,r); ptozp(r,1,(Q *)&c,pr);
                 }      }
         }    }
 }  }
   
 void pinva(p0, p, pr)  void pinva(p0, p, pr)
 P p,p0,*pr;  P p,p0,*pr;
 {  {
         V v;    V v;
         Q f,q,q1;    Q f,q,q1;
         Q u,t,a,b,s;    Q u,t,a,b,s;
         P c,c1;    P c,c1;
         P tg,th,tr;    P tg,th,tr;
         P z,zz,zzz;    P z,zz,zzz;
         UM wg,wh,ws,wt,wm;    UM wg,wh,ws,wt,wm;
         int n,m,modres,mod,index,i;    int n,m,modres,mod,index,i;
   
         v = VR(p); n=deg(v,p); m=deg(v,p0);    v = VR(p); n=deg(v,p); m=deg(v,p0);
         norm(p,&a); norm(p0,&b);    norm(p,&a); norm(p0,&b);
         STOQ(m,u); pwrq(a,u,&t); STOQ(n,u); pwrq(b,u,&s);    STOQ(m,u); pwrq(a,u,&t); STOQ(n,u); pwrq(b,u,&s);
         mulq(t,s,&u);    mulq(t,s,&u);
         factorial(n+m,&t); mulq(u,t,&s); addq(s,s,&f);    factorial(n+m,&t); mulq(u,t,&s); addq(s,s,&f);
         wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n);    wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n);
         wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n);    wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n);
         wm = W_UMALLOC(m+n);    wm = W_UMALLOC(m+n);
         for ( q = ONE, t = 0, c = 0, index = 0; ; ) {    for ( q = ONE, t = 0, c = 0, index = 0; ; ) {
                 mod = get_lprime(index++);      mod = get_lprime(index++);
                 if ( !rem(NM((Q)LC(p)),mod) || !rem(NM((Q)LC(p0)),mod) )      if ( !rem(NM((Q)LC(p)),mod) || !rem(NM((Q)LC(p0)),mod) )
                         continue;        continue;
                 ptomp(mod,p,&tg); ptomp(mod,p0,&th); srchump(mod,tg,th,&tr);      ptomp(mod,p,&tg); ptomp(mod,p0,&th); srchump(mod,tg,th,&tr);
                 if ( !tr )      if ( !tr )
                         continue;        continue;
                 else      else
                         modres = CONT((MQ)tr);        modres = CONT((MQ)tr);
                 mptoum(tg,wg); mptoum(th,wh);      mptoum(tg,wg); mptoum(th,wh);
                 eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */      eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */
                 DEG(wm) = 0; COEF(wm)[0] = modres;      DEG(wm) = 0; COEF(wm)[0] = modres;
                 mulum(mod,ws,wm,wt);      mulum(mod,ws,wm,wt);
                 for ( i = DEG(wt); i >= 0; i-- )      for ( i = DEG(wt); i >= 0; i-- )
                         if ( ( COEF(wt)[i] * 2 ) > mod )        if ( ( COEF(wt)[i] * 2 ) > mod )
                                 COEF(wt)[i] -= mod;          COEF(wt)[i] -= mod;
                 chnrem(mod,v,c,q,wt,&c1,&q1);      chnrem(mod,v,c,q,wt,&c1,&q1);
                 if ( !ucmpp(c,c1) ) {      if ( !ucmpp(c,c1) ) {
                         mulp(CO,c,p,&z); divsrp(CO,z,p0,&zz,&zzz);        mulp(CO,c,p,&z); divsrp(CO,z,p0,&zz,&zzz);
                         if ( NUM(zzz) ) {        if ( NUM(zzz) ) {
                                 q = q1; c = c1; break;          q = q1; c = c1; break;
                         }        }
                 }      }
                 q = q1; c = c1;      q = q1; c = c1;
                 if ( cmpq(f,q) < 0 )      if ( cmpq(f,q) < 0 )
                         break;        break;
         }    }
         ptozp(c,1,&s,pr);    ptozp(c,1,&s,pr);
   }
   
   /*
    * compute g s.t. h*p0+g*p = res(v,p0,p)
    */
   
   void inva_chrem(P p0, P p, NODE *pr)
   {
     V v;
     Q f,q,q1;
     Q u,t,a,b,s,dn;
     P c,c0,c1;
     P tg,th,tr;
     P z,zz,zzz;
     UM wg,wh,ws,wt,wm;
     int n,m,modres,mod,index,i;
   
     v = VR(p); n=deg(v,p); m=deg(v,p0);
     norm(p,&a); norm(p0,&b);
     STOQ(m,u); pwrq(a,u,&t); STOQ(n,u); pwrq(b,u,&s);
     mulq(t,s,&u);
     factorial(n+m,&t); mulq(u,t,&s); addq(s,s,&f);
     wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n);
     wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n);
     wm = W_UMALLOC(m+n);
     for ( q = ONE, t = 0, c = 0, index = 0; ; ) {
       mod = get_lprime(index++);
       if ( !rem(NM((Q)LC(p)),mod) || !rem(NM((Q)LC(p0)),mod) )
         continue;
       ptomp(mod,p,&tg); ptomp(mod,p0,&th); srchump(mod,tg,th,&tr);
       if ( !tr )
         continue;
       else
         modres = CONT((MQ)tr);
       mptoum(tg,wg); mptoum(th,wh);
       eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */
       DEG(wm) = 0; COEF(wm)[0] = modres;
       mulum(mod,ws,wm,wt);
       for ( i = DEG(wt); i >= 0; i-- )
         if ( ( COEF(wt)[i] * 2 ) > mod )
           COEF(wt)[i] -= mod;
       chnrem(mod,v,c,q,wt,&c1,&q1);
       if ( !ucmpp(c,c1) ) {
         mulp(CO,c,p,&z); divsrp(CO,z,p0,&zz,&zzz);
         if ( NUM(zzz) ) {
           q = (Q)zzz; break;
         }
       }
       q = q1; c = c1;
       if ( cmpq(f,q) < 0 ) {
         mulp(CO,c,p,&z); divsrp(CO,z,p0,&zz,&zzz);
         q = (Q)zzz;
         break;
       }
     }
     /* c*p = q mod p0 */
     /* c = s*c0 */
     ptozp(c,1,&s,&c0);
     /* c/q = c0/(q/s) */
     divq(q,s,&dn);
     *pr = (NODE)mknode(2,c0,dn);
 }  }
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.8
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>