text:
/**********************************************************\
Original Author: Richard Bateman (taxilian)

Created:    Oct 19, 2009
License:    Dual license model; choose one of two:
            New BSD License
            http://www.opensource.org/licenses/bsd-license.php
            - or -
            GNU Lesser General Public License, version 2.1
            http://www.gnu.org/licenses/lgpl-2.1.html

Copyright 2009 Richard Bateman, Firebreath development team
\**********************************************************/

#ifdef FB_MACOSX
#include <dlfcn.h>
#include <mach-o/dyld.h>
#endif
#include <cstdio>
#include "NpapiPlugin.h"
#include "FactoryBase.h"
#include "NpapiBrowserHost.h"
#include <boost/shared_ptr.hpp>
#include "precompiled_headers.h" // On windows, everything above this line in PCH
#include "AsyncFunctionCall.h"
#include "PluginInfo.h"
#include "SafeQueue.h"
#include "NpapiPluginModule.h"

#if FB_WIN
#  include "Win/NpapiBrowserHostAsyncWin.h"
#elif FB_MACOSX
#  include "Mac/OneShotManager.h"
#endif

using namespace FB::Npapi;

namespace
{
    bool needAsyncCallsWorkaround(NPP npp, NPNetscapeFuncs* funcs)
    {
        bool result(false);
        // work-around detection here
#if FB_WIN
        const char* const cstrUserAgent = funcs->uagent(npp);
        if(cstrUserAgent) {
            const std::string userAgent(cstrUserAgent);
            // Safari 5.1 NPN_PluginThreadAsyncCall doesn't seem to work anymore; use the workaround
            result = userAgent.find("Safari") != std::string::npos;
        }
#endif
        return result || (funcs->version < NPVERS_HAS_PLUGIN_THREAD_ASYNC_CALL);
    }

    bool asyncCallsWorkaround(NPP npp, NPNetscapeFuncs* funcs = 0)
    {
        static const bool useWorkaround = (funcs) ? needAsyncCallsWorkaround(npp, funcs) : false;
        //static const bool useWorkaround = false;
        return useWorkaround;
    }

    NpapiBrowserHostPtr createBrowserHost(NpapiPluginModule* module, NPP npp)
    {
        try {
            NpapiBrowserHostPtr host;
            NPNetscapeFuncs& npnFuncs = module->NPNFuncs;

            if(asyncCallsWorkaround(npp, &npnFuncs)) {
                npnFuncs.pluginthreadasynccall = NULL;
    #if FB_WIN
                NpapiBrowserHostPtr host(boost::make_shared<NpapiBrowserHostAsyncWin>(module, npp));
                return host;
    #else
                // no work-around for this platform
                NpapiBrowserHostPtr host(boost::make_shared<NpapiBrowserHost>(module, npp));
                return host;
    #endif
            } else {
                NpapiBrowserHostPtr host(boost::make_shared<NpapiBrowserHost>(module, npp));
                return host;
            }
        } catch (...) {
            // This function must not return an exception
            return NpapiBrowserHostPtr();
        }
    }

    class NpapiPDataHolder
    {
    private:
        NpapiBrowserHostPtr m_host;
        boost::shared_ptr<NpapiPlugin> m_plugin;

    public:
        NpapiPDataHolder(NpapiBrowserHostPtr host, boost::shared_ptr<NpapiPlugin> plugin)
          : m_host(host), m_plugin(plugin)
        {
#ifdef FB_MACOSX
            FB::OneShotManager::getInstance().npp_register(m_host->getContextID());
#endif
        }
        ~NpapiPDataHolder() {
#ifdef FB_MACOSX
            FB::OneShotManager::getInstance().npp_unregister(m_host->getContextID());
#endif
        }

        NpapiBrowserHostPtr getHost() const {
            return m_host;
        }
        NpapiPluginPtr getPlugin() const {
            return m_plugin;
        }
    };

    bool validInstance(NPP instance)
    {
        return instance != NULL && instance->pdata != NULL;
    }

    NpapiPDataHolder* getHolder(NPP instance)
    {
        if (validInstance(instance))
            return static_cast<NpapiPDataHolder*>(instance->pdata);
        return NULL;
    }
}

NpapiPluginPtr FB::Npapi::getPlugin(NPP instance)
{
    if (NpapiPDataHolder* holder = getHolder(instance))
        return holder->getPlugin();
    return NpapiPluginPtr();
}

#ifdef FB_MACOSX
// This is used on mac snow leopard safari
void NpapiPluginModule::scheduleAsyncCallback(NPP instance, void (*func)(void *), void *userData)
{
    FB::OneShotManager::getInstance().npp_scheduleAsyncCallback(instance, func, userData);
}
#endif

// These are the static NPP_ functions; NPP_New and NPP_Destroy create and destroy the
// plugin, the rest are wrappers that dereference NPP->pdata to get at the plugin object
// and proxy the call to there.
NPError NpapiPluginModule::NPP_New(NPMIMEType pluginType, NPP instance, uint16_t mode, int16_t argc,
                                   char* argn[], char* argv[], NPSavedData* saved)
{
    FBLOG_INFO("NPAPI", "NPP_New: " << (void*) instance);
    if (instance == NULL) {
        return NPERR_INVALID_INSTANCE_ERROR;
    }

    try
    {
#ifdef FB_MACOSX
        // Helps with certain weird embedding cases
		Dl_info info;
		dladdr(__builtin_return_address(0), &info);
        NpapiPluginModule *module = NpapiPluginModule::GetModule(info.dli_fbase);
#else
        NpapiPluginModule *module = NpapiPluginModule::GetModule(0);
#endif
        //printf("%p NpapiPluginModule::%s()\n", module, __func__);
        NPNetscapeFuncs& npnFuncs = module->NPNFuncs;

        NpapiBrowserHostPtr host(createBrowserHost(module, instance));

        host->setBrowserFuncs(&(npnFuncs));

        // TODO: We should probably change this and pass the MIMEType into _getNpapiPlugin instead
        // of into init later so that we can optionally return a different plugin type depending
        // on the specific mimetype
        NpapiPluginPtr plugin(getFactoryInstance()->createNpapiPlugin(host, pluginType));
        if (!plugin) {
            return NPERR_OUT_OF_MEMORY_ERROR;
        }

        NpapiPDataHolder* holder = new NpapiPDataHolder(host, plugin);
        instance->pdata = static_cast<void*>(holder);

        plugin->init(pluginType, argc, argn, argv);
    }
    catch (const PluginCreateError &e)
    {
        printf("%s\n", e.what());
        return NPERR_INCOMPATIBLE_VERSION_ERROR;
    }
    catch (const std::bad_alloc& e)
    {
        printf("%s\n", e.what());
        return NPERR_OUT_OF_MEMORY_ERROR;
    }
    catch (const std::exception& e)
    {
        printf("%s\n", e.what());
        return NPERR_GENERIC_ERROR;
    }

    return NPERR_NO_ERROR;
}

NPError NpapiPluginModule::NPP_Destroy(NPP instance, NPSavedData** save)
{
    FBLOG_INFO("NPAPI", "NPP_Destroy: " << (void*) instance);
    if (!validInstance(instance)) {
        return NPERR_INVALID_INSTANCE_ERROR;
    }
    NpapiBrowserHostWeakPtr weakHost;

    if (NpapiPDataHolder* holder = getHolder(instance)) {
        NpapiBrowserHostPtr host(holder->getHost());
        weakHost = host;
        if (host)
            host->shutdown();

        if (NpapiPluginPtr plugin = holder->getPlugin())
            plugin->shutdown();

        instance->pdata = NULL;
        delete holder; // Destroy plugin
        // host should be destroyed when it goes out of scope here
    } else {
        return NPERR_GENERIC_ERROR;
    }
    // If this assertion fails, you probably have a circular reference
    // to your BrowserHost object somewhere -- the host should be gone
    // by this point. This assertion is warning you of a bug.
    assert(weakHost.expired());

    return NPERR_NO_ERROR;
}

NPError NpapiPluginModule::NPP_SetWindow(NPP instance, NPWindow* window)
{
    FBLOG_TRACE("NPAPI", (void*) instance);
    if (!validInstance(instance)) {
        return NPERR_INVALID_INSTANCE_ERROR;
    }

    if (pluginGuiEnabled())
        if (NpapiPluginPtr plugin = getPlugin(instance)) {
            return plugin->SetWindow(window);
        }

    return NPERR_NO_ERROR;
}

NPError NpapiPluginModule::NPP_NewStream(NPP instance, NPMIMEType type, NPStream* stream,
                                       NPBool seekable, uint16_t* stype)
{
    FBLOG_INFO("NPAPI", (void*) instance);
    if (!validInstance(instance)) {
        return NPERR_INVALID_INSTANCE_ERROR;
    }

    if (NpapiPluginPtr plugin = getPlugin(instance)) {
        return plugin->NewStream(type, stream, seekable, stype);
    } else {
        return NPERR_GENERIC_ERROR;
    }
}

NPError NpapiPluginModule::NPP_DestroyStream(NPP instance, NPStream* stream, NPReason reason)
{
    FBLOG_INFO("NPAPI", (void*) instance);
    if (!validInstance(instance)) {
        return NPERR_INVALID_INSTANCE_ERROR;
    }

    if (NpapiPluginPtr plugin = getPlugin(instance)) {
        return plugin->DestroyStream(stream, reason);
    } else {
        return NPERR_GENERIC_ERROR;
    }
}

int32_t NpapiPluginModule::NPP_WriteReady(NPP instance, NPStream* stream)
{
    FBLOG_INFO("NPAPI",(void*) instance);
    if (!validInstance(instance)) {
        return NPERR_INVALID_INSTANCE_ERROR;
    }

    if (NpapiPluginPtr plugin = getPlugin(instance)) {
        return plugin->WriteReady(stream);
    } else {
        return NPERR_GENERIC_ERROR;
    }
}

int32_t NpapiPluginModule::NPP_Write(NPP instance, NPStream* stream, int32_t offset, int32_t len,
                                 void* buffer)
{
    FBLOG_INFO("NPAPI", (void*) instance);
    if (!validInstance(instance)) {
        return NPERR_INVALID_INSTANCE_ERROR;
    }

    if (NpapiPluginPtr plugin = getPlugin(instance)) {
        return plugin->Write(stream, offset, len, buffer);
    } else {
        return NPERR_GENERIC_ERROR;
    }
}

void NpapiPluginModule::NPP_StreamAsFile(NPP instance, NPStream* stream, const char* fname)
{
    FBLOG_INFO("NPAPI", (void*) instance);
    if (!validInstance(instance)) {
        return;
    }

    if (NpapiPluginPtr plugin = getPlugin(instance)) {
        plugin->StreamAsFile(stream, fname);
    }
}

void NpapiPluginModule::NPP_Print(NPP instance, NPPrint* platformPrint)
{
    FBLOG_INFO("NPAPI", (void*) instance);
    if (!validInstance(instance)) {
        return;
    }

    if (NpapiPluginPtr plugin = getPlugin(instance)) {
        plugin->Print(platformPrint);
    }
}

int16_t NpapiPluginModule::NPP_HandleEvent(NPP instance, void* event)
{
    FBLOG_TRACE("NPAPI", (void*) instance);
    if (!validInstance(instance)) {
        return 0;
    }

    if (NpapiPluginPtr plugin = getPlugin(instance)) {
        return plugin->HandleEvent(event);
    } else {
        return 0;
    }
}

void NpapiPluginModule::NPP_URLNotify(NPP instance, const char* url, NPReason reason,
                                    void* notifyData)
{
    FBLOG_INFO("NPAPI", (void*) instance);
    if (!validInstance(instance)) {
        return;
    }

    if (NpapiPluginPtr plugin = getPlugin(instance)) {
        plugin->URLNotify(url, reason, notifyData);
    }
}

