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

Diff for /OpenXM_contrib2/asir2000/engine/Q.c between version 1.9 and 1.10

version 1.9, 2002/08/14 04:49:52 version 1.10, 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/Q.c,v 1.8 2002/08/08 10:57:01 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/Q.c,v 1.9 2002/08/14 04:49:52 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "base.h"  #include "base.h"
Line 53 
Line 53 
   
 void addq(Q n1,Q n2,Q *nr)  void addq(Q n1,Q n2,Q *nr)
 {  {
         N nm,dn,nm1,nm2,nm3,dn0,dn1,dn2,g,g1,g0,m;    N nm,dn,nm1,nm2,nm3,dn0,dn1,dn2,g,g1,g0,m;
         int sgn;    int sgn;
         unsigned t,t1;    unsigned t,t1;
   
         if ( !n1 )    if ( !n1 )
                 *nr = n2;      *nr = n2;
         else if ( !n2 )    else if ( !n2 )
                 *nr = n1;      *nr = n1;
         else if ( INT(n1) )    else if ( INT(n1) )
                 if ( INT(n2) ) {      if ( INT(n2) ) {
                         nm1 = NM(n1); nm2 = NM(n2);        nm1 = NM(n1); nm2 = NM(n2);
                         if ( SGN(n1) == SGN(n2) ) {        if ( SGN(n1) == SGN(n2) ) {
                                 if ( PL(nm1) == 1 && PL(nm2) == 1 ) {          if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
                                         t1 = BD(nm1)[0]; t = t1+BD(nm2)[0];            t1 = BD(nm1)[0]; t = t1+BD(nm2)[0];
                                         if ( t < t1 ) {            if ( t < t1 ) {
                                                 m = NALLOC(2); PL(m) = 2; BD(m)[0] = t; BD(m)[1] = 1;              m = NALLOC(2); PL(m) = 2; BD(m)[0] = t; BD(m)[1] = 1;
                                         } else            } else
                                                 UTON(t,m);              UTON(t,m);
                                 } else          } else
                                         addn(NM(n1),NM(n2),&m);            addn(NM(n1),NM(n2),&m);
                                 NTOQ(m,SGN(n1),*nr);          NTOQ(m,SGN(n1),*nr);
                         } else {        } else {
                                 if ( PL(nm1) == 1 && PL(nm2) == 1 ) {          if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
                                         t1 = BD(nm1)[0]; t = t1-BD(nm2)[0];            t1 = BD(nm1)[0]; t = t1-BD(nm2)[0];
                                         if ( !t )            if ( !t )
                                                 sgn = 0;              sgn = 0;
                                         else if ( t > t1 ) {            else if ( t > t1 ) {
                                                 sgn = -1; t = -((int)t); UTON(t,m);              sgn = -1; t = -((int)t); UTON(t,m);
                                         } else {            } else {
                                                 sgn = 1; UTON(t,m);              sgn = 1; UTON(t,m);
                                         }            }
                                 } else          } else
                                         sgn = subn(NM(n1),NM(n2),&m);            sgn = subn(NM(n1),NM(n2),&m);
                                 if ( sgn ) {          if ( sgn ) {
                                         if ( SGN(n1) == sgn )            if ( SGN(n1) == sgn )
                                                 NTOQ(m,1,*nr);              NTOQ(m,1,*nr);
                                         else            else
                                                 NTOQ(m,-1,*nr);              NTOQ(m,-1,*nr);
                                 } else          } else
                                         *nr = 0;            *nr = 0;
                         }        }
                 } else {      } else {
                         kmuln(NM(n1),DN(n2),&m);        kmuln(NM(n1),DN(n2),&m);
                         if ( SGN(n1) == SGN(n2) ) {        if ( SGN(n1) == SGN(n2) ) {
                                 sgn = SGN(n1); addn(m,NM(n2),&nm);          sgn = SGN(n1); addn(m,NM(n2),&nm);
                         } else        } else
                                 sgn = SGN(n1)*subn(m,NM(n2),&nm);          sgn = SGN(n1)*subn(m,NM(n2),&nm);
                         if ( sgn ) {        if ( sgn ) {
                                 dn = DN(n2); NDTOQ(nm,dn,sgn,*nr);          dn = DN(n2); NDTOQ(nm,dn,sgn,*nr);
                         } else        } else
                                 *nr = 0;          *nr = 0;
                 }      }
         else if ( INT(n2) ) {    else if ( INT(n2) ) {
                 kmuln(NM(n2),DN(n1),&m);      kmuln(NM(n2),DN(n1),&m);
                 if ( SGN(n1) == SGN(n2) ) {      if ( SGN(n1) == SGN(n2) ) {
                         sgn = SGN(n1); addn(m,NM(n1),&nm);        sgn = SGN(n1); addn(m,NM(n1),&nm);
                 } else      } else
                         sgn = SGN(n1)*subn(NM(n1),m,&nm);        sgn = SGN(n1)*subn(NM(n1),m,&nm);
                 if ( sgn ) {      if ( sgn ) {
                         dn = DN(n1); NDTOQ(nm,dn,sgn,*nr);        dn = DN(n1); NDTOQ(nm,dn,sgn,*nr);
                 } else      } else
                         *nr = 0;        *nr = 0;
         } else {    } else {
                 gcdn(DN(n1),DN(n2),&g); divsn(DN(n1),g,&dn1); divsn(DN(n2),g,&dn2);      gcdn(DN(n1),DN(n2),&g); divsn(DN(n1),g,&dn1); divsn(DN(n2),g,&dn2);
                 kmuln(NM(n1),dn2,&nm1); kmuln(NM(n2),dn1,&nm2);      kmuln(NM(n1),dn2,&nm1); kmuln(NM(n2),dn1,&nm2);
                 if ( SGN(n1) == SGN(n2) ) {      if ( SGN(n1) == SGN(n2) ) {
                         sgn = SGN(n1); addn(nm1,nm2,&nm3);        sgn = SGN(n1); addn(nm1,nm2,&nm3);
                 } else      } else
                         sgn = SGN(n1)*subn(nm1,nm2,&nm3);        sgn = SGN(n1)*subn(nm1,nm2,&nm3);
                 if ( sgn ) {      if ( sgn ) {
                         gcdn(nm3,g,&g1); divsn(nm3,g1,&nm); divsn(g,g1,&g0);        gcdn(nm3,g,&g1); divsn(nm3,g1,&nm); divsn(g,g1,&g0);
                         kmuln(dn1,dn2,&dn0); kmuln(g0,dn0,&dn);        kmuln(dn1,dn2,&dn0); kmuln(g0,dn0,&dn);
                         if ( UNIN(dn) )        if ( UNIN(dn) )
                                 NTOQ(nm,sgn,*nr);          NTOQ(nm,sgn,*nr);
                         else        else
                                 NDTOQ(nm,dn,sgn,*nr);          NDTOQ(nm,dn,sgn,*nr);
                 } else      } else
                         *nr = 0;        *nr = 0;
         }    }
 }  }
   
 void subq(Q n1,Q n2,Q *nr)  void subq(Q n1,Q n2,Q *nr)
 {  {
         Q m;    Q m;
   
         if ( !n1 )    if ( !n1 )
                 if ( !n2 )      if ( !n2 )
                         *nr = 0;        *nr = 0;
                 else {      else {
                         DUPQ(n2,*nr); SGN(*nr) = -SGN(n2);        DUPQ(n2,*nr); SGN(*nr) = -SGN(n2);
                 }      }
         else if ( !n2 )    else if ( !n2 )
                 *nr = n1;      *nr = n1;
         else if ( n1 == n2 )    else if ( n1 == n2 )
                 *nr = 0;      *nr = 0;
         else {    else {
                 DUPQ(n2,m); SGN(m) = -SGN(n2); addq(n1,m,nr);      DUPQ(n2,m); SGN(m) = -SGN(n2); addq(n1,m,nr);
         }    }
 }  }
   
 void mulq(Q n1,Q n2,Q *nr)  void mulq(Q n1,Q n2,Q *nr)
 {  {
         N nm,nm1,nm2,dn,dn1,dn2,g;    N nm,nm1,nm2,dn,dn1,dn2,g;
         int sgn;    int sgn;
         unsigned u,l;    unsigned u,l;
   
         if ( !n1 || !n2 ) *nr = 0;    if ( !n1 || !n2 ) *nr = 0;
         else if ( INT(n1) )    else if ( INT(n1) )
                 if ( INT(n2) ) {      if ( INT(n2) ) {
                         nm1 = NM(n1); nm2 = NM(n2);        nm1 = NM(n1); nm2 = NM(n2);
                         if ( PL(nm1) == 1 && PL(nm2) == 1 ) {        if ( PL(nm1) == 1 && PL(nm2) == 1 ) {
                                 DM(BD(NM(n1))[0],BD(NM(n2))[0],u,l)          DM(BD(NM(n1))[0],BD(NM(n2))[0],u,l)
                                 if ( u ) {          if ( u ) {
                                         nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = l; BD(nm)[1] = u;            nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = l; BD(nm)[1] = u;
                                 } else {          } else {
                                         nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = l;            nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = l;
                                 }          }
                         } else        } else
                                 kmuln(nm1,nm2,&nm);          kmuln(nm1,nm2,&nm);
                         sgn = (SGN(n1)==SGN(n2)?1:-1); NTOQ(nm,sgn,*nr);        sgn = (SGN(n1)==SGN(n2)?1:-1); NTOQ(nm,sgn,*nr);
                 } else {      } else {
                         gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn);        gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn);
                         kmuln(nm1,NM(n2),&nm); sgn = SGN(n1)*SGN(n2);        kmuln(nm1,NM(n2),&nm); sgn = SGN(n1)*SGN(n2);
                         if ( UNIN(dn) )        if ( UNIN(dn) )
                                 NTOQ(nm,sgn,*nr);          NTOQ(nm,sgn,*nr);
                         else        else
                                 NDTOQ(nm,dn,sgn,*nr);          NDTOQ(nm,dn,sgn,*nr);
                 }      }
         else if ( INT(n2) ) {    else if ( INT(n2) ) {
                 gcdn(NM(n2),DN(n1),&g); divsn(NM(n2),g,&nm2); divsn(DN(n1),g,&dn);      gcdn(NM(n2),DN(n1),&g); divsn(NM(n2),g,&nm2); divsn(DN(n1),g,&dn);
                 kmuln(nm2,NM(n1),&nm); sgn = SGN(n1)*SGN(n2);      kmuln(nm2,NM(n1),&nm); sgn = SGN(n1)*SGN(n2);
                 if ( UNIN(dn) )      if ( UNIN(dn) )
                         NTOQ(nm,sgn,*nr);        NTOQ(nm,sgn,*nr);
                 else      else
                         NDTOQ(nm,dn,sgn,*nr);        NDTOQ(nm,dn,sgn,*nr);
         } else {    } else {
                 gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn2);      gcdn(NM(n1),DN(n2),&g); divsn(NM(n1),g,&nm1); divsn(DN(n2),g,&dn2);
                 gcdn(DN(n1),NM(n2),&g); divsn(DN(n1),g,&dn1); divsn(NM(n2),g,&nm2);      gcdn(DN(n1),NM(n2),&g); divsn(DN(n1),g,&dn1); divsn(NM(n2),g,&nm2);
                 kmuln(nm1,nm2,&nm); kmuln(dn1,dn2,&dn); sgn = SGN(n1) * SGN(n2);      kmuln(nm1,nm2,&nm); kmuln(dn1,dn2,&dn); sgn = SGN(n1) * SGN(n2);
                 if ( UNIN(dn) )      if ( UNIN(dn) )
                         NTOQ(nm,sgn,*nr);        NTOQ(nm,sgn,*nr);
                 else      else
                         NDTOQ(nm,dn,sgn,*nr);        NDTOQ(nm,dn,sgn,*nr);
         }    }
 }  }
   
 void divq(Q n1,Q n2,Q *nq)  void divq(Q n1,Q n2,Q *nq)
 {  {
         Q m;    Q m;
   
         if ( !n2 ) {    if ( !n2 ) {
                 error("division by 0");      error("division by 0");
                 *nq = 0;      *nq = 0;
                 return;      return;
         } else if ( !n1 )    } else if ( !n1 )
                 *nq = 0;      *nq = 0;
         else if ( n1 == n2 )    else if ( n1 == n2 )
                 *nq = ONE;      *nq = ONE;
         else {    else {
                 invq(n2,&m); mulq(n1,m,nq);      invq(n2,&m); mulq(n1,m,nq);
         }    }
 }  }
   
 void divsq(Q n1,Q n2,Q *nq)  void divsq(Q n1,Q n2,Q *nq)
 {  {
         int sgn;    int sgn;
         N tn;    N tn;
   
         if ( !n2 ) {    if ( !n2 ) {
                 error("division by 0");      error("division by 0");
                 *nq = 0;      *nq = 0;
                 return;      return;
         } else if ( !n1 )    } else if ( !n1 )
                 *nq = 0;      *nq = 0;
         else if ( n1 == n2 )    else if ( n1 == n2 )
                 *nq = ONE;      *nq = ONE;
         else {    else {
                 divsn(NM(n1),NM(n2),&tn);      divsn(NM(n1),NM(n2),&tn);
                 sgn = SGN(n1)*SGN(n2);      sgn = SGN(n1)*SGN(n2);
                 NTOQ(tn,sgn,*nq);      NTOQ(tn,sgn,*nq);
         }    }
 }  }
   
 void invq(Q n,Q *nr)  void invq(Q n,Q *nr)
 {  {
         N nm,dn;    N nm,dn;
   
         if ( INT(n) )    if ( INT(n) )
                 if ( UNIN(NM(n)) )      if ( UNIN(NM(n)) )
                         if ( SGN(n) > 0 )        if ( SGN(n) > 0 )
                                 *nr = ONE;          *nr = ONE;
                         else {        else {
                                 nm = ONEN; NTOQ(nm,SGN(n),*nr);          nm = ONEN; NTOQ(nm,SGN(n),*nr);
                         }        }
                 else {      else {
                         nm = ONEN; dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);        nm = ONEN; dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);
                 }      }
         else if ( UNIN(NM(n)) ) {    else if ( UNIN(NM(n)) ) {
                 nm = DN(n); NTOQ(nm,SGN(n),*nr);      nm = DN(n); NTOQ(nm,SGN(n),*nr);
         } else {    } else {
                 nm = DN(n); dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);      nm = DN(n); dn = NM(n); NDTOQ(nm,dn,SGN(n),*nr);
         }    }
 }  }
   
 void chsgnq(Q n,Q *nr)  void chsgnq(Q n,Q *nr)
 {  {
         Q t;    Q t;
   
         if ( !n )    if ( !n )
                 *nr = 0;      *nr = 0;
         else {    else {
                 DUPQ(n,t); SGN(t) = -SGN(t); *nr = t;      DUPQ(n,t); SGN(t) = -SGN(t); *nr = t;
         }    }
 }  }
   
 void pwrq(Q n1,Q n,Q *nr)  void pwrq(Q n1,Q n,Q *nr)
 {  {
         N nm,dn;    N nm,dn;
         int sgn;    int sgn;
   
         if ( !n || UNIQ(n1) )    if ( !n || UNIQ(n1) )
                 *nr = ONE;      *nr = ONE;
         else if ( !n1 )    else if ( !n1 )
                 *nr = 0;      *nr = 0;
         else if ( !INT(n) ) {    else if ( !INT(n) ) {
                 error("can't calculate fractional power."); *nr = 0;      error("can't calculate fractional power."); *nr = 0;
         } else if ( MUNIQ(n1) )    } else if ( MUNIQ(n1) )
                 *nr = BD(NM(n))[0]%2 ? n1 : ONE;      *nr = BD(NM(n))[0]%2 ? n1 : ONE;
         else if ( PL(NM(n)) > 1 ) {    else if ( PL(NM(n)) > 1 ) {
                 error("exponent too big."); *nr = 0;      error("exponent too big."); *nr = 0;
         } else {    } else {
                 sgn = ((SGN(n1)>0)||EVENN(NM(n))?1:-1);      sgn = ((SGN(n1)>0)||EVENN(NM(n))?1:-1);
                 pwrn(NM(n1),BD(NM(n))[0],&nm);      pwrn(NM(n1),BD(NM(n))[0],&nm);
                 if ( INT(n1) ) {      if ( INT(n1) ) {
                         if ( UNIN(nm) )        if ( UNIN(nm) )
                                 if ( sgn == 1 )          if ( sgn == 1 )
                                         *nr = ONE;            *nr = ONE;
                                 else          else
                                         NTOQ(ONEN,sgn,*nr);            NTOQ(ONEN,sgn,*nr);
                         else if ( SGN(n) > 0 )        else if ( SGN(n) > 0 )
                                 NTOQ(nm,sgn,*nr);          NTOQ(nm,sgn,*nr);
                         else        else
                                 NDTOQ(ONEN,nm,sgn,*nr);          NDTOQ(ONEN,nm,sgn,*nr);
                 } else {      } else {
                         pwrn(DN(n1),BD(NM(n))[0],&dn);        pwrn(DN(n1),BD(NM(n))[0],&dn);
                         if ( SGN(n) > 0 )        if ( SGN(n) > 0 )
                                 NDTOQ(nm,dn,sgn,*nr);          NDTOQ(nm,dn,sgn,*nr);
                         else if ( UNIN(nm) )        else if ( UNIN(nm) )
                                 NTOQ(dn,sgn,*nr);          NTOQ(dn,sgn,*nr);
                         else        else
                                 NDTOQ(dn,nm,sgn,*nr);          NDTOQ(dn,nm,sgn,*nr);
                 }      }
         }    }
 }  }
   
 int cmpq(Q q1,Q q2)  int cmpq(Q q1,Q q2)
 {  {
         int sgn;    int sgn;
         N t,s;    N t,s;
   
         if ( !q1 )    if ( !q1 )
                 if ( !q2 )      if ( !q2 )
                         return ( 0 );        return ( 0 );
                 else      else
                         return ( -SGN(q2) );        return ( -SGN(q2) );
         else if ( !q2 )    else if ( !q2 )
                 return ( SGN(q1) );      return ( SGN(q1) );
         else if ( SGN(q1) != SGN(q2) )    else if ( SGN(q1) != SGN(q2) )
                         return ( SGN(q1) );        return ( SGN(q1) );
         else if ( INT(q1) )    else if ( INT(q1) )
                         if ( INT(q2) ) {        if ( INT(q2) ) {
                                 sgn = cmpn(NM(q1),NM(q2));          sgn = cmpn(NM(q1),NM(q2));
                                 if ( !sgn )          if ( !sgn )
                                         return ( 0 );            return ( 0 );
                                 else          else
                                         return ( SGN(q1)==sgn?1:-1 );            return ( SGN(q1)==sgn?1:-1 );
                         } else {        } else {
                                 kmuln(NM(q1),DN(q2),&t); sgn = cmpn(t,NM(q2));          kmuln(NM(q1),DN(q2),&t); sgn = cmpn(t,NM(q2));
                                 return ( SGN(q1) * sgn );          return ( SGN(q1) * sgn );
                         }        }
         else if ( INT(q2) ) {    else if ( INT(q2) ) {
                 kmuln(NM(q2),DN(q1),&t); sgn = cmpn(NM(q1),t);      kmuln(NM(q2),DN(q1),&t); sgn = cmpn(NM(q1),t);
                 return ( SGN(q1) * sgn );      return ( SGN(q1) * sgn );
         } else {    } else {
                 kmuln(NM(q1),DN(q2),&t); kmuln(NM(q2),DN(q1),&s); sgn = cmpn(t,s);      kmuln(NM(q1),DN(q2),&t); kmuln(NM(q2),DN(q1),&s); sgn = cmpn(t,s);
                 return ( SGN(q1) * sgn );      return ( SGN(q1) * sgn );
         }    }
 }  }
   
 void remq(Q n,Q m,Q *nr)  void remq(Q n,Q m,Q *nr)
 {  {
         N q,r,s;    N q,r,s;
   
         if ( !n ) {    if ( !n ) {
                 *nr = 0;      *nr = 0;
                 return;      return;
         }    }
         divn(NM(n),NM(m),&q,&r);    divn(NM(n),NM(m),&q,&r);
         if ( !r )    if ( !r )
                 *nr = 0;      *nr = 0;
         else if ( SGN(n) > 0 )    else if ( SGN(n) > 0 )
                 NTOQ(r,1,*nr);      NTOQ(r,1,*nr);
         else {    else {
                 subn(NM(m),r,&s); NTOQ(s,1,*nr);      subn(NM(m),r,&s); NTOQ(s,1,*nr);
         }    }
 }  }
   
 /* t = [nC0 nC1 ... nCn] */  /* t = [nC0 nC1 ... nCn] */
   
 void mkbc(int n,Q *t)  void mkbc(int n,Q *t)
 {  {
         int i;    int i;
         N c,d;    N c,d;
   
         for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {    for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {
                 c = NALLOC(1); PL(c) = 1; BD(c)[0] = n-i+1;      c = NALLOC(1); PL(c) = 1; BD(c)[0] = n-i+1;
                 kmuln(NM(t[i-1]),c,&d); divin(d,i,&c); NTOQ(c,1,t[i]);      kmuln(NM(t[i-1]),c,&d); divin(d,i,&c); NTOQ(c,1,t[i]);
         }    }
         for ( ; i <= n; i++ )    for ( ; i <= n; i++ )
                 t[i] = t[n-i];      t[i] = t[n-i];
 }  }
   
 /*  /*
Line 384  void mkbc(int n,Q *t)
Line 384  void mkbc(int n,Q *t)
   
 void mkwc(int k,int l,Q *t)  void mkwc(int k,int l,Q *t)
 {  {
         int i,n,up,low;    int i,n,up,low;
         N nm,d,c;    N nm,d,c;
   
         n = MIN(k,l);    n = MIN(k,l);
         for ( t[0] = ONE, i = 1; i <= n; i++ ) {    for ( t[0] = ONE, i = 1; i <= n; i++ ) {
                 DM(k-i+1,l-i+1,up,low);      DM(k-i+1,l-i+1,up,low);
                 if ( up ) {      if ( up ) {
                         nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = low; BD(nm)[1] = up;        nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = low; BD(nm)[1] = up;
                 } else {      } else {
                         nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = low;        nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = low;
                 }      }
                 kmuln(NM(t[i-1]),nm,&d); divin(d,i,&c); NTOQ(c,1,t[i]);      kmuln(NM(t[i-1]),nm,&d); divin(d,i,&c); NTOQ(c,1,t[i]);
         }    }
 }  }
   
 /* mod m table */  /* mod m table */
