/*-
 * Copyright (c) 2005-2009 Ariff Abdullah <ariff@FreeBSD.org>
 * 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 AUTHOR 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 AUTHOR 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.
 */

/*
 * feeder_rate: (Codename: Z Resampler), which means any effort to create
 *              future replacement for this resampler are simply absurd unless
 *              the world decide to add new alphabet after Z.
 *
 * FreeBSD bandlimited sinc interpolator, technically based on
 * "Digital Audio Resampling" by Julius O. Smith III
 *  - http://ccrma.stanford.edu/~jos/resample/
 *
 * The Good:
 * + all out fixed point integer operations, no soft-float or anything like
 *   that.
 * + classic polyphase converters with high quality coefficient's polynomial
 *   interpolators.
 * + fast, faster, or the fastest of its kind.
 * + compile time configurable.
 * + etc etc..
 *
 * The Bad:
 * - The z, z_, and Z_ . Due to mental block (or maybe just 0x7a69), I
 *   couldn't think of anything simpler than that (feeder_rate_xxx is just
 *   too long). Expect possible clashes with other zitizens (any?).
 */

#ifdef _KERNEL
#ifdef HAVE_KERNEL_OPTION_HEADERS
#include "opt_snd.h"
#endif
#include <dev/sound/pcm/sound.h>
#include <dev/sound/pcm/pcm.h>
#include <dev/sound/pcm/intpcm.h>
#include "feeder_if.h"

#define SND_USE_FXDIV
#include "snd_fxdiv_gen.h"

SND_DECLARE_FILE("$FreeBSD: head/sys/dev/sound/pcm/feeder_rate.c 195689 2009-07-14 18:53:34Z ariff $");
#endif

#include "feeder_rate_gen.h"

#if !defined(_KERNEL) && defined(SND_DIAGNOSTIC)
#undef Z_DIAGNOSTIC
#define Z_DIAGNOSTIC            1
#elif defined(_KERNEL)
#undef Z_DIAGNOSTIC
#undef Z_STRESS_TEST
#endif

#ifndef Z_QUALITY_DEFAULT
#define Z_QUALITY_DEFAULT       Z_QUALITY_LINEAR
#endif

#define Z_RESERVOIR             2048
#define Z_RESERVOIR_MAX         131072

#define Z_DOWNMAX               48           /* 384000 / 8000 */

#define Z_SINC_MAX              0x3fffff
#define Z_SINC_DOWNMAX          Z_DOWNMAX

#ifdef _KERNEL
#define Z_POLYPHASE_MAX         183040         /* 286 taps, 640 phases */
#else
#define Z_POLYPHASE_MAX         1464320                /* 286 taps, 5120 phases */
#endif

#define Z_RATE_DEFAULT          48000

#define Z_RATE_MIN              FEEDRATE_RATEMIN
#define Z_RATE_MAX              FEEDRATE_RATEMAX
#define Z_ROUNDHZ               FEEDRATE_ROUNDHZ
#define Z_ROUNDHZ_MIN           FEEDRATE_ROUNDHZ_MIN
#define Z_ROUNDHZ_MAX           FEEDRATE_ROUNDHZ_MAX

#define Z_RATE_SRC              FEEDRATE_SRC
#define Z_RATE_DST              FEEDRATE_DST
#define Z_RATE_QUALITY          FEEDRATE_QUALITY
#define Z_RATE_CHANNELS         FEEDRATE_CHANNELS

#define Z_PARANOID              1

#define Z_MULTIFORMAT           0

#define Z_FACTOR_MIN            1
#define Z_FACTOR_MAX            Z_MASK
#define Z_FACTOR_SAFE(v)        (!((v) < Z_FACTOR_MIN || (v) > Z_FACTOR_MAX))

struct z_info;

typedef void (*z_resampler_t)(struct z_info *, intpcm_t *, int32_t);

struct z_info {
        int32_t rsrc, rdst;    /* original source / destination rates */
        int32_t src, dst;      /* rounded source / destination rates */
        int32_t channels;      /* total channels */
        int32_t quality;       /* resampling quality */

        int32_t z_gx, z_gy;    /* interpolation / decimation ratio */
        int32_t z_alpha;       /* output sample time phase / drift */
        intpcm_t *z_delay;     /* FIR delay line / linear buffer */
        int32_t *z_coeff;      /* FIR coefficients */
        int32_t *z_dcoeff;     /* FIR coefficients differences */
        int32_t *z_pcoeff;     /* FIR polyphase coefficients */
        int32_t z_scale;       /* output scaling */
        int32_t z_dx;          /* input sample drift increment */
        int32_t z_dy;          /* output sample drift increment */
        int32_t z_alphadrift;  /* alpha drift rate */
        int32_t z_startdrift;  /* buffer start position drift rate */
        int32_t z_size;                /* half width of FIR taps */
        int32_t z_full;                /* full size of delay line */
        int32_t z_alloc;       /* largest allocated full size of delay line */
        int32_t z_start;       /* buffer processing start position */
        int32_t z_pos;         /* current position for the next feed */
#ifdef Z_DIAGNOSTIC
        uint32_t z_cycle;      /* output cycle, purely for statistical */
#endif
        int32_t z_maxfeed;     /* maximum feed to avoid 32bit overflow */

        z_resampler_t z_resample;
};

int feeder_rate_min = Z_RATE_MIN;
int feeder_rate_max = Z_RATE_MAX;
int feeder_rate_round = Z_ROUNDHZ;
int feeder_rate_quality = Z_QUALITY_DEFAULT;

static int feeder_rate_polyphase_max = Z_POLYPHASE_MAX;

#ifdef _KERNEL
static const char feeder_rate_presets[] = FEEDER_RATE_PRESETS;
SYSCTL_STRING(_hw_snd, OID_AUTO, feeder_rate_presets, CTLFLAG_RD,
    &feeder_rate_presets, 0, "compile-time rate presets");

TUNABLE_INT("hw.snd.feeder_rate_min", &feeder_rate_min);
TUNABLE_INT("hw.snd.feeder_rate_max", &feeder_rate_max);
TUNABLE_INT("hw.snd.feeder_rate_round", &feeder_rate_round);
TUNABLE_INT("hw.snd.feeder_rate_quality", &feeder_rate_quality);

TUNABLE_INT("hw.snd.feeder_rate_polyphase_max", &feeder_rate_polyphase_max);
SYSCTL_INT(_hw_snd, OID_AUTO, feeder_rate_polyphase_max, CTLFLAG_RW,
    &feeder_rate_polyphase_max, 0, "maximum allowable polyphase entries");

static int
sysctl_hw_snd_feeder_rate_min(SYSCTL_HANDLER_ARGS)
{
        int err, val;

        val = feeder_rate_min;
        err = sysctl_handle_int(oidp, &val, 0, req);

        if (err != 0 || req->newptr == NULL || val == feeder_rate_min)
                return (err);

        if (!(Z_FACTOR_SAFE(val) && val < feeder_rate_max))
                return (EINVAL);

        feeder_rate_min = val;

        return (0);
}
SYSCTL_PROC(_hw_snd, OID_AUTO, feeder_rate_min, CTLTYPE_INT | CTLFLAG_RW,
    0, sizeof(int), sysctl_hw_snd_feeder_rate_min, "I",
    "minimum allowable rate");

static int
sysctl_hw_snd_feeder_rate_max(SYSCTL_HANDLER_ARGS)
{
        int err, val;

        val = feeder_rate_max;
        err = sysctl_handle_int(oidp, &val, 0, req);

        if (err != 0 || req->newptr == NULL || val == feeder_rate_max)
                return (err);

        if (!(Z_FACTOR_SAFE(val) && val > feeder_rate_min))
                return (EINVAL);

        feeder_rate_max = val;

        return (0);
}
SYSCTL_PROC(_hw_snd, OID_AUTO, feeder_rate_max, CTLTYPE_INT | CTLFLAG_RW,
    0, sizeof(int), sysctl_hw_snd_feeder_rate_max, "I",
    "maximum allowable rate");

static int
sysctl_hw_snd_feeder_rate_round(SYSCTL_HANDLER_ARGS)
{
        int err, val;

        val = feeder_rate_round;
        err = sysctl_handle_int(oidp, &val, 0, req);

        if (err != 0 || req->newptr == NULL || val == feeder_rate_round)
                return (err);

        if (val < Z_ROUNDHZ_MIN || val > Z_ROUNDHZ_MAX)
                return (EINVAL);

        feeder_rate_round = val - (val % Z_ROUNDHZ);

        return (0);
}
SYSCTL_PROC(_hw_snd, OID_AUTO, feeder_rate_round, CTLTYPE_INT | CTLFLAG_RW,
    0, sizeof(int), sysctl_hw_snd_feeder_rate_round, "I",
    "sample rate converter rounding threshold");

static int
sysctl_hw_snd_feeder_rate_quality(SYSCTL_HANDLER_ARGS)
{
        struct snddev_info *d;
        struct pcm_channel *c;
        struct pcm_feeder *f;
        int i, err, val;

        val = feeder_rate_quality;
        err = sysctl_handle_int(oidp, &val, 0, req);

        if (err != 0 || req->newptr == NULL || val == feeder_rate_quality)
                return (err);

        if (val < Z_QUALITY_MIN || val > Z_QUALITY_MAX)
                return (EINVAL);

        feeder_rate_quality = val;

        /*
         * Traverse all available channels on each device and try to
         * set resampler quality if and only if it is exist as
         * part of feeder chains and the channel is idle.
         */
        for (i = 0; pcm_devclass != NULL &&
            i < devclass_get_maxunit(pcm_devclass); i++) {
                d = devclass_get_softc(pcm_devclass, i);
                if (!PCM_REGISTERED(d))
                        continue;
                PCM_LOCK(d);
                PCM_WAIT(d);
                PCM_ACQUIRE(d);
                CHN_FOREACH(c, d, channels.pcm) {
                        CHN_LOCK(c);
                        f = chn_findfeeder(c, FEEDER_RATE);
                        if (f == NULL || f->data == NULL || CHN_STARTED(c)) {
                                CHN_UNLOCK(c);
                                continue;
                        }
                        (void)FEEDER_SET(f, FEEDRATE_QUALITY, val);
                        CHN_UNLOCK(c);
                }
                PCM_RELEASE(d);
                PCM_UNLOCK(d);
        }

        return (0);
}
SYSCTL_PROC(_hw_snd, OID_AUTO, feeder_rate_quality, CTLTYPE_INT | CTLFLAG_RW,
    0, sizeof(int), sysctl_hw_snd_feeder_rate_quality, "I",
    "sample rate converter quality ("__XSTRING(Z_QUALITY_MIN)"=low .. "
    __XSTRING(Z_QUALITY_MAX)"=high)");
