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

Diff for /OpenXM_contrib2/asir2000/fft/dft.c between version 1.4 and 1.5

version 1.4, 2003/02/14 22:29:10 version 1.5, 2018/03/29 01:32:53
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/fft/dft.c,v 1.3 2000/08/22 05:04:12 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/fft/dft.c,v 1.4 2003/02/14 22:29:10 ohara Exp $
 */  */
 /*  /*
  *              store                   arith           modulus   *    store      arith    modulus
  *      -------------------------------------------------------   *  -------------------------------------------------------
  *      IEEE:   float (23+1)    <->     double (52+1)   24 bits   *  IEEE:  float (23+1)  <->  double (52+1)  24 bits
  *              int   (32)      <->     double (52+1)   26 bits   *    int   (32)  <->  double (52+1)  26 bits
  *      370:    float (24)      <->     double (56)     24/25 bits   *  370:  float (24)  <->  double (56)  24/25 bits
  *              int   (32)      <->     double (56)     28 bits   *    int   (32)  <->  double (56)  28 bits
  */   */
   
   
Line 71 
Line 71 
   
 void C_DFT_FORE( in, nin, i1, K, powa,  void C_DFT_FORE( in, nin, i1, K, powa,
 #ifdef POWA_STRIDE  #ifdef POWA_STRIDE
                  a1,       a1,
 #endif  #endif
                  out, o1, P, Pinv, wk )       out, o1, P, Pinv, wk )
 ModNum in[], powa[], out[], wk[];  ModNum in[], powa[], out[], wk[];
 ModNum P;  ModNum P;
 int nin, i1, K, o1;  int nin, i1, K, o1;
Line 102  double Pinv;
Line 102  double Pinv;
   ModNum *qi, *qo, *qis, *qos, a, t0, tn;    ModNum *qi, *qo, *qis, *qos, a, t0, tn;
   
   /*    /*
    *    for k = K-1, ..., 0 do the following:     *  for k = K-1, ..., 0 do the following:
    *     *
    *        Out[2^{K-k}*i + j]             <= t0 + \alpha^{2^k * j} * tn;     *      Out[2^{K-k}*i + j]       <= t0 + \alpha^{2^k * j} * tn;
    *        Out[2^{K-k}*i + j + 2^{K-k-1}] <= t0 - \alpha^{2^k * j} * tn;     *      Out[2^{K-k}*i + j + 2^{K-k-1}] <= t0 - \alpha^{2^k * j} * tn;
    *     *
    *      where t0 = In[2^{K-k-1}*i + j] and     *      where t0 = In[2^{K-k-1}*i + j] and
    *            tn = In[2^{K-k-1}*i + j + 2^{K-1}],     *    tn = In[2^{K-k-1}*i + j + 2^{K-1}],
    *      for 0 <= i < 2^k and 0 <= j < 2^{K-k-1}.     *      for 0 <= i < 2^k and 0 <= j < 2^{K-k-1}.
    *     *
    */     */
