問題
ACコード

#include <bits/stdc++.h>
#define GET_REP(_1, _2, _3, NAME, ...) NAME
#define rep(...) GET_REP(__VA_ARGS__, irep, _rep)(__VA_ARGS__)
#define rep1(...) GET_REP(__VA_ARGS__, irep1, _rep1)(__VA_ARGS__)
#define _rep(i, n) irep (i, 0, n)
#define _rep1(i, n) irep1(i, 1, n)
#define irep(i, a, n) for (int i = a; i < (int)(n); ++i)
#define irep1(i, a, n) for (int i = a; i <= (int)(n); ++i)
#define rrep(i, n) for (int i = (int)(n) - 1; i >= 0; --i)
#define rrep1(i, n) for (int i = (int)(n); i >= 1; --i)
#define allrep(X, x) for (auto &&X : x)
#define all(x) begin(x), end(x)
#ifdef LOCAL
  #include "../../Lib/cout_container.hpp"
  #define debug(x) cout << #x " => " << (x) << endl
#else
  #define debug(x) 0
#endif
using lint = long long;
constexpr int    INF  = 1 << 29;
constexpr lint   INFL = 1LL << 61;
constexpr int    MOD  = (int)1e9 + 7;
constexpr double EPS  = 1e-9;
using namespace std;
namespace { struct INIT { INIT() { cin.tie(0); ios::sync_with_stdio(false); cout << fixed << setprecision(15); } } INIT; }

int main(void) {
  int n, m;
  cin >> n >> m;
  vector<int> a(m), b(m);
  rep (i, m) cin >> a[i];
  const vector<int> match{0, 2, 5, 5, 4, 5, 6, 3, 7, 6};
  rep (i, m) b[i] = match[a[i]];
  unordered_map<int, int> num;
  rep (i, m) num[b[i]] = max(num[b[i]], a[i]);
  sort(all(b));
  b.erase(unique(all(b)), b.end());
  const int bs = b.size();
  vector<map<int, int>> dp(n + 1);
  auto comp =
    [&](map<int, int> &x, map<int, int> &y) {
      if (x[0] != y[0]) {
        return x[0] < y[0];
      } else {
        for (auto xitr = x.rbegin(), yitr = y.rbegin(); xitr != x.rend() && yitr != y.rend(); ++xitr, ++yitr) {
          int xf = xitr->first, xs = xitr->second, yf = yitr->first, ys = yitr->second;
          if (xf != yf) {
            return xf < yf;
          } else {
            if (xs != ys) return xs < ys;
          }
        }
        return true;
      }
    };
  dp[0][0] = 0;
  rep (i, n + 1) {
    if (dp[i].empty()) continue;
    rep (j, bs) {
      if (i + b[j] > n) continue;
      auto tmp = dp[i];
      ++tmp[0];
      ++tmp[num[b[j]]];
      if (comp(dp[i + b[j]], tmp)) dp[i + b[j]] = tmp;
    }
  }
  string ans;
  for (auto itr = dp[n].rbegin(); itr != --dp[n].rend(); ++itr) {
    ans += string(itr->second, itr->first + '0');
  }
  cout << ans << endl;
  return 0;
}

青だけど自力で解けてうれC。
本番だとAからDまで通せていないくらいの時間はかかっているので、そろそろ速度も考えないと。

10分程解析的に解こうと色々考えたが、最終的にDPで良い感じになりそうなことが分かった。
その後色々試行錯誤して以下のようなDPに至った。
DP[使用したマッチの総数]=最大桁数 & map{使用する数字：その数字を使用する回数}
0はマッチで作る数字には存在しないので、0に桁数を格納することでpairがいらなくなる。
使用したマッチの総数iを1～Nまで、使用するj番目の数字を0からMまでループして値を更新する。
j番目の数字a[j]、a[j]を作るのに必要なマッチの数をb[j]とし、newDP[i]を、DP[i][0]を+1、DP[i][a[j]]を+1したものとする。
すると更新はDP[i+b[j]]=max(DP[i+b[j]], newDP[i])となる。

mapの比較は以下のように行う。
まず、桁数が異なれば大きい方を大とする。
次に、等しい場合は先頭要素のkeyを比較し、異なる場合は大きい方を大とする。
先頭要素のkeyも等しい場合は、先頭要素のvalueを比較し、異なる場合は大きい方を大とする。
それも同じ場合は次の要素を比較する。