NPError NpapiPluginModule::NPP_GetValue(NPP instance, NPPVariable variable, void *value)
{
    FBLOG_TRACE("NPAPI", (void*) instance);
    // These values may depend on the mimetype of the plugin
    if (!validInstance(instance)) {
        switch (variable) {
        case NPPVpluginNameString: {
            static const std::string pluginName = getFactoryInstance()->getPluginName("");
            *((const char **)value) = pluginName.c_str();
            break; }
        case NPPVpluginDescriptionString: {
            // If nothing is instantiated, use the first description found
            static const std::string pluginDesc = getFactoryInstance()->getPluginDescription("");
            *((const char **)value) = pluginDesc.c_str();
            break; }
        default:
            return NPERR_GENERIC_ERROR;
        }
        return NPERR_NO_ERROR;
    }

    if (NpapiPluginPtr plugin = getPlugin(instance)) {
        return plugin->GetValue(variable, value);
    }

    return NPERR_NO_ERROR;
}

NPError NpapiPluginModule::NPP_SetValue(NPP instance, NPNVariable variable, void *value)
{
    FBLOG_TRACE("NPAPI", (void*) instance << "variable: " << (void*)variable);
    if (!validInstance(instance)) {
        return NPERR_INVALID_INSTANCE_ERROR;
    }

    if (NpapiPluginPtr plugin = getPlugin(instance)) {
        return plugin->SetValue(variable, value);
    } else {
        return NPERR_GENERIC_ERROR;
    }
}

max_stars_repo_name: Faham/bric-n-brac

max_issues_repo_name: Faham/bric-n-brac

max_forks_repo_name: Faham/bric-n-brac


text:
/* ============================================================================

   recurrences_ntl.cpp:  recurrences solved via NTL polynomial arithmetic

   This file is part of hypellfrob (version 2.1.1).

   Copyright (C) 2007, 2008, David Harvey

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License along
   with this program; if not, write to the Free Software Foundation, Inc.,
   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

============================================================================ */

#include <NTL/ZZ_pX.h>
#include <NTL/mat_ZZ_p.h>
#include <NTL/lzz_pX.h>
#include <NTL/mat_lzz_p.h>
#include <cassert>
#include "recurrences_ntl.h"

NTL_CLIENT

namespace hypellfrob {

/* ============================================================================

   Some template stuff

The matrix evaluation code is templated so that it can work over either ZZ_p
or zz_p. There are several template parameters, which can have two settings:

   SCALAR:       ZZ_p           zz_p
   POLY:         ZZ_pX          zz_pX
   VECTOR:       vec_ZZ_p       vec_zz_p
   MATRIX:       mat_ZZ_p       mat_zz_p
   POLYMODULUS:  ZZ_pXModulus   zz_pXModulus
   FFTREP:       FFTRep         fftRep

For the most part NTL uses the same function names for both columns, which
makes life easy, but there's a few I needed to define explicitly:

   to_scalar() convert a ZZ or int into a SCALAR
   forward_fft() runs a forward FFT (either ToFFTRep() or TofftRep())
   inverse_fft() runs an inverse FFT (either FromFFTRep() or FromfftRep())

============================================================================ */

// to_scalar(ZZ)

template <typename SCALAR> SCALAR to_scalar(const ZZ& input);

template<> inline ZZ_p to_scalar<ZZ_p>(const ZZ& input)
{
   return to_ZZ_p(input);
}

template<> inline zz_p to_scalar<zz_p>(const ZZ& input)
{
   return to_zz_p(input);
}

// to_scalar(int)

template <typename SCALAR> SCALAR to_scalar(int input);

template<> inline ZZ_p to_scalar<ZZ_p>(int input)
{
   return to_ZZ_p(input);
}

template<> inline zz_p to_scalar<zz_p>(int input)
{
   return to_zz_p(input);
}

// forward_fft

template <typename POLY, typename FFTREP>
void forward_fft(FFTREP& y, const POLY& x, long k, long lo, long hi);

template<> inline void
forward_fft<ZZ_pX, FFTRep>(FFTRep& y, const ZZ_pX& x, long k, long lo, long hi)
{
   ToFFTRep(y, x, k, lo, hi);
}

template<> inline void
forward_fft<zz_pX, fftRep>(fftRep& y, const zz_pX& x, long k, long lo, long hi)
{
   TofftRep(y, x, k, lo, hi);
}

// inverse_fft

template <typename POLY, typename FFTREP>
void inverse_fft(POLY& x, FFTREP& y, long lo, long hi);

template<> inline void
inverse_fft<ZZ_pX, FFTRep>(ZZ_pX& x, FFTRep& y, long lo, long hi)
{
   FromFFTRep(x, y, lo, hi);
}

template<> inline void
inverse_fft<zz_pX, fftRep>(zz_pX& x, fftRep& y, long lo, long hi)
{
   FromfftRep(x, y, lo, hi);
}

/* ============================================================================

   Dyadic evaluation stuff

This section essentially implements Theorem 8 of [BGS], for the particular
case of the parameters that are used in Theorem 15.

============================================================================ */

/*
Assume that f has degree d = 2^n and that g has degree 2d.

This function computes a polynomial h whose x^d through x^{2d} coefficients
(inclusive) are the same as those of f*g. The bottom d coefficients of h
will be junk.

The parameter g_fft should be the precomputed length 2d FFT of g.

The algorithm in this function is based on the paper "The middle product
algorithm" by Hanrot, Quercia and Zimmermann. (Many thanks to Victor Shoup
for writing such wonderfully modular FFT code. It really made my day.)

*/
template <typename SCALAR, typename POLY, typename FFTREP>
void middle_product(POLY& h, const POLY& f,
                    const POLY& g, const FFTREP& g_fft, int n)
{
   int d = 1 << n;
   h.rep.SetLength(2*d + 1);

   FFTREP f_fft(INIT_SIZE, n+1);

   // Compute length 2d cyclic convolutions of f and g, letting the top
   // third of f*g "wrap around" to the bottom half of the output.
   forward_fft<POLY, FFTREP>(f_fft, f, n+1, 0, 2*d);
   mul(f_fft, f_fft, g_fft);
   inverse_fft<POLY, FFTREP>(h, f_fft, 0, 2*d);

   // Need to correct for the x^{2d} term of g which got wrapped around to the
   // constant term.
   h.rep[d] -= g.rep[2*d] * f.rep[d];

   // Now h contains terms x^d through x^{2d-1} of f*g.
   // To finish off, we just need the x^{2d} term.
   SCALAR temp;
   SCALAR& sum = h.rep[2*d];
   sum = 0;
   for (int i = 0; i <= d; i++)
   {
      mul(temp, f.rep[i], g.rep[2*d-i]);
      add(sum, sum, temp);
   }
}

/*
This struct stores precomputed information that can then be used to shift
evaluation values of a polynomial F(x) of degree d = 2^n.

Specifically, given the values
  F(0), F(b), F(2*b), ..., F(d*b),
the shift() method computes
  F(a), F(a + b), F(a + 2*b), ..., F(a + d*b).

PRECONDITIONS:
   n >= 1
   1, 2, ..., d + 1 are invertible
   a + i*b are invertible for -d <= i <= d

*/
template <typename SCALAR, typename POLY, typename VECTOR, typename FFTREP>
struct DyadicShifter
{
   int d, n;

   // input_twist is a vector of length d/2 + 1.
   // The i-th entry is \prod_{0 <= j <= d, j != i} (i-j)^(-1).
   VECTOR input_twist;

   // output_twist is a vector of length d+1.
   // The i-th entry is b^(-d) \prod_{0 <= j <= d} (a + (i-j)*b).
   VECTOR output_twist;

   // kernel is a polynomial of degree 2d.
   // The coefficients are (a + k*b)^(-1) for -d <= k <= d.
   // We also store its length 2d FFT.
   POLY kernel;
   FFTREP kernel_fft;

   // Polynomials for scratch space in shift()
   POLY scratch, scratch2;

   // Constructor (performs various precomputations)
   DyadicShifter(int n, const SCALAR& a, const SCALAR& b)
   {
      assert(n >= 1);
      this->n = n;
      d = 1 << n;

      // ------------------------ compute input_twist -------------------------

      input_twist.SetLength(d/2 + 1);

      // prod = (d!)^(-1)
      SCALAR prod;
      prod = 1;
      for (int i = 2; i <= d; i++)
         mul(prod, prod, i);
      prod = 1 / prod;

      // input_twist[i] = ((d-i)!)^(-1)
      input_twist[0] = prod;
      for (int i = 1; i <= d/2; i++)
         mul(input_twist[i], input_twist[i-1], d-(i-1));

      // input_twist[i] = ((d-i)!*i!)^(-1)
      prod = input_twist[d/2];
      for (int i = d/2; i >= 0; i--)
      {
         mul(input_twist[i], input_twist[i], prod);
         mul(prod, prod, i);
      }

      // input_twist[i] = \prod_{0 <= j <= d, j != i} (i-j)^(-1)    :-)
      for (int i = 1; i <= d/2; i += 2)
         NTL::negate(input_twist[i], input_twist[i]);

      // ----------------- compute output_twist and kernel --------------------

      output_twist.SetLength(d+1);

      // c[i] = c_i = a + (i-d)*b     for 0 <= i <= 2d
      VECTOR c;
      c.SetLength(2*d+1);
      c[0] = a - d*b;
      for (int i = 1; i <= 2*d; i++)
         add(c[i], c[i-1], b);

      // accum[i] = c_0 * c_1 * ... * c_i    for 0 <= i <= 2d
      VECTOR accum;
      accum.SetLength(2*d+1);
      accum[0] = c[0];
      for (int i = 1; i <= 2*d; i++)
         mul(accum[i], accum[i-1], c[i]);

      // accum_inv[i] = (c_0 * c_1 * ... * c_i)^(-1)    for 0 <= i <= 2d
      VECTOR accum_inv;
      accum_inv.SetLength(2*d+1);
      accum_inv[2*d] = 1 / accum[2*d];
      for (int i = 2*d-1; i >= 0; i--)
         mul(accum_inv[i], accum_inv[i+1], c[i+1]);

      // kernel[i] = (c_i)^(-1)    for 0 <= i <= 2d
      kernel.rep.SetLength(2*d+1);
      kernel.rep[0] = accum_inv[0];
      for (int i = 1; i <= 2*d; i++)
         mul(kernel.rep[i], accum_inv[i], accum[i-1]);

      // precompute transform of kernel
      forward_fft<POLY, FFTREP>(kernel_fft, kernel, n+1, 0, 2*d);

      // output_twist[i] = b^{-d} * c_i * c_{i+1} * ... * c_{i+d}
      // for 0 <= i <= d
      SCALAR factor = power(b, -d);
      SCALAR temp;
      output_twist.SetLength(d+1);
      output_twist[0] = factor * accum[d];
      for (int i = 1; i <= d; i++)
      {
         mul(temp, factor, accum[i+d]);
         mul(output_twist[i], temp, accum_inv[i-1]);
      }
   }

