TITLE 'P12DSH1.MLC - Calc BFP 128 Bit Standard Deviation'
*********************************************************************
* Program ID.  P12DSH1.MLC
* Author.      Don Higgins.
* Date.        06/20/07.
*********************************************************************
* 01/14/08 DSH1 convert previous test program for solution to #12
*********************************************************************
P10DSH1  ZMFACC CODE,START,NAME='DON HIGGINS'
         COPY  ASMMSP      STRUCTURED PROGRAMMING MACRO LIBRARY
         CALL  STDDEVLB,(LBX,N,LBSD)  CALCULATE STANDARD DEVIATION
*
* DISPLAY STD DEVIATION IN SCIENTIFIC NOTATION FOR DEBUGGIN
*
         CTD   CTD_LB,IN=LBSD,OUT=DLBSD
         ZMFACC CODE,END
N        DC    F'4'        INTEGER NUMBER OF ARRAY ELEMENTS
         DS     0D
         ZMFACC INPUT,START
LBX      DC    LB'1,2,3,6' BINARY EXTENDED PRECISION (BFP) ARRAY 
         ZMFACC INPUT,END
         ZMFACC OUTPUT,START
LBSD     DC    LB'0'       STANDARD DEVIATION BFP VALUE
DLBSD    DC    CL45' '     STANDARD DEVIATION IN SCIENTIFIC NOTATION 
         ZMFACC OUTPUT,END
         TITLE 'STDDEVLB - CALC BFP 128 BIT STANDARD DEVIATION'
STDDEVLB CSECT
         USING *,15
         SAVE  (14,12)
         BALR  12,0
         USING *,12
*
* CALCULATE BFP 128 BIT STANDARD DEVIATION FOR ARRAY OF BFP ELEMENTS
*
* CALL PARAMETERS:
*
*  1 = INPUT  ADDRESS OF BFP ARRAY ELEMENTS (16 BYTE BFP LB FORMAT)
*  2 - INPUT  FULL WORD ARRAY ELEMENT COUNT (MUST BE >= 1)
*  3 - OUTPUT STANDARD DEVIATION IN 16 BYTE BFP LB FORMAT
*
* RETURN CODES:
*
*   0 - STANDARD DEVIATION STORED SUCCESSFULLY IN PARAMETER 3
*  16 - ELEMENT COUNT LESS THAN 1 (STD.DEV. UNDEFINED)
*
*
* CALCULATE BFP LBXM MEAN = (SUM X(I))/N AND SAVE IN F4+F6
*
         LM    R3,R5,0(R1) R3=A(ARRAY) R4=A(N) R5=A(LBSD)
         L     R4,0(R4)    R4=N ELEMENT COUNT
         LR    R6,R3       R6=AARRY ELEMENT ADDRESS
         LR    R7,R4       R7=ELEMENTS TO PROCESS IN LOOP         
         IF    (CHI,R7,LT,1)
               RETURN (14,12),RC=16  EXIT WITH RC=16 IF COUNT < 1
         ENDIF
         STFPC  R8      SAVE FPC IN R8
         LR     R1,R8
FPC_MASK_UNF EQU X'01000000' BFP UNDERFLOW MASK BIT
         N      R1,=A(X'FFFFFFFF'-FPC_MASK_UNF)  TURN OFF UNDERFLOW
         LFPC   R1
         SXBR  F4,F4       F4+F6 = SUM X(I) TO CALC MEAN
         DO    WHILE=(CHI,R7,GT,0)
               LD  F1,0(R6) F1+F3 = ELEMENT
               LD  F3,8(R6)
               AXBR F4,F1  ADD ELEMENT TO SUM
               AHI R6,16   NEXT ELEMENT
               AHI R7,-1   REDUCE ELEMENTS REMAINING
         ENDDO
         CXFBR F1,R4       F1+F3 = N
         DXBR  F4,F1       F4+F6 = F4+F6 / F1+F3 = BFP LB MEAN
*
* CALCULATE SUM OF VARIANCE SQUARED IN F0+F2
*       
         LR    R6,R3       R6 = ADDRESS FIRST ARRAY ELEMENT
         LR    R7,R4       R7 = ELEMENTS TO PROCESS IN LOOP
         SXR   F0,F0       F0+F2 = (SUM X(I)-MEAN)**2
         DO    WHILE=(CHI,R7,GT,0)
               LD  F1,0(R6) F1+F3 = BFP LB ELEMENT
               LD  F3,8(R6)
               SXBR F1,F4  SUBTRACT MEAN
               MXBR F1,F1  SQUARE   DIFF
               AXBR F0,F1  ADD      DIFF * DIFF
               AHI R6,16   NEXT ELEMENT
               AHI R7,-1   REDUCE ELEMENTS REMAINING
         ENDDO
         CXFBR F1,R4       F1+F3 = N
         DXBR  F0,F1       F0+F2 = F0+F2 / F1+F3 = (SUM DIFF*DIFF)/N
*
* CALCULATE STANDARD DEVIATION
*
         SQXBR F0,F0     CALC SQRT SUM OF VAR**2 IN F0+F2
         STD   F0,0(R5)  STORE F0+F2 LB STANDARD DEVIATION RESULT
         STD   F2,8(R5)
         LFPC  R8        RESTORE FPC REG
         RETURN (14,12),RC=0
LBXM     DS    LB        VARIANCE FOR DEBUG DISPLAY
LBXV2    DS    LB        (SUM DIFF*DIFF)/N FOR DEBUG DISPLAY
         END