初期値はDP[0][0]=0とし、DP[i]が空だった場合は更新しない。
答えは、DP[N]から数字が大きい順に指定回数取り出したものになる。

この流れの実現に当たって少し工夫したのが、使用する数字の選別である。
同じ本数で作れる数字が複数あれば、必ず大きい方を使うべきなので小さい方は取り除いても良い。
mを減らせるので速度の改善が期待できる。
しかしこの作業を行ったがために、処理が色々と複雑になってしまった。
1発でACしたものの、実装中は色々とすんなり行かなかった。
実装中にDPに持つ情報を色々変えたので、当初は被りがあるとまずい予定だった。
最終的には被りがあろうと関係ないDPになったので、1～9のMを減らしてもせいぜい定数倍レベルしか改善さないようなこの操作はなくすべきだったかもしれない。

AC後に解法を調べてみると、どの解法も方針は大体同じものの、実装はもっと楽にできたらしい。
特に、DPの値に文字列を持つのはかなり頭良いと思った。
DPで持つ値が文字列と連想配列だとどちらが一般的かと言われれば間違いなく文字列だろうが、なぜか全く思いつかなかった。
DPは一言で言っても、本当に千差万別なのをまた実感した。
もっと慣れていきたい。

問題
ACコード

#include <bits/stdc++.h>
#define GET_REP(_1, _2, _3, NAME, ...) NAME
#define rep(...) GET_REP(__VA_ARGS__, irep, _rep)(__VA_ARGS__)
#define rep1(...) GET_REP(__VA_ARGS__, irep1, _rep1)(__VA_ARGS__)
#define _rep(i, n) irep (i, 0, n)
#define _rep1(i, n) irep1(i, 1, n)
#define irep(i, a, n) for (int i = a; i < (int)(n); ++i)
#define irep1(i, a, n) for (int i = a; i <= (int)(n); ++i)
#define rrep(i, n) for (int i = (int)(n) - 1; i >= 0; --i)
#define rrep1(i, n) for (int i = (int)(n); i >= 1; --i)
#define allrep(X, x) for (auto &&X : x)
#define all(x) (x).begin(), (x).end()
#ifdef LOCAL
  #include "../../Lib/cout_container.hpp"
  #define debug(x) cerr << #x " => " << (x) << endl
#else
  #define debug(x) 0
#endif
using lint = long long;
constexpr int    INF  = 1 << 29;
constexpr lint   INFL = 1LL << 61;
constexpr int    MOD  = (int)1e9 + 7;
constexpr double EPS  = 1e-9;
using namespace std;
namespace { struct INIT { INIT() { cin.tie(0); ios::sync_with_stdio(false); cout << fixed << setprecision(15); } } INIT; }

template <long long mod_num>
class ModInt {
  long long x;

  public:
  constexpr ModInt(const long long num = 0) : x(num % mod_num) { if (x < 0) x += mod_num; }
  constexpr bool   operator==(const ModInt &r) const { return this->x == r.x; }
  constexpr bool   operator!=(const ModInt &r) const { return this->x != r.x; }
  constexpr bool   operator<(const ModInt &r) const { return this->x < r.x; }
  constexpr bool   operator>(const ModInt &r) const { return this->x > r.x; }
  constexpr bool   operator<=(const ModInt &r) const { return this->x <= r.x; }
  constexpr bool   operator>=(const ModInt &r) const { return this->x >= r.x; }
  constexpr ModInt operator+(const ModInt &r) const { return ModInt(*this) += r; }
  constexpr ModInt operator-(const ModInt &r) const { return ModInt(*this) -= r; }
  constexpr ModInt operator*(const ModInt &r) const { return ModInt(*this) *= r; }
  constexpr ModInt operator/(const ModInt &r) const { return ModInt(*this) /= r; }
  friend ModInt operator+(const long long &l, const ModInt &r) { return ModInt(l) += r; }
  friend ModInt operator-(const long long &l, const ModInt &r) { return ModInt(l) -= r; }
  friend ModInt operator*(const long long &l, const ModInt &r) { return ModInt(l) *= r; }
  friend ModInt operator/(const long long &l, const ModInt &r) { return ModInt(l) /= r; }
  constexpr ModInt &operator++() { return *this += 1; }
  constexpr ModInt &operator--() { return *this -= 1; }
  constexpr ModInt &operator+=(const ModInt &r) {
    x += r.x;
    if (x >= mod_num) x -= mod_num;
    return *this;
  }
  constexpr ModInt &operator-=(const ModInt &r) {
    if (x < r.x) x += mod_num;
    x -= r.x;
    return *this;
  }
  constexpr ModInt &operator*=(const ModInt &r) {
    x = x * r.x % mod_num;
    return *this;
  }
  constexpr ModInt &operator/=(const ModInt &r) {
    *this *= r.inverse();
    return *this;
  }
  constexpr ModInt inverse(void) const {
    long long a = x, b = mod_num, u = 1, v = 0;
    while (b) {
      long long t = a / b;
      swap(a -= t * b, b);
      swap(u -= t * v, v);
    }
    return ModInt(u);
  }
  constexpr ModInt pow(long long n) const { return pow(x, n); }
  static ModInt pow(long long x, long long n) {
    ModInt ret(1), mul(x);
    while (n) {
      if (n & 1) ret *= mul;
      mul *= mul;
      n  >>= 1;
    }
    return ret;
  }
  friend istream &operator>>(istream &is, ModInt &r) {
    long long n;
    is >> n;
    r = ModInt(n);
    return is;
  }
  friend ostream &operator<<(ostream &os, const ModInt &r) { return os << r.x; }
  explicit operator int() const { return x; }
  explicit operator long long() const { return x; }
  explicit operator double() const { return x; }
};