   // Shifts evaluation values as described above.
   // Assumes both output and input have length d + 1.
   void shift(VECTOR& output, const VECTOR& input)
   {
      assert(input.length() == d+1);
      assert(output.length() == d+1);

      // multiply inputs pointwise by input_twist
      scratch.rep.SetLength(d+1);
      for (int i = 0; i <= d/2; i++)
         mul(scratch.rep[i], input[i], input_twist[i]);
      for (int i = 1; i <= d/2; i++)
         mul(scratch.rep[i+d/2], input[i+d/2], input_twist[d/2-i]);

      middle_product<SCALAR, POLY, FFTREP>(scratch2, scratch, kernel,
                                           kernel_fft, n);

      // multiply outputs pointwise by output_twist
      for (int i = 0; i <= d; i++)
         mul(output[i], scratch2.rep[i+d], output_twist[i]);
   }
};

/*
Let M0 and M1 be square matrices of size n*n. Let M(x) = M0 + x*M1; this is a
matrix of linear polys in x. Let P(x) = M(x+1) M(x+2) ... M(x+2^s); this is a
matrix of polynomials of degree 2^s. This function computes the values
    P(a), P(a + 2^t), P(a + 2*2^t), ..., P(a + 2^s*2^t).

The output array should have length n^2. Each entry should be a vector of
length 2^s+1, pre-initialised to all zeroes. The (y*n + x)-th vector will be
the values of the (y, x) entries of the above list of matrices. (This data
format is optimised for the case that 2^s+1 is much larger than n.)

PRECONDITIONS:
   0 <= s <= t
   2, 3, ..., 2^t + 1 must be invertible

*/
template <typename SCALAR, typename POLY, typename VECTOR,
          typename MATRIX, typename FFTREP>
void dyadic_evaluation(vector<VECTOR>& output,
                       const MATRIX& M0, const MATRIX& M1,
                       int s, int t, const SCALAR& a)
{
   int n = M0.NumRows();

   // base cases; just evaluate naively
   if (s <= 1)
   {
      MATRIX X[3];

      if (s == 0)
      {
         X[0] = M0 + (a+1) * M1;
         X[1] = M0 + (a+1 + (1 << t)) * M1;
      }
      else
      {
         for (int i = 0; i <= 2; i++)
            X[i] = (M0 + (a+1 + (i << t)) * M1) * (M0 + (a+2 + (i << t)) * M1);
      }

      for (int x = 0; x < n; x++)
      for (int y = 0; y < n; y++)
      for (int i = 0; i < output[0].length(); i++)
         output[y*n + x][i] = X[i][y][x];

      return;
   }

   // General case.
   // Let Q(x) = M(x+1) M(x+2) ... M(x+2^(s-1)).

   // Recursively compute Q(a), Q(a + 2^t), ..., Q(a + 2^(s-1)*2^t).
   vector<VECTOR> X(n*n);
   for (int i = 0; i < n*n; i++)
      X[i].SetLength((1 << (s-1)) + 1);
   dyadic_evaluation<SCALAR, POLY, VECTOR, MATRIX, FFTREP>
                           (X, M0, M1, s-1, t, a);

   // Do precomputations for shifting by 2^(s-1) and by (2^(s-1)+1)*2^t
   SCALAR c, b;
   c = 1 << (s-1);
   b = 1 << t;
   DyadicShifter<SCALAR, POLY, VECTOR, FFTREP> shifter1(s-1, c, b);
   DyadicShifter<SCALAR, POLY, VECTOR, FFTREP> shifter2(s-1, (c + 1) * b, b);

   // Shift by 2^(s-1) to obtain
   // Q(a + 2^(s-1)), Q(a + 2^t + 2^(s-1)), ..., Q(a + 2^(s-1)*2^t + 2^(s-1))
   vector<VECTOR> Y(n*n);
   for (int i = 0; i < n*n; i++)
   {
      Y[i].SetLength((1 << (s-1)) + 1);
      shifter1.shift(Y[i], X[i]);
   }

   // Multiply matrices to obtain
   // P(a), P(a + 2^t), ..., P(a + 2^(s-1)*2^t).
   SCALAR temp;
   for (int i = 0; i <= (1 << (s-1)); i++)
   for (int x = 0; x < n; x++)
   for (int y = 0; y < n; y++)
   for (int z = 0; z < n; z++)
   {
      mul(temp, X[y*n + z][i], Y[z*n + x][i]);
      output[y*n + x][i] += temp;
   }

   // Shift original sequence by (2^(s-1)+1)*2^t to obtain
   // Q(a + (2^(s-1)+1)*2^t), Q(a + (2^(s-1)+2)*2^t), ..., Q(a + (2^s+1)*2^t).
   for (int i = 0; i < n*n; i++)
      shifter2.shift(Y[i], X[i]);

   // Shift again by 2^(s-1) to obtain
   // Q(a + (2^(s-1)+1)*2^t + 2^(s-1)), Q(a + (2^(s-1)+2)*2^t + 2^(s-1)), ...,
   //                                             Q(a + (2^s+1)*2^t + 2^(s-1)).
   for (int i = 0; i < n*n; i++)
      shifter1.shift(X[i], Y[i]);

   // Multiply matrices to obtain
   // P(a + (2^(s-1)+1)*2^t), P(a + (2^(s-1)+2)*2^t), ..., P(a + (2^s+1)*2^t).
   // (we throw out the last one since it's surplus to requirements)
   for (int i = 0; i < (1 << (s-1)); i++)
   for (int x = 0; x < n; x++)
   for (int y = 0; y < n; y++)
   for (int z = 0; z < n; z++)
   {
      mul(temp, Y[y*n + z][i], X[z*n + x][i]);
      output[y*n + x][i + (1 << (s-1)) + 1] += temp;
   }
}

/* ============================================================================

   General evaluation stuff

This section essentially implements Corollary 10 of [BGS].

============================================================================ */

/*
This struct stores the product tree associated to a vector a[0], ..., a[n-1].

The top node stores the polynomial product
     (x - a[0]) ... (x - a[n-1]).
The two children nodes store
     (x - a[0]) ... (x - a[m-1])
and
     (x - a[m]) ... (x - a[n-1])
where m = floor(n/2). This continues recursively until we reach n = 1,
in which case just the polynomial x - a[0] is stored, and no children.

*/
template <typename SCALAR, typename POLY, typename VECTOR>
struct ProductTree
{
   // polynomial product stored at this node
   POLY poly;

   // children for left and right halves, if deg(poly) > 1
   ProductTree* child1;
   ProductTree* child2;

   // These are temp polys used by the Evaluator and Interpolator classes.
   // It's not very hygienic to keep them here... but it makes things more
   // efficient, because we need two temps for each node, and this prevent
   // unnecessary reallocations. (The lengths will be the same on repeated
   // calls to evaluate() and interpolate().)
   POLY scratch1, scratch2;

   // Constructs product tree for the supplied vector.
   ProductTree(const VECTOR& points)
   {
      build(points, 0, points.length());
   }

   ProductTree(const VECTOR& points, int start, int end)
   {
      build(points, start, end);
   }

   // Constructs product tree recursively for the subset [start, end) of
   // the supplied vector.
   void build(const VECTOR& points, int start, int end)
   {
      assert(end - start >= 1);
      assert(start >= 0);
      assert(end <= points.length());

      if (end - start == 1)
      {
         SetCoeff(poly, 1, 1);
         SetCoeff(poly, 0, -points[start]);
      }
      else
      {
         int m = (end - start) / 2;
         child1 = new ProductTree(points, start, start + m);
         child2 = new ProductTree(points, start + m, end);
         mul(poly, child1->poly, child2->poly);
      }
   }

   ~ProductTree()
   {
      if (deg(poly) > 1)
      {
         delete child1;
         delete child2;
      }
   }
};

/*
Given a list of evaluation points a[0], ..., a[n-1], this struct stores some
precomputed information to permit evaluating an arbitrary polynomial at those
points.
*/
template <typename SCALAR, typename POLY,
          typename POLYMODULUS, typename VECTOR>
struct Evaluator
{
   // The product tree for the evaluation points
   ProductTree<SCALAR, POLY, VECTOR>* tree;

   // A list of NTL ZZ_pXModulus/zz_pXModulus objects corresponding to the
   // polynomials in the product tree, in the order that they get used as the
   // tree is traversed in recursive_evaluate().
   vector<POLYMODULUS> moduli;

   // Constructs evaluator object for the given list of evaluation points
   Evaluator(const VECTOR& points)
   {
      assert(points.length() >= 1);
      tree = new ProductTree<SCALAR, POLY, VECTOR>(points);
      moduli.reserve(2*points.length());
      build(tree);
      assert(moduli.size() <= static_cast<size_t>(2*points.length()));
   }

   // Compute modulus objects for each polynomial under the supplied node of
   // the product tree; appends them in traversal order to "moduli".
   void build(const ProductTree<SCALAR, POLY, VECTOR>* node)
   {
      if (deg(node->poly) > 1)
      {
         moduli.push_back(POLYMODULUS(node->poly));
         build(node->child1);
         build(node->child2);
      }
   }

   ~Evaluator()
   {
      delete tree;
   }

   // Evaluates the input polynomial at the evaluation points, writes the
   // results to output. The output array must have the correct length.
   void evaluate(VECTOR& output, const POLY& input)
   {
      recursive_evaluate(output, input, tree, 0, 0);
   }

   // Evaluates the input polynomial at the subset [start, end) of the
   // evaluation points, which should correspond to the supplied product tree
   // node. (The length of the interval is implied by the degree of the poly
   // at that node.) Writes the output to the subset [start, end) of the
   // output array. The index parameter indicates which modulus in "moduli"
   // to use for this node of the tree. The return value is the index for
   // the modulus that should be used immediately after this call.
   int recursive_evaluate(VECTOR& output, const POLY& input,
                          ProductTree<SCALAR, POLY, VECTOR>* node,
                          int start, int index)
   {
      if (deg(node->poly) == 1)
      {
         eval(output[start], input, -coeff(node->poly, 0));
      }
      else
      {
         rem(node->scratch1, input, moduli[index++]);
         index = recursive_evaluate(output, node->scratch1, node->child1,
                                    start, index);
         index = recursive_evaluate(output, node->scratch1, node->child2,
                                    start + deg(node->child1->poly), index);
      }
      return index;
   }
};

/*
Given an integer L >= 1, this struct does some precomputations to permit
interpolating a polynomial whose values at 0, 1, ..., L are known.

PRECONDITIONS:
   1, 2, ..., L must be invertible.

*/
template <typename SCALAR, typename POLY, typename VECTOR>
struct Interpolator
{
   ProductTree<SCALAR, POLY, VECTOR>* tree;
   int L;

   // input_twist is a vector of length L+1.
   // The i-th entry is \prod_{0 <= j <= L, j != i} (i-j)^(-1).
   VECTOR input_twist;

   // vector of length L+1, used in interpolate()
   VECTOR temp;

   // Performs various precomputations for the given L.
   Interpolator(int L)
   {
      this->L = L;
      temp.SetLength(L+1);

      // Build a product tree for the evaluation points
      for (int i = 0; i <= L; i++)
         temp[i] = i;
      tree = new ProductTree<SCALAR, POLY, VECTOR>(temp);

      // prod = (L!)^(-1)
      SCALAR prod;
      prod = 1;
      for (int i = 2; i <= L; i++)
         mul(prod, prod, i);
      prod = 1 / prod;

      // input_twist[i] = (i!)^(-1),   0 <= i <= L
      input_twist.SetLength(L+1);
      input_twist[L] = prod;
      for (int i = L; i >= 1; i--)
         mul(input_twist[i-1], input_twist[i], i);

      // input_twist[i] = \prod_{0 <= j <= L, j != i} (i-j)^(-1).
      for (int i = 0; i <= L/2; i++)
      {
         mul(input_twist[i], input_twist[i], input_twist[L-i]);
         input_twist[L-i] = input_twist[i];
      }
      for (int i = L-1; i >= 0; i -= 2)
         NTL::negate(input_twist[i], input_twist[i]);
   }

   ~Interpolator()
   {
      delete tree;
   }

   // Returns the polynomial
   //    \sum_{i=start}^{end-1} values[i] * (x-start) (x-start+1) ... (x-end-1)
   // where [start, end) is the interval associated to the supplied product
   // tree node, and where the (x-i) term is omitted in each product.
   void combine(POLY& output, const VECTOR& values,
                ProductTree<SCALAR, POLY, VECTOR>* node, int start)
   {
      if (deg(node->poly) == 1)
      {
         // base case
         clear(output);
         SetCoeff(output, 0, values[start]);
      }
      else
      {
         // recursively build up from two halves
         // i.e. if f1, f2 are the results of "combine" for the two halves,
         // and if p1, p2 are the associated product tree polys, we compute
         // f1*p2 + f2*p1

         combine(node->scratch1, values, node->child1, start);
         mul(output, node->scratch1, node->child2->poly);

         combine(node->scratch1, values, node->child2,
                 start + deg(node->child1->poly));
         mul(node->scratch2, node->scratch1, node->child1->poly);

         add(output, output, node->scratch2);
      }
   }

