【模板】常用

#define ID isdigit(c = *next++)
#define IS isspace(c = *next++)

struct Istream {
    int size; char *next, buf[20030731];
    Istream & init(FILE *f = stdin) {fread(buf, 1, sizeof buf, f); next = buf; return *this;}
    Istream & operator >> (int &x) {
        int c; x = 0;
        for (; !ID; ) if (!~c) return *this;
        for (x = c & 15; ID; x = x * 10 + (c & 15));
        return *this;
    }
    Istream & operator >> (char *x) {
        int c;
        for (; IS; ) if (!~c) return *this;
        for (*x++ = c; !IS; *x++ = c) if (!~c) break;
        return *x = 0, *this;
    }
    char get() {return *next++;}
} cin;

struct Ostream {
    char *next, buf[20030731], _buf[34];
    Ostream () {next = buf;}
    void flush(FILE *f = stdout) {fwrite(buf, 1, next - buf, f); next = buf;}
    Ostream & operator << (int x) {
        if (!x) return put(48), *this;
        int i;
        for (i = 0; x; x /= 10) _buf[++i] = x % 10 | 48;
        for (; i; --i) put(_buf[i]);
        return *this;
    }
    Ostream & operator << (char c) {return put(c), *this;}
    Ostream & operator << (const char *s) {for (char *p = (char*)s; *p; ++p) put(*p); return *this;}
    void put(char c) {*next++ = c;}
} cout;

inline void up(int &x, const int y) {x < y ? x = y : 0;}
inline void down(int &x, const int y) {x > y ? x = y : 0;}
inline int min(const int x, const int y) {return x < y ? x : y;}
inline int max(const int x, const int y) {return x < y ? y : x;}
inline int & reduce(int &x) {return x += x >> 31 & mod;}
inline void add(int &x, const int y) {x += y - mod, x += x >> 31 & mod;}
inline void sub(int &x, const int y) {x -= y, x += x >> 31 & mod;}
inline int & half(int &x) {return x = (x >> 1) + (-(x & 1) & half_mod);}
inline int & neg(int &x) {return x = (!x - 1) & (mod - x);}

// another ver.
inline bool up(int &x, const int y) {return x < y ? x = y, 1 : 0;}
inline bool down(int &x, const int y) {return x > y ? x = y, 1 : 0;}
inline int & add(int &x, const int y) {return x += y - mod, x += x >> 31 & mod;}
inline int & sub(int &x, const int y) {return x -= y, x += x >> 31 & mod;}

// templated ver.
#define templated template <typename T>
inline bool up(T &x, const T y) {return x < y ? x = y, 1 : 0;}
inline bool down(T &x, const T y) {return x > y ? x = y, 1 : 0;}
inline T min(const T x, const T y) {return x < y ? x : y;}
inline T max(const T x, const T y) {return x < y ? y : x;}

inline ll MulMod(ll a, ll b, ll m){
    ll t = (a * b - (ll)((ld)a * b / m) * m) % m;
    return t + (t >> 63 & m);
}

void build_st_table() {
    int *f, *g = *st, i, j, k = cnt;
    for (j = 0; 1 << j + 1 <= cnt; ++j) {
        f = g; g = st[j + 1]; k -= 1 << j;
        for (i = 0; i < k; ++i)
            g[i] = min(f[i], f[i + (1 << j)]);
    }
}

inline int range(int L, int R) { // [L, R)
    int D = R - L, c = lg2(D);
    return min(st[c][L], st[c][R - (1 << c)]);
}

namespace DC {
    int F[N]; pr D[N];

    int Discretize(int n) {
        int i, cnt = 0;
        std::sort(D, D + n);
        for (i = 0; i < n; ++i)
            F[D[i].second] = (i && D[i].first == D[i - 1].first ? cnt - 1 : (D[cnt] = D[i], cnt++));
        return cnt;
    }
}

class hash_map{
public:
    static const int HASH_MAX = 0xffffff, N = 8000000;
    int cnt, first[HASH_MAX + 2], next[N]; data z[N];
    inline int getHash(int key) {return (key ^ key << 3 ^ key >> 2) & HASH_MAX;}

    void clear() {for(; cnt > 0; --cnt) first[z[cnt].hash] = 0;}

    data * find(int key, bool inserted){
        int x = getHash(key), i;
        for(i = first[x]; i; i = next[i]) if(z[i].key == key) return z + i;
        if(!inserted) return NULL;
        z[++cnt] = data(key, 0, x); next[cnt] = first[x]; first[x] = cnt;
        return z + cnt;
    }
};

// ‘Fast Number Theory Transform‘ is in memos/12.html.
namespace Poly {
    typedef std::complex <double> C;

    const int N = 530000;
    int l, n, rev[N];
    C x[N], y[N], B1[N], B2[N], B3[N];

    void FFT_init(int len) {
        if (l == len) return; n = 1 << (l = len);
        int i; double angle = M_PI;
        for (i = l - 1; i >= 0; angle *= .5, --i) x[1 << i] = C(cos(angle), sin(angle));
        for (i = 3; i < n; ++i) if (i & (i - 1)) x[i] = x[i & -i] * x[i ^ (i & -i)];
        *x = C(1.), *rev = 0;
        for (i = 1; i < n; ++i) rev[i] = rev[i >> 1] >> 1 | (i & 1) << (l - 1);
    }

