CPivot.cpp

/**
 * code generated by JHelper
 * More info: https://github.com/AlexeyDmitriev/JHelper
 * @author RiaD
 */

#include <iostream>
#include <fstream>

#include <iostream>
#include <cassert>
#include <vector>
#include <cmath>


#include <iterator>


#include <string>
#include <stdexcept>

#ifndef SPCPPL_ASSERT
  #ifdef SPCPPL_DEBUG
    #define SPCPPL_ASSERT(condition) \
    if(!(condition)) { \
      throw std::runtime_error(std::string() + #condition + " in line " + std::to_string(__LINE__) + " in " + __PRETTY_FUNCTION__); \
    }
  #else
    #define SPCPPL_ASSERT(condition)
  #endif
#endif


/**
* Support decrementing and multi-passing, but not declared bidirectional(or even forward) because
* it's reference type is not a reference.
*
* It doesn't return reference because
* 1. Anyway it'll not satisfy requirement [forward.iterators]/6
*   If a and b are both dereferenceable, then a == b if and only if *a and
*   b are bound to the same object.
* 2. It'll not work with reverse_iterator that returns operator * of temporary which is temporary for this iterator
*
* Note, reverse_iterator is not guaranteed to work  now too since it works only with bidirectional iterators,
* but it's seems to work at least on my implementation.
*
* It's not really useful anywhere except iterating anyway.
*/
template <typename T>
class IntegerIterator {
public:
  using value_type = T;
  using difference_type = std::ptrdiff_t;
  using pointer = T*;
  using reference = T;
  using iterator_category = std::input_iterator_tag;

  explicit IntegerIterator(T value): value(value) {

  }

  IntegerIterator& operator++() {
    ++value;
    return *this;
  }

  IntegerIterator operator++(int) {
    IntegerIterator copy = *this;
    ++value;
    return copy;
  }

  IntegerIterator& operator--() {
    --value;
    return *this;
  }

  IntegerIterator operator--(int) {
    IntegerIterator copy = *this;
    --value;
    return copy;
  }

  T operator*() const {
    return value;
  }

  bool operator==(IntegerIterator rhs) const {
    return value == rhs.value;
  }

  bool operator!=(IntegerIterator rhs) const {
    return !(*this == rhs);
  }

private:
  T value;
};

template <typename T>
class IntegerRange {
public:
  IntegerRange(T begin, T end): begin_(begin), end_(end) {
    SPCPPL_ASSERT(begin <= end);
  }

  IntegerIterator<T> begin() const {
    return IntegerIterator<T>(begin_);
  }

  IntegerIterator<T> end() const {
    return IntegerIterator<T>(end_);
  }

private:
  T begin_;
  T end_;
};

template <typename T>
class ReversedIntegerRange {
  using IteratorType = std::reverse_iterator<IntegerIterator<T>>;
public:
  ReversedIntegerRange(T begin, T end): begin_(begin), end_(end) {
    SPCPPL_ASSERT(begin >= end);
  }

  IteratorType begin() const {
    return IteratorType(IntegerIterator<T>(begin_));
  }

  IteratorType end() const {
    return IteratorType(IntegerIterator<T>(end_));
  }

private:
  T begin_;
  T end_;
};

template <typename T>
IntegerRange<T> range(T to) {
  return IntegerRange<T>(0, to);
}

template <typename T>
IntegerRange<T> range(T from, T to) {
  return IntegerRange<T>(from, to);
}

template <typename T>
IntegerRange<T> inclusiveRange(T to) {
  return IntegerRange<T>(0, to + 1);
}

template <typename T>
IntegerRange<T> inclusiveRange(T from, T to) {
  return IntegerRange<T>(from, to + 1);
}

template <typename T>
ReversedIntegerRange<T> downrange(T from) {
  return ReversedIntegerRange<T>(from, 0);
}

template <typename T>
ReversedIntegerRange<T> downrange(T from, T to) {
  return ReversedIntegerRange<T>(from, to);
}

template <typename T>
ReversedIntegerRange<T> inclusiveDownrange(T from) {
  return ReversedIntegerRange<T>(from + 1, 0);
}

template <typename T>
ReversedIntegerRange<T> inclusiveDownrange(T from, T to) {
  return ReversedIntegerRange<T>(from + 1, to);
}


using namespace std;

using Row = pair<vector<int64_t>, double>;

std::size_t reduceToRowEchelonForm(vector<Row>& matrix, size_t cols) {
  std::size_t rank = 0;

  for (std::size_t column: range(cols)) {
    bool found = false;
      /*
    for (std::size_t row: range(rank, matrix.size())) {
      if (matrix[row].first[column] != 0) {
        using std::swap;
        found = true;
        swap(matrix[row], matrix[rank]);
        break;
      }
    }
    */
    int64_t smallest = 100000000000;
    int64_t bestrow = -1;
    for (std::size_t row: range(rank, matrix.size())) {
      if (matrix[row].first[column] != 0) {
         if (found == false)
             bestrow = row;
         found = true;
         if (matrix[rank].first[column] != 0) {
             int64_t multiplier = matrix[row].first[column] / matrix[rank].first[column];
             if (abs(multiplier) < smallest) {
                 smallest = abs(multiplier);
                 bestrow = row; 
             }
         }
      }
    }
    if (!found) {
      continue;
    }
    {
        using std::swap;
        swap(matrix[bestrow], matrix[rank]);
    }

    for (std::size_t r: range(rank + 1, matrix.size())) {
      while (matrix[r].first[column] != 0) {
        int64_t multiplier = matrix[r].first[column] / matrix[rank].first[column];
        for (std::size_t c: range(column, cols)) {
          matrix[r].first[c] -= multiplier * matrix[rank].first[c];
        }
        matrix[r].second -= multiplier * matrix[rank].second;
        if (matrix[r].first[column] != 0) {
          swap(matrix[r], matrix[rank]);
        }
      }

    }
    ++rank;
  }
  return rank;
}


class TaskC {
public:
  static constexpr int kStressCount = 0;

  static void generateTest(std::ostream& test) {
  }

  void solve(std::istream& in, std::ostream& out) {
    //static int testnumber = 0;
    //out << "Case #" << ++testnumber << ": ";

    int y, c;
    in >> y >> c;
    //y = 10;
    //c = 100;
    --y;
    int q;
    in >> q;
    //q = 1000;
    int row_len = 2 * c + y;


    vector<Row> known_rows;

    for (int i: range(y)) {
      double x;
      in >> x;
      //x = 2;
      if (x > -0.5) {
        vector<int64_t> row(row_len);
        row[2 * c + i] = 1;
        known_rows.emplace_back(move(row), log(x));
      }
    }

    for (int cur_y: inclusiveRange(y)) {
      for (int id: range(c)) {
        double x;
        in >> x;
        //x = 2;
        if (x > -0.5) {
          vector<int64_t> row(row_len);
          row[id] = 1;
          row[id + c] = cur_y;
          for (int t: range(cur_y)) {
            row[2 * c + t] = 1;
          }
          known_rows.emplace_back(move(row), log(x));
        }
      }
    }

    size_t initial_rank = reduceToRowEchelonForm(known_rows, row_len);
    known_rows.resize(initial_rank);

    for (int i: range(q)) {
      int id, cur_y;
      in >> id >> cur_y;
      //id = 1;
      //cur_y = 1;
      --id;
      --cur_y;
      vector<int64_t> row(row_len);

      row[id] = 1;
      row[id + c] = cur_y;
      for (int t: range(cur_y)) {
        row[2 * c + t] = 1;
      }
      vector<double> drow(row.begin(), row.end());

      Row new_row{row, 0.0};

      auto copy = known_rows;
      copy.push_back(new_row);
      size_t new_rank = reduceToRowEchelonForm(copy, row_len);

      if (initial_rank == new_rank) {
        double ans = 0;
        for (int i: range(initial_rank)) {
          int first = 0;
          while (known_rows[i].first[first] == 0) {
            ++first;
          }

          auto mult = drow[first] / known_rows[i].first[first];
          for (int j: range(first, row_len)) {
            drow[j] -= mult * known_rows[i].first[j];
          }
          ans += mult * known_rows[i].second;

        }

        out << exp(ans) << "\n";
      }
      else {
        out << -1.0 << "\n";
      }
    }

  }
};


int main() {
  std::ios_base::sync_with_stdio(false);
  TaskC solver;
  std::istream& in(std::cin);
  std::ostream& out(std::cout);
  in.tie(nullptr);
  out << std::fixed;
  out.precision(20);
  solver.solve(in, out);
  return 0;
}