Branch data     Line data    Source code

       1                 :            : /******************************************************************************
       2                 :            :  * This file is part of the cvc5 project.
       3                 :            :  *
       4                 :            :  * Copyright (c) 2009-2026 by the authors listed in the file AUTHORS
       5                 :            :  * in the top-level source directory and their institutional affiliations.
       6                 :            :  * All rights reserved.  See the file COPYING in the top-level source
       7                 :            :  * directory for licensing information.
       8                 :            :  * ****************************************************************************
       9                 :            :  *
      10                 :            :  * Root finding for univariate polynomials in prime fields.
      11                 :            :  *
      12                 :            :  * Uses CoCoA for word-sized fields.
      13                 :            :  *
      14                 :            :  * Implements Rabin root-finding for larger ones.
      15                 :            :  *
      16                 :            :  * Reference: https://en.wikipedia.org/wiki/Berlekamp%E2%80%93Rabin_algorithm
      17                 :            :  */
      18                 :            : 
      19                 :            : #ifdef CVC5_USE_COCOA
      20                 :            : 
      21                 :            : #include "theory/ff/uni_roots.h"
      22                 :            : 
      23                 :            : #include <CoCoA/BigInt.H>
      24                 :            : #include <CoCoA/BigIntOps.H>
      25                 :            : #include <CoCoA/PolyRing.H>
      26                 :            : #include <CoCoA/RingZZ.H>
      27                 :            : #include <CoCoA/SmallFpImpl.H>
      28                 :            : #include <CoCoA/SparsePolyOps-RingElem.H>
      29                 :            : #include <CoCoA/factor.H>
      30                 :            : #include <CoCoA/factorization.H>
      31                 :            : #include <CoCoA/random.H>
      32                 :            : #include <CoCoA/ring.H>
      33                 :            : 
      34                 :            : #include <sstream>
      35                 :            : #include <unordered_map>
      36                 :            : #include <vector>
      37                 :            : 
      38                 :            : #include "base/output.h"
      39                 :            : #include "smt/assertions.h"
      40                 :            : 
      41                 :            : namespace cvc5::internal {
      42                 :            : namespace theory {
      43                 :            : namespace ff {
      44                 :            : 
      45                 :            : // Reduce b modulo m, for polynomials b and m.
      46                 :       8046 : CoCoA::RingElem redMod(CoCoA::RingElem b, CoCoA::RingElem m)
      47                 :            : {
      48                 :      24138 :   std::vector<CoCoA::RingElem> mm = {m};
      49                 :      16092 :   return CoCoA::NR(b, mm);
      50                 :       8046 : }
      51                 :            : 
      52                 :            : // Compute b^e modulo m.
      53                 :            : //
      54                 :            : // uses repeated squaring with reductions by m in each step
      55                 :         21 : CoCoA::RingElem powerMod(CoCoA::RingElem b, CoCoA::BigInt e, CoCoA::RingElem m)
      56                 :            : {
      57                 :         21 :   CoCoA::RingElem acc = CoCoA::owner(b)->myOne();
      58                 :         21 :   CoCoA::RingElem bPower = b;
      59         [ +  + ]:       4363 :   while (!CoCoA::IsZero(e))
      60                 :            :   {
      61         [ +  + ]:       4342 :     if (CoCoA::IsOdd(e))
      62                 :            :     {
      63                 :       3704 :       acc *= bPower;
      64                 :       3704 :       acc = redMod(acc, m);
      65                 :            :     }
      66                 :       4342 :     bPower *= bPower;
      67                 :       4342 :     bPower = redMod(bPower, m);
      68                 :       4342 :     e /= 2;
      69                 :            :   }
      70                 :         42 :   return acc;
      71                 :         21 : }
      72                 :            : 
      73                 :         20 : CoCoA::RingElem distinctRootsPoly(CoCoA::RingElem f)
      74                 :            : {
      75                 :         20 :   CoCoA::ring ring = CoCoA::owner(f);
      76                 :         20 :   CoCoA::ring field = CoCoA::owner(f)->myBaseRing();
      77                 :         20 :   int idx = CoCoA::UnivariateIndetIndex(f);
      78 [ -  + ][ -  + ]:         20 :   Assert(idx >= 0);
                 [ -  - ]
      79                 :         20 :   CoCoA::RingElem x = CoCoA::indet(ring, idx);
      80                 :            :   CoCoA::BigInt q =
      81                 :         20 :       CoCoA::power(CoCoA::characteristic(field), CoCoA::LogCardinality(field));
      82                 :         40 :   CoCoA::RingElem fieldPoly = powerMod(x, q, f) - x;
      83                 :         40 :   return gcd(f, fieldPoly);
      84                 :         20 : }
      85                 :            : 
      86                 :            : // get a string for `t` from its input
      87                 :            : template <typename T>
      88                 :        666 : std::string sstring(const T& t)
      89                 :            : {
      90                 :        666 :   std::ostringstream o;
      91                 :        666 :   o << t;
      92                 :       1332 :   return o.str();
      93                 :        666 : }
      94                 :            : 
      95                 :            : // sorting based on strings because CoCoA can't compare field elements:
      96                 :            : // it doesn't regard a integer quotient ring as an ordered domain.
      97                 :        607 : std::vector<CoCoA::RingElem> sortHack(
      98                 :            :     const std::vector<CoCoA::RingElem>& values)
      99                 :            : {
     100                 :        607 :   std::vector<std::string> strs;
     101                 :        607 :   std::unordered_map<std::string, size_t> origIndices;
     102         [ +  + ]:       1273 :   for (const auto& v : values)
     103                 :            :   {
     104                 :        666 :     std::string s = sstring(v);
     105                 :        666 :     origIndices.emplace(s, strs.size());
     106                 :        666 :     strs.push_back(s);
     107                 :        666 :   }
     108                 :        607 :   std::sort(strs.begin(), strs.end());
     109                 :        607 :   std::vector<CoCoA::RingElem> output;
     110         [ +  + ]:       1273 :   for (const auto& s : strs)
     111                 :            :   {
     112                 :        666 :     output.push_back(values[origIndices[s]]);
     113                 :            :   }
     114                 :       1214 :   return output;
     115                 :        607 : }
     116                 :            : 
     117                 :        607 : std::vector<CoCoA::RingElem> roots(CoCoA::RingElem f)
     118                 :            : {
     119                 :        607 :   CoCoA::ring ring = CoCoA::owner(f);
     120                 :        607 :   CoCoA::ring field = CoCoA::owner(f)->myBaseRing();
     121                 :        607 :   int idx = CoCoA::UnivariateIndetIndex(f);
     122 [ -  + ][ -  + ]:        607 :   Assert(idx >= 0);
                 [ -  - ]
     123                 :        607 :   CoCoA::RingElem x = CoCoA::indet(ring, idx);
     124                 :        607 :   CoCoA::BigInt q = CoCoA::characteristic(field);
     125                 :        607 :   std::vector<CoCoA::RingElem> output;
     126                 :            : 
     127                 :            :   // CoCoA has a good factorization routine, but it only works for small fields.
     128                 :        607 :   bool isSmall = false;
     129                 :            :   {
     130                 :            :     // I don't know how to check directly if their small field impl applies, so
     131                 :            :     // we try.
     132                 :            :     try
     133                 :            :     {
     134                 :        607 :       CoCoA::SmallFpImpl ModP(CoCoA::ConvertTo<long>(q));
     135                 :        595 :       isSmall = true;
     136                 :            :     }
     137         [ -  + ]:         12 :     catch (const CoCoA::ErrorInfo&)
     138                 :            :     {
     139                 :         12 :     }
     140                 :            :   }
     141         [ +  + ]:        607 :   if (isSmall)
     142                 :            :   {
     143                 :            :     // Use CoCoA
     144                 :        595 :     const auto factors = CoCoA::factor(f);
     145         [ +  + ]:       1271 :     for (const auto& factor : factors.myFactors())
     146                 :            :     {
     147         [ +  + ]:        676 :       if (CoCoA::deg(factor) == 1)
     148                 :            :       {
     149 [ -  + ][ -  + ]:        651 :         Assert(CoCoA::IsOne(CoCoA::LC(factor)));
                 [ -  - ]
     150                 :        651 :         output.push_back(-CoCoA::ConstantCoeff(factor));
     151                 :            :       }
     152                 :            :     }
     153                 :        595 :   }
     154                 :            :   else
     155                 :            :   {
     156                 :            :     // Rabin root finding
     157                 :            : 
     158                 :            :     // needed because of the random sampling below
     159 [ -  + ][ -  + ]:         12 :     Assert(CoCoA::LogCardinality(field) == 1);
                 [ -  - ]
     160 [ -  + ][ -  + ]:         12 :     Assert(CoCoA::IsOdd(q));
                 [ -  - ]
     161                 :         12 :     CoCoA::BigInt s = q / 2;
     162                 :            :     // Reduce the problem to factoring a product of linears.
     163                 :            :     // We need to factor everything in this queue.
     164                 :         48 :     std::vector<CoCoA::RingElem> toFactor{distinctRootsPoly(f)};
     165                 :            : 
     166                 :            :     // While there is more to factor.
     167         [ +  + ]:         32 :     while (toFactor.size())
     168                 :            :     {
     169                 :            :       // Grab a product of linears to factor
     170                 :         20 :       CoCoA::RingElem p = toFactor.back();
     171                 :         20 :       toFactor.pop_back();
     172         [ +  - ]:         20 :       Trace("ff::roots") << "toFactor " << p << std::endl;
     173         [ +  + ]:         20 :       if (CoCoA::deg(p) == 0)
     174                 :            :       {
     175                 :            :         // It's a constant: no factors
     176                 :            :       }
     177         [ +  + ]:         16 :       else if (CoCoA::ConstantCoeff(p) == 0)
     178                 :            :       {
     179                 :            :         // It has a zero root
     180                 :          6 :         output.push_back(CoCoA::ConstantCoeff(p));
     181                 :          6 :         toFactor.push_back(p / x);
     182                 :            :       }
     183         [ +  + ]:         10 :       else if (CoCoA::deg(p) == 1)
     184                 :            :       {
     185                 :            :         // It is linear
     186                 :          9 :         output.push_back(-CoCoA::ConstantCoeff(p));
     187                 :            :       }
     188                 :            :       else
     189                 :            :       {
     190                 :            :         // Its super-linear, without a zero-root
     191                 :            :         while (true)
     192                 :            :         {
     193                 :            :           // guess random delta
     194                 :          2 :           CoCoA::BigInt deltaInt = CoCoA::RandomBigInt(CoCoA::BigInt(0), q - 1);
     195                 :          1 :           CoCoA::RingElem delta = CoCoA::RingElem(field, deltaInt);
     196                 :            : 
     197                 :            :           // construct split(X) = (S - delta)^(q/2) - 1.
     198                 :            :           //
     199                 :            :           // the probability (over delta) that some fix field element is a root
     200                 :            :           // of split is exactly 1/2.
     201                 :          2 :           CoCoA::RingElem split = powerMod(x - delta, s, p) - 1;
     202                 :            : 
     203                 :            :           // Is the number of common roots between split and p at least 1 and
     204                 :            :           // less than p's degree?
     205                 :          1 :           CoCoA::RingElem h = gcd(p, split);
     206 [ +  - ][ +  - ]:          1 :           if (0 < CoCoA::deg(h) && CoCoA::deg(h) < CoCoA::deg(p))
                 [ +  - ]
     207                 :            :           {
     208                 :            :             // yes: replace p with (p / h) and h in the queue.
     209                 :          1 :             toFactor.push_back(h);
     210                 :          1 :             toFactor.push_back(p / h);
     211                 :          1 :             break;
     212                 :            :           }
     213                 :            :           // no: guess a new delta
     214 [ -  + ][ -  + ]:          4 :         }
         [ -  + ][ -  + ]
     215                 :            :       }
     216                 :         20 :     }
     217                 :         12 :   }
     218                 :       1214 :   return sortHack(output);
     219                 :        607 : }
     220                 :            : 
     221                 :            : }  // namespace ff
     222                 :            : }  // namespace theory
     223                 :            : }  // namespace cvc5::internal
     224                 :            : 
     225                 :            : #endif /* CVC5_USE_COCOA */