#endif  /* _KERNEL */


/*
 * Resampler type.
 */
#define Z_IS_ZOH(i)             ((i)->quality == Z_QUALITY_ZOH)
#define Z_IS_LINEAR(i)          ((i)->quality == Z_QUALITY_LINEAR)
#define Z_IS_SINC(i)            ((i)->quality > Z_QUALITY_LINEAR)

/*
 * Linear interpolation with downsampling ratio of 1 is practically zoh.
 */
#define Z_IS_LINEAR_ZOH(i)      (Z_IS_LINEAR(i) && (i)->z_gy == 1)

/*
 * Macroses for accurate sample time drift calculations.
 *
 * gy2gx : given the amount of output, return the _exact_ required amount of
 *         input.
 * gx2gy : given the amount of input, return the _maximum_ amount of output
 *         that will be generated.
 * drift : given the amount of input and output, return the elapsed
 *         sample-time.
 */
#define _Z_GCAST(x)             ((uint64_t)(x))

#if defined(__GNUCLIKE_ASM) && defined(__i386__)
/*
 * This is where i386 being beaten to a pulp. Fortunately this function is
 * rarely being called and if it is, it will decide the best (hopefully)
 * fastest way to do the division. If we can ensure that everything is dword
 * aligned, letting the compiler to call udivdi3 to do the division can be
 * faster compared to this.
 *
 * amd64 is the clear winner here, no question about it.
 */
static __inline uint32_t
Z_DIV(uint64_t v, uint32_t d)
{
        uint32_t hi, lo, quo, rem;

        hi = v >> 32;
        lo = v & 0xffffffff;

        /*
         * As much as we can, try to avoid long division like a plague.
         */
        if (hi == 0)
                quo = lo / d;
        else
                __asm("divl %2"
                    : "=a" (quo), "=d" (rem)
                    : "r" (d), "0" (lo), "1" (hi));

        return (quo);
}
#else
#define Z_DIV(x, y)             ((x) / (y))
#endif

#define _Z_GY2GX(i, a, v)                                               \
        Z_DIV(((_Z_GCAST((i)->z_gx) * (v)) + ((i)->z_gy - (a) - 1)),   \
        (i)->z_gy)

#define _Z_GX2GY(i, a, v)                                               \
        Z_DIV(((_Z_GCAST((i)->z_gy) * (v)) + (a)), (i)->z_gx)

#define _Z_DRIFT(i, x, y)                                               \
        ((_Z_GCAST((i)->z_gy) * (x)) - (_Z_GCAST((i)->z_gx) * (y)))

#define z_gy2gx(i, v)           _Z_GY2GX(i, (i)->z_alpha, v)
#define z_gx2gy(i, v)           _Z_GX2GY(i, (i)->z_alpha, v)
#define z_drift(i, x, y)        _Z_DRIFT(i, x, y)

#define Z_SRC_DOWN(i)           ((i)->z_gx > (i)->z_gy)
#define Z_SRC_UP(i)             (!Z_SRC_DOWN(i))

#define Z_FACTOR_DOWN(i)        (Z_SRC_DOWN(i) ? _Z_GY2GX(i, 0, 1) : 0)
#define Z_FACTOR_UP(i)          (Z_SRC_UP(i) ? _Z_GX2GY(i, 0, 1) : 0)

/*
 * Macroses for SINC coefficients table manipulations.. whatever.
 */
#define Z_SINC_COEFF_IDX(i)     ((i)->quality - Z_QUALITY_LINEAR - 1)

#define Z_SINC_LEN(i)                                                   \
        ((int32_t)((((uint64_t)z_coeff_tab[Z_SINC_COEFF_IDX(i)].len << \
            Z_SHIFT) - 1) / (i)->z_dy))

#define Z_SINC_BASE_LEN(i)                                              \
        ((z_coeff_tab[Z_SINC_COEFF_IDX(i)].len - 1) >> (Z_DRIFT_SHIFT - 1))

/*
 * Macroses for linear delay buffer operations.
 */
#define z_fetched(i)            ((i)->z_pos - (i)->z_start - 1)
#define z_free(i)               ((i)->z_full - (i)->z_pos)

/*
 * For !sinc, z_size represents full history.
 */
#define z_history(i)            (Z_IS_SINC(i) ? ((i)->z_size << 1) :      \
                                (i)->z_size)

#define z_start_drift(i)        _Z_GY2GX(i, 0, 1)

/*
 * Macroses for Bla Bla .. :)
 */
#define z_copy(src, dst, sz)    (void)memcpy(dst, src, sz)
#define z_feed(...)             FEEDER_FEED(__VA_ARGS__)

static __inline uint32_t
z_min(uint32_t x, uint32_t y)
{

        return ((x < y) ? x : y);
}

static int32_t
z_gcd(int32_t x, int32_t y)
{
        int32_t w;

        while (y != 0) {
                w = x % y;
                x = y;
                y = w;
        }

        return (x);
}

#ifndef Z_STRESS_TEST
static int32_t
z_roundpow2(int32_t v)
{
        int32_t i;

        i = 1;

        /*
         * Let it overflow at will..
         */
        while (i > 0 && i < v)
                i <<= 1;

        return (i);
}
#endif

/*
 * Zero Order Hold, the worst of the worst, an insult against quality,
 * but super fast.
 */
static void
z_feed_zoh(struct z_info *info, intpcm_t *b, int32_t count)
{
        int32_t ch, reqout, z_alpha, z_start;
        intpcm_t *src, *dst;

        ch = info->channels - 1;
        do {
                z_alpha = info->z_alpha;
                z_start = info->z_start;
                reqout = count;
                dst = b + ch;
                do {
                        z_alpha += info->z_alphadrift;
                        z_start += info->z_startdrift;
                        if (z_alpha >= info->z_gy) {
                                z_alpha -= info->z_gy;
                                z_start -= 1;
                        }
                        src = info->z_delay + (z_start * info->channels);
                        *dst = *src;
                        dst += info->channels;
                } while (--reqout);
        } while (ch-- != 0);

        info->z_alpha = z_alpha;
        info->z_start = z_start;
}

/*
 * Linear Interpolation. This at least sounds better (perceptually) and fast,
 * but without any proper filtering which means aliasing still exist and
 * could become worst with a right sample. Interpolation centered within
 * Z_LINEAR_ONE between the present and previous sample and everything is
 * done with simple 64bit or 32bit scaling arithmetic.
 *
 * 64bit linear interpolation - Full dynamic range, slower on 32bit arch
 *                              (though still quite fast).
 *
 * 32bit linear interpolation - Fraction and interpolation reduced to make
 *                              room for 32bit arithmetic at the cost
 *                              of conversion accuracy but lightning fast.
 *
 * Whether it is 64bit or 32bit, these are only interesting from analysis
 * point of view. Draging it further for the sake of quality is fighting
 * a losing cause.
 */
static void
z_feed_linear(struct z_info *info, intpcm_t *b, int32_t count)
{
        int32_t z, ch, reqout, z_alpha, z_start;
        intpcm_t *src, *dst;

        ch = info->channels - 1;

        do {
                z_alpha = info->z_alpha;
                z_start = info->z_start;
                reqout = count;
                dst = b + ch;
                do {
                        z_alpha += info->z_alphadrift;
                        z_start += info->z_startdrift;
                        if (z_alpha >= info->z_gy) {
                                z_alpha -= info->z_gy;
                                z_start -= 1;
                        }
                        z = Z_LINEAR_FRACTION_32(z_alpha * info->z_dx);
                        src = info->z_delay + (z_start * info->channels);
                        *dst = Z_LINEAR_INTERPOLATE_32(z, *src,
                            *(src - info->channels));
                        dst += info->channels;
                } while (--reqout);
        } while (ch-- != 0);

        info->z_alpha = z_alpha;
        info->z_start = z_start;
}

/*
 * Userland clipping diagnostic check, not enabled in kernel compilation.
 * While doing sinc interpolation, unrealistic samples like full scale sine
 * wav will clip, but for other things this will not make any noise at all.
 * Everybody should learn how to normalized perceived loudness of their own
 * music/sounds/samples (hint: ReplayGain).
 */
#ifdef Z_DIAGNOSTIC
#define Z_CLIP_CHECK(v) do {                                            \
        if ((v) > PCM_S32_MAX) {                                       \
                fprintf(stderr, "Overflow: v=%jd, max=%jd\n",         \
                    (intmax_t)(v), (intmax_t)PCM_S32_MAX);            \
        } else if ((v) < PCM_S32_MIN) {                                        \
                fprintf(stderr, "Underflow: v=%jd, min=%jd\n",                \
                    (intmax_t)(v), (intmax_t)PCM_S32_MIN);            \
        }                                                               \
} while (0)
#else
#define Z_CLIP_CHECK(...)
#endif

#define Z_CLAMP(v)                                                      \
        (((v) > PCM_S32_MAX) ? PCM_S32_MAX :                            \
        (((v) < PCM_S32_MIN) ? PCM_S32_MIN : (v)))

/*
 * Sine Cardinal (SINC) Interpolation. Scaling is done in 64 bit, so
 * there's no point to hold the plate any longer. All samples will be
 * shifted to a full 32 bit, scaled and restored during write for
 * maximum dynamic range (only for downsampling).
 */