template <typename T>
class SegmentTree {
  template <typename B>
  struct BinaryOperation {
    const static B unit = 0;
    B operator()(const B &a, const B &b) const { return a + b; }
  };
  BinaryOperation<T> bo;
  const T   unit;
  vector<T> tree;
  int v_size;

  public:
  SegmentTree(const int n = 0, const T unit = BinaryOperation<T>::unit) : unit(unit) { init(n); }
  SegmentTree(const vector<T> &v, const T unit = BinaryOperation<T>::unit) : unit(unit) { init(v); }
  void init(int n) {
    v_size = 1;
    while (v_size < n) v_size <<= 1;
    tree.assign(2 * v_size, unit);
  }
  void init(const vector<T> &v) {
    init(v.size());
    copy(v.begin(), v.end(), tree.end() - v_size);
    for (int i = v_size - 1; i > 0; --i) tree[i] = bo(tree[i * 2], tree[i * 2 + 1]);
  }
  void update(int i, const int x) {
    i += v_size;
    tree[i] = x;
    while (i >>= 1) tree[i] = bo(tree[i * 2], tree[i * 2 + 1]);
  }
  T query(int l, int r) const {
    T ret = unit;
    for (l += v_size, r += v_size; l < r; l = (l + 1) >> 1, r >>= 1) {
      if (l & 1) ret = bo(ret, tree[l]);
      if (r & 1) ret = bo(ret, tree[r - 1]);
    }
    return ret;
  }
  int size(void) const { return v_size; }
  T   operator[](const int i) const { return tree[i + v_size]; }
};

using mint = ModInt<998244353>;
int main(void) {
  int n;
  cin >> n;
  vector<pair<int, int>> xy(n);
  rep (i, n) cin >> xy[i].first >> xy[i].second;
  sort(all(xy));
  rep (i, n) xy[i].first = i;
  auto comp = [](pair<int, int> &a, pair<int, int> &b) { return a.second < b.second; };
  sort(all(xy), comp);
  rep (i, n) xy[i].second = i;
  vector<int> v(n, 1);
  SegmentTree<int> b(n), u(v);
  mint sum = mint::pow(2, n - 1) * n;
  rep (i, n) {
    u.update(xy[i].first, 0);
    const int  bl = b.query(0, xy[i].first), br = b.query(xy[i].first + 1, n);
    const int  ul = u.query(0, xy[i].first), ur = u.query(xy[i].first + 1, n);
    const mint bln = mint::pow(2, bl) - 1, brn = mint::pow(2, br) - 1;
    const mint uln = mint::pow(2, ul) - 1, urn = mint::pow(2, ur) - 1;
    sum += bln * urn * mint::pow(2, br + ul);
    sum += brn * uln * mint::pow(2, bl + ur);
    sum -= bln * brn * uln * urn;
    b.update(xy[i].first, 1);
  }
  cout << sum << endl;
  return 0;
}

