fixed64.c

wm PNE 3.3 source code, running at more than vxworks6.x version.

/* $Header: /usr/cvsroot/target/src/wrn/wm/common/lib/fixed64.c,v 1.3 2003/01/16 18:18:55 josh Exp $ */

/*
 * Copyright (C) 1999-2005 Wind River Systems, Inc.
 * All rights reserved.  Provided under license only.
 * Distribution or other use of this software is only
 * permitted pursuant to the terms of a license agreement
 * from Wind River Systems (and is otherwise prohibited).
 * Refer to that license agreement for terms of use.
 */


/****************************************************************************
 *  Copyright 1999 Integrated Systems, Inc.
 *  All rights reserved.
 ****************************************************************************/

/*
 * Fixed-point 64-bit math routines, intended primarily for time conversion
 * between Attache milliseconds and pSOSystem "ticks".  These routines
 * are in the common library because they may be useful for other
 * applications, but these functions are not a full-fleged 64-bit math
 * library.
 */

/*
 * $Log: fixed64.c,v $
 * Revision 1.3  2003/01/16 18:18:55  josh
 * directory structure shifting
 *
 * Revision 1.2  2001/11/06 22:15:50  tneale
 * Fixed for newest file layout
 *
 * Revision 1.1.1.1  2001/11/05 17:48:39  tneale
 * Tornado shuffle
 *
 * Revision 1.3  2001/01/19 22:21:30  paul
 * Update copyright.
 *
 * Revision 1.2  2000/03/17 00:16:53  meister
 * Update copyright message
 *
 * Revision 1.1  1999/03/01 20:08:56  sra
 * Add fixed-point 64-bit math routines; these are primarily intended for
 * converting between Attache and pSOSystem clock values, but the
 * functions themselves aren't specific to that.
 *
 */

/* [clearcase]
modification history
-------------------
01b,20apr05,job  update copyright notices
01a,11dec03,job  fix copyright statements
*/


#include <wrn/wm/common/config.h>
#include <wrn/wm/common/bug.h>
#include <wrn/wm/common/fixed64.h>
#include <wrn/wm/common/glue.h>

/*
 * Multiply two-64 bit unsigned fixed-point numbers (32 integer bits, 32
 * fractional bits) and return the integer part of the result, mod 2**32.
 *
 * The algorithm here is from Knuth, volume 2, 2nd edition, page 278.
 * We're trying to multiply two 2N-bit numbers u and v:
 *
 *	u = 2**N * u1 + u0
 *	v = 2**N * v1 + v0
 *
 * 	u * v = (2**(2*N) + 2**N    ) * (u1 * v1)
 *            + (           2**N    ) * ((u1 - u0) * (v0 - v1))
 *            + (           2**N + 1) * (u0 * v0)
 *
 * We use this algorithm twice, once with N=32, once with N=16.
 *
 * The final calculation yields a 128 bit result, but we can drop some
 * terms because we know they're just going to be outside the range of
 * the final result anyway.  All we really want are the bits of the
 * result from 2**32 through 2**63, inclusive.  So we can throw away
 * the 2**64 term, and we can use native arithmetic for all the 2**32
 * terms because we don't care about any high bits in those terms.
 * However, we have to calculate all 64 bits of the (u0 * v0) product
 * because of possible carries from the final term (1 * u0 * v0),
 * so we use this algorithm again, this time with N=16 (that is, using
 * the low 16 bits of 32-bit integers, so that we don't lose any bits).
 *
 * Most of the work here is just calculating the (x0 * y0) product,
 * and most of that is just the grunt work of summing up the 16-bit
 * partial products.  The only tricky bit is handling the sign of the
 * "b" term; the multiplication itself comes out correctly, but the
 * sign bit doesn't propegate correctly through the summation without
 * some help.
 */

bits32_t fixed64_mul(bits32_t x1,	/* high bits of X */
		     bits32_t x0,	/* low  bits of X */
		     bits32_t y1,	/* high bits of Y */
		     bits32_t y0)	/* low  bits of Y */
{
  bits32_t u1 = x0 >> 16, u0 = x0 & 0xFFFF;
  bits32_t v1 = y0 >> 16, v0 = y0 & 0xFFFF;
  bits32_t a  = u1 * v1;
  bits32_t b  = (u1 - u0) * (v0 - v1);
  bits32_t c  = u0 * v0;
  bits32_t p1 = a, p0 = c, w;

  w = p0, p0 += (a << 16), p1 += (a >> 16) + (p0 < w);
  w = p0, p0 += (c << 16), p1 += (c >> 16) + (p0 < w);

  w = p0; 
  if ((u1 < u0) != (v0 < v1))
    p0 -= (-b << 16), p1 -= (-b >> 16) + (p0 > w);
  else
    p0 += ( b << 16), p1 += ( b >> 16) + (p0 < w);

  return x1 * y1 + (x1 - x0) * (y0 - y1) + p0 + p1;
}

/*
 * Divide a 32-bit unsigned integer by another 32-bit unsigned integer
 * to form a 64-bit unsigned fixed point number (32 integer bits, 32
 * fractional bits).
 *
 * We only expect to use this function at boot time, so it's optimized
 * for small code space and readability rather than speed of the
 * algorithm.  Intended use is to compute things like the conversion
 * ratios between milliseconds (Attache) and "ticks" (pSOSystem).
 */

void fixed64_div(bits32_t numerator,
		 bits32_t denominator,
		 bits32_t *integer,
		 bits32_t *fractional)
{
  bits32_t dividend;
  bits32_t remainder;
  bits32_t overflow;

  BUG_ASSERT(denominator != 0 && integer != 0 && fractional != 0);

  *integer    = numerator / denominator;
  remainder   = numerator % denominator;
  *fractional = 0;

  for (dividend = 0x80000000;
       dividend && remainder;
       dividend >>= 1) {
    overflow = remainder & 0x80000000;
    remainder <<= 1;
    if (overflow || remainder >= denominator) {
      remainder -= denominator;
      *fractional += dividend;
    }
  }
}