Line 127  double Pinv;
Line 127  double Pinv;
     for ( i = nin; i > 0; i--, qi += i1, qo += ostep ) qo[0] = qo[odist] = qi[0];      for ( i = nin; i > 0; i--, qi += i1, qo += ostep ) qo[0] = qo[odist] = qi[0];
     for ( i = n - nin; i > 0; i--, qo += ostep ) qo[0] = qo[odist] = 0;      for ( i = n - nin; i > 0; i--, qo += ostep ) qo[0] = qo[odist] = 0;
   } else {    } else {
     idist = i1 << k;    /* distance = 2^{K-1} */      idist = i1 << k;  /* distance = 2^{K-1} */
     for ( j = i = nin - n; i > 0; i--, qi += i1, qo += ostep ) {      for ( j = i = nin - n; i > 0; i--, qi += i1, qo += ostep ) {
       t0 = qi[0],  tn = qi[idist];        t0 = qi[0],  tn = qi[idist];
       if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;        if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
       else {        else {
         qo[0] = AplusBmodP( t0, tn, P );    qo[0] = AplusBmodP( t0, tn, P );
         qo[odist] = A_BmodP( t0, tn, P );    qo[odist] = A_BmodP( t0, tn, P );
       }        }
     }      }
     for ( i = n - j; i > 0; i--, qi += i1, qo += ostep )      for ( i = n - j; i > 0; i--, qi += i1, qo += ostep )
Line 145  double Pinv;
Line 145  double Pinv;
 #ifndef OPT  #ifndef OPT
     K_k = (K_k_1 = K_k) + 1;      K_k = (K_k_1 = K_k) + 1;
 #endif  #endif
     n >>= 1;            /* == 2^k */      n >>= 1;    /* == 2^k */
     Ndiv2n <<= 1;       /* == 2^{K-k-1} */      Ndiv2n <<= 1;  /* == 2^{K-k-1} */
     if ( inwk )      if ( inwk )
       qi = wk, qo = out, iostep = 1, oostep = o1, idist = idistwk;        qi = wk, qo = out, iostep = 1, oostep = o1, idist = idistwk;
     else      else
       qi = out, qo = wk, iostep = o1, oostep = 1, idist = idistout;        qi = out, qo = wk, iostep = o1, oostep = 1, idist = idistout;
     qis = qi,  qos = qo;      qis = qi,  qos = qo;
     /**/      /**/
     istep = ostep;      /* = iostep << K_k_1; */      istep = ostep;  /* = iostep << K_k_1; */
 #ifndef OPT  #ifndef OPT
     ostep = oostep << K_k;      ostep = oostep << K_k;
     odist = oostep << K_k_1;      odist = oostep << K_k_1;
Line 167  double Pinv;
Line 167  double Pinv;
       t0 = qi[0],  tn = qi[idist];        t0 = qi[0],  tn = qi[idist];
       if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;        if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
       else {        else {
         qo[0] = AplusBmodP( t0, tn, P );    qo[0] = AplusBmodP( t0, tn, P );
         qo[odist] = A_BmodP( t0, tn, P );    qo[odist] = A_BmodP( t0, tn, P );
       }        }
     }      }
     /**/      /**/
Line 181  double Pinv;
Line 181  double Pinv;
 #endif  #endif
       qi = (qis += iostep),  qo = (qos += oostep);        qi = (qis += iostep),  qo = (qos += oostep);
       for ( i = 0; i < n; i++, qi += istep, qo += ostep ) {        for ( i = 0; i < n; i++, qi += istep, qo += ostep ) {
         t0 = qi[0],  tn = qi[idist];    t0 = qi[0],  tn = qi[idist];
         if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;    if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
         else {    else {
           AxBmodP( tn, ModNum, tn, a, P, Pinv );      AxBmodP( tn, ModNum, tn, a, P, Pinv );
           qo[0] = AplusBmodP( t0, tn, P );      qo[0] = AplusBmodP( t0, tn, P );
           qo[odist] = A_BmodP( t0, tn, P );      qo[odist] = A_BmodP( t0, tn, P );
         }    }
       }        }
     }      }
   }    }
   /*    /*
    *  When k = 0, for i = 0 (0 <= i < 2^0=1) and 0 <= j < 2^{K-1},     *  When k = 0, for i = 0 (0 <= i < 2^0=1) and 0 <= j < 2^{K-1},
    *    Out[j] <= t0 + \alpha^{j}*tn and Out[j+2^{K-1}] <= t0 - \alpha^{j}*tn,     *  Out[j] <= t0 + \alpha^{j}*tn and Out[j+2^{K-1}] <= t0 - \alpha^{j}*tn,
    *  where t0 = In[j] and tn = In[j + 2^{K-1}].     *  where t0 = In[j] and tn = In[j + 2^{K-1}].
    */     */
   qi = wk, qo = out, idist = idistwk, odist = idistout;    qi = wk, qo = out, idist = idistwk, odist = idistout;
