Branch data     Line data    Source code

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