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                 :            :  * Multivariate root finding. Implements "FindZero" from [OKTB23].
      11                 :            :  *
      12                 :            :  * [OKTB23]: https://doi.org/10.1007/978-3-031-37703-7_8
      13                 :            :  */
      14                 :            : 
      15                 :            : #ifdef CVC5_USE_COCOA
      16                 :            : 
      17                 :            : #include "theory/ff/multi_roots.h"
      18                 :            : 
      19                 :            : #include <CoCoA/BigIntOps.H>
      20                 :            : #include <CoCoA/RingFp.H>
      21                 :            : #include <CoCoA/SparsePolyOps-MinPoly.H>
      22                 :            : #include <CoCoA/SparsePolyOps-RingElem.H>
      23                 :            : #include <CoCoA/SparsePolyOps-ideal.H>
      24                 :            : #include <CoCoA/ring.H>
      25                 :            : 
      26                 :            : #include <algorithm>
      27                 :            : #include <memory>
      28                 :            : #include <sstream>
      29                 :            : 
      30                 :            : #include "theory/ff/cocoa_util.h"
      31                 :            : #include "theory/ff/uni_roots.h"
      32                 :            : #include "theory/ff/util.h"
      33                 :            : #include "util/resource_manager.h"
      34                 :            : 
      35                 :            : namespace cvc5::internal {
      36                 :            : namespace theory {
      37                 :            : namespace ff {
      38                 :            : 
      39                 :        950 : AssignmentEnumerator::~AssignmentEnumerator() = default;
      40                 :            : 
      41                 :        602 : ListEnumerator::ListEnumerator(std::vector<CoCoA::RingElem>&& options)
      42                 :        602 :     : d_remainingOptions(std::move(options))
      43                 :            : {
      44                 :        602 :   std::reverse(d_remainingOptions.begin(), d_remainingOptions.end());
      45                 :        602 : }
      46                 :            : 
      47                 :       1204 : ListEnumerator::~ListEnumerator() {};
      48                 :            : 
      49                 :        614 : std::optional<CoCoA::RingElem> ListEnumerator::next()
      50                 :            : {
      51         [ +  + ]:        614 :   if (d_remainingOptions.empty())
      52                 :            :   {
      53                 :         15 :     return {};
      54                 :            :   }
      55                 :            :   else
      56                 :            :   {
      57                 :        599 :     CoCoA::RingElem v = d_remainingOptions.back();
      58                 :        599 :     d_remainingOptions.pop_back();
      59                 :        599 :     return v;
      60                 :        599 :   }
      61                 :            : }
      62                 :            : 
      63                 :          0 : std::string ListEnumerator::name() { return "list"; }
      64                 :            : 
      65                 :        601 : std::unique_ptr<ListEnumerator> factorEnumerator(CoCoA::RingElem univariatePoly)
      66                 :            : {
      67                 :        601 :   long varIdx = CoCoA::UnivariateIndetIndex(univariatePoly);
      68 [ -  + ][ -  + ]:        601 :   Assert(varIdx >= 0);
                 [ -  - ]
      69         [ +  - ]:        601 :   Trace("ff::model::factor") << "roots for: " << univariatePoly << std::endl;
      70                 :        601 :   std::vector<CoCoA::RingElem> theRoots = roots(univariatePoly);
      71                 :        601 :   std::vector<CoCoA::RingElem> linears{};
      72                 :        601 :   CoCoA::RingElem var = CoCoA::indet(CoCoA::owner(univariatePoly), varIdx);
      73         [ +  + ]:       1259 :   for (const auto& r : theRoots)
      74                 :            :   {
      75                 :        658 :     linears.push_back(var - r);
      76                 :            :   }
      77                 :       1202 :   return std::make_unique<ListEnumerator>(std::move(linears));
      78                 :        601 : }
      79                 :            : 
      80                 :        348 : RoundRobinEnumerator::RoundRobinEnumerator(
      81                 :        348 :     const std::vector<CoCoA::RingElem>& vars, const CoCoA::ring& ring)
      82                 :        348 :     : d_vars(vars),
      83                 :        348 :       d_ring(ring),
      84                 :        348 :       d_idx(),
      85                 :        696 :       d_maxIdx(
      86                 :        696 :           CoCoA::power(CoCoA::characteristic(ring), CoCoA::LogCardinality(ring))
      87                 :        348 :           * vars.size())
      88                 :            : {
      89                 :        348 : }
      90                 :            : 
      91                 :        696 : RoundRobinEnumerator::~RoundRobinEnumerator() {}
      92                 :            : 
      93                 :        591 : std::optional<CoCoA::RingElem> RoundRobinEnumerator::next()
      94                 :            : {
      95                 :        591 :   std::optional<CoCoA::RingElem> ret{};
      96         [ +  + ]:        591 :   if (d_idx != d_maxIdx)
      97                 :            :   {
      98                 :        589 :     size_t whichVar = d_idx % d_vars.size();
      99                 :        589 :     CoCoA::BigInt whichVal = d_idx / d_vars.size();
     100                 :        589 :     CoCoA::RingElem val = d_ring->myZero();
     101                 :        589 :     val += whichVal;
     102                 :        589 :     ret = d_vars[whichVar] - val;
     103                 :        589 :     ++d_idx;
     104                 :        589 :   }
     105                 :        591 :   return ret;
     106                 :          0 : }
     107                 :            : 
     108                 :          0 : std::string RoundRobinEnumerator::name() { return "round-robin"; }
     109                 :            : 
     110                 :       1123 : bool isUnsat(const CoCoA::ideal& ideal)
     111                 :            : {
     112                 :       1123 :   const auto& gens = CoCoA::GBasis(ideal);
     113         [ +  - ]:       1317 :   return gens.size() == 1 && !CoCoA::IsZero(gens[0])
     114 [ +  + ][ +  + ]:       1317 :          && CoCoA::deg(gens[0]) <= 0;
     115                 :            : }
     116                 :            : 
     117                 :            : // True if a Groebner basis for `ideal` is available without (re)computation.
     118                 :            : //
     119                 :            : // Since CoCoALib 0.99850, computing the GBasis of the zero ideal no longer
     120                 :            : // marks it as stored (CoCoA::GBasis returns the trivially empty basis but
     121                 :            : // CoCoA::HasGBasis stays false); 0.99800 did mark it. The empty basis is always
     122                 :            : // available for the zero ideal, so treat that case as "has GBasis" too.
     123                 :       2589 : bool hasGBasis(const CoCoA::ideal& ideal)
     124                 :            : {
     125 [ +  + ][ +  - ]:       2589 :   return CoCoA::HasGBasis(ideal) || CoCoA::IsZero(ideal);
     126                 :            : }
     127                 :            : 
     128                 :            : template <typename T>
     129                 :          0 : std::string ostring(const T& t)
     130                 :            : {
     131                 :          0 :   std::ostringstream o;
     132                 :          0 :   o << t;
     133                 :          0 :   return o.str();
     134                 :          0 : }
     135                 :            : 
     136                 :        943 : std::pair<size_t, CoCoA::RingElem> extractAssignment(
     137                 :            :     const CoCoA::RingElem& elem)
     138                 :            : {
     139 [ -  + ][ -  + ]:        943 :   Assert(CoCoA::deg(elem) == 1);
                 [ -  - ]
     140 [ -  + ][ -  + ]:        943 :   Assert(CoCoA::NumTerms(elem) <= 2);
                 [ -  - ]
     141                 :        943 :   const CoCoA::RingElem m = CoCoA::monic(elem);
     142                 :        943 :   long varNumber = CoCoA::UnivariateIndetIndex(elem);
     143 [ -  + ][ -  + ]:        943 :   Assert(varNumber >= 0);
                 [ -  - ]
     144                 :       2829 :   return {varNumber, -CoCoA::ConstantCoeff(m)};
     145                 :        943 : }
     146                 :            : 
     147                 :        964 : std::unordered_set<std::string> assignedVars(const CoCoA::ideal& ideal)
     148                 :            : {
     149                 :        964 :   std::unordered_set<std::string> ret{};
     150 [ -  + ][ -  + ]:        964 :   Assert(hasGBasis(ideal));
                 [ -  - ]
     151         [ +  + ]:       5024 :   for (const auto& g : CoCoA::GBasis(ideal))
     152                 :            :   {
     153         [ +  + ]:       4060 :     if (CoCoA::deg(g) == 1)
     154                 :            :     {
     155                 :       3432 :       long varNumber = CoCoA::UnivariateIndetIndex(g);
     156         [ +  + ]:       3432 :       if (varNumber >= 0)
     157                 :            :       {
     158                 :       3221 :         ret.insert(ostring(CoCoA::indet(ideal->myRing(), varNumber)));
     159                 :            :       }
     160                 :            :     }
     161                 :            :   }
     162                 :        964 :   return ret;
     163                 :          0 : }
     164                 :            : 
     165                 :        788 : bool allVarsAssigned(const CoCoA::ideal& ideal)
     166                 :            : {
     167                 :       1576 :   return assignedVars(ideal).size()
     168                 :       1576 :          == (size_t)CoCoA::NumIndets(ideal->myRing());
     169                 :            : }
     170                 :            : 
     171                 :        266 : std::unique_ptr<AssignmentEnumerator> applyRule(const CoCoA::ideal& ideal,
     172                 :            :                                                 FfStatistics* stats = nullptr)
     173                 :            : {
     174                 :        266 :   CoCoA::PolyRing polyRing(ideal->myRing());
     175 [ -  + ][ -  + ]:        266 :   Assert(!isUnsat(ideal));
                 [ -  - ]
     176                 :            :   // first, we look for super-linear univariate polynomials.
     177 [ -  + ][ -  + ]:        266 :   Assert(hasGBasis(ideal));
                 [ -  - ]
     178                 :        266 :   const auto& gens = CoCoA::GBasis(ideal);
     179         [ +  + ]:       1147 :   for (const auto& p : gens)
     180                 :            :   {
     181                 :        971 :     long varNumber = CoCoA::UnivariateIndetIndex(p);
     182 [ +  + ][ +  + ]:        971 :     if (varNumber >= 0 && CoCoA::deg(p) > 1)
                 [ +  + ]
     183                 :            :     {
     184                 :         90 :       return factorEnumerator(p);
     185                 :            :     }
     186                 :            :   }
     187                 :            :   // now, we check the dimension
     188         [ -  + ]:        176 :   if (CoCoA::IsZeroDim(ideal))
     189                 :            :   {
     190         [ -  - ]:          0 :     if (stats) ++stats->d_idealMinPoly;
     191                 :            :     // If zero-dimensional, we compute a minimal polynomial in some unset
     192                 :            :     // variable.
     193                 :          0 :     std::unordered_set<std::string> alreadySet = assignedVars(ideal);
     194         [ -  - ]:          0 :     for (const auto& var : CoCoA::indets(polyRing))
     195                 :            :     {
     196                 :          0 :       std::string varName = ostring(var);
     197         [ -  - ]:          0 :       if (!alreadySet.count(ostring(var)))
     198                 :            :       {
     199                 :          0 :         CoCoA::RingElem minPoly = CoCoA::MinPolyQuot(var, ideal, var);
     200                 :          0 :         return factorEnumerator(minPoly);
     201                 :          0 :       }
     202         [ -  - ]:          0 :     }
     203                 :          0 :     Unreachable()
     204                 :          0 :         << "There should be no unset variables in zero-dimensional ideal";
     205                 :          0 :   }
     206                 :            :   else
     207                 :            :   {
     208         [ +  + ]:        176 :     if (stats) ++stats->d_idealPosDim;
     209                 :            :     // If positive dimensional, we make a list of unset variables and
     210                 :            :     // round-robin guess.
     211                 :            :     //
     212                 :            :     // TODO(aozdemir): better model construction (cvc5-wishues/issues/138)
     213                 :        176 :     std::unordered_set<std::string> alreadySet = assignedVars(ideal);
     214                 :        176 :     std::vector<CoCoA::RingElem> toGuess{};
     215         [ +  + ]:       1172 :     for (const auto& var : CoCoA::indets(polyRing))
     216                 :            :     {
     217                 :        996 :       std::string varName = ostring(var);
     218         [ +  + ]:        996 :       if (!alreadySet.count(ostring(var)))
     219                 :            :       {
     220                 :        364 :         toGuess.push_back(var);
     221                 :            :       }
     222                 :        996 :     }
     223                 :        352 :     return std::make_unique<RoundRobinEnumerator>(toGuess,
     224                 :        352 :                                                   polyRing->myBaseRing());
     225                 :        176 :   }
     226                 :        266 : }
     227                 :            : 
     228                 :        211 : std::vector<CoCoA::RingElem> findZero(const CoCoA::ideal& initialIdeal,
     229                 :            :                                       const Env& env,
     230                 :            :                                       FfStatistics* stats)
     231                 :            : {
     232                 :        211 :   CoCoA::ring polyRing = initialIdeal->myRing();
     233                 :            :   // We maintain two stacks:
     234                 :            :   // * one of ideals
     235                 :            :   // * one of branchers
     236                 :            :   //
     237                 :            :   // If brancher B has the same index as ideal I, then B represents possible
     238                 :            :   // expansions of ideal I (equivalently, restrictions of I's variety).
     239                 :            :   //
     240                 :            :   // NB: FindZero of [OKTB23] also takes a partial map M as input. GB(I)
     241                 :            :   // implicitly represents M: GB(I) contains a univariate linear polynomial
     242                 :            :   // Xi - k, if and only iff M[Xi] = k.
     243                 :            :   //
     244                 :            :   // NB: FindZero of [OKTB23] is recursive. That recursion is flattened here
     245                 :            :   // using the two stacks. The stack of ideals represents the input to
     246                 :            :   // recursive FindZero: GB(I). The stack of branchers represents the
     247                 :            :   // continuation context (which iteration of the for loop to return to).
     248                 :            : 
     249                 :            :   // goal: find a zero for any ideal in the stack.
     250                 :        633 :   std::vector<CoCoA::ideal> ideals{initialIdeal};
     251         [ -  + ]:        211 :   if (TraceIsOn("ff::model::branch"))
     252                 :            :   {
     253         [ -  - ]:          0 :     Trace("ff::model::branch") << "init polys: " << std::endl;
     254         [ -  - ]:          0 :     for (const auto& p : CoCoA::gens(initialIdeal))
     255                 :            :     {
     256         [ -  - ]:          0 :       Trace("ff::model::branch") << " * " << p << std::endl;
     257                 :            :     }
     258                 :            :   }
     259                 :            : 
     260                 :        211 :   std::vector<std::unique_ptr<AssignmentEnumerator>> branchers{};
     261                 :            :   // while some ideal might have a zero.
     262         [ +  + ]:        863 :   while (!ideals.empty())
     263                 :            :   {
     264                 :            :     // check for timeout
     265         [ -  + ]:        849 :     if (env.getResourceManager()->outOfTime())
     266                 :            :     {
     267                 :          0 :       throw FfTimeoutException("findZero");
     268                 :            :     }
     269                 :            :     // choose one ideal
     270                 :        849 :     const auto& ideal = ideals.back();
     271                 :            :     // make sure we have a GBasis:
     272                 :        849 :     GBasisTimeout(ideal, env.getResourceManager());
     273 [ -  + ][ -  + ]:        849 :     Assert(hasGBasis(ideal));
                 [ -  - ]
     274                 :            :     // If the ideal is UNSAT, drop it.
     275         [ +  + ]:        849 :     if (isUnsat(ideal))
     276                 :            :     {
     277                 :         61 :       ideals.pop_back();
     278                 :            :     }
     279                 :            :     // If the ideal has a linear polynomial in each variable, we've found a
     280                 :            :     // variety element (a model).
     281         [ +  + ]:        788 :     else if (allVarsAssigned(ideal))
     282                 :            :     {
     283                 :        197 :       std::unordered_map<size_t, CoCoA::RingElem> varNumToValue{};
     284 [ -  + ][ -  + ]:        197 :       Assert(hasGBasis(ideal));
                 [ -  - ]
     285                 :        197 :       const auto& gens = CoCoA::GBasis(ideal);
     286                 :        197 :       size_t numIndets = CoCoA::NumIndets(polyRing);
     287 [ -  + ][ -  + ]:        197 :       Assert(gens.size() == numIndets);
                 [ -  - ]
     288         [ +  + ]:       1136 :       for (const auto& g : gens)
     289                 :            :       {
     290                 :        939 :         varNumToValue.insert(extractAssignment(g));
     291                 :            :       }
     292                 :        197 :       std::vector<CoCoA::RingElem> values{};
     293         [ +  + ]:       1136 :       for (size_t i = 0; i < numIndets; ++i)
     294                 :            :       {
     295                 :        939 :         values.push_back(varNumToValue[i]);
     296                 :            :       }
     297                 :        197 :       return values;
     298                 :        197 :     }
     299                 :            :     // If there are more ideals than branchers, branch
     300         [ +  + ]:        591 :     else if (ideals.size() > branchers.size())
     301                 :            :     {
     302 [ -  + ][ -  + ]:        266 :       Assert(ideals.size() == branchers.size() + 1);
                 [ -  - ]
     303                 :        266 :       branchers.push_back(applyRule(ideal, stats));
     304         [ +  - ]:        532 :       Trace("ff::model::branch")
     305 [ -  + ][ -  - ]:        266 :           << "brancher: " << branchers.back()->name() << std::endl;
     306         [ -  + ]:        266 :       if (TraceIsOn("ff::model::branch"))
     307                 :            :       {
     308         [ -  - ]:          0 :         Trace("ff::model::branch") << "ideal polys: " << std::endl;
     309         [ -  - ]:          0 :         for (const auto& p : CoCoA::gens(ideal))
     310                 :            :         {
     311         [ -  - ]:          0 :           Trace("ff::model::branch") << " * " << p << std::endl;
     312                 :            :         }
     313                 :            :       }
     314                 :            :     }
     315                 :            :     // Otherwise, this ideal should have a brancher; get the next branch
     316                 :            :     else
     317                 :            :     {
     318 [ -  + ][ -  + ]:        325 :       Assert(ideals.size() == branchers.size());
                 [ -  - ]
     319                 :        325 :       std::optional<CoCoA::RingElem> choicePoly = branchers.back()->next();
     320                 :            :       // construct a new ideal from the branch
     321         [ +  + ]:        325 :       if (choicePoly.has_value())
     322                 :            :       {
     323         [ +  - ]:        626 :         Trace("ff::model::branch")
     324                 :          0 :             << "level: " << branchers.size()
     325 [ -  + ][ -  - ]:        313 :             << ", brancher: " << branchers.back()->name()
     326                 :        313 :             << ", branch: " << choicePoly.value() << std::endl;
     327 [ -  + ][ -  + ]:        313 :         Assert(hasGBasis(ideal));
                 [ -  - ]
     328                 :        313 :         std::vector<CoCoA::RingElem> newGens = CoCoA::GBasis(ideal);
     329                 :        313 :         newGens.push_back(choicePoly.value());
     330                 :        313 :         ideals.push_back(CoCoA::ideal(newGens));
     331                 :        313 :       }
     332                 :            :       // or drop this ideal & brancher if we're out of branches.
     333                 :            :       else
     334                 :            :       {
     335                 :         12 :         branchers.pop_back();
     336                 :         12 :         ideals.pop_back();
     337                 :            :       }
     338                 :        325 :     }
     339                 :            :   }
     340                 :            :   // Could not find any solution; return empty.
     341                 :         14 :   return {};
     342                 :        211 : }
     343                 :            : 
     344                 :            : }  // namespace ff
     345                 :            : }  // namespace theory
     346                 :            : }  // namespace cvc5::internal
     347                 :            : 
     348                 :            : #endif /* CVC5_USE_COCOA */