   // Returns a polynomial F(x) of degree at most L such that F(i) = values[i]
   // for each 0 <= i <= L.
   void interpolate(POLY& output, const VECTOR& values)
   {
      assert(values.length() == L+1);

      // multiply input values pointwise by input_twist; this corrects for
      // the factor (i-0) (i-1) ... (i-L) (where the i-i factor is omitted).
      for (int i = 0; i <= L; i++)
         mul(temp[i], values[i], input_twist[i]);

      // do the interpolation
      combine(output, temp, tree, 0);
   }
};

/* ============================================================================

   Matrix products over arbitrary, relatively short intervals

This section implements something similar to steps 1, 2, ... and the final
refining step of Theorem 15 of [BGS].

============================================================================ */

/*
Let M0 and M1 be matrices of constants. This function evaluates
   M(x) = M0 + x*M1
at x = a.

The output matrix must already have the correct dimensions.

*/
template <typename SCALAR, typename MATRIX>
void eval_matrix(MATRIX& output, const MATRIX& M0, const MATRIX& M1,
                 const SCALAR& a)
{
   int n = M0.NumRows();
   for (int x = 0; x < n; x++)
   for (int y = 0; y < n; y++)
   {
      mul(output[x][y], a, M1[x][y]);
      add(output[x][y], output[x][y], M0[x][y]);
   }
}

/*
Similar to ntl_interval_products. This is used as a subroutine of
ntl_interval_products() to handle the smaller "refining" subintervals.
Its asymptotic complexity theoretically has an extra logarithmic factor
over that of ntl_interval_products().

PRECONDITIONS:
   Let d = sum of lengths of intervals. Then 2, 3, ... 1 + floor(sqrt(d)) must
   all be invertible.

*/
template <typename SCALAR, typename POLY, typename POLYMODULUS,
          typename VECTOR, typename MATRIX>
void ntl_short_interval_products(vector<MATRIX>& output,
                                 const MATRIX& M0, const MATRIX& M1,
                                 const vector<ZZ>& target)
{
   output.clear();

   if (target.size() == 0)
      return;

   int dim = M0.NumRows();
   int num_intervals = target.size() / 2;

   // Determine maximum target interval length
   int max_length = -1;
   for (size_t i = 0; i < target.size(); i += 2)
   {
      int temp = to_ulong(target[i+1] - target[i]);
      if (temp > max_length)
         max_length = temp;
   }

   // Select an appropriate length for the matrix products we'll use
   int L, max_eval_points;
   if (max_length > 2*num_intervals)
   {
      // The intervals are still pretty long relative to the number of
      // intervals, so we're only going to do a single multipoint
      // evaluation.
      L = 1 + to_ulong(SqrRoot(num_intervals * to_ZZ(max_length)));
      max_eval_points = L;
   }
   else
   {
      // The intervals are getting pretty short, so we probably will need
      // to do several shorter multipoint evaluations.
      L = 1 + max_length/2;
      max_eval_points = num_intervals;
   }

   // =========================================================================
   // Step 1: compute entries of M(X, X+L) as polynomials in X.

   vector<POLY> polys(dim*dim);
   {
      // left_accum[i]  = M(L-i-1, L)  for 0 <= i <= L-1
      // right_accum[i] = M(L, L+i+1)  for 0 <= i <= L-1
      vector<MATRIX> left_accum(L), right_accum(L);

      MATRIX temp;
      temp.SetDims(dim, dim);

      left_accum[0].SetDims(dim, dim);
      eval_matrix<SCALAR, MATRIX>(left_accum[0], M0, M1, to_scalar<SCALAR>(L));
      for (int i = L-1; i >= 1; i--)
      {
         eval_matrix<SCALAR, MATRIX>(temp, M0, M1, to_scalar<SCALAR>(i));
         mul(left_accum[L-i], temp, left_accum[L-i-1]);
      }

      right_accum[0].SetDims(dim, dim);
      eval_matrix<SCALAR, MATRIX>(right_accum[0], M0, M1,
                                  to_scalar<SCALAR>(L+1));
      for (int i = 1; i <= L-1; i++)
      {
         eval_matrix<SCALAR, MATRIX>(temp, M0, M1, to_scalar<SCALAR>(L+1+i));
         mul(right_accum[i], right_accum[i-1], temp);
      }

      // Use left_accum and right_accum to compute:
      // initial[i] = M(i, L+i)    for 0 <= i <= L
      // i.e. initial[i] are the values of M(X, X+L) at X = 0, 1, ..., L.
      vector<MATRIX> initial(L+1);
      initial[0] = left_accum.back();
      initial[L] = right_accum.back();
      for (int i = 1; i <= L-1; i++)
         mul(initial[i], left_accum[L-1-i], right_accum[i-1]);

      // Now interpolate entries of initial[i] to get entries of M(X, X+L)
      // as polynomials of degree L.
      Interpolator<SCALAR, POLY, VECTOR> interpolator(L);
      VECTOR values;
      values.SetLength(L+1);
      for (int x = 0; x < dim; x++)
      for (int y = 0; y < dim; y++)
      {
         for (int j = 0; j <= L; j++)
            values[j] = initial[j][y][x];
         interpolator.interpolate(polys[y*dim + x], values);
      }
   }

   // =========================================================================
   // Step 2: decompose intervals into subintervals of length L which we'll
   // attack by direct multipoint evaluation, plus leftover pieces that we'll
   // handle with a recursive call to ntl_short_interval_products().

   // eval_points holds all the values of X for which we want to
   // evaluate M(X, X+L)
   VECTOR eval_points;
   eval_points.SetMaxLength(max_eval_points);

   // leftover_target is the list of leftover intervals that we're going to
   // later do recursively
   vector<ZZ> leftover_target;
   leftover_target.reserve(target.size());

   ZZ current, next;
   for (size_t i = 0; i < target.size(); i += 2)
   {
      current = target[i];
      next = current + L;
      while (next <= target[i+1])
      {
         // [current, next) fits inside this interval, so peel it off into
         // eval_points
         append(eval_points, to_scalar<SCALAR>(current));
         swap(current, next);
         next = current + L;
      }
      if (current < target[i+1])
      {
         // the rest of this interval is too short to handle with M(X, X+L),
         // so put it in the leftover bin
         leftover_target.push_back(current);
         leftover_target.push_back(target[i+1]);
      }
   }

   // =========================================================================
   // Step 3: recursively handle leftover pieces

   // leftover_matrices[i] holds the matrix for leftover interval #i
   vector<MATRIX> leftover_matrices;
   ntl_short_interval_products<SCALAR, POLY, POLYMODULUS, VECTOR, MATRIX>
                    (leftover_matrices, M0, M1, leftover_target);

   // =========================================================================
   // Step 4: evaluate M(X, X+L) at each of the evaluation points. We do this
   // by breaking up the list of evaluation points into blocks of length at
   // most L+1, and using multipoint evaluation on each block.

   // main_matrices[i] will hold M(X, X+L) for the i-th evaluation point X.
   vector<MATRIX> main_matrices(eval_points.length());
   for (size_t i = 0; i < main_matrices.size(); i++)
      main_matrices[i].SetDims(dim, dim);

   VECTOR block, values;
   block.SetMaxLength(L+1);
   values.SetMaxLength(L+1);

   // for each block...
   for (int i = 0; i < eval_points.length(); i += (L+1))
   {
      // determine length of this block, which is at most L+1
      int length = eval_points.length() - i;
      if (length >= (L+1))
         length = (L+1);
      block.SetLength(length);

      // construct Evaluator object for evaluating at these points
      for (int j = 0; j < length; j++)
         block[j] = eval_points[i+j];
      Evaluator<SCALAR, POLY, POLYMODULUS, VECTOR> evaluator(block);

      // evaluate each entry of M(X, X+L) at those points
      for (int x = 0; x < dim; x++)
      for (int y = 0; y < dim; y++)
      {
         evaluator.evaluate(values, polys[y*dim + x]);
         for (int k = 0; k < length; k++)
            main_matrices[i+k][y][x] = values[k];
      }
   }

   // =========================================================================
   // Step 5: merge together the matrices obtained from the multipoint
   // evaluation step and the recursive leftover interval step.

   output.clear();
   output.resize(target.size() / 2);
   for (size_t i = 0; i < target.size()/2; i++)
      output[i].SetDims(dim, dim);

   int main_index = 0;       // index into main_matrices
   int leftover_index = 0;   // index into leftover_matrices

   MATRIX temp;
   temp.SetDims(dim, dim);

   for (size_t i = 0; i < target.size(); i += 2)
   {
      current = target[i];
      next = current + L;
      ident(output[i/2], dim);

      while (next <= target[i+1])
      {
         // merge in a matrix from multipoint evaluation step
         mul(temp, output[i/2], main_matrices[main_index++]);
         swap(temp, output[i/2]);
         swap(current, next);
         next = current + L;
      }
      if (current < target[i+1])
      {
         // merge in a matrix from leftover interval step
         mul(temp, output[i/2], leftover_matrices[leftover_index++]);
         swap(temp, output[i/2]);
      }
   }
}

/* ============================================================================

   Matrix products over arbitrary, long intervals

This section implements an algorithm similar to Theorem 15 of [BGS].

============================================================================ */

/*
See interval_products_wrapper().

NOTE:
   This algorithm works best if the intervals are very long and don't have
   much space between them. The case where the gaps are relatively large is
   best handled by ntl_short_interval_products().

*/
template <typename SCALAR, typename POLY, typename POLYMODULUS,
          typename VECTOR, typename MATRIX, typename FFTREP>
void ntl_interval_products(vector<MATRIX>& output,
                           const MATRIX& M0, const MATRIX& M1,
                           const vector<ZZ>& target)
{
   assert(target.size() % 2 == 0);
   output.resize(target.size() / 2);

   int dim = M0.NumRows();
   assert(dim == M0.NumCols());
   assert(dim == M1.NumRows());
   assert(dim == M1.NumCols());

   // =========================================================================
   // Step 0: get as many intervals as possible using dyadic_evaluation().

   // step0_matrix[i] is the transition matrix between step0_index[2*i]
   // and step0_index[2*i+1].
   vector<MATRIX> step0_matrix;
   vector<ZZ> step0_index;
   // preallocate the maximum number of matrices that could arise (plus safety)
   size_t reserve_size = target.size() +
                      4*NumBits(target.back() - target.front());
   step0_matrix.reserve(reserve_size);
   step0_index.reserve(2 * reserve_size);

   ZZ current_index = target.front();
   size_t next_target = 0;   // index into "target" array

   // This flag indicates whether the last entry of step0_matrix is
   // still accumulating matrices (in which the right endpoint of the
   // corresponding interval hasn't been written to step0_index yet).
   int active = 0;

   MATRIX temp_mat;
   temp_mat.SetDims(dim, dim);

   while (current_index < target.back() - 3)
   {
      // find largest t such that 2^t*(2^t + 1) <= remaining distance to go
      ZZ remaining = target.back() - current_index;
      int t = 0;
      while ((to_ZZ(1) << (2*t)) + (1 << t) <= remaining)
         t++;
      t--;

      // evaluate matrices for 2^t+1 intervals of length 2^t
      vector<VECTOR> dyadic_output(dim*dim);
      for (int i = 0; i < dim*dim; i++)
         dyadic_output[i].SetLength((1 << t) + 1);
      dyadic_evaluation<SCALAR, POLY, VECTOR, MATRIX, FFTREP>
               (dyadic_output, M0, M1, t, t, to_scalar<SCALAR>(current_index));

      // Walk through the intervals we just computed. Find maximal subsequences
      // of intervals none of which contain any target endpoints. Merge them
      // together (by multiplying the appropriate matrices) and store results
      // in step0_matrix, step0_index.

      SCALAR scratch;

      for (int i = 0; i <= (1 << t); i++, current_index += (1 << t))
      {
         assert(next_target == target.size() ||
                target[next_target] >= current_index);

         // Skip over target endpoints which are exactly at the beginning
         // of this interval
         while ((next_target < target.size()) &&
                (target[next_target] == current_index))
         {
            // if there's an active matrix, don't forget to close it off
            if (active)
            {
               step0_index.push_back(current_index);
               active = 0;
            }
            next_target++;
         }

         // Test if any target endpoints are strictly within this interval.
         if ((next_target == target.size()) ||
             (target[next_target] >= current_index + (1 << t)))
         {
            // There are no target endpoints in this interval.
            if (active)
            {
               // Merge this matrix with the active one
               MATRIX& active_mat = step0_matrix.back();
               for (int y = 0; y < dim; y++)
               for (int x = 0; x < dim; x++)
               {
                  SCALAR& accum = temp_mat[y][x];
                  accum = 0;
                  for (int z = 0; z < dim; z++)
                  {
                     mul(scratch, active_mat[y][z],
                                  dyadic_output[z*dim + x][i]);
                     add(accum, accum, scratch);
                  }
               }

               swap(temp_mat, active_mat);
            }
            else
            {
               // Make this matrix into a new active one
               step0_index.push_back(current_index);
               step0_matrix.resize(step0_matrix.size() + 1);
               MATRIX& X = step0_matrix.back();
               X.SetDims(dim, dim);
               for (int y = 0; y < dim; y++)
               for (int x = 0; x < dim; x++)
                  X[y][x] = dyadic_output[y*dim + x][i];
               active = 1;
            }
         }
         else
         {
            // There are target endpoints in this interval.
            if (active)
            {
               // If there is still an active matrix, close it off.
               step0_index.push_back(current_index);
               active = 0;
            }

            // skip over any other endpoints in this interval
            while ((next_target < target.size()) &&
                   (target[next_target] < current_index + (1 << t)))
            {
               next_target++;
            }
         }
      }
   }

   // If there is still an active matrix, close it off.
   if (active)
      step0_index.push_back(current_index);

   assert(step0_index.size() == 2*step0_matrix.size());

   // =========================================================================
   // Step 1: Make a list of all subintervals that we are going to need in
   // the refining steps.

   size_t next_step0 = 0;        // index into step0_index
   vector<ZZ> step1_index;    // list of pairs of endpoints of needed intervals
   step1_index.reserve(2*target.size());

   // add sentinel endpoints to make the next loop simpler:
   step0_index.push_back(target.back() + 10);
   step0_index.push_back(target.back() + 20);

   for (next_target = 0; next_target < target.size(); next_target += 2)
   {
      // skip dyadic intervals that come before this target interval
      while (step0_index[next_step0+1] <= target[next_target])
         next_step0 += 2;

      if (step0_index[next_step0] < target[next_target+1])
      {
         // The next dyadic interval starts before the end of this target
         // interval.
         if (step0_index[next_step0] > target[next_target])
         {
            // The next dyadic interval starts strictly within this target
            // interval, so we need to create a refining subinterval for the
            // initial segment of this target interval.
            step1_index.push_back(target[next_target]);
            step1_index.push_back(step0_index[next_step0]);
         }

         // Skip over dyadic intervals to find the last one still contained
         // within this target interval.
         while (step0_index[next_step0+3] <= target[next_target+1])
            next_step0 += 2;

         if (step0_index[next_step0+1] < target[next_target+1])
         {
            // The next dyadic interval finishes strictly within this target
            // interval, so we need to create a refining subinterval for the
            // final segment of this target interval.
            step1_index.push_back(step0_index[next_step0+1]);
            step1_index.push_back(target[next_target+1]);
         }

         // Move on to next dyadic interval
         next_step0 += 2;
      }
      else
      {
         // The next dyadic interval starts beyond (or just at the end of)
         // this target interval, so we need to create a refining subinterval
         // for this *whole* target interval.
         step1_index.push_back(target[next_target]);
         step1_index.push_back(target[next_target+1]);
      }
   }

   // remove sentinels for my sanity
   step0_index.pop_back();
   step0_index.pop_back();

   // Step 1b: Compute matrix products over those refining subintervals.
   vector<MATRIX> step1_matrix;
   ntl_short_interval_products<SCALAR, POLY, POLYMODULUS, VECTOR, MATRIX>
                                   (step1_matrix, M0, M1, step1_index);

   assert(step1_index.size() == 2 * step1_matrix.size());

   // =========================================================================
   // Step 2: Merge together the dyadic intervals and refining intervals into
   // a single list, in sorted order.

   vector<MATRIX> step2_matrix(step0_matrix.size() + step1_matrix.size());
   vector<ZZ> step2_index(step0_index.size() + step1_index.size());

   // add sentinels to make the next loop simpler
   step0_index.push_back(target.back() + 10);
   step0_index.push_back(target.back() + 20);
   step1_index.push_back(target.back() + 10);
   step1_index.push_back(target.back() + 20);

   next_step0 = 0;       // index into step0_matrix
   size_t next_step1 = 0;   // index into step1_matrix

   for (size_t next_step2 = 0; next_step2 < step2_matrix.size(); next_step2++)
   {
      if (step0_index[2*next_step0] < step1_index[2*next_step1])
      {
         // grab a matrix and pair of indices from step0
         swap(step2_matrix[next_step2], step0_matrix[next_step0]);
         step2_index[2*next_step2] = step0_index[2*next_step0];
         step2_index[2*next_step2+1] = step0_index[2*next_step0+1];
         next_step0++;
      }
      else
      {
         // grab a matrix and pair of indices from step1
         swap(step2_matrix[next_step2], step1_matrix[next_step1]);
         step2_index[2*next_step2] = step1_index[2*next_step1];
         step2_index[2*next_step2+1] = step1_index[2*next_step1+1];
         next_step1++;
      }
   }

   // remove sentinels for my sanity
   step0_index.pop_back();
   step0_index.pop_back();
   step1_index.pop_back();
   step1_index.pop_back();

   assert(step2_index.size() == 2*step2_matrix.size());

   // =========================================================================
   // Step 3: Walk through target intervals, and merge together appropriate
   // intervals from step2 to get those target intervals.

   size_t next_step2 = 0;    // index into step2_matrix

   // add sentinels to make the next loop simpler
   step2_index.push_back(target.back() + 1);
   step2_index.push_back(target.back() + 2);

   for (size_t next_target = 0; next_target < target.size(); next_target += 2)
   {
      // search for step2 interval matching the start of this target interval
      while (step2_index[2*next_step2] < target[next_target])
         next_step2++;

      assert(step2_index[2*next_step2] == target[next_target]);

      // merge together matrices for step2 intervals contained in this target
      // interval
      swap(step2_matrix[next_step2++], output[next_target/2]);
      while (step2_index[2*next_step2+1] <= target[next_target+1])
      {
         mul(temp_mat, output[next_target/2], step2_matrix[next_step2]);
         swap(temp_mat, output[next_target/2]);
         next_step2++;
      }
   }
}

// explicit instantiations for zz_p and ZZ_p versions:

template void ntl_interval_products
      <ZZ_p, ZZ_pX, ZZ_pXModulus, vec_ZZ_p, mat_ZZ_p, FFTRep>
      (vector<mat_ZZ_p>& output, const mat_ZZ_p& M0, const mat_ZZ_p& M1,
                             const vector<ZZ>& target);

template void ntl_interval_products
      <zz_p, zz_pX, zz_pXModulus, vec_zz_p, mat_zz_p, fftRep>
      (vector<mat_zz_p>& output, const mat_zz_p& M0, const mat_zz_p& M1,
                             const vector<ZZ>& target);

};    // namespace hypellfrob

