/*
 * CDDL HEADER START
 *
 * The contents of this file are subject to the terms of the
 * Common Development and Distribution License, Version 1.0 only
 * (the "License").  You may not use this file except in compliance
 * with the License.
 *
 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
 * or http://www.opensolaris.org/os/licensing.
 * See the License for the specific language governing permissions
 * and limitations under the License.
 *
 * When distributing Covered Code, include this CDDL HEADER in each
 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
 * If applicable, add the following below this CDDL HEADER, with the
 * fields enclosed by brackets "[]" replaced with your own identifying
 * information: Portions Copyright [yyyy] [name of copyright owner]
 *
 * CDDL HEADER END
 */
/*
 * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
 * Use is subject to license terms.
 */

/*      Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T */
/*        All Rights Reserved   */

#include <unistd.h>
#include <stdlib.h>
#include <stdio.h>

#define BYTE 8
#define QW 1            /* width of bas-q digit in bits */

/*
 * this stuff should be local and hidden; it was made
 * accessible outside for dirty reasons: 20% faster spell
 */
#include "huff.h"
struct huff huffcode;

/*
 * Infinite Huffman code
 *
 * Let the messages be exponentially distributed with ratio r:
 *      P {message k} = r^k*(1-r),      k = 0, 1, ...
 * Let the messages be coded in base q, and suppose
 *      r^n = 1/q
 * If each decade(base q) contains n codes, then
 * the messages assigned to each decade will be q times
 * as probable as the next. Moreover the code for the tail of
 * the distribution after truncating one decade should look
 * just like the original, but longer by one leading digit q-1.
 *      q(z+n) = z + (q-1)q^w
 * where z is first code of decade, w is width of code, in shortest
 * full decade. Examples, base 2:
 *      r^1 = 1/2       r^5 = 1/2
 *      0               0110
 *      10              0111
 *      110             1000
 *      1110            1001
 *      ...             1010
 *                      10110
 *      w = 1, z = 0            w = 4, z = 0110
 * Rewriting slightly
 *      (q-1)z + q*n = (q-1)q^w
 * whence z is a multiple of q and n is a multiple of q-1. Let
 *      z = cq, n = d(q-1)
 * We pick w to be the least integer such that
 *      d = n/(q-1) <= q^(w-1)
 * Then solve for c
 *      c = q^(w-1) - d
 * If c is not zero, the first decade may be preceded by
 * even shorter (w-1)-digit codes 0, 1, ..., c-1. Thus
 * the example code with r^5 = 1/2 becomes
 *      000
 *      001
 *      010
 *      0110
 *      0111
 *      1000
 *      1001
 *      1010
 *      10110
 *      ...
 *      110110
 *      ...
 * The expected number of base-q digits in a codeword is then
 *      w - 1 + r^c/(1-r^n)
 * The present routines require q to be a power of 2
 */
/*
 * There is a lot of hanky-panky with left justification against
 * sign instead of simple left justification because
 * unsigned long is not available
 */
#define L (BYTE*(sizeof (long))-1)      /* length of signless long */
#define MASK (~((unsigned long)1<<L))   /* mask out sign */

/*
 * decode the prefix of word y (which is left justified against sign)
 * place mesage number into place pointed to by kp
 * return length (in bits) of decoded prefix or 0 if code is out of
 * range
 */
int
decode(long y, long *pk)
{
        int l;
        long v;
        if (y < cs) {
                *pk = y >> (long)(L+QW-w);
                return (w-QW);
        }
        for (l = w, v = v0; y >= qcs;
            y = ((unsigned long)y << QW) & MASK, v += n)
                if ((l += QW) > L)
                        return (0);
        *pk = v + (y>>(long)(L-w));
        return (l);
}

/*
 * encode message k and put result (right justified) into
 * place pointed to by py.
 * return length (in bits) of result,
 * or 0 if code is too long
 */

int
encode(long k, long *py)
{
        int l;
        long y;
        if (k < c) {
                *py = k;
                return (w-QW);
        }
        for (k -= c, y = 1, l = w; k >= n; k -= n, y <<= QW)
                if ((l += QW) > L)
                        return (0);
        *py = ((y-1)<<w) + cq + k;
        return (l);
}


/*
 * Initialization code, given expected value of k
 * E(k) = r/(1-r) = a
 * and given base width b
 * return expected length of coded messages
 */
static struct qlog {
        long p;
        double u;
} z;

static struct qlog
qlog(double x, double y, long p, double u)      /* find smallest p so x^p<=y */
{

        if (u/x <= y) {
                z.p = 0;
                z.u = 1;
        } else {
                z = qlog(x, y, p+p, u*u);
                if (u*z.u/x > y) {
                        z.p += p;
                        z.u *= u;
                }
        }
        return (z);
}

double
huff(float a)
{
        int i, q;
        long d, j;
        double r = a/(1.0 + a);
        double rc, rq;

        for (i = 0, q = 1, rq = r; i < QW; i++, q *= 2, rq *= rq)
                continue;
        rq /= r;        /* rq = r^(q-1) */
        (void) qlog(rq, 1./q, 1L, rq);
        d = z.p;
        n = d*(q-1);
        if (n != d * (q - 1))
                abort();        /* time to make n long */
        for (w = QW, j = 1; j < d; w += QW, j *= q)
                continue;
        c = j - d;
        cq = c*q;
        cs = cq<<(L-w);
        qcs = (((long)(q-1)<<w) + cq) << (L-QW-w);
        v0 = c - cq;
        for (i = 0, rc = 1; i < c; i++, rc *= r)        /* rc = r^c */
                continue;
        return (w + QW*(rc/(1-z.u) - 1));
}

void
whuff(void)
{
        (void) fwrite((char *) & huffcode, sizeof (huffcode), 1, stdout);
}

int
rhuff(FILE *f)
{
        return (read(fileno(f), (char *)&huffcode, sizeof (huffcode)) ==
            sizeof (huffcode));
}