Line 228  double Pinv;
Line 228  double Pinv;
 #ifndef OPT  #ifndef OPT
     K_k = (K_k_1 = K_k) + 1;      K_k = (K_k_1 = K_k) + 1;
 #endif  #endif
     n >>= 1;            /* == 2^k */      n >>= 1;    /* == 2^k */
     Ndiv2n <<= 1;       /* == 2^{K-k-1} */      Ndiv2n <<= 1;  /* == 2^{K-k-1} */
     if ( inwk ) qi = wk, qo = out, istep = 1, ostep = o1, idist = idistwk;      if ( inwk ) qi = wk, qo = out, istep = 1, ostep = o1, idist = idistwk;
     else qi = out, qo = wk, istep = o1, ostep = 1, idist = idistout;      else qi = out, qo = wk, istep = o1, ostep = 1, idist = idistout;
     qis = qi,  qos = qo;      qis = qi,  qos = qo;
     /**/      /**/
     iostep = oostep;    /* = istep << K_k_1; */      iostep = oostep;  /* = istep << K_k_1; */
 #ifndef OPT  #ifndef OPT
     oostep = ostep << K_k;      oostep = ostep << K_k;
     odist = ostep << K_k_1;      odist = ostep << K_k_1;
Line 248  double Pinv;
Line 248  double Pinv;
       t0 = qis[0],  tn = qis[idist];        t0 = qis[0],  tn = qis[idist];
       if ( tn == (ModNum)0 ) qos[0] = qos[odist] = t0;        if ( tn == (ModNum)0 ) qos[0] = qos[odist] = t0;
       else {        else {
         qos[0] = AplusBmodP( t0, tn, P );    qos[0] = AplusBmodP( t0, tn, P );
         qos[odist] = A_BmodP( t0, tn, P );    qos[odist] = A_BmodP( t0, tn, P );
       }        }
 #ifdef POWA_STRIDE  #ifdef POWA_STRIDE
       for ( aj = a1, j = 1, qi = qis, qo = qos; j < Ndiv2n; aj += a1, j++ ) {        for ( aj = a1, j = 1, qi = qis, qo = qos; j < Ndiv2n; aj += a1, j++ ) {
 #else  #else
       for ( j = 1, qi = qis, qo = qos; j < Ndiv2n; j++ ) {        for ( j = 1, qi = qis, qo = qos; j < Ndiv2n; j++ ) {
 #endif  #endif
         qi += istep,  qo += ostep;    qi += istep,  qo += ostep;
         t0 = qi[0],  tn = qi[idist];    t0 = qi[0],  tn = qi[idist];
         if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;    if ( tn == (ModNum)0 ) qo[0] = qo[odist] = t0;
         else {    else {
 #ifdef POWA_STRIDE  #ifdef POWA_STRIDE
           a = powa[aj << k];      a = powa[aj << k];
 #else  #else
           a = powa[j << k];      a = powa[j << k];
 #endif  #endif
           AxBmodP( tn, ModNum, tn, a, P, Pinv );      AxBmodP( tn, ModNum, tn, a, P, Pinv );
           qo[0] = AplusBmodP( t0, tn, P );      qo[0] = AplusBmodP( t0, tn, P );
           qo[odist] = A_BmodP( t0, tn, P );      qo[odist] = A_BmodP( t0, tn, P );
         }    }
       }        }
     }      }
   }    }
Line 280  double Pinv;
Line 280  double Pinv;
  *  NOTE:   *  NOTE:
  *  (1) Let \alpha be a primitive N-th root of unity modulo P,   *  (1) Let \alpha be a primitive N-th root of unity modulo P,
  *  where N is even,   *  where N is even,
  *      and t be an integer s.t. 0 <= t < N/2.   *  and t be an integer s.t. 0 <= t < N/2.
  *      Then, 1/\alpha^t is given by (- \alpha^s) with s = N/2 - t.   *  Then, 1/\alpha^t is given by (- \alpha^s) with s = N/2 - t.
  *   *
  *  (2) Let P be a prime s.t. P = 2^n*Q + 1, where Q is an odd integer.   *  (2) Let P be a prime s.t. P = 2^n*Q + 1, where Q is an odd integer.
  *      Then, for 0 < s <= n, 2^{-s} \equiv - 2^{n-s}*Q \bmod P,   *  Then, for 0 < s <= n, 2^{-s} \equiv - 2^{n-s}*Q \bmod P,
  *      because 1 \equiv - 2^n*Q \bmod P.   *  because 1 \equiv - 2^n*Q \bmod P.
  *   *
  **********************************************************************/   **********************************************************************/
   
 void C_DFT_BACK( in, N, i1, log2N, powa,  void C_DFT_BACK( in, N, i1, log2N, powa,
 #ifdef POWA_STRIDE  #ifdef POWA_STRIDE
                  a1,       a1,
 #endif  #endif
                  out, o1, osi, nout, Ninv, P, Pinv, wk )       out, o1, osi, nout, Ninv, P, Pinv, wk )
 ModNum in[], powa[], out[], wk[];  ModNum in[], powa[], out[], wk[];
 int N, log2N, osi, nout, i1, o1;  int N, log2N, osi, nout, i1, o1;
 #ifdef POWA_STRIDE  #ifdef POWA_STRIDE