#define _Z_SINC_ACCUMULATE(adv, acc)                                    \
        c = z >> Z_SHIFT;                                               \
        z &= Z_MASK;                                                    \
        z_coeff += c;                                                   \
        z_dcoeff += c;                                                  \
        coeff = Z_COEFF_INTERPOLATE(z, *z_coeff, *z_dcoeff);            \
        acc += Z_NORM_32((intpcm64_t)(*p) * coeff);                    \
        z += info->z_dy;                                                \
        p adv##= info->channels

#define _Z_SINC_FAST_ACCUMULATE(adv, acc)                               \
        z_coeff += z >> Z_SHIFT;                                        \
        z &= Z_MASK;                                                    \
        acc += Z_NORM_32((intpcm64_t)(*p) * *z_coeff);                 \
        z += info->z_dy;                                                \
        p adv##= info->channels

/* 
 * XXX GCC4 optimization is such a !@#$%, need manual unrolling.
 */
#if defined(__GNUC__) && __GNUC__ >= 4
#define Z_SINC_ACCUMULATE(...)  do {                                     \
        _Z_SINC_ACCUMULATE(__VA_ARGS__);                                \
        _Z_SINC_ACCUMULATE(__VA_ARGS__);                                \
} while (0)
#define Z_SINC_FAST_ACCUMULATE(...)     do {                                \
        _Z_SINC_FAST_ACCUMULATE(__VA_ARGS__);                           \
        _Z_SINC_FAST_ACCUMULATE(__VA_ARGS__);                           \
} while (0)
#define Z_SINC_ACCUMULATE_DECR          2
#else
#define Z_SINC_ACCUMULATE(...)  do {                                     \
        _Z_SINC_ACCUMULATE(__VA_ARGS__);                                \
} while (0)
#define Z_SINC_FAST_ACCUMULATE(...)     do {                                \
        _Z_SINC_FAST_ACCUMULATE(__VA_ARGS__);                           \
} while (0)
#define Z_SINC_ACCUMULATE_DECR          1
#endif

static void
z_feed_sinc(struct z_info *info, intpcm_t *b, int32_t count)
{
        intpcm64_t v, v0;
        intpcm_t *p, *dst;
        int32_t coeff, z, reqout, z_alpha, z_start, *z_coeff, *z_dcoeff;
        uint32_t c, ch, i;


        ch = info->channels - 1;

        do {
                z_alpha = info->z_alpha;
                z_start = info->z_start;
                reqout = count;
                dst = b + ch;
                do {
                        z_alpha += info->z_alphadrift;
                        z_start += info->z_startdrift;
                        if (z_alpha >= info->z_gy) {
                                z_alpha -= info->z_gy;
                                z_start -= 1;
                        }
                        v = v0 = 0;
                        z = z_alpha * info->z_dx;
                        p = info->z_delay + ((z_start - info->z_size + 1) *
                            info->channels) + ch;
                        z_coeff = info->z_coeff;
                        z_dcoeff = info->z_dcoeff;
                        for (i = info->z_size; i != 0;
                            i -= Z_SINC_ACCUMULATE_DECR)
                                Z_SINC_ACCUMULATE(+, v);
                        z = info->z_dy - (z_alpha * info->z_dx);
                        p = info->z_delay + ((z_start - info->z_size) *
                            info->channels) + ch;
                        z_coeff = info->z_coeff;
                        z_dcoeff = info->z_dcoeff;
                        for (i = info->z_size; i != 0;
                            i -= Z_SINC_ACCUMULATE_DECR)
                                Z_SINC_ACCUMULATE(-, v0);
                        v += v0;
                        v >>= Z_COEFF_SHIFT - Z_GUARD_BIT_32;
                        if (info->z_scale != Z_ONE)
                                v = (v * Z_SCALE_CAST(info->z_scale)) >>
                                    Z_SHIFT;
                        Z_CLIP_CHECK(v);
                        *dst = Z_CLAMP(v);
                        dst += info->channels;
                } while (--reqout);
        } while (ch-- != 0);

        info->z_alpha = z_alpha;
        info->z_start = z_start;
}

static void
z_feed_sinc_fast(struct z_info *info, intpcm_t *b, int32_t count)
{
        intpcm64_t v, v0;
        intpcm_t *p, *dst;
        int32_t z, reqout, z_alpha, z_start, *z_coeff;
        uint32_t ch, i;


        ch = info->channels - 1;

        do {
                z_alpha = info->z_alpha;
                z_start = info->z_start;
                reqout = count;
                dst = b + ch;
                do {
                        z_alpha += info->z_alphadrift;
                        z_start += info->z_startdrift;
                        if (z_alpha >= info->z_gy) {
                                z_alpha -= info->z_gy;
                                z_start -= 1;
                        }
                        v = v0 = 0;
                        z = z_alpha * info->z_dx;
                        p = info->z_delay + ((z_start - info->z_size + 1) *
                            info->channels) + ch;
                        z_coeff = info->z_coeff;
                        for (i = info->z_size; i != 0;
                            i -= Z_SINC_ACCUMULATE_DECR)
                                Z_SINC_FAST_ACCUMULATE(+, v);
                        z = info->z_dy - (z_alpha * info->z_dx);
                        p = info->z_delay + ((z_start - info->z_size) *
                            info->channels) + ch;
                        z_coeff = info->z_coeff;
                        for (i = info->z_size; i != 0;
                            i -= Z_SINC_ACCUMULATE_DECR)
                                Z_SINC_FAST_ACCUMULATE(-, v0);
                        v += v0;
                        v >>= Z_COEFF_SHIFT - Z_GUARD_BIT_32;
                        if (info->z_scale != Z_ONE)
                                v = (v * Z_SCALE_CAST(info->z_scale)) >>
                                    Z_SHIFT;
                        Z_CLIP_CHECK(v);
                        *dst = Z_CLAMP(v);
                        dst += info->channels;
                } while (--reqout);
        } while (ch-- != 0);

        info->z_alpha = z_alpha;
        info->z_start = z_start;
}

static void
z_feed_sinc_polyphase(struct z_info *info, intpcm_t *b, int32_t count)
{
        #if 0
        intpcm64_t v, v0;
        intpcm_t *p;
        int32_t *z_pcoeff;
        uint32_t ch, i, start;

        if (1) {
                z_feed_zoh(info, dst, count);
                return;
        }

        ch = info->channels;
        dst += ch;
        start = (info->z_start - (info->z_size << 1) + 1) * ch;

        do {
                v = v0 = 0;
                dst--;
                ch--;
                p = info->z_delay + start + ch;
                z_pcoeff = info->z_pcoeff +
                    ((info->z_alpha * info->z_size) << 1);
                for (i = info->z_size; i != 0; i--) {
                        v += Z_NORM_32((intpcm64_t)(*p) * *z_pcoeff);
                        z_pcoeff++;
                        p += info->channels;
                        v0 += Z_NORM_32((intpcm64_t)(*p) * *z_pcoeff);
                        z_pcoeff++;
                        p += info->channels;
                }
                v += v0;
                v >>= Z_COEFF_SHIFT - Z_GUARD_BIT_32;
                if (info->z_scale != Z_ONE)
                        v = (v * Z_SCALE_CAST(info->z_scale)) >> Z_SHIFT;
                Z_CLIP_CHECK(v);
                *dst = Z_CLAMP(v);
        } while (ch != 0);
        #endif
        intpcm64_t v, v0;
        intpcm_t *p, *dst;
        int32_t reqout, z_alpha, z_start, *z_pcoeff;
        uint32_t ch, i;

        ch = info->channels - 1;

        do {
                z_alpha = info->z_alpha;
                z_start = info->z_start;
                reqout = count;
                dst = b + ch;
                do {
                        z_alpha += info->z_alphadrift;
                        z_start += info->z_startdrift;
                        if (z_alpha >= info->z_gy) {
                                z_alpha -= info->z_gy;
                                z_start -= 1;
                        }
                        v = v0 = 0;
                        p = info->z_delay + ((z_start -
                            (info->z_size << 1) + 1) * info->channels) + ch;
                        z_pcoeff = info->z_pcoeff +
                            ((z_alpha * info->z_size) << 1);
                        for (i = info->z_size; i != 0; i--) {
                                v += Z_NORM_32((intpcm64_t)(*p) * *z_pcoeff);
                                z_pcoeff++;
                                p += info->channels;
                                v0 += Z_NORM_32((intpcm64_t)(*p) * *z_pcoeff);
                                z_pcoeff++;
                                p += info->channels;
                        }
                        v += v0;
                        v >>= Z_COEFF_SHIFT - Z_GUARD_BIT_32;
                        if (info->z_scale != Z_ONE)
                                v = (v * Z_SCALE_CAST(info->z_scale)) >>
                                    Z_SHIFT;
                        Z_CLIP_CHECK(v);
                        *dst = Z_CLAMP(v);
                        dst += info->channels;
                } while (--reqout);
        } while (ch-- != 0);

        info->z_alpha = z_alpha;
        info->z_start = z_start;
}

enum {
        Z_RESAMPLER_ZOH,
        Z_RESAMPLER_LINEAR,
        Z_RESAMPLER_SINC
};

#define Z_RESAMPLER_IDX(i)                                              \
        (Z_IS_SINC(i) ? Z_RESAMPLER_SINC : (i)->quality)

static void
z_resampler_reset(struct z_info *info)
{

        info->src = info->rsrc - (info->rsrc % ((feeder_rate_round > 0 &&
            info->rsrc > feeder_rate_round) ? feeder_rate_round : 1));
        info->dst = info->rdst - (info->rdst % ((feeder_rate_round > 0 &&
            info->rdst > feeder_rate_round) ? feeder_rate_round : 1));
        info->z_gx = 1;
        info->z_gy = 1;
        info->z_alpha = 0;
        info->z_resample = NULL;
        info->z_size = 1;
        info->z_coeff = NULL;
        info->z_dcoeff = NULL;
        if (info->z_pcoeff != NULL) {
                free(info->z_pcoeff, M_DEVBUF);
                info->z_pcoeff = NULL;
        }
        info->z_scale = Z_ONE;
        info->z_dx = Z_FULL_ONE;
        info->z_dy = Z_FULL_ONE;
#ifdef Z_DIAGNOSTIC
        info->z_cycle = 0;
#endif
        if (info->quality < Z_QUALITY_MIN)
                info->quality = Z_QUALITY_MIN;
        else if (info->quality > Z_QUALITY_MAX)
                info->quality = Z_QUALITY_MAX;
}

