90 lines
1.9 KiB
C
90 lines
1.9 KiB
C
|
/* void _SqrtDec80Bit(pDst,pSrc)
|
||
|
|
*
|
||
|
|
* ARGUMENT
|
||
|
|
* DEC *pDst;
|
||
|
|
* DEC *pSrc;
|
||
|
|
*
|
||
|
|
* DESCRIPTION
|
||
|
|
* Sets pDst = square root of pSrc;
|
||
|
|
*
|
||
|
|
* SIDE EFFECTS
|
||
|
|
* None.
|
||
|
|
*
|
||
|
|
* RETURNS
|
||
|
|
* None. Always succeeds because CALLER GUARANTEES a valid 64-bit
|
||
|
|
* value for square root operation.
|
||
|
|
*
|
||
|
|
* AUTHOR
|
||
|
|
* Jared Levy Aug 7, 1987
|
||
|
|
* Copyright (C) 1987-1990 Greenleaf Software Inc. All rights reserved.
|
||
|
|
*
|
||
|
|
* MODIFICATIONS
|
||
|
|
*
|
||
|
|
*/
|
||
|
|
|
||
|
|
#include <stdio.h>
|
||
|
|
#include "gm.h"
|
||
|
|
#include "gmsystem.h"
|
||
|
|
|
||
|
|
void _SqrtDec80Bit(pDst,pSrc)
|
||
|
|
DEC *pDst, *pSrc;
|
||
|
|
{
|
||
|
|
int m, precision;
|
||
|
|
DEC *temp, dtemp, *nsrc, dnsrc;
|
||
|
|
|
||
|
|
if(_MacIsDecZ(pSrc)) {
|
||
|
|
_MacDZero(pDst);
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
|
||
|
|
temp = &dtemp;
|
||
|
|
nsrc = &dnsrc;
|
||
|
|
|
||
|
|
/* adjust pSrc to range .01 < nsrc < 1 */
|
||
|
|
_MacDCopy(nsrc,pSrc);
|
||
|
|
m = MagnitudeOfDecimal(nsrc);
|
||
|
|
|
||
|
|
/* calculate first guess */
|
||
|
|
if ((m & 0X0001) == 0) {
|
||
|
|
m+=2;
|
||
|
|
nsrc->dc.id+= m;
|
||
|
|
if (nsrc->dc.id > GM_IMAXID) {
|
||
|
|
_DivUnsArrByPwrOf10(nsrc->dc.sl, 5, nsrc->dc.id-GM_IMAXID);
|
||
|
|
nsrc->dc.id = GM_IMAXID;
|
||
|
|
}
|
||
|
|
/* decTwoPoint59 = 2.59 */
|
||
|
|
(void) _MulDec80Bit(temp, &decTwoPoint59, nsrc);
|
||
|
|
/* decPoint0819 = .0819 */
|
||
|
|
(void) _AddDec80Bit(pDst, &decPoint0819, temp); }
|
||
|
|
else {
|
||
|
|
m++;
|
||
|
|
nsrc->dc.id+= m;
|
||
|
|
/* decPoint819 = .819 */
|
||
|
|
(void) _MulDec80Bit(temp, &decPoint819, nsrc);
|
||
|
|
/* decPoint259 = .259 */
|
||
|
|
(void) _AddDec80Bit(pDst, &decPoint259, temp); }
|
||
|
|
|
||
|
|
/* do iterative calculation */
|
||
|
|
precision = 3;
|
||
|
|
(void) _ScaleDec80Bit(pDst, pDst, 4);
|
||
|
|
while (precision < 18) {
|
||
|
|
precision = 2 * precision - 2;
|
||
|
|
(void) _DivRndDec80Bit(temp, nsrc, pDst, precision);
|
||
|
|
(void) _AddDec80Bit(pDst, pDst, temp);
|
||
|
|
_HalveUnsArr(pDst->dc.sl, 5);
|
||
|
|
}
|
||
|
|
|
||
|
|
/* last iteration */
|
||
|
|
(void) _DivDec80Bit(temp, nsrc, pDst);
|
||
|
|
(void) _AddDec80Bit(pDst, pDst, temp);
|
||
|
|
if (((pDst->dc.sl[0] & 1) == 0) || (pDst->dc.msd > 6553))
|
||
|
|
_HalveUnsArr(pDst->dc.sl, 5);
|
||
|
|
else
|
||
|
|
(void) _MulDec80Bit(pDst, pDst, &decPoint5);
|
||
|
|
|
||
|
|
/* recalibrate to correct magnitude */
|
||
|
|
pDst->dc.id -= (m/2);
|
||
|
|
|
||
|
|
}
|
||
|
|
|