現状での新ABC最終問題にして、初実装をいくつか経験させてもらえてとてもためになった問題。
考察部分では、各点ごとに何回含まれるかを考えれば良さそうというところまでは思いつけた。
頭が固くてそこから先には進めなかった。
解説を見たらあまりにそりゃそうだ感があって拍子抜けしてしまった。
座標圧縮の必要性はいまいち感じなかったのだが、セグメントツリーに乗せてなるほどとなった。

今回大きく経験になったポイントが3つ！

・包除原理
内容は知っていたが問題としては初。
ベン図的な表現なら分かりやすいのだが、何通りを求める問題だと何が論理積なのか一瞬迷った。
何通りかを考えるなら論理積は掛け算になるはずなので、今見ている点の左上に位置する点から1つ以上選ぶ選び方をA、右上をB、左下をC、右下をDとすると、A*C+B*D-A*B*C*Dとなるはずである。
しかしこれだとA*B*C*Dの方が大きいのでは？？？となったが、AとCの論理積ではAとC以外が二値のどちらで合っても良いので、A*CではなくA*C*(右上と左下の点の全ての選び方)が正しかった。

・座標圧縮
蟻本で眺めただけだが、内容のわりに実装難しくないか？という感想を持ったような覚えがある。
とりあえずその記憶は置いておいて、自分なりに実装してみるととても簡単だった。
被りがあったら少しややこしくなるのかなという気はするが本当に少しな気がする。

・セグメントツリーの活用
相当有名かつ汎用性の高そうなデータ構造なので、ライブラリとして実装だけはしていた。
確認用として、AOJのRMQの問題は解いているのだが、実際に問題で使うことになったのは初。
とは言っても過去にABC038 D プレゼントでBITを使っており、今回の問題もBITが使える以上セグメントツリーを使える機会としては初ではないのだが、今回はBITを使わずにセグメントツリーを使ったので事実上初である。
使い方も前回同様、数列の値のみでなくインデックスの値にも意味を持たせるもの。
やっぱりこの使い方は全然思いつけない。
調べたところ反転数を計算するときにも似たような使い方をするらしいので、反転数の問題が出てきた時はすぐに思いつきたい。


ここまで経験値になりそうな問題は初かもしれない。
なりそうだからなったかどうかは分からない。
なってると良いなぁ。

問題
ACコード

#include <bits/stdc++.h>
#define GET_REP(_1, _2, _3, NAME, ...) NAME
#define rep(...) GET_REP(__VA_ARGS__, irep, _rep)(__VA_ARGS__)
#define rep1(...) GET_REP(__VA_ARGS__, irep1, _rep1)(__VA_ARGS__)
#define _rep(i, n) irep (i, 0, n)
#define _rep1(i, n) irep1(i, 1, n)
#define irep(i, a, n) for (int i = a; i < (int)(n); ++i)
#define irep1(i, a, n) for (int i = a; i <= (int)(n); ++i)
#define rrep(i, n) for (int i = (int)(n) - 1; i >= 0; --i)
#define rrep1(i, n) for (int i = (int)(n); i >= 1; --i)
#define allrep(X, x) for (auto &&X : x)
#define all(x) (x).begin(), (x).end()
#ifdef LOCAL
  #include "../../Lib/cout_container.hpp"
  #define debug(x) cerr << #x " => " << (x) << endl
#else
  #define debug(x) 0
#endif
using lint = long long;
constexpr int    INF  = 1 << 29;
constexpr lint   INFL = 1LL << 61;
constexpr int    MOD  = (int)1e9 + 7;
constexpr double EPS  = 1e-9;
using namespace std;
namespace { struct INIT { INIT() { cin.tie(0); ios::sync_with_stdio(false); cout << fixed << setprecision(15); } } INIT; }

class VariableModInt {
  long long x;
  static long long &mod() {
    static long long mod_num = 1;
    return mod_num;
  }