Line 319  double Pinv;
Line 319  double Pinv;
   ModNum *qi, *qo, *qis, *qos, a, t0, tn, tt;    ModNum *qi, *qo, *qis, *qos, a, t0, tn, tt;
   
   /*    /*
    *    for k = K-1, ..., 0 do the following:     *  for k = K-1, ..., 0 do the following:
    *     *
    *        Out[2^{K-k}*i + j]             <= t0 + \alpha^{- 2^k * j} * tn;     *      Out[2^{K-k}*i + j]       <= t0 + \alpha^{- 2^k * j} * tn;
    *                                        = t0 - \alpha^{N/2 - 2^k * j}*tn;     *              = t0 - \alpha^{N/2 - 2^k * j}*tn;
    *        Out[2^{K-k}*i + j + 2^{K-k-1}] <= t0 - \alpha^{- 2^k * j} * tn;     *      Out[2^{K-k}*i + j + 2^{K-k-1}] <= t0 - \alpha^{- 2^k * j} * tn;
    *                                        = t0 + \alpha^{N/2 - 2^k * j}*tn;     *              = t0 + \alpha^{N/2 - 2^k * j}*tn;
    *     *
    *      where t0 = In[2^{K-k-1}*i + j] and     *      where t0 = In[2^{K-k-1}*i + j] and
    *            tn = In[2^{K-k-1}*i + j + 2^{K-1}],     *    tn = In[2^{K-k-1}*i + j + 2^{K-1}],
    *      for 0 <= i < 2^k and 0 <= j < 2^{K-k-1}.     *      for 0 <= i < 2^k and 0 <= j < 2^{K-k-1}.
    *     *
    */     */
   if ( log2N == 1 ) {    if ( log2N == 1 ) {
     /*  K = 1, k = 0, 0 <= i < 1, 0 <= j < 1.      /*  K = 1, k = 0, 0 <= i < 1, 0 <= j < 1.
      *  Out[0] <= t0 + tn, Out[1] <= t0 - tn,       *  Out[0] <= t0 + tn, Out[1] <= t0 - tn,
      *  where t0 = In[0] and Tn = In[1].       *  where t0 = In[0] and Tn = In[1].
      */       */
     t0 = in[0],  tn = in[i1];      t0 = in[0],  tn = in[i1];
Line 341  double Pinv;
Line 341  double Pinv;
       /*  out[0] = (t0 + tn)*Ninv; */        /*  out[0] = (t0 + tn)*Ninv; */
       tt = tn == (ModNum)0 ? t0 : AplusBmodP( t0, tn, P );        tt = tn == (ModNum)0 ? t0 : AplusBmodP( t0, tn, P );
       if ( tt != (ModNum)0 ) {        if ( tt != (ModNum)0 ) {
         AxBmodP( tt, ModNum, tt, Ninv, P, Pinv );    AxBmodP( tt, ModNum, tt, Ninv, P, Pinv );
       }        }
       out[0] = tt;        out[0] = tt;
       i = 1;        i = 1;
Line 351  double Pinv;
Line 351  double Pinv;
       /*  out[osi == 0 ? 1 : 0] = (t0 - tn)*Ninv; */        /*  out[osi == 0 ? 1 : 0] = (t0 - tn)*Ninv; */
       tt = tn == (ModNum)0 ? t0 : A_BmodP( t0, tn, P );        tt = tn == (ModNum)0 ? t0 : A_BmodP( t0, tn, P );
       if ( tt != (ModNum)0 ) {        if ( tt != (ModNum)0 ) {
         AxBmodP( tt, ModNum, tt, Ninv, P, Pinv );    AxBmodP( tt, ModNum, tt, Ninv, P, Pinv );
       }        }
       out[i] = tt;        out[i] = tt;
     }      }
     return;      return;
   }    }
   /****/    /****/
   halfN = n = 1 << (k = log2N - 1);     /* == 2^{K-1} */    halfN = n = 1 << (k = log2N - 1);  /* == 2^{K-1} */
 #ifdef POWA_STRIDE  #ifdef POWA_STRIDE
   halfNa = a1 << k;    halfNa = a1 << k;
 #endif  #endif
   /*    /*
    *  when k = K-1,     *  when k = K-1,
    *    Out[2*i] <= t0 + tn, and Out[2*i+1] <= t0 - tn,     *  Out[2*i] <= t0 + tn, and Out[2*i+1] <= t0 - tn,
    *  where t0 = In[i] and tn = In[i+2^{K-1}],     *  where t0 = In[i] and tn = In[i+2^{K-1}],
    *  for 0 <= i < 2^{K-1}.     *  for 0 <= i < 2^{K-1}.
    */     */
