/*
 * Copyright (c) 2013, Ambroz Bizjak <ambrop7@gmail.com>
 * All rights reserved.
 * 
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met: 
 * 
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer. 
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution. 
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
 
/*
 * Asymptotically efficient computation of a weighted sum of last N received numbers,
 * with worst case time complexity O(sqrt(N*log(N))) per received number.
 * See: http://c...content-available-to-author-only...e.com/questions/10612/weighted-sum-of-last-n-numbers/10670#10670
 * 
 * Initialization:
 *   WSum state;
 *   WSum_Init(&state, number_of_weights, pointer_to_weights);
 * The weights need to be specified in "reverse", meaning that the first weight
 * applies to the most recent number.
 * 
 * Operation:
 *   double current_wsum = WSum_AddNumber(&state, your_number);
 */
 
#include <stdio.h>
#include <complex.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <stdlib.h>
 
/*
 * Probably not the fastest FFT.
 */
static void fft_rec (int exp, double complex *data, double complex *temp, double complex rou)
{
    if (exp == 0) {
        return;
    }
    size_t half = ((size_t)1) << (exp - 1);
    for (size_t i = 0; i < half; i++) {
        temp[i] = data[2 * i];
        temp[half + i] = data[2 * i + 1];
    }
    double complex rou_2 = rou * rou;
    fft_rec(exp - 1, temp, data, rou_2);
    fft_rec(exp - 1, temp + half, data, rou_2);
    double complex f = 1;
    for (size_t i = 0; i < half; i++) {
        data[i] = temp[i] + f * temp[half + i];
        data[half + i] = temp[i] - f * temp[half + i];
        f *= rou;
    }
}
 
static void fft (int exp, double complex *data, double complex *temp)
{
    size_t n = ((size_t)1) << exp;
    double complex rou = cexp(I * (-2 * 3.14159265358979323846 / n));
    fft_rec(exp, data, temp, rou);
}
 
// Input: m weights, 2m-1 numbers.
// Output: m sums, sums[i] = sum (k=0 to m-1) weights[k]*numbers[i+m-1-k]
// This function implements everything FFT-related that the WSum code needs.
static void compute_sums (size_t m, const double *weights, const double *numbers, double *sums)
{
    int dft_exp = 0;
    size_t dft_n = 1;
    while (dft_n < 2 * m) {
        dft_exp++;
        dft_n *= 2;
    }
 
    double complex *a = malloc(dft_n * sizeof(*a));
    double complex *b = malloc(dft_n * sizeof(*b));
    double complex *c = malloc(dft_n * sizeof(*c));
 
    for (size_t i = 0; i < m; i++) {
        a[i] = weights[i];
    }
    for (size_t i = m; i < dft_n; i++) {
        a[i] = 0;
    }
 
    for (size_t i = 0; i < 2 * m - 1; i++) {
        b[i] = numbers[i];
    }
    for (size_t i = 2 * m - 1; i < dft_n; i++) {
        b[i] = 0;
    }
 
    fft(dft_exp, a, c);
    fft(dft_exp, b, c);
 
    for (size_t i = 0; i < dft_n; i++) {
        c[i] = conj(a[i] * b[i]);
    }
 
    fft(dft_exp, c, a);
 
    for (size_t i = 0; i < m; i++) {
        sums[i] = creal(c[m - 1 + i]) / dft_n;
    }
 
    free(c);
    free(b);
    free(a);
}
 
typedef struct {
    size_t n;
    size_t pos;
    double *weights;
    double *inputs;
} SimpleWSum;
 
void SimpleWSum_Init (SimpleWSum *o, size_t n, const double *weights)
{
    assert(n > 0);
 
    o->n = n;
    o->pos = n - 1;
    o->weights = malloc(n * sizeof(o->weights[0]));
    for (size_t i = 0; i < n; i++) {
        o->weights[i] = weights[i];
    }
    o->inputs = malloc(n * sizeof(o->inputs[0]));
    for (size_t i = 0; i < n; i++) {
        o->inputs[i] = 0;
    }
}
 
void SimpleWSum_Free (SimpleWSum *o)
{
    free(o->inputs);
    free(o->weights);
}
 
double SimpleWSum_AddNumber (SimpleWSum *o, double x)
{
    o->inputs[o->pos] = x;
 
    size_t n_minus_pos = o->n - o->pos;
    double sum = 0;
    for (size_t i = 0; i < n_minus_pos; i++) {
        sum += o->weights[i] * o->inputs[o->pos + i];
    }
    for (size_t i = n_minus_pos; i < o->n; i++) {
        sum += o->weights[i] * o->inputs[i - n_minus_pos];
    }
 
    if (o->pos == 0) {
        o->pos = o->n - 1;
    } else {
        o->pos--;
    }
 
    return sum;
}
 
typedef struct {
    size_t m;
    size_t n;
    size_t num_m;
    double *weights;
    double *inputs;
    double *sums;
    double *new_sums;
    size_t m_pos;
    size_t num_m_pos;
} WSum;
 