  public:
  VariableModInt(const long long num = 0) : x(num % get_mod()) { if (x < 0) x += get_mod(); }
  bool operator==(const VariableModInt &r) const { return this->x == r.x; }
  bool operator!=(const VariableModInt &r) const { return this->x != r.x; }
  bool operator<(const VariableModInt &r) const { return this->x < r.x; }
  bool operator>(const VariableModInt &r) const { return this->x > r.x; }
  bool operator<=(const VariableModInt &r) const { return this->x <= r.x; }
  bool operator>=(const VariableModInt &r) const { return this->x >= r.x; }
  VariableModInt operator+(const VariableModInt &r) const { return VariableModInt(*this) += r; }
  VariableModInt operator-(const VariableModInt &r) const { return VariableModInt(*this) -= r; }
  VariableModInt operator*(const VariableModInt &r) const { return VariableModInt(*this) *= r; }
  VariableModInt operator/(const VariableModInt &r) const { return VariableModInt(*this) /= r; }
  friend VariableModInt operator+(const long long &l, const VariableModInt &r) { return VariableModInt(l) += r; }
  friend VariableModInt operator-(const long long &l, const VariableModInt &r) { return VariableModInt(l) -= r; }
  friend VariableModInt operator*(const long long &l, const VariableModInt &r) { return VariableModInt(l) *= r; }
  friend VariableModInt operator/(const long long &l, const VariableModInt &r) { return VariableModInt(l) /= r; }
  VariableModInt &operator++() { return *this += 1; }
  VariableModInt &operator--() { return *this -= 1; }
  VariableModInt &operator+=(const VariableModInt &r) {
    x += r.x;
    if (x >= get_mod()) x -= get_mod();
    return *this;
  }
  VariableModInt &operator-=(const VariableModInt &r) {
    if (x < r.x) x += get_mod();
    x -= r.x;
    return *this;
  }
  VariableModInt &operator*=(const VariableModInt &r) {
    x = x * r.x % get_mod();
    return *this;
  }
  VariableModInt &operator/=(VariableModInt r) {
    *this *= r.inverse();
    return *this;
  }
  VariableModInt inverse(void) const {
    long long a = x, b = get_mod(), u = 1, v = 0;
    while (b) {
      long long t = a / b;
      swap(a -= t * b, b);
      swap(u -= t * v, v);
    }
    return VariableModInt(u);
  }
  VariableModInt pow(long long n) const {
    VariableModInt ret(1), mul(x);
    while (n) {
      if (n & 1) ret *= mul;
      mul *= mul;
      n  >>= 1;
    }
    return ret;
  }
  friend istream &operator>>(istream &is, VariableModInt &r) {
    long long n;
    is >> n;
    r = VariableModInt(n);
    return is;
  }
  friend ostream &operator<<(ostream &os, const VariableModInt &r) { return os << r.x; }
  static void set_mod(const long long mod_num) { mod() = mod_num; }
  static long long get_mod() { return mod(); }
  explicit operator int() const { return x; }
  explicit operator long long() const { return x; }
  explicit operator double() const { return x; }
};

template <typename T>
vector<T> lagrange_interpolation(const vector<T> &yvec) {
  const int n = yvec.size() - 1;
  vector<T> ret(n + 1), fac(n + 1), c(n + 2);
  fac[0] = 1;
  for (int i = 1; i <= n; ++i) fac[i] = fac[i - 1] * i;
  c[0] = 0, c[1] = 1;
  for (int i = 1; i <= n; ++i) {
    for (int j = i + 1; j > 0; --j) {
      c[j] = c[j - 1] - c[j] * i;
    }
  }
  for (int i = 0; i <= n; ++i) {
    T c2 = 0, qi = yvec[i] / (fac[n - i] * fac[i]);
    if ((n - i) & 1) qi *= -1;
    for (int j = n; j >= 0; --j) {
      c2 = c[j + 1] + c2 * i;
      ret[j] += c2 * qi;
    }
  }
  return ret;
}

using mint = VariableModInt;
int main(void) {
  int p;
  cin >> p;
  mint::set_mod(p);
  vector<mint> a(p);
  rep (i, p) cin >> a[i];
  auto ans = lagrange_interpolation(a);
  rep (i, ans.size()) cout << ans[i] << (i != ans.size() - 1 ? " " : "\n");
  return 0;
}

解説が少ない！！！
PDF解法のフェルマーの小定理を使う解法はかなり頭の良い解法らしいし、そもそもフェルマーの小定理なんて知らないので、ラグランジュ補間で頑張ろうと思ってライブラリの作成から始めた。
ラグランジュ補間自体は理解できるのだが、係数を復元するというのでかなり悩まされた。
おそらく丁寧に書かれているであろう解説記事はあったものの、解説を読んでもコードを読んでもとにかく理解できない。
数学力の無さが原因なのかと少し寂しい気持ちになる。
総じて大卒程度の数学力はあると思うが、実際は大卒ではないし、受けた数学の講義も数学Iや応用数学等で、幾何学や代数学等とジャンル分けして詳しく学べる講義もなかった。
数学とは違うところでFFTなんかは出てきたのだが、数学的な演習はしていないのでどういうところで使えるか等は分からない。
競プロやっていると少し大学行っておけば良かったなという気持ちになる。