#ifdef Z_PARANOID
static int32_t
z_resampler_sinc_len(struct z_info *info)
{
        int32_t c, z, len, lmax;

        if (!Z_IS_SINC(info))
                return (1);

        /*
         * A rather careful (or useless) way to calculate filter length.
         * Z_SINC_LEN() itself is accurate enough to do its job. Extra
         * sanity checking is not going to hurt though..
         */
        c = 0;
        z = info->z_dy;
        len = 0;
        lmax = z_coeff_tab[Z_SINC_COEFF_IDX(info)].len;

        do {
                c += z >> Z_SHIFT;
                z &= Z_MASK;
                z += info->z_dy;
        } while (c < lmax && ++len > 0);

        if (len != Z_SINC_LEN(info)) {
#ifdef _KERNEL
                printf("%s(): sinc l=%d != Z_SINC_LEN=%d\n",
                    __func__, len, Z_SINC_LEN(info));
#else
                fprintf(stderr, "%s(): sinc l=%d != Z_SINC_LEN=%d\n",
                    __func__, len, Z_SINC_LEN(info));
#endif
        }

        return (len);
}
#else
#define z_resampler_sinc_len(i)         (Z_IS_SINC(i) ? Z_SINC_LEN(i) : 1)
#endif

#define Z_POLYPHASE_COEFF_SHIFT         0

/*
 * Pick suitable polynomial interpolators based on filter oversampled ratio
 * (2 ^ Z_DRIFT_SHIFT).
 */
#if !(defined(Z_COEFF_INTERP_ZOH) || defined(Z_COEFF_INTERP_LINEAR) ||          \
    defined(Z_COEFF_INTERP_QUADRATIC) || defined(Z_COEFF_INTERP_HERMITE) ||     \
    defined(Z_COEFF_INTERP_BSPLINE) || defined(Z_COEFF_INTERP_BSPLINE_3) ||     \
    defined(Z_COEFF_INTERP_BSPLINE_5) || defined(Z_COEFF_INTERP_OPT32X) ||      \
    defined(Z_COEFF_INTERP_OPT16X) || defined(Z_COEFF_INTERP_OPT8X) ||          \
    defined(Z_COEFF_INTERP_OPT4X) || defined(Z_COEFF_INTERP_OPT2X))
#if Z_DRIFT_SHIFT >= 7
#define Z_COEFF_INTERP_BSPLINE_3        1
#elif Z_DRIFT_SHIFT == 6
#define Z_COEFF_INTERP_BSPLINE_5        1
#elif Z_DRIFT_SHIFT == 5
#define Z_COEFF_INTERP_OPT32X           1
#elif Z_DRIFT_SHIFT == 4
#define Z_COEFF_INTERP_OPT16X           1
#elif Z_DRIFT_SHIFT == 3
#define Z_COEFF_INTERP_OPT8X            1
#elif Z_DRIFT_SHIFT == 2
#define Z_COEFF_INTERP_OPT4X            1
#elif Z_DRIFT_SHIFT == 1
#define Z_COEFF_INTERP_OPT2X            1
#else
#error "Z_DRIFT_SHIFT screwed!"
#endif
#elif defined(Z_COEFF_INTERP_BSPLINE) && !(defined(Z_COEFF_INTERP_BSPLINE_3) || \
    defined(Z_COEFF_INTERP_BSPLINE_5))
#if Z_DRIFT_SHIFT >= 7
#define Z_COEFF_INTERP_BSPLINE_3        1
#else
#define Z_COEFF_INTERP_BSPLINE_5        1
#endif
#endif

/*
 * In classic polyphase mode, the actual coefficients for each phases need to
 * be calculated based on default prototype filters. For highly oversampled
 * filter, linear or quadradatic interpolator should be enough. Anything less
 * than that require 'special' interpolators to reduce interpolation errors.
 *
 * "Polynomial Interpolators for High-Quality Resampling of Oversampled Audio"
 *    by Olli Niemitalo
 *    - http://www.student.oulu.fi/~oniemita/dsp/deip.pdf
 *
 */
