From: S_JUFFA@iravcl.ira.uka.de (|S| Norbert Juffa)
Newsgroups: comp.sys.ibm.pc.hardware,comp.sys.intel
Subject: What you always wanted to know about math copros 5/5
Date: 28 Sep 1992 14:11:28 GMT
Organization: University of Karlsruhe (FRG) - Informatik Rechnerabt.
Distribution: world
NNTP-Posting-Host: irav1.ira.uka.de
X-News-Reader: VMS NEWS 1.23

         Appendix A


         {$N+,E+}
         PROGRAM PCtrl;

         VAR B,c: EXTENDED;
             Precision, L: WORD;

         PROCEDURE SetPrecisionControl (Precision: WORD);
         (* This procedure sets the internal precision of the NDP. Available *)
         (* precision values:  0  -  24 bits (SINGLE)                        *)
         (*                    1  -  n.a. (mapped to single)                 *)
         (*                    2  -  53 bits (DOUBLE)                        *)
         (*                    3  -  64 bits (EXTENDED)                      *)

         VAR CtrlWord: WORD;

         BEGIN {SetPrecisionCtrl}
            IF Precision = 1 THEN
               Precision := 0;
            Precision := Precision SHL 8; { make mask for PC field in ctrl word}
            ASM
               FSTCW    [CtrlWord]        { store NDP control word }
               MOV      AX, [CtrlWord]    { load control word into CPU }
               AND      AX, 0FCFFh        { mask out precision control field }
               OR       AX, [Precision]   { set desired precision in PC field }
               MOV      [CtrlWord], AX    { store new control word }
               FLDCW    [CtrlWord]        { set new precision control in NDP }
            END;
         END; {SetPrecisionCtrl}

         BEGIN {main}
            FOR Precision := 1 TO 3 DO BEGIN
               B := 1.2345678901234567890;
               SetPrecisionControl (Precision);
               FOR L := 1 TO 20 DO BEGIN
                  B := Sqrt (B);
               END;
               FOR L := 1 TO 20 DO BEGIN
                  B := B*B;
               END;
               SetPrecisionControl (3);   { full precision for printout }
               WriteLn (Precision, B:28);
            END;
         END.


         +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

         {$N+,E+}
         PROGRAM RCtrl;

         VAR B,c: EXTENDED;
             RoundingMode, L: WORD;


         PROCEDURE SetRoundingMode (RCMode: WORD);
         (* This procedure selects one of four available rounding modes *)
         (* 0  -  Round to nearest (default)                            *)
         (* 1  -  Round down (towards negative infinity)                *)
         (* 2  -  Round up (towards positive infinity)                  *)
         (* 3  -  Chop (truncate, round towards zero)                   *)

         VAR CtrlWord: WORD;

         BEGIN
            RCMode := RCMode SHL 10;      { make mask for RC field in control word}
            ASM
               FSTCW    [CtrlWord]        { store NDP control word }
               MOV      AX, [CtrlWord]    { load control word into CPU }
               AND      AX, 0F3FFh        { mask out rounding control field }
               OR       AX, [RCMode]      { set desired precision in RC field }
               MOV      [CtrlWord], AX    { store new control word }
               FLDCW    [CtrlWord]        { set new rounding control in NDP }
            END;
         END;

         BEGIN
            FOR RoundingMode := 0 TO 3 DO BEGIN
               B := 1.2345678901234567890e100;
               SetRoundingMode (RoundingMode);
               FOR L := 1 TO 51 DO BEGIN
                  B := Sqrt (B);
               END;
                  FOR L := 1 TO 51 DO BEGIN
                  B := -B*B;
               END;
               SetRoundingMode (0);        { round to nearest for printout }
               WriteLn (RoundingMode, B:28);
            END;
         END.


         +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

         {$N+,E+}

         PROGRAM DenormTs;

         VAR E: EXTENDED;
             D: DOUBLE;
             S: SINGLE;

         BEGIN
            WriteLn ('Testing support and printing of denormals');
            WriteLn;
            Write ('Coprocessor is: ');
            CASE Test8087 OF
               0: WriteLn ('Emulator');
               1: WriteLn ('8087 or compatible');
               2: WriteLn ('80287 or compatible');
               3: WriteLn ('80387 or compatible');
            END;
            WriteLn;
            S := 1.18e-38;
            S := S * 3.90625e-3;
            IF S = 0 THEN
               WriteLn ('SINGLE denormals not supported')
            ELSE BEGIN
               WriteLn ('SINGLE denormals supported');
               WriteLn ('SINGLE denormal prints as:   ', S);
               WriteLn ('Denormal should be printed as 4.60943...E-0041');
            END;
            WriteLn;
            D := 2.24e-308;
            D := D * 3.90625e-3;
            IF D = 0 THEN
               WriteLn ('DOUBLE denormals not supported')
            ELSE BEGIN
               WriteLn ('DOUBLE denormals supported');
               WriteLn ('DOUBLE denormal prints as:   ', D);
               WriteLn ('Denormal should be printed as 8.75...E-0311');
            END;
            WriteLn;
            E := 3.37e-4932;
            E := E * 3.90625e-3;
            IF E = 0 THEN
               WriteLn ('EXTENDED denormals not supported')
            ELSE BEGIN
               WriteLn ('EXTENDED denormals supported');
               WriteLn ('EXTENDED denormal prints as: ', E);
               WriteLn ('Denormal should be printed as 1.3164...E-4934');
            END;
         END.





         Appendix B


         ; FILE: APFELM4.ASM
         ; assemble with MASM /e APFELM4 or TASM /e APFELM4


         CODE        SEGMENT BYTE PUBLIC 'CODE'
                     ASSUME  CS: CODE

                     PAGE    ,120

                     PUBLIC  APPLE87;

         APPLE87     PROC    NEAR
                     PUSH    BP                  ; save caller's base pointer
                     MOV     BP, SP              ; make new frame pointer
                     PUSH    DS                  ; save caller's data segment
                     PUSH    SI                  ; save register
                     PUSH    DI                  ;  variables
                     LDS     BX, [BP+04]         ; pointer to parameter record
                     FINIT                       ; init 80x87          FSP->R0
                     FILD   WORD  PTR [BX+02]    ; maxrad              FSP->R7
                     FLD    QWORD PTR [BX+08]    ; qmax                FSP->R6
                     FSUB   QWORD PTR [BX+16]    ; qmax-qmin           FSP->R6
                     DEC    WORD  PTR [BX+04]    ; ymax-1
                     FIDIV  WORD  PTR [BX+04]    ; (qmax-qmin)/(ymax-1)FSP->R6
                     FSTP   QWORD PTR [BX+16]    ; save delta_q        FSP->R7
                     FLD    QWORD PTR [BX+24]    ; pmax                FSP->R6
                     FSUB   QWORD PTR [BX+32]    ; pmax-pmin           FSP->R6
                     DEC    WORD  PTR [BX+06]    ; xmax-1
                     FIDIV  WORD  PTR [BX+06]    ; delta_p             FSP->R6
                     MOV    AX, [BX]             ; save maxiter,[BX] needed for
                     MOV    [BX+2], AX           ;  80x87 status now
                     XOR    BP, BP               ; y=0
                     FLD    QWORD PTR [BX+08]    ; qmax                FSP->R5
                     CMP    WORD  PTR [BX+40], 0 ; fast mode on 8087 desired ?
                     JE     yloop                ; no, normal mode
                     FSTCW  [BX]                 ; save NDP control word
                     AND    WORD PTR [BX], 0FCFFh; set PCTRL = single precision
                     FLDCW  [BX]                 ; get back NDP control word
         yloop:      XOR    DI, DI               ; x=0
                     FLD    QWORD PTR [BX+32]    ; pmin                FSP->R4
         xloop:      FLDZ                        ; j**2= 0             FSP->R3
                     FLDZ                        ; 2ij = 0             FSP->R2
                     FLDZ                        ; i**2= 0             FSP->R1
                     MOV    CX, [BX+2]           ; maxiter
                     MOV    DL, 41h              ; mask for C0 and C3 cond.bits
         iteration:  FSUB   ST, ST(2)            ; i**2-j**2           FSP->R1
                     FADD   ST, ST(3)            ; i**2-j**2+p = i     FSP->R1
                     FLD    ST(0)                ; duplicate i         FSP->R0
                     FMUL   ST(1), ST            ; i**2                FSP->R0
                     FADD   ST, ST(0)            ; 2i                  FSP->R0
                     FXCH   ST(2)                ; 2*i*j               FSP->R0
                     FADD   ST, ST(5)            ; 2*i*j+q = j         FSP->R0
                     FMUL   ST(2), ST            ; 2*i*j               FSP->R0
                     FMUL   ST, ST(0)            ; j**2                FSP->R0
                     FST    ST(3)                ; save j**2           FSP->R0
                     FADD   ST, ST(1)            ; i**2+j**2           FSP->R0
                     FCOMP  ST(7)                ; i**2+j**2 > maxrad? FSP->R1
                     FSTSW  [BX]                 ; save 80x87 cond.codeFSP->R1
                     TEST   BYTE PTR [BX+1], DL  ; test carry and zero flags
                     LOOPNZ iteration            ; until maxiter if not diverg.
                     MOV    DX, CX               ; number of loops executed
                     NEG    CX                   ; carry set if CX <> 0
                     ADC    DX, 0                ; adjust DX if no. of loops<>0

                     ; plot point here (DI = X, BP = y, DX has the color)

                     FSTP   ST(0)                ; pop i**2            FSP->R2
                     FSTP   ST(0)                ; pop 2ij             FSP->R3
                     FSTP   ST(0)                ; pop j**2            FSP->R4
                     FADD   ST,ST(2)             ; p=p+delta_p         FSP->R4
                     INC    DI                   ; x:=x+1
                     CMP    DI, [BX+6]           ; x > xmax ?
                     JBE    xloop                ; no, continue on same line
                     FSTP   ST(0)                ; pop p               FSP->R5
                     FSUB   QWORD PTR [BX+16]    ; q=q-delta_q         FSP->R5
                     INC    BP                   ; y:=y+1
                     CMP    BP, [BX+4]           ; y > ymax ?
                     JBE    yloop                ; no, picture not done yet

         groesser:   POP    DI                   ; restore
                     POP    SI                   ;  register variables
                     POP    DS                   ; restore caller's data segm.
                     POP    BP                   ; save caller's base pointer
                     RET    4                    ; pop parameters and return
         APPLE87     ENDP

         CODE        ENDS

                     END

         ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

         UNIT Time;

         INTERFACE

         FUNCTION Clock: LONGINT;          { same as VMS; time in milliseconds }


         IMPLEMENTATION

         FUNCTION Clock: LONGINT; ASSEMBLER;
         ASM
                    PUSH    DS            { save caller's data segment }
                    XOR     DX, DX        { initialize data segment to }
                    MOV     DS, DX        {  access ticker counter }
                    MOV     BX, 46Ch      { offset of ticker counter in segm.}
                    MOV     DX, 43h       { timer chip control port }
                    MOV     AL, 4         { freeze timer 0 }
                    PUSHF                 { save caller's int flag setting }
                    STI                   { allow update of ticker counter }
                    LES     DI, DS:[BX]   { read BIOS ticker counter }
                    OUT     DX, AL        { latch timer 0 }
                    LDS     SI, DS:[BX]   { read BIOS ticker counter }
                    IN      AL, 40h       { read latched timer 0 lo-byte }
                    MOV     AH, AL        { save lo-byte }
                    IN      AL, 40h       { read latched timer 0 hi-byte }
                    POPF                  { restore caller's int flag }
                    XCHG    AL, AH        { correct order of hi and lo }
                    MOV     CX, ES        { ticker counter 1 in CX:DI:AX }
                    CMP     DI, SI        { ticker counter updated ? }
                    JE      @no_update    { no }
                    OR      AX, AX        { update before timer freeze ? }
                    JNS     @no_update    { no }
                    MOV     DI, SI        { use second }
                    MOV     CX, DS        {  ticker counter }
         @no_update:NOT     AX            { counter counts down }
                    MOV     BX, 36EDh     { load multiplier }
                    MUL     BX            { W1 * M }
                    MOV     SI, DX        { save W1 * M (hi) }
                    MOV     AX, BX        { get M }
                    MUL     DI            { W2 * M }
                    XCHG    BX, AX        { AX = M, BX = W2 * M (lo) }
                    MOV     DI, DX        { DI = W2 * M (hi) }
                    ADD     BX, SI        { accumulate }
                    ADC     DI, 0         {  result }
                    XOR     SI, SI        { load zero }
                    MUL     CX            { W3 * M }
                    ADD     AX, DI        { accumulate }
                    ADC     DX, SI        {  result in DX:AX:BX }
                    MOV     DH, DL        { move result }
                    MOV     DL, AH        {  from DL:AX:BX }
                    MOV     AH, AL        {   to }
                    MOV     AL, BH        {    DX:AX:BH }
                    MOV     DI, DX        { save result }
                    MOV     CX, AX        {  in DI:CX }
                    MOV     AX, 25110     { calculate correction }
                    MUL     DX            {  factor }
                    SUB     CX, DX        { subtract correction }
                    SBB     DI, SI        {  factor }
                    XCHG    AX, CX        { result back }
                    MOV     DX, DI        {  to DX:AX }
                    POP     DS            { restore caller's data segment }
         END;


         BEGIN
            Port [$43] := $34;           { need rate generator, not square wave}
            Port [$40] := 0;             { generator as prog. by some BIOSes }
            Port [$40] := 0;             { for timer 0 }
         END. { Time }


         ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

         {$A+,B-,R-,I-,V-,N+,E+}
         PROGRAM PeakFlop;

         USES Time;

         TYPE ParamRec = RECORD
                            MaxIter, MaxRad, YMax, XMax: WORD;
                            Qmax, Qmin, Pmax, Pmin: DOUBLE;
                            FastMod: WORD;
                            PlotFkt: POINTER;
                            FLOPS:LONGINT;
                         END;

         VAR Param: ParamRec;
             Start: LONGINT;


         {$L APFELM4.OBJ}

         PROCEDURE Apple87 (VAR Param: ParamRec);     EXTERNAL;


         BEGIN
            WITH Param DO BEGIN
               MaxIter:= 50;
               MaxRad := 30;
               YMax   := 30;
               XMax   := 30;
               Pmin   :=-2.1;
               Pmax   := 1.1;
               Qmin   :=-1.2;
               Qmax   := 1.2;
               FastMod:= Word (FALSE);
               PlotFkt:= NIL;
               Flops  := 0;
            END;
            Start := Clock;
            Apple87 (Param);         { executes 104002 FLOP }
            Start := Clock - Start;  { elapsed time in milliseconds }
            WriteLn ('Peak-MFLOPS: ', 104.002 / Start);
         END.

         ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

         ; FILE: M4X4.ASM
         ;
         ; assemble with TASM /e M4X4 or MASM /e M4X4

         CODE      SEGMENT BYTE PUBLIC 'CODE'

                   ASSUME  CS:CODE

                   PUBLIC  MUL_4x4
                   PUBLIC  IIT_MUL_4x4


         FSBP0     EQU     DB  0DBh, 0E8h        ; declare special IIT
         FSBP1     EQU     DB  0DBh, 0EBh        ;  instructions
         FSBP2     EQU     DB  0DBh, 0EAh
         F4X4      EQU     DB  0DBh, 0F1h


         ;---------------------------------------------------------------------
         ;
         ; MUL_4x4 multiplicates a four-by-four matrix by an array of four
         ; dimensional vectors. This operation is needed for 3D transformations
         ; in graphics data processing. There are arrays for each component of
         ; a vector. Thus there is an ; array containing all the x components,
         ; another containing all the y components and so on. Each component is
         ; an 8 byte IEEE floating point number. Two indices into the array of
         ; vectors are given. The first is the index of the vector that will be
         ; processed first, the second is the index of the vector processed
         ; last.
         ;
         ;---------------------------------------------------------------------

         MUL_4x4   PROC    NEAR

                   AddrX   EQU DWORD PTR [BP+24] ; address of X component array
                   AddrY   EQU DWORD PTR [BP+20] ; address of Y component array
                   AddrZ   EQU DWORD PTR [BP+16] ; address of Z component array
                   AddrW   EQU DWORD PTR [BP+12] ; address of W component array
                   AddrT   EQU DWORD PTR [BP+8]  ; addr. of 4x4 transform. mat.
                   F       EQU WORD  PTR [BP+6]  ; first vector to process
                   K       EQU WORD  PTR [BP+4]  ; last vector to process
                   RetAddr EQU WORD  PTR [BP+2]  ; return address saved by call
                   SavdBP  EQU WORD  PTR [BP+0]  ; saved frame pointer
                   SavdDS  EQU WORD  PTR [BP-2]  ; caller's data segment

                   PUSH    BP                    ; save TURBO-Pascal frame pointer
                   MOV     BP, SP                ; new frame pointer
                   PUSH    DS                    ; save TURBO-Pascal data segment

                   MOV     CX, K                 ; final index
                   SUB     CX, F                 ; final index - start index
                   JNC     $ok                   ; must not
                   JMP     $nothing              ;  be negative
         $ok:      INC     CX                    ; number of elements

                   MOV     SI, F                 ; init offset into arrays
                   SHL     SI, 1                 ; each
                   SHL     SI, 1                 ;  element
                   SHL     SI, 1                 ;   has 8 bytes

                   LDS     DI, AddrT             ; addr. of transformation mat.
                   FLD     QWORD PTR [DI]        ; load a[0,0]   = R7
                   FLD     QWORD PTR [DI+8]      ; load a[0,1]   = R6

         $mat_mul: LES     BX, AddrX             ; addr. of x component array
                   FLD     QWORD PTR ES:[BX+SI]  ; load x[a]     = R5
                   LES     BX, AddrY             ; addr. of y component array
                   FLD     QWORD PTR ES:[BX+SI]  ; load y[a]     = R4
                   LES     BX, AddrZ             ; addr. of z component array
                   FLD     QWORD PTR ES:[BX+SI]  ; load z[a]     = R3
                   LES     BX, AddrW             ; addr. of w component array
                   FLD     QWORD PTR ES:[BX+SI]  ; load w[a]     = R2

                   FLD     ST(5)                 ; load a[0,0]   = R1
                   FMUL    ST, ST(4)             ; a[0,0] * x[a] = R1
                   FLD     ST(5)                 ; load a[0,1]   = R0
                   FMUL    ST, ST(4)             ; a[0,1] * y[a] = R0
                   FADDP   ST(1), ST             ; a[0,0]*x[a]+a[0,1]*y[a]=R1
                   FLD     QWORD PTR [DI+16]     ; load a[0,2]   = R0
                   FMUL    ST, ST(3)             ; a[0,2] * z[a] = R0
                   FADDP   ST(1), ST             ; a[0,0]*x[a]...a[0,2]*z[a]=R1
                   FLD     QWORD PTR [DI+24]     ; load a[0,3]   = R0
                   FMUL    ST, ST(2)             ; a[0,3] * w[a] = R0
                   FADDP   ST(1), ST             ; a[0,0]*x[a]...a[0,3]*w[a]=R1
                   LES     BX, AddrX             ; get address of x vector
                   FSTP    QWORD PTR ES:[BX+SI]  ; write new x[a]

                   FLD     QWORD PTR [DI+32]     ; load a[1,0]   = R1
                   FMUL    ST, ST(4)             ; a[1,0] * x[a] = R1
                   FLD     QWORD PTR [DI+40]     ; load a[1,1]   = R0
                   FMUL    ST, ST(4)             ; a[1,1] * y[a] = R0
                   FADDP   ST(1), ST             ; a[1,0]*x[a]+a[1,1]*y[a]=R1
                   FLD     QWORD PTR [DI+48]     ; load a[1,2]   = R0
                   FMUL    ST, ST(3)             ; a[1,2] * z[a] = R0
                   FADDP   ST(1), ST             ; a[1,0]*x[a]...a[1,2]*z[a]=R1
                   FLD     QWORD PTR [DI+56]     ; load a[1,3]   = R0
                   FMUL    ST, ST(2)             ; a[1,3] * w[a] = R0
                   FADDP   ST(1), ST             ; a[1,0]*x[a]...a[1,3]*w[a]=R1
                   LES     BX, AddrY             ; get address of y vector
                   FSTP    QWORD PTR ES:[BX+SI]  ; write new y[a]

                   FLD     QWORD PTR [DI+64]     ; load a[2,0]   = R1
                   FMUL    ST, ST(4)             ; a[2,0] * x[a] = R1
                   FLD     QWORD PTR [DI+72]     ; load a[2,1]   = R0
                   FMUL    ST, ST(4)             ; a[2,1] * y[a] = R0
                   FADDP   ST(1), ST             ; a[2,0]*x[a]+a[2,1]*y[a]=R1
                   FLD     QWORD PTR [DI+80]     ; load a[2,2]   = R0
                   FMUL    ST, ST(3)             ; a[2,2] * z[a] = R0
                   FADDP   ST(1), ST             ; a[2,0]*x[a]...a[2,2]*z[a]=R1
                   FLD     QWORD PTR [DI+88]     ; load a[2,3]   = R0
                   FMUL    ST, ST(2)             ; a[2,3] * w[a] = R0
                   FADDP   ST(1), ST             ; a[2,0]*x[a]...a[2,3]*w[a]=R1
                   LES     BX, AddrZ             ; get address of z vector
                   FSTP    QWORD PTR ES:[BX+SI]  ; write new z[a]

                   FLD     QWORD PTR [DI+96]     ; load a[3,0]   = R1
                   FMULP   ST(4), ST             ; a[3,0] * x[a] = R5
                   FLD     QWORD PTR [DI+104]    ; load a[3,1]   = R1
                   FMULP   ST(3), ST             ; a[3,1] * y[a] = R4
                   FLD     QWORD PTR [DI+112]    ; load a[3,2]   = R1
                   FMULP   ST(2), ST             ; a[3,2] * z[a] = R3
                   FLD     QWORD PTR [DI+120]    ; load a[3,3]   = R1
                   FMULP   ST(1), ST             ; a[3,3] * w[a] = R2
                   FADDP   ST(1), ST             ; a[3,3]*w[a]+a[3,2]*z[a]=R3
                   FADDP   ST(1), ST             ; a[3,3]*w[a]...a[3,1]*y[a]=R4
                   FADDP   ST(1), ST             ; a[3,3]*w[a]...a[3,0]*x[a]=R5
                   LES     BX, AddrW             ; get address of w vector
                   FSTP    QWORD PTR ES:[BX+SI]  ; write new w[a]

                   ADD     SI, 8                 ; new offset into arrays
                   DEC     CX                    ; decrement element counter
                   JZ      $done                 ; no elements left, done
                   JMP     $mat_mul              ; transform next vector

         $done:    FSTP     ST(0)                ; clear
                   FSTP     ST(0)                ;  FPU stack
         $nothing: POP      DS                   ; restore TP data segment
                   POP      BP                   ; restore TP frame pointer
                   RET      24                   ; pop parameters and return

         MUL_4X4   ENDP


         ;---------------------------------------------------------------------
         ;
         ; IIT_MUL_4x4 multiplicates a four-by-four matrix by an array of four
         ; dimensional vectors. This operation is needed for 3D transformations
         ; in graphics data processing. There are arrays for each component of
         ; a vector.  Thus there is an array containing all the x components,
         ; another containing all the y components and so on. Each component is
         ; an 8 byte IEEE floating point number. Two indices into the array of
         ; vectors are given. The first is the index of the vector that will be
         ; processed first, the second is the index of the vector processed
         ; last. This subroutine uses the special instructions only available
         ; on IIT coprocessors to provide fast matrix multiply capabilities.
         ; So make sure to use it only on IIT coprocessors.
         ;
         ;---------------------------------------------------------------------

         IIT_MUL_4x4   PROC    NEAR

                   AddrX   EQU DWORD PTR [BP+24] ; address of X component array
                   AddrY   EQU DWORD PTR [BP+20] ; address of Y component array
                   AddrZ   EQU DWORD PTR [BP+16] ; address of Z component array
                   AddrW   EQU DWORD PTR [BP+12] ; address of W component array
                   AddrT   EQU DWORD PTR [BP+8]  ; addr. of 4x4 transf. matrix
                   F       EQU WORD  PTR [BP+6]  ; first vector to process
                   K       EQU WORD  PTR [BP+4]  ; last vector to process
                   RetAddr EQU WORD  PTR [BP+2]  ; return address saved by call
                   SavdBP  EQU WORD  PTR [BP+0]  ; saved frame pointer
                   SavdDS  EQU WORD  PTR [BP-2]  ; caller's data segment
                   Ctrl87  EQU WORD  PTR [BP-4]  ; caller's 80x87 control word

                   PUSH    BP                    ; save TURBO-Pascal frame ptr
                   MOV     BP, SP                ; new frame pointer
                   PUSH    DS                    ; save TURBO-Pascal data seg.
                   SUB     SP, 2                 ; make local variabe
                   FSTCW   [Ctrl87]              ; save 80x87 ctrl word
                   LES     SI, AddrT             ; ptr to transformation matrix
                   FINIT                         ; initialize coprocessor
                   FSBP2                         ; set register bank 2
                   FLD     QWORD PTR ES:[SI]     ; load a[0,0]
                   FLD     QWORD PTR ES:[SI+32]  ; load a[1,0]
                   FLD     QWORD PTR ES:[SI+64]  ; load a[2,0]
                   FLD     QWORD PTR ES:[SI+96]  ; load a[3,0]
                   FLD     QWORD PTR ES:[SI+8]   ; load a[0,1]
                   FLD     QWORD PTR ES:[SI+40]  ; load a[1,1]
                   FLD     QWORD PTR ES:[SI+72]  ; load a[2,1]
                   FLD     QWORD PTR ES:[SI+104] ; load a[3,1]
                   FINIT                         ; initialize coprocessor
                   FSBP1                         ; set register bank 1
                   FLD     QWORD PTR ES:[SI+16]  ; load a[0,2]
                   FLD     QWORD PTR ES:[SI+48]  ; load a[1,2]
                   FLD     QWORD PTR ES:[SI+80]  ; load a[2,2]
                   FLD     QWORD PTR ES:[SI+112] ; load a[3,2]
                   FLD     QWORD PTR ES:[SI+24]  ; load a[0,3]
                   FLD     QWORD PTR ES:[SI+56]  ; load a[1,3]
                   FLD     QWORD PTR ES:[SI+88]  ; load a[2,3]
                   FLD     QWORD PTR ES:[SI+120] ; load a[3,3]

                                                 ; transformation matrix loaded

                   MOV     AX, F                 ; index of first vector
                   MOV     DX, K                 ; index of last vector

                   MOV     BX, AX                ; index 1st vector to process
                   MOV     CL, 3                 ; component has 8 (2**3) bytes
                   SHL     BX, CL                ; compute offset into arrays

                   FINIT                         ; initialize coprocessor
                   FSBP0                         ; set register bank 0

         $mat_loop:LES     SI, AddrW             ; addr. of W component array
                   FLD     QWORD PTR ES:[SI+BX]  ; W component current vector
                   LES     SI, AddrZ             ; addr. of Z component array
                   FLD     QWORD PTR ES:[SI+BX]  ; Z component current vector
                   LES     SI, AddrY             ; addr. of Y component array
                   FLD     QWORD PTR ES:[SI+BX]  ; Y component current vector
                   LES     SI, AddrX             ; addr. of X component array
                   FLD     QWORD PTR ES:[SI+BX]  ; X component current vector
                   F4X4                          ; mul 4x4 matrix by 4x1 vector
                   INC     AX                    ; next vector
                   MOV     DI, AX                ; next vector
                   SHL     DI, CL                ; offset of vector into arrays

                   FSTP    QWORD PTR ES:[SI+BX]  ; store X comp. of curr. vect.
                   LES     SI, AddrY             ; address of Y component array
                   FSTP    QWORD PTR ES:[SI+BX]  ; store Y comp. of curr. vect.
                   LES     SI, AddrZ             ; address of Z component array
                   FSTP    QWORD PTR ES:[SI+BX]  ; store Z comp. of curr. vect.
                   LES     SI, AddrW             ; address of W component array
                   FSTP    QWORD PTR ES:[SI+BX]  ; store W comp. of curr. vect.

                   MOV     BX, DI                ; ofs nxt vect. in comp. arrays
                   CMP     AX, DX                ; nxt vector past upper bound?
                   JLE     $mat_loop             ; no, transform next vector
                   FLDCW   [Ctrl87]              ; restore orig 80x87 ctrl word

                   ADD      SP, 2                ; get rid of local variable
                   POP      DS                   ; restore TP data segment
                   POP      BP                   ; restore TP frame pointer
                   RET      24                   ; pop parameters and return
         IIT_MUL_4x4   ENDP

         CODE      ENDS

                   END

         ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

         {$N+,E+}

         PROGRAM Trnsform;

         USES Time;

         CONST VectorLen = 8190;

         TYPE  Vector    = ARRAY [0..VectorLen] OF DOUBLE;
               VectorPtr = ^Vector;
               Mat4      = ARRAY [1..4, 1..4] OF DOUBLE;

         VAR   X: VectorPtr;
               Y: VectorPtr;
               Z: VectorPtr;
               W: VectorPtr;
               T: Mat4;
               K: INTEGER;
               L: INTEGER;
               First: INTEGER;
               Last:  INTEGER;
               Start: LONGINT;
               Elapsed:LONGINT;

         PROCEDURE MUL_4X4     (X, Y, Z, W: VectorPtr;
                                VAR T: Mat4; First, Last: INTEGER); EXTERNAL;
         PROCEDURE IIT_MUL_4X4 (X, Y, Z, W: VectorPtr;
                                VAR T: Mat4; First, Last: INTEGER); EXTERNAL;

         {$L M4X4.OBJ}

         BEGIN
            WriteLn ('Test8087 = ', Test8087);
            New (X);
            New (Y);
            New (Z);
            New (W);
            FOR L := 1 TO VectorLen DO BEGIN
               X^ [L] := Random;
               Y^ [L] := Random;
               Z^ [L] := Random;
               W^ [L] := Random;
            END;
            X^ [0] := 1;
            Y^ [0] := 1;
            Z^ [0] := 1;
            W^ [0] := 1;
            FOR K := 1 TO 4 DO BEGIN
               FOR L := 1 TO 4 DO BEGIN
                  T [K, L] := (K-1)*4 + L;
               END;
            END;
            First := 0;
            Last  := 8190;
            Start := Clock;
            MUL_4X4 (X, Y, Z, W, T, First, Last);
            { IIT_MUL_4X4 (X, Y, Z, W, T, First, Last); }
            Elapsed := Clock - Start;
            WriteLn ('Number of vectors: ', Last-First+1);
            WriteLn ('Time: ', Elapsed, ' ms');
            WriteLn ('Equivalent to ', (28.0*(Last-First+1)/1e6)/
                     (Elapsed*1e-3):0:4, ' MFLOPS');
            WriteLn;
            WriteLn ('Last vector:');
            WriteLn;
            WriteLn (X^[Last]);
            WriteLn (Y^[Last]);
            WriteLn (Z^[Last]);
            WriteLn (W^[Last]);
         END.
Newsgroups: comp.sys.ibm.pc.hardware