void WSum_Init (WSum *o, size_t orig_n, const double *orig_weights)
{
    assert(orig_n > 0);
 
    // choose m
    o->m = sqrt(orig_n * log(orig_n));
    if (o->m == 0) {
        o->m = 1;
    } else if (o->m > orig_n) {
        o->m = orig_n;
    }
 
    // choose n, taking care it's a multiple of m
    o->n = orig_n;
    if (o->n % o->m) {
        o->n += o->m - (o->n % o->m);
    }
    o->num_m = o->n / o->m;
 
    // allocate and initialize weights
    o->weights = malloc((o->n + 1) * sizeof(o->weights[0]));
    for (size_t i = 0; i < orig_n; i++) {
        o->weights[i] = orig_weights[i];
    }
    for (size_t i = orig_n; i < o->n + 1; i++) {
        o->weights[i] = 0;
    }
 
    // allocate and initialize inputs
    o->inputs = malloc((o->n + o->m) * sizeof(o->inputs[0]));
    for (size_t i = 0; i < o->n; i++) {
        o->inputs[i] = 0;
    }
 
    // allocate and initialize sums
    o->sums = malloc(o->n * sizeof(o->sums[0]));
    o->new_sums = malloc(o->n * sizeof(o->new_sums[0]));
    for (size_t i = 0; i < o->n; i++) {
        o->sums[i] = 0;
        o->new_sums[i] = 0;
    }
 
    o->m_pos = 0;
    o->num_m_pos = 0;
}
 
void WSum_Free (WSum *o)
{
    free(o->new_sums);
    free(o->sums);
    free(o->inputs);
    free(o->weights);
}
 
double WSum_AddNumber (WSum *o, double x)
{
    assert(o->m_pos < o->m);
 
    double sum = 0;
 
    // add from sums
    for (size_t i = 0; i < o->num_m; i++) {
        sum += o->sums[i * o->m + o->m_pos];
    }
 
    // push this number to inputs
    o->inputs[o->n + o->m_pos] = x;
    o->m_pos++;
 
    // add from pending inputs (incl. this number)
    for (size_t i = 0; i < o->m_pos; i++) {
        sum += o->inputs[o->n + i] * o->weights[o->m_pos - 1 - i];
    }
 
    // decide up to which block we would like to update
    size_t end_update_pos = o->num_m * ((double)o->m_pos / o->m);
    if (o->m_pos == o->m) {
        end_update_pos = o->num_m; // don't assume any FP precision
    }
 
    // cannot compute last two new_sums blocks unless we have the new inputs
    if (o->m_pos < o->m) {
        while (end_update_pos > 0 && end_update_pos >= o->num_m - 1) {
            end_update_pos--;
        }
    }
 
    // compute new_sums blocks
    while (o->num_m_pos < end_update_pos) {
        if (o->num_m_pos == o->num_m - 1) {
            for (size_t j = 0; j < o->m; j++) {
                o->inputs[o->n + j] = 0;
            }
        } else {
            memcpy(o->inputs + (o->num_m_pos + 1) * o->m, o->inputs + (o->num_m_pos + 2) * o->m, o->m * sizeof(o->inputs[0]));
        }
        compute_sums(o->m, o->weights + (o->num_m - 1 - o->num_m_pos) * o->m + 1, o->inputs + o->num_m_pos * o->m, o->new_sums + o->num_m_pos * o->m);
        o->num_m_pos++;
    }
 
    // switch at end of input block
    if (o->m_pos == o->m) {
        assert(o->num_m_pos == o->num_m);
        memcpy(o->inputs, o->inputs + o->m, o->m * sizeof(o->inputs[0]));
        double *old_sums = o->sums;
        o->sums = o->new_sums;
        o->new_sums = old_sums;
        o->m_pos = 0;
        o->num_m_pos = 0;
    }
 
    return sum;
}
 
static void test (size_t n, int num_iter, int use_wsum, int use_simple_wsum)
{
    double *w = malloc(n * sizeof(w[0]));
    for (size_t i = 0; i < n; i++) {
        w[i] = (double)rand() / RAND_MAX;
    }
 
    double dummy = 0; // prevent optimizations
 
    WSum state;
    SimpleWSum simple_state;
    if (use_wsum) {
        WSum_Init(&state, n, w);
        printf("WSum m=%zu n=%zu num_m=%zu\n", state.m, state.n, state.num_m);
    }
    if (use_simple_wsum) {
        SimpleWSum_Init(&simple_state, n, w);
    }
 
    for (int k = 0; k < num_iter; k++) {
        double x = rand();
 
        double sum = 0;
        double simple_sum = 0;
        if (use_wsum) {
            sum = WSum_AddNumber(&state, x);
        }
        if (use_simple_wsum) {
            simple_sum = SimpleWSum_AddNumber(&simple_state, x);
        }
 
        if (use_wsum && use_simple_wsum) {
            int failed = (fabs(sum - simple_sum) > 0.02);
            printf("%f %f %f %s\n", x, sum, simple_sum, (failed ? "FAILED" : "OK"));
            if (failed) {
                exit(1);
            }
        }
 
        dummy += sum + simple_sum;
    }
 
    printf("dummy %f\n", dummy);
 
    if (use_simple_wsum) {
        SimpleWSum_Free(&simple_state);
    }
    if (use_wsum) {
        WSum_Free(&state);
    }
    free(w);
}
 
// Around here it WSum starts to beat SimpleWSum
//#define NUM_WEIGHTS 17000
//#define NUM_ITER 80000
#define NUM_WEIGHTS 50
#define NUM_ITER 200
 
