Changeset 6f62c3 in git


Ignore:
Timestamp:
Sep 26, 2007, 11:17:40 AM (16 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
Children:
8b388e1e07f2e848edfe629e622f372c9b16c0ea
Parents:
6c71caa8d1d8cd461d9ad0f2cece5d1fa1754f96
Message:
*hannes: chin.remainder/farey for gcd in char 0


git-svn-id: file:///usr/local/Singular/svn/trunk@10310 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
factory
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • factory/cf_algorithm.h

    r6c71ca r6f62c3  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: cf_algorithm.h,v 1.14 2006-05-16 13:43:18 Singular Exp $ */
     2/* $Id: cf_algorithm.h,v 1.15 2007-09-26 09:17:39 Singular Exp $ */
    33
    44#ifndef INCL_CF_ALGORITHM_H
    … …  
    5151
    5252void chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew );
     53
     54CanonicalForm Farey ( const CanonicalForm & f, const CanonicalForm & q );
    5355//}}}
    5456
    5557//{{{ function declarations from cf_factor.cc
    5658bool isPurePoly(const CanonicalForm & f);
     59
     60bool isPurePoly_m(const CanonicalForm & f);
    5761
    5862CFFList factorize ( const CanonicalForm & f, bool issqrfree = false );

  • factory/cf_chinese.cc

    r6c71ca r6f62c3  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: cf_chinese.cc,v 1.10 2005-02-08 10:28:46 Singular Exp $ */
     2/* $Id: cf_chinese.cc,v 1.11 2007-09-26 09:17:39 Singular Exp $ */
    33
    44//{{{ docu
    … …  
    1818
    1919#include "canonicalform.h"
     20#include "cf_iter.h"
     21
    2022
    2123//{{{ void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
    … …  
    7072    u = mod( v1, q2 );
    7173    d = mod( x2-u, q2 );
    72     if ( d.isZero() ) {
    73         xnew = v1;
    74         qnew = q1 * q2;
     74    if ( d.isZero() )
     75    {
     76        xnew = v1;
     77        qnew = q1 * q2;
    7578        DEBDECLEVEL( cerr, "chineseRemainder" );
    76         return;
     79        return;
    7780    }
    7881    (void)bextgcd( q1, q2, s, dummy );
    … …  
    120123    DEBOUTLN( cerr, "array size = " << n );
    121124
    122     while ( n != 1 ) {
    123         i = j = start;
    124         while ( i < start + n - 1 ) {
    125             // This is a little bit dangerous: X[i] and X[j] (and
    126             // Q[i] and Q[j]) may refer to the same object.  But
    127             // xnew and qnew in the above function are modified
    128             // at the very end of the function, so we do not
    129             // modify x1 and q1, resp., by accident.
    130             chineseRemainder( X[i], Q[i], X[i+1], Q[i+1], X[j], Q[j] );
    131             i += 2;
    132             j++;
    133         }
    134 
    135         if ( n & 1 ) {
    136             X[j] = X[i];
    137             Q[j] = Q[i];
    138         }
    139         // Maybe we would get some memory back at this point if
    140         // we would set X[j+1, ..., n] and Q[j+1, ..., n] to zero
    141         // at this point?
    142 
    143         n = ( n + 1) / 2;
     125    while ( n != 1 )
     126    {
     127        i = j = start;
     128        while ( i < start + n - 1 )
     129        {
     130            // This is a little bit dangerous: X[i] and X[j] (and
     131            // Q[i] and Q[j]) may refer to the same object.  But
     132            // xnew and qnew in the above function are modified
     133            // at the very end of the function, so we do not
     134            // modify x1 and q1, resp., by accident.
     135            chineseRemainder( X[i], Q[i], X[i+1], Q[i+1], X[j], Q[j] );
     136            i += 2;
     137            j++;
     138        }
     139
     140        if ( n & 1 )
     141        {
     142            X[j] = X[i];
     143            Q[j] = Q[i];
     144        }
     145        // Maybe we would get some memory back at this point if
     146        // we would set X[j+1, ..., n] and Q[j+1, ..., n] to zero
     147        // at this point?
     148
     149        n = ( n + 1) / 2;
    144150    }
    145151    xnew = X[start];
    … …  
    149155}
    150156//}}}
     157
     158CanonicalForm Farey_n (CanonicalForm N, const CanonicalForm P)
     159//"USAGE:  Farey_n (N,P); P, N number;
     160//RETURN:  a rational number a/b such that a/b=N mod P
     161//         and |a|,|b|<(P/2)^{1/2}
     162{
     163   //assume(P>0);
     164   if (N<0){N=N+P;}
     165   CanonicalForm A,B,C,D,E;
     166   E=P;
     167   B=1;
     168   while (N!=0)
     169   {
     170        if (2*N*N<P)
     171        {
     172           return(N/B);
     173        }
     174        D=E % N;
     175        C=A-(E-E % N)/N*B;
     176        E=N;
     177        N=D;
     178        A=B;
     179        B=C;
     180   }
     181   return(0);
     182}
     183
     184CanonicalForm Farey ( const CanonicalForm & f, const CanonicalForm & q )
     185{
     186    Variable x = f.mvar();
     187    CanonicalForm result = 0;
     188    CanonicalForm c;
     189    CFIterator i;
     190    for ( i = f; i.hasTerms(); i++ )
     191    {
     192        c = i.coeff();
     193        if ( c.inCoeffDomain())
     194        {
     195          result += power( x, i.exp() ) * Farey_n(c,q);
     196        }
     197        else
     198          result += power( x, i.exp() ) * Farey(c,q);
     199    }
     200    return result;
     201}
     202

  • factory/cf_defs.h

    r6c71ca r6f62c3  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: cf_defs.h,v 1.10 2005-05-03 09:35:33 Singular Exp $ */
     2/* $Id: cf_defs.h,v 1.11 2007-09-26 09:17:39 Singular Exp $ */
    33
    44#ifndef INCL_CF_DEFS_H
    … …  
    3939const int SW_USE_NTL_GCD_P=10;
    4040const int SW_USE_NTL_SORT=11;
     41const int SW_USE_CHINREM_GCD=12;
    4142//}}}
    4243

  • factory/cf_factor.cc

    r6c71ca r6f62c3  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: cf_factor.cc,v 1.33 2007-08-03 11:32:04 Singular Exp $ */
     2/* $Id: cf_factor.cc,v 1.34 2007-09-26 09:17:39 Singular Exp $ */
    33
    44//{{{ docu
    … …  
    8181
    8282#if 1
    83 #ifndef NOSTREAMIO
     83//#ifndef NOSTREAMIO
    8484void out_cf(char *s1,const CanonicalForm &f,char *s2)
    8585{
    … …  
    120120      printf("+%d",f.intval());
    121121    }
    122     else //printf("+...");
     122    else
     123    #ifdef NOSTREAMIO
     124      printf("+...");
     125    #else
    123126       std::cout << f;
     127    #endif
    124128    //if (f.inZ()) printf("(Z)");
    125129    //else if (f.inQ()) printf("(Q)");
    … …  
    161165  if ((f-t)!=0) { printf("problem:\n");out_cf("factor:",f," has problems\n");}
    162166}
     167//#endif
    163168#endif
    164 #endif
    165 
     169
     170bool isPurePoly_m(const CanonicalForm & f)
     171{
     172  if (f.inBaseDomain()) return true;
     173  if (f.level()<0) return false;
     174  for (CFIterator i=f;i.hasTerms();i++)
     175  {
     176    if (!isPurePoly_m(i.coeff())) return false;
     177  }
     178  return true;
     179}
    166180bool isPurePoly(const CanonicalForm & f)
    167181{
    … …  
    173187  return true;
    174188}
     189
    175190
    176191///////////////////////////////////////////////////////////////

  • factory/cf_gcd.cc

    r6c71ca r6f62c3  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: cf_gcd.cc,v 1.49 2007-04-26 08:22:48 Singular Exp $ */
     2/* $Id: cf_gcd.cc,v 1.50 2007-09-26 09:17:39 Singular Exp $ */
    33
    44#include <config.h>
    … …  
    2929static void cf_prepgcd( const CanonicalForm &, const CanonicalForm &, int &, int &, int & );
    3030
     31void out_cf(char *s1,const CanonicalForm &f,char *s2);
     32
     33CanonicalForm
     34chinrem_gcd ( const CanonicalForm & FF, const CanonicalForm & GG );
     35 
    3136bool
    3237gcd_test_one ( const CanonicalForm & f, const CanonicalForm & g, bool swap )
    … …  
    3843    delete sample;
    3944    CanonicalForm lcf, lcg;
    40     if ( swap ) {
     45    if ( swap )
     46    {
    4147        lcf = swapvar( LC( f ), Variable(1), f.mvar() );
    4248        lcg = swapvar( LC( g ), Variable(1), f.mvar() );
    4349    }
    44     else {
     50    else
     51    {
    4552        lcf = LC( f, Variable(1) );
    4653        lcg = LC( g, Variable(1) );
    … …  
    5562        return false;
    5663    CanonicalForm F, G;
    57     if ( swap ) {
     64    if ( swap )
     65    {
    5866        F=swapvar( f, Variable(1), f.mvar() );
    5967        G=swapvar( g, Variable(1), g.mvar() );
    6068    }
    61     else {
     69    else
     70    {
    6271        F = f;
    6372        G = g;
    … …  
    329338    if ( !( pi.isUnivariate() && pi1.isUnivariate() ) )
    330339    {
     340#if 0
     341        CanonicalForm newGCD(CanonicalForm A, CanonicalForm B);
     342        //out_cf("F:",f,"\n");
     343        //out_cf("G:",g,"\n");
     344        //out_cf("newGCD:",newGCD(f,g),"\n");
     345        return newGCD(f,g);
     346#endif
    331347        if ( gcd_test_one( pi1, pi, true ) )
    332             return C;
     348        {
     349          C=abs(C);
     350          //out_cf("GCD:",C,"\n");
     351          return C;
     352        }
    333353        bpure = false;
    334      
    335354    }
    336355    else
    … …  
    362381    }
    363382    if ( degree( pi1, v ) == 0 )
    364         return C;
     383    {
     384      C=abs(C);
     385      //out_cf("GCD:",C,"\n");
     386      return C;
     387    }
    365388    pi /= content( pi );
    366389    if ( bpure )
    367390        pi /= pi.lc();
    368     return C * pi;
     391    C=abs(C*pi);
     392    //out_cf("GCD:",C,"\n");
     393    return C;
    369394}
    370395
    … …  
    403428    else
    404429        bi = -1;
    405     while ( degree( pi1, v ) > 0 ) {
     430    while ( degree( pi1, v ) > 0 )
     431    {
    406432        pi2 = psr( pi, pi1, v );
    407433        pi2 = pi2 / bi;
    408434        pi = pi1; pi1 = pi2;
    409         if ( degree( pi1, v ) > 0 ) {
     435        if ( degree( pi1, v ) > 0 )
     436        {
    410437            delta = degree( pi, v ) - degree( pi1, v );
    411438            if ( (delta+1) % 2 )
    … …  
    474501        }
    475502    }
    476     else if ( isOn( SW_USE_EZGCD ) && !f.isUnivariate() )
    477     {
    478         if ( pe == 1 )
    479             fc = ezgcd( fc, gc );
    480         else if ( pe > 0 )// no variable at position 1
    481         {
    482             fc = replacevar( fc, Variable(pe), Variable(1) );
    483             gc = replacevar( gc, Variable(pe), Variable(1) );
    484             fc = replacevar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
    485         }
     503    else if (
     504    isOn(SW_USE_CHINREM_GCD)
     505    && (!fc.isUnivariate()) && (!gc.isUnivariate())
     506    && (isPurePoly_m(fc)) && (isPurePoly_m(gc))
     507    )
     508    {
     509        if ( p1 == fc.level() )
     510            fc = chinrem_gcd( fc, gc );
    486511        else
    487512        {
    488             pe = -pe;
    489             fc = swapvar( fc, Variable(pe), Variable(1) );
    490             gc = swapvar( gc, Variable(pe), Variable(1) );
    491             fc = swapvar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
    492         }
     513            fc = replacevar( fc, Variable(p1), Variable(mp) );
     514            gc = replacevar( gc, Variable(p1), Variable(mp) );
     515            fc = replacevar( chinrem_gcd( fc, gc ), Variable(mp), Variable(p1) );
     516        }
     517      //fc = chinrem_gcd( fc, gc);
     518    }
     519    else if ( isOn( SW_USE_EZGCD ) && !fc.isUnivariate() )
     520    {
     521      if ( pe == 1 )
     522          fc = ezgcd( fc, gc );
     523      else if ( pe > 0 )// no variable at position 1
     524      {
     525          fc = replacevar( fc, Variable(pe), Variable(1) );
     526          gc = replacevar( gc, Variable(pe), Variable(1) );
     527          fc = replacevar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
     528      }
     529      else
     530      {
     531          pe = -pe;
     532          fc = swapvar( fc, Variable(pe), Variable(1) );
     533          gc = swapvar( gc, Variable(pe), Variable(1) );
     534          fc = swapvar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
     535      }
    493536    }
    494537    else
    … …  
    515558cf_content ( const CanonicalForm & f, const CanonicalForm & g )
    516559{
    517     if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) ) {
     560    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
     561    {
    518562        CFIterator i = f;
    519563        CanonicalForm result = g;
    520         while ( i.hasTerms() && ! result.isOne() ) {
     564        while ( i.hasTerms() && ! result.isOne() )
     565        {
    521566            result = gcd( i.coeff(), result );
    522567            i++;
    … …  
    540585content ( const CanonicalForm & f )
    541586{
    542     if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) ) {
     587    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
     588    {
    543589        CFIterator i = f;
    544590        CanonicalForm result = abs( i.coeff() );
    545591        i++;
    546         while ( i.hasTerms() && ! result.isOne() ) {
     592        while ( i.hasTerms() && ! result.isOne() )
     593        {
    547594            result = gcd( i.coeff(), result );
    548595            i++;
    … …  
    910957    delete [] degsg;
    911958}
     959
     960
     961static CanonicalForm
     962balance_p ( const CanonicalForm & f, const CanonicalForm & q )
     963{
     964    Variable x = f.mvar();
     965    CanonicalForm result = 0, qh = q / 2;
     966    CanonicalForm c;
     967    CFIterator i;
     968    for ( i = f; i.hasTerms(); i++ )
     969    {
     970        c = i.coeff();
     971        if ( c.inCoeffDomain())
     972        {
     973          if ( c > qh )
     974            result += power( x, i.exp() ) * (c - q);
     975          else
     976            result += power( x, i.exp() ) * c;
     977        }
     978        else
     979          result += power( x, i.exp() ) * balance_p(c,q);
     980    }
     981    return result;
     982}
     983
     984#define GCD_CHINES_MIN_TRIES 3
     985#define GCD_CHINES_MAX_TRIES 7
     986CanonicalForm chinrem_gcd ( const CanonicalForm & FF, const CanonicalForm & GG )
     987{
     988  CanonicalForm f, g, cg, cl, q, Dp, newD, D, newq;
     989  int p, i, n, dp_deg, d_deg;;
     990
     991  CanonicalForm cd = bCommonDen( FF );
     992  f=cd*FF;
     993  f /=vcontent(f,Variable(1));
     994  cd = bCommonDen( f ); f *=cd;
     995  f /=vcontent(f,Variable(1));
     996
     997  cd = bCommonDen( GG );
     998  g=cd*GG;
     999  g /=vcontent(g,Variable(1));
     1000  cd = bCommonDen( g ); g *=cd;
     1001  g /=vcontent(g,Variable(1));
     1002
     1003  q = 0;
     1004  i = cf_getNumBigPrimes() - 1;
     1005  cl =  f.lc()* g.lc();
     1006
     1007  n=GCD_CHINES_MIN_TRIES;
     1008  while ( true )
     1009  {
     1010    p = cf_getBigPrime( i );
     1011    i--;
     1012    while ( i >= 0 && mod( cl, p ) == 0 )
     1013    {
     1014      p = cf_getBigPrime( i );
     1015      i--;
     1016    }
     1017    setCharacteristic( p );
     1018    n--;
     1019    Dp = gcd( mapinto( f ), mapinto( g ) );
     1020    setCharacteristic( 0 );
     1021    dp_deg=totaldegree(Dp);
     1022    if ( dp_deg == 0 )
     1023      return CanonicalForm(1);
     1024    if ( q.isZero() )
     1025    {
     1026      D = mapinto( Dp );
     1027      d_deg=dp_deg;
     1028      q = p;
     1029    }
     1030    else
     1031    {
     1032      if ( dp_deg == d_deg )
     1033      {
     1034        chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
     1035        q = newq;
     1036        D = newD;
     1037      }
     1038      else if ( dp_deg < d_deg )
     1039      {
     1040        n=GCD_CHINES_MIN_TRIES;
     1041        // all previous p's are bad primes
     1042        q = p;
     1043        D = mapinto( Dp );
     1044        d_deg=dp_deg;
     1045      }
     1046      //else
     1047      //  printf("bad prime\n");
     1048      // else p is a bad prime
     1049    }
     1050    if (( i >= 0 ) && (n >= GCD_CHINES_MIN_TRIES-GCD_CHINES_MAX_TRIES))
     1051    {
     1052      if (n<=0)
     1053      {
     1054      // now balance D mod q
     1055      //CanonicalForm Dn = balance_p( D, q );
     1056      //out_cf("test ",Dn,"\n");
     1057      //if ( fdivides( Dn, f ) && fdivides( Dn, g ) )
     1058      //{
     1059      //  //printf("does divide\n");
     1060      //  return Dn;
     1061      //}
     1062      CanonicalForm Dn= Farey(D,q);
     1063      CanonicalForm cd = bCommonDen( Dn );
     1064      Dn *=cd;
     1065      Dn /=vcontent(Dn,Variable(1));
     1066      if ( fdivides( Dn, f ) && fdivides( Dn, g ) )
     1067      {
     1068        //printf("does divide\n");
     1069        return Dn;
     1070      }
     1071      //else
     1072      // printf("does not divide\n");
     1073      }
     1074    }
     1075    else
     1076    {
     1077      #if 0
     1078      printf("primes left:%d, tries left:%d, gcd-deg:%d\n",i,n,d_deg);
     1079      out_cf("f= ",f,"\n");
     1080      out_cf("g= ",g,"\n");
     1081      CanonicalForm Dn= Farey(D,q);
     1082      CanonicalForm cd = bCommonDen( Dn );
     1083      Dn *=cd;
     1084      Dn /=vcontent(Dn,Variable(1));
     1085      out_cf("farey: ",Dn,"\n");
     1086      #endif
     1087      Off(SW_USE_CHINREM_GCD);
     1088      D=gcd_poly( f, g );
     1089      On(SW_USE_CHINREM_GCD);
     1090      #if 0
     1091      out_cf("gcd: ",D,"\n");
     1092      #endif
     1093      return D;
     1094    }
     1095  }
     1096}
     1097

  • factory/cf_linsys.cc

    r6c71ca r6f62c3  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: cf_linsys.cc,v 1.10 1997-08-29 07:26:40 schmidt Exp $ */
     2/* $Id: cf_linsys.cc,v 1.11 2007-09-26 09:17:40 Singular Exp $ */
    33
    44#include <config.h>
    … …  
    3535    int i, j;
    3636    for ( i = 1; i <= rows; i++ )
    37         for ( j = 1; j <= rows; j++ )
    38             if ( ! M(i,j).inZ() )
    39                 return false;
     37        for ( j = 1; j <= rows; j++ )
     38            if ( ! M(i,j).inZ() )
     39                return false;
    4040    return true;
    4141}
    … …  
    4646    int i, j, rows = M.rows(), cols = M.columns();
    4747    for ( i = 1; i <= rows; i++ )
    48         for ( j = 1; j <= cols; j++ )
    49             if ( ! M(i,j).inZ() )
    50                 return false;
     48        for ( j = 1; j <= cols; j++ )
     49            if ( ! M(i,j).inZ() )
     50                return false;
    5151    return true;
    5252}
    … …  
    5656{
    5757    if ( newpivot.isZero() )
    58         return false;
     58        return false;
    5959    else  if ( oldpivot.isZero() )
    60         return true;
     60        return true;
    6161    else  if ( level( oldpivot ) > level( newpivot ) )
    62         return true;
     62        return true;
    6363    else  if ( level( oldpivot ) < level( newpivot ) )
    64         return false;
     64        return false;
    6565    else
    66         return ( newpivot.lc() < oldpivot.lc() );
     66        return ( newpivot.lc() < oldpivot.lc() );
    6767}
    6868
    … …  
    7575
    7676    if ( ! matrix_in_Z( M ) ) {
    77         int nrows = M.rows(), ncols = M.columns();
    78         int i, j, k;
    79         CanonicalForm rowpivot, pivotrecip;
    80         // triangularization
    81         for ( i = 1; i <= nrows; i++ ) {
    82             //find "pivot"
    83             for (j = i; j <= nrows; j++ )
    84                 if ( M(j,i) != 0 ) break;
    85             if ( j > nrows ) return false;
    86             if ( j != i )
    87                 M.swapRow( i, j );
    88             pivotrecip = 1 / M(i,i);
    89             for ( j = 1; j <= ncols; j++ )
    90                 M(i,j) *= pivotrecip;
    91             for ( j = i+1; j <= nrows; j++ ) {
    92                 rowpivot = M(j,i);
    93                 if ( rowpivot == 0 ) continue;
    94                 for ( k = i; k <= ncols; k++ )
    95                     M(j,k) -= M(i,k) * rowpivot;
    96             }
    97         }
    98         // matrix is now upper triangular with 1s down the diagonal
    99         // back-substitute
    100         for ( i = nrows-1; i > 0; i-- ) {
    101             for ( j = nrows+1; j <= ncols; j++ ) {
    102                 for ( k = i+1; k <= nrows; k++ )
    103                     M(i,j) -= M(k,j) * M(i,k);
    104             }
    105         }
    106         return true;
     77        int nrows = M.rows(), ncols = M.columns();
     78        int i, j, k;
     79        CanonicalForm rowpivot, pivotrecip;
     80        // triangularization
     81        for ( i = 1; i <= nrows; i++ ) {
     82            //find "pivot"
     83            for (j = i; j <= nrows; j++ )
     84                if ( M(j,i) != 0 ) break;
     85            if ( j > nrows ) return false;
     86            if ( j != i )
     87                M.swapRow( i, j );
     88            pivotrecip = 1 / M(i,i);
     89            for ( j = 1; j <= ncols; j++ )
     90                M(i,j) *= pivotrecip;
     91            for ( j = i+1; j <= nrows; j++ ) {
     92                rowpivot = M(j,i);
     93                if ( rowpivot == 0 ) continue;
     94                for ( k = i; k <= ncols; k++ )
     95                    M(j,k) -= M(i,k) * rowpivot;
     96            }
     97        }
     98        // matrix is now upper triangular with 1s down the diagonal
     99        // back-substitute
     100        for ( i = nrows-1; i > 0; i-- ) {
     101            for ( j = nrows+1; j <= ncols; j++ ) {
     102                for ( k = i+1; k <= nrows; k++ )
     103                    M(i,j) -= M(k,j) * M(i,k);
     104            }
     105        }
     106        return true;
    107107    }
    108108    else {
    109         int rows = M.rows(), cols = M.columns();
    110         CFMatrix MM( rows, cols );
    111         int ** mm = new int_ptr[rows];
    112         CanonicalForm Q, Qhalf, mnew, qnew, B;
    113         int i, j, p, pno;
    114         bool ok;
    115 
    116         // initialize room to hold the result and the result mod p
    117         for ( i = 0; i < rows; i++ ) {
    118             mm[i] = new int[cols];
    119         }
    120 
    121         // calculate the bound for the result
    122         B = bound( M );
    123         DEBOUTLN( cerr, "bound = " <<  B );
    124 
    125         // find a first solution mod p
    126         pno = 0;
    127         do {
    128             DEBOUTSL( cerr );
    129             DEBOUT( cerr, "trying prime(" << pno << ") = " );
    130             p = cf_getBigPrime( pno );
    131             DEBOUT( cerr, p );
    132             DEBOUTENDL( cerr );
    133             setCharacteristic( p );
    134             // map matrix into char p
    135             for ( i = 1; i <= rows; i++ )
    136                 for ( j = 1; j <= cols; j++ )
    137                     mm[i-1][j-1] = mapinto( M(i,j) ).intval();
    138             // solve mod p
    139             ok = solve( mm, rows, cols );
    140             pno++;
    141         } while ( ! ok );
    142 
    143         // initialize the result matrix with first solution
    144         setCharacteristic( 0 );
    145         for ( i = 1; i <= rows; i++ )
    146             for ( j = rows+1; j <= cols; j++ )
    147                 MM(i,j) = mm[i-1][j-1];
    148 
    149         // Q so far
    150         Q = p;
    151         while ( Q < B && pno < cf_getNumBigPrimes() ) {
    152             do {
    153                 DEBOUTSL( cerr );
    154                 DEBOUT( cerr, "trying prime(" << pno << ") = " );
    155                 p = cf_getBigPrime( pno );
    156                 DEBOUT( cerr, p );
    157                 DEBOUTENDL( cerr );
    158                 setCharacteristic( p );
    159                 for ( i = 1; i <= rows; i++ )
    160                     for ( j = 1; j <= cols; j++ )
    161                         mm[i-1][j-1] = mapinto( M(i,j) ).intval();
    162                 // solve mod p
    163                 ok = solve( mm, rows, cols );
    164                 pno++;
    165             } while ( ! ok );
    166             // found a solution mod p
    167             // now chinese remainder it to a solution mod Q*p
    168             setCharacteristic( 0 );
    169             for ( i = 1; i <= rows; i++ )
    170                 for ( j = rows+1; j <= cols; j++ ) {
    171                     chineseRemainder( MM[i][j], Q, CanonicalForm(mm[i-1][j-1]), CanonicalForm(p), mnew, qnew );
    172                     MM(i, j) = mnew;
    173                 }
    174             Q = qnew;
    175         }
    176         if ( pno == cf_getNumBigPrimes() )
    177             fuzzy_result = true;
    178         else
    179             fuzzy_result = false;
    180         // store the result in M
    181         Qhalf = Q / 2;
    182         for ( i = 1; i <= rows; i++ ) {
    183             for ( j = rows+1; j <= cols; j++ )
    184                 if ( MM(i,j) > Qhalf )
    185                     M(i,j) = MM(i,j) - Q;
    186                 else
    187                     M(i,j) = MM(i,j);
    188             delete [] mm[i-1];
    189         }
    190         delete [] mm;
    191         return ! fuzzy_result;
     109        int rows = M.rows(), cols = M.columns();
     110        CFMatrix MM( rows, cols );
     111        int ** mm = new int_ptr[rows];
     112        CanonicalForm Q, Qhalf, mnew, qnew, B;
     113        int i, j, p, pno;
     114        bool ok;
     115
     116        // initialize room to hold the result and the result mod p
     117        for ( i = 0; i < rows; i++ ) {
     118            mm[i] = new int[cols];
     119        }
     120
     121        // calculate the bound for the result
     122        B = bound( M );
     123        DEBOUTLN( cerr, "bound = " <<  B );
     124
     125        // find a first solution mod p
     126        pno = 0;
     127        do {
     128            DEBOUTSL( cerr );
     129            DEBOUT( cerr, "trying prime(" << pno << ") = " );
     130            p = cf_getBigPrime( pno );
     131            DEBOUT( cerr, p );
     132            DEBOUTENDL( cerr );
     133            setCharacteristic( p );
     134            // map matrix into char p
     135            for ( i = 1; i <= rows; i++ )
     136                for ( j = 1; j <= cols; j++ )
     137                    mm[i-1][j-1] = mapinto( M(i,j) ).intval();
     138            // solve mod p
     139            ok = solve( mm, rows, cols );
     140            pno++;
     141        } while ( ! ok );
     142
     143        // initialize the result matrix with first solution
     144        setCharacteristic( 0 );
     145        for ( i = 1; i <= rows; i++ )
     146            for ( j = rows+1; j <= cols; j++ )
     147                MM(i,j) = mm[i-1][j-1];
     148
     149        // Q so far
     150        Q = p;
     151        while ( Q < B && pno < cf_getNumBigPrimes() ) {
     152            do {
     153                DEBOUTSL( cerr );
     154                DEBOUT( cerr, "trying prime(" << pno << ") = " );
     155                p = cf_getBigPrime( pno );
     156                DEBOUT( cerr, p );
     157                DEBOUTENDL( cerr );
     158                setCharacteristic( p );
     159                for ( i = 1; i <= rows; i++ )
     160                    for ( j = 1; j <= cols; j++ )
     161                        mm[i-1][j-1] = mapinto( M(i,j) ).intval();
     162                // solve mod p
     163                ok = solve( mm, rows, cols );
     164                pno++;
     165            } while ( ! ok );
     166            // found a solution mod p
     167            // now chinese remainder it to a solution mod Q*p
     168            setCharacteristic( 0 );
     169            for ( i = 1; i <= rows; i++ )
     170                for ( j = rows+1; j <= cols; j++ )
     171                {
     172                    chineseRemainder( MM[i][j], Q, CanonicalForm(mm[i-1][j-1]), CanonicalForm(p), mnew, qnew );
     173                    MM(i, j) = mnew;
     174                }
     175            Q = qnew;
     176        }
     177        if ( pno == cf_getNumBigPrimes() )
     178            fuzzy_result = true;
     179        else
     180            fuzzy_result = false;
     181        // store the result in M
     182        Qhalf = Q / 2;
     183        for ( i = 1; i <= rows; i++ ) {
     184            for ( j = rows+1; j <= cols; j++ )
     185                if ( MM(i,j) > Qhalf )
     186                    M(i,j) = MM(i,j) - Q;
     187                else
     188                    M(i,j) = MM(i,j);
     189            delete [] mm[i-1];
     190        }
     191        delete [] mm;
     192        return ! fuzzy_result;
    192193    }
    193194}
    … …  
    199200    bool ok = true;
    200201    for ( i = 0; i < rows && ok; i++ )
    201         for ( j = 0; j < rows && ok; j++ ) {
    202             if ( M(i+1,j+1).isZero() )
    203                 m[i][j] = 0;
    204             else {
    205                 m[i][j] = mapinto( M(i+1,j+1) ).intval();
    206 //              ok = m[i][j] != 0;
    207             }
    208         }
     202        for ( j = 0; j < rows && ok; j++ )
     203        {
     204            if ( M(i+1,j+1).isZero() )
     205                m[i][j] = 0;
     206            else
     207            {
     208                m[i][j] = mapinto( M(i+1,j+1) ).intval();
     209//                ok = m[i][j] != 0;
     210            }
     211        }
    209212    return ok;
    210213}
    … …  
    217220    ASSERT( rows <= M.rows() && rows <= M.columns() && rows > 0, "undefined determinant" );
    218221    if ( rows == 1 )
    219         return M(1,1);
     222        return M(1,1);
    220223    else  if ( rows == 2 )
    221         return M(1,1)*M(2,2)-M(2,1)*M(1,2);
    222     else  if ( matrix_in_Z( M, rows ) ) {
    223         int ** mm = new int_ptr[rows];
    224         CanonicalForm x, q, Qhalf, B;
    225         int n, i, intdet, p, pno;
    226         for ( i = 0; i < rows; i++ ) {
    227             mm[i] = new int[rows];
    228         }
    229         pno = 0; n = 0;
    230         TIMING_START(det_bound);
    231         B = detbound( M, rows );
    232         TIMING_END(det_bound);
    233         q = 1;
    234         TIMING_START(det_numprimes);
    235         while ( B > q && n < cf_getNumBigPrimes() ) {
    236             q *= cf_getBigPrime( n );
    237             n++;
    238         }
    239         TIMING_END(det_numprimes);
    240 
    241         CFArray X(1,n), Q(1,n);
    242 
    243         while ( pno < n ) {
    244             p = cf_getBigPrime( pno );
    245             setCharacteristic( p );
    246             // map matrix into char p
    247             TIMING_START(det_mapping);
    248             fill_int_mat( M, mm, rows );
    249             TIMING_END(det_mapping);
    250             pno++;
    251             DEBOUT( cerr, "." );
    252             TIMING_START(det_determinant);
    253             intdet = determinant( mm, rows );
    254             TIMING_END(det_determinant);
    255             setCharacteristic( 0 );
    256             X[pno] = intdet;
    257             Q[pno] = p;
    258         }
    259         TIMING_START(det_chinese);
    260         chineseRemainder( X, Q, x, q );
    261         TIMING_END(det_chinese);
    262         Qhalf = q / 2;
    263         if ( x > Qhalf )
    264             x = x - q;
    265         for ( i = 0; i < rows; i++ )
    266             delete [] mm[i];
    267         delete [] mm;
    268         return x;
    269     }
    270     else {
    271         CFMatrix m( M );
    272         CanonicalForm divisor = 1, pivot, mji;
    273         int i, j, k, sign = 1;
    274         for ( i = 1; i <= rows; i++ ) {
    275             pivot = m(i,i); k = i;
    276             for ( j = i+1; j <= rows; j++ ) {
    277                 if ( betterpivot( pivot, m(j,i) ) ) {
    278                     pivot = m(j,i);
    279                     k = j;
    280                 }
    281             }
    282             if ( pivot.isZero() )
    283                 return 0;
    284             if ( i != k ) {
    285                 m.swapRow( i, k );
    286                 sign = -sign;
    287             }
    288             for ( j = i+1; j <= rows; j++ ) {
    289                 if ( ! m(j,i).isZero() ) {
    290                     divisor *= pivot;
    291                     mji = m(j,i);
    292                     m(j,i) = 0;
    293                     for ( k = i+1; k <= rows; k++ )
    294                         m(j,k) = m(j,k) * pivot - m(i,k)*mji;
    295                 }
    296             }
    297         }
    298         pivot = sign;
    299         for ( i = 1; i <= rows; i++ )
    300             pivot *= m(i,i);
    301         return pivot / divisor;
     224        return M(1,1)*M(2,2)-M(2,1)*M(1,2);
     225    else  if ( matrix_in_Z( M, rows ) )
     226    {
     227        int ** mm = new int_ptr[rows];
     228        CanonicalForm x, q, Qhalf, B;
     229        int n, i, intdet, p, pno;
     230        for ( i = 0; i < rows; i++ )
     231        {
     232            mm[i] = new int[rows];
     233        }
     234        pno = 0; n = 0;
     235        TIMING_START(det_bound);
     236        B = detbound( M, rows );
     237        TIMING_END(det_bound);
     238        q = 1;
     239        TIMING_START(det_numprimes);
     240        while ( B > q && n < cf_getNumBigPrimes() )
     241        {
     242            q *= cf_getBigPrime( n );
     243            n++;
     244        }
     245        TIMING_END(det_numprimes);
     246
     247        CFArray X(1,n), Q(1,n);
     248
     249        while ( pno < n )
     250        {
     251            p = cf_getBigPrime( pno );
     252            setCharacteristic( p );
     253            // map matrix into char p
     254            TIMING_START(det_mapping);
     255            fill_int_mat( M, mm, rows );
     256            TIMING_END(det_mapping);
     257            pno++;
     258            DEBOUT( cerr, "." );
     259            TIMING_START(det_determinant);
     260            intdet = determinant( mm, rows );
     261            TIMING_END(det_determinant);
     262            setCharacteristic( 0 );
     263            X[pno] = intdet;
     264            Q[pno] = p;
     265        }
     266        TIMING_START(det_chinese);
     267        chineseRemainder( X, Q, x, q );
     268        TIMING_END(det_chinese);
     269        Qhalf = q / 2;
     270        if ( x > Qhalf )
     271            x = x - q;
     272        for ( i = 0; i < rows; i++ )
     273            delete [] mm[i];
     274        delete [] mm;
     275        return x;
     276    }
     277    else
     278    {
     279        CFMatrix m( M );
     280        CanonicalForm divisor = 1, pivot, mji;
     281        int i, j, k, sign = 1;
     282        for ( i = 1; i <= rows; i++ ) {
     283            pivot = m(i,i); k = i;
     284            for ( j = i+1; j <= rows; j++ ) {
     285                if ( betterpivot( pivot, m(j,i) ) ) {
     286                    pivot = m(j,i);
     287                    k = j;
     288                }
     289            }
     290            if ( pivot.isZero() )
     291                return 0;
     292            if ( i != k )
     293            {
     294                m.swapRow( i, k );
     295                sign = -sign;
     296            }
     297            for ( j = i+1; j <= rows; j++ )
     298            {
     299                if ( ! m(j,i).isZero() )
     300                {
     301                    divisor *= pivot;
     302                    mji = m(j,i);
     303                    m(j,i) = 0;
     304                    for ( k = i+1; k <= rows; k++ )
     305                        m(j,k) = m(j,k) * pivot - m(i,k)*mji;
     306                }
     307            }
     308        }
     309        pivot = sign;
     310        for ( i = 1; i <= rows; i++ )
     311            pivot *= m(i,i);
     312        return pivot / divisor;
    302313    }
    303314}
    … …  
    310321    ASSERT( rows <= M.rows() && rows <= M.columns() && rows > 0, "undefined determinant" );
    311322    if ( rows == 1 )
    312         return M(1,1);
     323        return M(1,1);
    313324    else  if ( rows == 2 )
    314         return M(1,1)*M(2,2)-M(2,1)*M(1,2);
     325        return M(1,1)*M(2,2)-M(2,1)*M(1,2);
    315326    else  if ( matrix_in_Z( M, rows ) ) {
    316         int ** mm = new int_ptr[rows];
    317         CanonicalForm QQ, Q, Qhalf, mnew, q, qnew, B;
    318         CanonicalForm det, detnew, qdet;
    319         int i, p, pcount, pno, intdet;
    320         bool ok;
    321 
    322         // initialize room to hold the result and the result mod p
    323         for ( i = 0; i < rows; i++ ) {
    324             mm[i] = new int[rows];
    325         }
    326 
    327         // calculate the bound for the result
    328         B = detbound( M, rows );
    329 
    330         // find a first solution mod p
    331         pno = 0;
    332         do {
    333             p = cf_getBigPrime( pno );
    334             setCharacteristic( p );
    335             // map matrix into char p
    336             ok = fill_int_mat( M, mm, rows );
    337             pno++;
    338         } while ( ! ok && pno < cf_getNumPrimes() );
    339         // initialize the result matrix with first solution
    340         // solve mod p
    341         DEBOUT( cerr, "." );
    342         intdet = determinant( mm, rows );
    343         setCharacteristic( 0 );
    344         det = intdet;
    345         // Q so far
    346         Q = p;
    347         QQ = p;
    348         while ( Q < B && cf_getNumPrimes() > pno ) {
    349             // find a first solution mod p
    350             do {
    351                 p = cf_getBigPrime( pno );
    352                 setCharacteristic( p );
    353                 // map matrix into char p
    354                 ok = fill_int_mat( M, mm, rows );
    355                 pno++;
    356             } while ( ! ok && pno < cf_getNumPrimes() );
    357             // initialize the result matrix with first solution
    358             // solve mod p
    359             DEBOUT( cerr, "." );
    360             intdet = determinant( mm, rows );
    361             setCharacteristic( 0 );
    362             qdet = intdet;
    363             // Q so far
    364             q = p;
    365             QQ *= p;
    366             pcount = 0;
    367             while ( QQ < B && cf_getNumPrimes() > pno && pcount < 500 ) {
    368                 do {
    369                     p = cf_getBigPrime( pno );
    370                     setCharacteristic( p );
    371                     ok = true;
    372                     // map matrix into char p
    373                     ok = fill_int_mat( M, mm, rows );
    374                     pno++;
    375                 } while ( ! ok && cf_getNumPrimes() > pno );
    376                 // solve mod p
    377                 DEBOUT( cerr, "." );
    378                 intdet = determinant( mm, rows );
    379                 // found a solution mod p
    380                 // now chinese remainder it to a solution mod Q*p
    381                 setCharacteristic( 0 );
    382                 chineseRemainder( qdet, q, intdet, p, detnew, qnew );
    383                 qdet = detnew;
    384                 q = qnew;
    385                 QQ *= p;
    386                 pcount++;
    387             }
    388             DEBOUT( cerr, "*" );
    389             chineseRemainder( det, Q, qdet, q, detnew, qnew );
    390             Q = qnew;
    391             QQ = Q;
    392             det = detnew;
    393         }
    394         if ( ! ok )
    395             fuzzy_result = true;
    396         else
    397             fuzzy_result = false;
    398         // store the result in M
    399         Qhalf = Q / 2;
    400         if ( det > Qhalf )
    401             det = det - Q;
    402         for ( i = 0; i < rows; i++ )
    403             delete [] mm[i];
    404         delete [] mm;
    405         return det;
     327        int ** mm = new int_ptr[rows];
     328        CanonicalForm QQ, Q, Qhalf, mnew, q, qnew, B;
     329        CanonicalForm det, detnew, qdet;
     330        int i, p, pcount, pno, intdet;
     331        bool ok;
     332
     333        // initialize room to hold the result and the result mod p
     334        for ( i = 0; i < rows; i++ ) {
     335            mm[i] = new int[rows];
     336        }
     337
     338        // calculate the bound for the result
     339        B = detbound( M, rows );
     340
     341        // find a first solution mod p
     342        pno = 0;
     343        do {
     344            p = cf_getBigPrime( pno );
     345            setCharacteristic( p );
     346            // map matrix into char p
     347            ok = fill_int_mat( M, mm, rows );
     348            pno++;
     349        } while ( ! ok && pno < cf_getNumPrimes() );
     350        // initialize the result matrix with first solution
     351        // solve mod p
     352        DEBOUT( cerr, "." );
     353        intdet = determinant( mm, rows );
     354        setCharacteristic( 0 );
     355        det = intdet;
     356        // Q so far
     357        Q = p;
     358        QQ = p;
     359        while ( Q < B && cf_getNumPrimes() > pno ) {
     360            // find a first solution mod p
     361            do {
     362                p = cf_getBigPrime( pno );
     363                setCharacteristic( p );
     364                // map matrix into char p
     365                ok = fill_int_mat( M, mm, rows );
     366                pno++;
     367            } while ( ! ok && pno < cf_getNumPrimes() );
     368            // initialize the result matrix with first solution
     369            // solve mod p
     370            DEBOUT( cerr, "." );
     371            intdet = determinant( mm, rows );
     372            setCharacteristic( 0 );
     373            qdet = intdet;
     374            // Q so far
     375            q = p;
     376            QQ *= p;
     377            pcount = 0;
     378            while ( QQ < B && cf_getNumPrimes() > pno && pcount < 500 ) {
     379                do {
     380                    p = cf_getBigPrime( pno );
     381                    setCharacteristic( p );
     382                    ok = true;
     383                    // map matrix into char p
     384                    ok = fill_int_mat( M, mm, rows );
     385                    pno++;
     386                } while ( ! ok && cf_getNumPrimes() > pno );
     387                // solve mod p
     388                DEBOUT( cerr, "." );
     389                intdet = determinant( mm, rows );
     390                // found a solution mod p
     391                // now chinese remainder it to a solution mod Q*p
     392                setCharacteristic( 0 );
     393                chineseRemainder( qdet, q, intdet, p, detnew, qnew );
     394                qdet = detnew;
     395                q = qnew;
     396                QQ *= p;
     397                pcount++;
     398            }
     399            DEBOUT( cerr, "*" );
     400            chineseRemainder( det, Q, qdet, q, detnew, qnew );
     401            Q = qnew;
     402            QQ = Q;
     403            det = detnew;
     404        }
     405        if ( ! ok )
     406            fuzzy_result = true;
     407        else
     408            fuzzy_result = false;
     409        // store the result in M
     410        Qhalf = Q / 2;
     411        if ( det > Qhalf )
     412            det = det - Q;
     413        for ( i = 0; i < rows; i++ )
     414            delete [] mm[i];
     415        delete [] mm;
     416        return det;
    406417    }
    407418    else {
    408         CFMatrix m( M );
    409         CanonicalForm divisor = 1, pivot, mji;
    410         int i, j, k, sign = 1;
    411         for ( i = 1; i <= rows; i++ ) {
    412             pivot = m(i,i); k = i;
    413             for ( j = i+1; j <= rows; j++ ) {
    414                 if ( betterpivot( pivot, m(j,i) ) ) {
    415                     pivot = m(j,i);
    416                     k = j;
    417                 }
    418             }
    419             if ( pivot.isZero() )
    420                 return 0;
    421             if ( i != k ) {
    422                 m.swapRow( i, k );
    423                 sign = -sign;
    424             }
    425             for ( j = i+1; j <= rows; j++ ) {
    426                 if ( ! m(j,i).isZero() ) {
    427                     divisor *= pivot;
    428                     mji = m(j,i);
    429                     m(j,i) = 0;
    430                     for ( k = i+1; k <= rows; k++ )
    431                         m(j,k) = m(j,k) * pivot - m(i,k)*mji;
    432                 }
    433             }
    434         }
    435         pivot = sign;
    436         for ( i = 1; i <= rows; i++ )
    437             pivot *= m(i,i);
    438         return pivot / divisor;
     419        CFMatrix m( M );
     420        CanonicalForm divisor = 1, pivot, mji;
     421        int i, j, k, sign = 1;
     422        for ( i = 1; i <= rows; i++ ) {
     423            pivot = m(i,i); k = i;
     424            for ( j = i+1; j <= rows; j++ ) {
     425                if ( betterpivot( pivot, m(j,i) ) ) {
     426                    pivot = m(j,i);
     427                    k = j;
     428                }
     429            }
     430            if ( pivot.isZero() )
     431                return 0;
     432            if ( i != k ) {
     433                m.swapRow( i, k );
     434                sign = -sign;
     435            }
     436            for ( j = i+1; j <= rows; j++ ) {
     437                if ( ! m(j,i).isZero() ) {
     438                    divisor *= pivot;
     439                    mji = m(j,i);
     440                    m(j,i) = 0;
     441                    for ( k = i+1; k <= rows; k++ )
     442                        m(j,k) = m(j,k) * pivot - m(i,k)*mji;
     443                }
     444            }
     445        }
     446        pivot = sign;
     447        for ( i = 1; i <= rows; i++ )
     448            pivot *= m(i,i);
     449        return pivot / divisor;
    439450    }
    440451}
    … …  
    448459    int i, j;
    449460    for ( i = 1; i <= rows; i++ )
    450         for ( j = 1; j <= rows; j++ )
    451             sum += M(i,j) * M(i,j);
     461        for ( j = 1; j <= rows; j++ )
     462            sum += M(i,j) * M(i,j);
    452463    DEBOUTLN( cerr, "bound(matrix)^2 = " << sum );
    453464    CanonicalForm vmax = 0, vsum;
    454465    for ( j = rows+1; j <= cols; j++ ) {
    455         vsum = 0;
    456         for ( i = 1; i <= rows; i++ )
    457             vsum += M(i,j) * M(i,j);
    458         if ( vsum > vmax ) vmax = vsum;
     466        vsum = 0;
     467        for ( i = 1; i <= rows; i++ )
     468            vsum += M(i,j) * M(i,j);
     469        if ( vsum > vmax ) vmax = vsum;
    459470    }
    460471    DEBOUTLN( cerr, "bound(lhs)^2 = " << vmax );
    … …  
    472483    int i, j;
    473484    for ( i = 1; i <= rows; i++ ) {
    474         sum = 0;
    475         for ( j = 1; j <= rows; j++ )
    476             sum += M(i,j) * M(i,j);
    477         prod *= 1 + sqrt(sum);
     485        sum = 0;
     486        for ( j = 1; j <= rows; j++ )
     487            sum += M(i,j) * M(i,j);
     488        prod *= 1 + sqrt(sum);
    478489    }
    479490    return prod;
    … …  
    495506    // triangularization
    496507    for ( i = 0; i < nrows; i++ ) {
    497         //find "pivot"
    498         for (j = i; j < nrows; j++ )
    499             if ( extmat[j][i] != 0 ) break;
    500         if ( j == nrows ) {
    501             DEBOUTLN( cerr, "solve failed" );
    502             DEBDECLEVEL( cerr, "solve" );
    503             return false;
    504         }
    505         if ( j != i ) {
    506             swap = extmat[i]; extmat[i] = extmat[j]; extmat[j] = swap;
    507         }
    508         pivotrecip = ff_inv( extmat[i][i] );
    509         rowi = extmat[i];
    510         for ( j = 0; j < ncols; j++ )
    511             rowi[j] = ff_mul( pivotrecip, rowi[j] );
    512         for ( j = i+1; j < nrows; j++ ) {
    513             rowj = extmat[j];
    514             rowpivot = rowj[i];
    515             if ( rowpivot == 0 ) continue;
    516             for ( k = i; k < ncols; k++ )
    517                 rowj[k] = ff_sub( rowj[k], ff_mul( rowpivot, rowi[k] ) );
    518         }
     508        //find "pivot"
     509        for (j = i; j < nrows; j++ )
     510            if ( extmat[j][i] != 0 ) break;
     511        if ( j == nrows ) {
     512            DEBOUTLN( cerr, "solve failed" );
     513            DEBDECLEVEL( cerr, "solve" );
     514            return false;
     515        }
     516        if ( j != i ) {
     517            swap = extmat[i]; extmat[i] = extmat[j]; extmat[j] = swap;
     518        }
     519        pivotrecip = ff_inv( extmat[i][i] );
     520        rowi = extmat[i];
     521        for ( j = 0; j < ncols; j++ )
     522            rowi[j] = ff_mul( pivotrecip, rowi[j] );
     523        for ( j = i+1; j < nrows; j++ ) {
     524            rowj = extmat[j];
     525            rowpivot = rowj[i];
     526            if ( rowpivot == 0 ) continue;
     527            for ( k = i; k < ncols; k++ )
     528                rowj[k] = ff_sub( rowj[k], ff_mul( rowpivot, rowi[k] ) );
     529        }
    519530    }
    520531    // matrix is now upper triangular with 1s down the diagonal
    521532    // back-substitute
    522533    for ( i = nrows-1; i >= 0; i-- ) {
    523         rowi = extmat[i];
    524         for ( j = 0; j < i; j++ ) {
    525             rowj = extmat[j];
    526             rowpivot = rowj[i];
    527             if ( rowpivot == 0 ) continue;
    528             for ( k = i; k < ncols; k++ )
    529                 rowj[k] = ff_sub( rowj[k], ff_mul( rowpivot, rowi[k] ) );
    530             // for (k=nrows; k<ncols; k++) rowj[k] = ff_sub(rowj[k], ff_mul(rowpivot, rowi[k]));
    531         }
     534        rowi = extmat[i];
     535        for ( j = 0; j < i; j++ ) {
     536            rowj = extmat[j];
     537            rowpivot = rowj[i];
     538            if ( rowpivot == 0 ) continue;
     539            for ( k = i; k < ncols; k++ )
     540                rowj[k] = ff_sub( rowj[k], ff_mul( rowpivot, rowi[k] ) );
     541            // for (k=nrows; k<ncols; k++) rowj[k] = ff_sub(rowj[k], ff_mul(rowpivot, rowi[k]));
     542        }
    532543    }
    533544    DEBOUTLN( cerr, "solve succeeded" );
    … …  
    549560
    550561    for ( i = 0; i < n; i++ ) {
    551         //find "pivot"
    552         for (j = i; j < n; j++ )
    553             if ( extmat[j][i] != 0 ) break;
    554         if ( j == n ) return 0;
    555         if ( j != i ) {
    556             multiplier = ff_neg( multiplier );
    557             swap = extmat[i]; extmat[i] = extmat[j]; extmat[j] = swap;
    558         }
    559         rowi = extmat[i];
    560         rowii = rowi[i];
    561         for ( j = i+1; j < n; j++ ) {
    562             rowj = extmat[j];
    563             rowji = rowj[i];
    564             if ( rowji == 0 ) continue;
    565             divisor = ff_mul( divisor, rowii );
    566             for ( k = i; k < n; k++ )
    567                 rowj[k] = ff_sub( ff_mul( rowj[k], rowii ), ff_mul( rowi[k], rowji ) );
    568         }
     562        //find "pivot"
     563        for (j = i; j < n; j++ )
     564            if ( extmat[j][i] != 0 ) break;
     565        if ( j == n ) return 0;
     566        if ( j != i ) {
     567            multiplier = ff_neg( multiplier );
     568            swap = extmat[i]; extmat[i] = extmat[j]; extmat[j] = swap;
     569        }
     570        rowi = extmat[i];
     571        rowii = rowi[i];
     572        for ( j = i+1; j < n; j++ ) {
     573            rowj = extmat[j];
     574            rowji = rowj[i];
     575            if ( rowji == 0 ) continue;
     576            divisor = ff_mul( divisor, rowii );
     577            for ( k = i; k < n; k++ )
     578                rowj[k] = ff_sub( ff_mul( rowj[k], rowii ), ff_mul( rowi[k], rowji ) );
     579        }
    569580    }
    570581    multiplier = ff_mul( multiplier, ff_inv( divisor ) );
    571582    for ( i = 0; i < n; i++ )
    572         multiplier = ff_mul( multiplier, extmat[i][i] );
     583        multiplier = ff_mul( multiplier, extmat[i][i] );
    573584    return multiplier;
    574585}
    … …  
    582593
    583594    for ( i = 1; i <= n; i++ )
    584         Q *= ( z - a[i] );
     595        Q *= ( z - a[i] );
    585596    for ( i = 1; i <= n; i++ ) {
    586         q = Q / ( z - a[i] );
    587         p = q / q( a[i], z );
    588         x[i] = 0;
    589         for ( j = p; j.hasTerms(); j++ )
    590             x[i] += w[j.exp()+1] * j.coeff();
    591     }
    592 }
     597        q = Q / ( z - a[i] );
     598        p = q / q( a[i], z );
     599        x[i] = 0;
     600        for ( j = p; j.hasTerms(); j++ )
     601            x[i] += w[j.exp()+1] * j.coeff();
     602    }
     603}

  • factory/cf_switches.cc

    r6c71ca r6f62c3  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: cf_switches.cc,v 1.7 2005-05-03 09:35:34 Singular Exp $ */
     2/* $Id: cf_switches.cc,v 1.8 2007-09-26 09:17:40 Singular Exp $ */
    33
    44//{{{ docu
    … …  
    3030#ifdef HAVE_NTL
    3131  On(SW_USE_NTL);
     32  On(SW_USE_CHINREM_GCD);
    3233  //Off(SW_USE_NTL_GCD_0);
    3334  //Off(SW_USE_NTL_GCD_P);

  • factory/cf_switches.h

    r6c71ca r6f62c3  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: cf_switches.h,v 1.9 2005-05-03 09:35:34 Singular Exp $ */
     2/* $Id: cf_switches.h,v 1.10 2007-09-26 09:17:40 Singular Exp $ */
    33
    44#ifndef INCL_CF_SWITCHES_H
    … …  
    1919//
    2020//}}}
    21 const int CFSwitchesMax = 12;
     21const int CFSwitchesMax = 13;
    2222//}}}
    2323
Note: See TracChangeset for help on using the changeset viewer.