//-------------pipercp.s-----------------------------------
//
// C Michael version 1.0 9-2-91
//           version 1.1 25-2-91 distributed loads, all loop in DIM.
//           version 1.2 pfld care reads off end of vector (10 cycles av)
//           version 1.3 12-3-91 fix so does not read off end of vector
//  tested: care result=source gave errors at DL. (hardware)
//------------------------------------------------------------------------
// Assembler routine to perform a reciprocal
//------------------------------------------------------------------------
// The source vector is single (32 bit) precision
// and the result is in single (32 bit) precision (can overwrite source).
// The result used two Newton Raphson iterations: accuracy 2^-23 = 10^-7
//    g1 = rcp(a1)  p1=2.0-a1*g1   g2=0.0+g1*p1
//                  p2=2.0-a1*g2   g3=0.0+g1*p2==> result
//    g1 to 2^-7,  g2 to 2^-14,  g3 to 2^-28 in theory
// but cancellation 2.0-ga gives 2^-23 in 0.999 error in worst case.
// This can give error in last 2 bits as Intel PRM says.
// The routine works for vectors of any length.  For a long vector
// it takes 34 cycles per 6 reciprocals in theory.
// version 1.1 n=100 takes 35 microsec (C-code 110 microsec)
// version 1.2 n=5000 * 10000 13 sec ==> 10.4 cycles (C-code 60 sec)
// with no loads or stores takes same time!
// ***** NEEDS result and vector on 64 bit boundary *****
//-------------------------------------------------------------------------
//
// ProtoType pipercp(float *result, float *vector, int n )
//
//      if( n <= 0 ) return;
//      for(i=0; i<n ; i++){
//        result[i] = 1.0 / vector[i];
//      }
//
//      n   is the number of real elements in the vector.
// n: 1 - 11  12-17   16-23  24-29  ,etc
//   scalar  + tail + loop  + loop  ,etc


        resptr   = r16
        ldptr    = r17
        n        = r18
        step     = r27
        count    = r28
        twoptr   = r31


        two       = f20
        one       = f21
        a1        = f22
        a2        = f23
        a3        = f24
        a4        = f25
        a5        = f26
        a6        = f27
        g1        = f28
        g2        = f29
        g3        = f30
        g4        = f31
        g5        = f12
        g6        = f13
        p1        = f14
        p2        = f15
        p3        = f16
        p4        = f17
        p5        = f18
        p6        = f19

.page
.align 4
_pipercp_::                         // fortran entry point
        ld.l        0(n),  n        // convert pointer to value.
_pipercp::                          // C entry

        adds  -64,sp,sp                  // space to save registers

        fst.d f12, 0(sp)
        fst.d f14, 8(sp)
        orh      ha%ctwo,r0,twoptr       // do as double
        fld.d    l%ctwo(twoptr),two      // 2.0,1.0 from low(high)
        adds        -6, r0,  step        // set up for bla
        adds        -18,  n, count       // bla decrements count by 6
        adds        -4,resptr,resptr     // for auto inc store
        adds        -4,ldptr,ldptr       // for auto inc load
        adds        -12,  n, r0
         bc          scalar              // n < 12 care ldptr,resptr
        adds        -4,resptr,resptr     // for auto inc store
        adds        -4,ldptr,ldptr       // for auto inc load
// Get pipe started -------------------------------------------
          pfld.d 8(ldptr)++,f0
          pfld.d 8(ldptr)++,f0
          pfld.d 8(ldptr)++,f0
          pfld.d 8(ldptr)++, a1
          pfld.d 8(ldptr)++, a3
          pfld.d 8(ldptr)++, a5
          frcp.ss         a1,g1
          frcp.ss         a2,g2
          frcp.ss         a3,g3
          frcp.ss         a4,g4
          frcp.ss         a5,g5
          frcp.ss         a6,g6
          m12tpm.ss       one,two, f0      // to set T = 2
          m12tpm.ss       a1,g1  , f0
          m12tpm.ss       a2,g2  , f0
.align 8
        d.m12ttpa.ss      a3,g3  , f0      // now T=2
        d.m12tsm.ss       a4,g4  , f0      // 2-a1*g1
        d.m12tsm.ss       a5,g5  , f0      // now in DIM
            bla step,count,dummy           // initialise LCC either branch
        d.m12tsm.ss       a6,g6  , f0
           nop
dummy:  d.m12tsm.ss       f0,f0  , p1      // to set T=0 ; p1 is 2-a1*g1
           adds -18,n,r0
        d.m12tsm.ss       g1,p1  , p2
            bc tail                        // skip loop if n < 12
