/*
 * Copyright 1/5/95 William Kahan.  All Rights Reserved.
 *
 * Coded for Professor William Kahan by Zhishun Alex Liu.
 * The coder's current electronic mail address as of December 1994 is
 * Alex.Liu@Eng.Sun.COM.
 *
 * Use of this software is granted with the understanding that all recipients
 * should regard themselves as participants in an ongoing research project and
 * hence should feel obligated to report their experiences (good or bad) with
 * it to the author.
 */
#ifndef lint
static char sccsid[] = "@(#)distd.c	1.1 (BeEF) 1/5/95";
#endif
/*
 * distd(n, x): Given n > 0 and x0, x1, ..., xn-1, we replace them by n'
 * and x0', x1', ..., xn-1' in such a way as nearly minimizes n' yet keeps
 * x0' + x1' + ... + xn-1' == x0 + x1 + ... + xn-1 exactly, using only
 * standard WP operations upon WP variables.
 *
 * distd() returns n'.
 *
 *	The first step is to sort {xj} by magnitude to ensure that
 *	|x0| >= |x1| >= ... >= |xn-1| > 0.  The sorting procedure used
 *	here should depend upon the provenance of the data; if {xj}
 *	is in no special order then a Heap-sort is appropriate;
 *	otherwise a Bubble-sort or Merge-sort may run faster.
 *
 *	The Distillation process proper consists of repeated summations
 *	alternating in direction from small-to-big, then big-to-small.
 *	Each summation leaves Sumj{xj} unchanged.  Each summation starts
 *	with x0 + x1 + ... + xkl and xkr + ... + xn-2 + xn-1 already
 *	distilled, and tracks changes in kl and kr.
 */
#include "sys.h"
			/* insertion sort; into decreasing order in magnitude */
#define	ISORTD(n, a)	{\
	register int i, j;\
	GENERIC gtmp;\
\
	for (i = 1; i <= n - 1; i++) {\
		gtmp = a[i];\
		j = i - 1;\
		while (j >= 0 && GABS(a[j]) < GABS(gtmp)) {\
			a[j + 1] = a[j];\
			j--;\
		}\
		a[j + 1] = gtmp;\
	}\
}
			/* sort a; into decreasing order in magnitude */
#define	GSORTD(n, a)	{\
	ISORTD(n, a)\
	while (n > 0 && a[n - 1] == ZERO)\
		n--;\
}
			/* Z := p + q = z + zeta where z = p + q rounded */
#define	DIST2(p, q, z, zeta)	{\
	GENERIC P = p, Q = q, gtmp;\
\
	if (GABS(P) < GABS(Q)) {\
		gtmp = P;\
		P = Q;\
		Q = gtmp;\
	}\
	z = P + Q;\
	zeta = P - z;\
	zeta += Q;\
}

int
#ifdef __STDC__
distd(int n, GENERIC *x)
#else
distd(n, x)
int n;
GENERIC *x;
#endif
{
	register int kl, kr, i, j, k;
	GENERIC p, q, z, zeta;

	GSORTD(n, x)
	kl = 0;
	kr = n - 1;
	while (kl < kr) {
		j = kr;				/* Backward pass */
		p = x[j];
		k = n;
		for (i = j - 1; i >= 0; i--) {
			if (i < kl && j + 1 < k) {
				x[j--] = p;
				p = x[i];
			}
			else {
				q = x[i];
				DIST2(p, q, z, zeta)
				if (zeta == ZERO)
					k = j;
				else {
					x[j] = zeta;
					if (z != q)
						k = j;
					j--;
					if (k == n)
						kr = j;
				}
				p = z;
			}
		}
		kl = k;
		if (kr < n - 1)
			kr++;
		if (p == ZERO)
			j++;
		else
			x[j] = p;
		if (j > 0) {
			n -= j;
			kl -= j;
			kr -= j;
			for (i = 0; i <= n - 1; i++)
				x[i] = x[i + j];
		}
		if (kl < 0)			/* forces valid kl */
			kl = 0;
		if (kl < kr) {			/* Forward pass */
			j = kl;
			p = x[j];
			k = -1;
			for (i = j + 1; i <= n - 1; i++) {
				if (i > kr && j - 1 > k) {
					x[j++] = p;
					p = x[i];
				}
				else {
					q = x[i];
					DIST2(p, q, z, zeta)
					if (zeta != q)
						k = j;
					if (zeta == ZERO)
						p = z;
					else {
						x[j++] = z;
						p = zeta;
						if (k <= -1)
							kl = j;
					}
				}
			}
			kr = k;
			if (kl > 0)
				kl--;
			if (p == ZERO)
				n = j;
			else {
				x[j] = p;
				n = j + 1;
			}
			if (kr >= n)		/* forces valid kr */
				kr = n - 1;
		}
	}
	return n;
}
#if 0
			/* shell sort; into decreasing order in magnitude */
#define	SSORTD(n, a)	{\
	register int i, j, gap = n;\
	GENERIC gtmp;\
\
	while ((gap >>= 1) > 0) {\
		for (i = gap; i <= n - 1; i++) {\
			gtmp = a[i];\
			j = i - gap;\
			while (j >= 0 && GABS(a[j]) < GABS(gtmp)) {\
				a[j + gap] = a[j];\
				j -= gap;\
			}\
			a[j + gap] = gtmp;\
		}\
	}\
}
			/* heap sort; into decreasing order in magnitude */
#define	HSORTD(n, a)	{\
	register int l, r, i, j;\
	GENERIC gtmp;\
\
	l = n >> 1;\
	r = n - 1;\
	while (1) {\
		if (l > 0)\
			gtmp = a[--l];\
		else {\
			gtmp = a[r];\
			a[r] = a[0];\
			if (--r == 0) {\
				a[0] = gtmp;\
				break;\
			}\
		}\
		i = l;\
		j = l ? l << 1 : 1;\
		while (j <= r) {\
			if (j < r && GABS(a[j]) > GABS(a[j + 1]))\
				j++;\
			if (gtmp < a[j]) {\
				a[i] = a[j];\
				j += (i = j);\
			}\
			else\
				j = r + 1;\
		}\
		a[i] = gtmp;\
	}\
}
#endif	/* 0 */