Line 390  double Pinv;
Line 390  double Pinv;
   idist = halfN;    idist = halfN;
   for ( ostep = 2; --k > 0; inwk = 1 - inwk ) {    for ( ostep = 2; --k > 0; inwk = 1 - inwk ) {
     n >>= 1;      n >>= 1;
     istep = ostep;      /* == 2^{K-k-1} */      istep = ostep;  /* == 2^{K-k-1} */
     ostep <<= 1;        /* == 2^{K-k} */      ostep <<= 1;  /* == 2^{K-k} */
     if ( inwk ) qi = &wk[N], qo = wk;      if ( inwk ) qi = &wk[N], qo = wk;
     else qi = wk, qo = &wk[N];      else qi = wk, qo = &wk[N];
     qis = qi,  qos = qo;      qis = qi,  qos = qo;
     /*      /*
      *  for j = 0,       *  for j = 0,
      *  Out[2^{K-k}*i] <= t0 + tn, and Out[2^{K-k}*i + 2^{K-k-1}] <= t0 - tn,       *  Out[2^{K-k}*i] <= t0 + tn, and Out[2^{K-k}*i + 2^{K-k-1}] <= t0 - tn,
      *  where t0 = In[2^{K-k-1}*i] and tn = In[2^{K-k-1}*i + 2^{K-1}].       *  where t0 = In[2^{K-k-1}*i] and tn = In[2^{K-k-1}*i + 2^{K-1}].
      */       */
     for ( i = n, is = os = 0; i-- > 0; qi += istep, qo += ostep ) {      for ( i = n, is = os = 0; i-- > 0; qi += istep, qo += ostep ) {
       t0 = qi[0],  tn = qi[idist];        t0 = qi[0],  tn = qi[idist];
       if ( tn == (ModNum)0 ) qo[0] = qo[istep] = t0;        if ( tn == (ModNum)0 ) qo[0] = qo[istep] = t0;
       else {        else {
         qo[0] = AplusBmodP( t0, tn, P );    qo[0] = AplusBmodP( t0, tn, P );
         qo[istep] = A_BmodP( t0, tn, P );    qo[istep] = A_BmodP( t0, tn, P );
       }        }
     }      }
 #ifdef POWA_STRIDE  #ifdef POWA_STRIDE
Line 419  double Pinv;
Line 419  double Pinv;
 #endif  #endif
       qi = ++qis,  qo = ++qos;        qi = ++qis,  qo = ++qos;
       for ( i = n; i-- > 0; qi += istep, qo += ostep ) {        for ( i = n; i-- > 0; qi += istep, qo += ostep ) {
         t0 = qi[0],  tn = qi[idist];    t0 = qi[0],  tn = qi[idist];
         if ( tn == (ModNum)0 ) qo[0] = qo[istep] = t0;    if ( tn == (ModNum)0 ) qo[0] = qo[istep] = t0;
         else {    else {
           AxBmodP( tn, ModNum, tn, a, P, Pinv );      AxBmodP( tn, ModNum, tn, a, P, Pinv );
 /***      qo[0] = AplusBmodP( t0, tn, P ); ***/  /***    qo[0] = AplusBmodP( t0, tn, P ); ***/
 /***      qo[istep] = A_BmodP( t0, tn, P ); ***/  /***    qo[istep] = A_BmodP( t0, tn, P ); ***/
           qo[0] = A_BmodP( t0, tn, P );      qo[0] = A_BmodP( t0, tn, P );
           qo[istep] = AplusBmodP( t0, tn, P );      qo[istep] = AplusBmodP( t0, tn, P );
         }    }
       }        }
     }      }
   }    }