//----------------------- main loop -------------------
loop:
        d.m12tsm.ss       g2,p2  , p3
           nop
        d.m12ttpa.ss      g3,p3  , p4      // now T=0
           nop
        d.m12tpm.ss       g4,p4  , p5      // 0+a1*p1
           nop
        d.m12tpm.ss       g5,p5  , p6
           nop
        d.m12tpm.ss       g6,p6  , f0      // output 0*0 => junk
           nop
        d.m12tpm.ss       one,two, g1      // to set T = 2; g1 is g1*p1
           nop
        d.m12tpm.ss       a1,g1  , g2
           nop
        d.m12tpm.ss       a2,g2  , g3
           nop
        d.m12ttpa.ss      a3,g3  , g4      // now T=2
           pfld.d 8(ldptr)++, a1            // dest not in next 2 insts
        d.m12tsm.ss       a4,g4  , g5      // 2-a1*g1
           nop
        d.m12tsm.ss       a5,g5  , g6
           nop
        d.m12tsm.ss       a6,g6  , f0      // output 1*2 ==> junk
           nop
        d.m12tsm.ss       f0,f0  , p1      // to set T=0 ; p1 is 2-a1*g1
           nop
        d.m12tsm.ss       g1,p1  , p2
           pfld.d 8(ldptr)++, a3
        d.m12tsm.ss       g2,p2  , p3
           nop
        d.m12ttpa.ss      g3,p3  , p4      // now T=0
           nop
        d.m12tpm.ss       g4,p4  , p5      // 0+g1*p1
           nop
        d.m12tpm.ss       g5,p5  , p6
           nop
        d.m12tpm.ss       g6,p6  , f0      // output 0*0 => junk
           nop
        d.m12tpm.ss       f0,f0  , p1      // p1 is g1*p1 = ANS
           pfld.d 8(ldptr)++, a5
        d.m12tpm.ss       f0,f0  , p2
           nop
        d.m12tpm.ss       f0,f0  , p3
           nop
        d.frcp.ss         a1,g1            // not pipelined
           nop
        d.frcp.ss         a2,g2            // frcp messes up multiplier pipe
           fst.d p1,8(resptr)++
        d.frcp.ss         a3,g3
           nop
        d.frcp.ss         a4,g4
           nop
        d.frcp.ss         a5,g5
           nop
        d.frcp.ss         a6,g6            // Adder pipe from above
           nop
        d.m12tpm.ss       one,two, p4      // to set T = 2
           nop
        d.m12tpm.ss       a1,g1  , p5
           fst.d p3,8(resptr)++
        d.m12tpm.ss       a2,g2  , p6
           nop
        d.m12ttpa.ss      a3,g3  , f0      // now T=2 ; junk from pipe
           nop
        d.m12tsm.ss       a4,g4  , f0      // 2-a1*g1
           nop
        d.m12tsm.ss       a5,g5  , f0
           nop
        d.m12tsm.ss       a6,g6  , f0      // output 1*2 ==> junk
           fst.d p5,8(resptr)++
        d.m12tsm.ss       f0,f0  , p1      // to set T=0 ; p1 is 2-a1*g1
            bla step,count,loop
        d.m12tsm.ss       g1,p1  , p2
           nop
 //------------- tidy up dual mode and pipeline -----------------
tail:
          m12tsm.ss       g2,p2  , p3
           nop
          m12ttpa.ss      g3,p3  , p4      // now T=0
           nop
          m12tpm.ss       g4,p4  , p5      // 0+a1*p1  ; left DIM
          m12tpm.ss       g5,p5  , p6
          m12tpm.ss       g6,p6  , f0      // output 0*0 => junk
          m12tpm.ss       one,two, g1      // to set T = 2
          m12tpm.ss       a1,g1  , g2
          m12tpm.ss       a2,g2  , g3
          m12ttpa.ss      a3,g3  , g4      // now T=2
          m12tsm.ss       a4,g4  , g5      // 2-a1*g1
          m12tsm.ss       a5,g5  , g6
          m12tsm.ss       a6,g6  , f0      // output 1*2 ==> junk

          m12tsm.ss       f0,f0  , p1      // to set T=0 ; p1 is 2-a1*g1
          m12tsm.ss       g1,p1  , p2
          m12tsm.ss       g2,p2  , p3
          m12ttpa.ss      g3,p3  , p4      // now T=0
          m12tpm.ss       g4,p4  , p5
          m12tpm.ss       g5,p5  , p6
          m12tpm.ss       g6,p6  , f0      // output 0*0 => junk
          m12tpm.ss       f0,f0  , g1      // flush results now
.align 8
        d.m12tpm.ss       f0,f0  , g2
        d.m12tpm.ss       f0,f0  , g3
        d.m12tpm.ss       f0,f0  , g4
           fst.d g1,8(resptr)++
          m12tpm.ss       f0,f0  , g5
           fst.d g3,8(resptr)++
          m12tpm.ss       f0,f0  , g6
// to fix pfld 4 to 4-24, scalar loads at least 24 bytes so no read off end
           adds        -20,ldptr,ldptr     // for auto load
          fst.d g5,8(resptr)++             // left DIM
          adds         4,resptr,resptr    // for auto store 32 bit

 scalar:                                  // if tail do at least 6 here
          adds 18,count,count             // count + 18 to do ( 0 to 11)
 scaloop: adds -1,count,count             // test at start so 0 pass is OK.
           bc Exit                        // exit if count <= 0
          fld.l    4(ldptr)++,a1          // load one float only
          frcp.ss  a1,g1                  // g ~ 1/a
          fmul.ss  a1,g1,p1               // g*a
          fsub.ss  two,p1,p1              // 2 -ga    = p
          fmul.ss  g1,p1,g1               // g*(2-ga) ~ g'
          fmul.ss  a1,g1,p1               // next iteration
          fsub.ss  two,p1,p1
          fmul.ss  g1,p1,g1
          fst.l    g1,4(resptr)++         // result 1/a
           br scaloop
          nop
Exit:
          fld.d  0(sp),f12                 // restore registers
          fld.d  8(sp),f14

        bri        r1
        adds 64,sp,sp

.align 8
         .data
ctwo:   .float  0f2.0
        .float  0f1.0