// ----------------------- end of file

ext: cpp

max_forks_repo_name: dimpase/sage


text:
//  Copyright (c) 2007-2014 Hartmut Kaiser
//
//  Distributed under the Boost Software License, Version 1.0. (See accompanying
//  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#if !defined(HPX_SEGMENTED_ITERATOR_TRAITS_OCT_14_2014_0229PM)
#define HPX_SEGMENTED_ITERATOR_TRAITS_OCT_14_2014_0229PM

#include <hpx/traits.hpp>

#include <boost/mpl/bool.hpp>

namespace hpx { namespace traits
{
    ///////////////////////////////////////////////////////////////////////////
    template <typename Iterator, typename Enable>
    struct segmented_iterator_traits
    {
        typedef boost::mpl::false_ is_segmented_iterator;
    };

    ///////////////////////////////////////////////////////////////////////////
    // traits allowing to distinguish iterators which have a purely local
    // representation
    template <typename Iterator, typename Enable>
    struct segmented_local_iterator_traits
    {
        typedef boost::mpl::false_ is_segmented_local_iterator;

        typedef Iterator iterator;
        typedef Iterator local_iterator;
        typedef Iterator local_raw_iterator;

        static local_raw_iterator const& local(local_iterator const& it)
        {
            return it;
        }

        static local_iterator const& remote(local_raw_iterator const& it)
        {
            return it;
        }

        static local_raw_iterator local(local_iterator&& it)
        {
            return std::move(it);
        }

        static local_iterator remote(local_raw_iterator&& it)
        {
            return std::move(it);
        }
    };
}}

#endif

ext: hpp


text:
#include <replay/planar_direction.hpp>
#include <boost/math/constants/constants.hpp>
#include <catch2/catch.hpp>

using namespace replay;
namespace
{
    constexpr auto pi = boost::math::constants::pi<float>();
    constexpr auto two_pi = boost::math::constants::two_pi<float>();
    constexpr auto half_pi = boost::math::constants::half_pi<float>();
}

TEST_CASE("can_move_around_singularities", "[planar_direction]")
{
    planar_direction current(-4.32833433f);
    planar_direction to(1.95485115f);

    auto next = planar_direction::move(current, to, 1000.f);
    REQUIRE(next == to);
}

TEST_CASE("can_move_toward_a_target", "[planar_direction]")
{
    planar_direction current(2.f);
    planar_direction to(1.f);

    auto next = planar_direction::move(current, to, 0.5f);
    REQUIRE(next.angle() == 1.5f);
}

TEST_CASE("move_picks_the_shorter_arc", "[planar_direction]")
{
    planar_direction current(0.1f);
    planar_direction to(two_pi - 0.1f);
    auto next = planar_direction::move(current, to, 0.1f);
    REQUIRE(next.angle() == 0.f);
}

TEST_CASE("normalizes_high_values_correctly", "[planar_direction]")
{
    REQUIRE(planar_direction(boost::math::constants::two_pi<float>()).normalized().angle() == Approx(0.f));
}

TEST_CASE("normalizes_low_values_correctly", "[planar_direction]")
{
    REQUIRE(planar_direction(-1.5f * pi).normalized().angle() == Approx(pi * 0.5f));
}

TEST_CASE("can_interpolate_between_positive_values", "[planar_direction]")
{
    REQUIRE(lerp(planar_direction(0.f), planar_direction(pi*2.5f), 0.25f).angle() == Approx(pi * 0.125f));
}

max_stars_repo_name: ltjax/replay

max_issues_repo_name: ltjax/replay

max_forks_repo_name: ltjax/replay


text: #include <boost/geometry/strategies/cartesian/point_in_box.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
#include "applicationoptionsconnection.hpp"

#include <boost/log/trivial.hpp>
#include <nlohmann/json.hpp>
#include <wx/taskbarbutton.h>

#include "../applicationoptions.hpp"
#include "../ui/mainframe.hpp"

using json = nlohmann::json;
using pt::IPC::ApplicationOptionsConnection;

ApplicationOptionsConnection::ApplicationOptionsConnection(pt::UI::MainFrame* frame)
    : m_frame(frame)
{
}

bool ApplicationOptionsConnection::OnExecute(const wxString&, const void *data, size_t size, wxIPCFormat format)
{
    std::string textData = GetTextFromData(data, size, format);

    m_frame->GetEventHandler()->CallAfter([this, textData]()
        {
            json j;
            pt::CommandLineOptions options;

            try
            {
                j = json::parse(textData);
                options.files = j["files"].get<std::vector<std::string>>();
                options.magnets = j["magnet_links"].get<std::vector<std::string>>();
                options.save_path = j["save_path"].get<std::string>();
                options.silent = j["silent"].get<bool>();
            }
            catch (std::exception const& e)
            {
                BOOST_LOG_TRIVIAL(error) << "Failed to parse application options JSON: " << e.what();
                return;
            }

            m_frame->MSWGetTaskBarButton()->Show();

            if (m_frame->IsIconized())
            {
                m_frame->Restore();
            }

            m_frame->Raise();
            m_frame->Show();
            m_frame->HandleParams(options);
        });

    return true;
}

max_stars_repo_name: vktr/PicoTorrent

max_issues_repo_name: vktr/PicoTorrent

max_forks_repo_name: vktr/PicoTorrent


