Morgan - Numerical Methods (523161), страница 37
Текст из файла (страница 37)
This routine takes a floating-pointnumber as an argument and replaces the exponent, that is, it raises the mantissa to anew power. It computes input_float * 2new-exponent Its operation is simple:Idxp: Algorithm1. Test the input floating point argument for zero.If it's zero, exit with zero as the result through step 6.2. Save the sign and replace the current exponent with 126D.3.
Add the new exponent and test for overflow.If there is an overflow, exit through the overflow error exit, step 7.4. Shift the sign back into place along with the exponent.5. Write the new float to the output and leave.6. Zero error exit; write zero out.7. Overflow-error exit; write infinite out.Idxp: Listing; *****;Ldxp is similar to ldexp in C, it is used for math functions.;Takes from the stack passed with it an input float (extended) and returns a;pointer to a value to the power of two.1dxpprocuses di, float:qword, power:word, exp:bytemovmovax, word ptr float [4]dx, word ptr float [2]subcx, cx;get upper word of float;extended bits are not;checked261NUMERICAL METHODSororjeshlrclmovaddjccx, axcx, dxreturn_zeroax, 1cl, 1ah, 7ehah, byte ptr expld_overflowshrrcrmovcl, 1word ptr ax, 1word ptr float[4], axldxp_exit:movmovleamovswrepretld_overflow:movsubmovmovjwreturn_zero:submovmovstoswrepjmpldxp endpcx, 4di, word ptr powersi, word ptr float;save the sign;add new exponent;return the sign;position exponent;write the result outword ptr float[4], 7f80hax, axword ptr float[2], axword ptr float[0], axldxp_exitax, axdi, word ptr powercx, 4ldxp_exitThe next three functions are all related.
The first, flr, implements the C functionfloor() and returns the largest floating-point mathematical integer not greater thanthe input.262THE ELEMENTARY FUNCTIONSflr: Algorithm1. Get the float, extract the exponent, and subtract 126D.If there is an underflow, the number must be less than .5; exitthrough step 5.2. Subtract the reduced exponent from 40D. This the mantissa portionplus extendedprecision.If the result is less than the reduced exponent, we already have thefloor (it's all integer); exit through step 3.Otherwise, save the number of shifts in shift.Shift the float right, shifting off the fraction bits, until theexponent is exhausted.
What remains are integer bits.3. Get the exponent back from shift.Shift the float back into its proper position, this time without thefractionbits. This is the floor of the argument.4. Leave, writing the result to the output.5. Exit with a result of zero.flr: Listing; ******;floor greatest integer less than or equal to x;single precisionflrprocuses bx dx di si, fp:qword, rptr:wordlocalshift :bytemovmovdi, word ptr rptrbx, word ptr fp[0]movmovmovandax,dx,cx,cx,word ptr fp[2]word ptr fp[4]dx7f80hshlmovsubsubcx,cl,ch,cl,1chch7eh;get float with extended;precision;get rid of sign and;mantissa portion;subtract bias (-1)from;exponent263NUMERICAL METHODSjbemovsubleave_with_zeroch, 40ch, cljbalready_floormovmovsubbyte ptr shift, chcl, chch, chshrdx, 1rcrrcrloopax, 1bx, 1fix;is it greater than the;mantissa portion?;there is no fractional;partfix:movre_position:shlrclrclloopbx, 1ax, 1dx, 1repositionalready-floor:movmovmovsubmovword ptr [di][4],word ptr [di][2],word ptr [di][0],ax, axwordptr [di][6],flr_exit:retleave_with_one:leamovmovmovswrepjmp;position as fixed pointcl, byte ptr shiftsi, word ptr onedi, word ptr rptrcx, 4;realign floatdxaxbx;write to outputax;floating-point oneflr_exitleave_with_zero:subax, axmovcx, 4264;shift the number the;number of times indicated;in the exponent-floating-point zeroTHE ELEMENTARY FUNCTIONSrepflrmovstoswjmpdi, word ptr rptrshort flr_exitendpThe complement to flr is flceil.
This routine is similar to the C function ceil() thatreturns the smallest floating-point mathematical integer not less than the inputargument.flceil: Algorithm1. Get the float and check for zero.If the input argument is zero, exit through step 6.If the input is not zero, continue.Extract the exponent and subtract 126D. If there is anunderflow,then the number must be less than .5; exit through step 5.2. Subtract the reduced exponent from 40D. This is the mantissa portionplus extended precision.If the result is less than the reduced exponent, we already have theceiling (it's all integer); exit through step 3.Otherwise, save the number of shifts in shift.Shift the float right, shifting the fractionbits into the MSW of thefloating-point data type until the exponent is exhausted. Whatremains are integer bits.Test the MSW of the floating-point data type.If it's zero, go to step 3.If it's anything else, round the integer portion up andcontinue with step 3.3.
Get the exponent back from shift.Shift the float back into its proper position, this time without thefraction bits. This is the floor of the argument.4. Leave, writing the result to the output.5. Exit with a result of one.6. Exit with a result of zero.265NUMERICAL METHODSflceil: Listing; *****;flceil least integer greater than or equal to x;single precision;flceil procuses bx dx di si, fp:qword, rptr:wordlocalshift:bytemovmovdi, word ptr rptrbx, word ptr fp[0]movmovsuborororjemovandax, word ptr fp[2]dx, word ptr fp[4]cx, cxcx, bxcx, axcx, dxleave_with_zerocx, dxcx, 7f80hshlmovsubsubcx,cl,ch,cl,jbmovsubleave-with-onech, 40ch, cljbalready_ceilmovmovsubbyte ptr shift, chcl, chch, chshrdx, 1rcrrcrrcrax, 1bx, 1word ptr [di] [6], 11chch7eh;get float with extended;precision;this is a zero;get rid of sign and;mantissa portion;subtract bias (-1) from;exponent;is it greater than the;mantissa portion?;there is no fractional;partfix:266;shift the number the;number of times indicated;in the exponent;put guard digits in MSW of;data typeTHE ELEMENTARY FUNCTIONSloopfix;position as fixed pointcmpjeaddadcadcword ptr [di][6],0hnot_quite_enoughbx, 1ax, 0dx, 0;roundupnot_quite_enough:movcl, byte ptr shiftreposition:shlbx, 1rclax, 1dx, 1rclloopre_positionalready_ceil:movmovmovsubmovwordwordwordax,wordptrptrptraxptr[di][4], dx[di][2], ax[di][0], bx;realign float;write to output[di][6], axceil-exit:retleave-with-one:leamovmovmovswrepjmpleave_with_zero:submovmovstoswrepjmpsi, word ptr onedi, word ptr rptrcx, 4;a floating-point oneceil-exitax, axcx, 4di, word ptr rptr;a floating-point zeroshort ceil_exitflceil endp267NUMERICAL METHODSFinally, intrnd rounds the input argument to its closest integer.
As used by Codyand Waite,9 this function returns an integer representing the mathematical integerclosest to the input float. It employs no rounding logic; if the mantissa portion of theinput float is greater than .5, the next higher whole integer is returned. In thisimplementation, however, it returns a floating-point number representing the mathematical integer closest to the input. It was written that way to accommodate otherroutines in the floating-point package.intrnd: Algorithm1.
Subtract the value returned by flr from the input and take theabsolute value of the result.2. Compare the result with .5.If it's greater, get the flceil of the input.If it's equal to or less than, go to step 3.3. Write the result to the output and return.intrnd: Listing; ******;intrnd is useful for the transcendental functions;it rounds to the nearest integer according to the logic;intrnd(x) =if((x-floor(x)) <.5) floor(x);else ceil(x);;intrnd procuses bx dx di si, fp:qword, rptr:wordlocal temp0:qword, temp1:qword,pushfcldsubmovleareprep268stoswmovleastoswmovax, axcx, 4di, word ptr temp0;prepare intermediate;registerscx, 4di, word ptr temp1di, word ptr rptr;clear the outputTHE ELEMENTARY FUNCTIONSrepmovstoswcx, 4invokeinvokeandinvokeflr, fp, addr temp0flsub, fp, tempo, addr templword ptr temp1[4], 7fffhflcomp, temp1, one-halfcmpjneax, 1intrnd_exit;greater than .5?flceil, p, addr temp0;get the ceiling of the;inputdo_ceil:invokeintrnd_exit:movmovmovmovmovpopfretintmd endpax, worddx, worddi, wordword ptrword ptr;cheap fabsptr temp0[2]ptr temp0[4]ptr rptr[di][2], ax[di][4], dxDealing with real numbers in a finite machine means we must deal withlimitations.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
