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

Diff for /OpenXM_contrib2/asir2000/engine/up.c between version 1.1.1.1 and 1.6

version 1.1.1.1, 1999/12/03 07:39:08 version 1.6, 2018/03/29 01:32:52
Line 1 
Line 1 
 /* $OpenXM: OpenXM/src/asir99/engine/up.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */  /*
    * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
    * All rights reserved.
    *
    * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
    * non-exclusive and royalty-free license to use, copy, modify and
    * redistribute, solely for non-commercial and non-profit purposes, the
    * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
    * conditions of this Agreement. For the avoidance of doubt, you acquire
    * only a limited right to use the SOFTWARE hereunder, and FLL or any
    * third party developer retains all rights, including but not limited to
    * copyrights, in and to the SOFTWARE.
    *
    * (1) FLL does not grant you a license in any way for commercial
    * purposes. You may use the SOFTWARE only for non-commercial and
    * non-profit purposes only, such as academic, research and internal
    * business use.
    * (2) The SOFTWARE is protected by the Copyright Law of Japan and
    * international copyright treaties. If you make copies of the SOFTWARE,
    * with or without modification, as permitted hereunder, you shall affix
    * to all such copies of the SOFTWARE the above copyright notice.
    * (3) An explicit reference to this SOFTWARE and its copyright owner
    * shall be made on your publication or presentation in any form of the
    * results obtained by use of the SOFTWARE.
    * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
    * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
    * for such modification or the source code of the modified part of the
    * SOFTWARE.
    *
    * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
    * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
    * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
    * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
    * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
    * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
    * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
    * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
    * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
    * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
    * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
    * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
    * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
    * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
    * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
    * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
    * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
    *
    * $OpenXM: OpenXM_contrib2/asir2000/engine/up.c,v 1.5 2001/10/09 01:36:13 noro Exp $
   */
 #include "ca.h"  #include "ca.h"
 #include <math.h>  #include <math.h>
   
