(* Release 3.10 *) (*-------------------------------------------------------------------------* * * * COREEMUL.A - 8087 emulator * * * * COPYRIGHT (C) 1989..1992 Clarion Software Corporation. * * All Rights Reserved * * * *--------------------------------------------------------------------------*) include "corelib.inc" module CoreEmul (********************************************************************) (* 8087 emulation support - part common to real and protected modes *) (********************************************************************) false = 0 (* Beware ! For incoming parameters, non-false = true.*) true = 1 (* For results, we generate the proper numbers.*) lesser = -1 (* Incoming, lesser < 0*) equal = 0 greater = 1 (* Incoming, greater > 0*) by0 = 0 by1 = 1 w0 = 0 w1 = 2 w2 = 4 w3 = 6 dd0 = 0 dd1 = 4 nearParam = 4 emInt = 34H (* allocated to 8087 by Microsoft Corp.*) (* The iAPX87 chips usually report exceptions when the next numeric *) (* instruction is attempted (since processing runs in parallel *) (* with the CPU). The emulator will normally do the same, but *) (* if the following symbol is defined some optional assembly is *) (* enabled which allows exceptions to be emulated in the failing *) (* instruction. The main benefit is that FWAIT may validly be *) (* converted to no-op, gaining speed. *) instantException = true (* report excep immediately *) (* error types *) fe_OK = 0 fe_Divide_By_Zero = 04H fe_Imaginary = 01H fe_Operand_Too_Small = 02H fe_Operand_Too_Big = 01H fe_Overflow = 08H fe_Underflow = 10H fe_Bad_Parameter = 01H fe_Invalid_Operation = 01H fe_No_Space = 01H fe_Precision = 20H (* Using a 15-bit signed exponent is enough to emulate the iNDP-87 *) (* and simplifies avoids overflow when combining exponents. *) zero_exponent = 0C001H infinite_exponent = 4001H positive_sign = 0 negative_sign = 1 (* 8087/287 status word values after comparisons *) e287_incomparable = 4500H e287_lesser = 0100H e287_equal = 4000H e287_greater = 0000H (* 8087/287 status word after XAM *) e287_plusInfinity = 0500H e287_plusNormal = 0400H e287_plusZero = 4000H e287_minusNormal = 0600H e287_empty = 4100H (* control word rounding control field *) e287_roundingControl = 0CH e287_round = 00H e287_floor = 04H e287_ceiling = 08H e287_chop = 0CH (* IEEE exponent bias constants, relative to temporary format *) iees_exponent_bias = 7EH ieel_exponent_bias = 3FEH temp_exponent_bias = 3FFEH (* The temporary format here is not identical to iNDP-87 format. The *) (* exponent is not biased and the sign is held in a separate *) (* byte. This format is faster for emulation software to use. *) emu_fraction = 0 emu_exponent = 8 emu_sign = 10 emu_spare = 11 emu_temp_size = 12 (*INCLUDE E086ENTR*) (* The emulation does not extend to some more exotic possibilities for *) (* using the iNDP registers. When stack overflow occurs, registers do *) (* not wrap around. When registers below stack bottom are addressed, the *) (* register actually used is undefined. In both cases, the emulator *) (* generates an "invalid operation" exception. If the exception is not *) (* masked and exception interrupts are properly supported, the program *) (* should behave as if the real chip was present. If exceptions are not *) (* processed, then the result will be unpredictable. *) fregs = 200 emws_limitSP = 200 (* org 8*emu_temp_size *) emws_initialSP = 296 (* org 9*emu_temp_size *) emws_saveVector = 404 (* org 4 *) emws_not_used = 408 (* org 6 *) emws_status = 414 (* org 2, aligned. Result of comparisons *) emws_control = 416 (* org 2 processing options and exceptions. *) emws_instrnPtrxx = 418 (* org 4 used with error recovery *) emws_dataPtrxx = 422 (* org 4 --------- " ------------ *) emws_instructionxx = 426 (* org 2 bytes swapped, used for error recovery *) emws_tos = 428 (* org 2 current level of e87 register stack *) emws_maxStack = 430 (* org 2 *) invLn2 = 432 (* org emu_temp_size, dw 0F0BBH, 05C17H, 03B29H, 0B8AAH, 1; db positive_sign, 0*) emws_resume_ip = 444 (* org 2 *) emws_resume_bp = 446 (* org 2 *) emuRecSize = 448 (* offsets from bp of saved registers *) oldFlags = 22 oldCS = 20 oldIP = 18 oldAX = 16 oldBX = 14 oldCX = 12 oldDX = 10 oldSI = 8 oldDI = 6 oldDS = 4 oldES = 2 oldBP = 0 segment EMU_TEXT (EMU_CODE,28H) (*%T _mthread *) extrn __FloatSave extrn __FloatRestore (*%E *) public __FloatInitStack : mov word ss:[InMathLib],0 (*%T RegParam *) push cx (* initialise fp registers except for stack top *) mov cx,7 push_zero: fldz st(0) loop push_zero pop cx (*%E *) (*%T _mthread *) push ds push es push ss pop ds push ss pop es call far __FloatSave call far __FloatRestore pop es pop ds (*%E *) ret far 0 segment EMU_TEXT(EMU_CODE,28H) public __EmulInitInstance: (* ax = resume address, bx = stack segment *) push cx push di push ds push es mov ds, bx mov es, bx mov [emws_resume_ip],ax mov di, emws_limitSP mov cx, emws_resume_ip - emws_limitSP sub ax, ax rep; stosb mov word [emws_status],0 mov word [emws_control],33FH + 0C00H and word [emws_control],~001DH (* enable invalid op/underflow/overflow/divide by zero exception *) mov word [emws_tos],(*offset*) emws_initialSP mov word [emws_maxStack],(*offset*) emws_initialSP mov word [emws_initialSP],0 mov word [emws_initialSP][2],0 mov word [emws_initialSP][4],0 mov word [emws_initialSP][6],0 mov word [emws_initialSP][8],infinite_exponent mov word [invLn2][0],0F0BBH mov word [invLn2][2],05C17H mov word [invLn2][4],03B29H mov word [invLn2][6],0B8AAH mov word [invLn2][8],00001H (*mov word [invLn2][10],00000H*) pop es pop ds pop di pop cx ret 0 other_jumps : dw e287_tos_Function dw e287_FPU_Control dw e287_Reserved_1 dw e287_Reserved_2 arith_push : dw iees_Push, emu_Float32, ieel_Push, emu_Float arith_actCase: dw arith_Add, arith_Multiply dw arith_Compare, arith_CompareP dw arith_Subtract, arith_SubtractR dw arith_Divide, arith_DivideR en_regCaseSwitch : dw en_regCase0 dw en_regCase1 dw en_regCase2 dw en_regCase3 dw en_regCase4 dw en_regCase5 dw en_regCase6 dw en_regCase7 public emulate8087: (* AL contains 0..7, the op code number, and AH contains mod/r/m DL indicates if there was a segment override ES contains the segment of the segment override (if any) *) mov ss:[emws_resume_bp], bp mov sp, ss:[emws_tos] mov bx, 0C007H and bl, ah (* bl = regMode *) and bh, ah (* bh = disp *) xchg cx, ax (* ch = mod/R/M, cl = op *) cmp bh, 0C0H jae en_NotMemoryOperand en_memoryOperand: sub ax, ax shl bl, 1 cmp bh, 40H mov bh, 0 ja en_dispCase_80 je en_dispCase_40 en_dispCase_00: cmp bl, 6*2 jne en_dispCaseEnd lodsw mov di, [bp][oldDS] jmp en_checkSegFix en_dispCase_40: lodsb cbw jmp cs:[en_regCaseSwitch][bx] en_dispCase_80: lodsw en_dispCaseEnd: jmp cs:[en_regCaseSwitch][bx] en_regCase0: add ax, [bp][oldBX] add ax, [bp][oldSI] mov di, [bp][oldDS] jmp en_regCaseEnd en_regCase1: add ax, [bp][oldBX] add ax, [bp][oldDI] mov di, [bp][oldDS] jmp en_regCaseEnd en_regCase2: add ax, [bp][oldBP] add ax, [bp][oldSI] mov di, ss jmp en_regCaseEnd en_regCase3: add ax, [bp][oldBP] add ax, [bp][oldDI] mov di, ss jmp en_regCaseEnd en_regCase4: add ax, [bp][oldSI] mov di, [bp][oldDS] jmp en_regCaseEnd en_regCase5: add ax, [bp][oldDI] mov di, [bp][oldDS] jmp en_regCaseEnd en_regCase6: add ax, [bp][oldBP] mov di, ss jmp en_regCaseEnd en_regCase7: add ax, [bp][oldBX] mov di, [bp][oldDS] en_regCaseEnd: en_checkSegFix: mov [bp][oldIP], si (* CSIP correct for next instruction *) xchg si, ax (* operandPtr := es:si *) cmp dl, true (* was segFix present ? *) je en_wasFixed mov es, di en_wasFixed: push ss pop ds (*mov [emws_dataPtr][0], si*) (*mov [emws_dataPtr][2], es*) jmp en_gotoLong en_NotMemoryOperand: mov [bp][oldIP], si (* CSIP correct for next instruction *) mov ax, ss mov ds, ax mov es, ax (* ss=cs=ds=es *) en_gotoLong: test cl, 1 jz e287_Arith cmp ch, 3 * 40H jb e287_Load_Store test ch, 20H jz e287_Register_Move other_case: mov bx, 6 and bl, cl jmp cs:[other_jumps][bx] (*public*) e287_Register_Move : (* which register ?*) mov al, emu_temp_size mul bl (* register number*) add ax, sp xchg si, ax (* operandPtr := es : si*) mov bx, 18H and bl, ch and cl, 06H or bl, cl cld mov cx, emu_temp_size/2 jmp near cs:[remo_case][bx] (*public*) e287_Arith : mov bp, cx (* save a copy of cx in BP*) cmp ch, 0C0H jnb arith_reg sub sp, emu_temp_size mov di, sp mov bx, 6 and bl, cl (* switch according to TYP field*) call cs:[arith_push][bx] mov bx, sp lea ax, [bx][emu_temp_size] mov dx, ax mov si, emu_temp_size (* si = popAfter*) jmp arith_action arith_reg: mov al, emu_temp_size mul bl (* register number*) add ax, sp xchg bx, ax (* op2 := ^register [bl];*) mov ax, sp test cl, 4 (* test .Reverse bit*) mov dx, ax jz arith_notReversed mov dx, bx arith_notReversed: arith_testPop: sub si, si (* si = popAfter*) test cl, 2 (* isolate the .Pop bit*) jz arith_action mov si, emu_temp_size arith_action: push ss pop es (* es = ds = ss*) and bp, 3800H mov cl, 6 rol bp, cl (* BP = word index [opb]*) mov cx, (*offset*) arith_end (* optimise return from cases*) jmp near cs:[arith_actCase][bp] (*public*) e287_Load_Store : (* these two statements are an optimised preparation for the case switches*) mov bp, 6 and bp, cx (* select CL.typ (same as CL.xtp)*) (* Note BP is conserved by subroutine calls*) test ch, 20H jnz ldst_extended test ch, 10H jnz ldst_store test ch, 08H jnz ldst_pushThenPop sub sp, emu_temp_size mov di, sp call cs:[ldst_loadCase][bp] (* jmp ldst_end*) ldst_pushThenPop: (*{ reserved, logically No-Op's }*) jmp ldst_end ldst_store: mov di, sp xchg si, di push cx call cs:[ldst_storeCase][bp] pop cx test ch, 08H (* pop required ?*) jz ldst_end add sp, emu_temp_size jmp ldst_end ldst_extended: mov ax, 08H and al, ch (* select CH.x*) or bp, ax (* BP = 2 * x|xtp*) test ch, 10H jnz ldst_xStore call cs:[ldst_xPushCase][bp] jmp ldst_end ldst_xStore: xchg si, di call cs:[ldst_xStoreCase][bp] ldst_end: jmp near e287_Exit (*public*) e287_Exit : mov ss:[emws_tos], sp (* record the e87 registers' new stack top *) ExceptionExit: mov bp, ss:[emws_resume_bp] jmp ss:[emws_resume_ip] (*INCLUDE E287ADSU.ASM*) (* algorithm notes*) (* subtraction is done as (a - b) = (a + (-b))*) (* align the lesser operand rightwards to match the greater.*) (* add them together, realign fraction rounded if overflow.*) (* The arithmetic is done using 5 bytes (dx:ax:CH), with the least significant*) (* byte (CH) holding the right-shifted part of the lesser of the two numbers.*) (* Actually only 2 bits are required, but use of a whole byte is easier.*) (* This 5th byte comes into play if the numbers are nearly equal. At the end*) (* of calculation the 5th byte is rounded off and a 4-byte result is*) (* delivered.*) x_p = 4+nearParam y_p = 2+nearParam z_p = 0+nearParam (* z := x +- y*) paramsize = 6 align = -1 (* used for shift counts*) excess = -2 (* bits 64..71 of calculation*) zExp = -4 (* effective exponent of z*) zSign = -5 (* effective sign of z*) op = -6 (* are we to add or to subtract*) (*public*) emu_Subtract : (*NEAR*) mov cl, negative_sign jmp @Shared_Add_Sub (*public*) emu_Add : (*NEAR*) mov cl, positive_sign @Shared_Add_Sub : (*NEAR*) push bp; mov bp,sp; lea sp,[bp][op]; push si; push di mov si, [bp][x_p] mov di, [bp][y_p] mov al, cl xor al, [di][emu_sign] xor al, [si][emu_sign] mov [bp][op], al (* op := negative if the operands are to*) (* be subtracted*) (* op := positive if same sign, addition.*) mov ax, [si][emu_exponent] mov bx, [di][emu_exponent] cmp ax, bx jge @share_x_ge_y @share_y_gt_x: xor cl, [di][emu_sign] xchg ax, bx xchg si, di jmp @share_limits @share_x_ge_y: mov cl, [si][emu_sign] @share_limits: (* si points to the larger, di to the smaller (or near-equal)*) (* ax exponent bx exponent*) mov [bp][zSign], cl (* larger's sign dominates result*) mov cx, [si][emu_exponent] (* the larger's exponent*) mov [bp][zExp], cx (* will dominate result*) cmp ax, infinite_exponent jge @larger_is_result cmp bx, zero_exponent jle @larger_is_result @share_align: (* now align the smaller number to match the larger.*) sub ax, bx cmp ax, 64 jle @share_normal @larger_is_result: mov ax, [si][emu_fraction][w0] mov bx, [si][emu_fraction][w1] mov cx, [si][emu_fraction][w2] mov dx, [si][emu_fraction][w3] jmp near @shared_end @share_normal: mov [bp][align], al mov ax, [di][emu_fraction][w0] mov bx, [di][emu_fraction][w1] mov cx, [di][emu_fraction][w2] mov dx, [di][emu_fraction][w3] mov byte [bp][excess], 0 sub byte [bp][align], 8 jl @align_bits @align_bytes: mov [bp][excess], al mov al, ah mov ah, bl mov bl, bh mov bh, cl mov cl, ch mov ch, dl mov dl, dh mov dh, 0 sub byte [bp][align], 8 jge @align_bytes @align_bits: and byte [bp][align], 7 jz @align_end @align_bit_shift: shr dx, 1 rcr cx, 1 rcr bx, 1 rcr ax, 1 rcr byte [bp][excess], 1 dec byte [bp][align] jnz @align_bit_shift @align_end: sub di, di (* convenient zero*) cmp byte [bp][op], positive_sign (* add or subtract*) jne @share_do_subtract (* start of addition*) @share_do_add: add ax, [si][emu_fraction][w0] adc bx, [si][emu_fraction][w1] adc cx, [si][emu_fraction][w2] adc dx, [si][emu_fraction][w3] jnc @shared_round (* normalise if overflowed*) rcr dx, 1 rcr cx, 1 rcr bx, 1 rcr ax, 1 rcr byte [bp][excess], 1 inc word [bp][zExp] jmp @shared_round (* end of addition*) (* start of subtraction*) @share_do_subtract: xor byte [bp][zSign], negative_sign sub ax, [si][emu_fraction][w0] (* smaller - larger*) sbb bx, [si][emu_fraction][w1] (* is the*) sbb cx, [si][emu_fraction][w2] (* negative of*) sbb dx, [si][emu_fraction][w3] (* the expected result*) (* Since we compared only exponents when checking which is "larger", in*) (* the case where exponents were equal there is a chance the fraction is*) (* smaller. Since the smaller is supposed to be in registers, the subtraction*) (* would give a negative result (carry=borrow bit flagged) if they are the*) (* correct way around.*) jnb @sub_align (* "smaller" actually larger*) (* "larger" was larger, so negate to get correct result.*) xor byte [bp][zSign], negative_sign not dx not cx (* use the identity -x = not(x) + 1*) not bx not ax neg byte [bp][excess] cmc adc ax, di adc bx, di adc cx, di adc dx, di @sub_align: mov si, 64 or dh, dh js @shared_round (* Bit shifting appears inefficient. However, note that only when x and y*) (* are nearly equal will their difference require large realignments, and*) (* such cases are unlikely. It does not seem worthwhile to provide byte*) (* or word shifts.*) @sub_bit_shift: dec si jz @sub_result_zero shl byte [bp][excess], 1 rcl ax, 1 rcl bx, 1 rcl cx, 1 adc dx, dx jns @sub_bit_shift sub si, 64 add [bp][zExp], si jmp @shared_round @sub_result_zero: mov word [bp][zExp], zero_exponent mov byte [bp][zSign], positive_sign jmp @shared_end (* end of subtraction*) (* shared code for both add and subtract*) @shared_round: shl byte [bp][excess], 1 adc ax, di adc bx, di adc cx, di adc dx, di jnc @roundEnd rcr dx, 1 (* note dx:cx:bx:ax must be 8000:0:0:0*) inc word [bp][zExp] (* if you get to here*) @roundEnd: cmp word [bp][zExp], infinite_exponent jge @overflow cmp word [bp][zExp], zero_exponent jle @underflow @shared_end: cld mov di, [bp][z_p] stosw xchg ax, bx stosw xchg ax, cx stosw xchg ax, dx stosw mov ax, [bp][zExp] stosw mov al, [bp][zSign] stosb pop di; pop si; mov sp,bp; pop bp ret paramsize @overflow: mov ch, fe_Overflow mov word [bp][zExp], infinite_exponent jmp @extreme @underflow: mov ch, fe_Underflow mov word [bp][zExp], zero_exponent @extreme: call e287_Exception mov byte [bp][zSign], positive_sign (* extremes always positive*) sub ax, ax mov bx, ax mov cx, ax mov dx, ax jmp @shared_end (*INCLUDE E287ATAN.ASM*) pi_by_2 : dw 0C235H, 2168H, 0DAA2H, 0C90FH, 1; db positive_sign, 0 (* Polynomial coefficients for Chebyshev approx. to ArcTan (x) / x,*) (* in terms of every second power.*) poly_terms : dw 8 poly_coeffs : dw 0H, 0H, 0H, 0H, 0H; db positive_sign, 0 (* 1.0*) dw 05557H, 05555H, 05555H, 0155H, 0H; db negative_sign, 0 dw 032BEH, 03333H, 03333H, 03H, 0H; db positive_sign, 0 dw 01E7DH, 09249H, 0924H, 0H, 0H; db negative_sign, 0 dw 0FEBCH, 071C6H, 01CH, 0H, 0H; db positive_sign, 0 dw 0FF5CH, 05D16H, 0H, 0H, 0H; db negative_sign, 0 dw 0BBD8H, 013AH, 0H, 0H, 0H; db positive_sign, 0 dw 00C0AH, 04H, 0H, 0H, 0H; db negative_sign, 0 numOfIntervals = 4 topInterval = 3 topsOf : dw 00000H, 00000H, 00000H, 0E700H, -3; db positive_sign, 0 dw 0A4BDH, 07BD6H, 064EEH, 0B35CH, -1; db positive_sign, 0 dw 085B5H, 0FC47H, 03074H, 0A111H, 0; db positive_sign, 0 one : dw 00000H, 00000H, 00000H, 08000H, 1; db positive_sign, 0 arcTans : dw 0FA9CH, 0B064H, 01DB2H, 0E607H, -2; db positive_sign, 0 dw 0FA9CH, 0B064H, 01DB2H, 0E607H, -1; db positive_sign, 0 dw 0BBF5H, 0044BH, 05646H, 0AC85H, 0; db positive_sign, 0 middles : dw 06EE6H, 01FD9H, 009BDH, 0E9FAH, -2; db positive_sign, 0 dw 07B8DH, 0BD35H, 0845BH, 0F6DDH, -1; db positive_sign, 0 dw 0D57FH, 00235H, 080DBH, 0CC73H, 0; db positive_sign, 0 (* algorithm notes*) (* Firstly we take the ratio z = x / y, with x and y the parameters. We only*) (* need z.*) (* A Chebyshev polynomial expansion is used for a very narrow domain*) (* of arguments, -1/N .. +1/N, where N is on the order of 8. The remaining*) (* range +1/N .. 1.0 is mapped onto the smaller interval.*) (* Beginning with high-school trigonometry:*) (* Tan (a + b) = (Tan (a) + Tan (b)) / (1 - Tan(a) * Tan(b))*) (* We seek to find:*) (* psi = ArcTan (z)*) (* let z = (centre + v) / (1 - centre * v)*) (* where the closest value of centre is chosen by comparing z to a table*) (* of bounds of intervals from 0 to infinity, and v is a number in the*) (* range -1/N .. +1/N. The intervals and centres are pre-computed to*) (* ensure that v will be in this range.*) (* We can calculate v by re-arranging the algebra:*) (* v = (z - centre) / (1 + centre * z)*) (* The importance of this formula is that if*) (* theta = ArcTan (centre) -- pre-calculated table of theta values*) (* phi = ArcTan (v) -- calculated by rapid polynomial*) (* then*) (* ArcTan (z) = theta + phi*) (* which is the answer we seek.*) (* VAR*) (* sign : (positive_sign, negative_sign);*) (* theta : 0 .. topInterval;*) (* result : emu_temp;*) (* BEGIN*) (* z := x / y;*) (* theta := 0;*) (* WHILE (theta < top_interval) AND (x > topsOf [theta]) DO theta := theta + 1;*) (* IF theta = 0 THEN result := ArcTan (z)*) (* ELSE*) (* theta := theta - 1;*) (* z := (z - middles [theta]) / (1 + z * middles [theta]);*) (* result := ArcTan (z) + arcTans [theta];*) (* emu_PAtan := result;*) (* END;*) (*public*) emu_PAtan : (*NEAR*) _theta = -2 _temp = _theta - emu_temp_size push bp; mov bp,sp; lea sp,[bp][_temp]; push si; push di push si push di push si call emu_Divide mov byte [bp][_theta], 0 mov di, (*offset*) topsOf (* di := ^ topsOf [theta]*) at@while: cmp byte [bp][_theta], topInterval jae at@whileEnd push di (*1*) push si call copy_di_to_temp push di call emu_Compare pop di (*1*) cmp ax, e287_greater jne at@whileEnd inc byte [bp][_theta] add di, emu_temp_size jmp at@while at@whileEnd: cmp byte [bp][_theta], 0 jne at@notSmall at@small: call ArcTan (* _result := ArcTan (si)*) jmp at@normalEnd at@notSmall: dec byte [bp][_theta] (* theta := theta - 1;*) (* z := (z - middles [theta]) / (1 + z * middles [theta]);*) mov al, emu_temp_size mul byte [bp][_theta] add ax, (*offset*) middles xchg di, ax sub sp, emu_temp_size mov ax, sp push si call copy_di_to_temp push di push ax call emu_Subtract push si push di push si call emu_Multiply mov di, one call copy_di_to_temp push di push si push si call emu_Add mov ax, sp push ax push si push si call emu_Divide add sp, emu_temp_size call ArcTan mov al, emu_temp_size mul byte [bp][_theta] add ax, (*offset*) arcTans xchg di, ax call copy_di_to_temp push di push si push si call emu_Add at@normalEnd: at@end: pop di; pop si; mov sp,bp; pop bp ret 0 copy_di_to_temp: push ax mov ax, cs:[di] mov [bp][_temp],ax mov ax, cs:[di][2] mov [bp][_temp][2],ax mov ax, cs:[di][4] mov [bp][_temp][4],ax mov ax, cs:[di][6] mov [bp][_temp][6],ax mov ax, cs:[di][8] mov [bp][_temp][8],ax mov ax, cs:[di][0AH] mov [bp][_temp][0AH],ax lea di, [bp][_temp] pop ax ret 0 (* FUNCTION ArcTan ( _z : emu_temp ) _result : emu_temp;*) (* { Polynomial evaluation of ArcTan of a number within a small domain*) (* centred upon zero.*) (* BEGIN*) (* _result := _z * Poly ( _z * _z, poly_terms, poly_coeffs);*) (* END;*) ArcTan : (* parameter and result both at [si]*) linear_bound = -32 push si push di call emup_Push cmp word [si][emu_exponent], linear_bound jg art@nonLinear mov di, si call emu_Pop (* ArcTan (x) ~= x for very small x*) jmp art@end art@nonLinear: mov di, sp add word [di][emu_exponent], 3 (* the polynomial is designed*) (* for a *8 magnified parameter.*) call emup_Square_Fix push cs:[poly_terms] (* number of terms to evaluate*) mov ax, poly_coeffs push cs push ax call emup_Polynomial push di push si push si call emu_Multiply add sp, emu_temp_size (* pop Poly value off stack*) art@end: pop di pop si ret 0 (*INCLUDE E287CNST.ASM*) (* one: dw 0, 0, 0, 8000H, 1; db positive_sign, 0 *) Log2of10: dw 08AFEH, 0CD1BH, 0784BH, 0D49AH, 2; db positive_sign, 0 Log2ofE: dw 0F0BBH, 05C17H, 03B29H, 0B8AAH, 1; db positive_sign, 0 Pi: dw 0C235H, 2168H, 0DAA2H, 0C90FH, 2; db positive_sign, 0 Log10of2: dw 0F799H, 0FBCFH, 09A84H, 09A20H, -1; db positive_sign, 0 LogEof2: dw 079ACH, 0D1CFH, 017F7H, 0B172H, 0; db positive_sign, 0 zero: dw 0, 0, 0, 0, zero_exponent; db positive_sign, 0 (* All constants' procedures are in0-out1 pattern (see e287long.asm) and*) (* thus take no input and place their output at ss=ds=es:[di]*) (*public*) emu_1 : (*NEAR*) mov ax, (*offset*) one jmp MoveConstant (*public*) emu_Log2of10 : (*NEAR*) mov ax, (*offset*) Log2of10 jmp MoveConstant (*public*) emu_Log2ofE : (*NEAR*) mov ax, (*offset*) Log2ofE jmp MoveConstant (*public*) emu_Pi : (*NEAR*) mov ax, (*offset*) Pi jmp MoveConstant (*public*) emu_Log10of2 : (*NEAR*) mov ax, (*offset*) Log10of2 jmp MoveConstant (*public*) emu_LogEof2 : (*NEAR*) mov ax, (*offset*) LogEof2 jmp MoveConstant (*public*) emu_0 : (*NEAR*) mov ax, (*offset*) zero (* jmp MoveConstant*) MoveConstant: push ds push cs pop ds xchg si, ax cld mov cx, emu_temp_size / 2 rep; movsw sub di, emu_temp_size (* restore original di*) xchg si, ax (* and si values*) pop ds ret 0 (*public*) emu_Extreme : (*NEAR*) push di mov es:[di][emu_exponent], ax (* infinity or zero, chosen by caller*) mov byte es:[di][emu_sign], positive_sign sub ax, ax cld stosw (* zero fraction*) stosw stosw stosw pop di ret 0 (*INCLUDE E287CNVT.ASM*) (* procedure to float an integer to yield a long temporary real*) (* input is an integer at es : [si]*) (* result is a temporary real at ds : [di]*) (* no exceptions.*) (*public*) emu_Float : (*NEAR*) mov ax, es : [si][w0] sub cx, cx (* convenient zero*) mov dx, positive_sign or ax, ax jl @fl_negative jg @fl_signed @fl_zero: mov bx, zero_exponent jmp @fl_end @fl_negative: neg ax mov dl, negative_sign @fl_signed: sub bx, bx xchg cx, ax @fl_shift: inc bx shr cx, 1 rcr ax, 1 jcxz @fl_aligned jmp @fl_shift @fl_aligned: @fl_end: mov [di][emu_sign][w0], dx mov [di][emu_exponent], bx mov [di][emu_fraction][w3], ax mov [di][emu_fraction][w2], cx (* cx = zero*) mov [di][emu_fraction][w1], cx mov [di][emu_fraction][w0], cx ret 0 (*public*) emu_Round : (*NEAR*) (* Convert a long temporary floating point number to a 16 bit integer*) (* with fraction rounded according to current rounding control.*) (* Input operand at ds:[si]*) (* Result at es:[di]*) (* Side effects: may exit via exception*) mov cx, [si][emu_exponent] cmp cx, 15 jg @ro_overflow cmp cx, zero_exponent jg @ro_standard @ro_zero: sub ax, ax jmp @ro_end @ro_overflow: mov ch, fe_Operand_Too_Big call e287_Exception mov ax, 8000H jmp @ro_end @ro_standard: mov bx, [si][emu_fraction][w3] (* ax:bx holds the*) sub ax, ax (* integer:fraction*) mov dx, ax (* dx collects the sticky bits*) cmp cx, 0 jge @ro_getSticky shr bx, 1 (* ensure tiny numbers are*) rcr dx, 1 (* aligned <= 1/4*) @ro_getSticky: or dx, [si][emu_fraction][w0] (* stick the lesser*) or dx, [si][emu_fraction][w1] (* words all*) or dx, [si][emu_fraction][w2] (* together*) cmp cx, 0 jle @ro_bitAligned @ro_bitLoop: shl bx, 1 rcl ax, 1 loop @ro_bitLoop @ro_bitAligned: or bl, dh (* BL is completely sticky*) or bl, dl (* bx is fraction byte + sticky byte*) mov cl, e287_roundingControl and cl, [emws_control][by1] (* get rounding control*) cmp cl, e287_chop je @ro_chop cmp cl, e287_round je @ro_roundNear add cl, [si][emu_sign][by0] cmp cl, e287_floor + positive_sign je @ro_chop cmp cl, e287_ceiling + negative_sign je @ro_chop @ro_roundUp: (* this is what sticky bits are for !*) neg bx (* sets carry if any bit not zero*) adc ax, 0 js @ro_overflow jmp @ro_setSign @ro_roundNear: mov dl, 1 (* banker's rounding needs test*) and dl, al (* of ls bit of integer*) or bl, dl add bx, 7FFFH adc ax, 0 js @ro_overflow @ro_chop: @ro_setSign: cmp byte [si][emu_sign], negative_sign jne @ro_end neg ax @ro_end: mov es : [di][w0], ax ret 0 (* algorithm notes*) (* The scaling operation is quite simple, just adding an integer to*) (* the exponent of the floating point number. The complications are*) (* that the "integer" is a real in tos(1) which must first be truncated,*) (* then to check against scaling zeroes and infinities, and to report*) (* on underflows or overflows created.*) (*public*) emu_Scale : (*NEAR*) push bp; mov bp,sp; push si; push di (* first we truncate sc_n*) mov cx, [si][emu_exponent] cmp cx, 15 jg sc@overScale or cx, cx jg sc@non_zero sub ax, ax jmp sc@truncated sc@overScale: mov ch, fe_Operand_Too_Big call e287_Exception mov ax, 7FFFH (* largest signable magnitude*) jmp sc@signScale sc@non_zero: mov ax, [si][emu_fraction][w3] neg cl add cl, 16 shr ax, cl sc@signScale: cmp byte [si][emu_sign], negative_sign jne sc@truncated neg ax sc@truncated: (* now apply scale to tos' exponent*) mov cx, [di][emu_exponent] cmp cx, zero_exponent jle sc@end (* zero * anything = zero*) cmp cx, infinite_exponent jge sc@end (* infinity cannot be changed*) add cx, ax jo sc@under_or_overflow cmp cx, zero_exponent jle sc@underflow cmp cx, infinite_exponent jge sc@overflow mov [di][emu_exponent], cx sc@end: pop di; pop si; pop bp ret 0 sc@under_or_overflow: test ax,ax jl @underflow sc@overflow: mov ch, fe_Overflow call e287_Exception mov ax, infinite_exponent jmp sc@extreme sc@underflow: mov ch, fe_Underflow call e287_Exception mov ax, zero_exponent sc@extreme: call emu_Extreme jmp sc@end (* (*INCLUDE E287COMP1.ASM*) (* second parameter in cs *) (*public*) emu_Compare1 : (*NEAR*) co_x_p = 2+nearParam (* based on ds*) co_y_p = 0+nearParam (* based on cs*) co_paramsize = 4 (* result in [emws_status]*) push bp; mov bp,sp; push si; push di mov si, [bp][co_x_p] mov di, [bp][co_y_p] mov ax, [si][emu_exponent] cmp ax, cs:[di][emu_exponent] jge @1ax_is_greater mov ax, cs:[di][emu_exponent] @1ax_is_greater: cmp ax, zero_exponent jle @1both_equal (* both must be zero*) cmp ax, infinite_exponent jge @1incomparable (* infinities are incomparable*) mov cl, [si][emu_sign] cmp cl, cs:[di][emu_sign] jl @1x_greater (* positive_sign < negative_sign*) jg @1x_lesser mov ax, [si][emu_exponent] cmp ax, cs:[di][emu_exponent] jl @1x_smaller jg @1x_larger mov ax, [si][emu_fraction][w3] cmp ax, cs:[di][emu_fraction][w3] jne co@1notEqual mov ax, [si][emu_fraction][w2] cmp ax, cs:[di][emu_fraction][w2] jne co@1notEqual mov ax, [si][emu_fraction][w1] cmp ax, cs:[di][emu_fraction][w1] jne co@1notEqual mov ax, [si][emu_fraction][w0] cmp ax, cs:[di][emu_fraction][w0] jne co@1notEqual @1both_equal: mov ax, e287_equal jmp @1co_end co@1notEqual: ja @1x_larger @1x_smaller: cmp cl, positive_sign jne @1x_greater @1x_lesser: mov ax, e287_lesser jmp @1co_end @1x_larger: cmp cl, positive_sign jne @1x_lesser @1x_greater: mov ax, e287_greater jmp @1co_end @1incomparable: mov ax, e287_incomparable @1co_end: mov [emws_status][by1], ah pop di; pop si; pop bp ret co_paramsize *) (*INCLUDE E287COMP.ASM*) (* algorithm notes:*) (* first check for infinities. Infinities are incomparable.*) (* Next check for sign. If signs differ, and not both numbers are zero,*) (* then the positive number is the greater.*) (* Check for the greater absolute value, beginning with the exponent then*) (* working down through the fraction. If not both are equal, then if both*) (* are positive the greater maginitude is the greater, else if negative*) (* the greater magnitude is lesser.*) (*public*) emu_Compare : (*NEAR*) co_x_p = 2+nearParam (* based on ds*) co_y_p = 0+nearParam (* based on ds*) co_paramsize = 4 (* result in [emws_status]*) push bp; mov bp,sp; push si; push di mov si, [bp][co_x_p] mov di, [bp][co_y_p] mov ax, [si][emu_exponent] cmp ax, [di][emu_exponent] jge @ax_is_greater mov ax, [di][emu_exponent] @ax_is_greater: cmp ax, zero_exponent jle @both_equal (* both must be zero*) cmp ax, infinite_exponent jge @incomparable (* infinities are incomparable*) mov cl, [si][emu_sign] cmp cl, [di][emu_sign] jl @x_greater (* positive_sign < negative_sign*) jg @x_lesser mov ax, [si][emu_exponent] cmp ax, [di][emu_exponent] jl @x_smaller jg @x_larger mov ax, [si][emu_fraction][w3] cmp ax, [di][emu_fraction][w3] jne co@notEqual mov ax, [si][emu_fraction][w2] cmp ax, [di][emu_fraction][w2] jne co@notEqual mov ax, [si][emu_fraction][w1] cmp ax, [di][emu_fraction][w1] jne co@notEqual mov ax, [si][emu_fraction][w0] cmp ax, [di][emu_fraction][w0] jne co@notEqual @both_equal: mov ax, e287_equal jmp @co_end co@notEqual: ja @x_larger @x_smaller: cmp cl, positive_sign jne @x_greater @x_lesser: mov ax, e287_lesser jmp @co_end @x_larger: cmp cl, positive_sign jne @x_lesser @x_greater: mov ax, e287_greater jmp @co_end @incomparable: mov ax, e287_incomparable @co_end: mov [emws_status][by1], ah pop di; pop si; pop bp ret co_paramsize (* A very simple procedure.*) (* Compare_Zero :=*) (* IF x = 0 THEN equal*) (* ELSE x = infinity THEN incomparable*) (* ELSE x < 0 THEN lesser*) (* ELSE greater;*) (* with apologies to Algol-68, Prof. Dijkstra, and P. B. Hansen for the syntax.*) (*public*) emu_Compare_Zero : (*NEAR*) (* result in ax and [emws_status]*) mov ax, e287_equal cmp word [si][emu_exponent], zero_exponent jle cz@end mov ax, e287_incomparable cmp word [si][emu_exponent], infinite_exponent jge cz@end mov ax, e287_lesser cmp byte [si][emu_sign], negative_sign je cz@end mov ax, e287_greater cz@end: mov [emws_status][by1], ah ret 0 (* Another simple procedure.*) (* XAM :=*) (* IF x = 0 THEN plusZero*) (* ELSE x = infinity THEN plusInfinity*) (* ELSE x < 0 THEN minusNormal*) (* ELSE plusNormal;*) (*public*) emu_XAM : (*NEAR*) (* result in ax and [emws_status]*) mov ax, e287_plusZero cmp word [si][emu_exponent], zero_exponent jle xa@end mov ax, e287_plusInfinity cmp word [si][emu_exponent], infinite_exponent jge xa@end mov ax, e287_minusNormal cmp byte [si][emu_sign], negative_sign je xa@end mov ax, e287_plusNormal xa@end: mov [emws_status][by1], ah ret 0 (*INCLUDE E287DCML.ASM*) (* Packed Decimal Format*) (* The format used by the 8087/287 coprocessors uses 10 bytes to hold a*) (* signed 18-digit decimal integer. The sign occupies the 79th (most*) (* significant) bit (1 = negative), the bits 72..78 are don't cares (always*) (* generated as zeroes by this software), and the 4-bit groups 68..71,*) (* 64..67, .., 0..3 are the decimal digits in decreasing significance.*) (* The number is an integer, so the least significant digit is always*) (* interpreted as the units' digit. The Intel co-processors do not*) (* directly handle decimal fractions or floating points.*) (*packedDecimal STRUC*) pdec_digits = 0 (*db 9 dup (?)*) pdec_sign = 9 (*db ?*) (*packedDecimal ENDS*) (* FUNCTION emu_Push_Decimal ( VAR dec : packedDecimal) : emu_temp;*) (* VAR*) (* res : emu_temp;*) (* FUNCTION WordDecToBin ( posn : cardinal );*) (* BEGIN*) (* WordDecToBin := ((dec.digit [posn + [3] * 10 +*) (* dec.digit [posn + [2]) * 10 +*) (* dec.digit [posn + [1]) * 10 + dec.digit [posn];*) (* END;*) WordDecToBin : (* posn supplied in si*) (* result in ax*) push bx push cx push dx mov cl, 4 mov ch, 10 mov bx, es:[si][w0] mov al, bh shr al, cl mul ch (* ax in range 0..90*) mov dl, 0FH and dl, bh add al, dl mul ch (* ax in range 0..990*) mov dx, 0F0H and dl, bl shr dx, cl add ax, dx mov cx, 10 mul cx (* ax in range 0..9990*) and bx, 0FH add ax, bx pop dx pop cx pop bx ret 0 (* BEGIN*) (* res.sign := dec.sign;*) (* res.fraction := dec.digit [16] + tens [dec.digit [17];*) (* res.fraction := res.fraction * 10000 + WordDecToBin (12);*) (* res.fraction := res.fraction * 10000 + WordDecToBin (8);*) (* res.fraction := res.fraction * 10000 + WordDecToBin (4);*) (* res.fraction := res.fraction * 10000 + WordDecToBin (0);*) (* res.exponent := 64;*) (* Normalise (res);*) (* END;*) (*public*) emu_Push_Decimal : (*NEAR*) (* inputs: packed decimal at es:[si]*) (* outputs: emu_temp at ds:[di]*) (* side effects: no exceptions reported.*) push bp mov bp, sp push si pu?decPtr = -2 (* saved value of si*) push di pu?resPtr = -4 (* saved value of di*) (* res.sign := dec.sign;*) mov al, es : [si][pdec_sign] and ax, 80H (* separate the sign bit*) rol al, 1 (* into bit 0*) mov [di][emu_sign][w0], ax (* res.fraction := dec.digit [16] + tens [dec.digit [17];*) mov cl, 4 (* mov ah, 0 ; done above*) mov al, es:[si][8][by0] shl ax, cl shr al, cl aad (* AL := AL + 10 * AH, AH := 0*) (* res.fraction := res.fraction * 10000 + WordDecToBin (12);*) mov di, 10000 mul di xchg bx, ax (* dx:bx in range 0..990000*) lea si, [si][w3] call WordDecToBin add ax, bx adc dl, dh (* dh=0 dx:ax in range 0..999,999*) (* res.fraction := res.fraction * 10000 + WordDecToBin (8);*) mov bx, dx mul di xchg ax, bx mov cx, dx mul di add cx, ax adc dl, dh (* dh=0 dx:cx:bx in range 0..9,999,990,000*) sub si, 2 call WordDecToBin add ax, bx adc cx, 0 adc dl, dh (* dh=0 dx:cx:ax in range 0..9,999,999,999*) (* res.fraction := res.fraction * 10000 + WordDecToBin (4);*) push si mov bx, dx (* bx:cx:ax*) mul di xchg ax, cx mov si, dx (* si:cx*) mul di xchg ax, bx xchg di, dx (* di:bx*) mul dx (* dx:ax*) add bx, si adc di, ax (* dx = 0 di:bx:cx range 0..10^14 < 2^47*) pop si sub si, 2 call WordDecToBin add ax, cx adc bx, dx (* dx = 0*) adc di, dx (* di:bx:ax range 0..10^14*) (* res.fraction := res.fraction * 10000 + WordDecToBin (0);*) mov si, 10000 mul si xchg ax, bx mov cx, dx (* cx:bx*) mul si xchg ax, si xchg di, dx (* di:si*) mul dx (* dx:ax*) add cx, si adc di, ax adc dx, 0 (* dx:di:cx:bx range 0..10^18 < 2^60*) mov si, [bp][pu?decPtr] call WordDecToBin add bx, ax adc cx, 0 adc di, 0 adc dx, 0 (* res.exponent := 64;*) mov ax, 64 (* Normalise (res);*) pu@alignWord: or dx, dx jnz pu@alignBit sub ax, 16 jz pu@zero xchg dx, di xchg di, cx xchg cx, bx jmp pu@alignWord pu@alignBit: js pu@finish pu@anotherBit: dec ax shl bx, 1 rcl cx, 1 rcl di, 1 adc dx, dx jns pu@anotherBit pu@finish: mov si, [bp][pu?resPtr] mov [si][emu_fraction][w0], bx mov [si][emu_fraction][w1], cx mov [si][emu_fraction][w2], di mov [si][emu_fraction][w3], dx mov [si][emu_exponent], ax pu@end: pop di pop si mov sp, bp pop bp ret 0 pu@zero: mov ax, zero_exponent jmp pu@finish public __emul_fbstp : TempReal = -10 locals = 10 push bp mov bp, sp sub sp, locals push bx push cx push dx push di push si push ds fld st(0), st(0) fstp tbyte [bp][TempReal], st(0) push ss pop ds push ss pop es lea si, [bp][TempReal] mov di, ax call emu_Pop_Decimal pop ds pop si pop di pop dx pop cx pop bx mov sp, bp pop bp ret far 0 (* PROCEDURE emu_Pop_Decimal ( z : emu_temp ;*) (* VAR dec : packedDecimal );*) (* Convert a floating point temporary to a packed decimal integer.*) (* VAR*) (* exp : integer;*) (* frac : cardinal64;*) (* res ; packedDecimal;*) (* PROCEDURE PutFourDigits ( rem, posn : cardinal );*) (* BEGIN*) (* res.digits [posn] := rem MOD 10;*) (* rem := rem DIV 10;*) (* res.digits [posn+[1] := rem MOD 10;*) (* rem := rem DIV 10;*) (* res.digits [posn+[2] := rem MOD 10;*) (* res.digits [posn+[3] := rem DIV 10;*) (* END;*) PutFourDigits : (* rem in dx, posn in di*) push ax push cx mov al, 100 mov cl, 4 xchg ax, dx div dl mov dl, ah aam (* AH, AL := AL DIV 10, AL MOD 10*) shl ah, cl or ah, al xchg ax, dx aam shl ah, cl or al, ah mov ah, dh stosw pop cx pop ax ret 0 (* BEGIN*) (* exp := z.exp;*) (* frac := z.frac;*) (* IF exp < 1 THEN res := 0*) (* ELSE exp > 60 THEN res := +infinity*) (* ELSE*) (* WHILE exp < 64 DO*) (* frac := frac DIV 2;*) (* exp := exp + 1;*) (* frac := frac + carry; { rounding }*) (* res := 0;*) (* PutFourDigits (frac MOD 10000, 0);*) (* frac := frac DIV 10000;*) (* PutFourDigits (frac MOD 10000, 4);*) (* frac := frac DIV 10000;*) (* PutFourDigits (frac MOD 10000, 8);*) (* frac := frac DIV 10000;*) (* PutFourDigits (frac MOD 10000, 12);*) (* frac := frac DIV 10000;*) (* res.digits [16] := frac MOD 10;*) (* res.digits [17] := frac DIV 10;*) (* res.sign := z.sign;*) (* emu_Pop_Decimal := res;*) (* END;*) (*public*) emu_Pop_Decimal : (*NEAR*) push bp mov bp, sp push si (* ds:si points to input*) push di (* es:di used to point to result*) cld (* forward string direction*) (* exp := z.exp;*) (* frac := z.frac;*) mov ax, [si][emu_exponent] mov bx, [si][emu_fraction][w0] mov cx, [si][emu_fraction][w1] mov dx, [si][emu_fraction][w3] mov si, [si][emu_fraction][w2] (* IF exp < 1 THEN res := 0*) (* { will the result exceed 999,999,999,999,999,999 ? }*) (* ELSE z > [3C, DE0B 6B3A 763F [FFF0] THEN res := +infinity*) cmp ax, 0 jl po@underflowJmp sub ax, 3CH jl po@inRange jg po@overflowJmp cmp dx, 0DE0BH jb po@inRange ja po@overflowJmp cmp si, 06B3AH jb po@inRange ja po@overflowJmp cmp cx, 0763FH jb po@inRange ja po@overflowJmp cmp bx, 0FFF0H ja po@overflowJmp (* ELSE*) (* WHILE exp < 64 DO*) (* frac := frac DIV 2;*) (* exp := exp + 1;*) po@inRange: mov ah, 0 (* excess precision for rounding*) sub al, 4 po@wordAlign: add al, 16 jg po@bitAlign mov ah, bh mov bx, cx mov cx, si mov si, dx sub dx, dx jmp po@wordAlign (*----*) po@underflowJmp: jmp po@underflow po@overflowJmp: jmp po@overflow (*----*) po@bitAlign: sub al, 16 jnl po@aligned po@anotherBit: shr dx, 1 rcr si, 1 rcr cx, 1 rcr bx, 1 rcr ah, 1 inc al jl po@anotherBit po@aligned: (* frac := frac + carry; { rounding }*) add ah, ah adc bx, 0 adc cx, 0 adc si, 0 adc dx, 0 (* dx:si:cx:bx is now an integer < 2^60*) (* PutFourDigits (frac MOD 10000, 0);*) (* frac := frac DIV 10000;*) (* dx < 4096 since 2^60 = 4096 * 2^48*) xchg ax, si mov si, 10000 div si xchg ax, cx div si xchg ax, bx div si (* frac = cx:bx:ax < 2^47, dx = mod 10000*) call PutFourDigits (* dx is parameter*) (* PutFourDigits (frac MOD 10000, 4);*) (* frac := frac DIV 10000;*) sub dx, dx xchg ax, cx div si xchg ax, bx div si xchg ax, cx div si (* frac = bx:cx:ax < 2^34, dx = mod 10000*) call PutFourDigits (* dx is parameter*) (* PutFourDigits (frac MOD 10000, 8);*) (* frac := frac DIV 10000;*) mov dx, bx xchg ax, cx div si xchg ax, cx div si (* frac = cx:ax < 2^20, dx = mod 10000*) call PutFourDigits (* dx is parameter*) (* PutFourDigits (frac MOD 10000, 12);*) (* frac := frac DIV 10000;*) mov dx, cx div si (* frac = ax < 2^7, dx = mod 10000*) call PutFourDigits (* dx is parameter*) (* res.digits [16] := frac MOD 10;*) (* res.digits [17] := frac DIV 10;*) aam (* AH, AL := AL DIV 10, AL MOD 10*) mov cl, 4 shl ah, cl or al, ah stosb po@end: (* res.sign := z.sign;*) mov si, [bp][-2][w0] mov al, [si][emu_sign] ror al, 1 stosb pop di pop si pop bp ret 0 (* ----*) po@underflow: mov al, 0 po@extreme: mov cx, 9 rep; stosb jmp po@end po@overflow: mov ch, fe_Overflow call e287_Exception mov al, 99H jmp po@extreme (* ----*) (*INCLUDE E287DIV.ASM*) (* algorithm notes.*) (* We are looking for the solution to the ratio:*) (* aNNN + bNN + cN + d*) (* ---------------------- = (vNNN +wNN + xN + y) / K*) (* pNNN + qNN + rN + s*) (* Where a, p, v are in the range N/2 .. N-1 and the other lower-case*) (* symbols range from 0..N-1 , and K is either NNNN or NNNN/2.*) (* The machine provides a primitive for NN / N division.*) (* Knuth proposed a method for double length divide which uses an*) (* approximation:*) (* AN + B AN + B PN AN + B*) (* ------ = ------ x ------ = ------ x ( 1 - Q/PN + (Q/PN)^2 ...)*) (* PN + Q PN PN + Q PN*) (* We can use the same algebra, noting that B = bcd, Q = qrs.*) (* The algebra is better changed into a finite equation:*) (* AN + B = ((PN + Q) * V) - error*) (* = (PN*V + QV) - error*) (* Now if V = (AN + B) / PN then we have*) (* error = QV*) (* Unfortunately, calculating QV is not likely to be fast enough, as*) (* in our case it involves roughly 12 multiplications plus some additions.*) (* A faster method is to be more crude, and accept only a rounded*) (* version of V such that*) (* V + Ev = (AN + B) / PN*) (* This gives the equation*) (* error = (PN + Q)*V - (AN + B)*) (* It will be convenient for the error to have a predictable sign, so*) (* V is always rounded up. The truncation PN also causes an error in*) (* the same direction. It saves code to know which direction to correct*) (* at each stage, and will not actually change speed.*) (* If another stage is applied, so that*) (* W = error / P*) (* is the next estimate, then a new error will remain*) (* error(2) = error(1) - (PN + Q) * W*) (* and that can be substituted to give*) (* (PN + Q)*W + error(2) = (PN + Q)*V - (AN + B)*) (* error(2) = (PN + Q)(V - W) - (AN + B)*) (* Generalised to 5 stages the result is*) (* error(5) = (PN + Q)(V - W + X - Y + Z) - (AN + B)*) (* The crucial question is, how large are the error terms at each stage ?*) (* If we look at the first stage, it can be seen that the error term is*) (* made up of two parts. Let us introduce the symbol*) (* V' = (AN + B) / P accurately*) (* V = RoundUp(V')*) (* Then the error has two components:*) (* error(1) = V'*Q + (V - V') * (PN + Q)*) (* For the floating point fractions we have:*) (* NNNN/2 <= AN + B < NNNN and N/2 <= P < N*) (* (don't forget to normalise by dividing products by NNNN).*) (* so the first component of V has bounds*) (* NNNN/2 < V < 2NNNN*) (* and that leads to a first error term of*) (* 0 < V * Q < 2NNN*) (* The second component depends on (V - V'), which ranges over*) (* 0 <= V' - V < NNN*) (* so the second component of the error ranges*) (* 0 <= (PN + Q) * (V' - V) < NNN*) (* The two terms are almost independent, so the total error range will*) (* be*) (* 0 < error < 3NNN*) (* More succinctly, the error is up to 3 times the least bit.*) (* In fact, that will be our error at every stage since the code uses*) (* a double length estimate for V, W, X, Y, Z, where the most significant*) (* word acts as a guard so that later stages are not complicated by over-*) (* flow from early stages. The 4'th stage will thus have an error of*) (* up to 3 LSB and biased positive, which is acceptable if intended for*) (* use only to support an HLL which does not provide temporary-precision*) (* variables accessible to the user.*) numP = 4+nearParam (* based on ds*) dvrP = 2+nearParam (* based on ds*) quoP = 0+nearParam (* based on es*) div_paramsize = 6 quoExp = -2 quoSign = -3 excess = -4 result = -12 temp = -16 guard = -18 (* note that temp and result, both quad-word objects, overlap by two words.*) (* This is allowed because the most significant words of temp are discarded*) (* as the least significant words of result are generated.*) (*public*) emu_Divide : (*NEAR*) push bp; mov bp,sp; lea sp, [bp][temp][w3][w1] mov bx, [bp][numP] push [bx][emu_fraction][w3] (* copy fraction to*) push [bx][emu_fraction][w2] push [bx][emu_fraction][w1] (* temp work area*) push [bx][emu_fraction][w0] sub ax, ax push ax (* zero guard*) push si push di cld mov di, [bp][dvrP] mov si, [di][emu_exponent] mov ax, [bx][emu_exponent] (* first check for all the extreme cases*) cmp si, zero_exponent jle @div_by_zero cmp ax, infinite_exponent jge @div_infinite_result cmp ax, zero_exponent jle @div_result_zero cmp si, infinite_exponent jge @div_underflow sub ax, si cmp ax, infinite_exponent jge @div_overflow cmp ax, zero_exponent jg @ordinary_division @div_underflow: mov ch, fe_Underflow call e287_Exception @div_result_zero: mov si, zero_exponent jmp @div_extreme_result @div_overflow: mov ch, fe_Overflow jmp @div_infinite_excep @div_by_zero: mov ch, fe_Divide_By_Zero @div_infinite_excep: call e287_Exception @div_infinite_result: mov si, infinite_exponent @div_extreme_result: sub ax, ax mov di, [bp][quoP] stosw stosw stosw stosw xchg ax, si stosw mov ax, positive_sign (* zero/infinity always positive*) stosw jmp near @div_emu_end (* here begin the non-extreme cases*) @ordinary_division: mov [bp][quoExp], ax mov al, [bx][emu_sign] xor al, [di][emu_sign] mov [bp][quoSign], al sub cx, cx mov dx, [bp][bp][temp][w3] mov si, [di][emu_fraction][w3] cmp si, dx ja div_noExcess sub dx, si inc cx div_noExcess: mov ax, [bp][temp][w2] div si adc ax, 1 (* round up always !*) adc cl, ch mov [bp][excess], cl mov [bp][bp][result][w3], ax push ax (* Now we have V = excess:[bp][result][w3], up to 18 bits in size. Next*) (* calculate V * (pqrs)*) mov ax, [di][emu_fraction][w0] mov bx, [di][emu_fraction][w1] jcxz w3excessSubtracted mov dx, [di][emu_fraction][w2] w3subExcess: sub [bp][temp][w0], ax (* anticipate calculation*) sbb [bp][temp][w1], bx (* of (abcd) - V * (pqrs)*) sbb [bp][temp][w2], dx (* by subtracting*) sbb [bp][temp][w3], si (* excess * (pqrs)*) loop w3subExcess (* V <= 2:0000H*) (* result may be negative*) w3excessSubtracted: (* Calculate [bp][result][w3] (= si) * (pqrs)*) pop si mul si mov [bp][guard], ax xchg ax, dx xchg ax, bx mul si add bx, ax adc cx, dx (* cx = 0 since excessSubtracted*) mov ax, [di][emu_fraction][w2] mul si add cx, ax adc dx, 0 xchg ax, dx xchg ax, si mul word [di][emu_fraction][w3] add ax, si adc dx, 0 (* We know that V * (pqrs>= (abcd) since V = Roundup ( (abcd) / p )*) (* Therefore, subtract abcd from (V * (pqrs)) so that we can work*) (* with a conveniently positive remainder.*) sub bx, [bp][temp][w0] sbb cx, [bp][temp][w1] sbb ax, [bp][temp][w2] sbb dx, [bp][temp][w3] mov [bp][temp][w0], bx mov [bp][temp][w1], cx mov [bp][temp][w2], ax (* mov [bp][temp][w3], dx not needed, use only dx copy*) (* Now the next stage begins. Calculate W = Roundup (temp / p).*) mov si, [di][emu_fraction][w3] div si xchg cx, ax div si sub dx, dx stc adc ax, dx (* rounded up always !*) adc cx, dx (* cx:ax = W*) sub dx, ax mov [bp][result][w2], dx (* result = V - W*) sbb [bp][result][w3], cx sbb byte [bp][excess], 0 (* The error from the previous stage was up to 3NNN, so W <= 6:0000H.*) (* However, W is typically < 3:0000, so we will get the fastest times*) (* by using subtraction loops to do the effect of*) (* temp - W[w1] * pqrs*) push ax (* save W[w0]*) mov ax, [di][emu_fraction][w0] mov bx, [di][emu_fraction][w1] jcxz w2excessSubtracted mov dx, [di][emu_fraction][w2] w2subExcess: sub [bp][guard], ax sbb [bp][temp][w0], bx sbb [bp][temp][w1], dx sbb [bp][temp][w2], si (* sbb [bp][temp][w3], 0 ignore [bp][temp][w3], it will always be 0*) loop w2subExcess (* Now calculate W[w0] * pqrs.*) w2excessSubtracted: pop si (* retrieve W[w0]*) mul si xchg ax, dx (* forget contents of ax, too insignificant*) xchg ax, bx mul si add bx, ax adc cx, dx (* ASSERT cx = 0*) mov ax, [di][emu_fraction][w2] mul si add cx, ax adc dx, 0 xchg ax, dx xchg si, ax mul word [di][emu_fraction][w3] add ax, si adc dx, 0 (* Again, we know that temp < W[w0] * pqrs, and we want to continue*) (* with a positive remainder, so temp := W[w0] * pqrs - temp*) (* The approximation is now accurate enough that [bp][temp][w3] is always 0.*) sub bx, [bp][guard] sbb cx, [bp][temp][w0] sbb ax, [bp][temp][w1] sbb dx, [bp][temp][w2] mov [bp][guard], bx mov [bp][temp][w0], cx mov [bp][temp][w1], ax (* mov [bp][temp][w2], dx not needed, use only the dx copy*) (* Next we calculate X = RoundUp ( temp / p );*) mov si, [di][emu_fraction][w3] div si xchg cx, ax div si sub dx, dx stc adc ax, dx (* round up always !*) adc cx, dx (* cx:ax = X*) mov [bp][result][w1], ax add [bp][result][w2], cx (* result = V - W + X*) adc [bp][result][w3], dx adc [bp][excess], dl (* Again, V <= 6:0000H but typically V < 3:0000, so we will*) (* use a subtraction loop to do the effect of*) (* temp - V[w1] * pqrs*) (* A new feature at this stage is that we ignore the least word of*) (* pqrs, because it is 16 bits below significance.*) xchg si, ax mov dx, [di][emu_fraction][w1] xchg ax, dx mov bx, [di][emu_fraction][w2] jcxz w1excessSubtracted w1subExcess: sub [bp][guard], ax sbb [bp][temp][w0], bx sbb [bp][temp][w1], dx (* sbb [bp][temp][w2], 0 ignore [bp][temp][w2], it will always be 0*) loop w1subExcess (* Now calculate V[w0] * pqr.*) w1excessSubtracted: mul si xchg ax, dx (* forget contents of ax, too insignificant*) xchg ax, bx mul si add bx, ax adc cx, dx (* ASSERT cx = 0 from excessSubtracted*) mov ax, [di][emu_fraction][w3] mul si add ax, cx adc dx, 0 (* As usual, calculate temp := V[w0] * pqr - temp.*) (* Temp[w2] is now always zero.*) sub bx, [bp][guard] sbb ax, [bp][temp][w0] sbb dx, [bp][temp][w1] mov [bp][guard], bx mov [bp][temp][w0], ax (* mov [bp][temp][w1], dx not needed, use only the dx copy*) (* Lastly we calculate Y = temp / p, rounded to nearest.*) mov si, [di][emu_fraction][w2] add si, si mov si, [di][emu_fraction][w3] adc si, 0 jnc calcY (* if carry occured, then P rounded to 1:0000, so Y = temp. Rearrange*) (* the registers to look as if division has occured:*) xchg bx, dx (* temp = Y = bx:ax:DH*) jmp finalCombine calcY: shl bx, 1 rcl ax, 1 rcl dx, 1 div si xchg bx, ax div si shr bx, 1 rcr ax, 1 rcr dh, 1 (* least 7 bits of DH are don't cares*) finalCombine: mov cx, -1 mov si, cx mov dl, cl not bx not ax neg dh cmc (* negate Y*) adc ax, 0 adc bx, [bp][result][w1] adc cx, [bp][result][w2] (* result = V - W + X - Y*) adc si, [bp][result][w3] adc dl, [bp][excess] (* result is now in DL:si:cx:bx:ax:DH*) div_align: shr dl, 1 jnc div_round rcr si, 1 rcr cx, 1 rcr bx, 1 rcr ax, 1 rcr dh, 1 inc word [bp][quoExp] div_round: sub di, di add dh, dh adc ax, di adc bx, di adc cx, di adc si, di adc dl, 0 mov dx, [bp][quoExp] jz div_checkExp mov si, 8000H (* result must be 8000:0:0:0H*) inc dx div_checkExp: cmp dx, infinite_exponent jl div_store jmp near @div_infinite_result div_store: mov di, [bp][quoP] mov [di][emu_fraction][w0], ax mov [di][emu_fraction][w1], bx mov [di][emu_fraction][w2], cx mov [di][emu_fraction][w3], si mov [di][emu_exponent], dx mov al, [bp][quoSign] mov [di][emu_sign], al @div_emu_end: pop di; pop si; mov sp,bp; pop bp ret div_paramsize (*INCLUDE E287EXPS.ASM*) (* CONST*) (* { coeffs for range -Ln2/2 .. +Ln2/2, magnified twice. }*) (* exp_polycount = 13;*) (* exp_polycoeffs = ARRAY [0..[14] OF emu_temp;*) exp_polycount : dw 13 exp_polycoeffs : dw 0, 0, 0, 0, 0; db positive_sign, 0 dw 00003H, 00000H, 00000H, 4000H, 0; db positive_sign, 0 dw 0AAA9H, 0AAAAH, 0AAAAH, 0AAAH, 0; db positive_sign, 0 dw 05578H, 05555H, 05555H, 0155H, 0; db positive_sign, 0 dw 02326H, 02222H, 02222H, 0022H, 0; db positive_sign, 0 dw 02B1CH, 082D8H, 0D82DH, 0002H, 0; db positive_sign, 0 dw 0F9FCH, 04033H, 03403H, 0000H, 0; db positive_sign, 0 dw 05214H, 03403H, 00340H, 0000H, 0; db positive_sign, 0 dw 06C1EH, 03BC7H, 0002EH, 0000H, 0; db positive_sign, 0 dw 0B04CH, 04FC9H, 00002H, 0000H, 0; db positive_sign, 0 dw 00E91H, 01AE6H, 00000H, 0000H, 0; db positive_sign, 0 dw 07651H, 0011FH, 00000H, 0000H, 0; db positive_sign, 0 dw 02CC5H, 0000BH, 00000H, 0000H, 0; db positive_sign, 0 (* FUNCTION emu_2XM1 ( x : emu_temp ) : emu_temp;*) (* BEGIN*) (* x := x * LogEof2;*) (* emu_2XM1 := x * Polynomial ( FixedPoint (x), exp_polycount, exp_polycoeffs);*) (* END;*) (*public*) emu_2XM1 : (*NEAR*) push bp; mov bp,sp; push si; push di sub sp, emu_temp_size mov di, sp call emu_LogEof2 push di push si push si call emu_Multiply add sp, emu_temp_size (* discard LogEof2*) mov di, [si][emu_exponent] cmp di, -32 jle e@end (* linear domain, answer = x * LogEof2*) mov ax, [si][emu_fraction][w0] mov bx, [si][emu_fraction][w1] mov cx, [si][emu_fraction][w2] mov dx, [si][emu_fraction][w3] e@notLinear: inc di jnl e@fixed (* fix point for Poly*) e@fix: shr dx, 1 rcr cx, 1 rcr bx, 1 rcr ax, 1 inc di jl e@fix adc ax, 0 (* round after shift*) adc bx, 0 adc cx, 0 adc dx, 0 e@fixed: mov di, positive_sign push di sub di, di (* exponent has been fixed at zero*) push di push dx push cx push bx push ax (* fixed ?x now tos*) push cs:[exp_polycount] push cs mov ax, (*offset*) exp_polycoeffs push ax call emup_Polynomial mov ax, sp push ax push si push si call emu_Multiply add sp, emu_temp_size (* discard polynomial result*) e@end: pop di; pop si; pop bp ret 0 (*INCLUDE E287LOGS.ASM*) (* CONST*) (* invLn2 = 1 / Ln (2);*) (* polycount = 9;*) (* log_polycoeffs = ARRAY [0..(log_polycount - [1)] OF emu_temp;*) log_polycount : dw 9 log_polycoeffs : dw 0,0,0,0, 0; db positive_sign, 0 dw 05568H, 05555H, 05555H, 0555H, 0; db positive_sign, 0 dw 034BAH, 03333H, 03333H, 033H, 0; db positive_sign, 0 dw 0C3A7H, 09248H, 04924H, 02H, 0; db positive_sign, 0 dw 05D4DH, 0C722H, 01C71H, 0H, 0; db positive_sign, 0 dw 05624H, 05CEBH, 0174H, 0H, 0; db positive_sign, 0 dw 0AD39H, 0B1ECH, 013H, 0H, 0; db positive_sign, 0 dw 0D6FDH, 00F80H, 01H, 0H, 0; db positive_sign, 0 dw 07AB5H, 010E4H, 0H, 0H, 0; db positive_sign, 0 (* FUNCTION emu_YL2X ( y, x : emu_temp ) : emu_temp;*) (* VAR*) (* frac : emu_temp;*) (* BEGIN*) (* frac := x.fraction;*) (* { We convert frac to a domain symmetric around 0, as required for the domain*) (* of LnXP1. }*) (* IF frac <= Sqrt (1/2) THEN*) (* frac := (2 * frac) - 1;*) (* x.exponent := x.exponent - 1;*) (* ELSE*) (* frac := frac - 1;*) (* emu_YL2X := ( LnXP1 (frac) * invLn2 + x.exponent ) * Y;*) (* END;*) (*public*) emu_YL2X : (*NEAR*) _exponent = -2 _yPtr = -4 _xPtr = -6 push bp mov bp, sp push [di][emu_exponent] push si push di mov ax, [di][emu_fraction][w0] mov bx, [di][emu_fraction][w1] mov cx, [di][emu_fraction][w2] mov dx, [di][emu_fraction][w3] cmp dx, 0B505H (* approx. Sqrt (1/2)*) ja x@else shl ax, 1 rcl bx, 1 rcl cx, 1 rcl dx, 1 (* ignoring carry achieves subtraction of 1*) mov si, positive_sign dec word [bp][_exponent] (* frac * 2^exp = 2 * frac * 2^(exp-1)*) jmp x@endIf x@else: not dx not cx not bx neg ax cmc adc bx, 0 adc cx, 0 adc dx, 0 mov si, negative_sign x@endIf: sub di, di (* di will be exponent for frac, begin = 0*) x@wordNorm: or dx, dx jnz x@bitNorm xchg ax, bx xchg ax, cx xchg ax, dx sub di, 16 cmp di, -64 jg x@wordNorm (* drop through to here if fraction is very tiny*) sub sp, emu_temp_size mov di, sp call emu_0 (* if frac <= 2^-64 then LnXP1 = 0*) jmp x@logDone (* so take a -cut*) x@bitNorm: js x@normalised x@bitNormLoop: dec di shl ax, 1 rcl bx, 1 rcl cx, 1 adc dx, dx jns x@bitNormLoop x@normalised: push si (* sign*) push di (* exponent*) push dx push cx push bx push ax (* pushed frac*) mov si, sp call LnXP1 (* result at [si], which is stack*) mov ax, (*offset*) invLn2 push ax push si push si call emu_Multiply (* tos := LnXP1 (frac) * invLn2*) x@logDone: lea si, [bp][_exponent] sub sp, emu_temp_size (* temp on stack*) mov di, sp call emu_Float lea si, [di][emu_temp_size] push si push di push di call emu_Add (* (temp) := LnXP1 (frac) * invLn2 + x.exponent*) push di push [bp][_yPtr] push [bp][_yPtr] (* result overwrites _y*) call emu_Multiply (* emu_YL2X := temp * Y*) add sp, 2 * emu_temp_size (* clear operands from stack*) x@end: pop di pop si mov sp, bp pop bp ret 0 (* FUNCTION emu_YL2XP1 ( y, x : emu_temp ) : emu_temp;*) (* X is limited by 8087/287 rules to 0 <= x < (1 - Sqrt (1/2))*) (* which is comfortably within (1 - Sqrt (2)) <= x <= (Sqrt(2) - 1),*) (* the domain of LnXP1.*) (* BEGIN*) (* emu_YL2XP1 := LnXP1 (x) * invLn2 * Y;*) (* END;*) (*public*) emu_YL2XP1 : (*NEAR*) (* parameters arranged same as YLn2X*) p1_yPtr = -2 p1_xPtr = -4 push bp mov bp, sp push si push di mov si, di call emup_Push (* Do not overwrite _x -- the 8087 doesn't.*) mov si, sp call LnXP1 (* parameter at [si], result overwrites param.*) mov ax, (*offset*) invLn2 push ax push si push si call emu_Multiply push si push [bp][p1_yPtr] push [bp][p1_yPtr] (* result overwrites _y*) call emu_Multiply add sp, emu_temp_size (* pop temporary from stack*) xp1@end: pop di pop si mov sp, bp pop bp ret 0 (* FUNCTION LnXP1 ( frac : emu_temp ) : emu_temp;*) (* Input domain is (1 - Sqrt(2)) <= frac <= (Sqrt(2) - 1)*) (* We first transform frac into a more linear domain*) (* which is easier to use with polynomials.*) (* ln (frac) := 2 * y * (1 + y^2/3 + y^4/5 + ...)*) (* The polynomial is prepared for * 4 scaling.*) (* VAR*) (* interval, step : cardinal;*) (* BEGIN*) (* y := frac / (2 + frac);*) (* LnXP1 := 2 * y * Poly (Square_Fix (4 * y), log_polycount, log_polycoeffs);*) (* END;*) LnXP1 : push bp mov bp, sp push si push di cmp word [si][emu_exponent], -32 (* in linear domain of small fracs,*) jle LnXP1@end (* LnXP1 (frac) = frac*) sub sp, emu_temp_size mov di, sp call emu_1 inc word ss : [di][emu_exponent] (* tos now = 2.0*) push si push di push di call emu_Add push si push di push di call emu_Divide (* [di] = tos = y = frac / (2 + frac)*) call emu_Duplicate_tos mov bx, sp add word ss : [bx][emu_exponent], 2 (* * 4*) call emup_Square_Fix (* uses and overwrites tos*) push cs:[log_polycount] push cs mov ax, (*offset*) log_polycoeffs push ax call emup_Polynomial (* Poly (Sqfx (y * 4), count, coeffs*) (* uses and overwrites tos*) mov ax, sp (* push sp is incompatible 8086/80286 !!*) push ax push di (* points to duplicate of y*) push si (* result overwrites _x*) call emu_Multiply inc word ss : [si][emu_exponent] (* * 2*) add sp, 2 * emu_temp_size (* clear intermediate values from stack*) LnXP1@end: pop di pop si mov sp, bp pop bp ret 0 arith_Add: push ax push bx push dx push cx (* return address = arith_end*) jmp near emu_Add arith_Multiply: push ax push bx push dx push cx (* return address = arith_end*) jmp near emu_Multiply arith_CompareP: add si, emu_temp_size arith_Compare: push ax push bx push cx (* return address = arith_end*) jmp near emu_Compare arith_SubtractR: xchg ax, bx arith_Subtract: push ax push bx push dx push cx (* return address = arith_end*) jmp near emu_Subtract arith_DivideR: xchg ax, bx arith_Divide: push ax push bx push dx call emu_Divide arith_end: add sp, si (* popAfter*) jmp near e287_Exit (* format name op number modr/m*) (* - - - - - - - - - - -*) (* load / store typ.1 mod.0.S.P. r/m .*) (* IF (CH & 20H) = 0 THEN*) (* IF (CH.S) = 0 THEN*) (* IF CH.P = 0 THEN*) (* CASE CL.typ OF*) (* 0 : iees_Push (es:si);*) (* 1 : emu_Float32 (es:si);*) (* 2 : ieel_Push (es:si);*) (* 3 : emu_Float (es:si);*) (* ELSE*) (* { reserved, logically these are No-Op's }*) (* ELSE*) (* IF (CH.P) = 0 THEN emu_Duplicate_tos;*) (* CASE CL.typ OF*) (* 0 : iees_Pop (es:si);*) (* 1 : emu_Round32 (es:si);*) (* 2 : ieel_Pop (es:si);*) (* 3 : emu_Round (es:si);*) (* ELSE {*) (* - - - - - - - - - - -*) (* -- " -- extended xtp.1 mod.1.S.x. r/m .*) (* IF CH.S = 0 THEN*) (* CASE (x|xtp) OF { load NDP from memory }*) (* 0 : FLDENV;*) (* 1 : { reserved }*) (* 2 : { Restore }*) (* 3 : emu_Push_Decimal (es:si);*) (* 4 : FLDCW;*) (* 5 : emu_Push_Temp87 (es:si);*) (* 6 : { reserved, would be "FLDSW" }*) (* 7 : emu_Float64 (es:si);*) (* ELSE*) (* di := si;*) (* CASE (x|xtp) OF { store NDP to memory }*) (* 0 : FSTENV;*) (* 1 : { reserved };*) (* 2 : { Save }*) (* 3 : emu_Pop_Decimal (es:di);*) (* 4 : [es:[di] := [emws_control];*) (* 5 : emu_Pop_Temp87 (es:di);*) (* 6 : [es:[di] := [emws_status];*) (* 7 : emu_Round64 (es:di);*) (* END;*) ldst_loadCase : dw iees_Push, emu_Float32, ieel_Push, emu_Float ldst_storeCase : dw iees_Pop, emu_Round32, ieel_Pop, emu_Round ldst_notEmulated: ldst_reserved : mov ch, fe_Invalid_Operation call e287_Exception ret 0 (*environment STRUC*) env_control = 0 (*dw ?*) env_status = 2 (*dw ?*) env_tag = 4 (*dw ?*) env_instrnPtr = 6 (*dd ?*) env_dataPtr = 10 (*dd ?*) environment_size = 14 (*environment ENDS*) ldst_ldEnv : mov cl, 4 mov ax, es:[si][env_status] mov [emws_status], ax (* instrnPtr not supported mov ax, env_instrnPtr+w0 (*!!!???*) (*george*) mov [emws_instrnPtr][w0], ax mov ax, env_instrnPtr+w1 (*!!!???*) mov bx, 7FFH and bx, ax xor ax, bx xchg bh, bl (* emws_instruction bytes are swapped*) mov [emws_instruction], bx (* for speed of normal operations*) rol ax, cl mov [emws_instrnPtr][w1], ax mov ax, env_dataPtr+w0 mov [emws_dataPtr][w0], ax mov ax, env_dataPtr+w1 rol ax, cl mov [emws_dataPtr][w1], ax *) (* jmp ldst_loadControl*) ldst_loadControl : (* note keep ldEnv just above here !!!*) mov ax, es : [si][env_control] mov [emws_control], ax mov ch, 0 jmp near e287_Exception (* any exceptions newly unmasked ?*) (* Exception returns to our caller.*) ldst_Float64 : pop dx (* return address*) sub sp, emu_temp_size (* make space on stack*) mov di, sp push dx jmp near emu_Float64 (* will return to this routine's caller*) ldst_Push_Decimal : pop dx (* return address*) sub sp, emu_temp_size (* make space on stack*) mov di, sp push dx jmp near emu_Push_Decimal (* will return to this routine's caller*) ldst_Push_Temp87 : pop dx (* return address*) sub sp, emu_temp_size (* make space on stack*) mov di, sp push dx jmp near emu_Push_Temp87 (* will return to this routine's caller*) ldst_xPushCase : dw ldst_ldEnv, ldst_reserved dw ldst_Restore, ldst_Push_Decimal dw ldst_loadControl, ldst_Push_Temp87 dw ldst_reserved, ldst_Float64 ldst_stEnv : call ldst_storeControl call ldst_storeStatus mov ax, (*offset*) emws_initialSP sub ax, sp mov dl, emu_temp_size div dl (* AL = count of regs in use*) mov cx, -1 xchg cx, ax add cx, cx (* 2 bits per tag*) shr ax, cl (* simulate the tag word*) stosw (*rest not supported mov ax, 16 mul word [emws_instrnPtr][w1] add ax, [emws_instrnPtr][w0] adc dl, dh stosw mov cl, 4 ror dx, cl mov ax, [emws_instruction] xchg ah, al (* bytes are swapped as Emu287 does decode*) or ax, dx stosw mov ax, 16 mul word [emws_dataPtr][w1] add ax, [emws_dataPtr][w0] adc dl, dh stosw ror dx, cl xchg ax, dx stosw *) add di,8 (* Skip *) ret 0 ldst_Save : call ldst_stEnv mov si, sp add si, 2 mov bp, es:[di][-10][w0] (* retrieve the tags*) mov cx, 8 or bp, bp jz ldst_sv_Full mov cl, -1 decLoop: inc cl shl bp, 1 shl bp, 1 jnc decLoop ldst_sv_Full: mov bp, 8 sub bp, cx jcxz ldst_sv_loop2 (* The first loop copies non-empty registers.*) ldst_sv_loop: push cx call emu_Pop_Temp87 pop cx add si, emu_temp_size add di, 10 loop ldst_sv_loop (* The second loop writes out zeroes in place of unused registers.*) or bp, bp jz ldst_sv_done sub ax, ax ldst_sv_loop2: mov cx, 5 rep; stosw dec bp jnz ldst_sv_loop2 ldst_sv_done: ret 0 ldst_Restore : pop [emws_tos] (* remember return address*) (* in a convenient location*) mov sp, (*offset*) emws_initialSP (* stack wiped !!*) mov bp, es:[si][env_tag] add si, (7 * 10) + environment_size (* last stack element*) or bp, 5555H (* The first loop skips empty registers.*) ldst_rs_loop: shr bp, 1 shr bp, 1 jnc ldst_rs_loop2 sub si, 10 or bp, bp jnz ldst_rs_loop jmp ldst_rs_allPushed (* The second loop reads in the defined registers.*) ldst_rs_moreRegs: shr bp, 1 ldst_rs_loop2: sub sp, emu_temp_size mov di, sp call emu_Push_Temp87 (* from es:[si] to ds:[di]*) sub si, 10 shr bp, 1 jc ldst_rs_moreRegs (* when reach here si is pointing one temp87 below the register save vector.*) ldst_rs_allPushed: mov ax, sp xchg ax, [emws_tos] (* fix the stack depth*) push ax (* and restore return address*) add si, (10) - environment_size (* si = & environment*) jmp near ldst_ldEnv (* will return to our caller*) ldst_storeControl : mov ax, [emws_control] stosw ret 0 ldst_storeStatus : mov ax, (*offset*) emws_initialSP sub ax, sp mov cl, emu_temp_size div cl neg al and al, 7 mov cl, 3 shl al, cl mov cx, [emws_status] and ch, 0C7H or ch, al xchg ax, cx stosw ret 0 ldst_Round64 : mov si, sp add si, 2 call emu_Round64 ret emu_temp_size ldst_Pop_Decimal : mov si, sp add si, 2 call emu_Pop_Decimal ret emu_temp_size ldst_Pop_Temp87 : mov si, sp add si, 2 call emu_Pop_Temp87 ret emu_temp_size ldst_xStoreCase : dw ldst_stEnv, ldst_reserved dw ldst_Save, ldst_Pop_Decimal dw ldst_storeControl, ldst_Pop_Temp87 dw ldst_storeStatus, ldst_Round64 (* format name op number modr/m*) (* - - - - - - - - - - -*) (* register moves opa.1 1 1 0.opc. reg .*) (* CASE opc|opa OF*) (* 0 : emu_Push; { FLD reg }*) (* 1, 5, 9, 13 : reserved;*) (* 2 : FFREE;*) (* 3 : FFREEP;*) (* 4, 6, 7 : emu_Push; emu_Swap_tos; emu_Pop; { FXCHG tos, reg }*) (* 8 : IF reg = 0 THEN FNOP*) (* ELSE reserved;*) (* 10 : emu_Duplicate_tos; emu_Pop; { FST reg }*) (* 11, 12, 14, 15 : emu_Pop; { FSTP reg }*) (* END;*) remo_case : dw remo_push, remo_res, remo_free, remo_freep dw remo_xch, remo_res, remo_xch, remo_xch dw remo_nop, remo_res, remo_st, remo_stp dw remo_stp, remo_res, remo_stp, remo_stp remo_push: sub sp, emu_temp_size mov di, sp rep; movsw jmp remo_end remo_freep: mov di, emu_temp_size jmp remo_freeMayPop remo_free: xchg di, si mov ax, infinite_exponent call emu_Extreme jmp remo_end remo_freeMayPop: xchg di, si mov ax, infinite_exponent call emu_Extreme add sp, di jmp remo_end remo_notEm: remo_res: mov ch, fe_Invalid_Operation call e287_Exception jmp remo_end remo_xch: mov di, sp remo_xchgLoop: mov ax, [di] xchg ax, [si] stosw add si, 2 loop remo_xchgLoop jmp remo_end remo_nop: test ch, 7 jnz remo_res jmp remo_end remo_st: mov di, si mov si, sp rep; movsw jmp remo_end remo_stp: mov di, si mov si, sp rep; movsw mov sp, si remo_end: jmp near e287_Exit (* format name op number modr/m*) (* - - - - - - - - - - -*) (* tos functions 0 0 1 1 1 1. fun .*) (* CASE CH.fun OF*) (* 0 : [tos][emu_sign] := -[tos][emu_sign] {=FCHS}*) (* 1 : [tos][emu_sign] := positive_sign; {=FABS}*) (* 2..3 : reserved;*) (* 4 : emu_Compare_Zero; {=FTST}*) (* 5 : emu_XAM;*) (* 6..7 : reserved;*) (* 8 : emu_1;*) (* 9 : emu_Log2of10;*) (* 10 : emu_Log2ofE;*) (* 11 : emu_PI;*) (* 12 : emu_Log10of2;*) (* 13 : emu_LogEof2;*) (* 14 : emu_0;*) (* 15 : reserved;*) (* 16 : emu_2XM1;*) (* 17 : emu_YL2X;*) (* 18 : emu_PTan;*) (* 19 : emu_PAtan;*) (* 20 : tosf_Xtract;*) (* 21 : reserved;*) (* 22 : emu_DecStP;*) (* 23 : emu_IncStP;*) (* 24 : emu_PRem;*) (* 25 : emu_YL2XP1;*) (* 26 : emu_Square_Root;*) (* 27 : reserved;*) (* 28 : emu_RndInt;*) (* 29 : emu_Scale;*) (* 30..31 : reserved;*) (* END;*) tosf_pat01 = 0 tosf_pat10 = 2 tosf_pat11 = 4 tosf_pat12 = 6 tosf_pat21 = 8 tosf_pat22 = 10 tosf_pattern : db tosf_pat11,tosf_pat11 (* tosf_Negate, tosf_Absolute*) db tosf_pat11,tosf_pat11 (* tosf_Reserved, tosf_Reserved*) db tosf_pat11 (* emu_Compare_Zero*) db tosf_pat11 (* emu_XAM*) db tosf_pat11,tosf_pat11 (* tosf_Reserved, tosf_Reserved*) db tosf_pat01,tosf_pat01,tosf_pat01 (* emu_1, emu_Log2of10, emu_Log2ofE*) db tosf_pat01,tosf_pat01,tosf_pat01 (* emu_Pi, emu_Log10of2, emu_LogEof2*) db tosf_pat01 (* emu_0*) db tosf_pat11 (* tosf_Reserved*) db tosf_pat11 (* emu_2XM1*) db tosf_pat21 (* emu_YL2X*) db tosf_pat12 (* emu_PTan*) db tosf_pat21 (* emu_PAtan*) db tosf_pat12 (* tosf_Xtract*) db tosf_pat11 (* tosf_Reserved*) db tosf_pat10 (* emu_DecStP*) db tosf_pat01 (* emu_IncStP*) db tosf_pat22 (* emu_PRem*) db tosf_pat21 (* emu_YL2XP1*) db tosf_pat11 (* emu_Square_Root*) db tosf_pat11 (* tosf_Reserved*) db tosf_pat11 (* emu_RndInt*) db tosf_pat22 (* emu_Scale*) db tosf_pat11 (* tosf_Reserved*) db tosf_pat11 (* tosf_Reserved*) tosf_case : dw tosf_Negate, tosf_Absolute dw tosf_Reserved, tosf_Reserved dw emu_Compare_Zero dw emu_XAM dw tosf_Reserved, tosf_Reserved dw emu_1, emu_Log2of10, emu_Log2ofE dw emu_Pi, emu_Log10of2, emu_LogEof2 dw emu_0, tosf_Reserved dw emu_2XM1 dw emu_YL2X dw emu_PTan dw emu_PAtan dw tosf_Xtract dw tosf_Reserved dw tosf_NoOp (* emu_DecStP*) dw tosf_NoOp (* emu_IncStP*) dw emu_PRem dw emu_YL2XP1 dw emu_Square_Root dw tosf_Reserved dw emu_RndInt dw emu_Scale dw tosf_Reserved dw tosf_Reserved tosf_Negate : xor byte [si][emu_sign], negative_sign ret 0 tosf_Absolute : mov byte [si][emu_sign], positive_sign ret 0 tosf_Xtract : push si push di push si push di cld mov cx, emu_temp_size/2 rep; movsw pop di pop si xchg si, di (* swap so [si] = fraction, [di] = exp*) lea si, [si][emu_exponent] cmp word [si][w0], zero_exponent jle toxt_end cmp word [si][w0], infinite_exponent jge toxt_infinite dec word [si][w0] (* e287 binary point aligned*) call emu_Float mov word [si][w0], 1 (* differently to 80287*) toxt_end: pop di pop si ret 0 toxt_infinite: mov word [si][w0], zero_exponent (* infinity.fraction = 0.0*) mov word [di][emu_exponent], 0DH mov byte [di][emu_fraction][w3][by1], 80H (* infinity.exponent = 4000H*) jmp toxt_end tosf_Reserved : mov ch, fe_Invalid_Operation jmp near e287_Exception (* ret*) tosf_NoOp : ret 0 tosf_patternSwitch : dw tosf_in0out1, tosf_in1out0 dw tosf_in1out1, tosf_in1out2, tosf_in2out1 dw tosf_in2out2 (*public*) e287_tos_Function : mov bx, 1FH and bl, ch mov cl, cs:[tosf_pattern][bx] mov ch, 0 shl bx, 1 (* ready for indexing a word table*) mov di, cx jmp near cs:[tosf_patternSwitch][di] tosf_in0out1: sub sp, emu_temp_size mov di, sp (* no input, output pushed to tos*) jmp tosf_viaCase tosf_in1out1: mov si, sp (* input and output are both tos element*) jmp tosf_viaCase tosf_in1out2: mov si, sp (* tos input*) sub sp, emu_temp_size (* push...*) mov di, sp (* tos1 and tos output*) jmp tosf_viaCase tosf_in1out0: (* only used for DecStP*) add sp, emu_temp_size jmp near e287_Exit tosf_in2out1: mov di, sp lea si, [di][emu_temp_size] call cs:[tosf_case][bx] add sp, emu_temp_size jmp near e287_Exit tosf_in2out2: mov di, sp lea si, [di][emu_temp_size] (* jmp tosf_via_case*) tosf_viaCase: call cs:[tosf_case][bx] jmp near e287_Exit (* format name op number modr/m*) (* - - - - - - - - - - -*) (* FPU control 0 1 1 1 1 1. ctl .*) (* CASE CH.ctl OF*) (* 0 : { FENI -- no-op on 80287 };*) (* 1 : { FDISI -- no-op on 80287 };*) (* 2 : FCLEX;*) (* 3 : FINIT;*) (* 4 : { FSETPM };*) (* 5..31 : reserved;*) (* END;*) (*public*) e287_FPU_Control : and ch, 31 cmp ch, 3 (* FINIT ?*) jne fpuc_notInit (* Set infinite values into 9 registers, corresponding to 1 underflow*) (* guard and 8 normal.*) mov sp, (*offset*) emws_initialSP mov di, sp mov ax, infinite_exponent call emu_Extreme (* set the underflow guard*) lea si, [di][emu_temp_size][-2] sub di, 2 std (* reverse string direction*) sub sp, 8 * emu_temp_size mov cx, (8 * emu_temp_size) / 2 rep; movsw (* make 8 infinite (undefined) temps*) mov sp, (*offset*) emws_initialSP mov word [emws_status], 0 mov word [emws_control], 033FH (* precision = temp real,*) (* all exceptions masked.*) jmp fpuc_end fpuc_notInit: cmp ch, 2 (* FCLEX ?*) jne Reserved mov byte [emws_status][by0], 0 (* clear exceptions*) jmp fpuc_end (*------- This label is called from several nearby groups*) e287_Reserved_1: e287_Reserved_2: Reserved: mov ch, fe_Invalid_Operation call e287_Exception (*-------*) fpuc_end: jmp near e287_Exit (* format name op number modr/m*) (* - - - - - - - - - - -*) (* reserved 1 0 1 1 1 1.- - - - -.*) (* reserved;*) (* END;*) (* (*public*) e287_Reserved_1 : jmp Reserved (* format name op number modr/m*) (* - - - - - - - - - - -*) (* reserved 1 1 1 1 1 1.- - - - -.*) (* IF (CH & 1Fh) = 0 THEN { FSTSW ax }*) (* ELSE reserved;*) (* END;*) (*public*) e287_Reserved_2 : test ch, 1FH jnz Reserved mov ax, (*offset*) emws_initialSP sub ax, sp mov cl, emu_temp_size div cl and al, 7 mov cl, 3 shl al, cl mov cx, [emws_status] and ch, 0C7H or ch, al mov [bp][oldAX], cx jmp near e287_Exit *) (*INCLUDE E287MULT][ASM*) (* Multiplication first involves checks for the special cases of zero*) (* and infinity][ If these are encountered, an exception may be raised*) (* and a special value is returned][ Otherwise, a subroutine is called*) (* to perform a fixed point 64-bit multiplication][ This is then*) (* normalised and rounded, and the result and flags set][*) (*public*) emu_Multiply : x_p = 4+nearParam y_p = 2+nearParam z_p = 0+nearParam mul_paramsize = 6 push bp; mov bp,sp; push si; push di cld mov si, [bp][x_p] mov di, [bp][y_p] mov bx, [di][emu_exponent] mov ax, [si][emu_exponent] cmp ax, bx jl @exponents_ordered xchg ax, bx @exponents_ordered: cmp bx, infinite_exponent jge @mult_infinite cmp ax, zero_exponent jle @mult_zero add ax, bx cmp ax, infinite_exponent jge @mult_overflow cmp ax, zero_exponent jg @mult_ordinary @mult_underflow: mov ch, fe_Underflow call e287_Exception @mult_zero: mov bx, zero_exponent jmp @mult_extreme_result @mult_overflow: mov ch, fe_Overflow call e287_Exception @mult_infinite: mov bx, infinite_exponent @mult_extreme_result: sub ax, ax mov di, [bp][z_p] stosw stosw stosw stosw xchg ax, bx stosw mov al, positive_sign stosb jmp @mult_end @mult_ordinary: mov bl, [si][emu_sign] xor bl, [di][emu_sign] push bx (* save result's sign*) push ax (* save result's exponent*) call emup_Push xchg si, di call emup_Push mov bx, sp push bp lea bp, [bx][-nearParam] (* Inner@Mult caller sets BP*) call Inner@Multiply (* result in DX:CX:BX:AX:SI*) pop bp add sp, 2 * emu_temp_size (* Inner routine does not pop parameters*) pop di (* retrieve result's exponent*) or dx, dx (* result must be normalised*) js @round dec di shl si, 1 rcl ax, 1 rcl bx, 1 rcl cx, 1 rcl dx, 1 (* 1 bit shift is the most ever needed*) @round: shl si, 1 adc ax, 0 adc bx, 0 adc cx, 0 adc dx, 0 jnc @ready rcr dx, 1 (* here, result is 8000:0:0:0H*) inc di @ready: mov si, di mov di, [bp][z_p] stosw xchg ax, bx stosw xchg ax, cx stosw xchg ax, dx stosw xchg ax, si stosw pop ax mov ah, 0 stosw @mult_end: pop di; pop si; pop bp ret mul_paramsize (* Fraction multiplication ignores the exponent and sign and returns the*) (* non-aligned 5-word result of fraction multiplication, the 5th word*) (* being only significant in the top two bits, used for rounding][*) (*public*) emup_Fraction_Multiply : ?x = emu_temp_size+nearParam (* emu_temp*) ?y = 0+nearParam (* emu_temp*) ?paramsize = 2 * emu_temp_size ?res = emu_temp_size+nearParam (* result overwrites [bp][?x]*) push bp; mov bp,sp; push si; push di mov al, [bp][?y][emu_sign] xor [bp][bp][?res][emu_sign], al (* multiply the signs together*) call Inner@Multiply (* result in DX:CX:BX:AX:SI*) mov [bp][?res][emu_fraction][w0], ax mov [bp][?res][emu_fraction][w1], bx mov [bp][?res][emu_fraction][w2], cx mov [bp][?res][emu_fraction][w3], dx mov ax, si pop di; pop si; pop bp ret ?paramsize - emu_temp_size (* result on stack and in AX*) Inner@Multiply : (* Algorithm:*) (* let N = 2^32*) (* the problem to multiply abcd * uvwx, where each letter represents 16 bits][*) (* use the transforms:*) (* abcd = (ab + i)N + (cd - iN), where i = Ord (cd >= N/2)*) (* uvwx = (uv + j)N + (wx - jN), j = Ord (wx >= N/2)*) (* The useful aspect of this transform is that the trailing terms*) (* |(cd - iN)| <= N/2, and |(wx - jN)| <= N/2*) (* are limited in magnitude, and the product of those terms <= NN/4*) (* On the other hand, both leading terms*) (* (ab + i) >= N/2, and (uv + j) >= N/2, since both are left-normalised,*) (* and the product of the leading terms is therefore >= NNNN/4][*) (* The consequence of ignoring the lesser terms' mutual product is to cause*) (* a maximum error of 1/NN, and a probable error of 1/4NN, unbiased][ The*) (* gain is to avoid four 16-bit multiplications][*) (* It turns out that the "i" and "j" terms can easily be added to the result,*) (* given a bit more algebraic juggling:*) (* result = (ab+i)N * (uv+j)N + (ab+i)N * (wx-jN) + (uv+j)N * (cd-iN)*) (* = (ab)(uv)NN + abjNN + uviNN + ijNN*) (* (ab)(wx)N - abjNN - ijNN + wxiN*) (* (uv)(cd)N - uviNN - ijNN + cdjN*) (* -----------------------------------------------------------*) (* = (ab)(uv)NN + (ab)(wx)N + (uv)(cd)N + wxiN + cdjN - ijNN*) (* and so, fortunately, the calculation can be done in terms of the*) (* unsigned 16-bit integers a,b,c,d,u,v,w,x, with the terms using i or j*) (* calculated conditionally if i or j <> 0][*) (* A further optimisation relates to the evaluation of polynomials][*) (* This procedure is used for the multiplication of fixed point*) (* coefficients of polynomials, and many of these have 1,2, or even 3*) (* zeroed most significant words][ The method of coefficient evaluation*) (* ensures that the first of the two multiplicands will be the coeffs,*) (* so checks for leading zeroes can look at ?x and ignore ?y][*) (* The accumulation of products begins with the least significant][ As*) (* more significant products are introduced, the least significant bits*) (* of the accumulation are truncated, so that only the most significant*) (* 80 bits are kept][ These are then normalised and rounded to 64 bits][*) (* the multiplication terms are aligned as follows:*) (* |w x -- * iN*) (* + |c d -- * jN*) (* - 1| -- * ijNN*) (* + |d v*) (* + |b x*) (* + d|u*) (* + c|v*) (* + b|w*) (* + a|x*) (* + c u|*) (* + b v|*) (* + a w|*) (* + b u |*) (* + a v |*) (* + a u |*) (* - - - -*) (* (64 bits)*) (* The four registers DI:SI:CX:BX act as a 64-bit accumulator][ However,*) (* since the terms produce results spanning 96 bits, the alignment of the*) (* "accumulator" changes by 16 bits between each group of similarily*) (* aligned terms][ The 6th, least significant word is discarded when no*) (* longer of use][ The 5th word is kept as it is significant for rounding][*) sub di, di mov cx, di mov si, di test byte [bp][?x][emu_fraction][w1][by1], 80H jz @i_skip mov cx, [bp][?y][emu_fraction][w0] (* + wxiN*) mov si, [bp][?y][emu_fraction][w1] @i_skip: test byte [bp][?y][emu_fraction][w1][by1], 80H jz @j_skip add cx, [bp][?x][emu_fraction][w0] (* + cdjN*) adc si, [bp][?x][emu_fraction][w1] adc di, di test byte [bp][?x][emu_fraction][w1][by1], 80H jz @ij_skip dec di (* - ijNN*) @ij_skip: @j_skip: sub bx, bx (*mulAdd [w0], [w2], cx, si, di*) (* + d*v N*) mov ax, [bp][?x][emu_fraction][w0]; mul word [bp][?y][emu_fraction][w2] add cx, ax; adc si, dx; adc di, 0 (*mulAdd [w2], [w0], cx, si, di, @skip1*) (* + b*x N*) mov ax, [bp][?x][emu_fraction][w2]; or ax, ax jz @skip1 mul word [bp][?y][emu_fraction][w0] add cx, ax; adc si, dx; adc di, 0 (*mulAdd [w2], [w1], si, di, bx*) (* + b*w N*) mov ax, [bp][?x][emu_fraction][w2]; mul word [bp][?y][emu_fraction][w1] add si, ax; adc di, dx; adc bx, 0 @skip1: (*mulAdd [w0], [w3], si, di, bx*) (* + d*u N*) mov ax, [bp][?x][emu_fraction][w0]; mul word [bp][?y][emu_fraction][w3] add si, ax; adc di, dx; adc bx, 0 (*mulAdd [w1], [w2], si, di, bx*) (* + c*v N*) mov ax, [bp][?x][emu_fraction][w1]; mul word [bp][?y][emu_fraction][w2] add si, ax; adc di, dx; adc bx, 0 (*mulAdd [w3], [w0], si, di, bx, @skip2*) (* + a*x N*) mov ax, [bp][?x][emu_fraction][w3]; or ax, ax jz @skip2 mul word [bp][?y][emu_fraction][w0] add si, ax; adc di, dx; adc bx, 0 @skip2: sub cx, cx (* discard bits 80 downwards*) push si (* save bits 64][][79 for later rounding*) mov si, cx (*mulAdd [w1], [w3], di, bx, cx*) (* + c*u N*) mov ax, [bp][?x][emu_fraction][w1]; mul word [bp][?y][emu_fraction][w3] add di, ax; adc bx, dx; adc cx, 0 (*mulAdd [w2], [w2], di, bx, cx, @skip3*) (* + b*v NN*) mov ax, [bp][?x][emu_fraction][w2]; or ax, ax jz @skip3 mul word [bp][?y][emu_fraction][w2] add di, ax; adc bx, dx; adc cx, 0 (*mulAdd [w2], [w3], bx, cx, si*) (* + b*u NN*) mov ax, [bp][?x][emu_fraction][w2]; mul word [bp][?y][emu_fraction][w3] add bx, ax; adc cx, dx; adc si, 0 @skip3: (*mulAdd [w3], [w1], di, bx, cx, @skip4*) (* + a*w N*) mov ax, [bp][?x][emu_fraction][w3]; or ax, ax jz @skip4 mul word [bp][?y][emu_fraction][w1] add di, ax; adc bx, dx; adc cx, 0 adc si, 0 (*mulAdd [w3], [w2], bx, cx, si*) (* + a*v NN*) mov ax, [bp][?x][emu_fraction][w3]; mul word [bp][?y][emu_fraction][w2] add bx, ax; adc cx, dx; adc si, 0 (* final product added in as a special case*) mov ax, [bp][?x][emu_fraction][w3] (* + a*u NN*) mul word [bp][?y][emu_fraction][w3] add cx, ax adc si, dx (* never a carry-out at the end, maximum*) @skip4: (* possible is ][FFFF:FFFF:FFFF:FFFFh*) mov dx, si xchg ax, di pop si (* result now in dx:cx:bx:ax:si*) ret 0 (*INCLUDE E287POLY.ASM*) (* algorithm notes:*) (* The polynomial is rearranged in the usual way (due to Horner ?):*) (* poly := (..(c [N] * x + c [N-[1]) * x + c [N-[2])..) * x + c [0];*) (* which minimises the number of calculations needed.*) (* The procedure is a loop of the form:*) (* poly = coeff [terms-[1];*) (* FOR i := terms-2 DOWNTO 0 DO poly := poly * x + coef [i]*) (* result := 1.0 + poly;*) (* The calculations are fixed point. The largest error is that introduced*) (* by the final stage of multiplication and addition, since x is usually*) (* small and multiplication by x therefore shrinks the prior error in*) (* the calculation.*) (*public*) emup_Polynomial : (*NEAR*) _x = 6+nearParam (* emu_temp*) _terms = 4+nearParam _coef_p = 0+nearParam _result = 6+nearParam (* overwrites parameter*) _paramsize = 18 push bp; mov bp,sp; push si; push di push ds push es mov al, emu_temp_size mul byte [bp][_terms] sub ax, emu_temp_size lds si, [bp][_coef_p] add si, ax push ss pop es (* es = ss until @end*) sub sp, emu_temp_size mov di, sp cld mov cx, emu_temp_size rep; movsb (* poly now on Top Of Stack*) mov di, si dec byte [bp][_terms] sub di, emu_temp_size @for: dec byte [bp][_terms] jl @forEnd lea si, [bp][_x] call emu_Push call emup_Fraction_Multiply sub di, emu_temp_size (* ds:di points to coeff [terms - [2]*) push bp mov bp, sp _poly = 2 (* access to tos*) mov bx, [di][emu_fraction][w0] mov cx, [di][emu_fraction][w1] mov dx, [di][emu_fraction][w2] mov si, [di][emu_fraction][w3] (* coeff now in si:dx:cx:bx*) mov al, [di][emu_sign] cmp al, [bp][_poly][emu_sign] je @add @subtract: (* note that the coeff to be subtracted is almost*) (* always larger than the accumulated lesser*) (* terms of the polynomial. Hence subtract*) (* poly from coeff to keep a positive fraction.*) sub bx, [bp][_poly][emu_fraction][w0] sbb cx, [bp][_poly][emu_fraction][w1] sbb dx, [bp][_poly][emu_fraction][w2] sbb si, [bp][_poly][emu_fraction][w3] jnc @storePoly @negate: (* carry flag = borrow when subtracting,*) (* so we have now a negative fraction. This*) (* must be restored to positive.*) not si not dx not cx neg bx (* -x = (not x) + 1*) cmc adc cx, 0 adc dx, 0 adc si, 0 xor al, negative_sign @storePoly: mov [bp][_poly][emu_fraction][w0], bx mov [bp][_poly][emu_fraction][w1], cx mov [bp][_poly][emu_fraction][w2], dx mov [bp][_poly][emu_fraction][w3], si mov [bp][_poly][emu_sign], al pop bp jmp @for @add: add [bp][_poly][emu_fraction][w0], bx adc [bp][_poly][emu_fraction][w1], cx adc [bp][_poly][emu_fraction][w2], dx adc [bp][_poly][emu_fraction][w3], si (* it is a requirement that*) (* the polynomial partial sum*) (* never exceeds 1, so there*) (* is no carry.*) pop bp jmp @for @forEnd: (* pop [bp][_poly] from tos into registers*) pop bx (* w0*) pop cx (* w1*) pop dx (* w2*) pop si (* w3*) pop ax (* ignore exponent*) pop ax (* sign in al*) sub di, di (* convenient zero*) cmp al, negative_sign (* result above or below 1.0 _*) jne @above @below: not si (* res := 1.0 - poly;*) not dx not cx neg bx cmc (* -poly = (not poly) + 1*) adc cx, di adc dx, di adc si, di (* di = 0 is the exponent, conveniently*) jmp @log_polyend @above: stc (* res := 1.0 + poly;*) rcr si, 1 rcr dx, 1 rcr cx, 1 rcr bx, 1 adc bx, di (* round off*) adc cx, di adc dx, di adc si, di inc di (* di = 1 is the exponent*) @log_polyend: mov [bp][_result][emu_exponent], di mov byte [bp][_result][emu_sign], positive_sign mov [bp][_result][emu_fraction][w0], bx mov [bp][_result][emu_fraction][w1], cx mov [bp][_result][emu_fraction][w2], dx mov [bp][_result][emu_fraction][w3], si pop es pop ds pop di; pop si; pop bp ret _paramsize - emu_temp_size (*INCLUDE E287PTAN.ASM*) (* Polynomial coefficients for Chebyshev approx. to Sine (x) / x,*) (* in terms of every second power.*) sine_terms : dw 8 sine_coeffs : dw 0H, 0H, 0H, 0H, 0H; db positive_sign, 0 dw 0AAB6H, 0AAAAH, 0AAAAH, 2AAAH, 0H; db negative_sign, 0 dw 2403H, 2222H, 2222H, 0222H, 0H; db positive_sign, 0 dw 0E4B3H, 0D00H, 00D0H, 0DH, 0H; db negative_sign, 0 dw 086EH, 0C74BH, 2E3BH, 0H, 0H; db positive_sign, 0 dw 40C6H, 9916H, 6BH, 0H, 0H; db negative_sign, 0 dw 450CH, 0B092H, 0H, 0H, 0H; db positive_sign, 0 dw 45D5H, 0D6H, 0H, 0H, 0H; db negative_sign, 0 (* Coeffs converted into hexadecimals, Cosine (-pi/4 .. pi/4)*) cosine_terms : dw 8 cosine_coeffs : dw 0H, 0H, 0H, 0H, 0H; db positive_sign, 0 (* 1.0*) dw 0FF88H, 0FFFFH, 0FFFFH, 07FFFH, 0H; db negative_sign, 0 dw 9AA6H, 0AAAAH, 0AAAAH, 0AAAH, 0H; db positive_sign, 0 dw 0E67FH, 5B04H, 05B0H, 05BH, 0H; db negative_sign, 0 dw 26EFH, 019BH, 0A01AH, 01H, 0H; db positive_sign, 0 dw 0CE1DH, 93DCH, 49FH, 0H, 0H; db negative_sign, 0 dw 0B10FH, 0F74BH, 8H, 0H, 0H; db positive_sign, 0 dw 0D803H, 0C7BH, 0H, 0H, 0H; db negative_sign, 0 (* algorithm notes*) (* Two Chebyshev series expansions are used. The first series calculates*) (* y = Sine (x) for 0 <= x < pi/4*) (* and the second:*) (* y = Cosine (x) for 0 <= x <= pi/4*) (* Tangent is represented as a Sine, Cosine pair left on top of the stack.*) (* BEGIN*) (* fix := Square_Fix (x);*) (* cos := Polynomial (fix, cosine_terms, cosine_coeffs);*) (* sin := x * Polynomial (fix, sine_terms, sine_coeffs);*) (* emu_PTan := (sin, cos);*) (* END;*) (*public*) emu_PTan : (*NEAR*) push bp; mov bp,sp; push si; push di call emup_Push call emup_Square_Fix (* tos = fix*) call emu_Duplicate_tos (* save a copy for sine*) push cs:[cosine_terms] (* number of terms to evaluate*) mov ax, (*offset*) cosine_coeffs push cs push ax call emup_Polynomial (* cosine now on tos*) call emu_Pop (* put cos in its place in result.*) (* "fix" is now tos*) push cs:[sine_terms] (* number of terms to evaluate*) mov ax, (*offset*) sine_coeffs push cs push ax call emup_Polynomial (* tos = polynomial result*) mov ax, sp (* caution : push sp 8086/286 incompatible*) push ax (* * poly (fix)*) push si (* x*) push si (* x := *) call emu_Multiply (* sine now in tan?x*) add sp, emu_temp_size (* discard Sine poly result*) tan@end: pop di; pop si; pop bp ret 0 (*INCLUDE E287REM.ASM*) (* algorithm notes*) (* The method is similar to division, but the process stops when the*) (* numerator is or becomes less than the divisor. Only the least significant*) (* 16 bits of the quotient are retained.*) (* The procedure ignores operand signs and does not guard against 0 divisor,*) (* trusting the caller, so is not suitable for carefree use.*) (* In detail:*) (* We seek to find (Remainder, Quotient) such that*) (* Q * Y + R = X exactly*) (* X = j.x.2^m 1/2 <= x < 1*) (* Y = y.2^n 1/2 <= y < 1*) (* R = j.r.2^p 1/2 <= r < 1*) (* j is either +1 or -1, and m, n, p and Q are integers.*) (* Initially, R = X, Q = 0.*) (* The eventual goal is for 0 <= R < Y if 0 < Y,*) (* ELSE Y < R <= 0.*) (* STAGE 1 :*) (* IF p > n THEN GOTO end; { X = R is the answer }*) (* ELSE*) (* introduce Y' = y.2^p*) (* transform Y' -> Y , maintaining invariant Q * Y' + R = X :*) (* WHILE*) (* (p > n) AND (r < y)*) (* r, p, Q := 2r, p-1, 2Q { A }*) (* (p >= n) AND (r >= y)*) (* r, Q := (r - y), Q + 1 { B }*) (* Transform { A } does not change the value of R or of (Q * Y'), it merely*) (* re-aligns their representation.*) (* Transform { B } increases the magnitude of (Q * Y') and decreases that*) (* of R. We can further prove that R will eventually become less than Y:*) (* If {A} is used, then p is reduced towards n, while r < 2*y*) (* IF {B} is used, then p is constant but since r < 2*y (initial*) (* condition invariant under A and B), r is reduced to less than*) (* half its present magnitude.*) (* Now, either R < Y, or there are a finite number of applications of {A}*) (* possible before r >= y (a fundamental property of floating point*) (* representation). We can set an upper limit on the number of uses of*) (* {A}, N{A} <= (p-n). Therefore no more than (p-n) uses of {A} can*) (* occur.*) (* Similarily, there are a finite number of applications of {B} possible*) (* before R < Y (given R is not infinite), since X < Y * 2^(p+1-n).*) (* Each application of {B} at least halves R from its initial value X,*) (* so N{B} <= (p+1-n).*) (* Since only {A} or {B} is possible at each turn of the WHILE loop, at*) (* most N{A} + N{B} loops are possible, and we therefore have proof*) (* that R < Y obtains within a strictly bounded time.*) (* Implementation note:*) (* It is crucial that there be no rounding error in this algorithm. The*) (* transform {A} can produce an carry of one bit. The implementation*) (* will use the carry flag to give the extra bit of precision. This is*) (* possible because if an carry occurs during {A} we know that r > y,*) (* and also p >= n, so we know that the guard of {B} will be satisfied and*) (* can immediately apply {B}. After {B} we know that r < y, so we are*) (* again safely within the bounds of normal precision.*) (* STAGE 2 : { normalisation }*) (* After stage 1, we have the correct value of R but it is not properly*) (* normalised. Therefore normalise in the usual way, shifiting left until*) (* the most significant bit is in the MS bit of the MS word.*) (* IF r <> 0*) (* WHILE |r| < .5*) (* r, p := 2r, p-1*) (* After normalisation, we are finished.*) (* The 8087/287 delivers the least significant three bits of the remainder,*) (* and scrambles them into the Condition bits of the status word.*) conditions : db 0H, 02H, 40H, 42H, 01H, 03H, 41H, 43H (*public*) emu_PRem : (*NEAR*) ?sign = -2 Q_ = -4 ?incr = -6 dvrP = -8 resP = -10 push bp mov bp, sp push [di][emu_sign][w0] lea sp, [bp][?incr] push si push di mov ax, [di][emu_fraction][w0] (* fetch the numerator*) mov bx, [di][emu_fraction][w1] mov cx, [di][emu_fraction][w2] mov dx, [di][emu_fraction][w3] mov word [bp][Q_], 0 (* initial quotient := 0*) (* NOTE p_ overwrites numerator pointer, which we no longer need.*) mov di, [di][emu_exponent] sub di, [si][emu_exponent] (* compare the exponents*) jge @while jmp @finish (* The WHILE statement is re-ordered to minimise the number of (slow) jumps*) (* needed. It is optimal to place the WHILE following clause {A}.*) @isA: (* assert p >= n, r < y*) dec di jl @overshootA (* quit if was p = n, maintain p now >= n.*) shl word [bp][Q_], 1 shl ax, 1 rcl bx, 1 rcl cx, 1 adc dx, dx jc @isB (* carry implies r >= 1.0 > y*) jns @isA (* optimise if r < .5 <= y*) @while: (* assert p >= n, 1.0 > r >= .5*) cmp dx, [si][emu_fraction][w3] ja @isB jb @isA cmp cx, [si][emu_fraction][w2] ja @isB jb @isA cmp bx, [si][emu_fraction][w1] ja @isB jb @isA cmp ax, [si][emu_fraction][w0] ja @isB jb @isA jmp @isZero @overshootA: inc di (* assert p=n and r < y, therefore finished*) jmp @normalise @isB: (* assert y < r < 2y, p >= n*) inc word [bp][Q_] sub ax, [si][emu_fraction][w0] sbb bx, [si][emu_fraction][w1] sbb cx, [si][emu_fraction][w2] sbb dx, [si][emu_fraction][w3] or di, di jg @isA (* optimise, since r < y, if p > n*) (* ELSE r < y, p = n, so finished.*) @normalise: or dx, dx js @finish jz @check_zero @align_loop: dec di shl ax, 1 rcl bx, 1 rcl cx, 1 adc dx, dx jns @align_loop @finish: add di, [si][emu_exponent] @end: (* result overwrites ?dvr*) push di (* p_*) mov di, [bp][resP] stosw (* ax*) xchg ax, bx stosw (* bx*) xchg ax, cx stosw (* cx*) xchg ax, dx stosw (* dx*) pop ax stosw (* di*) mov al, [bp][?sign] (* result sign = numerator sign*) cbw stosw (* sign*) and byte [emws_status][by1], ~47H (* ggfb - fixed 31/5/88 *) mov di, 7 mov ax, [bp][Q_] cmp byte [bp][?sign], negative_sign jne @scramble neg ax @scramble: and di, ax (* iNDP delivers only three quotient*) mov dl, cs:[conditions][di] (* bits, scrambled into*) or [emws_status][by1], dl (* a condition code*) pop di pop si mov sp, bp pop bp ret 0 @check_zero: xchg dx, cx xchg cx, bx xchg bx, ax sub di, 16 or dx, dx jnz @normalise xchg dx, cx xchg cx, bx sub di, 16 xchg dx, cx sub di, 16 or dx, dx jnz @normalise @isZero: inc word [bp][Q_] (* ggfb added 15/2/90 *) mov di, zero_exponent sub dx, dx sub cx, cx sub bx, bx sub ax, ax jmp @end (*INCLUDE E287SQFX.ASM*) (* Squaring can be almost twice as fast as multiplying two different*) (* numbers together, since most terms are duplicated.*) (* Fixing the point is done at the start, since the number of bit-shifts*) (* necessary will be doubled after the multiplication. It might seem that*) (* this causes a loss of accuracy. However, consider that the least bit*) (* is to be multiplied by the overall fraction. If the fraction is more*) (* than .5 , in other words if there are no shifts anyway, the least bit*) (* will account for up to one bit of the result. If the fraction is less*) (* than .5, then the least bit will count for less than 1/2 bit of error.*) (* If this is combined with rounding after the shift, then the error due*) (* to shifting first is at most 1/4 bit, which isworth sacrificing for the*) (* speed.*) (*public*) emup_Square_Fix : (*NEAR*) sf_x = 0+nearParam sf_result = 0+nearParam (* result overwrites parameter*) push bp; mov bp,sp; push si; push di mov ax, [bp][sf_x][emu_fraction][w0] mov bx, [bp][sf_x][emu_fraction][w1] mov cx, [bp][sf_x][emu_fraction][w2] mov dx, [bp][sf_x][emu_fraction][w3] sub di, di (* hold the rounding bits*) mov si, [bp][sf_x][emu_exponent] cmp si, -16 jg sf@bitShift mov di, ax mov ax, bx mov bx, cx mov cx, dx sub dx, dx add si, 16 sf@bitShift: inc si jg sf@aligned shr dx, 1 rcr cx, 1 rcr bx, 1 rcr ax, 1 rcr di, 1 jmp sf@bitShift sf@aligned: shl di, 1 mov di, 0 (* di is ZERO throughout.*) adc ax, di adc bx, di adc cx, di adc dx, di (* let the factors be (a:b:c:d) ^2*) (* We apply the usual cut to rounding (see emu_Multiply), which*) (* works even better for squares:*) (* Let i = Ord (c > 8000H)*) (* then a:b:c:d = (ab + i)N : (cd - iN)*) (* square = (ab+i)N ^2 + 2 * (ab+i)N * (cd-iN) + (cd-iN)^2*) (* If the final term is ignored, on the grounds that it has a maximum*) (* error of 1/4 least bit (biased positive unfortunately, but the average*) (* value is 1/12 which can be tolerated)*) (* approx = (ab)^2 NN + 2 abiNN + iiNN*) (* 2(ab)(cd)N - 2 abiNN + 2 cdiN - 2 iiNN*) (* -----------------------------------------------------------*) (* = (ab)^2NN + 2(ab)(cd)N + 2 cdiN - iiNN*) (* the multiplication terms are aligned as follows:*) (* | 2c d -- * iN*) (* - 1| -- * iN*) (* + 2|b d*) (* + 2b |c*) (* + 2a |d *) (* + b b |*) (* + 2a c |*) (* + 2a b |*) (* + a a |*) (* - - - - |- - - -*) (* put the shifted value back into the parameter area, clearing the*) (* registers for work space:*) mov [bp][sf_x][emu_fraction][w0], ax mov [bp][sf_x][emu_fraction][w1], bx mov [bp][sf_x][emu_fraction][w2], cx mov [bp][sf_x][emu_fraction][w3], dx (* si:cx:bx will serve as a 48-bit accumulator which scans from right to left*) (* building up the product. It accumulates the sum of partial products of*) (* the same alignment, then shifts the accumulation right and adds in the sum*) (* of the next, more significant group. Since only a small mumber of 32-bit*) (* numbers are being added to the accumulator, si will be a small number and*) (* carry never overflows si. This greatly simplifies matters. Rounding can*) (* be done at an early stage, as soon as all remaining products are to be*) (* more significant than the rounding position.*) mov bx, di mov cx, di mov si, di test byte [bp][sf_x][emu_fraction][w1][by1], 80H jz @sf_i_skip mov bx, [bp][sf_x][emu_fraction][w0] (* + c D*) mov cx, [bp][sf_x][emu_fraction][w1] shl bx, 1 (* * 2*) rcl cx, 1 (* ignore the carry, equivalent*) (* to -iN*) @sf_i_skip: mov ax, [bp][sf_x][emu_fraction][w0] mul word [bp][sf_x][emu_fraction][w2] add bx, ax adc cx, dx adc si, di (* + b d*) add bx, ax adc cx, dx adc si, di (* * 2*) mov bx, cx (* discard bx since below 80 bit level*) mov cx, si (* and move the rest "right" one word*) sub si, si (* to match next group's alignment*) mov ax, [bp][sf_x][emu_fraction][w1] mul word [bp][sf_x][emu_fraction][w2] add bx, ax adc cx, dx adc si, di (* + b c*) add bx, ax adc cx, dx adc si, di (* * 2*) mov ax, [bp][sf_x][emu_fraction][w0] mul word [bp][sf_x][emu_fraction][w3] add bx, ax adc cx, dx adc si, di (* + a d*) add bx, ax adc cx, dx adc si, di (* * 2*) shl bx, 1 (* round off at bit 64*) adc cx, di adc si, di mov bx, cx mov cx, si (* and move the rest "right" one word*) sub si, si (* to match next group's alignment*) mov ax, [bp][sf_x][emu_fraction][w2] mul word ax add bx, ax adc cx, dx adc si, di (* + b b*) mov ax, [bp][sf_x][emu_fraction][w3] mul word [bp][sf_x][emu_fraction][w1] add bx, ax adc cx, dx adc si, di (* + a c*) add bx, ax adc cx, dx adc si, di (* * 2*) (* henceforth di:si:cx:bx = sum*) mov ax, [bp][sf_x][emu_fraction][w3] mul word [bp][sf_x][emu_fraction][w2] add cx, ax adc si, dx adc di, di (* + a b*) (* zero = di, so "zero" is no longer available.*) add cx, ax adc si, dx adc di, 0 (* * 2*) (* final product added in as a special case*) mov ax, [bp][sf_x][emu_fraction][w3] mul ax (* + a a*) add si, ax adc di, dx (* never a carry-out at the end.*) (* max result is .FFFF:FFFF:..h*) (* result now in di:si:cx:bx*) mov [bp][sf_result][emu_fraction][w0], bx mov [bp][sf_result][emu_fraction][w1], cx mov [bp][sf_result][emu_fraction][w2], si mov [bp][sf_result][emu_fraction][w3], di mov byte [bp][sf_result][emu_sign], positive_sign mov word [bp][sf_result][emu_exponent], 0 @sf_end: pop di; pop si; pop bp ret 0 (*INCLUDE E287SQRT.ASM*) (* Algorithm notes.*) (* 25 May 84 -- Preliminary comments.*) (* I was rather partial to the earlier algorithm, but it really didn't*) (* suit the iAPX-86. It required rather too much shifting and ideally*) (* needed two extended shift registers to hold the intermediate values.*) (* Given a machine such as the Z-8000 or AMD-2903 (or, presumably, the*) (* i8087) it was an algorithm of about the same speed as division.*) (* However, for the iAPX-86 there is a quicker way.*) (* It was shown in the algorithm notes for square root that the*) (* formula*) (* x [n] := (y / x [n-[1] + x [n-[1]) / 2*) (* converges rapidly to Sqrt (y), indeed more than doubling the number*) (* of significant bits if x [n-[1] is already a fair estimate. Using*) (* this method has yielded a very fast precision square root,*) (* typically faster than single precision divide.*) (* We will use a precision Square root to provide the estimate*) (* for the long precision version. Applying the above formula once*) (* should double the accuracy. This is about 10% faster than the bit-*) (* shift algorithm was, but requires only 1/3 as much code.*) (*public*) emu_Square_Root : (*NEAR*) push bp; mov bp,sp; push si; push di mov ax, [si][emu_exponent] (* Square roots of zero or infinity are zero or infinity, respectively.*) cmp ax, zero_exponent jng @noChange cmp ax, infinite_exponent jnl @noChange (* Square root only valid on positive numbers*) cmp byte [si][emu_sign], positive_sign je @isPositive mov ch, fe_Imaginary call e287_Exception mov di, si call emu_0 mov word [di][emu_exponent], infinite_exponent (* mark bad result*) jmp @sq_end @isPositive: lea si, [si][by0] call emup_Push call Inner_Square_Root (* x [n-[1] now on tos*) mov di, sp (* [di] := x [n-[1]*) push si push di push si call emu_Divide (* (y / x [n-[1]*) push di push si push si call emu_Add (* + x [n-[1])*) dec word [si][emu_exponent] (* / 2*) add sp, emu_temp_size (* discard x [n-[1]*) @noChange: @sq_end: pop di; pop si; pop bp ret 0 (* Algorithm Notes*) (* The method is simple first-order successive approximation,*) (* based on the refinement:*) (* x' = (x + z/x)/2*) (* where the aim is to generate x = Sqrt (z).*) (* let y = Sqrt (z) exactly, and x = (1 + a) y { approximately }*) (* then x' = y (1 + a + 1 - a + a**2 - a**3) / 2*) (* = y (1 + a*a/2 - a**3/2...)*) (* which, if abs (a) < 1, gives*) (* a' < a*a/2*) (* Actually, the convergence is excellent even when 1 = a:*) (* let z = 1/4, y = 1/2, x = 1, a = 1*) (* x' = (1 + 1/4)/2= 5/8, a' = 1/4*) (* x'' = (5/8 + 8/20)/2 = 41/80, a'' = 1/40*) (* and we can predict a''' < 1/3200, a'''' < 1/(2**24)*) (* Incidentally, that was the worst-case, since the range we*) (* will cover is 1/4 <= z < 1. It converged to better than*) (* 17 bits in 4 iterations. This is important, because we*) (* break the 32-bit process into two stages. The first stage*) (* finds the most significant 16 bits, then one final iteration*) (* extends that to 32 bits.*) (* In order to stop the iteration after the optimal number of*) (* stages, we simply check the difference between x and z/x.*) (* When the difference is one or less, it is time to finish the*) (* calculation to 32 bits.*) Inner_Square_Root : (*NEAR*) s?x = 0+nearParam s?res = 0+nearParam (* result overlays [bp][s?x]*) push bp; mov bp,sp; push di mov bx, [bp][s?x][emu_exponent] mov dx, [bp][s?x][emu_fraction][w3] mov ax, [bp][s?x][emu_fraction][w2] mov cx, [bp][s?x][emu_fraction][w1] (* cx will act as guard bits.*) sar bx, 1 jnc s@aligned inc bx shr dx, 1 rcr ax, 1 rcr cx, 1 s@aligned: cmp dx, 0FFFEH jae s@nearly_one push bx (* save exponent*) push cx (* save guard bits*) mov cx, dx mov bx, ax (* save fraction in cs:bx*) mov di, dx stc rcr di, 1 (* di = (1 + z) / 2, first approximation.*) s@loop: div di dec di cmp di, ax jna s@stage_2 inc di add di, ax rcr di, 1 mov dx, cx mov ax, bx jmp s@loop (* When the input aligned fraction is greater than .FFFE0001 it must be*) (* a special case, since the square root is greater than .FFFF and thus*) (* stage 1 (16 bit precision) would overflow. Fortunately, no 16 bit*) (* iterations are necessary since the estimate (1 + z) / 2 is accurate*) (* to 32 bits within this range.*) s@nearly_one: stc rcr dx, 1 rcr ax, 1 jmp s@end s@stage_2: inc di cmp di, ax jae s@balanced xchg ax, di (* minimise the difference. The remainder*) (* is not alterred by the swap.*) s@balanced: pop cx (* retrieve guard bits*) pop bx (* retrieve exponent*) xchg ax, cx div di sub dx, dx (* dx insignificant, use as convenient zero*) add di, cx rcr di, 1 rcr ax, 1 adc ax, dx (* .. dx is conveniently zero.*) adc dx, di s@end: mov word [bp][s?res][emu_fraction][w0], 0 mov word [bp][s?res][emu_fraction][w1], 0 mov [bp][s?res][emu_fraction][w2], ax mov [bp][s?res][emu_fraction][w3], dx mov [bp][s?res][emu_exponent], bx mov byte [bp][s?res][emu_sign], positive_sign pop di; pop bp ret 0 (*INCLUDE E287STAK.ASM*) (*public*) emu_Push : (*NEAR*) (* push ES:[si].emu_temp*) pop ax push es : [si][emu_sign][w0] push es : [si][emu_exponent] push es : [si][emu_fraction][w3] push es : [si][emu_fraction][w2] push es : [si][emu_fraction][w1] push es : [si][emu_fraction][w0] jmp near ax (*public*) emup_Push : (*NEAR*) (* push ds:[si][emu_temp*) pop ax push [si][emu_sign][w0] push [si][emu_exponent] push [si][emu_fraction][w3] push [si][emu_fraction][w2] push [si][emu_fraction][w1] push [si][emu_fraction][w0] jmp near ax (*public*) emu_Pop : (*NEAR*) (* pop emu_temp to es:[di]*) pop dx (* save the return address*) cld pop ax stosw (* fraction[w0]*) pop ax stosw (* fraction[w1]*) pop ax stosw (* fraction[w2]*) pop ax stosw (* fraction[w3]*) pop ax stosw (* exponent*) pop ax stosb (* sign*) sub di, 11 jmp near dx (* return*) (*public*) emu_Duplicate_tos : (*NEAR*) pop ax (* save the return address*) mov dx, bp (* save old BP*) mov bp, sp push [bp][emu_sign][w0] push [bp][emu_exponent] push [bp][emu_fraction][w3] push [bp][emu_fraction][w2] push [bp][emu_fraction][w1] push [bp][emu_fraction][w0] mov bp, dx jmp near ax (* Swap the two emu_temps currently in tos and tos[1]*) (*public*) emu_Swap_tos : (*NEAR*) pop cx (* save the return address*) mov dx, bp (* save old BP*) mov bp, sp mov al, [bp][emu_sign] xchg al, [bp][emu_temp_size][emu_sign] mov [bp][emu_sign], al mov ax, [bp][emu_exponent] xchg ax, [bp][emu_temp_size][emu_exponent] mov [bp][emu_exponent], ax mov ax, [bp][emu_fraction][w3] xchg ax, [bp][emu_temp_size][emu_fraction][w3] mov [bp][emu_fraction][w3], ax mov ax, [bp][emu_fraction][w2] xchg ax, [bp][emu_temp_size][emu_fraction][w2] mov [bp][emu_fraction][w2], ax mov ax, [bp][emu_fraction][w1] xchg ax, [bp][emu_temp_size][emu_fraction][w1] mov [bp][emu_fraction][w1], ax mov ax, [bp][emu_fraction][w0] xchg ax, [bp][emu_temp_size][emu_fraction][w0] mov [bp][emu_fraction][w0], ax mov bp, dx jmp near cx (*INCLUDE E287TEMP.ASM*) (* input Intel 8087 temp at es:[si]*) (* output emu_temp at ds:[di]*) (* no exceptions.*) (*public*) emu_Push_Temp87 : (*NEAR*) push si push di push ds push es mov ax, ds mov bx, es mov es, ax mov ds, bx std (* reverse directions*) lea si, [si][8][w0] lea di, [di][emu_sign] lodsw xchg bx, ax (* bx = sign bit : exponent*) sub ax, ax shl bx, 1 jz put_zeroExponent rcl ax, 1 stosw (* write the sign*) shr bx, 1 sub bx, temp_exponent_bias xchg ax, bx stosw (* write the exponent*) movsw (* copy 4 fraction words*) movsw movsw movsw put_end: cld (* default is forward strings*) pop es pop ds pop di pop si ret 0 put_zeroExponent: (* interpret denormals as zeroes.*) stosw (* sign = 0, zero always positive*) mov ax, zero_exponent stosw sub ax, ax stosw stosw stosw stosw (* fraction = 0*) jmp put_end (*public*) emu_Pop_Temp87 : (*NEAR*) (* input : Emu temp at ds:[si]*) (* output: Intel 8087 Temp format at es:[di]*) push si push di cld movsw (* fractions are the same*) movsw movsw movsw lodsw mov cx, [si] cmp ax, infinite_exponent jge pot_infinite add ax, temp_exponent_bias jl pot_zero shl ax, 1 shr cx, 1 rcr ax, 1 pot_exit: stosw (* exponent and sign together*) pop di pop si ret 0 pot_infinite: mov bx, 7FFFH (* biased infinite exponent*) jmp pot_extreme pot_zero: sub bx, bx (* biased zero exponent*) pot_extreme: sub ax, ax (* make sure a zero fraction is written*) sub di, 8 stosw stosw stosw stosw xchg ax, bx jmp pot_exit (*INCLUDE E287XCVT.ASM*) (* procedure to float an integer to yield a long temporary real*) (* input is a 32-bit integer at es:[si]*) (* result is an emu_temp real at ds:[di]*) (* No exceptions.*) (*public*) emu_Float32 : (*NEAR*) mov bx, es : [si][w0] mov dx, es : [si][w1] sub ax, ax mov cx, positive_sign or dx, dx jns @fl32_signed @fl32_negative: not dx neg bx sbb dx, -1 (* negate dx:bx*) mov cl, negative_sign @fl32_signed: mov [di][emu_sign][w0], cx mov cx, 16 or dx, dx jnz @fl32_shift xchg bx, dx mov cl, 0 @fl32_shift: or dx, dx jz @fl32_aligned inc cx shr dx, 1 rcr bx, 1 rcr ax, 1 jmp @fl32_shift @fl32_zero: mov cx, zero_exponent (* jmp @fl32_end*) @fl32_aligned: jcxz @fl32_zero @fl32_end: mov [di][emu_exponent], cx mov [di][emu_fraction][w3], bx mov [di][emu_fraction][w2], ax mov [di][emu_fraction][w1], dx (* dx is zero*) mov [di][emu_fraction][w0], dx ret 0 (*public*) emu_Round32 : (*NEAR*) (* Convert a long temporary floating point number to a 32 bit integer*) (* with fraction rounded according to current rounding mode.*) (* Input operand an emu_temp at ds:[si]*) (* Result at es:[di]*) (* Side effects may exit via exception*) mov cx, [si][emu_exponent] cmp cx, 31 jg @ro32_overflow or cx, cx jge @ro32_non_zero (* -.5..+.5 may round to -1, 0, or 1, according to rounding control*) cmp cx, zero_exponent jle @ro32_zero mov bl, e287_roundingControl and bl, [emws_control][by1] (* get rounding control*) add bl, [si][emu_sign][by0] cmp bl, e287_floor + negative_sign je @ro32_one cmp bl, e287_ceiling + positive_sign je @ro32_one @ro32_zero: sub dx, dx mov ax, dx jmp @ro32_endJmp @ro32_one: sub dx, dx mov ax, 1 cmp bl, e287_floor + negative_sign jne @ro32_endJmp neg ax not dx jmp @ro32_endJmp @ro32_overflow: mov ch, fe_Operand_Too_Big call e287_Exception mov dx, 8000H sub ax, ax @ro32_endJmp: jmp near @ro32_end @ro32_non_zero: mov bx, [si][emu_fraction][w1] (* bx will be sticky,*) or bl, [si][emu_fraction][w0][by0] (* used for rounding up*) or bl, [si][emu_fraction][w0][by1] mov ax, [si][emu_fraction][w2] mov dx, [si][emu_fraction][w3] sub cl, 16 ja @ro32_wordAligned @ro32_below_1: or al, bl or al, bh xchg bx, ax (* sticky bx*) xchg ax, dx sub dx, dx add cl, 16 jng @ro32_below_1 (* if 0.5 <= x < 1 go round twice*) and cl, 0FH @ro32_wordAligned: jcxz @ro32_noBitShift push si mov si, 0FFFFH rol dx, cl rol ax, cl (* align the fraction*) shl si, cl mov cx, si (* and two copies of a mask*) and cx, ax xor ax, cx (* mask out the low-end shift-out*) and si, dx xor dx, si or ax, si (* transfer bits from msw to lsw*) or bl, bh or bl, cl (* bl is a "sticky" byte for rounding up.*) mov bh, ch (* bh is sticky + 1/2 ls bit.*) pop si @ro32_noBitShift: mov cl, e287_roundingControl and cl, [emws_control][by1] (* get rounding control*) cmp cl, e287_chop je @ro32_chop cmp cl, e287_round je @ro32_roundNear add cl, [si][emu_sign][by0] cmp cl, e287_floor + positive_sign je @ro32_chop cmp cl, e287_ceiling + negative_sign je @ro32_chop @ro32_roundUp: (* this is what sticky bits are for !*) neg bx (* sets carry if any bit not zero*) jc @ro32_addCarry jmp @ro32_setSign @ro32_roundNear: mov cl, 1 and cl, al or bl, cl (* apply banker's rounding if sticky was 8000H*) add bx, 7FFFH @ro32_addCarry: adc ax, 0 adc dx, 0 jns @ro32_setSign jmp near @ro32_overflow @ro32_chop: @ro32_setSign: cmp byte [si][emu_sign], negative_sign jne @ro32_end not dx neg ax sbb dx, -1 (* negate dx:ax*) @ro32_end: mov es:[di][w0], ax mov es:[di][w1], dx ret 0 (* procedure to float an integer to yield a long temporary real*) (* input is a 64-bit signed integer at es:[si]*) (* result is an emu_temp at ds:[di]*) (* No exceptions.*) (*public*) emu_Float64 : (*NEAR*) push bp push si mov ax, es : [si][w0] mov bx, es : [si][w1] mov cx, es : [si][w2] mov dx, es : [si][w3] mov bp, positive_sign or dx, dx jl @fl64_negative jg @fl64_signed or cx, cx jnz @fl64_signed or bx, bx jnz @fl64_signed or ax, ax jnz @fl64_signed @fl64_zero: mov si, zero_exponent jmp @fl64_end @fl64_negative: not dx not cx not bx not ax add ax, 1 adc bx, 0 adc cx, 0 adc dx, 0 (* negate dx:cx:bx:ax*) mov bp, negative_sign @fl64_signed: mov si, 64 @fl64_maybeTiny: or dx, dx jnz @fl64_notTiny xchg dx, cx xchg cx, bx xchg bx, ax sub si, 16 jmp @fl64_maybeTiny @fl64_notTiny: js @fl64_aligned @fl64_shift: dec si add ax, ax adc bx, bx adc cx, cx adc dx, dx jns @fl64_shift @fl64_aligned: @fl64_end: mov [di][emu_sign][w0], bp mov [di][emu_exponent], si mov [di][emu_fraction][w3], dx mov [di][emu_fraction][w2], cx mov [di][emu_fraction][w1], bx mov [di][emu_fraction][w0], ax pop si pop bp ret 0 (*public*) emu_Round64 : (*NEAR*) (* Convert a long temporary floating point number to a 64 bit integer*) (* with fraction rounded as specified in [emws_control].*) (* Input operand an emu_temp at ds:[si]*) (* Result at es:[di]*) (* May generate an exception.*) push bp push di mov cx, [si][emu_exponent] cmp cx, 63 jg @ro64_overflow or cx, cx jge @ro64_non_zero (* small fractions may round to 1 or -1, according to rounding control*) cmp cx, zero_exponent jle @ro64_zero (* true zero rounds to zero*) mov bl, e287_roundingControl and bl, [emws_control][by1] (* get rounding control*) add bl, [si][emu_sign][by0] cmp bl, e287_floor + negative_sign je @ro64_one cmp bl, e287_ceiling + positive_sign je @ro64_one @ro64_zero: sub bp, bp jmp @ro64_extreme @ro64_one: sub bp, bp mov dx, bp mov ax, 1 cmp bl, e287_floor + negative_sign mov bx, bp jne @ro64_endJmp neg ax not bx not dx not bp jmp @ro64_endJmp @ro64_overflow: mov ch, fe_Operand_Too_Big call e287_Exception mov bp, 8000H @ro64_extreme: sub dx, dx mov bx, dx mov ax, bx @ro64_endJmp: jmp near @ro64_end @ro64_non_zero: mov bp, [si][emu_fraction][w3] mov dx, [si][emu_fraction][w2] mov bx, [si][emu_fraction][w1] mov di, [si][emu_fraction][w0] sub ax, ax (* ax will act as sticky collector*) sub cl, 48 ja @ro64_wordAligned @ro64_wordShift: or al, ah (* AL is completely sticky*) mov ah, 0 or ax, di (* AH is 1/2 ls bit + sticky*) mov di, bx mov bx, dx mov dx, bp sub bp, bp add cl, 16 jng @ro64_wordShift and cl, 0FH @ro64_wordAligned: neg cl jz @ro64_bitAligned add cl, 16 @ro64_bitShift: or al, ah shr bp, 1 rcr dx, 1 rcr bx, 1 rcr di, 1 rcr ah, 1 loop @ro64_bitShift @ro64_bitAligned: mov cl, e287_roundingControl and cl, [emws_control][by1] cmp cl, e287_chop je @ro64_chop cmp cl, e287_round je @ro64_roundNear add cl, [si][emu_sign][by0] cmp cl, e287_floor + positive_sign je @ro64_chop cmp cl, e287_ceiling + negative_sign je @ro64_chop @ro64_roundUp: neg ax (* set carry if any sticky or rounding bit = 1*) jc @ro64_addCarry jmp @ro64_setsign @ro64_overflowJmp: jmp @ro64_overflow @ro64_roundNear: mov cl, 1 and cx, di (* bankers' rounding, add the odd integer bit*) or al, cl (* to sticky, force 8000H to 8001H*) add ax, 7FFFH @ro64_addCarry: mov cx, 0 (* convenient zero*) adc di, cx adc bx, cx adc dx, cx adc bp, cx js @ro64_overflowJmp @ro64_chop: @ro64_setsign: xchg ax, di (* ax is now more convenient than di*) cmp byte [si][emu_sign], negative_sign jne @ro64_end not bp not dx not bx neg ax cmc (* Fixed 12/11/91 *) adc bx, 0 adc dx, 0 adc bp, 0 @ro64_end: pop di mov es:[di][w0], ax mov es:[di][w1], bx mov es:[di][w2], dx mov es:[di][w3], bp pop bp ret 0 (*public*) emu_RndInt : (*NEAR*) (* Round the tos contents to the nearest integer, although the number remains*) (* in emu_temp format on tos.*) (* Input operand at [si], but we change to more convenient [bp]*) (* Output overwrites input.*) (* may report exceptions.*) push bp push si push di xchg bp, si mov cx, [bp][emu_exponent] cmp cx, 64 (* note this is wider range than int64*) jg rndi_overflow or cx, cx jge rndi_normalRange cmp cx, zero_exponent jle rndi_zero mov cl, e287_roundingControl and cl, [emws_control][by1] add cl, [bp][emu_sign][by0] cmp cl, e287_floor + negative_sign je rndi_one cmp cl, e287_ceiling + positive_sign je rndi_one rndi_zero: sub bx, bx mov si, zero_exponent jmp rndi_special rndi_one: mov si, 1 mov bx, 8000H jmp rndi_special rndi_overflow: mov ch, fe_Precision call e287_Exception jmp near rndi_end rndi_special: sub ax, ax mov [bp][emu_fraction][w0], ax mov [bp][emu_fraction][w1], ax mov [bp][emu_fraction][w2], ax mov [bp][emu_fraction][w3], bx mov [bp][emu_exponent], si jmp near rndi_end rndi_normalRange: (* locate the byte where the most significant fraction bits lie.*) mov si, 7 rndi_nextByte: sub cl, 8 jl rndi_split dec si jge rndi_nextByte (* if you fall through here, then there are no fraction bits. The*) (* value in X is effectively an integer.*) jmp near rndi_end rndi_split: (* byte [si] of the fraction has at least one MS bit of the fraction.*) (* load registers AH = least sig all integer byte*) (* AL = split integer/fraction byte*) (* BL = sticky fraction bits*) cmp si, 7 jne rndi_notMSbyte mov ah, 0 mov al, [bp][si][by0] jmp rndi_AXloaded rndi_notMSbyte: mov ax, [bp][si][w0] rndi_AXloaded: sub bx, bx mov di, si rndi_nextSticky: dec di jl rndi_stickyLoaded or bl, [bp][di][by0] mov [bp][di][by0], bh (* clear x.fraction*) jmp rndi_nextSticky rndi_stickyLoaded: mov dx, 0FFH and cl, 07H shr dx, cl (* zeroes = integers, ones = fractions*) mov di, dx inc di (* di = fraction carry / least integer bit*) and dx, ax (* dx = fraction MS bits*) xor ax, dx (* ax = integer LS bits*) or bl, dl jz rndi_maskAligned push cx mov ch, fe_Precision call e287_Exception pop cx rndi_maskAligned: mov cl, e287_roundingControl and cl, [emws_control][by1] cmp cl, e287_chop je rndi_chop cmp cl, e287_round je rndi_roundNear add cl, [bp][emu_sign][by0] cmp cl, e287_floor + positive_sign je rndi_chop cmp cl, e287_ceiling + negative_sign je rndi_chop rndi_roundUp: or bl, dl (* were there any fraction bits ?*) jnz rndi_oneUp jmp rndi_end (* if not, number already is an integer*) rndi_chop: mov [bp][si][by0], al cmp byte [bp][emu_fraction][w3][by1], 0 (* did we chop to zero ?*) jne rndi_end jmp near rndi_zero rndi_roundNear: or bl, bl jnz rndi_bankersUp (* Banker's rounding aims to balance out the rounding of exact halves by*) (* rounding always towards the nearest even integer. Thus half the*) (* halves round up, and half down, on average cancelling.*) rndi_bankersBit: test ax, di (* test for even/odd integer*) jnz rndi_bankersUp rndi_bankersDown: add dx, dx cmp dx, di jna rndi_chop jmp rndi_oneUp rndi_bankersUp: add dx, dx (* was fraction at least half ?*) cmp dx, di (* if so, carry will be found*) jb rndi_chop rndi_oneUp: sub si, 7 jl rndi_severalBytes add ax, di mov [bp][si][7], al neg ah (* set carry if not zero*) jmp rndi_oneUpped rndi_severalBytes: add ax, di mov [bp][si][7], ax inc si rndi_oneUpLoop: inc si jg rndi_oneUpped adc byte [bp][si][7][by0], 0 jc rndi_oneUpLoop rndi_oneUpped: jnc rndi_end (* is there an overflowed carry ?*) stc (* if carry, then result has form 8000:0:0:0H*) rcr word [bp][emu_fraction][w3], 1 inc word [bp][emu_exponent] rndi_end: pop di pop si pop bp ret 0 (*INCLUDE IEELSTAK.ASM*) (* An IEEE long real is packed into 64 bits.*) (* -------- ---- ---- -------- -------- -------- -------- -------- --------*) (* |seeeeeee eeee.ffff ffffffff ffffffff ffffffff ffffffff ffffffff ffffffff|*) (* -------- ----.---- -------- -------- -------- -------- -------- --------*) (* 63 . 0*) (* 1.*) (* The three parts of the number are its sign, exponent, and fraction.*) (* If the number is denoted by the triple (sign, exponent, fraction) then*) (* its value is:*) (* (-1)^sign * 2^(exponent-1023) * 1.fraction*) (* with the exception of the special cases with all-zeroes exponents.*) (* The sign is in the most significant bit and is "1" if negative.*) (* The exponent occupies the next 11 bits. It represents powers of 2, biased*) (* so that an exponent of 3FFH is used for numbers in the range 1.0 to 1.999...*) (* If the exponent is all ones, then the number is either "infinity" or*) (* "Not A Number" (NAN), described below. If the exponent is all zeroes, then*) (* the number is either zero or a "denormal", as described below.*) (* The fraction is the remaining 52 bits. Excepting "denormal" special cases,*) (* numbers are normalised into the triple form shown above. Since the final*) (* term is a number in the range 1.0 to 1.999... it is predictable that the*) (* first bit of the final term shoul be "1" (the "integer bit"), and this*) (* need not be packed into the 64 bits. However, this bit is necessary for*) (* calculation, and so whenever the number is unpacked into temporary format*) (* the integer bit is restored. The 52 bits of packed IEEE fraction are*) (* therefore equivalent to 53 bits of temporary fraction.*) (* If the exponent is all ones and the fraction is all zeroes, then the number*) (* represented is infinite. IEEE makes allowance for both positive and*) (* negative infinities, but both are treated as positive by the Combinatorics*) (* algorithms. If the exponent is all ones but the fraction not zeroes, then*) (* the value is a "NAN", an illegal combination. A huge range of NAN values*) (* is possible, and one suggested use is as distinct space-fill values for use*) (* in detecting uninitialised variables in programs. When a NAN is found*) (* an overflow status is generated and an infinite temporary used.*) (* If the exponent and fraction are all zeroes then the number represented is*) (* zero. Zero can have either sign, but Combinatorics treats either as if*) (* positive. If the exponent is zeroes but the fraction is not, then the*) (* number is a denormal, normally a value with an exponent in the range from*) (* -1024 to -1075. Combinatorics algorithms do not use denormals, instead*) (* flagging them with underflow status and substituting zero. It is an*) (* extremely insignificant range.*) (* convert a temporary floating-point number to an IEEE 64-bit real.*) (* inputs*) (* - a long temporary floating point (emu format) at ds : [si]*) (* outputs*) (* - a long IEEE real at es:[di]*) (* side effects*) (* - instruction may quit in case of exception*) (*public*) ieel_Pop : (*NEAR*) push di (* move fraction to DL:cx:di:ax:DH. The strange choice of registers*) (* is the result of forward planning to minimise register*) (* shuffling and shifting at the pack-and-align stage.*) mov ax, [si][emu_fraction][by1][w0] mov di, [si][emu_fraction][by1][w1] mov cx, [si][emu_fraction][by1][w2] mov dl, [si][emu_fraction][w3][by1] mov bx, [si][emu_exponent] mov dh, 3 and dh, al or dh, [si][emu_fraction][w0][by0] (* sticky rounding..*) jnz @lpop_isSticky test al, 8 (* bankers' rounding to nearest-even*) jz @lpop_rounded (* means alternate round up/down of exact .5*) @lpop_isSticky: add ax, 4 (* round off at the 53rd bit*) adc di, 0 adc cx, 0 adc dl, 0 jnc @lpop_rounded rcr dl, 1 (* carry overflow implies 8000:0:0:0 result*) inc bx @lpop_rounded: add bx, ieel_exponent_bias jle @lpop_small cmp bx, 7FFH jl @lpop_normal_range (* All numbers too large for IEEE range will be given as positive infinity.*) (* If the number was already infinity, then ok, but if it becomes an*) (* infinity through limitations in representation, give overflow status.*) @lpop_large: cmp word [si][emu_exponent], infinite_exponent jge @lpop_infinite mov ch, fe_Overflow call e287_Exception @lpop_infinite: mov dx, 07FF0H (* note clear sign bit = positive*) jmp @lpop_extreme (* All number too small for IEEE normal range will be given as zeroes.*) (* If the number was already zero, this is no problem. However, if the*) (* number is zeroed because it falls outside IEEE range, an underflow*) (* status is reported.*) @lpop_small: cmp word [si][emu_exponent], zero_exponent jle @lpop_zero mov ch, fe_Underflow call e287_Exception @lpop_zero: sub dx, dx @lpop_extreme: sub cx, cx (* both zero and infinity have zero fraction*) mov bx, cx mov ax, bx jmp @lpop_end (* For normal range numbers, we must do some rotations and pack the three*) (* parts into their correct positions.*) @lpop_normal_range: and al, 0F8H (* Discard least bits of fraction.*) shl dl, 1 (* Discard integer bit*) shr bx, 1 rcr dl, 1 or al, bh (* Replace with upper bits of exponent.*) mov dh, bl (* Move lower byte of exponent into place.*) mov bx, di (* di will be consumed to prime the carry flag.*) (* Number is now in dx:cx:bx:ax, but lacks*) (* sign and needs rotation.*) shr di, 1 rcr ax, 1 rcr dx, 1 rcr cx, 1 rcr bx, 1 shr di, 1 rcr ax, 1 rcr dx, 1 rcr cx, 1 rcr bx, 1 or al, [si][emu_sign] (* Put the sign bit into place.*) shr di, 1 rcr ax, 1 rcr dx, 1 rcr cx, 1 rcr bx, 1 @lpop_end: pop di cld (* store dx:cx:bx:ax at es:[di]*) stosw xchg ax, bx stosw xchg ax, cx stosw xchg ax, dx stosw sub di, 8 (* preserve original di value*) ret 0 (* procedure to unpack a long IEEE real yielding a long temporary*) (* inputs*) (* - a long IEEE real in four words at es : [si]*) (* outputs*) (* - a long temporary floating point in Emu287 format at ds:[di]*) (* side effects*) (* - may cause exit via exception for NANs and Denormals*) (*public*) ieel_Push : (*NEAR*) push di (* more convenient to have ds and es the other way around, so swap*) mov ax, ds mov bx, es mov ds, bx mov es, ax mov dx, [si][w3] (* word 3 contains sign,*) sub ax, ax (* exponent, and*) shl dx, 1 (* top 4 bits of fraction*) rcl ax, 1 (* ax is sign*) mov cl, 5 shr dx, cl (* dx now holds only the exponent*) jz @lpush_small cmp dx, 7FFH je @lpush_large @lpush_normal_range: std (* reverse string direction*) lea di, [di][emu_sign] (* result "pushed" in reverse order*) stosw (* push the sign*) sub dx, ieel_exponent_bias xchg ax, dx stosw (* push the exponent*) mov dx, [si][1][w2] mov cx, [si][1][w1] mov bx, [si][1][w0] mov ah, [si][by0] mov al, 0 (* fraction now in dx:cx:bx:ax*) and dh, 0FH (* but needs alignment*) shl ax, 1 rcl bx, 1 rcl cx, 1 rcl dx, 1 shl ax, 1 rcl bx, 1 rcl cx, 1 rcl dx, 1 shl ax, 1 rcl bx, 1 rcl cx, 1 rcl dx, 1 or dh, 80H (* set the integer bit*) xchg ax, dx stosw xchg ax, cx stosw xchg ax, bx stosw xchg ax, dx stosw (* the fraction is now pushed*) @lpush_end: mov ax, ds (* swap segments back as they were*) mov bx, es mov ds, bx mov es, ax pop di ret 0 (* handle extreme values here.*) @lpush_small: mov dx, zero_exponent mov ch, fe_Operand_Too_Small (* anticipate a Denormal*) jmp @lpush_extreme @lpush_large: mov dx, infinite_exponent mov ch, fe_Operand_Too_Big (* anticipate a NAN*) @lpush_extreme: mov al, [si][w3][by0] (* Check for NAN's or Denormals*) and ax, 000FH (* which will not have*) or ax, [si][w2] (* all-zeroes*) or ax, [si][w1] (* fractions.*) or ax, [si][w0] jz @lpush_extremeCorrect call e287_Exception @lpush_extremeCorrect: xchg ax, dx call emu_Extreme jmp @lpush_end (*INCLUDE IEESSTAK.ASM*) (* An IEEE real is packed into 32 bits.*) (* byte 3 byte 0*) (* - ------- - ------- -------- --------*) (* |s eeeeeee|e.fffffff|ffffffff|ffffffff|*) (* - ------- -.------- -------- --------*) (* bit 32 . 0 bit*) (* 1.*) (* The most significant bit is the sign bit, one if negative.*) (* The next 8 bits are the exponent, "biased" so that a value of 7Fh*) (* is used for numbers from 1.0 to 1.999..., and the extreme values*) (* 00H and FFH represent zero and infinity respectively. The exponent*) (* is in powers of 2.*) (* The remaining 23 bits are the normalised fraction. When the fraction*) (* is normalised its leftmost bit is always "1", so it can be ignored*) (* and only the 23 next, unpredictable, bits are kept. The 24th bit is*) (* restored before any calculations. For both infinities and zero the*) (* fraction part is all zeroes.*) (* If the number is denoted by the triple (sign, exponent, fraction) then*) (* its value is*) (* (-1)^sign * 2^(exponent-127) * 1.fraction*) (* There are in addition some special values. NAN's (Not A Number) are*) (* variations on infinity, having an exponent of FFH but fraction parts*) (* which are not all zeroes. NAN's have various uses the most likely of*) (* which is use as memory-fill values when checking for non-initialisation*) (* errors. When a NAN is encountered, an overflow status is generated.*) (* The other special values are denormals, which are underflows. These*) (* have zero exponent but non-zero fraction, while true zero has an*) (* all-zeroes fraction. They represent numbers in the marginal range*) (* too small to be represented normally but close enough to the edge*) (* of normal to allow partial representation. Here no attempt is made*) (* to interpret de-normals, which are of marginal utility, and they*) (* give rise to underflow status if they are encountered.*) (*public*) iees_Pop : (*NEAR*) (* FUNCTION iees_Pop (x : emu_temp) : IEEE_short;*) (* convert a temporary floating-point number to an IEEE 32-bit real.*) (* input:*) (* - a temporary floating point number at ds:[si]*) (* output:*) (* - an IEEE real at es:[di]*) (* side effects*) (* - may exit via exception*) mov ax, [si][emu_fraction][w2] mov dx, [si][emu_fraction][w3] mov bx, [si][emu_exponent] mov cl, [si][emu_sign] cmp word [si][emu_fraction][w0], 0 jnz @spop_isSticky cmp word [si][emu_fraction][w1], 0 jnz @spop_isSticky test al, 7FH jnz @spop_isSticky test ah, 1 (* bankers' rounding causes alternate*) jz @spop_carried (* round up/down of exact .5 fraction.*) @spop_isSticky: add al, al (* round off least 8 bits*) adc ah, 0 adc dx, 0 jnc @spop_carried (* rounding can force renormalisation*) rcr dx, 1 rcr ax, 1 inc bx @spop_carried: add bx, iees_exponent_bias (* IEEE exponent is biased*) jle @spop_underflow cmp bx, 255 jnl @spop_overflow (* If the number is not to be zero or infinite, then pack the three*) (* parts into 32 bits as per IEEE.*) @spop_normal: shl dx, 1 (* throw away the integer bit (always = 1)*) shr cl, 1 (* get the sign bit*) rcr bl, 1 (* which becomes most sig bit of result*) rcr dx, 1 (* .. and exponent shuffles up fraction*) mov al, ah mov ah, dl mov dl, dh mov dh, bl @spop_end: stosw (* store result at es:di*) xchg ax, dx stosw sub di, 4 (* restore original value of di*) ret 0 (* Note that status warning need be generated only if the number has*) (* just become infinite or zero due to change in representation. If*) (* the input was already zero or infinite, the status remains OK.*) @spop_overflow: (* all overflows become positive infinity*) cmp word [si][emu_exponent], infinite_exponent jge @spop_prior_infinity mov ch, fe_Overflow call e287_Exception @spop_prior_infinity: mov dx, 7F80H sub ax, ax jmp @spop_end @spop_underflow: (* all underflows become positive zero*) cmp word [si][emu_exponent], zero_exponent jle @spop_prior_zero mov ch, fe_Underflow call e287_Exception @spop_prior_zero: sub dx, dx mov ax, dx jmp @spop_end (*public*) iees_Push : (*NEAR*) (* FUNCTION iees_Push (x : IEEE_short) : emu_temp;*) (* Procedure to take a IEEE real and push the equivalent*) (* temporary format number onto stack.*) (* inputs*) (* - short IEEE real at es:si*) (* outputs*) (* - temp floating point at ds:di*) (* side effects*) (* - can set exception exit*) push si push di push es cld (* forward string direction*) mov ax, es : [si][w0] mov dx, es : [si][w1] sub si, si shl dx, 1 (* carry := sign, DH := exponent*) rcl si, 1 (* si is sign*) or dh, dh jz @spush_small (* check for extreme exponents*) cmp dh, 0FFH je @spush_large @spush_normal: mov bl, dh (* get the exponent byte*) mov bh, 0 sub bx, iees_exponent_bias (* bx is exponent*) stc (* set carry = integer bit*) rcr dl, 1 mov dh, dl mov dl, ah mov ch, al mov cl, 0 (* dx:cx is now the fraction*) @spush_result: push ds pop es (* arrange so es:[di] points to result*) sub ax, ax stosw (* extended to 64 bits with zeroes*) stosw xchg ax, cx stosw xchg ax, dx stosw xchg ax, bx stosw xchg ax, si stosw @spush_end: pop es pop di pop si ret 0 (* special treatment for extreme values.*) @spush_large: inc dh (* convert FFH to 00H*) mov bx, infinite_exponent mov ch, fe_Operand_Too_Big (* anticipate a NAN*) jmp @spush_extreme @spush_small: mov bx, zero_exponent mov ch, fe_Operand_Too_Small (* anticipate a Denormal*) @spush_extreme: or dx, ax (* check for de-normals*) jz @spush_continue call e287_Exception @spush_continue: mov si, positive_sign (* extremes always positive*) sub dx, dx sub cx, cx (* zero or infinite, fraction = 0*) jmp @spush_result (*INCLUDE EMUINIT.ASM*) (*public*) e287_Exception : (* called when exception detected. Exception status in CH. Process *) (* CH against mask. If masked, then return to caller to continue with *) (* default exception processing, otherwise do exception exit. If *) (* allowed, then stop the current instruction and forget all partial *) (* calculations, just return. *) (* Precision Errors do not stop completion of the current operation, *) (* being considered a "mild" error, and often being signalled with the *) (* default action of a harsher masked error, though this if unmasked *) (* does cause error interrupt when the next instruction is attempted. *) (* The routine is designed also to be called with CH=0 after the Control *) (* Word (and hence the mask) is changed, and can thus clear Error Summary *) (* as well as set it. *) (* The error is actually detected when the program tries to execute the *) (* next numerics instruction. Presumably the iNDP-87 is designed like *) (* that because the iAPX-?86 may have switched into unrelated code by *) (* the time the exception occurs, and the iNDP-87 can't be sure the *) (* CPU will be in the correct context until the next numeric instruction. *) push ax (* must not upset caller's registers *) push cx push ds push ss pop ds (* some callers may use other ds values *) mov al, [emws_status] mov cl, [emws_control] and cl, 3FH (* select the exception mask bits *) xor cl, 3FH (* invert the mask for convenience *) or al, ch (* add in new exceptions *) mov ah, al and ah, cl xor ah, al (* AH shows masked exceptions. *) test ah, fe_Overflow (* if Overflow but masked, *) jz notMaskedOver or al, 20H (* then convert into Precision error *) notMaskedOver: test al, cl (* any unmasked exceptions ? *) jz resumeDefault or al, 80H (* set Error Summary *) mov [emws_status], al and al, cl cmp al, 20H (* if only a precision error, resume. *) (*ifndef instantException jne ExceptionExit (* else abort current instruction. *) else*) jne QuitException (* emulate exception interrupt now. *) (*endif*) pop ds pop cx pop ax ret 0 (* resume if precision error. *) resumeDefault: and al, 7FH (* clear Error Summary *) mov [emws_status], al pop ds pop cx pop ax ret 0 public QuitException : mov sp, ss:[emws_resume_bp] pop bp pop es pop ds pop di pop si pop dx pop cx pop bx pop ax extrn __FloatEmulNMI jmp far __FloatEmulNMI end