Line 404  void mkwc(int k,int l,Q *t)
Line 404  void mkwc(int k,int l,Q *t)
   
 void mkwcm(int k,int l,int m,int *t)  void mkwcm(int k,int l,int m,int *t)
 {  {
         int i,n;    int i,n;
         Q *s;    Q *s;
   
         n = MIN(k,l);    n = MIN(k,l);
         s = (Q *)ALLOCA((n+1)*sizeof(Q));    s = (Q *)ALLOCA((n+1)*sizeof(Q));
         mkwc(k,l,s);    mkwc(k,l,s);
         for ( i = 0; i <= n; i++ ) {    for ( i = 0; i <= n; i++ ) {
                 t[i] = rem(NM(s[i]),m);      t[i] = rem(NM(s[i]),m);
         }    }
 }  }
   
 #if 0  #if 0
 void factorial(int n,Q *r)  void factorial(int n,Q *r)
 {  {
         Q t,iq,s;    Q t,iq,s;
         unsigned int i,m,m0;    unsigned int i,m,m0;
         unsigned int c;    unsigned int c;
   
         if ( !n )    if ( !n )
                 *r = ONE;      *r = ONE;
         else if ( n < 0 )    else if ( n < 0 )
                 *r = 0;      *r = 0;
         else {    else {
                 for ( i = 1, t = ONE; (int)i <= n; ) {      for ( i = 1, t = ONE; (int)i <= n; ) {
                         for ( m0 = m = 1, c = 0; ((int)i <= n) && !c; i++ ) {        for ( m0 = m = 1, c = 0; ((int)i <= n) && !c; i++ ) {
                                 m0 = m; DM(m0,i,c,m)          m0 = m; DM(m0,i,c,m)
                         }        }
                         if ( c ) {        if ( c ) {
                                 m = m0; i--;          m = m0; i--;
                         }        }
                         UTOQ(m,iq); mulq(t,iq,&s); t = s;        UTOQ(m,iq); mulq(t,iq,&s); t = s;
                 }      }
                 *r = t;      *r = t;
         }    }
 }  }
 #endif  #endif
   