Line 28  extern int lm_lazy;
Line 76  extern int lm_lazy;
 extern int current_ff;  extern int current_ff;
 extern int GC_dont_gc;  extern int GC_dont_gc;
   
 void monicup(a,b)  void monicup(UP a,UP *b)
 UP a;  
 UP *b;  
 {  {
         UP w;    UP w;
   
         if ( !a )    if ( !a )
                 *b = 0;      *b = 0;
         else {    else {
                 w = W_UPALLOC(0); w->d = 0;      w = W_UPALLOC(0); w->d = 0;
                 divnum(0,(Num)ONE,(Num)a->c[a->d],(Num *)&w->c[0]);      divnum(0,(Num)ONE,(Num)a->c[a->d],(Num *)&w->c[0]);
                 mulup(a,w,b);      mulup(a,w,b);
         }    }
 }  }
   
 void simpup(a,b)  void simpup(UP a,UP *b)
 UP a;  
 UP *b;  
 {  {
         int i,d;    int i,d;
         UP c;    UP c;
   
         if ( !a )    if ( !a )
                 *b = 0;      *b = 0;
         else {    else {
                 d = a->d;      d = a->d;
                 c = UPALLOC(d);      c = UPALLOC(d);
   
                 for ( i = 0; i <= d; i++ )      for ( i = 0; i <= d; i++ )
                         simpnum(a->c[i],&c->c[i]);        simpnum(a->c[i],&c->c[i]);
                 for ( i = d; i >= 0 && !c->c[i]; i-- );      for ( i = d; i >= 0 && !c->c[i]; i-- );
                 if ( i < 0 )      if ( i < 0 )
                         *b = 0;        *b = 0;
                 else {      else {
                         c->d = i;        c->d = i;
                         *b = c;        *b = c;
                 }      }
         }    }
 }  }
   
 void simpnum(a,b)  void simpnum(Num a,Num *b)
 Num a;  
 Num *b;  
 {  {
         LM lm;    LM lm;
         GF2N gf;    GF2N gf;
         GFPN gfpn;    GFPN gfpn;
   
         if ( !a )    if ( !a )
                 *b = 0;      *b = 0;
         else    else
                 switch ( NID(a) ) {      switch ( NID(a) ) {
                         case N_LM:        case N_LM:
                                 simplm((LM)a,&lm); *b = (Num)lm; break;          simplm((LM)a,&lm); *b = (Num)lm; break;
                         case N_GF2N:        case N_GF2N:
                                 simpgf2n((GF2N)a,&gf); *b = (Num)gf; break;          simpgf2n((GF2N)a,&gf); *b = (Num)gf; break;
                         case N_GFPN:        case N_GFPN:
                                 simpgfpn((GFPN)a,&gfpn); *b = (Num)gfpn; break;          simpgfpn((GFPN)a,&gfpn); *b = (Num)gfpn; break;
                         default:        default:
                                 *b = a; break;          *b = a; break;
                 }      }
 }  }
   
 void uremp(p1,p2,rp)  void uremp(P p1,P p2,P *rp)
 P p1,p2;  
 P *rp;  
 {  {
         UP n1,n2,r;    UP n1,n2,r;
   
         if ( !p1 || NUM(p1) )    if ( !p1 || NUM(p1) )
                 *rp = p1;      *rp = p1;
         else {    else {
                 ptoup(p1,&n1); ptoup(p2,&n2);      ptoup(p1,&n1); ptoup(p2,&n2);
                 remup(n1,n2,&r);      remup(n1,n2,&r);
                 uptop(r,rp);      uptop(r,rp);
         }    }
 }  }
   
 void ugcdp(p1,p2,rp)  void ugcdp(P p1,P p2,P *rp)
 P p1,p2;  
 P *rp;  
 {  {
         UP n1,n2,r;    UP n1,n2,r;
   
         ptoup(p1,&n1); ptoup(p2,&n2);    ptoup(p1,&n1); ptoup(p2,&n2);
         gcdup(n1,n2,&r);    gcdup(n1,n2,&r);
         uptop(r,rp);    uptop(r,rp);
 }  }
   
 void reversep(p1,d,rp)  void reversep(P p1,Q d,P *rp)
 P p1;  
 Q d;  
 P *rp;  
 {  {
         UP n1,r;    UP n1,r;
   
         if ( !p1 )    if ( !p1 )
                 *rp = 0;      *rp = 0;
         else {    else {
                 ptoup(p1,&n1);      ptoup(p1,&n1);
                 reverseup(n1,QTOS(d),&r);      reverseup(n1,QTOS(d),&r);
                 uptop(r,rp);      uptop(r,rp);
         }    }
 }  }
   
 void invmodp(p1,d,rp)  void invmodp(P p1,Q d,P *rp)
 P p1;  
 Q d;  
 P *rp;  
 {  {
         UP n1,r;    UP n1,r;
   
         if ( !p1 )    if ( !p1 )
                 error("invmodp : invalid input");      error("invmodp : invalid input");
         else {    else {
                 ptoup(p1,&n1);      ptoup(p1,&n1);
                 invmodup(n1,QTOS(d),&r);      invmodup(n1,QTOS(d),&r);
                 uptop(r,rp);      uptop(r,rp);
         }    }
 }  }
   
 void addup(n1,n2,nr)  void addup(UP n1,UP n2,UP *nr)
 UP n1,n2;  
 UP *nr;  
 {  {
         UP r,t;    UP r,t;
         int i,d1,d2;    int i,d1,d2;
         Num *c,*c1,*c2;    Num *c,*c1,*c2;
   
         if ( !n1 ) {    if ( !n1 ) {
                 *nr = n2;      *nr = n2;
         } else if ( !n2 ) {    } else if ( !n2 ) {
                 *nr = n1;      *nr = n1;
         } else {    } else {
                 if ( n2->d > n1->d ) {      if ( n2->d > n1->d ) {
                         t = n1; n1 = n2; n2 = t;        t = n1; n1 = n2; n2 = t;
                 }      }
                 d1 = n1->d; d2 = n2->d;      d1 = n1->d; d2 = n2->d;
                 *nr = r = UPALLOC(d1);      *nr = r = UPALLOC(d1);
                 c1 = n1->c; c2 = n2->c; c = r->c;      c1 = n1->c; c2 = n2->c; c = r->c;
                 for ( i = 0; i <= d2; i++ )      for ( i = 0; i <= d2; i++ )
                         addnum(0,*c1++,*c2++,c++);        addnum(0,*c1++,*c2++,c++);
                 for ( ; i <= d1; i++ )      for ( ; i <= d1; i++ )
                         *c++ = *c1++;        *c++ = *c1++;
                 for ( i = d1, c = r->c; i >=0 && !c[i]; i-- );      for ( i = d1, c = r->c; i >=0 && !c[i]; i-- );
                 if ( i < 0 )      if ( i < 0 )
                         *nr = 0;        *nr = 0;
                 else      else
                         r->d = i;        r->d = i;
         }    }
 }  }
   
 void subup(n1,n2,nr)  void subup(UP n1,UP n2,UP *nr)
 UP n1,n2;  
 UP *nr;  
 {  {
         UP r;    UP r;
         int i,d1,d2,d;    int i,d1,d2,d;
         Num *c,*c1,*c2;    Num *c,*c1,*c2;
   
         if ( !n1 ) {    if ( !n1 ) {
                 chsgnup(n2,nr);      chsgnup(n2,nr);
         } else if ( !n2 ) {    } else if ( !n2 ) {
                 *nr = n1;      *nr = n1;
         } else {    } else {
                 d1 = n1->d; d2 = n2->d; d = MAX(d1,d2);      d1 = n1->d; d2 = n2->d; d = MAX(d1,d2);
                 *nr = r = UPALLOC(d);      *nr = r = UPALLOC(d);
                 c1 = n1->c; c2 = n2->c; c = r->c;      c1 = n1->c; c2 = n2->c; c = r->c;
                 if ( d1 >= d2 ) {      if ( d1 >= d2 ) {
                         for ( i = 0; i <= d2; i++ )        for ( i = 0; i <= d2; i++ )
                                 subnum(0,*c1++,*c2++,c++);          subnum(0,*c1++,*c2++,c++);
                         for ( ; i <= d1; i++ )        for ( ; i <= d1; i++ )
                                 *c++ = *c1++;          *c++ = *c1++;
                 } else {      } else {
                         for ( i = 0; i <= d1; i++ )        for ( i = 0; i <= d1; i++ )
                                 subnum(0,*c1++,*c2++,c++);          subnum(0,*c1++,*c2++,c++);
                         for ( ; i <= d2; i++ )        for ( ; i <= d2; i++ )
                                 chsgnnum(*c2++,c++);          chsgnnum(*c2++,c++);
                 }      }
                 for ( i = d, c = r->c; i >=0 && !c[i]; i-- );      for ( i = d, c = r->c; i >=0 && !c[i]; i-- );
                 if ( i < 0 )      if ( i < 0 )
                         *nr = 0;        *nr = 0;
                 else      else
                         r->d = i;        r->d = i;
         }    }
 }  }
   
 void chsgnup(n1,nr)  void chsgnup(UP n1,UP *nr)
 UP n1;  
 UP *nr;  
 {  {
         UP r;    UP r;
         int d1,i;    int d1,i;
         Num *c1,*c;    Num *c1,*c;
   
         if ( !n1 ) {    if ( !n1 ) {
                 *nr = 0;      *nr = 0;
         } else {    } else {
                 d1 = n1->d;      d1 = n1->d;
                 *nr = r = UPALLOC(d1); r->d = d1;      *nr = r = UPALLOC(d1); r->d = d1;
                 c1 = n1->c; c = r->c;      c1 = n1->c; c = r->c;
                 for ( i = 0; i <= d1; i++ )      for ( i = 0; i <= d1; i++ )
                         chsgnnum(*c1++,c++);        chsgnnum(*c1++,c++);
         }    }
 }  }
   
 void hybrid_mulup(ff,n1,n2,nr)  void hybrid_mulup(int ff,UP n1,UP n2,UP *nr)
 int ff;  
 UP n1,n2;  
 UP *nr;  
 {  {
         if ( !n1 || !n2 )    if ( !n1 || !n2 )
                 *nr = 0;      *nr = 0;
         else if ( MAX(n1->d,n2->d) < up_fft_mag )    else if ( MAX(n1->d,n2->d) < up_fft_mag )
                 kmulup(n1,n2,nr);      kmulup(n1,n2,nr);
         else    else
                 switch ( ff ) {      switch ( ff ) {
                         case FF_GFP:        case FF_GFP:
                                 fft_mulup_lm(n1,n2,nr); break;          fft_mulup_lm(n1,n2,nr); break;
                         case FF_GF2N:        case FF_GF2N:
                                 kmulup(n1,n2,nr); break;          kmulup(n1,n2,nr); break;
                         default:        default:
                                 fft_mulup(n1,n2,nr); break;          fft_mulup(n1,n2,nr); break;
                 }      }
 }  }
   
 void hybrid_squareup(ff,n1,nr)  void hybrid_squareup(int ff,UP n1,UP *nr)
 int ff;  
 UP n1;  
 UP *nr;  
 {  {
         if ( !n1 )    if ( !n1 )
                 *nr = 0;      *nr = 0;
         else if ( n1->d < up_fft_mag )    else if ( n1->d < up_fft_mag )
                 ksquareup(n1,nr);      ksquareup(n1,nr);
         else    else
                 switch ( ff ) {      switch ( ff ) {
                         case FF_GFP:        case FF_GFP:
                                 fft_squareup_lm(n1,nr); break;          fft_squareup_lm(n1,nr); break;
                         case FF_GF2N:        case FF_GF2N:
                                 ksquareup(n1,nr); break;          ksquareup(n1,nr); break;
                         default:        default:
                                 fft_squareup(n1,nr); break;          fft_squareup(n1,nr); break;
                 }      }
 }  }
   
 void hybrid_tmulup(ff,n1,n2,d,nr)  void hybrid_tmulup(int ff,UP n1,UP n2,int d,UP *nr)
 int ff;  
 UP n1,n2;  
 int d;  
 UP *nr;  
 {  {
         if ( !n1 || !n2 )    if ( !n1 || !n2 )
                 *nr = 0;      *nr = 0;
         else if ( MAX(n1->d,n2->d) < up_fft_mag )    else if ( MAX(n1->d,n2->d) < up_fft_mag )
                 tkmulup(n1,n2,d,nr);      tkmulup(n1,n2,d,nr);
         else    else
                 switch ( ff ) {      switch ( ff ) {
                         case FF_GFP:        case FF_GFP:
                                 trunc_fft_mulup_lm(n1,n2,d,nr); break;          trunc_fft_mulup_lm(n1,n2,d,nr); break;
                         case FF_GF2N:        case FF_GF2N:
                                 tkmulup(n1,n2,d,nr); break;          tkmulup(n1,n2,d,nr); break;
                         default:        default:
                                 trunc_fft_mulup(n1,n2,d,nr); break;          trunc_fft_mulup(n1,n2,d,nr); break;
                 }      }
 }  }
   
 void mulup(n1,n2,nr)  void mulup(UP n1,UP n2,UP *nr)
 UP n1,n2;  
 UP *nr;  
 {  {
         UP r;    UP r;
         Num *pc1,*pc,*c1,*c2,*c;    Num *pc1,*pc,*c1,*c2,*c;
         Num mul,t,s;    Num mul,t,s;
         int i,j,d1,d2;    int i,j,d1,d2;
   
         if ( !n1 || !n2 )    if ( !n1 || !n2 )
                 *nr = 0;      *nr = 0;
         else {    else {
                 d1 = n1->d; d2 = n2->d;      d1 = n1->d; d2 = n2->d;
                 *nr = r = UPALLOC(d1+d2); r->d = d1+d2;      *nr = r = UPALLOC(d1+d2); r->d = d1+d2;
                 c1 = n1->c; c2 = n2->c; c = r->c;      c1 = n1->c; c2 = n2->c; c = r->c;
                 lm_lazy = 1;      lm_lazy = 1;
                 for ( i = 0; i <= d2; i++, c++ )      for ( i = 0; i <= d2; i++, c++ )
                         if ( mul = *c2++ )        if ( mul = *c2++ )
                                 for ( j = 0, pc1 = c1, pc = c; j <= d1; j++ ) {          for ( j = 0, pc1 = c1, pc = c; j <= d1; j++ ) {
                                         mulnum(0,*pc1++,mul,&t); addnum(0,*pc,t,&s); *pc++ = s;            mulnum(0,*pc1++,mul,&t); addnum(0,*pc,t,&s); *pc++ = s;
                         }        }
                 lm_lazy = 0;      lm_lazy = 0;
                 if ( !up_lazy )      if ( !up_lazy )
                         for ( i = 0, c = r->c; i <= r->d; i++ ) {        for ( i = 0, c = r->c; i <= r->d; i++ ) {
                                 simpnum(c[i],&t); c[i] = t;          simpnum(c[i],&t); c[i] = t;
                         }        }
         }    }
 }  }
   
 /* nr = c*n1 */  /* nr = c*n1 */
   
 void mulcup(c,n1,nr)  void mulcup(Num c,UP n1,UP *nr)
 Num c;  
 UP n1;  
 UP *nr;  
 {  {
         int d;    int d;
         UP r;    UP r;
         Num t;    Num t;
         Num *c1,*cr;    Num *c1,*cr;
         int i;    int i;
   
         if ( !c || !n1 )    if ( !c || !n1 )
                 *nr = 0;      *nr = 0;
         else {    else {
                 d = n1->d;      d = n1->d;
                 *nr = r = UPALLOC(d); r->d = d;      *nr = r = UPALLOC(d); r->d = d;
                 c1 = n1->c; cr = r->c;      c1 = n1->c; cr = r->c;
                 lm_lazy = 1;      lm_lazy = 1;
                 for ( i = 0; i <= d; i++ )      for ( i = 0; i <= d; i++ )
                         mulnum(0,c1[i],c,&cr[i]);        mulnum(0,c1[i],c,&cr[i]);
                 lm_lazy = 0;      lm_lazy = 0;
                 if ( !up_lazy )      if ( !up_lazy )
                         for ( i = 0; i <= d; i++ ) {        for ( i = 0; i <= d; i++ ) {
                                 simpnum(cr[i],&t); cr[i] = t;          simpnum(cr[i],&t); cr[i] = t;
                         }        }
         }    }
 }  }
   
 void tmulup(n1,n2,d,nr)  void tmulup(UP n1,UP n2,int d,UP *nr)
 UP n1,n2;  
 int d;  
 UP *nr;  
 {  {
         UP r;    UP r;
         Num *pc1,*pc,*c1,*c2,*c;    Num *pc1,*pc,*c1,*c2,*c;
         Num mul,t,s;    Num mul,t,s;
         int i,j,d1,d2,dr;    int i,j,d1,d2,dr;
   
         if ( !n1 || !n2 )    if ( !n1 || !n2 )
                 *nr = 0;      *nr = 0;
         else {    else {
                 d1 = n1->d; d2 = n2->d;      d1 = n1->d; d2 = n2->d;
                 dr = MAX(d-1,d1+d2);      dr = MAX(d-1,d1+d2);
                 *nr = r = UPALLOC(dr);      *nr = r = UPALLOC(dr);
                 c1 = n1->c; c2 = n2->c; c = r->c;      c1 = n1->c; c2 = n2->c; c = r->c;
                 lm_lazy = 1;      lm_lazy = 1;
                 for ( i = 0; i <= d2; i++, c++ )      for ( i = 0; i <= d2; i++, c++ )
                         if ( mul = *c2++ )        if ( mul = *c2++ )
                                 for ( j = 0, pc1 = c1, pc = c; (j <= d1) && (i+j < d); j++ ) {          for ( j = 0, pc1 = c1, pc = c; (j <= d1) && (i+j < d); j++ ) {
                                         mulnum(0,*pc1++,mul,&t); addnum(0,*pc,t,&s); *pc++ = s;            mulnum(0,*pc1++,mul,&t); addnum(0,*pc,t,&s); *pc++ = s;
                         }        }
                 lm_lazy = 0;      lm_lazy = 0;
                 c = r->c;      c = r->c;
                 if ( !up_lazy )      if ( !up_lazy )
                         for ( i = 0; i < d; i++ ) {        for ( i = 0; i < d; i++ ) {
                                 simpnum(c[i],&t); c[i] = t;          simpnum(c[i],&t); c[i] = t;
                         }        }
                 for ( i = d-1; i >= 0 && !c[i]; i-- );      for ( i = d-1; i >= 0 && !c[i]; i-- );
                 if ( i < 0 )      if ( i < 0 )
                         *nr = 0;        *nr = 0;
                 else      else
                         r->d = i;        r->d = i;
         }    }
 }  }
   
 void squareup(n1,nr)  void squareup(UP n1,UP *nr)
 UP n1;  
 UP *nr;  
 {  {
         UP r;    UP r;
         Num *c1,*c;    Num *c1,*c;
         Num t,s,u;    Num t,s,u;
         int i,j,h,d1,d;    int i,j,h,d1,d;
   
         if ( !n1 )    if ( !n1 )
                 *nr = 0;      *nr = 0;
         else if ( !n1->d ) {    else if ( !n1->d ) {
                 *nr = r = UPALLOC(0); r->d = 0;      *nr = r = UPALLOC(0); r->d = 0;
                 mulnum(0,n1->c[0],n1->c[0],&r->c[0]);      mulnum(0,n1->c[0],n1->c[0],&r->c[0]);
         } else {    } else {
                 d1 = n1->d;      d1 = n1->d;
                 d = 2*d1;      d = 2*d1;
                 *nr = r = UPALLOC(2*d); r->d = d;      *nr = r = UPALLOC(2*d); r->d = d;
                 c1 = n1->c; c = r->c;      c1 = n1->c; c = r->c;
                 lm_lazy = 1;      lm_lazy = 1;
                 for ( i = 0; i <= d1; i++ )      for ( i = 0; i <= d1; i++ )
                         mulnum(0,c1[i],c1[i],&c[2*i]);        mulnum(0,c1[i],c1[i],&c[2*i]);
                 for ( j = 1; j < d; j++ ) {      for ( j = 1; j < d; j++ ) {
                         h = (j-1)/2;        h = (j-1)/2;
                         for ( i = MAX(j-d1,0), s = 0; i <= h; i++ ) {        for ( i = MAX(j-d1,0), s = 0; i <= h; i++ ) {
                                 mulnum(0,c1[i],c1[j-i],&t); addnum(0,t,s,&u); s = u;          mulnum(0,c1[i],c1[j-i],&t); addnum(0,t,s,&u); s = u;
                         }        }
                         addnum(0,s,s,&t); addnum(0,t,c[j],&u); c[j] = u;        addnum(0,s,s,&t); addnum(0,t,c[j],&u); c[j] = u;
                 }      }
                 lm_lazy = 0;      lm_lazy = 0;
                 if ( !up_lazy )      if ( !up_lazy )
                         for ( i = 0, c = r->c; i <= d; i++ ) {        for ( i = 0, c = r->c; i <= d; i++ ) {
                                 simpnum(c[i],&t); c[i] = t;          simpnum(c[i],&t); c[i] = t;
                         }        }
         }    }
 }  }
   
 void remup(n1,n2,nr)  void remup(UP n1,UP n2,UP *nr)
 UP n1,n2;  
 UP *nr;  
 {  {
         UP w,r;    UP w,r;
   
         if ( !n2 )    if ( !n2 )
                 error("remup : division by 0");      error("remup : division by 0");
         else if ( !n1 || !n2->d )    else if ( !n1 || !n2->d )
                 *nr = 0;      *nr = 0;
         else if ( n1->d < n2->d )    else if ( n1->d < n2->d )
                 *nr = n1;      *nr = n1;
         else {    else {
                 w = W_UPALLOC(n1->d); copyup(n1,w);      w = W_UPALLOC(n1->d); copyup(n1,w);
                 remup_destructive(w,n2);      remup_destructive(w,n2);
                 if ( w->d < 0 )      if ( w->d < 0 )
                         *nr = 0;        *nr = 0;
                 else {      else {
                         *nr = r = UPALLOC(w->d); copyup(w,r);        *nr = r = UPALLOC(w->d); copyup(w,r);
                 }      }
         }    }
 }  }
   
 void remup_destructive(n1,n2)  void remup_destructive(UP n1,UP n2)
 UP n1,n2;  
 {  {
         Num *c1,*c2;    Num *c1,*c2;
         int d1,d2,i,j;    int d1,d2,i,j;
         Num m,t,s,mhc;    Num m,t,s,mhc;
   
         c1 = n1->c; c2 = n2->c;    c1 = n1->c; c2 = n2->c;
         d1 = n1->d; d2 = n2->d;    d1 = n1->d; d2 = n2->d;
         chsgnnum(c2[d2],&mhc);    chsgnnum(c2[d2],&mhc);
         for ( i = d1; i >= d2; i-- ) {    for ( i = d1; i >= d2; i-- ) {
                 simpnum(c1[i],&t); c1[i] = t;      simpnum(c1[i],&t); c1[i] = t;
                 if ( !c1[i] )      if ( !c1[i] )
                         continue;        continue;
                 divnum(0,c1[i],mhc,&m);      divnum(0,c1[i],mhc,&m);
                 lm_lazy = 1;      lm_lazy = 1;
                 for ( j = d2; j >= 0; j-- ) {      for ( j = d2; j >= 0; j-- ) {
                         mulnum(0,c2[j],m,&t); addnum(0,c1[i+j-d2],t,&s);        mulnum(0,c2[j],m,&t); addnum(0,c1[i+j-d2],t,&s);
                         c1[i+j-d2] = s;        c1[i+j-d2] = s;
                 }      }
                 lm_lazy = 0;      lm_lazy = 0;
         }    }
         for ( i = 0; i < d2; i++ ) {    for ( i = 0; i < d2; i++ ) {
                 simpnum(c1[i],&t); c1[i] = t;      simpnum(c1[i],&t); c1[i] = t;
         }    }
         for ( i = d2-1; i >= 0 && !c1[i]; i-- );    for ( i = d2-1; i >= 0 && !c1[i]; i-- );
         n1->d = i;    n1->d = i;
 }  }
   
 void qrup(n1,n2,nq,nr)  void qrup(UP n1,UP n2,UP *nq,UP *nr)
 UP n1,n2;  
 UP *nq,*nr;  
 {  {
         UP w,r,q;    UP w,r,q;
         struct oUP t;    struct oUP t;
         Num m;    Num m;
         int d;    int d;
   
         if ( !n2 )    if ( !n2 )
                 error("qrup : division by 0");      error("qrup : division by 0");
         else if ( !n1 ) {    else if ( !n1 ) {
                 *nq = 0; *nr = 0;      *nq = 0; *nr = 0;
         } else if ( !n2->d )    } else if ( !n2->d )
                 if ( !n1->d ) {      if ( !n1->d ) {
                         divnum(0,n1->c[0],n2->c[0],(Num *)nq); *nr = 0;        divnum(0,n1->c[0],n2->c[0],(Num *)nq); *nr = 0;
                 } else {      } else {
                         divnum(0,(Num)ONE,n2->c[0],&m);        divnum(0,(Num)ONE,n2->c[0],&m);
                         t.d = 0; t.c[0] = m;        t.d = 0; t.c[0] = m;
                         mulup(n1,&t,nq); *nr = 0;        mulup(n1,&t,nq); *nr = 0;
                 }      }
         else if ( n1->d < n2->d ) {    else if ( n1->d < n2->d ) {
                 *nq = 0; *nr = n1;      *nq = 0; *nr = n1;
         } else {    } else {
                 w = W_UPALLOC(n1->d); copyup(n1,w);      w = W_UPALLOC(n1->d); copyup(n1,w);
                 qrup_destructive(w,n2);      qrup_destructive(w,n2);
                 d = n1->d-n2->d;      d = n1->d-n2->d;
                 *nq = q = UPALLOC(d); q->d = d;      *nq = q = UPALLOC(d); q->d = d;
                 bcopy(w->c+n2->d,q->c,(d+1)*sizeof(Num));      bcopy(w->c+n2->d,q->c,(d+1)*sizeof(Num));
                 if ( w->d < 0 )      if ( w->d < 0 )
                         *nr = 0;        *nr = 0;
                 else {      else {
                         *nr = r = UPALLOC(w->d); copyup(w,r);        *nr = r = UPALLOC(w->d); copyup(w,r);
                 }      }
         }    }
 }  }
   
 void qrup_destructive(n1,n2)  void qrup_destructive(UP n1,UP n2)
 UP n1,n2;  
 {  {
         Num *c1,*c2;    Num *c1,*c2;
         int d1,d2,i,j;    int d1,d2,i,j;
         Num q,t,s,hc;    Num q,t,s,hc;
   
         c1 = n1->c; c2 = n2->c;    c1 = n1->c; c2 = n2->c;
         d1 = n1->d; d2 = n2->d;    d1 = n1->d; d2 = n2->d;
         hc = c2[d2];    hc = c2[d2];
         for ( i = d1; i >= d2; i-- ) {    for ( i = d1; i >= d2; i-- ) {
                 simpnum(c1[i],&t); c1[i] = t;      simpnum(c1[i],&t); c1[i] = t;
                 if ( c1[i] ) {      if ( c1[i] ) {
                         divnum(0,c1[i],hc,&q);        divnum(0,c1[i],hc,&q);
                         lm_lazy = 1;        lm_lazy = 1;
                         for ( j = d2; j >= 0; j-- ) {        for ( j = d2; j >= 0; j-- ) {
                                 mulnum(0,c2[j],q,&t); subnum(0,c1[i+j-d2],t,&s);          mulnum(0,c2[j],q,&t); subnum(0,c1[i+j-d2],t,&s);
                                 c1[i+j-d2] = s;          c1[i+j-d2] = s;
                         }        }
                         lm_lazy = 0;        lm_lazy = 0;
                 } else      } else
                         q = 0;        q = 0;
                 c1[i] = q;      c1[i] = q;
         }    }
         for ( i = 0; i < d2; i++ ) {    for ( i = 0; i < d2; i++ ) {
                 simpnum(c1[i],&t); c1[i] = t;      simpnum(c1[i],&t); c1[i] = t;
         }    }
         for ( i = d2-1; i >= 0 && !c1[i]; i-- );    for ( i = d2-1; i >= 0 && !c1[i]; i-- );
         n1->d = i;    n1->d = i;
 }  }
   
 void gcdup(n1,n2,nr)  void gcdup(UP n1,UP n2,UP *nr)
 UP n1,n2;  
 UP *nr;  
 {  {
         UP w1,w2,t,r;    UP w1,w2,t,r;
         int d1,d2;    int d1,d2;
   
         if ( !n1 )    if ( !n1 )
                 *nr = n2;      *nr = n2;
         else if ( !n2 )    else if ( !n2 )
                 *nr = n1;      *nr = n1;
         else if ( !n1->d || !n2->d ) {    else if ( !n1->d || !n2->d ) {
                 *nr = r = UPALLOC(0); r->d = 0;      *nr = r = UPALLOC(0); r->d = 0;
                 divnum(0,n1->c[n1->d],n1->c[n1->d],&r->c[0]);      divnum(0,n1->c[n1->d],n1->c[n1->d],&r->c[0]);
         } else {    } else {
                 d1 = n1->d; d2 = n2->d;      d1 = n1->d; d2 = n2->d;
                 if ( d2 > d1 ) {      if ( d2 > d1 ) {
                         w1 = W_UPALLOC(d2); copyup(n2,w1);        w1 = W_UPALLOC(d2); copyup(n2,w1);
                         w2 = W_UPALLOC(d1); copyup(n1,w2);        w2 = W_UPALLOC(d1); copyup(n1,w2);
                 } else {      } else {
                         w1 = W_UPALLOC(d1); copyup(n1,w1);        w1 = W_UPALLOC(d1); copyup(n1,w1);
                         w2 = W_UPALLOC(d2); copyup(n2,w2);        w2 = W_UPALLOC(d2); copyup(n2,w2);
                 }      }
                 d1 = w1->d; d2 = w2->d;      d1 = w1->d; d2 = w2->d;
                 while ( 1 ) {      while ( 1 ) {
                         remup_destructive(w1,w2);        remup_destructive(w1,w2);
                         if ( w1->d < 0 ) {        if ( w1->d < 0 ) {
                                 *nr = r = UPALLOC(w2->d); copyup(w2,r); return;          *nr = r = UPALLOC(w2->d); copyup(w2,r); return;
                         } else if ( !w1->d ) {        } else if ( !w1->d ) {
                                 *nr = r = UPALLOC(0); r->d = 0;          *nr = r = UPALLOC(0); r->d = 0;
                                 divnum(0,n1->c[n1->d],n1->c[n1->d],&r->c[0]); return;          divnum(0,n1->c[n1->d],n1->c[n1->d],&r->c[0]); return;
                         } else {        } else {
                                 t = w1; w1 = w2; w2 = t;          t = w1; w1 = w2; w2 = t;
                         }        }
                 }      }
         }    }
 }  }
   
 /* compute r s.t. a*r = 1 mod m */  /* compute r s.t. a*r = 1 mod m */
   
 void extended_gcdup(a,m,r)  void extended_gcdup(UP a,UP m,UP *r)
 UP a,m;  
 UP *r;  
 {  {
         UP one,g1,g2,a1,a2,a3,b1,b2,b3,inv,quo,rem,t;    UP one,g1,g2,a1,a2,a3,b1,b2,b3,inv,quo,rem,t;
         Num i;    Num i;
   
         one = UPALLOC(0); one->d = 0; one->c[0] = (Num)ONELM;    one = UPALLOC(0); one->d = 0; one->c[0] = (Num)ONELM;
         g1 = m; g2 = a;    g1 = m; g2 = a;
         a1 = one; a2 = 0;    a1 = one; a2 = 0;
         b1 = 0; b2 = one;    b1 = 0; b2 = one;
         /* a2*m+b2*a = g2 always holds */    /* a2*m+b2*a = g2 always holds */
         while ( g2->d ) {    while ( g2->d ) {
                 qrup(g1,g2,&quo,&rem); g1 = g2; g2 = rem;      qrup(g1,g2,&quo,&rem); g1 = g2; g2 = rem;
                 mulup(a2,quo,&t); subup(a1,t,&a3); a1 = a2; a2 = a3;      mulup(a2,quo,&t); subup(a1,t,&a3); a1 = a2; a2 = a3;
                 mulup(b2,quo,&t); subup(b1,t,&b3); b1 = b2; b2 = b3;      mulup(b2,quo,&t); subup(b1,t,&b3); b1 = b2; b2 = b3;
         }    }
         /* now a2*m+b2*a = g2 (const) */    /* now a2*m+b2*a = g2 (const) */
         divnum(0,(Num)ONE,g2->c[0],&i);    divnum(0,(Num)ONE,g2->c[0],&i);
         inv = UPALLOC(0); inv->d = 0; inv->c[0] = i;    inv = UPALLOC(0); inv->d = 0; inv->c[0] = i;
         mulup(b2,inv,r);    mulup(b2,inv,r);
 }  }
   
 void reverseup(n1,d,nr)  void reverseup(UP n1,int d,UP *nr)
 UP n1;  
 int d;  
 UP *nr;  
 {  {
         Num *c1,*c;    Num *c1,*c;
         int i,d1;    int i,d1;
         UP r;    UP r;
   
         if ( !n1 )    if ( !n1 )
                 *nr = 0;      *nr = 0;
         else if ( d < (d1 = n1->d) )    else if ( d < (d1 = n1->d) )
                 error("reverseup : invalid input");      error("reverseup : invalid input");
         else {    else {
                 *nr = r = UPALLOC(d);      *nr = r = UPALLOC(d);
                 c = r->c;      c = r->c;
                 bzero((char *)c,(d+1)*sizeof(Num));      bzero((char *)c,(d+1)*sizeof(Num));
                 c1 = n1->c;      c1 = n1->c;
                 for ( i = 0; i <= d1; i++ )      for ( i = 0; i <= d1; i++ )
                         c[d-i] = c1[i];        c[d-i] = c1[i];
                 for ( i = d; i >= 0 && !c[i]; i-- );      for ( i = d; i >= 0 && !c[i]; i-- );
                 r->d = i;      r->d = i;
                 if ( i < 0 )      if ( i < 0 )
                         *nr = 0;        *nr = 0;
         }    }
 }  }
   
 void invmodup(n1,d,nr)  void invmodup(UP n1,int d,UP *nr)
 UP n1;  
 int d;  
 UP *nr;  
 {  {
         UP r;    UP r;
         Num s,t,u,hinv;    Num s,t,u,hinv;
         int i,j,dmin;    int i,j,dmin;
         Num *w,*c,*cr;    Num *w,*c,*cr;
   
         if ( !n1 || !n1->c[0] )    if ( !n1 || !n1->c[0] )
                 error("invmodup : invalid input");      error("invmodup : invalid input");
         else if ( !n1->d ) {    else if ( !n1->d ) {
                 *nr = r = UPALLOC(0); r->d = 0;      *nr = r = UPALLOC(0); r->d = 0;
                 divnum(0,(Num)ONE,n1->c[0],&r->c[0]);      divnum(0,(Num)ONE,n1->c[0],&r->c[0]);
         } else {    } else {
                 *nr = r = UPALLOC(d);      *nr = r = UPALLOC(d);
                 w = (Num *)ALLOCA((d+1)*sizeof(Q));      w = (Num *)ALLOCA((d+1)*sizeof(Q));
                 bzero((char *)w,(d+1)*sizeof(Q));      bzero((char *)w,(d+1)*sizeof(Q));
                 dmin = MIN(d,n1->d);      dmin = MIN(d,n1->d);
                 divnum(0,(Num)ONE,n1->c[0],&hinv);      divnum(0,(Num)ONE,n1->c[0],&hinv);
                 for ( i = 0, c = n1->c; i <= dmin; i++ )      for ( i = 0, c = n1->c; i <= dmin; i++ )
                         mulnum(0,c[i],hinv,&w[i]);        mulnum(0,c[i],hinv,&w[i]);
                 cr = r->c;      cr = r->c;
                 cr[0] = w[0];      cr[0] = w[0];
                 for ( i = 1; i <= d; i++ ) {      for ( i = 1; i <= d; i++ ) {
                         for ( j = 1, s = w[i]; j < i; j++ ) {        for ( j = 1, s = w[i]; j < i; j++ ) {
                                 mulnum(0,cr[j],w[i-j],&t);          mulnum(0,cr[j],w[i-j],&t);
                                 addnum(0,s,t,&u);          addnum(0,s,t,&u);
                                 s = u;          s = u;
                         }        }
                         chsgnnum(s,&cr[i]);        chsgnnum(s,&cr[i]);
                 }      }
                 for ( i = 0; i <= d; i++ ) {      for ( i = 0; i <= d; i++ ) {
                         mulnum(0,cr[i],hinv,&t); cr[i] = t;        mulnum(0,cr[i],hinv,&t); cr[i] = t;
                 }      }
                 for ( i = d; i >= 0 && !cr[i]; i-- );      for ( i = d; i >= 0 && !cr[i]; i-- );
                 r->d = i;      r->d = i;
         }    }
 }  }
   
 void pwrup(n,e,nr)  void pwrup(UP n,Q e,UP *nr)
 UP n;  
 Q e;  
 UP *nr;  
 {  {
         UP y,z,t;    UP y,z,t;
         N en,qn;    N en,qn;
         int r;    int r;
   
         y = UPALLOC(0); y->d = 0; y->c[0] = (Num)ONE;    y = UPALLOC(0); y->d = 0; y->c[0] = (Num)ONE;
         z = n;    z = n;
         if ( !e )    if ( !e )
                 *nr = y;      *nr = y;
         else if ( UNIQ(e) )    else if ( UNIQ(e) )
                 *nr = n;      *nr = n;
         else {    else {
                 en = NM(e);      en = NM(e);
                 while ( 1 ) {      while ( 1 ) {
                         r = divin(en,2,&qn); en = qn;        r = divin(en,2,&qn); en = qn;
                         if ( r ) {        if ( r ) {
                                 mulup(z,y,&t); y = t;          mulup(z,y,&t); y = t;
                                 if ( !en ) {          if ( !en ) {
                                         *nr = y;            *nr = y;
                                         return;            return;
                                 }          }
                         }        }
                         mulup(z,z,&t); z = t;        mulup(z,z,&t); z = t;
                 }      }
         }    }
 }  }
   
 int compup(n1,n2)  int compup(UP n1,UP n2)
 UP n1,n2;  
 {  {
         int i,r;    int i,r;
   
         if ( !n1 )    if ( !n1 )
                 return n2 ? -1 : 0;      return n2 ? -1 : 0;
         else if ( !n2 )    else if ( !n2 )
                 return 1;      return 1;
         else if ( n1->d > n2->d )    else if ( n1->d > n2->d )
                 return 1;      return 1;
         else if ( n1->d < n2->d )    else if ( n1->d < n2->d )
                 return -1;      return -1;
         else {    else {
                 for ( i = n1->d; i >= 0; i-- ) {      for ( i = n1->d; i >= 0; i-- ) {
                         r = compnum(CO,n1->c[i],n2->c[i]);        r = compnum(CO,n1->c[i],n2->c[i]);
                         if ( r )        if ( r )
                                 return r;          return r;
                 }      }
                 return 0;      return 0;
         }    }
 }  }
   
 void kmulp(vl,n1,n2,nr)  void kmulp(VL vl,P n1,P n2,P *nr)
 VL vl;  
 P n1,n2;  
 P *nr;  
 {  {
         UP b1,b2,br;    UP b1,b2,br;
   
         if ( !n1 || !n2 )    if ( !n1 || !n2 )
                 *nr = 0;      *nr = 0;
         else if ( OID(n1) == O_N || OID(n2) == O_N )    else if ( OID(n1) == O_N || OID(n2) == O_N )
                 mulp(vl,n1,n2,nr);      mulp(vl,n1,n2,nr);
         else {    else {
                 ptoup(n1,&b1); ptoup(n2,&b2);      ptoup(n1,&b1); ptoup(n2,&b2);
                 kmulup(b1,b2,&br);      kmulup(b1,b2,&br);
                 uptop(br,nr);      uptop(br,nr);
         }    }
 }  }
   
 void ksquarep(vl,n1,nr)  void ksquarep(VL vl,P n1,P *nr)
 VL vl;  
 P n1;  
 P *nr;  
 {  {
         UP b1,br;    UP b1,br;
   
         if ( !n1 )    if ( !n1 )
                 *nr = 0;      *nr = 0;
         else if ( OID(n1) == O_N )    else if ( OID(n1) == O_N )
                 mulp(vl,n1,n1,nr);      mulp(vl,n1,n1,nr);
         else {    else {
                 ptoup(n1,&b1);      ptoup(n1,&b1);
                 ksquareup(b1,&br);      ksquareup(b1,&br);
                 uptop(br,nr);      uptop(br,nr);
         }    }
 }  }
   
 void kmulup(n1,n2,nr)  void kmulup(UP n1,UP n2,UP *nr)
 UP n1,n2,*nr;  
 {  {
         UP n,t,s,m,carry;    int d1,d2;
         int d,d1,d2,len,i,l;  
         pointer *r,*r0;  
   
         if ( !n1 || !n2 ) {    if ( !n1 || !n2 ) {
                 *nr = 0; return;      *nr = 0; return;
         }    }
         d1 = DEG(n1)+1; d2 = DEG(n2)+1;    d1 = DEG(n1)+1; d2 = DEG(n2)+1;
         if ( (d1 < up_kara_mag) && (d2 < up_kara_mag) )    if ( (d1 < up_kara_mag) && (d2 < up_kara_mag) )
                 mulup(n1,n2,nr);      mulup(n1,n2,nr);
         else    else
                 kmulupmain(n1,n2,nr);      kmulupmain(n1,n2,nr);
 }  }
   
 void ksquareup(n1,nr)  void ksquareup(UP n1,UP *nr)
 UP n1,*nr;  
 {  {
         int d1;    int d1;
         extern Q TWO;    extern Q TWO;
   
         if ( !n1 ) {    if ( !n1 ) {
                 *nr = 0; return;      *nr = 0; return;
         }    }
         d1 = DEG(n1)+1;    d1 = DEG(n1)+1;
         if ( (d1 < up_kara_mag) )    if ( (d1 < up_kara_mag) )
                 squareup(n1,nr);      squareup(n1,nr);
         else    else
                 ksquareupmain(n1,nr);      ksquareupmain(n1,nr);
 }  }
   
 void copyup(n1,n2)  void copyup(UP n1,UP n2)
 UP n1,n2;  
 {  {
         n2->d = n1->d;    n2->d = n1->d;
         bcopy((char *)n1->c,(char *)n2->c,(n1->d+1)*sizeof(Q));    bcopy((char *)n1->c,(char *)n2->c,(n1->d+1)*sizeof(Q));
 }  }
   
 void c_copyup(n,len,p)  void c_copyup(UP n,int len,pointer *p)
 UP n;  
 int len;  
 pointer *p;  
 {  {
         if ( n )    if ( n )
                 bcopy((char *)COEF(n),(char *)p,MIN((DEG(n)+1),len)*sizeof(pointer));      bcopy((char *)COEF(n),(char *)p,MIN((DEG(n)+1),len)*sizeof(pointer));
 }  }
   
 void kmulupmain(n1,n2,nr)  void kmulupmain(UP n1,UP n2,UP *nr)
 UP n1,n2,*nr;  
 {  {
         int d1,d2,h,len;    int d1,d2,h;
         UP n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;    UP n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
   
         d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (MAX(d1,d2)+1)/2;    d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (MAX(d1,d2)+1)/2;
         decompup(n1,h,&n1lo,&n1hi);    decompup(n1,h,&n1lo,&n1hi);
         decompup(n2,h,&n2lo,&n2hi);    decompup(n2,h,&n2lo,&n2hi);
         kmulup(n1lo,n2lo,&lo);    kmulup(n1lo,n2lo,&lo);
         kmulup(n1hi,n2hi,&hi);    kmulup(n1hi,n2hi,&hi);
         shiftup(hi,2*h,&t2);    shiftup(hi,2*h,&t2);
         addup(lo,t2,&t1);    addup(lo,t2,&t1);
         addup(hi,lo,&mid1);    addup(hi,lo,&mid1);
         subup(n1hi,n1lo,&s1); subup(n2lo,n2hi,&s2); kmulup(s1,s2,&mid2);    subup(n1hi,n1lo,&s1); subup(n2lo,n2hi,&s2); kmulup(s1,s2,&mid2);
         addup(mid1,mid2,&mid);    addup(mid1,mid2,&mid);
         shiftup(mid,h,&t2);    shiftup(mid,h,&t2);
         addup(t1,t2,nr);    addup(t1,t2,nr);
 }  }
   
 void ksquareupmain(n1,nr)  void ksquareupmain(UP n1,UP *nr)
 UP n1,*nr;  
 {  {
         int d1,h,len;    int d1,h;
         UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;    UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;
   
         d1 = DEG(n1)+1; h = (d1+1)/2;    d1 = DEG(n1)+1; h = (d1+1)/2;
         decompup(n1,h,&n1lo,&n1hi);    decompup(n1,h,&n1lo,&n1hi);
         ksquareup(n1hi,&hi); ksquareup(n1lo,&lo);    ksquareup(n1hi,&hi); ksquareup(n1lo,&lo);
         shiftup(hi,2*h,&t2);    shiftup(hi,2*h,&t2);
         addup(lo,t2,&t1);    addup(lo,t2,&t1);
         addup(hi,lo,&mid1);    addup(hi,lo,&mid1);
         subup(n1hi,n1lo,&s1); ksquareup(s1,&mid2);    subup(n1hi,n1lo,&s1); ksquareup(s1,&mid2);
         subup(mid1,mid2,&mid);    subup(mid1,mid2,&mid);
         shiftup(mid,h,&t2);    shiftup(mid,h,&t2);
         addup(t1,t2,nr);    addup(t1,t2,nr);
 }  }
   
 void rembymulup(n1,n2,nr)  void rembymulup(UP n1,UP n2,UP *nr)
 UP n1,n2;  
 UP *nr;  
 {  {
         int d1,d2,d;    int d1,d2,d;
         UP r1,r2,inv2,t,s,q;    UP r1,r2,inv2,t,s,q;
   
         if ( !n2 )    if ( !n2 )
                 error("rembymulup : division by 0");      error("rembymulup : division by 0");
         else if ( !n1 || !n2->d )    else if ( !n1 || !n2->d )
                 *nr = 0;      *nr = 0;
         else if ( (d1 = n1->d) < (d2 = n2->d) )    else if ( (d1 = n1->d) < (d2 = n2->d) )
                 *nr = n1;      *nr = n1;
         else {    else {
                 d = d1-d2;      d = d1-d2;
                 reverseup(n1,d1,&t); truncup(t,d+1,&r1);      reverseup(n1,d1,&t); truncup(t,d+1,&r1);
                 reverseup(n2,d2,&r2);      reverseup(n2,d2,&r2);
                 invmodup(r2,d,&inv2);      invmodup(r2,d,&inv2);
                 kmulup(r1,inv2,&t);      kmulup(r1,inv2,&t);
                 truncup(t,d+1,&s);      truncup(t,d+1,&s);
                 reverseup(s,d,&q);      reverseup(s,d,&q);
                 kmulup(q,n2,&t); subup(n1,t,nr);      kmulup(q,n2,&t); subup(n1,t,nr);
         }    }
 }  }
   
 /*  /*
         deg n2 = d    deg n2 = d
         deg n1 <= 2*d-1    deg n1 <= 2*d-1
         inv2 = inverse of reversep(n2) mod x^d    inv2 = inverse of reversep(n2) mod x^d
 */  */
   
 void hybrid_rembymulup_special(ff,n1,n2,inv2,nr)  void hybrid_rembymulup_special(int ff,UP n1,UP n2,UP inv2,UP *nr)
 int ff;  
 UP n1,n2,inv2;  
 UP *nr;  
 {  {
         int d1,d2,d;    int d1,d2,d;
         UP r1,t,s,q,u;    UP r1,t,s,q;
   
         if ( !n2 )    if ( !n2 )
                 error("hybrid_rembymulup : division by 0");      error("hybrid_rembymulup : division by 0");
         else if ( !n1 || !n2->d )    else if ( !n1 || !n2->d )
                 *nr = 0;      *nr = 0;
         else if ( (d1 = n1->d) < (d2 = n2->d) )    else if ( (d1 = n1->d) < (d2 = n2->d) )
                 *nr = n1;      *nr = n1;
         else {    else {
                 d = d1-d2;      d = d1-d2;
                 reverseup(n1,d1,&t); truncup(t,d+1,&r1);      reverseup(n1,d1,&t); truncup(t,d+1,&r1);
                 hybrid_tmulup(ff,r1,inv2,d+1,&t);      hybrid_tmulup(ff,r1,inv2,d+1,&t);
                 truncup(t,d+1,&s);      truncup(t,d+1,&s);
                 reverseup(s,d,&q);      reverseup(s,d,&q);
   
                 hybrid_tmulup(ff,q,n2,d2,&t);      hybrid_tmulup(ff,q,n2,d2,&t);
                 truncup(n1,d2,&s);      truncup(n1,d2,&s);
                 subup(s,t,nr);      subup(s,t,nr);
         }    }
 }  }
   
 void rembymulup_special(n1,n2,inv2,nr)  void rembymulup_special(UP n1,UP n2,UP inv2,UP *nr)
 UP n1,n2,inv2;  
 UP *nr;  
 {  {
         int d1,d2,d;    int d1,d2,d;
         UP r1,t,s,q,u;    UP r1,t,s,q;
   
         if ( !n2 )    if ( !n2 )
                 error("rembymulup : division by 0");      error("rembymulup : division by 0");
         else if ( !n1 || !n2->d )    else if ( !n1 || !n2->d )
                 *nr = 0;      *nr = 0;
         else if ( (d1 = n1->d) < (d2 = n2->d) )    else if ( (d1 = n1->d) < (d2 = n2->d) )
                 *nr = n1;      *nr = n1;
         else {    else {
                 d = d1-d2;      d = d1-d2;
                 reverseup(n1,d1,&t); truncup(t,d+1,&r1);      reverseup(n1,d1,&t); truncup(t,d+1,&r1);
                 tkmulup(r1,inv2,d+1,&t);      tkmulup(r1,inv2,d+1,&t);
                 truncup(t,d+1,&s);      truncup(t,d+1,&s);
                 reverseup(s,d,&q);      reverseup(s,d,&q);
   
                 tkmulup(q,n2,d2,&t);      tkmulup(q,n2,d2,&t);
                 truncup(n1,d2,&s);      truncup(n1,d2,&s);
                 subup(s,t,nr);      subup(s,t,nr);
         }    }
 }  }
   
 /* *nr = n1*n2 mod x^d */  /* *nr = n1*n2 mod x^d */
   
 void tkmulup(n1,n2,d,nr)  void tkmulup(UP n1,UP n2,int d,UP *nr)
 UP n1,n2;  
 int d;  
 UP *nr;  
 {  {
         int m;    int m;
         UP n1l,n1h,n2l,n2h,l,h,t,s,u,afo;    UP n1l,n1h,n2l,n2h,l,h,t,s,u;
   
         if ( d < 0 )    if ( d < 0 )
                 error("tkmulup : invalid argument");      error("tkmulup : invalid argument");
         else if ( d == 0 )    else if ( d == 0 )
                 *nr = 0;      *nr = 0;
         else {    else {
                 truncup(n1,d,&t); n1 = t;      truncup(n1,d,&t); n1 = t;
                 truncup(n2,d,&t); n2 = t;      truncup(n2,d,&t); n2 = t;
                 if ( !n1 || !n2 )      if ( !n1 || !n2 )
                         *nr = 0;        *nr = 0;
                 else if ( n1->d < up_tkara_mag || n2->d < up_tkara_mag )      else if ( n1->d < up_tkara_mag || n2->d < up_tkara_mag )
                         tmulup(n1,n2,d,nr);        tmulup(n1,n2,d,nr);
                 else {      else {
                         m = (d+1)/2;        m = (d+1)/2;
                         decompup(n1,m,&n1l,&n1h);        decompup(n1,m,&n1l,&n1h);
                         decompup(n2,m,&n2l,&n2h);        decompup(n2,m,&n2l,&n2h);
                         kmulup(n1l,n2l,&l);        kmulup(n1l,n2l,&l);
                         tkmulup(n1l,n2h,d-m,&t);        tkmulup(n1l,n2h,d-m,&t);
                         tkmulup(n2l,n1h,d-m,&s);        tkmulup(n2l,n1h,d-m,&s);
                         addup(t,s,&u);        addup(t,s,&u);
                         shiftup(u,m,&h);        shiftup(u,m,&h);
                         addup(l,h,&t);        addup(l,h,&t);
                         truncup(t,d,nr);        truncup(t,d,nr);
                 }      }
         }    }
 }  }
   
 /* n->n*x^d */  /* n->n*x^d */
   
 void shiftup(n,d,nr)  void shiftup(UP n,int d,UP *nr)
 UP n;  
 int d;  
 UP *nr;  
 {  {
         int dr;    int dr;
         UP r;    UP r;
   
         if ( !n )    if ( !n )
                 *nr = 0;      *nr = 0;
         else {    else {
                 dr = n->d+d;      dr = n->d+d;
                 *nr = r = UPALLOC(dr); r->d = dr;      *nr = r = UPALLOC(dr); r->d = dr;
                 bzero(r->c,(dr+1)*sizeof(Num));      bzero(r->c,(dr+1)*sizeof(Num));
                 bcopy(n->c,r->c+d,(n->d+1)*sizeof(Num));      bcopy(n->c,r->c+d,(n->d+1)*sizeof(Num));
         }    }
 }  }
   
 void fft_rembymulup_special(n1,n2,inv2,nr)  void fft_rembymulup_special(UP n1,UP n2,UP inv2,UP *nr)
 UP n1,n2,inv2;  
 UP *nr;  
 {  {
         int d1,d2,d;    int d1,d2,d;
         UP r1,t,s,q,u;    UP r1,t,s,q,u;
   
         if ( !n2 )    if ( !n2 )
                 error("rembymulup : division by 0");      error("rembymulup : division by 0");
         else if ( !n1 || !n2->d )    else if ( !n1 || !n2->d )
                 *nr = 0;      *nr = 0;
         else if ( (d1 = n1->d) < (d2 = n2->d) )    else if ( (d1 = n1->d) < (d2 = n2->d) )
                 *nr = n1;      *nr = n1;
         else {    else {
                 d = d1-d2;      d = d1-d2;
                 reverseup(n1,d1,&t); truncup(t,d+1,&r1);      reverseup(n1,d1,&t); truncup(t,d+1,&r1);
                 trunc_fft_mulup(r1,inv2,d+1,&t);      trunc_fft_mulup(r1,inv2,d+1,&t);
                 truncup(t,d+1,&s);      truncup(t,d+1,&s);
                 reverseup(s,d,&q);      reverseup(s,d,&q);
                 trunc_fft_mulup(q,n2,d2,&t); truncup(t,d2,&u);      trunc_fft_mulup(q,n2,d2,&t); truncup(t,d2,&u);
                 truncup(n1,d2,&s);      truncup(n1,d2,&s);
                 subup(s,u,nr);      subup(s,u,nr);
         }    }
 }  }
   
 void set_degreeup(n,d)  void set_degreeup(UP n,int d)
 UP n;  
 int d;  
 {  {
         int i;    int i;
   
         for ( i = d; i >= 0 && !n->c[i]; i-- );    for ( i = d; i >= 0 && !n->c[i]; i-- );
         n->d = i;    n->d = i;
 }  }
   
 /* n -> n0 + x^d n1 */  /* n -> n0 + x^d n1 */
   
 void decompup(n,d,n0,n1)  void decompup(UP n,int d,UP *n0,UP *n1)
 UP n;  
 int d;  
 UP *n0,*n1;  
 {  {
         int dn;    int dn;
         UP r0,r1;    UP r0,r1;
   
         if ( !n || d > n->d ) {    if ( !n || d > n->d ) {
                 *n0 = n; *n1 = 0;      *n0 = n; *n1 = 0;
         } else if ( d < 0 )    } else if ( d < 0 )
                 error("decompup : invalid argument");      error("decompup : invalid argument");
         else if ( d == 0 ) {    else if ( d == 0 ) {
                 *n0 = 0; *n1 = n;      *n0 = 0; *n1 = n;
         } else {    } else {
                 dn = n->d;      dn = n->d;
                 *n0 = r0 = UPALLOC(d-1);      *n0 = r0 = UPALLOC(d-1);
                 *n1 = r1 = UPALLOC(dn-d);      *n1 = r1 = UPALLOC(dn-d);
                 bcopy(n->c,r0->c,d*sizeof(Num));      bcopy(n->c,r0->c,d*sizeof(Num));
                 bcopy(n->c+d,r1->c,(dn-d+1)*sizeof(Num));      bcopy(n->c+d,r1->c,(dn-d+1)*sizeof(Num));
                 set_degreeup(r0,d-1);      set_degreeup(r0,d-1);
                 if ( r0->d < 0 )      if ( r0->d < 0 )
                         *n0 = 0;        *n0 = 0;
                 set_degreeup(r1,dn-d);      set_degreeup(r1,dn-d);
                 if ( r1->d < 0 )      if ( r1->d < 0 )
                         *n1 = 0;        *n1 = 0;
         }    }
 }  }
   
 /* n -> n mod x^d */  /* n -> n mod x^d */
   
 void truncup(n1,d,nr)  void truncup(UP n1,int d,UP *nr)
 UP n1;  
 int d;  
 UP *nr;  
 {  {
         int i;    int i;
         UP r;    UP r;
   
         if ( !n1 || d > n1->d )    if ( !n1 || d > n1->d )
                 *nr = n1;      *nr = n1;
         else if ( d < 0 )    else if ( d < 0 )
                 error("truncup : invalid argument");      error("truncup : invalid argument");
         else if ( d == 0 )    else if ( d == 0 )
                 *nr = 0;      *nr = 0;
         else {    else {
                 *nr = r = UPALLOC(d-1);       *nr = r = UPALLOC(d-1);
                 bcopy(n1->c,r->c,d*sizeof(Num));      bcopy(n1->c,r->c,d*sizeof(Num));
                 for ( i = d-1; i >= 0 && !r->c[i]; i-- );      for ( i = d-1; i >= 0 && !r->c[i]; i-- );
                 if ( i < 0 )      if ( i < 0 )
                         *nr = 0;        *nr = 0;
                 else      else
                         r->d = i;        r->d = i;
         }    }
 }  }
   
 int int_bits(t)  int int_bits(int t)
 int t;  
 {  {
         int k;    int k;
   
         for ( k = 0; t; t>>=1, k++);    for ( k = 0; t; t>>=1, k++);
         return k;    return k;
 }  }
   
 /* n is assumed to be LM or integer coefficient */  /* n is assumed to be LM or integer coefficient */
   
 int maxblenup(n)  int maxblenup(UP n)
 UP n;  
 {  {
         int m,r,i,d;    int m,r,i,d;
         Num *c;    Num *c;
   
         if ( !n )    if ( !n )
                 return 0;      return 0;
         d = n->d; c = (Num *)n->c;    d = n->d; c = (Num *)n->c;
         for ( r = 0, i = 0; i <= d; i++ )    for ( r = 0, i = 0; i <= d; i++ )
                 if ( c[i] ) {      if ( c[i] ) {
                         switch ( NID(c[i]) ) {        switch ( NID(c[i]) ) {
                                 case N_LM:          case N_LM:
                                         m = n_bits(((LM)c[i])->body);            m = n_bits(((LM)c[i])->body);
                                         break;            break;
                                 case N_Q:          case N_Q:
                                         m = n_bits(((Q)c[i])->nm);            m = n_bits(((Q)c[i])->nm);
                                         break;            break;
                                 default:          default:
                                         error("maxblenup : invalid coefficient");            error("maxblenup : invalid coefficient");
                         }        }
                         r = MAX(m,r);        r = MAX(m,r);
                 }      }
         return r;    return r;
 }  }
   
 void uptofmarray(mod,n,f)  void uptofmarray(int mod,UP n,ModNum *f)
 int mod;  
 UP n;  
 ModNum *f;  
 {  {
         int d,i;    int d,i;
         unsigned int r;    unsigned int r;
         Num *c;    Num *c;
   
         if ( !n )    if ( !n )
                 return;      return;
         else {    else {
                 d = n->d; c = (Num *)n->c;      d = n->d; c = (Num *)n->c;
                 for ( i = 0; i <= d; i++ ) {      for ( i = 0; i <= d; i++ ) {
                         if ( c[i] ) {        if ( c[i] ) {
                                 switch ( NID(c[i]) ) {          switch ( NID(c[i]) ) {
                                         case N_LM:            case N_LM:
                                                 f[i] = (ModNum)rem(((LM)c[i])->body,mod);              f[i] = (ModNum)rem(((LM)c[i])->body,mod);
                                                 break;              break;
                                         case N_Q:            case N_Q:
                                                 r = rem(NM((Q)c[i]),mod);              r = rem(NM((Q)c[i]),mod);
                                                 if ( r && (SGN((Q)c[i])<0) )              if ( r && (SGN((Q)c[i])<0) )
                                                         r = mod-r;                r = mod-r;
                                                 f[i] = (ModNum)r;              f[i] = (ModNum)r;
                                                 break;              break;
                                         default:            default:
                                                 error("uptofmarray : invalid coefficient");              error("uptofmarray : invalid coefficient");
                                 }          }
                         } else        } else
                                 f[i] = 0;          f[i] = 0;
                 }      }
         }    }
 }  }
   
 void fmarraytoup(f,d,nr)  void fmarraytoup(ModNum *f,int d,UP *nr)
 ModNum *f;  
 int d;  
 UP *nr;  
 {  {
         int i;    int i;
         Q *c;    Q *c;
         UP r;    UP r;
   
         if ( d < 0 ) {    if ( d < 0 ) {
                 *nr = 0;      *nr = 0;
         } else {    } else {
                 *nr = r = UPALLOC(d); c = (Q *)r->c;      *nr = r = UPALLOC(d); c = (Q *)r->c;
                 for ( i = 0; i <= d; i++ ) {      for ( i = 0; i <= d; i++ ) {
                         UTOQ((unsigned int)f[i],c[i]);        UTOQ((unsigned int)f[i],c[i]);
                 }      }
                 for ( i = d; i >= 0 && !c[i]; i-- );      for ( i = d; i >= 0 && !c[i]; i-- );
                 if ( i < 0 )      if ( i < 0 )
                         *nr = 0;        *nr = 0;
                 else      else
                         r->d = i;        r->d = i;
         }    }
 }  }
   
 /* f[i]: an array of length n */  /* f[i]: an array of length n */
   
 void uiarraytoup(f,n,d,nr)  void uiarraytoup(unsigned int **f,int n,int d,UP *nr)
 unsigned int **f;  
 int n,d;  
 UP *nr;  
 {  {
         int i,j;    int i,j;
         unsigned int *fi;    unsigned int *fi;
         N ci;    N ci;
         Q *c;    Q *c;
         UP r;    UP r;
   
         if ( d < 0 ) {    if ( d < 0 ) {
                 *nr = 0;      *nr = 0;
         } else {    } else {
                 *nr = r = UPALLOC(d); c = (Q *)r->c;      *nr = r = UPALLOC(d); c = (Q *)r->c;
                 for ( i = 0; i <= d; i++ ) {      for ( i = 0; i <= d; i++ ) {
                         fi = f[i];        fi = f[i];
                         for ( j = n-1; j >= 0 && !fi[j]; j-- );        for ( j = n-1; j >= 0 && !fi[j]; j-- );
                         j++;        j++;
                         if ( j ) {        if ( j ) {
                                 ci = NALLOC(j);          ci = NALLOC(j);
                                 PL(ci) = j;          PL(ci) = j;
                                 bcopy(fi,BD(ci),j*sizeof(unsigned int));          bcopy(fi,BD(ci),j*sizeof(unsigned int));
                                 NTOQ(ci,1,c[i]);          NTOQ(ci,1,c[i]);
                         } else        } else
                                 c[i] = 0;          c[i] = 0;
                 }      }
                 for ( i = d; i >= 0 && !c[i]; i-- );      for ( i = d; i >= 0 && !c[i]; i-- );
                 if ( i < 0 )      if ( i < 0 )
                         *nr = 0;        *nr = 0;
                 else      else
                         r->d = i;        r->d = i;
         }    }
 }  }
   
 void adj_coefup(n,m,m2,nr)  void adj_coefup(UP n,N m,N m2,UP *nr)
 UP n;  
 N m,m2;  
 UP *nr;  
 {  {
         int d;    int d;
         Q *c,*cr;    Q *c,*cr;
         Q mq;    Q mq;
         int i;    int i;
         UP r;    UP r;
   
         if ( !n )    if ( !n )
                 *nr = 0;      *nr = 0;
         else {    else {
                 d = n->d; c = (Q *)n->c;      d = n->d; c = (Q *)n->c;
                 *nr = r = UPALLOC(d); cr = (Q *)r->c;      *nr = r = UPALLOC(d); cr = (Q *)r->c;
                 NTOQ(m,1,mq);      NTOQ(m,1,mq);
                 for ( i = 0; i <= d; i++ ) {      for ( i = 0; i <= d; i++ ) {
                         if ( !c[i] )        if ( !c[i] )
                                 continue;          continue;
                         if ( cmpn(NM(c[i]),m2) > 0 )        if ( cmpn(NM(c[i]),m2) > 0 )
                                 subq(c[i],mq,&cr[i]);          subq(c[i],mq,&cr[i]);
                         else        else
                                 cr[i] = c[i];          cr[i] = c[i];
                 }      }
                 for ( i = d; i >= 0 && !cr[i]; i-- );      for ( i = d; i >= 0 && !cr[i]; i-- );
                 if ( i < 0 )      if ( i < 0 )
                         *nr = 0;        *nr = 0;
                 else      else
                         r->d = i;        r->d = i;
         }    }
 }  }
   
 /* n is assumed to have positive integer coefficients. */  /* n is assumed to have positive integer coefficients. */
   
 void remcup(n,mod,nr)  void remcup(UP n,N mod,UP *nr)
 UP n;  
 N mod;  
 UP *nr;  
 {  {
         int i,d;    int i,d;
         Q *c,*cr;    Q *c,*cr;
         UP r;    UP r;
         N t;    N t;
   
         if ( !n )    if ( !n )
                 *nr = 0;      *nr = 0;
         else {    else {
                 d = n->d; c = (Q *)n->c;      d = n->d; c = (Q *)n->c;
                 *nr = r = UPALLOC(d); cr = (Q *)r->c;      *nr = r = UPALLOC(d); cr = (Q *)r->c;
                 for ( i = 0; i <= d; i++ )      for ( i = 0; i <= d; i++ )
                         if ( c[i] ) {        if ( c[i] ) {
                                 remn(NM(c[i]),mod,&t);          remn(NM(c[i]),mod,&t);
                                 if ( t )          if ( t )
                                         NTOQ(t,1,cr[i]);            NTOQ(t,1,cr[i]);
                                 else          else
                                         cr[i] = 0;            cr[i] = 0;
                         } else        } else
                                 cr[i] = 0;          cr[i] = 0;
                 for ( i = d; i >= 0 && !cr[i]; i-- );      for ( i = d; i >= 0 && !cr[i]; i-- );
                 if ( i < 0 )      if ( i < 0 )
                         *nr = 0;        *nr = 0;
                 else      else
                         r->d = i;        r->d = i;
         }    }
 }  }
   
 void fft_mulup_main(UP,UP,int,UP *);  void fft_mulup_main(UP,UP,int,UP *);
   
 void fft_mulup(n1,n2,nr)  void fft_mulup(UP n1,UP n2,UP *nr)
 UP n1,n2;  
 UP *nr;  
 {  {
         int d1,d2,d,b1,b2,h;    int d1,d2,d,b1,b2,h;
         UP n1lo,n1hi,n2lo,n2hi,lo,hi,t1,t2,mid1,mid2,mid,s1,s2;    UP n1lo,n1hi,n2lo,n2hi,lo,hi,t1,t2,mid1,mid2,mid,s1,s2;
   
         if ( !n1 || !n2 )    if ( !n1 || !n2 )
                 *nr = 0;      *nr = 0;
         else {    else {
                 d1 = n1->d; d2 = n2->d;      d1 = n1->d; d2 = n2->d;
                 if ( MAX(d1,d2) < up_fft_mag )      if ( MAX(d1,d2) < up_fft_mag )
                         kmulup(n1,n2,nr);        kmulup(n1,n2,nr);
                 else {      else {
                         b1 = maxblenup(n1); b2 = maxblenup(n2);        b1 = maxblenup(n1); b2 = maxblenup(n2);
                         if ( fft_available(d1,b1,d2,b2) )        if ( fft_available(d1,b1,d2,b2) )
                                 fft_mulup_main(n1,n2,0,nr);          fft_mulup_main(n1,n2,0,nr);
                         else {        else {
                                 d = MAX(d1,d2)+1; h = (d+1)/2;          d = MAX(d1,d2)+1; h = (d+1)/2;
                                 decompup(n1,h,&n1lo,&n1hi);          decompup(n1,h,&n1lo,&n1hi);
                                 decompup(n2,h,&n2lo,&n2hi);          decompup(n2,h,&n2lo,&n2hi);
                                 fft_mulup(n1lo,n2lo,&lo);          fft_mulup(n1lo,n2lo,&lo);
                                 fft_mulup(n1hi,n2hi,&hi);          fft_mulup(n1hi,n2hi,&hi);
                                 shiftup(hi,2*h,&t2);          shiftup(hi,2*h,&t2);
                                 addup(lo,t2,&t1);          addup(lo,t2,&t1);
                                 addup(hi,lo,&mid1);          addup(hi,lo,&mid1);
                                 subup(n1hi,n1lo,&s1); subup(n2lo,n2hi,&s2);          subup(n1hi,n1lo,&s1); subup(n2lo,n2hi,&s2);
                                 fft_mulup(s1,s2,&mid2);          fft_mulup(s1,s2,&mid2);
                                 addup(mid1,mid2,&mid);          addup(mid1,mid2,&mid);
                                 shiftup(mid,h,&t2);          shiftup(mid,h,&t2);
                                 addup(t1,t2,nr);          addup(t1,t2,nr);
                         }        }
                 }      }
         }    }
 }  }
   
 void trunc_fft_mulup(n1,n2,dbd,nr)  void trunc_fft_mulup(UP n1,UP n2,int dbd,UP *nr)
 UP n1,n2;  
 int dbd;  
 UP *nr;  
 {  {
         int d1,d2,b1,b2,m;    int d1,d2,b1,b2,m;
         UP n1l,n1h,n2l,n2h,l,h,t,s,u;    UP n1l,n1h,n2l,n2h,l,h,t,s,u;
   
         if ( !n1 || !n2 )    if ( !n1 || !n2 )
                 *nr = 0;      *nr = 0;
         else if ( dbd < 0 )    else if ( dbd < 0 )
                 error("trunc_fft_mulup : invalid argument");      error("trunc_fft_mulup : invalid argument");
         else if ( dbd == 0 )    else if ( dbd == 0 )
                 *nr = 0;      *nr = 0;
         else {    else {
                 truncup(n1,dbd,&t); n1 = t;      truncup(n1,dbd,&t); n1 = t;
                 truncup(n2,dbd,&t); n2 = t;      truncup(n2,dbd,&t); n2 = t;
                 d1 = n1->d; d2 = n2->d;      d1 = n1->d; d2 = n2->d;
                 if ( MAX(d1,d2) < up_fft_mag )      if ( MAX(d1,d2) < up_fft_mag )
                         tkmulup(n1,n2,dbd,nr);        tkmulup(n1,n2,dbd,nr);
                 else {      else {
                         b1 = maxblenup(n1); b2 = maxblenup(n2);        b1 = maxblenup(n1); b2 = maxblenup(n2);
                         if ( fft_available(d1,b1,d2,b2) )        if ( fft_available(d1,b1,d2,b2) )
                                 fft_mulup_main(n1,n2,dbd,nr);          fft_mulup_main(n1,n2,dbd,nr);
                         else {        else {
                                 m = (dbd+1)/2;          m = (dbd+1)/2;
                                 decompup(n1,m,&n1l,&n1h);          decompup(n1,m,&n1l,&n1h);
                                 decompup(n2,m,&n2l,&n2h);          decompup(n2,m,&n2l,&n2h);
                                 fft_mulup(n1l,n2l,&l);          fft_mulup(n1l,n2l,&l);
                                 trunc_fft_mulup(n1l,n2h,dbd-m,&t);          trunc_fft_mulup(n1l,n2h,dbd-m,&t);
                                 trunc_fft_mulup(n2l,n1h,dbd-m,&s);          trunc_fft_mulup(n2l,n1h,dbd-m,&s);
                                 addup(t,s,&u);          addup(t,s,&u);
                                 shiftup(u,m,&h);          shiftup(u,m,&h);
                                 addup(l,h,&t);          addup(l,h,&t);
                                 truncup(t,dbd,nr);          truncup(t,dbd,nr);
                         }        }
                 }      }
         }    }
 }  }
   
 void fft_squareup(n1,nr)  void fft_squareup(UP n1,UP *nr)
 UP n1;  
 UP *nr;  
 {  {
         int d1,d,h,len,b1;    int d1,d,h,b1;
         UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;    UP n1lo,n1hi,hi,lo,mid1,mid2,mid,s1,t1,t2;
   
         if ( !n1 )    if ( !n1 )
                 *nr = 0;      *nr = 0;
         else {    else {
                 d1 = n1->d;      d1 = n1->d;
                 if ( d1 < up_fft_mag )      if ( d1 < up_fft_mag )
                         ksquareup(n1,nr);        ksquareup(n1,nr);
                 else {      else {
                         b1 = maxblenup(n1);        b1 = maxblenup(n1);
                         if ( fft_available(d1,b1,d1,b1) )        if ( fft_available(d1,b1,d1,b1) )
                                 fft_mulup_main(n1,n1,0,nr);          fft_mulup_main(n1,n1,0,nr);
                         else {        else {
                                 d = d1+1; h = (d1+1)/2;          d = d1+1; h = (d1+1)/2;
                                 decompup(n1,h,&n1lo,&n1hi);          decompup(n1,h,&n1lo,&n1hi);
                                 fft_squareup(n1hi,&hi);          fft_squareup(n1hi,&hi);
                                 fft_squareup(n1lo,&lo);          fft_squareup(n1lo,&lo);
                                 shiftup(hi,2*h,&t2);          shiftup(hi,2*h,&t2);
                                 addup(lo,t2,&t1);          addup(lo,t2,&t1);
                                 addup(hi,lo,&mid1);          addup(hi,lo,&mid1);
                                 subup(n1hi,n1lo,&s1);          subup(n1hi,n1lo,&s1);
                                 fft_squareup(s1,&mid2);          fft_squareup(s1,&mid2);
                                 subup(mid1,mid2,&mid);          subup(mid1,mid2,&mid);
                                 shiftup(mid,h,&t2);          shiftup(mid,h,&t2);
                                 addup(t1,t2,nr);          addup(t1,t2,nr);
                         }        }
                 }      }
         }    }
 }  }
   
 /*  /*
Line 1407  UP *nr;
Line 1330  UP *nr;
  * n1 == n2 => squaring   * n1 == n2 => squaring
  */   */
   
 void fft_mulup_main(n1,n2,dbd,nr)  void fft_mulup_main(UP n1,UP n2,int dbd,UP *nr)
 UP n1,n2;  
 UP *nr;  
 {  {
         ModNum *f1,*f2,*w,*fr;    ModNum *f1,*f2,*w,*fr;
         ModNum **frarray,**fa;    ModNum **frarray,**fa;
         unsigned int *modarray,*ma;    unsigned int *modarray,*ma;
         int modarray_len;    int modarray_len;
         int frarray_index = 0;    int frarray_index = 0;
         N m,m1,m2;    N m,m1,m2;
         int d1,d2,dmin,i,mod,root,d,cond,bound;    int d1,d2,dmin,i,mod,root,d,cond,bound;
         UP r;    UP r;
   
         if ( !n1 || !n2 ) {    if ( !n1 || !n2 ) {
                 *nr = 0; return;      *nr = 0; return;
         }    }
         d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);    d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);
         if ( !d1 || !d2 ) {    if ( !d1 || !d2 ) {
                 mulup(n1,n2,nr); return;      mulup(n1,n2,nr); return;
         }    }
         m = ONEN;    m = ONEN;
         bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+1;    bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+1;
         f1 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));    f1 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
         if ( n1 == n2 )    if ( n1 == n2 )
                 f2 = 0;      f2 = 0;
         else    else
                 f2 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));      f2 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
         w = (ModNum *)MALLOC_ATOMIC(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));    w = (ModNum *)MALLOC_ATOMIC(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));
         modarray_len = 1024;    modarray_len = 1024;
         modarray = (unsigned int *)MALLOC_ATOMIC(modarray_len*sizeof(unsigned int));    modarray = (unsigned int *)MALLOC_ATOMIC(modarray_len*sizeof(unsigned int));
         frarray = (ModNum **)MALLOC(modarray_len*sizeof(ModNum *));    frarray = (ModNum **)MALLOC(modarray_len*sizeof(ModNum *));
         for ( i = 0; i < NPrimes; i++ ) {    for ( i = 0; i < NPrimes; i++ ) {
                 FFT_primes(i,&mod,&root,&d);      FFT_primes(i,&mod,&root,&d);
                 if ( (1<<d) < d1+d2+1 )      if ( (1<<d) < d1+d2+1 )
                         continue;        continue;
                 if ( frarray_index == modarray_len ) {      if ( frarray_index == modarray_len ) {
                         ma = (unsigned int *)MALLOC_ATOMIC(2*modarray_len*sizeof(unsigned int));        ma = (unsigned int *)MALLOC_ATOMIC(2*modarray_len*sizeof(unsigned int));
                         bcopy(modarray,ma,modarray_len*sizeof(unsigned int));        bcopy(modarray,ma,modarray_len*sizeof(unsigned int));
                         modarray = ma;        modarray = ma;
                         fa = (ModNum **)MALLOC(2*modarray_len*sizeof(ModNum *));        fa = (ModNum **)MALLOC(2*modarray_len*sizeof(ModNum *));
                         bcopy(frarray,fa,modarray_len*sizeof(ModNum *));        bcopy(frarray,fa,modarray_len*sizeof(ModNum *));
                         frarray = fa;        frarray = fa;
                         modarray_len *= 2;        modarray_len *= 2;
                 }      }
                 modarray[frarray_index] = mod;      modarray[frarray_index] = mod;
                 frarray[frarray_index++] = fr      frarray[frarray_index++] = fr
                         = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));        = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
                 uptofmarray(mod,n1,f1);      uptofmarray(mod,n1,f1);
                 if ( !f2 )      if ( !f2 )
                         cond = FFT_pol_square(d1,f1,fr,i,w);        cond = FFT_pol_square(d1,f1,fr,i,w);
                 else {      else {
                         uptofmarray(mod,n2,f2);        uptofmarray(mod,n2,f2);
                         cond = FFT_pol_product(d1,f1,d2,f2,fr,i,w);        cond = FFT_pol_product(d1,f1,d2,f2,fr,i,w);
                 }      }
                 if ( cond )      if ( cond )
                         error("fft_mulup : error in FFT_pol_product");        error("fft_mulup : error in FFT_pol_product");
                 STON(mod,m1); muln(m,m1,&m2); m = m2;      STON(mod,m1); muln(m,m1,&m2); m = m2;
                 if ( n_bits(m) > bound ) {      if ( n_bits(m) > bound ) {
                         if ( !dbd )        if ( !dbd )
                                 dbd = d1+d2+1;          dbd = d1+d2+1;
                         crup(frarray,MIN(d1+d2,dbd-1),modarray,frarray_index,m,&r);        crup(frarray,MIN(d1+d2,dbd-1),modarray,frarray_index,m,&r);
                         divin(m,2,&m2);        divin(m,2,&m2);
                         adj_coefup(r,m,m2,nr);        adj_coefup(r,m,m2,nr);
                         return;        return;
                 }      }
         }    }
         error("fft_mulup : FFT_primes exhausted");    error("fft_mulup : FFT_primes exhausted");
 }  }
 #if 0  #if 0
 /* inefficient version */  /* inefficient version */
   
 void crup(f,d,mod,index,m,r)  void crup(ModNum **f,int d,int *mod,int index,N m,UP *r)
 ModNum **f;  
 int d;  
 int *mod;  
 int index;  
 N m;  
 UP *r;  
 {  {
         N *cof,*c;    N *cof,*c;
         int *inv;    int *inv;
         struct oUP w;    struct oUP w;
         int i;    int i;
         UP s,s1,s2;    UP s,s1,s2;
         N t;    N t;
         UP *g;    UP *g;
         Q q;    Q q;
         struct oEGT eg0,eg1;    struct oEGT eg0,eg1;
   
         get_eg(&eg0);    get_eg(&eg0);
         w.d = 0;    w.d = 0;
         inv = (int *)ALLOCA(index*sizeof(int));    inv = (int *)ALLOCA(index*sizeof(int));
         cof = (N *)ALLOCA(index*sizeof(N));    cof = (N *)ALLOCA(index*sizeof(N));
         c = (N *)ALLOCA(index*sizeof(N));    c = (N *)ALLOCA(index*sizeof(N));
         g = (UP *)ALLOCA(index*sizeof(UP));    g = (UP *)ALLOCA(index*sizeof(UP));
         up_lazy = 1;    up_lazy = 1;
         for ( i = 0, s = 0; i < index; i++ ) {    for ( i = 0, s = 0; i < index; i++ ) {
                 divin(m,mod[i],&cof[i]);      divin(m,mod[i],&cof[i]);
                 inv[i] = invm(rem(cof[i],mod[i]),mod[i]);      inv[i] = invm(rem(cof[i],mod[i]),mod[i]);
                 STON(inv[i],t);      STON(inv[i],t);
                 muln(cof[i],t,&c[i]);      muln(cof[i],t,&c[i]);
                 NTOQ(c[i],1,q); w.c[0] = (Num)q;      NTOQ(c[i],1,q); w.c[0] = (Num)q;
                 fmarraytoup(f[i],d,&g[i]);      fmarraytoup(f[i],d,&g[i]);
                 mulup(g[i],&w,&s1);      mulup(g[i],&w,&s1);
                 addup(s,s1,&s2);      addup(s,s1,&s2);
                 s = s2;      s = s2;
         }    }
         up_lazy = 0;    up_lazy = 0;
         remcup(s,m,r);    remcup(s,m,r);
         get_eg(&eg1);    get_eg(&eg1);
         add_eg(&eg_chrem,&eg0,&eg1);    add_eg(&eg_chrem,&eg0,&eg1);
 }  }
   
 #else  #else
 /* improved version */  /* improved version */
   
 void crup(f,d,mod,index,m,r)  void crup(ModNum **f,int d,int *mod,int index,N m,UP *r)
 ModNum **f;  
 int d;  
 int *mod;  
 int index;  
 N m;  
 UP *r;  
 {  {
         N cof,c,t,w,w1;    N cof,c,t,w,w1;
         struct oN fc;    struct oN fc;
         N *s;    N *s;
         UP u;    UP u;
         Q q;    Q q;
         int inv,i,j,rlen;    int inv,i,j,rlen;
   
         rlen = PL(m)+10; /* XXX */    rlen = PL(m)+10; /* XXX */
         PL(&fc) = 1;    PL(&fc) = 1;
         c = NALLOC(rlen);    c = NALLOC(rlen);
         w = NALLOC(rlen);    w = NALLOC(rlen);
         w1 = NALLOC(rlen);    w1 = NALLOC(rlen);
         s = (N *)MALLOC((d+1)*sizeof(N));    s = (N *)MALLOC((d+1)*sizeof(N));
         for ( i = 0; i <= d; i++ ) {    for ( i = 0; i <= d; i++ ) {
                 s[i] = NALLOC(rlen);      s[i] = NALLOC(rlen);
                 PL(s[i]) = 0;      PL(s[i]) = 0;
                 bzero(BD(s[i]),rlen*sizeof(unsigned int));      bzero(BD(s[i]),rlen*sizeof(unsigned int));
         }    }
         for ( i = 0; i < index; i++ ) {    for ( i = 0; i < index; i++ ) {
                 divin(m,mod[i],&cof);      divin(m,mod[i],&cof);
                 inv = invm(rem(cof,mod[i]),mod[i]); STON(inv,t);      inv = invm(rem(cof,mod[i]),mod[i]); STON(inv,t);
                 _muln(cof,t,c);      _muln(cof,t,c);
                 /* s += c*f[i] */      /* s += c*f[i] */
                 for ( j = 0; j <= d; j++ )      for ( j = 0; j <= d; j++ )
                         if ( f[i][j] ) {        if ( f[i][j] ) {
                                 BD(&fc)[0] = f[i][j];          BD(&fc)[0] = f[i][j];
                                 _muln(c,&fc,w);          _muln(c,&fc,w);
                                 _addn(s[j],w,w1);          _addn(s[j],w,w1);
                                 dupn(w1,s[j]);          dupn(w1,s[j]);
                         }        }
         }    }
         for ( i = d; i >= 0; i-- )    for ( i = d; i >= 0; i-- )
                 if ( s[i] )      if ( s[i] )
                         break;        break;
         if ( i < 0 )    if ( i < 0 )
                 *r = 0;      *r = 0;
         else {    else {
                 u = UPALLOC(i);      u = UPALLOC(i);
                 DEG(u) = i;      DEG(u) = i;
                 for ( j = 0; j <= i; j++ ) {      for ( j = 0; j <= i; j++ ) {
                         if ( PL(s[j]) )        if ( PL(s[j]) )
                                 NTOQ(s[j],1,q);          NTOQ(s[j],1,q);
                         else        else
                                 q = 0;          q = 0;
                         COEF(u)[j] = (Num)q;        COEF(u)[j] = (Num)q;
                 }      }
                 remcup(u,m,r);      remcup(u,m,r);
         }    }
 }  }
 #endif  #endif
   
   /*
    * dbd == 0 => n1 * n2
    * dbd > 0  => n1 * n2 mod x^dbd
    * n1 == n2 => squaring
    * return: n1*n2 mod Primes[modind[0]]*.prime...*Primes[modind[nmod-1]].prime
    */
   
   void fft_mulup_specialmod_main(UP n1,UP n2,int dbd,int *modind,int nmod,UP *nr)
   {
     ModNum *f1,*f2,*w,*fr;
     ModNum **frarray;
     N m,m1,m2;
     unsigned int *modarray;
     int d1,d2,dmin,i,root,d,cond,bound;
   
     if ( !n1 || !n2 ) {
       *nr = 0; return;
     }
     d1 = n1->d; d2 = n2->d; dmin = MIN(d1,d2);
     if ( !d1 || !d2 ) {
       mulup(n1,n2,nr); return;
     }
     m = ONEN;
     bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+1;
     f1 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
     if ( n1 == n2 )
       f2 = 0;
     else
       f2 = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
     w = (ModNum *)MALLOC_ATOMIC(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));
     frarray = (ModNum **)MALLOC(nmod*sizeof(ModNum *));
     modarray = (unsigned int *)MALLOC_ATOMIC(nmod*sizeof(unsigned int *));
   
     for ( i = 0; i < nmod; i++ ) {
       FFT_primes(modind[i],&modarray[i],&root,&d);
         if ( (1<<d) < d1+d2+1 )
           error("fft_mulup_specialmod_main : invalid modulus");
       frarray[i] = fr
         = (ModNum *)MALLOC_ATOMIC((d1+d2+1)*sizeof(ModNum));
       uptofmarray(modarray[i],n1,f1);
       if ( !f2 )
         cond = FFT_pol_square(d1,f1,fr,modind[i],w);
       else {
         uptofmarray(modarray[i],n2,f2);
         cond = FFT_pol_product(d1,f1,d2,f2,fr,modind[i],w);
       }
       if ( cond )
         error("fft_mulup_specialmod_main : error in FFT_pol_product");
       STON(modarray[i],m1); muln(m,m1,&m2); m = m2;
     }
     if ( !dbd )
       dbd = d1+d2+1;
     crup(frarray,MIN(d1+d2,dbd-1),modarray,nmod,m,nr);
   }
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.6
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>