arch/m68k/math-emu/fp_log.c

Source file repositories/reference/linux-study-clean/arch/m68k/math-emu/fp_log.c

File Facts

System
Linux kernel
Corpus path
arch/m68k/math-emu/fp_log.c
Extension
.c
Size
3816 bytes
Lines
208
Domain
Architecture Layer
Bucket
arch/m68k
Inferred role
Architecture Layer: implementation source
Status
source implementation candidate

Why This File Exists

CPU and platform-specific kernel glue: boot entry, traps, syscall entry, interrupts, page tables, context switch, and low-level barriers.

Dependency Surface

Detected Declarations

Annotated Snippet

fp_log.c: floating-point math routines for the Linux-m68k
  floating point emulator.

  Copyright (c) 1998-1999 David Huggins-Daines / Roman Zippel.

  I hereby give permission, free of charge, to copy, modify, and
  redistribute this software, in source or binary form, provided that
  the above copyright notice and the following disclaimer are included
  in all such copies.

  THIS SOFTWARE IS PROVIDED "AS IS", WITH ABSOLUTELY NO WARRANTY, REAL
  OR IMPLIED.

*/

#include "fp_arith.h"
#include "fp_emu.h"
#include "fp_log.h"

static const struct fp_ext fp_one = {
	.exp = 0x3fff,
};

struct fp_ext *fp_fsqrt(struct fp_ext *dest, struct fp_ext *src)
{
	struct fp_ext tmp, src2;
	int i, exp;

	dprint(PINSTR, "fsqrt\n");

	fp_monadic_check(dest, src);

	if (IS_ZERO(dest))
		return dest;

	if (dest->sign) {
		fp_set_nan(dest);
		return dest;
	}
	if (IS_INF(dest))
		return dest;

	/*
	 *		 sqrt(m) * 2^(p)	, if e = 2*p
	 * sqrt(m*2^e) =
	 *		 sqrt(2*m) * 2^(p)	, if e = 2*p + 1
	 *
	 * So we use the last bit of the exponent to decide whether to
	 * use the m or 2*m.
	 *
	 * Since only the fractional part of the mantissa is stored and
	 * the integer part is assumed to be one, we place a 1 or 2 into
	 * the fixed point representation.
	 */
	exp = dest->exp;
	dest->exp = 0x3FFF;
	if (!(exp & 1))		/* lowest bit of exponent is set */
		dest->exp++;
	fp_copy_ext(&src2, dest);

	/*
	 * The taylor row around a for sqrt(x) is:
	 *	sqrt(x) = sqrt(a) + 1/(2*sqrt(a))*(x-a) + R
	 * With a=1 this gives:
	 *	sqrt(x) = 1 + 1/2*(x-1)
	 *		= 1/2*(1+x)
	 */
	/* It is safe to cast away the constness, as fp_one is normalized */
	fp_fadd(dest, (struct fp_ext *)&fp_one);
	dest->exp--;		/* * 1/2 */

	/*
	 * We now apply the newton rule to the function
	 *	f(x) := x^2 - r
	 * which has a null point on x = sqrt(r).
	 *
	 * It gives:
	 *	x' := x - f(x)/f'(x)
	 *	    = x - (x^2 -r)/(2*x)
	 *	    = x - (x - r/x)/2
	 *          = (2*x - x + r/x)/2
	 *	    = (x + r/x)/2
	 */
	for (i = 0; i < 9; i++) {
		fp_copy_ext(&tmp, &src2);

		fp_fdiv(&tmp, dest);
		fp_fadd(dest, &tmp);
		dest->exp--;
	}

Annotation

Implementation Notes