アルゴリズムは理解してACできたものの、解説記事に関しては未だに理解していない。
分からない部分は、結局のところ(x-xi)の総乗の係数の計算をやっているだけだと思うので、解説記事の方法は置いておいて自分流になんとか完成させた。
オーダーもO(N^2)に収まっているので問題ないだろう。
こういう他の解説は完全に理解できなかったが、自分なりに実装できた問題のためにブログを書いていると言っても過言ではない。
後から理解できなかったところが理解できて、過去の自分の解法と比較できるようになりたい。

コードで計算した順番に内容を書いておく。
求める多項式をN次多項式として、入力の点をN+1点とする。
・0!~N!まで階乗を計算する。
DP的に計算するだけでO(N)。
階乗の計算は高速化のために事前計算として行っているのだが、オーダー的には改善されない。
係数の復元ではなく、任意のaに対してf(a)を求める場合にオーダーレベルで高速化される。

・N+1次式{(x-xi)の総乗}の係数(C0~CN+1)を計算する。
解説記事はおそらくここを解説しているであろう部分が理解できなかった。
自分なりに愚直に計算した。
各桁ごとに計算し、(x-xi)をかけるごとに係数の数が1つ増えていくので、全体としてはO(N^2)。

・ラグランジュ補間多項式の各iについて直接計算できる部分Qiを計算する。
xの係数に関係なく、xiに代入することで直接計算できる部分をQiとおく。
具体的には、ラグランジュ補間多項式
\[L(x) = \sum_{i=0}^{N+1} y_i \frac{\prod_{j=0, j \neq i}^{N+1} (x - x_j)}{\prod_{j=0, j \neq i}^{N+1} (x_i - x_j)}\]のxがかからない部分。
\[Q(x) = \frac{y_i} {\prod_{j=0, j \neq i}^{N+1} (x_i - x_j)}\]
xiが0~Nまで連続しているので、(N-i)!*i!*-1^(N-i)で計算できる。
事前に階乗を計算しているのでO(1)。

・i=jを含まないN次式{(xi-xj)の総乗}の係数(C20~C2N)を各iについて計算する。
C2=C/(xi-xi)なので、c2はcを使って上位の桁から各桁ごとにO(1)で求められ、全桁でO(N)。
各iについて計算するのでO(N^2)。

・C2*Qiを各iについて全て足し合わせる。
ラグランジュ補間の公式に従って全て足し合わせる。
今回は係数について考えているので、係数ごとにQiをかけて足し合わせる。
QiがO(N)、C2がO(N^2)なので、全体としてO(N^2)

ラグランジュ補間の公式の見た目があまりにも複雑で、通常の書き方ではとても理解できそうになかったので、急遽数式の表示スクリプトを導入した。
しかしLATEXの書き方微塵も覚えていなくて笑った。
5,6年前？に学校で少し習ったのだが、それから自分含めクラスメイトが本当に誰一人として使っていないんじゃないかというくらい使わなかったLATEX形式。
論文は知る限りみんなWord。
当時は難しく感じたが、かなり簡単に書けてこんなところでも成長を感じてしまう。
最も当時は段落や章等も含んだ1つの文書として書けるようになる講義だったので、さすがに今回ほど簡単とはいかないだろうが、プログラムを触ったこともなかった時代と比べたらさすがに成長しているだろう。

ついでに今回、シンタックスハイライトで実体参照していないがためにtemplateあたりが崩れているのに気付いた。
変換しなくても正常に見えるじゃーんと思っていたが実はダメだった・・・。
わりと問題ない部分ではあるが、一応今回から実体参照にして修正することにした。
どんどんブログ書くのが面倒臭くなっていく。

問題
ACコード