text:
/*
 * Copyright 2017 Fraunhofer Institute for Manufacturing Engineering and Automation (IPA)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0

 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include <ros/ros.h>
#include <geometry_msgs/Pose.h>
#include <moveit_msgs/CollisionObject.h>
#include <shape_msgs/SolidPrimitive.h>
#include <geometric_shapes/shape_operations.h>
#include <geometric_shapes/mesh_operations.h>

#include <boost/variant.hpp>
#include <boost/shared_ptr.hpp>

int main(int argc, char** argv)
{
  ros::init(argc, argv, "test_obstacle_publisher_node");
  ros::NodeHandle n;
  ros::Publisher pub = n.advertise<moveit_msgs::CollisionObject>("/collision_object", 1);
  ros::Duration(1.0).sleep();

  // primitive sphere
  moveit_msgs::CollisionObject col_sphere_msg;
  col_sphere_msg.id = "sphere_primitive";
  col_sphere_msg.header.frame_id = "odom_combined";
  col_sphere_msg.operation = moveit_msgs::CollisionObject::ADD;

  shape_msgs::SolidPrimitive sphere;
  sphere.type = shape_msgs::SolidPrimitive::SPHERE;
  sphere.dimensions.push_back(0.1);

  geometry_msgs::Pose pose;
  pose.position.x = 1.0;
  pose.position.y = 0.0;
  pose.position.z = 0.0;
  pose.orientation.w = 1.0;

  col_sphere_msg.primitives.push_back(sphere);
  col_sphere_msg.primitive_poses.push_back(pose);

  pub.publish(col_sphere_msg);
  ros::Duration(0.1).sleep();

  // mesh sphere
  moveit_msgs::CollisionObject col_sphere_mesh_msg;
  col_sphere_mesh_msg.id = "sphere_mesh";
  col_sphere_mesh_msg.header.frame_id = "odom_combined";
  col_sphere_mesh_msg.operation = moveit_msgs::CollisionObject::ADD;

  shapes::Mesh* sphere_shape = shapes::createMeshFromResource("package://cob_obstacle_distance_moveit/files/sphere.stl");
  shapes::ShapeMsg sphere_mesh_msg;
  shapes::constructMsgFromShape(sphere_shape, sphere_mesh_msg);
  shape_msgs::Mesh sphere_mesh = boost::get<shape_msgs::Mesh>(sphere_mesh_msg);

  pose.position.x = 1.0;
  // pose.position.y = 1.0;
  pose.position.z = 0.0;
  pose.orientation.w = 1.0;

  col_sphere_mesh_msg.meshes.push_back(sphere_mesh);
  col_sphere_mesh_msg.mesh_poses.push_back(pose);

  pub.publish(col_sphere_mesh_msg);
  ros::Duration(0.1).sleep();

  // primitive box
  moveit_msgs::CollisionObject col_box_msg;
  col_box_msg.id = "box_primitive";
  col_box_msg.header.frame_id = "odom_combined";
  col_box_msg.operation = moveit_msgs::CollisionObject::ADD;

  shape_msgs::SolidPrimitive box;
  box.type = shape_msgs::SolidPrimitive::BOX;
  box.dimensions.push_back(0.1);
  box.dimensions.push_back(0.1);
  box.dimensions.push_back(0.1);

  pose.position.x = 1.0;
  pose.position.y = 0.0;
  pose.position.z = 1.0;
  pose.orientation.w = 1.0;

  col_box_msg.primitives.push_back(box);
  col_box_msg.primitive_poses.push_back(pose);

  pub.publish(col_box_msg);
  ros::Duration(0.1).sleep();

  // mesh box
  moveit_msgs::CollisionObject col_box_mesh_msg;
  col_box_mesh_msg.id = "box_mesh";
  col_box_mesh_msg.header.frame_id = "odom_combined";
  col_box_mesh_msg.operation = moveit_msgs::CollisionObject::ADD;

  shapes::Mesh* box_shape = shapes::createMeshFromResource("package://cob_obstacle_distance_moveit/files/box.stl");
  shapes::ShapeMsg box_mesh_msg;
  shapes::constructMsgFromShape(box_shape, box_mesh_msg);
  shape_msgs::Mesh box_mesh = boost::get<shape_msgs::Mesh>(box_mesh_msg);

  pose.position.x = 1.0;
  // pose.position.y = 1.0;
  pose.position.z = 1.0;
  pose.orientation.w = 1.0;

  col_box_mesh_msg.meshes.push_back(box_mesh);
  col_box_mesh_msg.mesh_poses.push_back(pose);

  pub.publish(col_box_mesh_msg);
  ros::Duration(0.1).sleep();

  // primitive cylinder
  moveit_msgs::CollisionObject col_cylinder_msg;
  col_cylinder_msg.id = "cylinder_primitive";
  col_cylinder_msg.header.frame_id = "odom_combined";
  col_cylinder_msg.operation = moveit_msgs::CollisionObject::ADD;

  shape_msgs::SolidPrimitive cylinder;
  cylinder.type = shape_msgs::SolidPrimitive::CYLINDER;
  cylinder.dimensions.push_back(0.1);
  cylinder.dimensions.push_back(0.1);

  pose.position.x = 1.0;
  pose.position.y = 0.0;
  pose.position.z = 2.0;
  pose.orientation.w = 1.0;

  col_cylinder_msg.primitives.push_back(cylinder);
  col_cylinder_msg.primitive_poses.push_back(pose);

  pub.publish(col_cylinder_msg);
  ros::Duration(0.1).sleep();

  // mesh cylinder
  moveit_msgs::CollisionObject col_cylinder_mesh_msg;
  col_cylinder_mesh_msg.id = "cylinder_mesh";
  col_cylinder_mesh_msg.header.frame_id = "odom_combined";
  col_cylinder_mesh_msg.operation = moveit_msgs::CollisionObject::ADD;

  shapes::Mesh* cylinder_shape = shapes::createMeshFromResource("package://cob_obstacle_distance_moveit/files/cylinder.stl");
  shapes::ShapeMsg cylinder_mesh_msg;
  shapes::constructMsgFromShape(cylinder_shape, cylinder_mesh_msg);
  shape_msgs::Mesh cylinder_mesh = boost::get<shape_msgs::Mesh>(cylinder_mesh_msg);

  pose.position.x = 1.0;
  // pose.position.y = 1.0;
  pose.position.z = 2.0;
  pose.orientation.w = 1.0;

  col_cylinder_mesh_msg.meshes.push_back(cylinder_mesh);
  col_cylinder_mesh_msg.mesh_poses.push_back(pose);

  pub.publish(col_cylinder_mesh_msg);
  ros::Duration(0.1).sleep();

  /// test robot objects
  // test primitive
  moveit_msgs::CollisionObject test_col_prim_msg;
  test_col_prim_msg.id = "test_primitive";
  test_col_prim_msg.header.frame_id = "odom_combined";
  test_col_prim_msg.operation = moveit_msgs::CollisionObject::ADD;

  shape_msgs::SolidPrimitive test_prim;
  //test_prim.type = shape_msgs::SolidPrimitive::SPHERE;
  //test_prim.type = shape_msgs::SolidPrimitive::BOX;
  test_prim.type = shape_msgs::SolidPrimitive::CYLINDER;
  test_prim.dimensions.push_back(0.1);
  test_prim.dimensions.push_back(0.1);
  //test_prim.dimensions.push_back(0.1);

  geometry_msgs::Pose test_pose;
  test_pose.position.x = 2.0;
  test_pose.position.y = 0.0;
  test_pose.position.z = 1.0;
  test_pose.orientation.w = 1.0;

  test_col_prim_msg.primitives.push_back(test_prim);
  test_col_prim_msg.primitive_poses.push_back(test_pose);

  pub.publish(test_col_prim_msg);
  ros::Duration(0.1).sleep();

  // test mesh
  moveit_msgs::CollisionObject test_col_mesh_msg;
  test_col_mesh_msg.id = "test_mesh";
  test_col_mesh_msg.header.frame_id = "odom_combined";
  test_col_mesh_msg.operation = moveit_msgs::CollisionObject::ADD;

  //shapes::Mesh* test_shape = shapes::createMeshFromResource("package://cob_obstacle_distance_moveit/files/sphere.stl");
  //shapes::Mesh* test_shape = shapes::createMeshFromResource("package://cob_obstacle_distance_moveit/files/box.stl");
  shapes::Mesh* test_shape = shapes::createMeshFromResource("package://cob_obstacle_distance_moveit/files/cylinder.stl");
  shapes::ShapeMsg test_mesh_msg;
  shapes::constructMsgFromShape(test_shape, test_mesh_msg);
  shape_msgs::Mesh test_mesh = boost::get<shape_msgs::Mesh>(test_mesh_msg);

  test_pose.position.x = 2.0;
  // test_pose.position.y = 1.0;
  test_pose.position.z = 1.0;
  test_pose.orientation.w = 1.0;

  test_col_mesh_msg.meshes.push_back(test_mesh);
  test_col_mesh_msg.mesh_poses.push_back(test_pose);

  pub.publish(test_col_mesh_msg);
  ros::Duration(0.1).sleep();

  ROS_INFO("Done");

  return 0;
}


text:
//
// Copyright (c) 2016 CNRS
//

#include "pinocchio/algorithm/energy.hpp"
#include "pinocchio/algorithm/crba.hpp"

#include "pinocchio/parsers/sample-models.hpp"

#include <boost/test/unit_test.hpp>
#include <boost/utility/binary.hpp>
#include <boost/test/floating_point_comparison.hpp>

BOOST_AUTO_TEST_SUITE ( BOOST_TEST_MODULE )

BOOST_AUTO_TEST_CASE(test_kinetic_energy)
{
  using namespace Eigen;
  using namespace pinocchio;

  pinocchio::Model model;
  pinocchio::buildModels::humanoidRandom(model);
  pinocchio::Data data(model);

  VectorXd q = VectorXd::Zero(model.nq);
  VectorXd v = VectorXd::Ones(model.nv);
  VectorXd a = VectorXd::Ones(model.nv);

  data.M.fill(0);  crba(model,data,q);
  data.M.triangularView<Eigen::StrictlyLower>()
  = data.M.transpose().triangularView<Eigen::StrictlyLower>();

  double kinetic_energy_ref = 0.5*v.transpose() * data.M * v;
  double kinetic_energy = kineticEnergy(model, data, q, v);

  BOOST_CHECK_SMALL(kinetic_energy_ref - kinetic_energy, 1e-12);
}

BOOST_AUTO_TEST_SUITE_END()

max_stars_repo_name: andreadelprete/pinocchio

max_issues_repo_name: andreadelprete/pinocchio

max_forks_repo_name: andreadelprete/pinocchio


text:
#ifndef MTGP_CALC_JUMP_HPP
#define MTGP_CALC_JUMP_HPP
/**
 * @file mtgp-calc-jump.hpp
 *
 * @brief functions for calculating jump polynomial.
 *
 * calculate jump polynomial from jump step and the
 * characteristic polynomial using PowerXMod.

 * jump polynomial p is calculated

 *     p = Xs mod q,

 * where s is jump step and q is characteristic polynomial.
 *
 * @author Mutsuo Saito (Hiroshima University)
 * @author Makoto Matsumoto (The University of Tokyo)
 *
 * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto,
 * Hiroshima University and The University of Tokyo.
 * All rights reserved.
 *
 * The 3-clause BSD License is applied to this software, see
 * LICENSE.txt
 */
#include <iostream>
#include <iomanip>
#include <sstream>
#include <NTL/GF2X.h>

/**
 * convert polynomial to string
 * @param[out] x output string
 * @param[in] polynomial polynomial
 */
static inline void polytostring(std::string& x, NTL::GF2X& polynomial)
{
    using namespace NTL;
    using namespace std;

    long degree = deg(polynomial);
    int buff = 0;
    int index = 0;
    stringstream ss;
    for (int i = 0; i <= degree; i++) {
	if (index < i / 4) {
	    ss << hex << buff;
	    buff = 0;
	    index = i / 4;
	}
	if (IsOne(coeff(polynomial, i))) {
	    buff |= 1 << (i % 4);
	}
    }
    if (buff != 0) {
	ss << hex << buff;
    }
    ss << flush;
    x = ss.str();
}

/**
 * Convert polynomial to 32-bit unsigned array.
 * Resulted array is used in kernel functions of cuda and OpenCL.
 * @param[out] array array format polynomial
 * @param[in] size max size of array
 * @param polynomial polynomial
 */
static inline void polytoarray(uint32_t array[],
			       int size,
			       NTL::GF2X& polynomial)
{
    using namespace NTL;
    using namespace std;

    long degree = deg(polynomial);
    if (size == 0) {
	return;
    } else if (size < degree / 32 + 1) {
	array[0] = 0;
	return;
    }
    for (int i = 0; i < size; i++) {
	array[i] = 0;
    }
    for (int i = 0; i <= degree; i++) {
	int index = i / 32;
	int pos = i % 32;
	if (IsOne(coeff(polynomial, i))) {
	    array[index] |= 1 << pos;
	}
    }
}

