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