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

Diff for /OpenXM_contrib2/asir2000/engine/up_lm.c between version 1.2 and 1.12

version 1.2, 2000/08/21 08:31:28 version 1.12, 2018/03/29 01:32:52
Line 23 
Line 23 
  * shall be made on your publication or presentation in any form of the   * shall be made on your publication or presentation in any form of the
  * results obtained by use of the SOFTWARE.   * results obtained by use of the SOFTWARE.
  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by   * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
  * e-mail at risa-admin@flab.fujitsu.co.jp of the detailed specification   * 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   * for such modification or the source code of the modified part of the
  * SOFTWARE.   * SOFTWARE.
  *   *
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/up_lm.c,v 1.1.1.1 1999/12/03 07:39:08 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/up_lm.c,v 1.11 2015/08/14 13:51:55 fujimoto Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include <math.h>  #include <math.h>
Line 54  extern GC_dont_gc;
Line 54  extern GC_dont_gc;
 extern struct oEGT eg_chrem,eg_fore,eg_back;  extern struct oEGT eg_chrem,eg_fore,eg_back;
 extern int debug_up;  extern int debug_up;
 extern int up_lazy;  extern int up_lazy;
   extern int current_ff;
   
 void crup_lm(ModNum **,int,int *,int,N,N,UP *);  void fft_mulup_lm(UP n1,UP n2,UP *nr)
   
 void fft_mulup_lm(n1,n2,nr)  
 UP n1,n2;  
 UP *nr;  
 {  {
         ModNum *f1,*f2,*w,*fr;    ModNum *f1,*f2,*w,*fr;
         ModNum *frarray[1024];    ModNum *frarray[1024];
         int modarray[1024];    int modarray[1024];
         int frarray_index = 0;    int frarray_index = 0;
         N m,m1,m2,lm_mod;    N m,m1,m2,lm_mod;
         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;
         }    }
         getmod_lm(&lm_mod);    getmod_lm(&lm_mod);
         if ( !lm_mod )    if ( !lm_mod )
                 error("fft_mulup_lm : current_mod_lm is not set");      error("fft_mulup_lm : current_mod_lm is not set");
         m = ONEN;    m = ONEN;
         bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+2;    bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+2;
         f1 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));    f1 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
         f2 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));    f2 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
         w = (ModNum *)ALLOCA(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));    w = (ModNum *)ALLOCA(6*(1<<int_bits(d1+d2+1))*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;
                 modarray[frarray_index] = mod;      modarray[frarray_index] = mod;
                 frarray[frarray_index++] = fr      frarray[frarray_index++] = fr
                         = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));        = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
                 uptofmarray(mod,n1,f1);      uptofmarray(mod,n1,f1);
                 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 ) {
 /*                      GC_dont_gc = 1; */  /*      GC_dont_gc = 1; */
                         crup_lm(frarray,d1+d2,modarray,frarray_index,m,lm_mod,&r);        crup_lm(frarray,d1+d2,modarray,frarray_index,m,lm_mod,&r);
                         uptolmup(r,nr);        uptolmup(r,nr);
                         GC_dont_gc = 0;        GC_dont_gc = 0;
                         return;        return;
                 }      }
         }    }
         error("fft_mulup : FFT_primes exhausted");    error("fft_mulup : FFT_primes exhausted");
 }  }
   
 void fft_squareup_lm(n1,nr)  void fft_squareup_lm(UP n1,UP *nr)
 UP n1;  
 UP *nr;  
 {  {
         ModNum *f1,*w,*fr;    ModNum *f1,*w,*fr;
         ModNum *frarray[1024];    ModNum *frarray[1024];
         int modarray[1024];    int modarray[1024];
         int frarray_index = 0;    int frarray_index = 0;
         N m,m1,m2,lm_mod;    N m,m1,m2,lm_mod;
         int d1,dmin,i,mod,root,d,cond,bound;    int d1,dmin,i,mod,root,d,cond,bound;
         UP r;    UP r;
   
         if ( !n1 ) {    if ( !n1 ) {
                 *nr = 0; return;      *nr = 0; return;
         }    }
         d1 = n1->d; dmin = d1;    d1 = n1->d; dmin = d1;
         if ( !d1 ) {    if ( !d1 ) {
                 mulup(n1,n1,nr); return;      mulup(n1,n1,nr); return;
         }    }
         getmod_lm(&lm_mod);    getmod_lm(&lm_mod);
         if ( !lm_mod )    if ( !lm_mod )
                 error("fft_squareup_lm : current_mod_lm is not set");      error("fft_squareup_lm : current_mod_lm is not set");
         m = ONEN;    m = ONEN;
         bound = 2*maxblenup(n1)+int_bits(d1)+2;    bound = 2*maxblenup(n1)+int_bits(d1)+2;
         f1 = (ModNum *)ALLOCA((2*d1+1)*sizeof(ModNum));    f1 = (ModNum *)ALLOCA((2*d1+1)*sizeof(ModNum));
         w = (ModNum *)ALLOCA(6*(1<<int_bits(2*d1+1))*sizeof(ModNum));    w = (ModNum *)ALLOCA(6*(1<<int_bits(2*d1+1))*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) < 2*d1+1 )      if ( (1<<d) < 2*d1+1 )
                         continue;        continue;
                 modarray[frarray_index] = mod;      modarray[frarray_index] = mod;
                 frarray[frarray_index++] = fr      frarray[frarray_index++] = fr
                         = (ModNum *)ALLOCA((2*d1+1)*sizeof(ModNum));        = (ModNum *)ALLOCA((2*d1+1)*sizeof(ModNum));
                 uptofmarray(mod,n1,f1);      uptofmarray(mod,n1,f1);
                 cond = FFT_pol_square(d1,f1,fr,i,w);      cond = FFT_pol_square(d1,f1,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 ) {
 /*                      GC_dont_gc = 1; */  /*      GC_dont_gc = 1; */
                         crup_lm(frarray,2*d1,modarray,frarray_index,m,lm_mod,&r);        crup_lm(frarray,2*d1,modarray,frarray_index,m,lm_mod,&r);
                         uptolmup(r,nr);        uptolmup(r,nr);
                         GC_dont_gc = 0;        GC_dont_gc = 0;
                         return;        return;
                 }      }
         }    }
         error("fft_squareup : FFT_primes exhausted");    error("fft_squareup : FFT_primes exhausted");
 }  }
   
 void trunc_fft_mulup_lm(n1,n2,dbd,nr)  void trunc_fft_mulup_lm(UP n1,UP n2,int dbd,UP *nr)
 UP n1,n2;  
 int dbd;  
 UP *nr;  
 {  {
         ModNum *f1,*f2,*fr,*w;    ModNum *f1,*f2,*fr,*w;
         ModNum *frarray[1024];    ModNum *frarray[1024];
         int modarray[1024];    int modarray[1024];
         int frarray_index = 0;    int frarray_index = 0;
         N m,m1,m2,lm_mod;    N m,m1,m2,lm_mod;
         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 ) {
                 tmulup(n1,n2,dbd,nr); return;      tmulup(n1,n2,dbd,nr); return;
         }    }
         getmod_lm(&lm_mod);    getmod_lm(&lm_mod);
         if ( !lm_mod )    if ( !lm_mod )
                 error("trunc_fft_mulup_lm : current_mod_lm is not set");      error("trunc_fft_mulup_lm : current_mod_lm is not set");
         m = ONEN;    m = ONEN;
         bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+2;    bound = maxblenup(n1)+maxblenup(n2)+int_bits(dmin)+2;
         f1 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));    f1 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
         f2 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));    f2 = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
         w = (ModNum *)ALLOCA(6*(1<<int_bits(d1+d2+1))*sizeof(ModNum));    w = (ModNum *)ALLOCA(6*(1<<int_bits(d1+d2+1))*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;
   
                 modarray[frarray_index] = mod;      modarray[frarray_index] = mod;
                 frarray[frarray_index++] = fr      frarray[frarray_index++] = fr
                         = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));        = (ModNum *)ALLOCA((d1+d2+1)*sizeof(ModNum));
                 uptofmarray(mod,n1,f1);      uptofmarray(mod,n1,f1);
                 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 ) {
 /*                      GC_dont_gc = 1; */  /*      GC_dont_gc = 1; */
                         crup_lm(frarray,MIN(dbd-1,d1+d2),modarray,frarray_index,m,lm_mod,&r);        crup_lm(frarray,MIN(dbd-1,d1+d2),modarray,frarray_index,m,lm_mod,&r);
                         uptolmup(r,nr);        uptolmup(r,nr);
                         GC_dont_gc = 0;        GC_dont_gc = 0;
                         return;        return;
                 }      }
         }    }
         error("trunc_fft_mulup : FFT_primes exhausted");    error("trunc_fft_mulup : FFT_primes exhausted");
 }  }
   
 void crup_lm(f,d,mod,index,m,lm_mod,r)  void crup_lm(ModNum **f,int d,int *mod,int index,N m,N lm_mod,UP *r)
 ModNum **f;  
 int d;  
 int *mod;  
 int index;  
 N m;  
 N lm_mod;  
 UP *r;  
 {  {
         double *k;    double *k;
         double c2;    double c2;
         unsigned int *w;    unsigned int *w;
         unsigned int zi,au,al;    unsigned int zi,au,al;
         UL a;    UL a;
         int i,j,l,len;    int i,j,l,len;
         UP s,s1,u;    UP s,s1,u;
         struct oUP c;    struct oUP c;
         N t,ci,mi,qn;    N t,ci,mi,qn;
         unsigned int **sum;    unsigned int **sum;
         unsigned int *sum_b;    unsigned int *sum_b;
         Q q;    Q q;
         struct oEGT eg0,eg1;  
   
         if ( !lm_mod )    if ( !lm_mod )
                 error("crup_lm : current_mod_lm is not set");      error("crup_lm : current_mod_lm is not set");
         k = (double *)ALLOCA((d+1)*sizeof(double));    k = (double *)ALLOCA((d+1)*sizeof(double));
         for ( j = 0; j <= d; j++ )    for ( j = 0; j <= d; j++ )
                 k[j] = 0.5;      k[j] = 0.5;
         up_lazy = 1;    up_lazy = 1;
         sum = (unsigned int **)ALLOCA((d+1)*sizeof(unsigned int *));    sum = (unsigned int **)ALLOCA((d+1)*sizeof(unsigned int *));
         len = (int_bits(index)+n_bits(lm_mod)+31)/32+1;    len = (int_bits(index)+n_bits(lm_mod)+31)/32+1;
         w = (unsigned int *)ALLOCA((len+1)*sizeof(unsigned int));    w = (unsigned int *)ALLOCA((len+1)*sizeof(unsigned int));
         sum_b = (unsigned int *)MALLOC_ATOMIC((d+1)*len*sizeof(unsigned int));    sum_b = (unsigned int *)MALLOC_ATOMIC((d+1)*len*sizeof(unsigned int));
         for ( j = 0; j <= d; j++ ) {    for ( j = 0; j <= d; j++ ) {
                 sum[j] = sum_b+len*j;      sum[j] = sum_b+len*j;
                 bzero((char *)sum[j],len*sizeof(unsigned int));      bzero((char *)sum[j],len*sizeof(unsigned int));
         }    }
         for ( i = 0, s = 0; i < index; i++ ) {    for ( i = 0, s = 0; i < index; i++ ) {
                 divin(m,mod[i],&ci);      divin(m,mod[i],&ci);
                 zi = (unsigned int)invm((unsigned int)rem(ci,mod[i]),mod[i]);      zi = (unsigned int)invm((unsigned int)rem(ci,mod[i]),mod[i]);
   
                 STON(zi,t); muln(ci,t,&mi);      STON(zi,t); muln(ci,t,&mi);
                 divn(mi,lm_mod,&qn,&t);      divn(mi,lm_mod,&qn,&t);
                 if ( t )      if ( t )
                         for ( j = 0; j <= d; j++ ) {        for ( j = 0; j <= d; j++ ) {
                                 bzero((char *)w,(len+1)*sizeof(unsigned int));          bzero((char *)w,(len+1)*sizeof(unsigned int));
                                 muln_1(BD(t),PL(t),(unsigned int)f[i][j],w);          muln_1(BD(t),PL(t),(unsigned int)f[i][j],w);
                                 for ( l = PL(t); l >= 0 && !w[l]; l--);          for ( l = PL(t); l >= 0 && !w[l]; l--);
                                 l++;          l++;
                                 if ( l )          if ( l )
                                         addarray_to(w,l,sum[j],len);            addarray_to(w,l,sum[j],len);
                         }        }
                 c2 = (double)zi/(double)mod[i];      c2 = (double)zi/(double)mod[i];
                 for ( j = 0; j <= d; j++ )      for ( j = 0; j <= d; j++ )
                         k[j] += c2*f[i][j];        k[j] += c2*f[i][j];
         }    }
         uiarraytoup(sum,len,d,&s);    uiarraytoup(sum,len,d,&s);
         GC_free(sum_b);    GCFREE(sum_b);
   
         u = UPALLOC(d);    u = UPALLOC(d);
         for ( j = 0; j <= d; j++ ) {    for ( j = 0; j <= d; j++ ) {
 #if 1  #if 1
                 a = (UL)floor(k[j]);      a = (UL)floor(k[j]);
 #if defined(i386) || defined(__alpha) || defined(VISUAL)  #if defined(i386) || defined(__alpha) || defined(VISUAL) || defined(__MINGW32__) || defined(__x86_64)
                 au = ((unsigned int *)&a)[1];      au = ((unsigned int *)&a)[1];
                 al = ((unsigned int *)&a)[0];      al = ((unsigned int *)&a)[0];
 #else  #else
                 al = ((unsigned int *)&a)[1];      al = ((unsigned int *)&a)[1];
                 au = ((unsigned int *)&a)[0];      au = ((unsigned int *)&a)[0];
 #endif  #endif
                 if ( au ) {      if ( au ) {
                         NEWQ(q); SGN(q) = 1; NM(q)=NALLOC(2); DN(q)=0;        NEWQ(q); SGN(q) = 1; NM(q)=NALLOC(2); DN(q)=0;
                         PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;        PL(NM(q))=2; BD(NM(q))[0]=al; BD(NM(q))[1] = au;
                 } else      } else
                         UTOQ(al,q);        UTOQ(al,q);
 #else  #else
                 al = (int)floor(k[j]); STOQ(al,q);      al = (int)floor(k[j]); STOQ(al,q);
 #endif  #endif
                  u->c[j] = (Num)q;       u->c[j] = (Num)q;
         }    }
         for ( j = d; j >= 0 && !u->c[j]; j-- );    for ( j = d; j >= 0 && !u->c[j]; j-- );
         if ( j < 0 )    if ( j < 0 )
                 u = 0;      u = 0;
         else    else
                 u->d = j;      u->d = j;
         divn(m,lm_mod,&qn,&t); NTOQ(t,-1,q);    divn(m,lm_mod,&qn,&t); NTOQ(t,-1,q);
         c.d = 0; c.c[0] = (Num)q;    c.d = 0; c.c[0] = (Num)q;
         mulup(u,&c,&s1);    mulup(u,&c,&s1);
         up_lazy = 0;    up_lazy = 0;
   
         addup(s,s1,r);    addup(s,s1,r);
 }  }
   
 void fft_rembymulup_special_lm(n1,n2,inv2,nr)  void fft_rembymulup_special_lm(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_lm(r1,inv2,d+1,&t);      trunc_fft_mulup_lm(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_lm(q,n2,d2,&t); truncup(t,d2,&u);      trunc_fft_mulup_lm(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 uptolmup(n,nr)  void uptolmup(UP n,UP *nr)
 UP n;  
 UP *nr;  
 {  {
         int i,d;    int i,d;
         Q *c;    Q *c;
         LM *cr;    LM *cr;
         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 = (LM *)r->c;      *nr = r = UPALLOC(d); cr = (LM *)r->c;
                 for ( i = 0; i <= d; i++ )      for ( i = 0; i <= d; i++ )
                         qtolm(c[i],&cr[i]);        qtolm(c[i],&cr[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;
         }    }
 }  }
   
 save_up(obj,name)  void save_up(UP obj,char *name)
 UP obj;  
 char *name;  
 {  {
         P p;    P p;
         Obj ret;    Obj ret;
         NODE n0,n1;    NODE n0,n1;
         STRING s;    STRING s;
     void Pbsave();
   
         uptop(obj,&p);    uptop(obj,&p);
         MKSTR(s,name);    MKSTR(s,name);
         MKNODE(n1,s,0); MKNODE(n0,p,n1);    MKNODE(n1,s,0); MKNODE(n0,p,n1);
         Pbsave(n0,&ret);    Pbsave(n0,&ret);
 }  }
   
 void hybrid_powermodup(f,xp)  void hybrid_powermodup(UP f,UP *xp)
 UP f;  
 UP *xp;  
 {  {
         N n;    N n;
         UP x,y,t,invf,s;    UP x,y,t,invf,s;
         int k;    int k;
         LM lm;    LM lm;
         struct oEGT eg_sq,eg_rem,eg_mul,eg_inv,eg0,eg1,eg2,eg3;  
         char name[BUFSIZ];  
   
         getmod_lm(&n);    getmod_lm(&n);
         if ( !n )    if ( !n )
                 error("hybrid_powermodup : current_mod_lm is not set");      error("hybrid_powermodup : current_mod_lm is not set");
         MKLM(ONEN,lm);    MKLM(ONEN,lm);
         x = UPALLOC(1); x->d = 1; x->c[1] = (Num)lm;    x = UPALLOC(1); x->d = 1; x->c[1] = (Num)lm;
         y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;    y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;
   
         reverseup(f,f->d,&t);    reverseup(f,f->d,&t);
         invmodup(t,f->d,&s); uptolmup(s,&invf);    invmodup(t,f->d,&s); uptolmup(s,&invf);
         for ( k = n_bits(n)-1; k >= 0; k-- ) {    for ( k = n_bits(n)-1; k >= 0; k-- ) {
                 hybrid_squareup(FF_GFP,y,&t);      hybrid_squareup(FF_GFP,y,&t);
                 hybrid_rembymulup_special(FF_GFP,t,f,invf,&s);      hybrid_rembymulup_special(FF_GFP,t,f,invf,&s);
                 y = s;      y = s;
                 if ( n->b[k/32] & (1<<(k%32)) ) {      if ( n->b[k/32] & (1<<(k%32)) ) {
                         mulup(y,x,&t);        mulup(y,x,&t);
                         remup(t,f,&s);        remup(t,f,&s);
                         y = s;        y = s;
                 }      }
         }    }
         *xp = y;    *xp = y;
 }  }
   
 void powermodup(f,xp)  void powermodup(UP f,UP *xp)
 UP f;  
 UP *xp;  
 {  {
         N n;    N n;
         UP x,y,t,invf,s;    UP x,y,t,invf,s;
         int k;    int k;
         LM lm;    Num c;
         struct oEGT eg_sq,eg_rem,eg_mul,eg_inv,eg0,eg1,eg2,eg3;  
   
         field_order_ff(&n);    if ( !current_ff )
         if ( !n )      error("powermodup : current_ff is not set");
                 error("powermodup : current_mod_lm is not set");    field_order_ff(&n);
         MKLM(ONEN,lm);    one_ff(&c);
         x = UPALLOC(1); x->d = 1; x->c[1] = (Num)lm;    x = UPALLOC(1); x->d = 1; x->c[1] = c;
         y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;    y = UPALLOC(0); y->d = 0; y->c[0] = c;
   
         reverseup(f,f->d,&t);    reverseup(f,f->d,&t);
         invmodup(t,f->d,&s); uptolmup(s,&invf);    invmodup(t,f->d,&s);
         for ( k = n_bits(n)-1; k >= 0; k-- ) {    switch ( current_ff ) {
                 ksquareup(y,&t);      case FF_GFP:
                 rembymulup_special(t,f,invf,&s);      case FF_GFPN:
                 y = s;        uptolmup(s,&invf);
                 if ( n->b[k/32] & (1<<(k%32)) ) {        break;
                         mulup(y,x,&t);      case FF_GFS:
                         remup(t,f,&s);      case FF_GFSN:
                         y = s;        invf = s; /* XXX */
                 }        break;
         }      default:
         *xp = y;        error("powermodup : not implemented yet");
     }
     for ( k = n_bits(n)-1; k >= 0; k-- ) {
       ksquareup(y,&t);
       rembymulup_special(t,f,invf,&s);
       y = s;
       if ( n->b[k/32] & (1<<(k%32)) ) {
         mulup(y,x,&t);
         remup(t,f,&s);
         y = s;
       }
     }
     *xp = y;
 }  }
   
 /* g^d mod f */  /* g^d mod f */
   
 void hybrid_generic_powermodup(g,f,d,xp)  void hybrid_generic_powermodup(UP g,UP f,Q d,UP *xp)
 UP g,f;  
 Q d;  
 UP *xp;  
 {  {
         N e;    N e;
         UP x,y,t,invf,s;    UP x,y,t,invf,s;
         int k;    int k;
         LM lm;    LM lm;
         struct oEGT eg_sq,eg_rem,eg_mul,eg_inv,eg0,eg1,eg2,eg3;  
   
         e = NM(d);    e = NM(d);
         MKLM(ONEN,lm);    MKLM(ONEN,lm);
         y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;    y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;
         remup(g,f,&x);    remup(g,f,&x);
         if ( !x ) {    if ( !x ) {
                 *xp = !d ? y : 0;      *xp = !d ? y : 0;
                 return;      return;
         } else if ( !x->d ) {    } else if ( !x->d ) {
                 pwrup(x,d,xp);      pwrup(x,d,xp);
                 return;      return;
         }    }
         reverseup(f,f->d,&t);    reverseup(f,f->d,&t);
         invmodup(t,f->d,&invf);    invmodup(t,f->d,&invf);
         for ( k = n_bits(e)-1; k >= 0; k-- ) {    for ( k = n_bits(e)-1; k >= 0; k-- ) {
                 hybrid_squareup(FF_GFP,y,&t);      hybrid_squareup(FF_GFP,y,&t);
                 hybrid_rembymulup_special(FF_GFP,t,f,invf,&s);      hybrid_rembymulup_special(FF_GFP,t,f,invf,&s);
                 y = s;      y = s;
                 if ( e->b[k/32] & (1<<(k%32)) ) {      if ( e->b[k/32] & (1<<(k%32)) ) {
                         mulup(y,x,&t);        mulup(y,x,&t);
                         remup(t,f,&s);        remup(t,f,&s);
                         y = s;        y = s;
                 }      }
         }    }
         *xp = y;    *xp = y;
 }  }
   
 void generic_powermodup(g,f,d,xp)  void generic_powermodup(UP g,UP f,Q d,UP *xp)
 UP g,f;  
 Q d;  
 UP *xp;  
 {  {
         N e;    N e;
         UP x,y,t,invf,s;    UP x,y,t,invf,s;
         int k;    int k;
         LM lm;    Num c;
         struct oEGT eg_sq,eg_rem,eg_mul,eg_inv,eg0,eg1,eg2,eg3;  
   
         e = NM(d);    e = NM(d);
         MKLM(ONEN,lm);    one_ff(&c);
         y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;    y = UPALLOC(0); y->d = 0; y->c[0] = c;
         remup(g,f,&x);    remup(g,f,&x);
         if ( !x ) {    if ( !x ) {
                 *xp = !d ? y : 0;      *xp = !d ? y : 0;
                 return;      return;
         } else if ( !x->d ) {    } else if ( !x->d ) {
                 pwrup(x,d,xp);      pwrup(x,d,xp);
                 return;      return;
         }    }
         reverseup(f,f->d,&t);    reverseup(f,f->d,&t);
         invmodup(t,f->d,&invf);    invmodup(t,f->d,&invf);
         for ( k = n_bits(e)-1; k >= 0; k-- ) {    for ( k = n_bits(e)-1; k >= 0; k-- ) {
                 ksquareup(y,&t);      ksquareup(y,&t);
                 rembymulup_special(t,f,invf,&s);      rembymulup_special(t,f,invf,&s);
                 y = s;      y = s;
                 if ( e->b[k/32] & (1<<(k%32)) ) {      if ( e->b[k/32] & (1<<(k%32)) ) {
                         mulup(y,x,&t);        mulup(y,x,&t);
                         remup(t,f,&s);        remup(t,f,&s);
                         y = s;        y = s;
                 }      }
         }    }
         *xp = y;    *xp = y;
 }  }
   
 void hybrid_powertabup(f,xp,tab)  void hybrid_powertabup(UP f,UP xp,UP *tab)
 UP f;  
 UP xp;  
 UP *tab;  
 {  {
         UP y,t,invf;    UP y,t,invf;
         int i,d;    int i,d;
         LM lm;    LM lm;
         struct oEGT eg_rem,eg_mul,eg0,eg1,eg2;  
   
         d = f->d;    d = f->d;
         MKLM(ONEN,lm);    MKLM(ONEN,lm);
         y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;    y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;
         tab[0] = y;    tab[0] = y;
         tab[1] = xp;    tab[1] = xp;
   
         reverseup(f,f->d,&t);    reverseup(f,f->d,&t);
         invmodup(t,f->d,&invf);    invmodup(t,f->d,&invf);
   
         for ( i = 2; i < d; i++ ) {    for ( i = 2; i < d; i++ ) {
                 if ( debug_up )      if ( debug_up ){
                         fprintf(stderr,".");        fprintf(stderr,".");
                 hybrid_mulup(FF_GFP,tab[i-1],xp,&t);      }
                 hybrid_rembymulup_special(FF_GFP,t,f,invf,&tab[i]);      hybrid_mulup(FF_GFP,tab[i-1],xp,&t);
         }      hybrid_rembymulup_special(FF_GFP,t,f,invf,&tab[i]);
     }
 }  }
   
 void powertabup(f,xp,tab)  void powertabup(UP f,UP xp,UP *tab)
 UP f;  
 UP xp;  
 UP *tab;  
 {  {
         UP y,t,invf;    UP y,t,invf;
         int i,d;    int i,d;
         LM lm;    Num c;
         struct oEGT eg_rem,eg_mul,eg0,eg1,eg2;  
   
         d = f->d;    d = f->d;
         MKLM(ONEN,lm);    one_ff(&c);
         y = UPALLOC(0); y->d = 0; y->c[0] = (Num)lm;    y = UPALLOC(0); y->d = 0; y->c[0] = c;
         tab[0] = y;    tab[0] = y;
         tab[1] = xp;    tab[1] = xp;
   
         reverseup(f,f->d,&t);    reverseup(f,f->d,&t);
         invmodup(t,f->d,&invf);    invmodup(t,f->d,&invf);
   
         for ( i = 2; i < d; i++ ) {    for ( i = 2; i < d; i++ ) {
                 if ( debug_up )      if ( debug_up ){
                         fprintf(stderr,".");        fprintf(stderr,".");
                 kmulup(tab[i-1],xp,&t);      }
                 rembymulup_special(t,f,invf,&tab[i]);      kmulup(tab[i-1],xp,&t);
         }      rembymulup_special(t,f,invf,&tab[i]);
     }
 }  }
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.12
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>