static int32_t
z_coeff_interpolate(int32_t z, int32_t *z_coeff)
{
        int32_t coeff;
#if defined(Z_COEFF_INTERP_ZOH)

        /* 1-point, 0th-order (Zero Order Hold) */
        z = z;
        coeff = z_coeff[0];
#elif defined(Z_COEFF_INTERP_LINEAR)
        int32_t zl0, zl1;

        /* 2-point, 1st-order Linear */
        zl0 = z_coeff[0];
        zl1 = z_coeff[1] - z_coeff[0];

        coeff = Z_RSHIFT((int64_t)zl1 * z, Z_SHIFT) + zl0;
#elif defined(Z_COEFF_INTERP_QUADRATIC)
        int32_t zq0, zq1, zq2;

        /* 3-point, 2nd-order Quadratic */
        zq0 = z_coeff[0];
        zq1 = z_coeff[1] - z_coeff[-1];
        zq2 = z_coeff[1] + z_coeff[-1] - (z_coeff[0] << 1);

        coeff = Z_RSHIFT((Z_RSHIFT((int64_t)zq2 * z, Z_SHIFT) +
            zq1) * z, Z_SHIFT + 1) + zq0;
#elif defined(Z_COEFF_INTERP_HERMITE)
        int32_t zh0, zh1, zh2, zh3;

        /* 4-point, 3rd-order Hermite */
        zh0 = z_coeff[0];
        zh1 = z_coeff[1] - z_coeff[-1];
        zh2 = ((int64_t)z_coeff[1] << 2) - ((int64_t)z_coeff[0] * 5) +
            (z_coeff[-1] << 1) - z_coeff[2];
        zh3 = ((int64_t)(z_coeff[0] - z_coeff[1]) * 3) + z_coeff[2] -
            z_coeff[-1];

        coeff = Z_RSHIFT((Z_RSHIFT((Z_RSHIFT((int64_t)zh3 * z, Z_SHIFT) +
            zh2) * z, Z_SHIFT) + zh1) * z, Z_SHIFT + 1) + zh0;
#elif defined(Z_COEFF_INTERP_BSPLINE_3)
        int32_t zb0, zb1, zb2, zb3;

        /* 4-point, 3rd-order B-Spline */
        zb0 = Z_RSHIFT(0x15555555LL * (((int64_t)z_coeff[0] << 2) +
            z_coeff[-1] + z_coeff[1]), 30);
        zb1 = z_coeff[1] - z_coeff[-1];
        zb2 = z_coeff[-1] + z_coeff[1] - (z_coeff[0] << 1);
        zb3 = Z_RSHIFT(0x15555555LL *
            (((int64_t)(z_coeff[0] - z_coeff[1]) * 3) + z_coeff[2] -
            z_coeff[-1]), 30);

        coeff = Z_RSHIFT(Z_RSHIFT((Z_RSHIFT((Z_RSHIFT(
            (int64_t)zb3 * z, Z_SHIFT) + zb2) * z, Z_SHIFT) +
            zb1) * z, Z_SHIFT) + zb0, 1);
#elif defined(Z_COEFF_INTERP_BSPLINE_5)
        int32_t zb0, zb1, zb2, zb3, zb4, zb5;

        /* 6-point, 5th-order B-Spline */
        zb0 = Z_RSHIFT((0x01111111LL * (z_coeff[-2] + z_coeff[2])) +
            (0x1bbbbbbcLL * (z_coeff[-1] + z_coeff[1])) +
            (0x46666666LL * z_coeff[0]), 30);
        zb1 = Z_RSHIFT((0x05555555LL * (z_coeff[2] - z_coeff[-2])) +
            (0x35555555LL * (z_coeff[1] - z_coeff[-1])), 30);
        zb2 = Z_RSHIFT((0x0aaaaaabLL * (z_coeff[-2] + z_coeff[2])) +
            (0x15555555LL * (z_coeff[-1] + z_coeff[1])) -
            (0x40000000LL * z_coeff[0]), 30);
        zb3 = Z_RSHIFT((0x0aaaaaabLL * (z_coeff[2] - z_coeff[-2])) -
            (0x15555555LL * (z_coeff[1] - z_coeff[-1])), 30);
        zb4 = Z_RSHIFT((0x05555555LL * (z_coeff[-2] + z_coeff[2])) -
            (0x15555555LL * (z_coeff[-1] + z_coeff[1])) +
            (0x20000000LL * z_coeff[0]), 30);
        zb5 = Z_RSHIFT((0x01111111LL * (z_coeff[3] - z_coeff[-2])) +
            (0x05555555LL * (z_coeff[-1] - z_coeff[2])) +
            (0x0aaaaaabLL * (z_coeff[1] - z_coeff[0])), 30);

        coeff = Z_RSHIFT(Z_RSHIFT((Z_RSHIFT((Z_RSHIFT((Z_RSHIFT((Z_RSHIFT(
            (int64_t)zb5 * z, Z_SHIFT) + zb4) * z, Z_SHIFT) +
            zb3) * z, Z_SHIFT) + zb2) * z, Z_SHIFT) + zb1) * z, Z_SHIFT) +
            zb0, 1);
#elif defined(Z_COEFF_INTERP_OPT32X)
        int32_t zoz, zoe1, zoe2, zoe3, zoo1, zoo2, zoo3;
        int32_t zoc0, zoc1, zoc2, zoc3, zoc4, zoc5;

        /* 6-point, 5th-order Optimal 32x */
        zoz = z - (Z_ONE >> 1);
        zoe1 = z_coeff[1] + z_coeff[0];
        zoe2 = z_coeff[2] + z_coeff[-1];
        zoe3 = z_coeff[3] + z_coeff[-2];
        zoo1 = z_coeff[1] - z_coeff[0];
        zoo2 = z_coeff[2] - z_coeff[-1];
        zoo3 = z_coeff[3] - z_coeff[-2];

        zoc0 = Z_RSHIFT((0x36a357d2LL * zoe1) + (0x0943c9cfLL * zoe2) +
            (0x0018de5fLL * zoe3), 30);
        zoc1 = Z_RSHIFT((0x2ddd5aecLL * zoo1) + (0x1a2d984bLL * zoo2) +
            (0x00b85f3eLL * zoo3), 30);
        zoc2 = Z_RSHIFT((-0x1bc6f4ecLL * zoe1) + (0x19aa6f62LL * zoe2) +
            (0x021c858aLL * zoe3), 30);
        zoc3 = Z_RSHIFT((-0x2024ef40LL * zoo1) + (0x0567cd19LL * zoo2) +
            (0x032f8198LL * zoo3), 30);
        zoc4 = Z_RSHIFT((0x05556cd2LL * zoe1) + (-0x0800233eLL * zoe2) +
            (0x02aab66bLL * zoe3), 30);
        zoc5 = Z_RSHIFT((0x0ab00fedLL * zoo1) + (-0x05580913LL * zoo2) +
            (0x01119bdcLL * zoo3), 30);

        coeff = Z_RSHIFT(Z_RSHIFT((Z_RSHIFT((Z_RSHIFT((Z_RSHIFT((Z_RSHIFT(
            (int64_t)zoc5 * zoz, Z_SHIFT) +
            zoc4) * zoz, Z_SHIFT) + zoc3) * zoz, Z_SHIFT) +
            zoc2) * zoz, Z_SHIFT) + zoc1) * zoz, Z_SHIFT) + zoc0, 1);
#elif defined(Z_COEFF_INTERP_OPT16X)
        int32_t zoz, zoe1, zoe2, zoe3, zoo1, zoo2, zoo3;
        int32_t zoc0, zoc1, zoc2, zoc3, zoc4, zoc5;

        /* 6-point, 5th-order Optimal 16x */
        zoz = z - (Z_ONE >> 1);
        zoe1 = z_coeff[1] + z_coeff[0];
        zoe2 = z_coeff[2] + z_coeff[-1];
        zoe3 = z_coeff[3] + z_coeff[-2];
        zoo1 = z_coeff[1] - z_coeff[0];
        zoo2 = z_coeff[2] - z_coeff[-1];
        zoo3 = z_coeff[3] - z_coeff[-2];

        zoc0 = Z_RSHIFT((0x35844c1aLL * zoe1) + (0x0a4d9b93LL * zoe2) +
            (0x002e1852LL * zoe3), 30);
        zoc1 = Z_RSHIFT((0x29f14933LL * zoo1) + (0x1ada2213LL * zoo2) +
            (0x0119a9b7LL * zoo3), 30);
        zoc2 = Z_RSHIFT((-0x1a7d2947LL * zoe1) + (0x17bbbda7LL * zoe2) +
            (0x02c16ba0LL * zoe3), 30);
        zoc3 = Z_RSHIFT((-0x1bc21989LL * zoo1) + (0x033654faLL * zoo2) +
            (0x039fd228LL * zoo3), 30);
        zoc4 = Z_RSHIFT((0x055425aeLL * zoe1) + (-0x07fe3767LL * zoe2) +
            (0x02aa11b9LL * zoe3), 30);
        zoc5 = Z_RSHIFT((0x0a3a53caLL * zoo1) + (-0x051cec8eLL * zoo2) +
            (0x0105b035LL * zoo3), 30);

        coeff = Z_RSHIFT(Z_RSHIFT((Z_RSHIFT((Z_RSHIFT((Z_RSHIFT((Z_RSHIFT(
            (int64_t)zoc5 * zoz, Z_SHIFT) +
            zoc4) * zoz, Z_SHIFT) + zoc3) * zoz, Z_SHIFT) +
            zoc2) * zoz, Z_SHIFT) + zoc1) * zoz, Z_SHIFT) + zoc0, 1);
#elif defined(Z_COEFF_INTERP_OPT8X)
        int32_t zoz, zoe1, zoe2, zoe3, zoo1, zoo2, zoo3;
        int32_t zoc0, zoc1, zoc2, zoc3, zoc4, zoc5;

        /* 6-point, 5th-order Optimal 8x */
        zoz = z - (Z_ONE >> 1);
        zoe1 = z_coeff[1] + z_coeff[0];
        zoe2 = z_coeff[2] + z_coeff[-1];
        zoe3 = z_coeff[3] + z_coeff[-2];
        zoo1 = z_coeff[1] - z_coeff[0];
        zoo2 = z_coeff[2] - z_coeff[-1];
        zoo3 = z_coeff[3] - z_coeff[-2];

        zoc0 = Z_RSHIFT((0x355368f9LL * zoe1) + (0x0a7b3288LL * zoe2) +
            (0x0031647fLL * zoe3), 30);
        zoc1 = Z_RSHIFT((0x294209a1LL * zoo1) + (0x1afa4a08LL * zoo2) +
            (0x01296b33LL * zoo3), 30);
        zoc2 = Z_RSHIFT((-0x1a44a615LL * zoe1) + (0x1766f457LL * zoe2) +
            (0x02ddb1bfLL * zoe3), 30);
        zoc3 = Z_RSHIFT((-0x1ae89637LL * zoo1) + (0x02c92b21LL * zoo2) +
            (0x03b5d273LL * zoo3), 30);
        zoc4 = Z_RSHIFT((0x054fdc36LL * zoe1) + (-0x07f7b647LL * zoe2) +
            (0x02a7da0eLL * zoe3), 30);
        zoc5 = Z_RSHIFT((0x099f36d8LL * zoo1) + (-0x04cd66efLL * zoo2) +
            (0x00f4f84dLL * zoo3), 30);

        coeff = Z_RSHIFT(Z_RSHIFT((Z_RSHIFT((Z_RSHIFT((Z_RSHIFT((Z_RSHIFT(
            (int64_t)zoc5 * zoz, Z_SHIFT) +
            zoc4) * zoz, Z_SHIFT) + zoc3) * zoz, Z_SHIFT) +
            zoc2) * zoz, Z_SHIFT) + zoc1) * zoz, Z_SHIFT) + zoc0, 1);
#elif defined(Z_COEFF_INTERP_OPT4X)
        int32_t zoz, zoe1, zoe2, zoe3, zoo1, zoo2, zoo3;
        int32_t zoc0, zoc1, zoc2, zoc3, zoc4, zoc5;

        /* 6-point, 5th-order Optimal 4x */
        zoz = z - (Z_ONE >> 1);
        zoe1 = z_coeff[1] + z_coeff[0];
        zoe2 = z_coeff[2] + z_coeff[-1];
        zoe3 = z_coeff[3] + z_coeff[-2];
        zoo1 = z_coeff[1] - z_coeff[0];
        zoo2 = z_coeff[2] - z_coeff[-1];
        zoo3 = z_coeff[3] - z_coeff[-2];

        zoc0 = Z_RSHIFT((0x351db485LL * zoe1) + (0x0aaddc70LL * zoe2) +
            (0x00346f09LL * zoe3), 30);
        zoc1 = Z_RSHIFT((0x287b0c7bLL * zoo1) + (0x1b221c6cLL * zoo2) +
            (0x01395113LL * zoo3), 30);
        zoc2 = Z_RSHIFT((-0x1a04d042LL * zoe1) + (0x1706eee6LL * zoe2) +
            (0x02fde18bLL * zoe3), 30);
        zoc3 = Z_RSHIFT((-0x19de2a03LL * zoo1) + (0x0240f51cLL * zoo2) +
            (0x03d26db6LL * zoo3), 30);
        zoc4 = Z_RSHIFT((0x053fcc86LL * zoe1) + (-0x07de7f90LL * zoe2) +
            (0x029eb247LL * zoe3), 30);
        zoc5 = Z_RSHIFT((0x08753a11LL * zoo1) + (-0x042a9fd6LL * zoo2) +
            (0x00ce1b7aLL * zoo3), 30);

        coeff = Z_RSHIFT(Z_RSHIFT((Z_RSHIFT((Z_RSHIFT((Z_RSHIFT((Z_RSHIFT(
            (int64_t)zoc5 * zoz, Z_SHIFT) +
            zoc4) * zoz, Z_SHIFT) + zoc3) * zoz, Z_SHIFT) +
            zoc2) * zoz, Z_SHIFT) + zoc1) * zoz, Z_SHIFT) + zoc0, 1);
#elif defined(Z_COEFF_INTERP_OPT2X)
        int32_t zoz, zoe1, zoe2, zoe3, zoo1, zoo2, zoo3;
        int32_t zoc0, zoc1, zoc2, zoc3, zoc4, zoc5;

        /* 6-point, 5th-order Optimal 2x */
        zoz = z - (Z_ONE >> 1);
        zoe1 = z_coeff[1] + z_coeff[0];
        zoe2 = z_coeff[2] + z_coeff[-1];
        zoe3 = z_coeff[3] + z_coeff[-2];
        zoo1 = z_coeff[1] - z_coeff[0];
        zoo2 = z_coeff[2] - z_coeff[-1];
        zoo3 = z_coeff[3] - z_coeff[-2];

        zoc0 = Z_RSHIFT((0x33db6dfbLL * zoe1) + (0x0bd7a0c5LL * zoe2) +
            (0x004cf101LL * zoe3), 30);
        zoc1 = Z_RSHIFT((0x24475eecLL * zoo1) + (0x1bc7bad5LL * zoo2) +
            (0x01ad079bLL * zoo3), 30);
        zoc2 = Z_RSHIFT((-0x187dc0d1LL * zoe1) + (0x14b86ed1LL * zoe2) +
            (0x03c559d4LL * zoe3), 30);
        zoc3 = Z_RSHIFT((-0x15156c27LL * zoo1) + (-0x0032a45cLL * zoo2) +
            (0x0459df8eLL * zoo3), 30);
        zoc4 = Z_RSHIFT((0x04ec30fbLL * zoe1) + (-0x075003d1LL * zoe2) +
            (0x0263b26bLL * zoe3), 30);
        zoc5 = Z_RSHIFT((0x0586e7eaLL * zoo1) + (-0x024ebf05LL * zoo2) +
            (0x0031dcf2LL * zoo3), 30);

        coeff = Z_RSHIFT(Z_RSHIFT((Z_RSHIFT((Z_RSHIFT((Z_RSHIFT((Z_RSHIFT(
            (int64_t)zoc5 * zoz, Z_SHIFT) +
            zoc4) * zoz, Z_SHIFT) + zoc3) * zoz, Z_SHIFT) +
            zoc2) * zoz, Z_SHIFT) + zoc1) * zoz, Z_SHIFT) + zoc0, 1);
#else
#error "Interpolation type screwed!"
#endif

#if Z_POLYPHASE_COEFF_SHIFT > 0
        coeff = Z_RSHIFT(coeff, Z_POLYPHASE_COEFF_SHIFT);
#endif
        return (coeff);
}

