nullspac.cpp
上传用户:jtjnyq9001
上传日期:2014-11-21
资源大小:3974k
文件大小:2k
源码类别:

3G开发

开发平台:

Visual C++

  1. //
  2. // nullspac.cpp
  3. //
  4. #include <iostream>
  5. #include "nullspac.h"
  6. #include <fstream>
  7. extern ofstream DebugFile;
  8. matrix_pf * NullSpace(BerlekampMatrix* b_mtx, int *r_ret)
  9. {
  10.   int deg = b_mtx->GetDegree();
  11.   int prime_char = b_mtx->PrimeBase();
  12.   matrix_pf *v = new matrix_pf(0,deg,0,deg);
  13.   int *c = new int[deg];
  14.   PrimeFieldElem b_val;
  15.   int i, j, k, k1, r, s;
  16.   bool dependence_found;
  17.   for(i=0; i<deg; i++) c[i]=-1;
  18.   r=0;
  19.   for( k=0; k<deg; k++)
  20.     {
  21.     dependence_found = false;
  22.     for( j=0; j<deg; j++)
  23.       {
  24.       if( (*b_mtx)[k][j] == 0 ) continue;
  25.         if( c[j] >= 0 ) continue;
  26.           //---------------------------------
  27.           //  multiply column by -1/B[k,j]
  28.           b_val = neg(recip( (*b_mtx)[k][j]));
  29.           for( k1=0; k1<deg; k1++)
  30.             (*b_mtx)[k1][j] *= b_val;
  31.           //---------------------------------------
  32.           // add B[k,i] times column j to column i
  33.           // for all i .ne. j
  34.           for( i=0; i<deg; i++)
  35.             {
  36.             if( (i==j) || ((*b_mtx)[k][i] == 0) ) continue;
  37.               b_val = (*b_mtx)[k][i];
  38.               for( k1=0; k1<deg; k1++)
  39.                 {
  40.                 (*b_mtx)[k1][i] += b_val * (*b_mtx)[k1][j];
  41.                 }
  42.             }
  43.           c[j] = k;
  44.           dependence_found = true;
  45.           break;
  46.       } // end of loop over j
  47.     if( !dependence_found )
  48.       {
  49.       for( j=0; j<deg; j++)
  50.         {
  51.         (*v)[r][j] = PrimeFieldElem(prime_char,0);
  52.         for( s=0; s<deg; s++)
  53.           {
  54.           if( c[s] == j ) (*v)[r][j] = (*b_mtx)[k][s];
  55.           }
  56.         if( j==k ) (*v)[r][j] = PrimeFieldElem(prime_char,1);
  57.         }
  58.       r++;
  59.       }
  60.     } // end of loop over k
  61.   *r_ret = r;
  62.   return(v);
  63. };