Line 445  void partial_factorial(int s,int e,N *r);
Line 445  void partial_factorial(int s,int e,N *r);
   
 void factorial(int n,Q *r)  void factorial(int n,Q *r)
 {  {
         N nm;    N nm;
   
         if ( !n )    if ( !n )
                 *r = ONE;      *r = ONE;
         else if ( n < 0 )    else if ( n < 0 )
                 *r = 0;      *r = 0;
         else {    else {
                 partial_factorial(1,n,&nm);      partial_factorial(1,n,&nm);
                 NTOQ(nm,1,*r);      NTOQ(nm,1,*r);
         }    }
 }  }
   
 /* s*(s+1)*...*e */  /* s*(s+1)*...*e */
 void partial_factorial(int s,int e,N *r)  void partial_factorial(int s,int e,N *r)
 {  {
         unsigned int len,i,m,m0,c;    unsigned int len,i,m,m0,c;
         unsigned int *p,*p1,*p2;    unsigned int *p,*p1,*p2;
         N nm,r1,r2;    N nm,r1,r2;
   
         /* XXX */    /* XXX */
         if ( e-s > 2000 ) {    if ( e-s > 2000 ) {
                 m = (e+s)/2;      m = (e+s)/2;
                 partial_factorial(s,m,&r1);      partial_factorial(s,m,&r1);
                 partial_factorial(m+1,e,&r2);      partial_factorial(m+1,e,&r2);
                 kmuln(r1,r2,r);      kmuln(r1,r2,r);
                 return;      return;
         }    }
         /* find i s.t. 2^(i-1) < m <= 2^i */    /* find i s.t. 2^(i-1) < m <= 2^i */
         for ( i = 0, m = 1; m < e; m <<=1, i++ );    for ( i = 0, m = 1; m < e; m <<=1, i++ );
          /* XXX estimate of word length of the result */     /* XXX estimate of word length of the result */
         len = (i*(e-s+1)+1+31)/32;    len = (i*(e-s+1)+1+31)/32;
         p = ALLOCA((len+1)*sizeof(int));    p = ALLOCA((len+1)*sizeof(int));
         p1 = ALLOCA((len+1)*sizeof(int));    p1 = ALLOCA((len+1)*sizeof(int));
         p[0] = s;    p[0] = s;
         len = 1;    len = 1;
         for ( i = s+1; (int)i <= e; ) {    for ( i = s+1; (int)i <= e; ) {
                 for ( m0 = m = 1, c = 0; ((int)i <= e) && !c; i++ ) {      for ( m0 = m = 1, c = 0; ((int)i <= e) && !c; i++ ) {
                         m0 = m; DM(m0,i,c,m)        m0 = m; DM(m0,i,c,m)
                 }      }
                 if ( c ) {      if ( c ) {
                         m = m0; i--;        m = m0; i--;
                 }      }
                 bzero(p1,(len+1)*sizeof(int));      bzero(p1,(len+1)*sizeof(int));
                 muln_1(p,len,m,p1);      muln_1(p,len,m,p1);
                 if ( p1[len] )      if ( p1[len] )
                         len++;        len++;
                 p2 = p; p = p1; p1 = p2;      p2 = p; p = p1; p1 = p2;
         }    }
         *r = nm = NALLOC(len);    *r = nm = NALLOC(len);
         bcopy(p,BD(nm),sizeof(int)*len);    bcopy(p,BD(nm),sizeof(int)*len);
         PL(nm) = len;    PL(nm) = len;
 }  }
   
 void invl(Q a,Q mod,Q *ar)  void invl(Q a,Q mod,Q *ar)
 {  {
         Q f1,f2,a1,a2,q,m,s;    Q f1,f2,a1,a2,q,m,s;
         N qn,rn;    N qn,rn;
   
         a1 = ONE; a2 = 0;    a1 = ONE; a2 = 0;
         DUPQ(a,f1); SGN(f1) = 1; f2 = mod;    DUPQ(a,f1); SGN(f1) = 1; f2 = mod;
   
         while ( 1 ) {    while ( 1 ) {
                 divn(NM(f1),NM(f2),&qn,&rn);      divn(NM(f1),NM(f2),&qn,&rn);
                 if ( !qn )      if ( !qn )
                         q = 0;        q = 0;
                 else      else
                         NTOQ(qn,1,q);        NTOQ(qn,1,q);
   
                 f1 = f2;      f1 = f2;
                 if ( !rn )      if ( !rn )
                         break;        break;
                 else      else
                         NTOQ(rn,1,f2);        NTOQ(rn,1,f2);
   
                 mulq(a2,q,&m); subq(a1,m,&s); a1 = a2;      mulq(a2,q,&m); subq(a1,m,&s); a1 = a2;
                 if ( !s )      if ( !s )
                         a2 = 0;        a2 = 0;
                 else      else
                         remq(s,mod,&a2);        remq(s,mod,&a2);
         }    }
         if ( SGN(a) < 0 )    if ( SGN(a) < 0 )
                 chsgnq(a2,&m);      chsgnq(a2,&m);
         else    else
                 m = a2;      m = a2;
   
         if ( SGN(m) >= 0 )    if ( SGN(m) >= 0 )
                 *ar = m;      *ar = m;
         else    else
                 addq(m,mod,ar);      addq(m,mod,ar);
 }  }
   
 int kara_mag=100;  int kara_mag=100;
   
 void kmuln(N n1,N n2,N *nr)  void kmuln(N n1,N n2,N *nr)
 {  {
         N n,t,s,m,carry;    N n,t,s,m,carry;
         int d,d1,d2,len,i,l;    int d,d1,d2,len,i,l;
         int *r,*r0;    int *r,*r0;
   
         if ( !n1 || !n2 ) {    if ( !n1 || !n2 ) {
                 *nr = 0; return;      *nr = 0; return;
         }    }
         d1 = PL(n1); d2 = PL(n2);    d1 = PL(n1); d2 = PL(n2);
         if ( (d1 < kara_mag) || (d2 < kara_mag) ) {    if ( (d1 < kara_mag) || (d2 < kara_mag) ) {
                 muln(n1,n2,nr); return;      muln(n1,n2,nr); return;
         }    }
         if ( d1 < d2 ) {    if ( d1 < d2 ) {
                 n = n1; n1 = n2; n2 = n;      n = n1; n1 = n2; n2 = n;
                 d = d1; d1 = d2; d2 = d;      d = d1; d1 = d2; d2 = d;
         }    }
         if ( d2 > (d1+1)/2 ) {    if ( d2 > (d1+1)/2 ) {
                 kmulnmain(n1,n2,nr); return;      kmulnmain(n1,n2,nr); return;
         }    }
         d = (d1/d2)+((d1%d2)!=0);    d = (d1/d2)+((d1%d2)!=0);
         len = (d+1)*d2;    len = (d+1)*d2;
         r0 = (int *)ALLOCA(len*sizeof(int));    r0 = (int *)ALLOCA(len*sizeof(int));
         bzero((char *)r0,len*sizeof(int));    bzero((char *)r0,len*sizeof(int));
         for ( carry = 0, i = 0, r = r0; i < d; i++, r += d2 ) {    for ( carry = 0, i = 0, r = r0; i < d; i++, r += d2 ) {
                 extractn(n1,i*d2,d2,&m);      extractn(n1,i*d2,d2,&m);
                 if ( m ) {      if ( m ) {
                         kmuln(m,n2,&t);        kmuln(m,n2,&t);
                         addn(t,carry,&s);        addn(t,carry,&s);
                         copyn(s,d2,r);        copyn(s,d2,r);
                         extractn(s,d2,d2,&carry);        extractn(s,d2,d2,&carry);
                 } else {      } else {
                         copyn(carry,d2,r);        copyn(carry,d2,r);
                         carry = 0;        carry = 0;
                 }      }
         }    }
         copyn(carry,d2,r);    copyn(carry,d2,r);
         for ( l = len - 1; !r0[l]; l-- );    for ( l = len - 1; !r0[l]; l-- );
         l++;    l++;
         *nr = t = NALLOC(l); PL(t) = l;    *nr = t = NALLOC(l); PL(t) = l;
         bcopy((char *)r0,(char *)BD(t),l*sizeof(int));    bcopy((char *)r0,(char *)BD(t),l*sizeof(int));
 }  }
   
 void extractn(N n,int index,int len,N *nr)  void extractn(N n,int index,int len,N *nr)
 {  {
         unsigned int *m;    unsigned int *m;
         int l;    int l;
         N t;    N t;
   
         if ( !n ) {    if ( !n ) {
                 *nr = 0; return;      *nr = 0; return;
         }    }
         m = BD(n)+index;    m = BD(n)+index;
         if ( (l = PL(n)-index) >= len ) {    if ( (l = PL(n)-index) >= len ) {
                 for ( l = len - 1; (l >= 0) && !m[l]; l-- );      for ( l = len - 1; (l >= 0) && !m[l]; l-- );
                 l++;      l++;
         }    }
         if ( l <= 0 ) {    if ( l <= 0 ) {
                 *nr = 0; return;      *nr = 0; return;
         } else {    } else {
                 *nr = t = NALLOC(l); PL(t) = l;      *nr = t = NALLOC(l); PL(t) = l;
                 bcopy((char *)m,(char *)BD(t),l*sizeof(int));      bcopy((char *)m,(char *)BD(t),l*sizeof(int));
         }    }
 }  }
   
 void copyn(N n,int len,int *p)  void copyn(N n,int len,int *p)
 {  {
         if ( n )    if ( n )
                 bcopy((char *)BD(n),(char *)p,MIN(PL(n),len)*sizeof(int));      bcopy((char *)BD(n),(char *)p,MIN(PL(n),len)*sizeof(int));
 }  }
   
 void dupn(N n,N p)  void dupn(N n,N p)
 {  {
         if ( n )    if ( n )
                 bcopy((char *)n,(char *)p,(1+PL(n))*sizeof(int));      bcopy((char *)n,(char *)p,(1+PL(n))*sizeof(int));
 }  }
   
 void kmulnmain(N n1,N n2,N *nr)  void kmulnmain(N n1,N n2,N *nr)
 {  {
         int d1,d2,h,sgn,sgn1,sgn2,len;    int d1,d2,h,sgn,sgn1,sgn2,len;
         N n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;    N n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
   
         d1 = PL(n1); d2 = PL(n2); h = (d1+1)/2;    d1 = PL(n1); d2 = PL(n2); h = (d1+1)/2;
         extractn(n1,0,h,&n1lo); extractn(n1,h,d1-h,&n1hi);    extractn(n1,0,h,&n1lo); extractn(n1,h,d1-h,&n1hi);
         extractn(n2,0,h,&n2lo); extractn(n2,h,d2-h,&n2hi);    extractn(n2,0,h,&n2lo); extractn(n2,h,d2-h,&n2hi);
         kmuln(n1hi,n2hi,&hi); kmuln(n1lo,n2lo,&lo);    kmuln(n1hi,n2hi,&hi); kmuln(n1lo,n2lo,&lo);
         len = PL(hi)+2*h; t1 = NALLOC(len); PL(t1) = len;    len = PL(hi)+2*h; t1 = NALLOC(len); PL(t1) = len;
         bzero((char *)BD(t1),len*sizeof(int));    bzero((char *)BD(t1),len*sizeof(int));
         if ( lo )    if ( lo )
                 bcopy((char *)BD(lo),(char *)BD(t1),PL(lo)*sizeof(int));      bcopy((char *)BD(lo),(char *)BD(t1),PL(lo)*sizeof(int));
         if ( hi )    if ( hi )
                 bcopy((char *)BD(hi),(char *)(BD(t1)+2*h),PL(hi)*sizeof(int));      bcopy((char *)BD(hi),(char *)(BD(t1)+2*h),PL(hi)*sizeof(int));
   
         addn(hi,lo,&mid1);    addn(hi,lo,&mid1);
         sgn1 = subn(n1hi,n1lo,&s1); sgn2 = subn(n2lo,n2hi,&s2);    sgn1 = subn(n1hi,n1lo,&s1); sgn2 = subn(n2lo,n2hi,&s2);
         kmuln(s1,s2,&mid2); sgn = sgn1*sgn2;    kmuln(s1,s2,&mid2); sgn = sgn1*sgn2;
         if ( sgn > 0 )    if ( sgn > 0 )
                 addn(mid1,mid2,&mid);      addn(mid1,mid2,&mid);
         else    else
                 subn(mid1,mid2,&mid);      subn(mid1,mid2,&mid);
         if ( mid ) {    if ( mid ) {
                 len = PL(mid)+h; t2 = NALLOC(len); PL(t2) = len;      len = PL(mid)+h; t2 = NALLOC(len); PL(t2) = len;
                 bzero((char *)BD(t2),len*sizeof(int));      bzero((char *)BD(t2),len*sizeof(int));
                 bcopy((char *)BD(mid),(char *)(BD(t2)+h),PL(mid)*sizeof(int));      bcopy((char *)BD(mid),(char *)(BD(t2)+h),PL(mid)*sizeof(int));
                 addn(t1,t2,nr);      addn(t1,t2,nr);
         } else    } else
                 *nr = t1;      *nr = t1;
 }  }
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.10
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>