static int
z_resampler_build_polyphase(struct z_info *info)
{
        int32_t *z_coeff, *z_pcoeff, alpha, i, z;

        /* Let this be here first. */
        if (info->z_pcoeff != NULL) {
                free(info->z_pcoeff, M_DEVBUF);
                info->z_pcoeff = NULL;
        }

        if (feeder_rate_polyphase_max < 1)
                return (ENOTSUP);

        if (((int64_t)info->z_size * info->z_gy * 2) >
            feeder_rate_polyphase_max) {
#ifndef _KERNEL
                fprintf(stderr, "Polyphase entries exceed: [%d/%d] %jd > %d\n",
                    info->z_gx, info->z_gy,
                    (intmax_t)info->z_size * info->z_gy * 2,
                    feeder_rate_polyphase_max);
#endif
                return (E2BIG);
        }

        info->z_pcoeff = malloc(sizeof(int32_t) *
            info->z_size * info->z_gy * 2, M_DEVBUF, M_NOWAIT | M_ZERO);
        if (info->z_pcoeff == NULL)
                return (ENOMEM);

        for (alpha = 0; alpha < info->z_gy; alpha++) {
                z = alpha * info->z_dx;
                z_coeff = info->z_coeff;
                z_pcoeff = info->z_pcoeff + (alpha * info->z_size * 2) +
                    info->z_size;
                for (i = info->z_size; i != 0; i--) {
                        z_coeff += z >> Z_SHIFT;
                        z &= Z_MASK;
                        *z_pcoeff++ = z_coeff_interpolate(z, z_coeff);
                        z += info->z_dy;
                }
                z = info->z_dy - (alpha * info->z_dx);
                z_coeff = info->z_coeff;
                z_pcoeff = info->z_pcoeff + (alpha * info->z_size * 2) +
                    info->z_size - 1;
                for (i = info->z_size; i != 0; i--) {
                        z_coeff += z >> Z_SHIFT;
                        z &= Z_MASK;
                        *z_pcoeff-- = z_coeff_interpolate(z, z_coeff);
                        z += info->z_dy;
                }
        }
        
#ifndef _KERNEL
        fprintf(stderr, "Polyphase: [%d/%d] %d entries\n",
            info->z_gx, info->z_gy, info->z_size * info->z_gy * 2);
#endif

        return (0);
}