    void DFT(C *d, C *t) {
        int i, len = 1, delta = n; C *j, *k, R;
        for (i = 0; i < n; ++i) t[rev[i]] = d[i];
        for (i = 0; i < l; ++i) {
            delta >>= 1;
            for (k = x, j = y; j < y + len; k += delta, ++j) *j = *k;
            for (j = t; j < t + n; j += len << 1)
                for (k = j; k < j + len; ++k)
                    R = y[k - j] * k[len], k[len] = *k - R, *k += R;
            len <<= 1;
        }
    }

    void Mul(int degA, int degB, double *a, double *b, double *c) {
        if (!(degA || degB)) {*c = *a * *b; return;}
        FFT_init(lg2(degA + degB) + 1);
        int i; double iv = 1.0 / n;
        for (i = 0; i <= degA; ++i) B1[i] = C(a[i]); std::fill(B1 + i, B1 + n, C());
        for (i = 0; i <= degB; ++i) B2[i] = C(b[i]); std::fill(B2 + i, B2 + n, C());
        DFT(B1, B3), DFT(B2, B1);
        for (i = 0; i < n; ++i) B1[i] *= B3[i];
        DFT(B1, B3), std::reverse(B3 + 1, B3 + n);
        for (i = 0; i <= degA + degB; ++i) c[i] = B3[i].real() * iv;
    }
}

// Gauss elimination with type ‘double‘
struct LnEqn{
    int sz;
    double **m, *b;
    LnEqn (): sz(0) {m = NULL; b = NULL;}
    void resize(int size){
        sz = size; m = new double *[sz];
        for(int i = 0; i < sz; i++){
            m[i] = new double[sz];
            memset(m[i], 0, sz << 3);
        }
        b = new double[sz];
        memset(b, 0, sz << 3);
    }
    ~LnEqn (){
        if(m) {for(int i = 0; i < sz; i++) delete [] (m[i]); delete [] (m);}
        if(b) delete [] (b);
    }
    bool solve(){
        int i, j, k, maxi; double coe;
        for(k = 0; k < sz; k++){
            maxi = k;
            for(i = k + 1; i < sz; i++)
                if(fabs(m[i][k]) > fabs(m[maxi][k]))
                    maxi = i;
            if(fabs(m[maxi][k]) < 1e-8) return false;
            if(maxi != k){
                swap(m[maxi], m[k]);
                swap(b[maxi], b[k]);
            }
            coe = 1.0 / m[k][k];
            for(j = 0; j < sz; j++)
                m[k][j] *= coe;
            m[k][k] = 1.0;
            b[k] *= coe;
            for(i = 0; i < sz; i++){
                if(i == k) continue;
                coe = m[i][k];
                for(j = 0; j < sz; j++)
                    m[i][j] -= coe * m[k][j];
                m[i][k] = 0.0;
                b[i] -= coe * b[k];
            }
        }
        return true;
    }
};

const double eps = 1e-8;

int id[N + E];

double m[E][N], b[E], *c = m[0], ans[N + E];

void pivot(int toFree, int toBasic); // basic(1,m) -> free, free(1,n) -> basic

void pivot(int r, int c){
    int i, j; double coe = 1.0 / m[r][c];
    swap(id[n + r], id[c]);
    m[r][c] = 1.0;
    for(j = 1; j <= n; ++j)
        m[r][j] *= coe;
    b[r] *= coe;
    for(i = 0; i <= e; ++i){
        if(i == r) continue;
        coe = m[i][c]; m[i][c] = 0.0;
        for(j = 1; j <= n; ++j)
            m[i][j] -= coe * m[r][j];
        b[i] -= coe * b[r];
    }
}

bool simplex(){
    int i, j, basic, free; double G;
    for(; ; ){
        basic = free = 0; G = INFINITY;
        for(i = 1; i <= n; ++i) // free (nonbasic) variable
            if(c[i] > c[free]) free = i;
        if(!free) return true;
        for(j = 1; j <= e; ++j) // basic variable
            if(m[j][free] > eps && b[j] < G * m[j][free]){
                G = b[j] / m[j][free]; basic = j;
            }
        if(!basic) return puts("Unbounded"), false;
        pivot(basic, free);
    }
}

// initialize :
for(j = 1; j <= n + e; ++j) id[j] = j;
c[0] = eps;

// output:
if(simplex()){
    printf("%.8lg\n", -*b);
    for(i = 1; i <= e; ++i) ans[id[n + i]] = b[i];
    for(j = 1; j <= n; ++j) printf("%.8lg%c", ans[j], j == n ? 10 : 32);
}

int gauss(int n) {
    int i, j, k, det = 1; ll coe;
    static int *m[N];
    for (i = 0; i < n; ++i) m[i] = mat[i];
    for (i = 0; i < n; ++i) {
        for (j = i; j < n && !m[j][i]; ++j);
        if (j == n) return 0;
        if (i != j) det = mod - det, std::swap(m[i], m[j]);
        det = (ll)det * m[i][i] % mod;
        coe = PowerMod(m[i][i], mod - 2);
        for (j = 0; j < n; ++j) m[i][j] = m[i][j] * coe % mod;
        for (k = i + 1; k < n; ++k)
            for (coe = mod - m[k][i], j = i; j < n; ++j) m[k][j] = (m[k][j] + coe * m[i][j]) % mod;
    }
    return det;
}