#include <bits/stdc++.h>
#define GET_REP(_1, _2, _3, NAME, ...) NAME
#define rep(...) GET_REP(__VA_ARGS__, irep, _rep)(__VA_ARGS__)
#define rep1(...) GET_REP(__VA_ARGS__, irep1, _rep1)(__VA_ARGS__)
#define _rep(i, n) irep (i, 0, n)
#define _rep1(i, n) irep1(i, 1, n)
#define irep(i, a, n) for (int i = a; i < (int)(n); ++i)
#define irep1(i, a, n) for (int i = a; i <= (int)(n); ++i)
#define rrep(i, n) for (int i = (int)(n) - 1; i >= 0; --i)
#define rrep1(i, n) for (int i = (int)(n); i >= 1; --i)
#define allrep(X, x) for (auto &&X : x)
#define all(x) (x).begin(), (x).end()
#ifdef LOCAL
  #include "../../Lib/cout_container.hpp"
  #define debug(x) cerr << #x " => " << (x) << endl
#else
  #define debug(x) 0
#endif
using lint = long long;
constexpr int    INF  = 1 << 29;
constexpr lint   INFL = 1LL << 61;
constexpr int    MOD  = (int)1e9 + 7;
constexpr double EPS  = 1e-9;
using namespace std;
namespace { struct INIT { INIT() { cin.tie(0); ios::sync_with_stdio(false); cout << fixed << setprecision(15); } } INIT; }

template 
class ModInt {
  long long n_;

  public:
  constexpr ModInt(const long long num = 0) : n_(num % mod_num) {}
  constexpr const long long&value() const { return n_; }
  constexpr bool   operator==(const ModInt &r) const { return this->n_ == r.n_; }
  constexpr bool   operator==(const long long &r) const { return this->n_ == r; }
  constexpr bool   operator!=(const ModInt &r) const { return this->n_ != r.n_; }
  constexpr bool   operator!=(const long long &r) const { return this->n_ != r; }
  constexpr bool   operator<(const ModInt &r) const { return this->n_ < r.n_; }
  constexpr bool   operator<(const long long &r) const { return this->n_ < r; }
  constexpr bool   operator>(const ModInt &r) const { return this->n_ > r.n_; }
  constexpr bool   operator>(const long long &r) const { return this->n_ > r; }
  constexpr bool   operator<=(const ModInt &r) const { return this->n_ <= r.n_; }
  constexpr bool   operator<=(const long long &r) const { return this->n_ <= r; }
  constexpr bool   operator>=(const ModInt &r) const { return this->n_ >= r.n_; }
  constexpr bool   operator>=(const long long &r) const { return this->n_ >= r; }
  constexpr ModInt operator+(const ModInt &r) const { return ModInt(*this) += r; }
  constexpr ModInt operator-(const ModInt &r) const { return ModInt(*this) -= r; }
  constexpr ModInt operator*(const ModInt &r) const { return ModInt(*this) *= r; }
  constexpr ModInt operator/(const ModInt &r) const { return ModInt(*this) /= r; }
  constexpr ModInt operator+(const long long &r) const { return ModInt(*this) += ModInt(r); }
  constexpr ModInt operator-(const long long &r) const { return ModInt(*this) -= ModInt(r); }
  constexpr ModInt operator*(const long long &r) const { return ModInt(*this) *= ModInt(r); }
  constexpr ModInt operator/(const long long &r) const { return ModInt(*this) /= ModInt(r); }
  friend ModInt operator+(const long long &l, const ModInt &r) { return ModInt(l) += r; }
  friend ModInt operator-(const long long &l, const ModInt &r) { return ModInt(l) -= r; }
  friend ModInt operator*(const long long &l, const ModInt &r) { return ModInt(l) *= r; }
  friend ModInt operator/(const long long &l, const ModInt &r) { return ModInt(l) /= r; }
  constexpr ModInt &operator++() { return *this += 1; }
  constexpr ModInt &operator--() { return *this -= 1; }
  constexpr ModInt &operator+=(const ModInt &r) {
    n_ += r.n_;
    if (n_ >= mod_num) n_ -= mod_num;
    return *this;
  }
  constexpr ModInt &operator-=(const ModInt &r) {
    if (n_ < r.n_) n_ += mod_num;
    n_ -= r.n_;
    return *this;
  }
  constexpr ModInt &operator*=(const ModInt &r) {
    n_ = n_ * r.n_ % mod_num;
    return *this;
  }
  constexpr ModInt &operator/=(ModInt r) {
    long long exp = mod_num - 2;
    while (exp) {
      if (exp & 2) *this *= r;
      r   *= r;
      exp /= 2;
    }
    return *this;
  }
  friend istream &operator>>(istream &is, ModInt<mod_num> &r) {
    is >> r.n_;
    r.n_ %= r.mod_num;
    return is;
  }
  friend ostream &operator<<(ostream &os, const ModInt<mod_num> &r) { return os << r.n_; }
  explicit operator int() const { return n_; }
  explicit operator long long() const { return n_; }
  explicit operator double() const { return n_; }
};