static int
z_resampler_setup(struct pcm_feeder *f)
{
        struct z_info *info;
        int64_t gy2gx_max, gx2gy_max;
        uint32_t format;
        int32_t align, bufsz, i, startdrift, z_scale;
        int adaptive;

        info = f->data;
        z_resampler_reset(info);

        if (info->src == info->dst)
                return (0);

        /* Shrink by greatest common divisor. */
        i = z_gcd(info->src, info->dst);
        info->z_gx = info->src / i;
        info->z_gy = info->dst / i;

        /* Too big, or too small. Bail out. */
        if (!(Z_FACTOR_SAFE(info->z_gx) && Z_FACTOR_SAFE(info->z_gy)))
                return (EINVAL);

        format = f->desc->in;
        adaptive = 0;
        z_scale = 0;

        /*
         * Setup everything: filter length, conversion factor, etc.
         */
        if (Z_IS_SINC(info)) {
                /*
                 * Downsampling, or upsampling scaling factor. As long as the
                 * factor can be represented by a fraction of 1 << Z_SHIFT,
                 * we're pretty much in business. Scaling is not needed for
                 * upsampling, so we just slap Z_ONE there.
                 */
                if (Z_SRC_DOWN(info)) {
                        /*
                         * If the downsampling ratio is beyond sanity,
                         * enable semi-adaptive mode. Although handling
                         * extreme ratio is possible, the result of the
                         * conversion is just pointless, unworthy,
                         * nonsensical noises, etc.
                         */
                        if (Z_FACTOR_DOWN(info) > Z_SINC_DOWNMAX)
                                z_scale = Z_ONE / Z_SINC_DOWNMAX;
                        else
                                z_scale = ((uint64_t)info->z_gy << Z_SHIFT) /
                                    info->z_gx;
                } else
                        z_scale = Z_ONE;

                /*
                 * This is actually impossible, unless anything above
                 * overflow.
                 */
                if (z_scale < 1)
                        return (E2BIG);

                /*
                 * Calculate sample time/coefficients index drift. It is
                 * a constant for upsampling, but downsampling require
                 * heavy duty filtering with possible too long filters.
                 * If anything goes wrong, revisit again and enable
                 * adaptive mode.
                 */
z_setup_adaptive_sinc:
                if (info->z_pcoeff != NULL) {
                        free(info->z_pcoeff, M_DEVBUF);
                        info->z_pcoeff = NULL;
                }

                if (adaptive == 0) {
                        info->z_dy = z_scale << Z_DRIFT_SHIFT;
                        if (info->z_dy < 1)
                                return (E2BIG);
                        info->z_scale = z_scale;
                } else {
                        info->z_dy = Z_FULL_ONE;
                        info->z_scale = Z_ONE;
                }

#if 0
#define Z_SCALE_DIV     10000
#define Z_SCALE_LIMIT(s, v)                                             \
        ((((uint64_t)(s) * (v)) + (Z_SCALE_DIV >> 1)) / Z_SCALE_DIV)

                info->z_scale = Z_SCALE_LIMIT(info->z_scale, 9780);
#endif

                /* Smallest drift increment. */
                info->z_dx = info->z_dy / info->z_gy;

                /*
                 * Overflow or underflow. Try adaptive, let it continue and
                 * retry.
                 */
                if (info->z_dx < 1) {
                        if (adaptive == 0) {
                                adaptive = 1;
                                goto z_setup_adaptive_sinc;
                        }
                        return (E2BIG);
                }

                /*
                 * Round back output drift.
                 */
                info->z_dy = info->z_dx * info->z_gy;

                for (i = 0; i < Z_COEFF_TAB_SIZE; i++) {
                        if (Z_SINC_COEFF_IDX(info) != i)
                                continue;
                        /*
                         * Calculate required filter length and guard
                         * against possible abusive result. Note that
                         * this represents only 1/2 of the entire filter
                         * length.
                         */
                        info->z_size = z_resampler_sinc_len(info);

                        /*
                         * Multiple of 2 rounding, for better accumulator
                         * performance.
                         */
                        info->z_size &= ~1;

                        if (info->z_size < 2 || info->z_size > Z_SINC_MAX) {
                                if (adaptive == 0) {
                                        adaptive = 1;
                                        goto z_setup_adaptive_sinc;
                                }
                                return (E2BIG);
                        }
                        info->z_coeff = z_coeff_tab[i].coeff + Z_COEFF_OFFSET;
                        info->z_dcoeff = z_coeff_tab[i].dcoeff;
                        break;
                }

                if (info->z_coeff == NULL || info->z_dcoeff == NULL)
                        return (EINVAL);
        } else if (Z_IS_LINEAR(info) && !Z_IS_LINEAR_ZOH(info)) {
                /*
                 * Don't put much effort if we're doing linear interpolation.
                 * Just center the interpolation distance within Z_LINEAR_ONE,
                 * and be happy about it.
                 */
                info->z_dx = Z_LINEAR_ONE / info->z_gy;
                info->z_size = 2;
        }

        /*
         * We're safe for now, lets continue.. Look for our resampler
         * depending on configured format and quality.
         */
        switch (Z_RESAMPLER_IDX(info)) {
        case Z_RESAMPLER_ZOH:
                info->z_resample = z_feed_zoh;
                break;
        case Z_RESAMPLER_LINEAR:
                info->z_resample = Z_IS_LINEAR_ZOH(info) ? z_feed_zoh :
                    z_feed_linear;
                break;
        case Z_RESAMPLER_SINC:
                if (adaptive == 0 && z_resampler_build_polyphase(info) == 0)
                        info->z_resample = z_feed_sinc_polyphase;
                else if (feeder_rate_polyphase_max < 0)
                        info->z_resample = z_feed_sinc_fast;
                else
                        info->z_resample = z_feed_sinc;
                break;
        default:
                return (EINVAL);
                break;
        }

        align = info->channels * PCM_32_BPS;

        /*
         * Calculate largest value that can be fed into z_gy2gx() and
         * z_gx2gy() without causing (signed) 32bit overflow. z_gy2gx() will
         * be called early during feeding process to determine how much input
         * samples that is required to generate requested output, while
         * z_gx2gy() will be called just before samples filtering /
         * accumulation process based on available samples that has been
         * calculated using z_gx2gy().
         *
         * Now that is damn confusing, I guess ;-) .
         */
        gy2gx_max = (((uint64_t)info->z_gy * INT32_MAX) - info->z_gy + 1) /
            info->z_gx;

        if ((gy2gx_max * align) > SND_FXDIV_MAX)
                gy2gx_max = SND_FXDIV_MAX / align;

        if (gy2gx_max < 1)
                return (E2BIG);

        gx2gy_max = (((uint64_t)info->z_gx * INT32_MAX) - info->z_gy) /
            info->z_gy;

        if (gx2gy_max > INT32_MAX)
                gx2gy_max = INT32_MAX;

        if (gx2gy_max < 1)
                return (E2BIG);

        /*
         * Ensure that z_gy2gx() at its largest possible calculated value
         * (alpha = 0) will not cause overflow further late during z_gx2gy()
         * stage.
         */
        if (z_gy2gx(info, gy2gx_max) > _Z_GCAST(gx2gy_max))
                return (E2BIG);

        info->z_maxfeed = gy2gx_max * align;

        startdrift = z_start_drift(info);
        info->z_startdrift = startdrift;
        info->z_alphadrift = z_drift(info, startdrift, 1);

#ifdef Z_STRESS_TEST
        info->z_full = z_history(info);
#else
        /*
         * Extreme downsampling does not deserve buffering optimization
         * since startdrift might become too large.
         */
        if (Z_FACTOR_DOWN(info) > Z_DOWNMAX)
                info->z_full = 0;
        else {
                /*
                 * First, optimize delay buffer size by allocating enough
                 * samples so that at least it can accomodate multiple
                 * passes of resampling routine in a single feed.
                 */
                info->z_full = z_roundpow2(z_history(info) + startdrift);

                /* Test against multi-pass without ^2 rounding. */
                bufsz = z_history(info) + startdrift;
                if ((bufsz * align) > (info->z_full * align))
                        info->z_full = bufsz;
        }

        /* Test against single-pass with ^2 rounding. */
        bufsz = z_roundpow2(z_history(info));
        if ((bufsz * align) > (info->z_full * align))
                info->z_full = bufsz;

        /* Test against single-pass without ^2 rounding */
        bufsz = z_history(info);
        if ((bufsz * align) > (info->z_full * align))
                info->z_full = bufsz;

        /*
         * Too big to be true, and overflowing left and right like mad ..
         */
        if ((info->z_full * align) < 1) {
                if (adaptive == 0 && Z_IS_SINC(info)) {
                        adaptive = 1;
                        goto z_setup_adaptive_sinc;
                }
                return (E2BIG);
        }

        /*
         * Increase full buffer size if its too small to reduce cyclic
         * buffer shifting in main conversion/feeder loop.
         */
        while (info->z_full < Z_RESERVOIR_MAX &&
            (info->z_full - z_history(info)) < Z_RESERVOIR)
                info->z_full <<= 1;
#endif

        /* Initialize buffer position. */
        info->z_pos = z_history(info) - 1;
        info->z_start = info->z_pos - startdrift;

        /*
         * Allocate or reuse delay line buffer, whichever makes sense.
         */
        bufsz = info->z_full * align;

        if (info->z_delay == NULL || info->z_alloc < bufsz ||
            bufsz <= (info->z_alloc >> 1)) {
                if (info->z_delay != NULL)
                        free(info->z_delay, M_DEVBUF);
                info->z_delay = malloc(bufsz, M_DEVBUF, M_NOWAIT | M_ZERO);
                if (info->z_delay == NULL)
                        return (ENOMEM);
                info->z_alloc = bufsz;
        }

        /*
         * Zero out head of buffer to avoid pops and clicks.
         */
        memset(info->z_delay, 0, info->z_pos * align);

#ifdef Z_DIAGNOSTIC
        /*
         * XXX Debuging mess !@#$%^
         */
#define dumpz(x)        fprintf(stderr, "\t%12s = %10u : %-11d\n",     \
                            "z_"__STRING(x), (uint32_t)info->z_##x,  \
                            (int32_t)info->z_##x)
        fprintf(stderr, "\n%s():\n", __func__);
        fprintf(stderr, "\tchannels=%d, bps=%d, format=0x%08x, quality=%d\n",
            info->channels, info->bps, format, info->quality);
        fprintf(stderr, "\t%d (%d) -> %d (%d), ",
            info->src, info->rsrc, info->dst, info->rdst);
        fprintf(stderr, "[%d/%d]\n", info->z_gx, info->z_gy);
        fprintf(stderr, "\tminreq=%d, ", z_gy2gx(info, 1));
        if (adaptive != 0)
                z_scale = Z_ONE;
        fprintf(stderr, "factor=0x%08x/0x%08x (%f)\n",
            z_scale, Z_ONE, (double)z_scale / Z_ONE);
        fprintf(stderr, "\tbase_length=%d, ", Z_SINC_BASE_LEN(info));
        fprintf(stderr, "adaptive=%s\n", (adaptive != 0) ? "YES" : "NO");
        dumpz(size);
        dumpz(alloc);
        if (info->z_alloc < 1024)
                fprintf(stderr, "\t%15s%10d Bytes\n",
                    "", info->z_alloc);
        else if (info->z_alloc < (1024 << 10))
                fprintf(stderr, "\t%15s%10d KBytes\n",
                    "", info->z_alloc >> 10);
        else if (info->z_alloc < (1024 << 20))
                fprintf(stderr, "\t%15s%10d MBytes\n",
                    "", info->z_alloc >> 20);
        else
                fprintf(stderr, "\t%15s%10d GBytes\n",
                    "", info->z_alloc >> 30);
        fprintf(stderr, "\t%12s   %10d (min output samples)\n",
            "",
            (int32_t)z_gx2gy(info, info->z_full - z_history(info)));
        fprintf(stderr, "\t%12s   %10d (min allocated output samples)\n",
            "",
            (int32_t)z_gx2gy(info, (info->z_alloc / align) - z_history(info)));
        fprintf(stderr, "\t%12s = %10d\n",
            "z_gy2gx()", (int32_t)z_gy2gx(info, 1));
        fprintf(stderr, "\t%12s = %10d -> z_gy2gx() -> %d\n",
            "Max", (int32_t)gy2gx_max, (int32_t)z_gy2gx(info, gy2gx_max));
        fprintf(stderr, "\t%12s = %10d\n",
            "z_gx2gy()", (int32_t)z_gx2gy(info, 1));
        fprintf(stderr, "\t%12s = %10d -> z_gx2gy() -> %d\n",
            "Max", (int32_t)gx2gy_max, (int32_t)z_gx2gy(info, gx2gy_max));
        dumpz(maxfeed);
        dumpz(full);
        dumpz(start);
        dumpz(pos);
        dumpz(alphadrift);
        dumpz(startdrift);
        dumpz(scale);
        fprintf(stderr, "\t%12s   %10f\n", "",
            (double)info->z_scale / Z_ONE);
        dumpz(dx);
        fprintf(stderr, "\t%12s   %10f\n", "",
            (double)info->z_dx / info->z_dy);
        dumpz(dy);
        fprintf(stderr, "\t%12s   %10d (drift step)\n", "",
            info->z_dy >> Z_SHIFT);
        fprintf(stderr, "\t%12s   %10d (scaling differences)\n", "",
            (z_scale << Z_DRIFT_SHIFT) - info->z_dy);
        fprintf(stderr, "\t%12s = %u bytes\n",
            "intpcm32_t", sizeof(intpcm32_t));
        fprintf(stderr, "\t%12s = 0x%08x, smallest=%.16lf\n",
            "Z_ONE", Z_ONE, (double)1.0 / (double)Z_ONE);
#endif

        return (0);
}

static int
z_resampler_set(struct pcm_feeder *f, int what, int32_t value)
{
        struct z_info *info;
        int32_t oquality;

        info = f->data;

        switch (what) {
        case Z_RATE_SRC:
                if (value < feeder_rate_min || value > feeder_rate_max)
                        return (E2BIG);
                if (value == info->rsrc)
                        return (0);
                info->rsrc = value;
                break;
        case Z_RATE_DST:
                if (value < feeder_rate_min || value > feeder_rate_max)
                        return (E2BIG);
                if (value == info->rdst)
                        return (0);
                info->rdst = value;
                break;
        case Z_RATE_QUALITY:
                if (value < Z_QUALITY_MIN || value > Z_QUALITY_MAX)
                        return (EINVAL);
                if (value == info->quality)
                        return (0);
                /*
                 * If we failed to set the requested quality, restore
                 * the old one. We cannot afford leaving it broken since
                 * passive feeder chains like vchans never reinitialize
                 * itself.
                 */
                oquality = info->quality;
                info->quality = value;
                if (z_resampler_setup(f) == 0)
                        return (0);
                info->quality = oquality;
                break;
        case Z_RATE_CHANNELS:
                if (value < SND_CHN_MIN || value > SND_CHN_MAX)
                        return (EINVAL);
                if (value == info->channels)
                        return (0);
                info->channels = value;
                break;
        default:
                return (EINVAL);
                break;
        }

        return (z_resampler_setup(f));
}

static int
z_resampler_get(struct pcm_feeder *f, int what)
{
        struct z_info *info;

        info = f->data;

        switch (what) {
        case Z_RATE_SRC:
                return (info->rsrc);
                break;
        case Z_RATE_DST:
                return (info->rdst);
                break;
        case Z_RATE_QUALITY:
                return (info->quality);
                break;
        case Z_RATE_CHANNELS:
                return (info->channels);
                break;
        default:
                break;
        }

        return (-1);
}

static int
z_resampler_init(struct pcm_feeder *f)
{
        struct z_info *info;
        int ret;

        if (f->desc->in != f->desc->out ||
            AFMT_ENCODING(f->desc->in) != AFMT_S32_NE)
                return (EINVAL);

        info = malloc(sizeof(*info), M_DEVBUF, M_NOWAIT | M_ZERO);
        if (info == NULL)
                return (ENOMEM);

        info->rsrc = Z_RATE_DEFAULT;
        info->rdst = Z_RATE_DEFAULT;
        info->quality = feeder_rate_quality;
        info->channels = AFMT_CHANNEL(f->desc->in);

        f->data = info;

        ret = z_resampler_setup(f);
        if (ret != 0) {
                if (info->z_pcoeff != NULL)
                        free(info->z_pcoeff, M_DEVBUF);
                if (info->z_delay != NULL)
                        free(info->z_delay, M_DEVBUF);
                free(info, M_DEVBUF);
                f->data = NULL;
        }

        return (ret);
}