/**
 * converts string to polynomial
 * @param[in] str string
 * @param[out] poly output polynomial
 */
static inline void stringtopoly(NTL::GF2X& poly, std::string& str)
{
    using namespace NTL;
    using namespace std;

    stringstream ss(str);
    char c;
    long p = 0;
    clear(poly);
    while(ss) {
	ss >> c;
	if (!ss) {
	    break;
	}
	if (c >= 'a') {
	    c = c - 'a' + 10;
	} else {
	    c = c - '0';
	}
	for (int j = 0; j < 4; j++) {
	    if (c & (1 << j)) {
		SetCoeff(poly, p, 1);
	    } else {
		SetCoeff(poly, p, 0);
	    }
	    p++;
	}
    }
}

/**
 * calculate the jump polynomial.
 * @param[out] array output string which represents jump polynomial.
 * @param[in] size max size of array
 * @param[in] step jump step of internal state
 * @param[in] characteristic characteristic polynomial
 */
static inline void calc_jump(uint32_t array[],
			     int size,
			     NTL::ZZ& step,
			     NTL::GF2X& characteristic)
{
    using namespace NTL;
    using namespace std;
    GF2X jump;
    PowerXMod(jump, step, characteristic);
    polytoarray(array, size, jump);
}

/**
 * calculate the jump polynomial.
 * @param[out] jump_str output string which represents jump polynomial.
 * @param[in] step jump step of internal state
 * @param[in] characteristic characteristic polynomial
 */
static inline void calc_jump(std::string& jump_str,
			     NTL::ZZ& step,
			     NTL::GF2X& characteristic)
{
    using namespace NTL;
    using namespace std;
    GF2X jump;
    PowerXMod(jump, step, characteristic);
    polytostring(jump_str, jump);
}

#endif

max_stars_repo_name: mkt-matsumoto-lab/MTGP

max_issues_repo_name: mkt-matsumoto-lab/MTGP

max_forks_repo_name: mkt-matsumoto-lab/MTGP


text:
// Copyright (c) 2018 Baidu, Inc. All Rights Reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#include "transaction.h"
#include "transaction_pool.h"
#include "tuple_record.h"
#include <boost/scoped_array.hpp>
#include <gflags/gflags.h>