bool check(const lint n, const int dig) { return n & (1LL << (59 - dig)); }

using mint = ModInt;
int main(void) {
  lint l, r;
  cin >> l >> r;
  vector<vector<vector<vector>>> dp(61, vector<vector<vector>>(2, vector<vector>(2, vector(2))));
  // 60桁目(0-index)以上が1の数は必ずRより大きくなってしまうので0のみから成る1通りのみ
  dp[0][0][0][0] = 1;
  rep (d, 60) {
    bool cl = check(l, d), cr = check(r, d);
    rep (i, 2) {
      rep (j, 2) {
        rep (k, 2) {
          auto val = dp[d][i][j][k];
          rep (x, 2) {
            rep (y, 2) {
              // 条件に当てはまるかチェック
              // xの方が大きい
              if (x && !y) continue;
              // 最上位ビットが異なる
              if (!k && (x ^ y)) continue;
              // L以上が未確定な状態でLよりも小さい
              if (!i && cl && (!x || !y)) continue;
              // R以下が未確定な状態でRよりも大きい
              if (!j && !cr && (x || y)) continue;
              //
              // 遷移先が変わるかチェック
              int ni = i, nj = j, nk = k;
              // L以上が確定
              if (!cl && x) ni = 1;
              // R以下が確定
              if (cr && !y) nj = 1;
              // 最上位ビットが確定
              if (!k && x) nk = 1;
              //
              // 遷移
              dp[d + 1][ni][nj][nk] += val;
            }
          }
        }
      }
    }
  }
  mint ans = 0;
  rep (i, 2) {
    rep (j, 2) {
      rep (k, 2) {
        ans += dp[60][i][j][k];
      }
    }
  }
  cout << ans << endl;
  return 0;
}

過去最高に混乱したDP。
分からなさすぎて珍しくコードにコメントを書いた。
分からなかったのは遷移条件ではなく遷移方法(DPの値全てに(x,y)の値を持って遷移すること)なのだが、1から記録していかないとつらいと感じてコメントを残すことにした。
解法も当然のごとく思いつなかったのだが、解法を見た後でも自力で実装できなかった。
添字4つの時点で初体験なのに、それに加えて別途遷移先を判断する数が必要とは。
各所で「2つの数字を管理する」というようなことが書かれていたが、なるほどこういうことかと実装できて初めて意味を理解する。
確かにこのDPで遷移先を全て網羅できている。
網羅できていることは分かるし、遷移条件も分かる。
だが、全体として見ると自然に捉えられないなぁという感想。
経験が足りない分余計な思考が入ってしまうので、それら全てに間違っているという判断を下しながら、全体を結合して考えても合っていると判断するには脳のメモリが足りない。
最近ますますDP is 経験を認識してきた。

最初考えたのが、(x,y)の値を持つのではなく、以下のように何回足し合わせられるか、というような遷移である。
dp[d][0][0][0] = dp[d+1][0][0][0] * A + dp[d+1][0][0][1] * B + ...
dp[d][1][0][0] = dp[d+1][0][0][0] * X + dp[d+1][0][0][1] * Y + ...
確かにこれでもできそうだが、(x,y)それぞれについて考えてこの遷移量なので、この方針だとさらにややこしいことになりそう・・・。
そもそもこの方針を思いついたのは、「DPは余計にループせずに全遷移を網羅するもの」という考えで組もうとしているから。
確かにループの数は小さくなるのだが、結局四則演算の数はあまり変わらなさそうなので全体としてもあまり変わらないかもしれない。
(x, y)の状態が大きかったら変わる？かもしれない？
だとしたら難易度のレベルがいくつか上がりそうなので考えるべきではないと判断することにした。

大事なのは、全ての状態・遷移を網羅できて、なおかつ実行可能時間内に収めること。
余計なこと考えるのはせめて黄色くらいになってからにしようと思い直した。
余計なことではないだろうが、慣れていない種類の難しい問題の詳細はそりゃ相当難しい。