static int
z_resampler_free(struct pcm_feeder *f)
{
        struct z_info *info;

        info = f->data;
        if (info != NULL) {
                if (info->z_pcoeff != NULL)
                        free(info->z_pcoeff, M_DEVBUF);
                if (info->z_delay != NULL)
                        free(info->z_delay, M_DEVBUF);
                free(info, M_DEVBUF);
        }

        f->data = NULL;

        return (0);
}

/*
 * The math, though pretty much elementary, looks scary enough.  But just
 * how accurate is this entire conversion and buffering routine?
 *
 * 1) It is accurate to the point that you don't need more than the total
 *    number of history samples that need to be kept, which means at least:
 *
 *      sinc = total number of taps (z_size * 2, including 1 current sample)
 *    linear = current and previous sample (1 + 1 = 2)
 *       zoh = current sample (1)
 *
 *    However, small buffer means that a lot of buffer shifting need to be
 *    done and will effect the entire conversion performance. It is also the
 *    best way to stress things out.
 *
 * 2) For zero-order-hold (ZOH), symmetric conversion of up and down is
 *    lossless. Up followed by down symmetrically, not the other way round.
 *    Regardless of whatever ratios involved. Lossless. No kidding.
 *
 * ...and thus, it is the very definition of accurate.
 */
static uint32_t
z_resampler_feed_internal(struct pcm_feeder *f, struct pcm_channel *c,
    intpcm_t *b, uint32_t count, void *source)
{
        struct z_info *info;
        int32_t reqout, ocount, reqin, align, fetch, fetched;
        intpcm_t *dst;

        info = f->data;
        if (info->z_resample == NULL)
                return (z_feed(f->source, c, (uint8_t *)b, count, source));

        /*
         * Calculate sample size alignment and amount of sample output.
         * We will do everything in sample domain, but at the end we
         * will jump back to byte domain.
         */
        align = info->channels * PCM_32_BPS;
        ocount = SND_FXDIV(count, align);
        if (ocount == 0)
                return (0);

        /*
         * Calculate amount of input samples that is needed to generate
         * exact amount of output.
         */
        reqin = z_gy2gx(info, ocount) - z_fetched(info);

        dst = b;

        do {
                if (reqin != 0) {
                        fetch = z_min(z_free(info), reqin);
                        if (fetch == 0) {
#define z_future_startdrift()   (((info->z_alpha + info->z_alphadrift) <        \
                                info->z_gy) ? info->z_startdrift :              \
                                (info->z_startdrift - 1))
#define z_future_start()        (info->z_start + z_future_startdrift())
#define z_future_begin()        (z_future_start() - z_history(info) + 1)
                                /*
                                 * No more free spaces, so wind required
                                 * samples back to the head of delay line
                                 * in byte domain. The winding process is
                                 * done in a such way that the future begin
                                 * offset will fall exactly at index 0 of
                                 * delay/history buffer.
                                 *
                                 *  zfb = future z_start begin offset
                                 *        (constant = 0)
                                 *  zfs = future z_start
                                 *  zfp = future z_pos
                                 * zfsd = future startdrift
                                 *  zh  = z_history
                                 *    f = current fetched buffer (constant)
                                 *   zs = current z_start
                                 *
                                 *          zfb = 0
                                 * zfs - zh + 1 = 0
                                 *          zfs = zh - 1
                                 *    zs + zfsd = zh - 1
                                 *           zs = zh - zfsd - 1
                                 *  zfp - f - 1 = zh - zfsd - 1
                                 *          zfp = zh - zfsd - 1 + f + 1
                                 *          zfp = zh - zfsd + f
                                 */

                                /* Get the amount of fetched samples. */
                                fetched = z_fetched(info);

                                /* Calculate future z_pos. */
                                info->z_pos = z_history(info) -
                                    z_future_startdrift() + fetched;

                                /*
                                 *      z_pos <= 0 : The current fetched
                                 *                   samples will fall far
                                 *                   behind the delay buffer.
                                 *                   Dispose it entirely.
                                 *                   Set z_pos to 0.
                                 *
                                 *  z_pos < z_full : Early part of fetched
                                 *                   samples will fall behind
                                 *                   the delay buffer. Move
                                 *                   remaining to the head.
                                 *
                                 * z_pos == z_full : Do nothing. The entire
                                 *                   fetched samples might be
                                 *                   reused.
                                 *
                                 *  z_pos > z_full : Mathematically impossible.
                                 *
                                 * It is obvious that z_start might become
                                 * negative at certain point, but the drift
                                 * process will reposition it within z_full
                                 * boundary just before the resampling routine
                                 * begin.
                                 */
                                if (info->z_pos != info->z_full) {
#if defined(Z_DIAGNOSTIC) || defined(Z_STRESS_TEST)
                                        /* Mathematically impossible. */
                                        if (info->z_pos > info->z_full)
                                                fprintf(stderr, "OVERFLOW: "
                                                    "z_pos=%d >= z_full\n",
                                                    info->z_pos, info->z_full);
#endif
                                        if (info->z_pos > 0)
                                                z_copy(
                                                    (uint8_t *)info->z_delay +
                                                    ((info->z_full -
                                                    info->z_pos) * align),
                                                    info->z_delay,
                                                    info->z_pos * align);
                                        else
                                                info->z_pos = 0;
                                        info->z_start = info->z_pos -
                                            fetched - 1;
                                        fetch = z_min(z_free(info), reqin);
                                }
#ifdef Z_DIAGNOSTIC
                                if (1) {
                                        static uint32_t kk = 0;
                                        fprintf(stderr,
                                            "Buffer Move: "
                                            "start=%d fetched=%d cp=%d "
                                            "cycle=%u [%u]\r",
                                            info->z_full - info->z_pos,
                                            fetched, info->z_pos,
                                            info->z_cycle, ++kk);
                                }
                                info->z_cycle = 0;
#endif
                        }
                        if (fetch != 0) {
                                /*
                                 * Fetch in byte domain and jump back
                                 * to sample domain.
                                 */
#if defined(Z_DIAGNOSTIC) || defined(Z_STRESS_TEST)
                                /* Mathematically impossible. */
                                if (fetch < 0)
                                        fprintf(stderr, "NEGATIVE FETCH: %d\n",
                                            fetch);
                                if (info->z_pos >= info->z_full)
                                        fprintf(stderr, "OVERFLOW: "
                                            "z_pos=%d >= z_full=%d\n",
                                            info->z_pos, info->z_full);
                                if ((info->z_full - info->z_pos) < fetch)
                                        fprintf(stderr, "IMPOSSIBLE FETCH: "
                                            "fetch=%d, z_pos=%d, z_full=%d\n",
                                            fetch, info->z_pos, info->z_full);
#endif
                                fetched = SND_FXDIV(z_feed(f->source, c,
                                    (uint8_t *)info->z_delay +
                                    (info->z_pos * align),
                                    fetch * align, source), align);
                                /*
                                 * Prepare to convert fetched buffer,
                                 * or mark us done if we cannot fulfill
                                 * the request.
                                 */
                                reqin -= fetched;
                                info->z_pos += fetched;
                                if (fetched != fetch)
                                        reqin = 0;
                        }
                }

                reqout = z_min(z_gx2gy(info, z_fetched(info)), ocount);
                if (reqout != 0) {
#ifdef Z_DIAGNOSTIC
                        info->z_cycle += reqout;
#endif
                        ocount -= reqout;

                        /*
                         * Drift.. drift.. drift..
                         *
                         * Notice that there are 2 methods of doing the drift
                         * operations: The former is much cleaner (in a sense
                         * of mathematical readings of my eyes), but slower
                         * due to integer division in z_gy2gx(). Nevertheless,
                         * both should give the same exact accurate drifting
                         * results, so the later is favourable.
                         *
                         * startdrift = z_gy2gx(info, 1);
                         * alphadrift = z_drift(info, startdrift, 1);
                         * info->z_start += startdrift;
                         * info->z_alpha += alphadrift;
                         */
                        info->z_resample(info, dst, reqout);
                        dst += info->channels * reqout;
                }
        } while (reqin != 0 && ocount != 0);

        /*
         * Back to byte domain..
         */
        return (PCM_32_BPS * (dst - b));
}

static int
z_resampler_feed(struct pcm_feeder *f, struct pcm_channel *c, uint8_t *b,
    uint32_t count, void *source)
{
        uint32_t feed, maxfeed, left;

        /*
         * Split count to smaller chunks to avoid possible 32bit overflow.
         */
        maxfeed = ((struct z_info *)(f->data))->z_maxfeed;
        left = count;

        do {
                feed = z_resampler_feed_internal(f, c, (intpcm_t *)b,
                    z_min(maxfeed, left), source);
                b += feed;
                left -= feed;
        } while (left != 0 && feed != 0);

        return (count - left);
}

static struct pcm_feederdesc feeder_rate_desc[] = {
        { FEEDER_RATE, AFMT_S32_NE, AFMT_S32_NE, 0, 0 },
        { 0, 0, 0, 0, 0 },
};

static kobj_method_t feeder_rate_methods[] = {
        KOBJMETHOD(feeder_init,         z_resampler_init),
        KOBJMETHOD(feeder_free,         z_resampler_free),
        KOBJMETHOD(feeder_set,          z_resampler_set),
        KOBJMETHOD(feeder_get,          z_resampler_get),
        KOBJMETHOD(feeder_feed,         z_resampler_feed),
        KOBJMETHOD_END
};

FEEDER_DECLARE(feeder_rate, NULL);