Can't Copy and Paste this?

Click here for a copy-and-paste friendly version of this code!

//**************************************

//     

//INCLUDE files for :Integer square root

//     

//**************************************

//     

/* +++Date last modified: 05-Jul-1997 */

/*

** SNIPMATH.H - Header file for SNIPPETS math functions and macros

*/

#ifndef SNIPMATH__H

#define SNIPMATH__H

#include 

#include "sniptype.h"

#include "round.h"

/*

** Callable library functions begin here

*/

voidSetBCDLen(int n); /* Bcdl.C */

longBCDtoLong(char *BCDNum); /* Bcdl.C */

voidLongtoBCD(long num, char BCDNum[]);/* Bcdl.C */

double bcd_to_double(void *buf, size_t len, /* Bcdd.C */

int digits);

int double_to_bcd(double arg, char *buf, /* Bcdd.C */

size_t length, size_t digits );

DWORDncomb1 (int n, int m);/* Combin.C*/

DWORDncomb2 (int n, int m);/* Combin.C*/

voidSolveCubic(double a, double b, double c, /* Cubic.C*/

double d, int *solutions,

double *x);

DWORDdbl2ulong(double t); /* Dbl2Long.C */

longdbl2long(double t);/* Dbl2Long.C */

double dround(double x); /* Dblround.C */

/* Use #defines for Permutations and Combinations -- Factoryl.C */

#define log10P(n,r) (log10factorial(n)-log10factorial((n)-(r)))

#define log10C(n,r) (log10P((n),(r))-log10factorial(r))

double log10factorial(double N); /* Factoryl.C */

double fibo(unsigned short term);/* Fibo.C */

double frandom(int n);/* Frand.C*/

double ipow(double x, int n);/* Ipow.C */

int ispow2(int x);/* Ispow2.C*/

longdouble ldfloor(long double a);/* Ldfloor.C */

int initlogscale(long dmax, long rmax);/* Logscale.C */

longlogscale(long d); /* Logscale.C */

floatMSBINToIEEE(float f); /* Msb2Ieee.C */

floatIEEEToMSBIN(float f); /* Msb2Ieee.C */

int perm_index (char pit[], int size);/* Perm_Idx.C */

int round_div(int n, int d); /* Rnd_Div.C */

longround_ldiv(long n, long d);/* Rnd_Div.C */

double rad2deg(double rad); /* Rad2Deg.C */

double deg2rad(double deg); /* Rad2Deg.C */

#include "pi.h"

#ifndef PHI

#define PHI ((1.0+sqrt(5.0))/2.0) /* the golden number*/

#define INV_PHI (1.0/PHI) /* the golden ratio */

#endif

/*

** File: ISQRT.C

*/

    struct int_sqrt {

    unsigned sqrt,

    frac;

};

void usqrt(unsigned long x, struct int_sqrt *q);

#endif /* SNIPMATH__H */

code:

Can't Copy and Paste this?

Click here for a copy-and-paste friendly version of this code!

Terms of Agreement:   

By using this code, you agree to the following terms...   

1) You may use this code in your own programs (and may compile it into a program and distribute it in compiled format for langauges that allow it) freely and with no charge.   

2) You MAY NOT redistribute this code (for example to a web site) without written permission from the original author. Failure to do so is a violation of copyright laws.   

3) You may link to this code from another website, but ONLY if it is not wrapped in a frame. 

4) You will abide by any additional copyright restrictions which the author may have placed in the code or code's description.  

/* +++Date last modified: 05-Jul-1997 */

#include 

#include "snipmath.h"

#define BITSPERLONG 32

#define TOP2BITS(x) ((x & (3L << (BITSPERLONG-2))) >> (BITSPERLONG-2))

/* usqrt:

ENTRY x: unsigned long

EXIT returns floor(sqrt(x) * pow(2, BITSPERLONG/2))

Since the square root never uses more than half the bits

of the input, we use the other half of the bits to contain

extra bits of precision after the binary point.

EXAMPLE

suppose BITSPERLONG = 32

thenusqrt(144) = 786432 = 12 * 65536

usqrt(32) = 370727 = 5.66 * 65536

NOTES

(1) change BITSPERLONG to BITSPERLONG/2 if you do not want

the answer scaled. Indeed, if you want n bits of

precision after the binary point, use BITSPERLONG/2+n.

The code assumes that BITSPERLONG is even.

(2) this is really better off being written in assembly.

The line marked below is really a "arithmetic shift left"

on the double-long value with r in the upper half

and x in the lower half. this operation is typically

expressible in only one or two assembly instructions.

(3) Unrolling this loop is probably not a bad idea.

ALGORITHM

The calculations are the base-two analogue of the square

root algorithm we all learned in grammar school. Since we're

in base 2, there is only one nontrivial trial multiplier.

Notice that absolutely no multiplications or divisions are performed.

this means it'll be fast on a wide range of processors.

*/

void usqrt(unsigned long x, struct int_sqrt *q)

    {

    unsigned long a = 0L;/* accumulator */

    unsigned long r = 0L;/* remainder*/

    unsigned long e = 0L;/* trial product*/

    int i;

    for (i = 0; i < BITSPERLONG; i++)/* NOTE 1 */

        {

        r = (r << 2) + TOP2BITS(x); x <<= 2; /* NOTE 2 */

        a <<= 1;

        e = (a << 1) + 1;

        if (r >= e)

            {

            r -= e;

            a++;

        }

    }

    memcpy(q, &a, sizeof(long));

}

#ifdef TEST

#include 

#include 

main(void)

    {

    int i;

    unsigned long l = 0x3fed0169L;

    struct int_sqrt q;

    for (i = 0; i < 101; ++i)

        {

        usqrt(i, &q);

        printf("sqrt(%3d) = %2d, remainder = %2d\n",

        i, q.sqrt, q.frac);

    }

    usqrt(l, &q);

    printf("\nsqrt(%lX) = %X, remainder = %X\n", l, q.sqrt, q.frac);

    return 0;

}

#endif /* TEST */