int main ()
{
    // Run both, comparing the results.
    test(NUM_WEIGHTS, NUM_ITER, 1, 1);
 
    // Run just WSum or just SimpleWSum for benchmarking
    //test(NUM_WEIGHTS, NUM_ITER, 1, 0); // WSum
    //test(NUM_WEIGHTS, NUM_ITER, 0, 1); // SimpleWSum
}
 

LyoKICogQ29weXJpZ2h0IChjKSAyMDEzLCBBbWJyb3ogQml6amFrIDxhbWJyb3A3QGdtYWlsLmNvbT4KICogQWxsIHJpZ2h0cyByZXNlcnZlZC4KICogCiAqIFJlZGlzdHJpYnV0aW9uIGFuZCB1c2UgaW4gc291cmNlIGFuZCBiaW5hcnkgZm9ybXMsIHdpdGggb3Igd2l0aG91dAogKiBtb2RpZmljYXRpb24sIGFyZSBwZXJtaXR0ZWQgcHJvdmlkZWQgdGhhdCB0aGUgZm9sbG93aW5nIGNvbmRpdGlvbnMgYXJlIG1ldDogCiAqIAogKiAxLiBSZWRpc3RyaWJ1dGlvbnMgb2Ygc291cmNlIGNvZGUgbXVzdCByZXRhaW4gdGhlIGFib3ZlIGNvcHlyaWdodCBub3RpY2UsIHRoaXMKICogICAgbGlzdCBvZiBjb25kaXRpb25zIGFuZCB0aGUgZm9sbG93aW5nIGRpc2NsYWltZXIuIAogKiAyLiBSZWRpc3RyaWJ1dGlvbnMgaW4gYmluYXJ5IGZvcm0gbXVzdCByZXByb2R1Y2UgdGhlIGFib3ZlIGNvcHlyaWdodCBub3RpY2UsCiAqICAgIHRoaXMgbGlzdCBvZiBjb25kaXRpb25zIGFuZCB0aGUgZm9sbG93aW5nIGRpc2NsYWltZXIgaW4gdGhlIGRvY3VtZW50YXRpb24KICogICAgYW5kL29yIG90aGVyIG1hdGVyaWFscyBwcm92aWRlZCB3aXRoIHRoZSBkaXN0cmlidXRpb24uIAogKiAKICogVEhJUyBTT0ZUV0FSRSBJUyBQUk9WSURFRCBCWSBUSEUgQ09QWVJJR0hUIEhPTERFUlMgQU5EIENPTlRSSUJVVE9SUyAiQVMgSVMiIEFORAogKiBBTlkgRVhQUkVTUyBPUiBJTVBMSUVEIFdBUlJBTlRJRVMsIElOQ0xVRElORywgQlVUIE5PVCBMSU1JVEVEIFRPLCBUSEUgSU1QTElFRAogKiBXQVJSQU5USUVTIE9GIE1FUkNIQU5UQUJJTElUWSBBTkQgRklUTkVTUyBGT1IgQSBQQVJUSUNVTEFSIFBVUlBPU0UgQVJFCiAqIERJU0NMQUlNRUQuIElOIE5PIEVWRU5UIFNIQUxMIFRIRSBDT1BZUklHSFQgT1dORVIgT1IgQ09OVFJJQlVUT1JTIEJFIExJQUJMRSBGT1IKICogQU5ZIERJUkVDVCwgSU5ESVJFQ1QsIElOQ0lERU5UQUwsIFNQRUNJQUwsIEVYRU1QTEFSWSwgT1IgQ09OU0VRVUVOVElBTCBEQU1BR0VTCiAqIChJTkNMVURJTkcsIEJVVCBOT1QgTElNSVRFRCBUTywgUFJPQ1VSRU1FTlQgT0YgU1VCU1RJVFVURSBHT09EUyBPUiBTRVJWSUNFUzsKICogTE9TUyBPRiBVU0UsIERBVEEsIE9SIFBST0ZJVFM7IE9SIEJVU0lORVNTIElOVEVSUlVQVElPTikgSE9XRVZFUiBDQVVTRUQgQU5ECiAqIE9OIEFOWSBUSEVPUlkgT0YgTElBQklMSVRZLCBXSEVUSEVSIElOIENPTlRSQUNULCBTVFJJQ1QgTElBQklMSVRZLCBPUiBUT1JUCiAqIChJTkNMVURJTkcgTkVHTElHRU5DRSBPUiBPVEhFUldJU0UpIEFSSVNJTkcgSU4gQU5ZIFdBWSBPVVQgT0YgVEhFIFVTRSBPRiBUSElTCiAqIFNPRlRXQVJFLCBFVkVOIElGIEFEVklTRUQgT0YgVEhFIFBPU1NJQklMSVRZIE9GIFNVQ0ggREFNQUdFLgogKi8KCi8qCiAqIEFzeW1wdG90aWNhbGx5IGVmZmljaWVudCBjb21wdXRhdGlvbiBvZiBhIHdlaWdodGVkIHN1bSBvZiBsYXN0IE4gcmVjZWl2ZWQgbnVtYmVycywKICogd2l0aCB3b3JzdCBjYXNlIHRpbWUgY29tcGxleGl0eSBPKHNxcnQoTipsb2coTikpKSBwZXIgcmVjZWl2ZWQgbnVtYmVyLgogKiBTZWU6IGh0dHA6Ly9jLi4uY29udGVudC1hdmFpbGFibGUtdG8tYXV0aG9yLW9ubHkuLi5lLmNvbS9xdWVzdGlvbnMvMTA2MTIvd2VpZ2h0ZWQtc3VtLW9mLWxhc3Qtbi1udW1iZXJzLzEwNjcwIzEwNjcwCiAqIAogKiBJbml0aWFsaXphdGlvbjoKICogICBXU3VtIHN0YXRlOwogKiAgIFdTdW1fSW5pdCgmc3RhdGUsIG51bWJlcl9vZl93ZWlnaHRzLCBwb2ludGVyX3RvX3dlaWdodHMpOwogKiBUaGUgd2VpZ2h0cyBuZWVkIHRvIGJlIHNwZWNpZmllZCBpbiAicmV2ZXJzZSIsIG1lYW5pbmcgdGhhdCB0aGUgZmlyc3Qgd2VpZ2h0CiAqIGFwcGxpZXMgdG8gdGhlIG1vc3QgcmVjZW50IG51bWJlci4KICogCiAqIE9wZXJhdGlvbjoKICogICBkb3VibGUgY3VycmVudF93c3VtID0gV1N1bV9BZGROdW1iZXIoJnN0YXRlLCB5b3VyX251bWJlcik7CiAqLwoKI2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxjb21wbGV4Lmg+CiNpbmNsdWRlIDxzdHJpbmcuaD4KI2luY2x1ZGUgPG1hdGguaD4KI2luY2x1ZGUgPGFzc2VydC5oPgojaW5jbHVkZSA8c3RkbGliLmg+CgovKgogKiBQcm9iYWJseSBub3QgdGhlIGZhc3Rlc3QgRkZULgogKi8Kc3RhdGljIHZvaWQgZmZ0X3JlYyAoaW50IGV4cCwgZG91YmxlIGNvbXBsZXggKmRhdGEsIGRvdWJsZSBjb21wbGV4ICp0ZW1wLCBkb3VibGUgY29tcGxleCByb3UpCnsKICAgIGlmIChleHAgPT0gMCkgewogICAgICAgIHJldHVybjsKICAgIH0KICAgIHNpemVfdCBoYWxmID0gKChzaXplX3QpMSkgPDwgKGV4cCAtIDEpOwogICAgZm9yIChzaXplX3QgaSA9IDA7IGkgPCBoYWxmOyBpKyspIHsKICAgICAgICB0ZW1wW2ldID0gZGF0YVsyICogaV07CiAgICAgICAgdGVtcFtoYWxmICsgaV0gPSBkYXRhWzIgKiBpICsgMV07CiAgICB9CiAgICBkb3VibGUgY29tcGxleCByb3VfMiA9IHJvdSAqIHJvdTsKICAgIGZmdF9yZWMoZXhwIC0gMSwgdGVtcCwgZGF0YSwgcm91XzIpOwogICAgZmZ0X3JlYyhleHAgLSAxLCB0ZW1wICsgaGFsZiwgZGF0YSwgcm91XzIpOwogICAgZG91YmxlIGNvbXBsZXggZiA9IDE7CiAgICBmb3IgKHNpemVfdCBpID0gMDsgaSA8IGhhbGY7IGkrKykgewogICAgICAgIGRhdGFbaV0gPSB0ZW1wW2ldICsgZiAqIHRlbXBbaGFsZiArIGldOwogICAgICAgIGRhdGFbaGFsZiArIGldID0gdGVtcFtpXSAtIGYgKiB0ZW1wW2hhbGYgKyBpXTsKICAgICAgICBmICo9IHJvdTsKICAgIH0KfQoKc3RhdGljIHZvaWQgZmZ0IChpbnQgZXhwLCBkb3VibGUgY29tcGxleCAqZGF0YSwgZG91YmxlIGNvbXBsZXggKnRlbXApCnsKICAgIHNpemVfdCBuID0gKChzaXplX3QpMSkgPDwgZXhwOwogICAgZG91YmxlIGNvbXBsZXggcm91ID0gY2V4cChJICogKC0yICogMy4xNDE1OTI2NTM1ODk3OTMyMzg0NiAvIG4pKTsKICAgIGZmdF9yZWMoZXhwLCBkYXRhLCB0ZW1wLCByb3UpOwp9CgovLyBJbnB1dDogbSB3ZWlnaHRzLCAybS0xIG51bWJlcnMuCi8vIE91dHB1dDogbSBzdW1zLCBzdW1zW2ldID0gc3VtIChrPTAgdG8gbS0xKSB3ZWlnaHRzW2tdKm51bWJlcnNbaSttLTEta10KLy8gVGhpcyBmdW5jdGlvbiBpbXBsZW1lbnRzIGV2ZXJ5dGhpbmcgRkZULXJlbGF0ZWQgdGhhdCB0aGUgV1N1bSBjb2RlIG5lZWRzLgpzdGF0aWMgdm9pZCBjb21wdXRlX3N1bXMgKHNpemVfdCBtLCBjb25zdCBkb3VibGUgKndlaWdodHMsIGNvbnN0IGRvdWJsZSAqbnVtYmVycywgZG91YmxlICpzdW1zKQp7CiAgICBpbnQgZGZ0X2V4cCA9IDA7CiAgICBzaXplX3QgZGZ0X24gPSAxOwogICAgd2hpbGUgKGRmdF9uIDwgMiAqIG0pIHsKICAgICAgICBkZnRfZXhwKys7CiAgICAgICAgZGZ0X24gKj0gMjsKICAgIH0KICAgIAogICAgZG91YmxlIGNvbXBsZXggKmEgPSBtYWxsb2MoZGZ0X24gKiBzaXplb2YoKmEpKTsKICAgIGRvdWJsZSBjb21wbGV4ICpiID0gbWFsbG9jKGRmdF9uICogc2l6ZW9mKCpiKSk7CiAgICBkb3VibGUgY29tcGxleCAqYyA9IG1hbGxvYyhkZnRfbiAqIHNpemVvZigqYykpOwogICAgCiAgICBmb3IgKHNpemVfdCBpID0gMDsgaSA8IG07IGkrKykgewogICAgICAgIGFbaV0gPSB3ZWlnaHRzW2ldOwogICAgfQogICAgZm9yIChzaXplX3QgaSA9IG07IGkgPCBkZnRfbjsgaSsrKSB7CiAgICAgICAgYVtpXSA9IDA7CiAgICB9CiAgICAKICAgIGZvciAoc2l6ZV90IGkgPSAwOyBpIDwgMiAqIG0gLSAxOyBpKyspIHsKICAgICAgICBiW2ldID0gbnVtYmVyc1tpXTsKICAgIH0KICAgIGZvciAoc2l6ZV90IGkgPSAyICogbSAtIDE7IGkgPCBkZnRfbjsgaSsrKSB7CiAgICAgICAgYltpXSA9IDA7CiAgICB9CiAgICAKICAgIGZmdChkZnRfZXhwLCBhLCBjKTsKICAgIGZmdChkZnRfZXhwLCBiLCBjKTsKICAgIAogICAgZm9yIChzaXplX3QgaSA9IDA7IGkgPCBkZnRfbjsgaSsrKSB7CiAgICAgICAgY1tpXSA9IGNvbmooYVtpXSAqIGJbaV0pOwogICAgfQogICAgCiAgICBmZnQoZGZ0X2V4cCwgYywgYSk7CiAgICAKICAgIGZvciAoc2l6ZV90IGkgPSAwOyBpIDwgbTsgaSsrKSB7CiAgICAgICAgc3Vtc1tpXSA9IGNyZWFsKGNbbSAtIDEgKyBpXSkgLyBkZnRfbjsKICAgIH0KICAgIAogICAgZnJlZShjKTsKICAgIGZyZWUoYik7CiAgICBmcmVlKGEpOwp9Cgp0eXBlZGVmIHN0cnVjdCB7CiAgICBzaXplX3QgbjsKICAgIHNpemVfdCBwb3M7CiAgICBkb3VibGUgKndlaWdodHM7CiAgICBkb3VibGUgKmlucHV0czsKfSBTaW1wbGVXU3VtOwoKdm9pZCBTaW1wbGVXU3VtX0luaXQgKFNpbXBsZVdTdW0gKm8sIHNpemVfdCBuLCBjb25zdCBkb3VibGUgKndlaWdodHMpCnsKICAgIGFzc2VydChuID4gMCk7CiAgICAKICAgIG8tPm4gPSBuOwogICAgby0+cG9zID0gbiAtIDE7CiAgICBvLT53ZWlnaHRzID0gbWFsbG9jKG4gKiBzaXplb2Yoby0+d2VpZ2h0c1swXSkpOwogICAgZm9yIChzaXplX3QgaSA9IDA7IGkgPCBuOyBpKyspIHsKICAgICAgICBvLT53ZWlnaHRzW2ldID0gd2VpZ2h0c1tpXTsKICAgIH0KICAgIG8tPmlucHV0cyA9IG1hbGxvYyhuICogc2l6ZW9mKG8tPmlucHV0c1swXSkpOwogICAgZm9yIChzaXplX3QgaSA9IDA7IGkgPCBuOyBpKyspIHsKICAgICAgICBvLT5pbnB1dHNbaV0gPSAwOwogICAgfQp9Cgp2b2lkIFNpbXBsZVdTdW1fRnJlZSAoU2ltcGxlV1N1bSAqbykKewogICAgZnJlZShvLT5pbnB1dHMpOwogICAgZnJlZShvLT53ZWlnaHRzKTsKfQoKZG91YmxlIFNpbXBsZVdTdW1fQWRkTnVtYmVyIChTaW1wbGVXU3VtICpvLCBkb3VibGUgeCkKewogICAgby0+aW5wdXRzW28tPnBvc10gPSB4OwogICAgCiAgICBzaXplX3Qgbl9taW51c19wb3MgPSBvLT5uIC0gby0+cG9zOwogICAgZG91YmxlIHN1bSA9IDA7CiAgICBmb3IgKHNpemVfdCBpID0gMDsgaSA8IG5fbWludXNfcG9zOyBpKyspIHsKICAgICAgICBzdW0gKz0gby0+d2VpZ2h0c1tpXSAqIG8tPmlucHV0c1tvLT5wb3MgKyBpXTsKICAgIH0KICAgIGZvciAoc2l6ZV90IGkgPSBuX21pbnVzX3BvczsgaSA8IG8tPm47IGkrKykgewogICAgICAgIHN1bSArPSBvLT53ZWlnaHRzW2ldICogby0+aW5wdXRzW2kgLSBuX21pbnVzX3Bvc107CiAgICB9CiAgICAKICAgIGlmIChvLT5wb3MgPT0gMCkgewogICAgICAgIG8tPnBvcyA9IG8tPm4gLSAxOwogICAgfSBlbHNlIHsKICAgICAgICBvLT5wb3MtLTsKICAgIH0KICAgIAogICAgcmV0dXJuIHN1bTsKfQoKdHlwZWRlZiBzdHJ1Y3QgewogICAgc2l6ZV90IG07CiAgICBzaXplX3QgbjsKICAgIHNpemVfdCBudW1fbTsKICAgIGRvdWJsZSAqd2VpZ2h0czsKICAgIGRvdWJsZSAqaW5wdXRzOwogICAgZG91YmxlICpzdW1zOwogICAgZG91YmxlICpuZXdfc3VtczsKICAgIHNpemVfdCBtX3BvczsKICAgIHNpemVfdCBudW1fbV9wb3M7Cn0gV1N1bTsKCnZvaWQgV1N1bV9Jbml0IChXU3VtICpvLCBzaXplX3Qgb3JpZ19uLCBjb25zdCBkb3VibGUgKm9yaWdfd2VpZ2h0cykKewogICAgYXNzZXJ0KG9yaWdfbiA+IDApOwogICAgCiAgICAvLyBjaG9vc2UgbQogICAgby0+bSA9IHNxcnQob3JpZ19uICogbG9nKG9yaWdfbikpOwogICAgaWYgKG8tPm0gPT0gMCkgewogICAgICAgIG8tPm0gPSAxOwogICAgfSBlbHNlIGlmIChvLT5tID4gb3JpZ19uKSB7CiAgICAgICAgby0+bSA9IG9yaWdfbjsKICAgIH0KICAgIAogICAgLy8gY2hvb3NlIG4sIHRha2luZyBjYXJlIGl0J3MgYSBtdWx0aXBsZSBvZiBtCiAgICBvLT5uID0gb3JpZ19uOwogICAgaWYgKG8tPm4gJSBvLT5tKSB7CiAgICAgICAgby0+biArPSBvLT5tIC0gKG8tPm4gJSBvLT5tKTsKICAgIH0KICAgIG8tPm51bV9tID0gby0+biAvIG8tPm07CiAgICAKICAgIC8vIGFsbG9jYXRlIGFuZCBpbml0aWFsaXplIHdlaWdodHMKICAgIG8tPndlaWdodHMgPSBtYWxsb2MoKG8tPm4gKyAxKSAqIHNpemVvZihvLT53ZWlnaHRzWzBdKSk7CiAgICBmb3IgKHNpemVfdCBpID0gMDsgaSA8IG9yaWdfbjsgaSsrKSB7CiAgICAgICAgby0+d2VpZ2h0c1tpXSA9IG9yaWdfd2VpZ2h0c1tpXTsKICAgIH0KICAgIGZvciAoc2l6ZV90IGkgPSBvcmlnX247IGkgPCBvLT5uICsgMTsgaSsrKSB7CiAgICAgICAgby0+d2VpZ2h0c1tpXSA9IDA7CiAgICB9CiAgICAKICAgIC8vIGFsbG9jYXRlIGFuZCBpbml0aWFsaXplIGlucHV0cwogICAgby0+aW5wdXRzID0gbWFsbG9jKChvLT5uICsgby0+bSkgKiBzaXplb2Yoby0+aW5wdXRzWzBdKSk7CiAgICBmb3IgKHNpemVfdCBpID0gMDsgaSA8IG8tPm47IGkrKykgewogICAgICAgIG8tPmlucHV0c1tpXSA9IDA7CiAgICB9CiAgICAKICAgIC8vIGFsbG9jYXRlIGFuZCBpbml0aWFsaXplIHN1bXMKICAgIG8tPnN1bXMgPSBtYWxsb2Moby0+biAqIHNpemVvZihvLT5zdW1zWzBdKSk7CiAgICBvLT5uZXdfc3VtcyA9IG1hbGxvYyhvLT5uICogc2l6ZW9mKG8tPm5ld19zdW1zWzBdKSk7CiAgICBmb3IgKHNpemVfdCBpID0gMDsgaSA8IG8tPm47IGkrKykgewogICAgICAgIG8tPnN1bXNbaV0gPSAwOwogICAgICAgIG8tPm5ld19zdW1zW2ldID0gMDsKICAgIH0KICAgIAogICAgby0+bV9wb3MgPSAwOwogICAgby0+bnVtX21fcG9zID0gMDsKfQoKdm9pZCBXU3VtX0ZyZWUgKFdTdW0gKm8pCnsKICAgIGZyZWUoby0+bmV3X3N1bXMpOwogICAgZnJlZShvLT5zdW1zKTsKICAgIGZyZWUoby0+aW5wdXRzKTsKICAgIGZyZWUoby0+d2VpZ2h0cyk7Cn0KCmRvdWJsZSBXU3VtX0FkZE51bWJlciAoV1N1bSAqbywgZG91YmxlIHgpCnsKICAgIGFzc2VydChvLT5tX3BvcyA8IG8tPm0pOwogICAgCiAgICBkb3VibGUgc3VtID0gMDsKICAgIAogICAgLy8gYWRkIGZyb20gc3VtcwogICAgZm9yIChzaXplX3QgaSA9IDA7IGkgPCBvLT5udW1fbTsgaSsrKSB7CiAgICAgICAgc3VtICs9IG8tPnN1bXNbaSAqIG8tPm0gKyBvLT5tX3Bvc107CiAgICB9CiAgICAKICAgIC8vIHB1c2ggdGhpcyBudW1iZXIgdG8gaW5wdXRzCiAgICBvLT5pbnB1dHNbby0+biArIG8tPm1fcG9zXSA9IHg7CiAgICBvLT5tX3BvcysrOwogICAgCiAgICAvLyBhZGQgZnJvbSBwZW5kaW5nIGlucHV0cyAoaW5jbC4gdGhpcyBudW1iZXIpCiAgICBmb3IgKHNpemVfdCBpID0gMDsgaSA8IG8tPm1fcG9zOyBpKyspIHsKICAgICAgICBzdW0gKz0gby0+aW5wdXRzW28tPm4gKyBpXSAqIG8tPndlaWdodHNbby0+bV9wb3MgLSAxIC0gaV07CiAgICB9CiAgICAKICAgIC8vIGRlY2lkZSB1cCB0byB3aGljaCBibG9jayB3ZSB3b3VsZCBsaWtlIHRvIHVwZGF0ZQogICAgc2l6ZV90IGVuZF91cGRhdGVfcG9zID0gby0+bnVtX20gKiAoKGRvdWJsZSlvLT5tX3BvcyAvIG8tPm0pOwogICAgaWYgKG8tPm1fcG9zID09IG8tPm0pIHsKICAgICAgICBlbmRfdXBkYXRlX3BvcyA9IG8tPm51bV9tOyAvLyBkb24ndCBhc3N1bWUgYW55IEZQIHByZWNpc2lvbgogICAgfQogICAgCiAgICAvLyBjYW5ub3QgY29tcHV0ZSBsYXN0IHR3byBuZXdfc3VtcyBibG9ja3MgdW5sZXNzIHdlIGhhdmUgdGhlIG5ldyBpbnB1dHMKICAgIGlmIChvLT5tX3BvcyA8IG8tPm0pIHsKICAgICAgICB3aGlsZSAoZW5kX3VwZGF0ZV9wb3MgPiAwICYmIGVuZF91cGRhdGVfcG9zID49IG8tPm51bV9tIC0gMSkgewogICAgICAgICAgICBlbmRfdXBkYXRlX3Bvcy0tOwogICAgICAgIH0KICAgIH0KICAgIAogICAgLy8gY29tcHV0ZSBuZXdfc3VtcyBibG9ja3MKICAgIHdoaWxlIChvLT5udW1fbV9wb3MgPCBlbmRfdXBkYXRlX3BvcykgewogICAgICAgIGlmIChvLT5udW1fbV9wb3MgPT0gby0+bnVtX20gLSAxKSB7CiAgICAgICAgICAgIGZvciAoc2l6ZV90IGogPSAwOyBqIDwgby0+bTsgaisrKSB7CiAgICAgICAgICAgICAgICBvLT5pbnB1dHNbby0+biArIGpdID0gMDsKICAgICAgICAgICAgfQogICAgICAgIH0gZWxzZSB7CiAgICAgICAgICAgIG1lbWNweShvLT5pbnB1dHMgKyAoby0+bnVtX21fcG9zICsgMSkgKiBvLT5tLCBvLT5pbnB1dHMgKyAoby0+bnVtX21fcG9zICsgMikgKiBvLT5tLCBvLT5tICogc2l6ZW9mKG8tPmlucHV0c1swXSkpOwogICAgICAgIH0KICAgICAgICBjb21wdXRlX3N1bXMoby0+bSwgby0+d2VpZ2h0cyArIChvLT5udW1fbSAtIDEgLSBvLT5udW1fbV9wb3MpICogby0+bSArIDEsIG8tPmlucHV0cyArIG8tPm51bV9tX3BvcyAqIG8tPm0sIG8tPm5ld19zdW1zICsgby0+bnVtX21fcG9zICogby0+bSk7CiAgICAgICAgby0+bnVtX21fcG9zKys7CiAgICB9CiAgICAKICAgIC8vIHN3aXRjaCBhdCBlbmQgb2YgaW5wdXQgYmxvY2sKICAgIGlmIChvLT5tX3BvcyA9PSBvLT5tKSB7CiAgICAgICAgYXNzZXJ0KG8tPm51bV9tX3BvcyA9PSBvLT5udW1fbSk7CiAgICAgICAgbWVtY3B5KG8tPmlucHV0cywgby0+aW5wdXRzICsgby0+bSwgby0+bSAqIHNpemVvZihvLT5pbnB1dHNbMF0pKTsKICAgICAgICBkb3VibGUgKm9sZF9zdW1zID0gby0+c3VtczsKICAgICAgICBvLT5zdW1zID0gby0+bmV3X3N1bXM7CiAgICAgICAgby0+bmV3X3N1bXMgPSBvbGRfc3VtczsKICAgICAgICBvLT5tX3BvcyA9IDA7CiAgICAgICAgby0+bnVtX21fcG9zID0gMDsKICAgIH0KICAgIAogICAgcmV0dXJuIHN1bTsKfQoKc3RhdGljIHZvaWQgdGVzdCAoc2l6ZV90IG4sIGludCBudW1faXRlciwgaW50IHVzZV93c3VtLCBpbnQgdXNlX3NpbXBsZV93c3VtKQp7CiAgICBkb3VibGUgKncgPSBtYWxsb2MobiAqIHNpemVvZih3WzBdKSk7CiAgICBmb3IgKHNpemVfdCBpID0gMDsgaSA8IG47IGkrKykgewogICAgICAgIHdbaV0gPSAoZG91YmxlKXJhbmQoKSAvIFJBTkRfTUFYOwogICAgfQogICAgCiAgICBkb3VibGUgZHVtbXkgPSAwOyAvLyBwcmV2ZW50IG9wdGltaXphdGlvbnMKICAgIAogICAgV1N1bSBzdGF0ZTsKICAgIFNpbXBsZVdTdW0gc2ltcGxlX3N0YXRlOwogICAgaWYgKHVzZV93c3VtKSB7CiAgICAgICAgV1N1bV9Jbml0KCZzdGF0ZSwgbiwgdyk7CiAgICAgICAgcHJpbnRmKCJXU3VtIG09JXp1IG49JXp1IG51bV9tPSV6dVxuIiwgc3RhdGUubSwgc3RhdGUubiwgc3RhdGUubnVtX20pOwogICAgfQogICAgaWYgKHVzZV9zaW1wbGVfd3N1bSkgewogICAgICAgIFNpbXBsZVdTdW1fSW5pdCgmc2ltcGxlX3N0YXRlLCBuLCB3KTsKICAgIH0KICAgIAogICAgZm9yIChpbnQgayA9IDA7IGsgPCBudW1faXRlcjsgaysrKSB7CiAgICAgICAgZG91YmxlIHggPSByYW5kKCk7CiAgICAgICAgCiAgICAgICAgZG91YmxlIHN1bSA9IDA7CiAgICAgICAgZG91YmxlIHNpbXBsZV9zdW0gPSAwOwogICAgICAgIGlmICh1c2Vfd3N1bSkgewogICAgICAgICAgICBzdW0gPSBXU3VtX0FkZE51bWJlcigmc3RhdGUsIHgpOwogICAgICAgIH0KICAgICAgICBpZiAodXNlX3NpbXBsZV93c3VtKSB7CiAgICAgICAgICAgIHNpbXBsZV9zdW0gPSBTaW1wbGVXU3VtX0FkZE51bWJlcigmc2ltcGxlX3N0YXRlLCB4KTsKICAgICAgICB9CiAgICAgICAgCiAgICAgICAgaWYgKHVzZV93c3VtICYmIHVzZV9zaW1wbGVfd3N1bSkgewogICAgICAgICAgICBpbnQgZmFpbGVkID0gKGZhYnMoc3VtIC0gc2ltcGxlX3N1bSkgPiAwLjAyKTsKICAgICAgICAgICAgcHJpbnRmKCIlZiAlZiAlZiAlc1xuIiwgeCwgc3VtLCBzaW1wbGVfc3VtLCAoZmFpbGVkID8gIkZBSUxFRCIgOiAiT0siKSk7CiAgICAgICAgICAgIGlmIChmYWlsZWQpIHsKICAgICAgICAgICAgICAgIGV4aXQoMSk7CiAgICAgICAgICAgIH0KICAgICAgICB9CiAgICAgICAgCiAgICAgICAgZHVtbXkgKz0gc3VtICsgc2ltcGxlX3N1bTsKICAgIH0KICAgIAogICAgcHJpbnRmKCJkdW1teSAlZlxuIiwgZHVtbXkpOwogICAgCiAgICBpZiAodXNlX3NpbXBsZV93c3VtKSB7CiAgICAgICAgU2ltcGxlV1N1bV9GcmVlKCZzaW1wbGVfc3RhdGUpOwogICAgfQogICAgaWYgKHVzZV93c3VtKSB7CiAgICAgICAgV1N1bV9GcmVlKCZzdGF0ZSk7CiAgICB9CiAgICBmcmVlKHcpOwp9CgovLyBBcm91bmQgaGVyZSBpdCBXU3VtIHN0YXJ0cyB0byBiZWF0IFNpbXBsZVdTdW0KLy8jZGVmaW5lIE5VTV9XRUlHSFRTIDE3MDAwCi8vI2RlZmluZSBOVU1fSVRFUiA4MDAwMAojZGVmaW5lIE5VTV9XRUlHSFRTIDUwCiNkZWZpbmUgTlVNX0lURVIgMjAwCgppbnQgbWFpbiAoKQp7CiAgICAvLyBSdW4gYm90aCwgY29tcGFyaW5nIHRoZSByZXN1bHRzLgogICAgdGVzdChOVU1fV0VJR0hUUywgTlVNX0lURVIsIDEsIDEpOwogICAgCiAgICAvLyBSdW4ganVzdCBXU3VtIG9yIGp1c3QgU2ltcGxlV1N1bSBmb3IgYmVuY2htYXJraW5nCiAgICAvL3Rlc3QoTlVNX1dFSUdIVFMsIE5VTV9JVEVSLCAxLCAwKTsgLy8gV1N1bQogICAgLy90ZXN0KE5VTV9XRUlHSFRTLCBOVU1fSVRFUiwgMCwgMSk7IC8vIFNpbXBsZVdTdW0KfQo=