namespace baikaldb {
DEFINE_bool(disable_wal, true, "disable rocksdb interanal WAL log, only use raft log");
// DEFINE_int32(rocks_transaction_expiration_ms, 600 * 1000,
//         "rocksdb transaction_expiration timeout(us)");

int Transaction::begin() {
    if (nullptr == (_db = RocksWrapper::get_instance())) {
        DB_WARNING("get rocksdb instance failed");
        return -1;
    }
    if (nullptr == (_data_cf = _db->get_data_handle())) {
        DB_WARNING("get rocksdb data column family failed");
        return -1;
    }
    //_txn_opt.lock_timeout = 2000;
    //_txn_opt.expiration = FLAGS_rocks_transaction_expiration_ms;
    if (nullptr == (_txn = _db->begin_transaction(_write_opt, _txn_opt))) {
        DB_WARNING("start_trananction failed");
        return -1;
    }
    last_active_time = butil::gettimeofday_us();
    return 0;
}

int Transaction::begin(rocksdb::Transaction* txn) {
    if (txn == nullptr) {
        DB_WARNING("txn is nullptr");
        return -1;
    }
    if (nullptr == (_db = RocksWrapper::get_instance())) {
        DB_WARNING("get rocksdb instance failed");
        return -1;
    }
    if (nullptr == (_data_cf = _db->get_data_handle())) {
        DB_WARNING("get rocksdb data column family failed");
        return -1;
    }
    _txn = txn;
    last_active_time = butil::gettimeofday_us();
    _is_prepared = true;
    _prepare_time_us = butil::gettimeofday_us();
    _pool->increase_prepared();
    return 0;
}

int Transaction::get_full_primary_key(
        rocksdb::Slice index_bytes,
        rocksdb::Slice pk_bytes,
        IndexInfo& pk_index,
        IndexInfo& index_info,
        char* pk_buf) {
    //remove null flag byte
    //index_bytes.remove_prefix(1);
    int pos = 0;
    for (uint32_t i = 0; i < index_info.pk_pos.size(); ++i) {
        if (index_info.pk_pos[i].first == 1) {
            memcpy(pk_buf + pos, index_bytes.data() + index_info.pk_pos[i].second, pk_index.fields[i].size);
            pos += pk_index.fields[i].size;
        } else if (index_info.pk_pos[i].first == -1) {
            memcpy(pk_buf + pos, pk_bytes.data() + index_info.pk_pos[i].second, pk_index.fields[i].size);
            pos += pk_index.fields[i].size;
        } else {
            DB_WARNING("error.");
            return -1;
        }
    }
    return 0;
}

// start指针为空表示不用判断region start_key
// end指针为空表示不用判断region end_key
// 传入待测试的key和value不包含regionid+table前缀
bool Transaction::fits_region_range(rocksdb::Slice key, rocksdb::Slice value,
        const std::string* start, const std::string* end,
        IndexInfo& pk_index, IndexInfo& index_info) {
    if (!start && !end) {
        return true;
    }
    int ret1 = 1;
    int ret2 = -1;
    if (index_info.type == pb::I_PRIMARY) {
        if (start) {
            ret1 = key.compare(*start);
            if (ret1 < 0) {
                return false;
            }
        }
        if (end && !end->empty()) {
            ret2 = key.compare(*end);
        }
    } else if (index_info.type == pb::I_UNIQ) {
        if (pk_index.length > 0 && index_info.length > 0 && index_info.overlap) {
            boost::scoped_array<char> pk_buf(new(std::nothrow)char[pk_index.length]);
            int ret = get_full_primary_key(key, value, pk_index, index_info, pk_buf.get());
            if (ret != 0) {
                return false;
            }
            rocksdb::Slice primary_key(pk_buf.get(), pk_index.length);
            if (start) {
                ret1 = primary_key.compare(*start);
                if (ret1 < 0) {
                    return false;
                }
            }
            if (end && !end->empty()) {
                ret2 = primary_key.compare(*end);
            }
        } else {
            if (start) {
                ret1 = value.compare(*start);
                if (ret1 < 0) {
                    return false;
                }
            }
            if (end && !end->empty()) {
                ret2 = value.compare(*end);
            }
        }
    } else if (index_info.type == pb::I_KEY) {
        if (pk_index.length > 0 && index_info.length > 0 && index_info.overlap) {
            rocksdb::Slice _value(key);
            _value.remove_prefix(index_info.length);
            boost::scoped_array<char> pk_buf(new(std::nothrow)char[pk_index.length]);
            int ret = get_full_primary_key(key, _value, pk_index, index_info, pk_buf.get());
            if (ret != 0) {
                return false;
            }
            rocksdb::Slice primary_key(pk_buf.get(), pk_index.length);
            if (start) {
                ret1 = primary_key.compare(*start);
                if (ret1 < 0) {
                    return false;
                }
            }
            if (end && !end->empty()) {
                ret2 = primary_key.compare(*end);
            }
        } else {
            if (index_info.length > 0) {
                // index_info为定长
                key.remove_prefix(index_info.length);
            } else {
                // todo, index_info为变长或有Null字段
                TableKey table_key(key);
                uint8_t null_flag = table_key.extract_u8(0);
                int pos = 1;
                for (uint32_t idx = 0; idx < index_info.fields.size(); ++idx) {
                    //flagbit为1且can_null为true同时成立，该字段才真正不存储
                    if (((null_flag >> (7 - idx)) & 0x01)
                            && index_info.fields[idx].can_null) {
                        continue;
                    }
                    if (index_info.fields[idx].type == pb::STRING) {
                        // string 编码后面有\0
                        pos += strlen(key.data() + pos) + 1;
                    } else {
                        pos += index_info.fields[idx].size;
                    }
                }
                key.remove_prefix(pos);
            }
            if (start) {
                ret1 = key.compare(*start);
                if (ret1 < 0) {
                    return false;
                }
            }
            if (end && !end->empty()) {
                ret2 = key.compare(*end);
            }
        }
    }
    if (ret1 >= 0 && ret2 < 0) {
        return true;
    }
    return false;
}

//TODO: finer return status
//return -3 when region not match
int Transaction::put_primary(int64_t region, IndexInfo& pk_index, SmartRecord record) {
    BAIDU_SCOPED_LOCK(_txn_mutex);
    last_active_time = butil::gettimeofday_us();
    MutTableKey key;
    int ret = -1;
    key.append_i64(region).append_i64(pk_index.id);
    //encode key (not allowing prefix, do post field_clear)
    if (0 != key.append_index(pk_index, record.get(), -1, true)) {
        DB_FATAL("Fail to append_index, reg=%ld, tab=%ld", region, pk_index.id);
        return -1;
    }
    std::string value;
    ret = record->encode(value);
    if (ret != 0) {
        DB_WARNING("encode record failed: reg=%ld, tab=%ld", region, pk_index.id);
        return -1;
    }
    auto res = _txn->Put(_data_cf, key.data(), value);
    if (!res.ok()) {
        return -1;
    }
    return 0;
}

//TODO: finer return status
//txt->Put always return OK when using OptimisticTransactionDB
int Transaction::put_secondary(int64_t region, IndexInfo& index, SmartRecord record) {
    BAIDU_SCOPED_LOCK(_txn_mutex);
    last_active_time = butil::gettimeofday_us();
    //IndexInfo index_info = _factory->get_index_info(index);
    if (index.type != pb::I_KEY && index.type != pb::I_UNIQ) {
        DB_WARNING("invalid index type, region_id: %ld, table_id: %ld, index_type:%d", region, index.id, index.type);
        return -1;
    }
    MutTableKey key;
    key.append_i64(region).append_i64(index.id);

    //encode key (not allowing prefix, no clear)
    if(0 != key.append_index(index, record.get(), -1, false)) {
        DB_FATAL("Fail to append_index, reg:%ld, tab:%ld", region, index.id);
        return -1;
    }
    rocksdb::Status res;
    if (index.type == pb::I_KEY) {
        if (0 != record->encode_primary_key(index, key, -1)) {
            DB_FATAL("Fail to append_index, reg:%ld, tab:%ld", region, index.pk);
            return -1;
        }
        res = _txn->Put(_data_cf, key.data(), "");
        //DB_FATAL("data:%s", str_to_hex(key.data()).c_str());
    } else if (index.type == pb::I_UNIQ) {
        MutTableKey pk;
        if (0 != record->encode_primary_key(index, pk, -1)) {
            DB_FATAL("Fail to append_index, reg:%ld, tab:%ld", region, index.pk);
            return -1;
        }
        res = _txn->Put(_data_cf, key.data(), pk.data());
    }
    if (!res.ok()) {
        return -1;
    }
    return 0;
}

//val should be nullptr (create inside the func if the key is found)
//TODO: update return status
int Transaction::get_update_primary(
        int64_t             region,
        IndexInfo&          pk_index,
        SmartRecord         key,
        std::vector<int32_t>& fields,
        GetMode             mode,
        bool                check_region) {
    MutTableKey         _key;
    //full key, no prefix allowed
    if (0 != _key.append_index(pk_index, key.get(), -1, false)) {
        DB_WARNING("Fail to append_index, reg:%ld, tab:%ld", region, pk_index.id);
        return -1;
    }
    return get_update_primary(region, pk_index, TableKey(_key),
            mode, key, fields, false, check_region);
}

//TODO: update return status
int Transaction::get_update_primary(
        int64_t         region,
        IndexInfo&      pk_index,
        const TableKey& key,
        GetMode         mode,
        SmartRecord     val,
        std::vector<int32_t>& fields,
        bool            check_region) {
    return get_update_primary(region, pk_index, key,
            mode, val, fields, true, check_region);
}

//TODO: update return status
//if val == null, don't read value
//return -3 when region not match
int Transaction::get_update_primary(
        int64_t         region,
        IndexInfo&      pk_index,
        const TableKey& key,
        GetMode         mode,
        SmartRecord     val,
        std::vector<int32_t>& fields,
        bool            parse_key,
        bool            check_region) {
    BAIDU_SCOPED_LOCK(_txn_mutex);
    if (_region_info == nullptr) {
        DB_WARNING("no region_info");
        return -1;
    }
    last_active_time = butil::gettimeofday_us();
    int ret = -1;
    //IndexInfo index_info = _factory->get_index_info(table);
    if (pk_index.type != pb::I_PRIMARY) {
        DB_WARNING("invalid index type: %d", pk_index.type);
        return -1;
    }
    // region_info.end_key() == "" 是最后一个region
    if (/*_need_check_region &&*/check_region) {
        rocksdb::Slice pure_key(key.data());
        rocksdb::Slice value;
        if (!fits_region_range(pure_key, value,
            &_region_info->start_key(), &_region_info->end_key(), pk_index, pk_index)) {
            //DB_WARNING("fail to fit: %s", pure_key.ToString(true).c_str());
            return -3;
        }
    }
    MutTableKey _key;
    _key.append_i64(region).append_i64(pk_index.id).append_index(key);

    std::string _value;
    std::string* val_ptr = nullptr;
    if (mode == GET_ONLY || mode == GET_LOCK) {
        val_ptr = &_value;
    }
    rocksdb::Status res;
    if (mode == GET_ONLY) {
        //TimeCost cost;
        res = _txn->Get(_read_opt, _data_cf, _key.data(), val_ptr);
        //DB_NOTICE("txn get time:%ld", cost.get_time());
    } else if (mode == LOCK_ONLY || mode == GET_LOCK) {
        res = _txn->GetForUpdate(_read_opt, _data_cf, _key.data(), val_ptr);
        //DB_WARNING("data: %s %d", _value.c_str(), _value.size());
    } else {
        DB_WARNING("invalid GetMode: %d", mode);
        return -1;
    }

    if (res.ok()) {
        DB_DEBUG("lock ok and key exist");
        if (mode == GET_ONLY || mode == GET_LOCK) {
            //TimeCost cost;
            TupleRecord tuple_record(_value);
            // only decode the required field (field_ids stored in fields)
            if (0 != tuple_record.decode_fields(fields, val)) {
                DB_WARNING("decode value failed: %d", pk_index.id);
                return -1;
            }
            // // 外部传来的record可能包含一些额外的信息需要保留
            // SmartRecord tmp_val = val->clone();
            // ret = tmp_val->decode(_value);
            // if (ret != 0) {
            //     DB_WARNING("decode value failed: %ld", pk_index.id);
            //     return -1;
            // }
            if (parse_key) {
                ret = val->decode_key(pk_index, key);
                if (ret != 0) {
                    DB_WARNING("decode primary index failed: %ld", pk_index.id);
                    return -1;
                }
            }
            //DB_NOTICE("decode time:%ld", cost.get_time());
            //val->merge_from(tmp_val);
        }
    } else if (res.IsNotFound()) {
        DB_DEBUG("lock ok but key not exist");
        return -2;
    } else if (res.IsBusy()) {
        DB_WARNING("lock failed, busy: %s", res.ToString().c_str());
        return -1;
    } else if (res.IsTimedOut()) {
        uint32_t cf_id = _data_cf->GetID();
        std::vector<uint64_t> txn_ids = _txn->GetWaitingTxns(&cf_id, &_key.data());
        if (txn_ids.size() > 0) {
            for (auto txn_id : txn_ids) {
                DB_WARNING("locked by id: %lu", txn_id);
            }
        }
        DB_WARNING("lock failed, timedout: %s", res.ToString().c_str());
        return -1;
    } else {
        DB_WARNING("unknown error: %d, %s", res.code(), res.ToString().c_str());
        return -1;
    }
    return 0;
}

//TODO: update return status
int Transaction::get_update_secondary(
        int64_t             region,
        IndexInfo&          pk_index,
        IndexInfo&          index,
        SmartRecord         key,
        GetMode             mode,
        bool                check_region) {
    last_active_time = butil::gettimeofday_us();

    MutTableKey pk_val;
    int ret = -1;
    if (index.type != pb::I_UNIQ) {
        //DB_WARNING("invalid index type: %d", index.type);
        return -2;
    }
    int res = get_update_secondary(region, pk_index, index, key, mode, pk_val, check_region);
    if (res == 0 && index.type == pb::I_UNIQ && (mode == GET_ONLY || mode == GET_LOCK)) {
        int pos = 0;
        ret = key->decode_primary_key(index, pk_val, pos);
        if (ret != 0) {
            DB_WARNING("decode value failed: %ld", index.pk);
            return -1;
        }
    }
    return res;
}

int Transaction::get_update_secondary(
        int64_t             region,
        IndexInfo&          pk_index,
        IndexInfo&          index,
        const SmartRecord   key,
        GetMode             mode,
        MutTableKey&        pk,
        bool                check_region) {
    BAIDU_SCOPED_LOCK(_txn_mutex);
    if (_region_info == nullptr) {
        DB_WARNING("no region_info");
        return -1;
    }
    last_active_time = butil::gettimeofday_us();
    std::string& pk_val = pk.data();
    if (index.type != pb::I_UNIQ) {
        //DB_WARNING("invalid index type: %d", index.type);
        return -2;
    }
    MutTableKey _key;
    _key.append_i64(region).append_i64(index.id);

    //full key, no prefix allowed
    if (0 != _key.append_index(index, key.get(), -1, false)) {
        DB_FATAL("Fail to append_index, region:%ld,tab:%ld", region, index.id);
        return -1;
    }
    // if (index.type == pb::I_KEY) {
    //     if (0 != _key.append_index(pk_index, key.get(), -1, false)) {
    //         DB_FATAL("Fail to append_pk_index, region:%ld,tab:%ld", region, index.id);
    //         return -1;
    //     }
    // }
    std::string* val_ptr = nullptr;
    if (mode == GET_ONLY || mode == GET_LOCK) {
        val_ptr = &pk_val;
    }
    rocksdb::Status res;
    if (mode == GET_ONLY) {
        res = _txn->Get(_read_opt, _data_cf, _key.data(), &pk_val);
    } else if (mode == LOCK_ONLY || mode == GET_LOCK) {
        res = _txn->GetForUpdate(_read_opt, _data_cf, _key.data(), val_ptr);
    } else {
        DB_WARNING("invalid GetMode: %d", mode);
        return -1;
    }
    if (res.ok()) {
        DB_DEBUG("lock ok and key exist");
    } else if (res.IsNotFound()) {
        DB_DEBUG("lock ok but key not exist");
        return -2;
    } else if (res.IsBusy()) {
        DB_WARNING("lock failed, busy: %s", res.ToString().c_str());
        return -1;
    } else if (res.IsTimedOut()) {
        DB_WARNING("lock failed, timedout: %s", res.ToString().c_str());
        return -1;
    } else {
        DB_WARNING("unknown error: %d, %s", res.code(), res.ToString().c_str());
        return -1;
    }

    if (/*_need_check_region &&*/check_region && (mode == GET_ONLY || mode == GET_LOCK)) {
        rocksdb::Slice pure_key(_key.data());
        pure_key.remove_prefix(2 * sizeof(int64_t));
        rocksdb::Slice value(pk_val);
        if (!fits_region_range(pure_key, value,
            &_region_info->start_key(), &_region_info->end_key(), pk_index, index)) {
            return -3;
        }
    }
    return 0;
}

int Transaction::remove(int64_t region, IndexInfo& index, /*IndexInfo& pk_index,*/
            const SmartRecord key) {
    BAIDU_SCOPED_LOCK(_txn_mutex);
    MutTableKey _key;
    last_active_time = butil::gettimeofday_us();
    _key.append_i64(region).append_i64(index.id);
    if (0 != _key.append_index(index, key.get(), -1, false)) {
        DB_FATAL("Fail to append_index, reg:%ld,tab:%ld", region, index.id);
        return -1;
    }
    if (index.type == pb::I_KEY) {
        if (0 != key->encode_primary_key(index, _key, -1)) {
            DB_FATAL("Fail to append_pk_index, reg:%ld,tab:%ld", region, index.id);
            return -1;
        }
    }
    auto res = _txn->Delete(_data_cf, _key.data());
    if (!res.ok()) {
        DB_WARNING("delete error: code=%d, msg=%s", res.code(), res.ToString().c_str());
        return -1;
    }
    return 0;
}

int Transaction::remove(int64_t region, IndexInfo& index, const TableKey& key) {
    BAIDU_SCOPED_LOCK(_txn_mutex);
    last_active_time = butil::gettimeofday_us();
    MutTableKey _key;
    _key.append_i64(region).append_i64(index.id).append_index(key);
    if (index.type == pb::I_KEY) {
        // cannot append primary index
        DB_WARNING("cannot delete type KEY index");
        return -1;
    }
    auto res = _txn->Delete(_data_cf, _key.data());
    if (!res.ok()) {
        DB_WARNING("delete error: code=%d, msg=%s", res.code(), res.ToString().c_str());
        return -1;
    }
    return 0;
}

rocksdb::Status Transaction::prepare() {
    BAIDU_SCOPED_LOCK(_txn_mutex);
    if (_is_prepared) {
        return rocksdb::Status();
    }
    last_active_time = butil::gettimeofday_us();
    auto res = _txn->Prepare();
    if (res.ok()) {
        if (_pool && !_is_prepared) {
            _pool->increase_prepared();
            DB_WARNING("increase_prepared: %d", _pool->num_prepared());
        }
        _is_prepared = true;
        _prepare_time_us = butil::gettimeofday_us();
    }
    return res;
}

rocksdb::Status Transaction::commit() {
    BAIDU_SCOPED_LOCK(_txn_mutex);
    last_active_time = butil::gettimeofday_us();
    if (_is_finished) {
        DB_WARNING("TransactionWarn: commit a finished txn: %lu", _txn_id);
        return rocksdb::Status();
    }
    if (_txn->GetName().size() != 0 && !_is_prepared) {
        DB_FATAL("TransactionError: commit a un-prepare txn: %lu", _txn_id);
        return rocksdb::Status::Aborted("commit a un-prepare txn");
    }
    auto res = _txn->Commit();
    if (res.ok()) {
        if (_pool && _is_prepared && !_is_finished) {
            _pool->decrease_prepared();
        }
        _is_finished = true;
    }
    return res;
}

rocksdb::Status Transaction::rollback() {
    BAIDU_SCOPED_LOCK(_txn_mutex);
    last_active_time = butil::gettimeofday_us();
    if (_is_finished) {
        DB_WARNING("rollback a finished txn: %lu", _txn_id);
        return rocksdb::Status();
    }
    auto res = _txn->Rollback();
    if (res.ok()) {
        if (_pool && _is_prepared && !_is_finished) {
            _pool->decrease_prepared();
        }
        _is_finished = true;
        _is_rolledback = true;
    }
    return res;
}

int Transaction::set_save_point() {
    BAIDU_SCOPED_LOCK(_txn_mutex);
    last_active_time = butil::gettimeofday_us();
    if (_save_point_seq.empty() || _save_point_seq.top() < _seq_id) {
        _txn->SetSavePoint();
        _save_point_seq.push(_seq_id);
        _save_point_increase_rows.push(num_increase_rows);
    }
    return _seq_id;
}

void Transaction::rollback_to_point(int seq_id) {
    BAIDU_SCOPED_LOCK(_txn_mutex);
    last_active_time = butil::gettimeofday_us();
    if (!_save_point_seq.empty() && _save_point_seq.top() >= seq_id) {
        num_increase_rows = _save_point_increase_rows.top();
        _save_point_seq.pop();
        _save_point_increase_rows.pop();
        _txn->RollbackToSavePoint();
        auto iter = _cache_plan_map.rbegin();
        DB_WARNING("rollback cmd seq_id: %d, num_increase_rows: %ld",
            iter->first, num_increase_rows);
        _cache_plan_map.erase(iter->first);
    }
    if (!_save_point_seq.empty() && _save_point_seq.top() >= seq_id) {
        DB_FATAL("TransactionError: need to rollback to earlier point: top: %d, seq_id: %d",
            _save_point_seq.top(), seq_id);
    }
}

} //nanespace baikaldb