From ffca5a5837d1d2c579d5c7a109fe631e284451da Mon Sep 17 00:00:00 2001 From: Bruno Haible Date: Tue, 9 Mar 2004 12:52:13 +0000 Subject: [PATCH] Fix a bug occurring with extremely high exponents. --- ChangeLog | 2 ++ src/float/lfloat/algebraic/cl_LF_sqrt.cc | 29 ++++++++++++------------ 2 files changed, 17 insertions(+), 14 deletions(-) diff --git a/ChangeLog b/ChangeLog index e986632..67dc8ba 100644 --- a/ChangeLog +++ b/ChangeLog @@ -2,6 +2,8 @@ * src/float/lfloat/elem/cl_LF_mul.cc (operator*): Fix the second underflow condition. + * src/float/lfloat/algebraic/cl_LF_sqrt.cc (sqrt): Fix a bug with large + uexp whereby SQRT of MOST-POSITIVE-LONG-FLOAT was less than 1. 2004-03-04 Richard B. Kreckel diff --git a/src/float/lfloat/algebraic/cl_LF_sqrt.cc b/src/float/lfloat/algebraic/cl_LF_sqrt.cc index c551ca4..613e4fb 100644 --- a/src/float/lfloat/algebraic/cl_LF_sqrt.cc +++ b/src/float/lfloat/algebraic/cl_LF_sqrt.cc @@ -28,7 +28,7 @@ const cl_LF sqrt (const cl_LF& x) // Bei ungeradem e schiebe dies (oder nur die ersten n+1 Digits davon) // um 1 Bit nach rechts. // Bilde daraus die Ganzzahl-Wurzel, eine n+1-Digit-Zahl mit einer -// führenden 1. +// fhrenden 1. // Runde das letzte Digit weg: // Bit 15 = 0 -> abrunden, // Bit 15 = 1, Rest =0 und Wurzel exakt -> round-to-even, @@ -42,25 +42,26 @@ const cl_LF sqrt (const cl_LF& x) CL_ALLOCA_STACK; var uintD* r_MSDptr; var uintD* r_LSDptr; - var uintL r_len = 2*(uintL)len+2; // Länge des Radikanden + var uintL r_len = 2*(uintL)len+2; // L�ge des Radikanden num_stack_alloc(r_len, r_MSDptr=,r_LSDptr=); - uexp = uexp - LF_exp_mid + 1; - if (uexp & bit(0)) + if ((uexp & bit(0)) == (LF_exp_mid & bit(0))) // Exponent gerade {var uintD* ptr = copy_loop_msp(arrayMSDptr(TheLfloat(x)->data,len),r_MSDptr,len); // n Digits kopieren - clear_loop_msp(ptr,len+2); // n+2 Nulldigits anhängen + clear_loop_msp(ptr,len+2); // n+2 Nulldigits anh�gen } else // Exponent ungerade {var uintD carry_rechts = // n Digits kopieren und um 1 Bit rechts shiften shiftrightcopy_loop_msp(arrayMSDptr(TheLfloat(x)->data,len),r_MSDptr,len,1,0); var uintD* ptr = r_MSDptr mspop len; - msprefnext(ptr) = carry_rechts; // Übertrag und - clear_loop_msp(ptr,len+1); // n+1 Nulldigits anhängen + msprefnext(ptr) = carry_rechts; // �ertrag und + clear_loop_msp(ptr,len+1); // n+1 Nulldigits anh�gen } - uexp = (sintL)((sintL)uexp >> 1); // Exponent halbieren - uexp = uexp + LF_exp_mid; + // Compute ((uexp - LF_exp_mid + 1) >> 1) + LF_exp_mid without risking + // uintL overflow. + uexp = ((uexp - ((LF_exp_mid - 1) & 1)) >> 1) - ((LF_exp_mid - 1) >> 1) + + LF_exp_mid; // Ergebnis allozieren: var Lfloat y = allocate_lfloat(len,uexp,0); var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len); @@ -80,7 +81,7 @@ const cl_LF sqrt (const cl_LF& x) // Ablegen und runden: copy_loop_msp(p_MSDptr mspop 1,y_mantMSDptr,len); // NUDS nach y kopieren if (mspref(p_MSDptr,0) == 0) - { if ( ((sintD)mspref(p_MSDptr,len+1) >= 0) // nächstes Bit =0 -> abrunden + { if ( ((sintD)mspref(p_MSDptr,len+1) >= 0) // n�hstes Bit =0 -> abrunden || ( ((mspref(p_MSDptr,len+1) & ((uintD)bit(intDsize-1)-1)) ==0) // =1 und weitere Bits >0 -> aufrunden && !test_loop_msp(p_MSDptr mspop (len+2),len+1) // round-to-even (etwas witzlos, da eh alles ungenau ist) @@ -91,13 +92,13 @@ const cl_LF sqrt (const cl_LF& x) else // aufrunden { if ( inc_loop_lsp(y_mantMSDptr mspop len,len) ) - // Übertrag durchs Aufrunden + // �ertrag durchs Aufrunden { mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0 (TheLfloat(y)->expo)++; // Exponenten incrementieren } } } else - // Übertrag durch Rundungsfehler + // �ertrag durch Rundungsfehler { if (test_loop_msp(y_mantMSDptr,len)) cl_abort(); mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0 (TheLfloat(y)->expo)++; // Exponenten incrementieren @@ -111,7 +112,7 @@ const cl_LF sqrt (const cl_LF& x) // w ist die Ganzzahl-Wurzel, eine n+1-Digit-Zahl. copy_loop_msp(w.MSDptr,y_mantMSDptr,len); // NUDS nach y kopieren // Runden: - if ( ((sintD)lspref(w.LSDptr,0) >= 0) // nächstes Bit =0 -> abrunden + if ( ((sintD)lspref(w.LSDptr,0) >= 0) // n�hstes Bit =0 -> abrunden || ( ((lspref(w.LSDptr,0) & ((uintD)bit(intDsize-1)-1)) ==0) // =1 und weitere Bits >0 oder Rest >0 -> aufrunden && exactp // round-to-even @@ -122,7 +123,7 @@ const cl_LF sqrt (const cl_LF& x) else // aufrunden { if ( inc_loop_lsp(y_mantMSDptr mspop len,len) ) - // Übertrag durchs Aufrunden + // �ertrag durchs Aufrunden { mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0 (TheLfloat(y)->expo)++; // Exponenten incrementieren } }