Line 437  double Pinv;
Line 437  double Pinv;
   /*    /*
    *  The final step of k = 0.  i can take only the value 0.     *  The final step of k = 0.  i can take only the value 0.
    *     *
    *        Out[j]           <= t0 + \alpha^{-j} * tn;     *      Out[j]       <= t0 + \alpha^{-j} * tn;
    *        Out[j + 2^{K-1}] <= t0 - \alpha^{-j} * tn;     *      Out[j + 2^{K-1}] <= t0 - \alpha^{-j} * tn;
    *     *
    *      where t0 = In[j] and tn = In[j + 2^{K-1}],     *      where t0 = In[j] and tn = In[j + 2^{K-1}],
    *      for 0 <= j < 2^{K-1}.     *      for 0 <= j < 2^{K-1}.
    *     *
    */     */
   qo = out,  qi = &wk[inwk ? N : 0];    qo = out,  qi = &wk[inwk ? N : 0];
   if ( (n = osi + nout) > N ) nout = (n = N) - osi;     /* true nout */    if ( (n = osi + nout) > N ) nout = (n = N) - osi;  /* true nout */
   if ( (k = nout - halfN) <= 0 ) {      /* no overlap, i.e., no +- */    if ( (k = nout - halfN) <= 0 ) {  /* no overlap, i.e., no +- */
     if ( osi <= 0 ) {      if ( osi <= 0 ) {
       t0 = qi[0],  tn = qi[idist];        t0 = qi[0],  tn = qi[idist];
       if ( tn != (ModNum)0 ) t0 = AplusBmodP( t0, tn, P );        if ( tn != (ModNum)0 ) t0 = AplusBmodP( t0, tn, P );
       if ( t0 == (ModNum)0 ) qo[0] = t0;        if ( t0 == (ModNum)0 ) qo[0] = t0;
       else {        else {
         AxBmodP( qo[0], ModNum, t0, Ninv, P, Pinv );    AxBmodP( qo[0], ModNum, t0, Ninv, P, Pinv );
       }        }
 #ifdef POWA_STRIDE  #ifdef POWA_STRIDE
       aj = a1;        aj = a1;
Line 463  double Pinv;
Line 463  double Pinv;
 #ifdef POWA_STRIDE  #ifdef POWA_STRIDE
       aj = osi*a1;        aj = osi*a1;
 #endif  #endif
       if ( n <= halfN ) i = 0, k = n;           /* lower only */        if ( n <= halfN ) i = 0, k = n;    /* lower only */
       else if ( j == halfN ) goto L_halfN;      /* start with j = N/2 */        else if ( j == halfN ) goto L_halfN;  /* start with j = N/2 */
       else if ( j > halfN ) goto L_upper;       /* upper only */        else if ( j > halfN ) goto L_upper;  /* upper only */
       else i = n - halfN + 1, k = halfN;        /* we have lower, and i upper */        else i = n - halfN + 1, k = halfN;  /* we have lower, and i upper */
     }      }
   } else {    } else {
     os = o1*halfN;      os = o1*halfN;
Line 476  double Pinv;
Line 476  double Pinv;
     if ( osi <= 0 ) {      if ( osi <= 0 ) {
       t0 = qi[0],  tn = qi[idist];        t0 = qi[0],  tn = qi[idist];
       if ( tn == (ModNum)0 ) {        if ( tn == (ModNum)0 ) {
         if ( t0 == (ModNum)0 ) tt = t0;    if ( t0 == (ModNum)0 ) tt = t0;
         else { AxBmodP( tt, ModNum, t0, Ninv, P, Pinv ); }    else { AxBmodP( tt, ModNum, t0, Ninv, P, Pinv ); }
         qo[0] = qo[os] = tt;    qo[0] = qo[os] = tt;
       } else {        } else {
         tt = AplusBmodP( t0, tn, P );    tt = AplusBmodP( t0, tn, P );
         if ( tt == (ModNum)0 ) qo[0] = tt;    if ( tt == (ModNum)0 ) qo[0] = tt;
         else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }    else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
         tt = A_BmodP( t0, tn, P );    tt = A_BmodP( t0, tn, P );
         if ( tt == (ModNum)0 ) qo[os] = tt;    if ( tt == (ModNum)0 ) qo[os] = tt;
         else { AxBmodP( qo[os], ModNum, tt, Ninv, P, Pinv ); }    else { AxBmodP( qo[os], ModNum, tt, Ninv, P, Pinv ); }
       }        }
 #ifdef EBUG  #ifdef EBUG
     fprintf( stderr, "::::: DFT^{-1} ::: [0]=%d, [%d]=%d\n", (int)qo[0], os, (int)qo[os] );      fprintf( stderr, "::::: DFT^{-1} ::: [0]=%d, [%d]=%d\n", (int)qo[0], os, (int)qo[os] );
Line 503  double Pinv;
Line 503  double Pinv;
 #endif  #endif
       t0 = qi[j],  tn = qi[j+idist];        t0 = qi[j],  tn = qi[j+idist];
       if ( tn == (ModNum)0 ) {        if ( tn == (ModNum)0 ) {
         if ( t0 != (ModNum)0 ) {    if ( t0 != (ModNum)0 ) {
           AxBmodP( t0, ModNum, t0, Ninv, P, Pinv );      AxBmodP( t0, ModNum, t0, Ninv, P, Pinv );
         }    }
         qo[0] = qo[os] = t0;    qo[0] = qo[os] = t0;
       } else {        } else {
 #ifdef POWA_STRIDE  #ifdef POWA_STRIDE
 /***    a = P - powa[halfNa - aj]; ***/  /***  a = P - powa[halfNa - aj]; ***/
         a = powa[halfNa - aj];    a = powa[halfNa - aj];
 #else  #else
 /***    a = P - powa[halfN - j]; ***/  /***  a = P - powa[halfN - j]; ***/
         a = powa[halfN - j];    a = powa[halfN - j];
 #endif  #endif
         AxBmodP( tn, ModNum, tn, a, P, Pinv );    AxBmodP( tn, ModNum, tn, a, P, Pinv );
 /***    tt = AplusBmodP( t0, tn, P ); ***/  /***  tt = AplusBmodP( t0, tn, P ); ***/
         tt = A_BmodP( t0, tn, P );    tt = A_BmodP( t0, tn, P );
         if ( tt == (ModNum)0 ) qo[0] = tt;    if ( tt == (ModNum)0 ) qo[0] = tt;
         else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }    else { AxBmodP( qo[0], ModNum, tt, Ninv, P, Pinv ); }
 /***    tt = A_BmodP( t0, tn, P ); ***/  /***  tt = A_BmodP( t0, tn, P ); ***/
         tt = AplusBmodP( t0, tn, P );    tt = AplusBmodP( t0, tn, P );
         if ( tt == (ModNum)0 ) qo[os] = tt;    if ( tt == (ModNum)0 ) qo[os] = tt;
         else { AxBmodP( qo[os], ModNum, tt, Ninv, P, Pinv ); }    else { AxBmodP( qo[os], ModNum, tt, Ninv, P, Pinv ); }
       }        }
     }      }
     k = halfN;      k = halfN;
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>