campo-sirio/gfm/dsqx.c
alex ba237a9d91 Patch level : no patch
Files correlati     :
Ricompilazione Demo : [ ]
Commento            :
Aggiunti i sorgenti per Greenleaf Math Library (gfm.dll)


git-svn-id: svn://10.65.10.50/trunk@10079 c028cbd2-c16b-5b4b-a496-9718f37d4682
2002-02-26 12:19:02 +00:00

90 lines
1.9 KiB
C
Executable File

/* 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);
}