//---pgauss.s-----------------------------------------------
//    24-4-91
// 25 -4-91 tested DL: ave,var,var of non-tail,etc all OK
// 26-4-91 fix nr < size of temstr
// 14-5-91 tested Maxwell - reproducible with large stats
// Timing: 37 cycles/pass (n=5000), disable tail&not_short: 30 cycles
//
// Gaussian random number generator : exp(-x*x*0.5)
// using Ahrens Dieter Math Comp 27 (1973) 927 method.
// This is suitable for i860 since no division, uses look-up tables
// and has no subroutine calls.  Version FL5 is used which has
// 1.23156 flat random number calls on average. These follow the
// algorithm of pranmar.s: integers (0..2^32-1) are obtained from
// a call to pranmari which puts n values into resptr, the (variable)
// number of extra floating point random numbers come from calls to pranmar
// and are held in temstr.

// pgauss(resptr,n,uptr,temstr,nr) calls
// pranmari(resptr,n,uptr)once and then repeatedly (if necessary) calls
//  pranmar(temstr,nr,uptr)  and
//  pgauss1(resptr',n',uptr,temstr,nr) - nleft returned as n'-

// temstr should be nr (approx=n/4+1) words long. For uptr see pranmar.
// increasing nr helps to avoid "retries" but costs more to do. - too small
// a value of nr could give some bias because retries will come from large
// float rand parts of code - though rand int (first random decision) is reused
// in retries.

// pgauss1 tries to evaluate n gaussian random deviates
// from n pranmari random integers  and  a further nr+n/4+1 random floats
// after calling pranmar. If this is not enough then pranmar and pgauss1
// are called again with n=nleft to make up the deficit,etc.

// integer registers
        resptr  =r16
        n       =r17        // number gaussians left to do
        uptr    =r18
        temstr  =r19         // start of float rand numbers
        nr      =r20
        nleft   =r21
        nword   =r22
        ndone   =r23
        ntem    =r24



_pgauss_::                            // fortran load point
        ld.l    0(n),   n               // pointer to value
        ld.l    0(nr),  nr              // pointer to value

_pgauss::                             // C entry

     adds -1,n,r0         // test zero-pass
      bnc  go_on
     bri r1
     nop

// call pranmari(resptr,n,uptr)  use resptr as store
// call pranmar(temstr,nr,uptr)  use temstr as store
// call pgauss1(resptr,n,uptr,temstr,nr)   (nleft returned as r17)
//               r16  r17 r18  r19   r20

go_on:
      adds -32,sp,sp
      st.l n,0(sp)     // save  arguments
      st.l resptr,4(sp)
      st.l temstr,8(sp)
      st.l nr,12(sp)
      st.l r1,16(sp) // return address

        call _pranmari   // preserves uptr, regs 2..16
         nop

      ld.l  12(sp),nr
      ld.l  8(sp),temstr
      ld.l  4(sp),resptr
      ld.l  0(sp),n
      st.l  nr,20(sp)    // size of temstr

redo:
      st.l n,0(sp)      // save  n,resptr  (temstr not changed)
      st.l resptr,4(sp)
      mov  nr,r17      // in r17 :  number of float32 randoms estimated
      st.l nr,12(sp)
      mov temstr,r16   // reuse temstr for any extras needed

        call _pranmar    // preserves uptr, regs 2..16
         nop

      ld.l  12(sp),nr
      ld.l  8(sp),temstr
      ld.l  4(sp),resptr
      ld.l  0(sp),n

        call _pgauss1    // preserves uptr, regs 2..16
          nop

      mov  r17,nleft
      andnot  r0,nleft,r0
       bc EXITG          //CC set if zero

      ld.l  8(sp),temstr
      ld.l  4(sp),resptr
      ld.l  0(sp),n


      shr 2,nleft,nr
      adds 4,nr,nr       // new estimate of float rands needed (was 1)
      ld.l 20(sp),ntem  // size of temstr
      subs  ntem,nr,r0
       bc.t  too_big
        mov ntem,nr   // skip if nr <= size of temstr
too_big:

      subs n,nleft,ndone
      shl 2,ndone,nword     // number words done
      adds nword,resptr,resptr   // new address of results & rand_int

      mov nleft,n      // new number to do
      br redo     // or resptr=resptr+4
        nop

EXITG:

     ld.l  16(sp),r1
     bri r1
       adds  32,sp,sp

//----------------pgauss1.s------------------------------------------

// pgauss1(resptr,n,uptr,temstr,nr)
// resptr has n int32 random numbers in it, temstr is of size nr
// and has float32 random numbers in it. nleft is returned as n(r17)
// if temstr did not have enough numbers. uptr is unchanged.

// STRATEGY: (note fst increments resptr for use by ld.l rand-int )
//           bla dummy
//
//  NEXT_ONE:ldi
//           bla (not) left_overs             not_short_cut: refresh?
//                                                           fst
//  TOP:     bc tail       ldi                               test NEXT_ONE
//           bc not_short cut                                Quit
//            fst
//           bla  TOP                                tail:   refresh?
//                                                           fst
//  left_overs:            bc tail                           test NEXT_ONE
//                         bc not_short_cut                  Quit
//                         fst



// integer registers  (some defined in _pgauss )
  //    resptr  =r16
 //     n       =r17        // number gaussians left to do
  //    uptr    =r18
  //    temstr  =r19         // start of float rand numbers
  //    nr      =r20         // number float rands left
        nttable =r21         // a/h -t  -1.0  ; t+1.0    0..31  (0 tail)
        hstable =r22         // h (0..63) inc sign ; s (+/-1.0) 0..63
        datable =r23         // a-a ,  a (32 needed tail) 0..32  (0 tail)
        dtable  =r24         // d tail                   0..24
        count   =r17        // == n number gaussians left to do
        step    =r25        // bla decrement step
        ru      =r26
        ri      =r27
        ris     =r28
        rnd1    =r29
        ruf     =r30      // r31 pointer

// floating registers
     fn   = f16                     // double
     ft   = f17
     fd   = f18                     // double
     fa   = f19
     fh   = f20                     // double
     fs   = f21
     one  = f22
     half = f23
     ff   = f24
     fr   = f25
     fr1  = f26
     fu   = f27
     fw   = f28
     fv   = f29


_pgauss1_::                             // fortran load point
        ld.l    0(n),   n               // pointer to value
        ld.l    0(nr),   nr             // pointer to value

_pgauss1::                             // C entry


      orh   ha%const,r0,r31
      fld.d l%const(r31),one // one, half
      orh   ha%ntlable,r0,r31
      adds  l%ntlable,r31,nttable
      orh   ha%dalable,r0,r31
      adds  l%dalable,r31,datable
      orh   ha%hslable,r0,r31
      adds  l%hslable,r31,hstable
      orh   ha%dlable,r0,r31
      adds  l%dlable,r31,dtable

        adds    -4,resptr,resptr        // for auto increment load/store
        adds    -4,temstr,temstr        // for auto increment load


      adds   -1,n,count        // for bla n=1 no TOP loop, etc
      adds   -1,r0,step     // count -1 => -2 for bla fall thru.
.align 8
      bla    step,count,NEXT_ONE
        nop

.align 8
NEXT_ONE:       // gets here from not_short_cut or tail after bla test
    d.pfadd.ss f0,f0,f0
      ld.l   4(resptr),rnd1             // from pranmari call in pgauss
    d.pfadd.ss f0,f0,f0
      andnoth  64512,rnd1,ru            // mask not 111111000000..
    d.fnop                        // in DIM
      shr      3,ru,ru                  // for fp use: msb 9 to right
    d.fnop
      xorh     16256,ru,ruf             // 0 01111111 ru = 1.ru
    d.fnop
      ixfr     ruf,fu                   // fu is 1.0 to 2.0 random flat
    d.fnop
      andh     64512,rnd1,ris           // mask 11111100..=> int 0..63
    d.fnop
      shr      23,ris,ris               // to 2word address 0..63 *8  (32-6-3)
   d.fnop
        bla step,count,TOP           // test if >= 2 left to finish
    d.fnop
      and      248,ris,ri                // mask 0000 11111000: ri 0..31*8

                                // here if last go only
    d.fnop
      br left_overs            // count = -2 here
    d.fnop
      nop

//------int rand number--to-----fast-gaussian--code---------------

TOP:              // here at least 2 still to complete, count >= -1
    d.fnop
      fld.d    nttable(ri),fn           // fn,ft loaded (n,t+1)
    d.fnop
      ld.l   8(resptr),rnd1             // from pranmari call in pgauss
    d.fnop                              // 8 since fst not yet done.
        bc  tail                        // 3% only,  ris/8 = 0 or 32, ru OK
    d.fnop
                                    // (ai + (rnd - ti)*hi) => (ni+rnd)*hi
      fld.d    hstable(ris),fh          // h contains sign (ris/8= 1..63,not32)
    d.pfadd.ss fn,fu,f0                 // n+1.ru  (so just take 1 from n )
      andh     64512,rnd1,ris           // mask 11111100..=> int 0..63
    d.pfadd.ss f0,f0,f0
      andnoth  64512,rnd1,ru            // mask not 111111000000..
    d.pfle.ss  fu,ft,f0                 // in adder pipe
      shr      3,ru,ru                  // for fp use: msb 9 to right
    d.m12tpm.ss fh,ff,ff                // h*(n+rnd)  [care T=0]
      shr      23,ris,ris               // to word address 0..63 *8  (32-6-3)
    d.pfmul.ss f0,f0,f0
        bnc   not_short_cut             // u < t CC clear: not_short_cut
    d.pfmul.ss f0,f0,f0                 //               ri fu fs ft fn fh used
      xorh     16256,ru,ruf             // 0 01111111 ru = 1.ru
    d.pfmul.ss f0,f0,ff
      ixfr     ruf,fu                   // fu is 1.0 to 2.0 random flat
    d.fnop
      fst.l    ff,4(resptr)++     // resptr is same as rndptr in effect
    d.fnop
        bla step,count,TOP        // test if >= 2 left to finish
    d.fnop
      and      248,ris,ri                // mask 0000 11111000: ri 0..31*8
//-----------------------------------------------------------------

left_overs:   // end of fast loop for last case. count = -2
    d.fnop
      fld.d    nttable(ri),fn           // fn,ft loaded (n,t+1)
    d.fnop
        bc  tail                        // 3% only,  ris/8 = 0 or 32, ru OK
    d.fnop
      fld.d    hstable(ris),fh          // h contains sign (ris/8= 1..63,not32)
    d.pfadd.ss fn,fu,f0                 // n+1.ru  (so just take 1 from n )
      nop
    d.pfle.ss  fu,ft,f0                 // keep adder pipe going
      nop
    d.pfadd.ss f0,f0,f0                 // (ai + (rnd - ti)*hi) => (ni+rnd)*hi
      nop
    d.m12tpm.ss fh,ff,ff                // h*(n+rnd)  [care T=0]
        bnc   not_short_cut             // u < t CC clear: not_short_cut
    d.pfmul.ss f0,f0,f0
      nop
      pfmul.ss f0,f0,f0
      nop
      pfmul.ss f0,f0,ff
      fst.l    ff,4(resptr)++     // resptr is same as rndptr in effect

                       // all done in orderly fashion
       bri r1
      adds r0,r0,n     // none left to do

//------------------------------------------------------------------
refresh:         // nr ran out so return with n gaussians left to do in r17
      bri r1      // count=n-1, then 2 bla before store
      adds 3,count,n   // n  left to do:  CARE to reuse random int

//------------------------------------------------------------------

.align 8
             // 7.5% only, use ri(double word),fu (1.ru),fs,fn,ft,fh
not_short_cut:         // one rnd call (fp); called DIM and not DIM.
   d.pfsub.ss  fu,one,f0                // 1.ru - 1 => rnd
      adds     -1,nr,nr
   d.pfsub.ss  f0,f0,f0
      fld.d    datable(ri),fd            // fd,fa loaded
     pfsub.ss  f0,f0,f0
      nop
     pfsub.ss  f0,f0,fu       // u < tmax so test if u < v ( v < tmax)
      bc   refresh                       // refresh taken rarely
      fld.l    4(temstr)++,fr       // left DIM
      fmul.ss  fd,fr,fw            //  w=rnd2*d, new trial x value
      fmul.ss  half,fw,ff
      fadd.ss  fa,ff,ff
      fmul.ss  fw,ff,fv                 // v=(w/2+a)*w (estimate ft above)
odd:  pfle.ss  fu,fv,f0                 // u < v  CC clear:  even
        bnc  even                       // < 50% jumps here
      fadd.ss  fw,fa,ff                 //  (w+a)
      adds 1,count,r0               // if count = -2 QUIT
      fmul.ss  fs,ff,ff                 // mult by sign: loaded with h(ris)
      bnc.t  NEXT_ONE                // mostly taken
     fst.l    ff,4(resptr)++

    fst.l    ff,4(resptr)++
        bri r1
      adds  r0,r0,n
//-----------------------------------------------------------------

// -- 2 rnd calls per odd/even iteration-------
 even:
      adds     -2,nr,nr
       bc    refresh                   // if nr -1 or -2 refresh
      fld.l    4(temstr)++,fr1
      pfle.ss  fr1,fu,f0              // fr1 < fu  CC clear: next_odd
      fld.l    4(temstr)++,fu
        bnc.t odd
      fmov.ss  fr1,fv                  // skipped if restart

// even termination so re do from start (same ri though)
      fadd.ss one,fu,fu               // to match 1.ru convention in fastloop
      fadd.ss fn,fu,ff                // code as fast part:fn,ri,ft,fh used
      pfle.ss fu,ft,f0                // 1.ru < ti+1 CC clear: not_short_cut
        bnc  not_short_cut             //  use ri, fu (1+rnd), etc
        adds 1,count,r0               // QUIT if count = -2
      fmul.ss fh,ff,ff                // h*(n+1.rnd)
        bnc.t  NEXT_ONE
    fst.l    ff,4(resptr)++

    fst.l    ff,4(resptr)++
        bri r1
      adds  r0,r0,n


//------------------------------------------------------------------

.align 8
tail:    // assume ru,ris exist (ris/8 = 0 or 32)  in DIM
                              // ru now clear ms 9 bits
     fnop
       adds  r0,r0,ri
     fnop
       fld.d 256(datable),fd   // fd,fa 8*32 up - a of end of centre
test_0:  // find first zero bit of ru - so always terminates
       shl  1,ru,ru
       andh 128,ru,r0       // test msb   mask 00000000 10000000 000..
         bc found_0          // CC set if msb zero
       fld.l dtable(ri),fd
       adds 4,ri,ri         // ri is look-up word address
        br test_0
       fadd.ss fd,fa,fa

found_0:
       fld.l  dtable(ri),fd
       andnoth 65408,ru,ru    // 9 msb bits cleared: mask 11111111 10000000
       xorh    16256,ru,ru    // 1.ru
       ixfr   ru,fr
       fsub.ss fr,one,fr      // fr 0..1 random

tail_even:
       fmul.ss fd,fr,fw       // w = di*rand
       adds   -1,nr,nr
       fmul.ss fw,half,ff
        bc    refresh
       fadd.ss fa,ff,ff
       fld.l  4(temstr)++,fr
       fmul.ss fw,ff,ft


tail_odd:
       pfle.ss fr,ft,f0      // r < t CC clear: keep going
        bc done_tail
       adds -2,nr,nr
        bc  refresh
       fld.l  4(temstr)++,fr1
       pfle.ss fr1,fr,f0    // r1 < r CC clear: keep going- tail_odd
       fld.l  4(temstr)++,fr
         bc     tail_even      // restart random choice of w
         br     tail_odd
       fmov.ss fr1,ft


done_tail:
       fld.d   hstable(ris),fh   // 0 or 32 only   +/- 1.   fh,fs
       fadd.ss fa,fw,ff
       fmul.ss fs,ff,ff
         adds   1,count,r0        // QUIT if count = -2
      fst.l    ff,4(resptr)++
        bnc  NEXT_ONE

        bri r1
      adds  r0,r0,n


.data
.align 8
const:
     .float   1.000000000   // one
     .float   0.500000000   // half
//     a1..32: READ in 22-4-91 by CM from Ahrens Dieter
//  d-a 0..32        gauss.num
dalable:
  .float 0f0.0              // tail i=0
 .float 0f0.0
  .float 0f0.0391760850     // d[1] = a[2]-a[1]
 .float 0f0.0000000000      //a[1]
  .float 0f0.0392363280
 .float 0f0.0391760850
  .float 0f0.0393574615
 .float 0f0.0784124130
  .float 0f0.0395408101
 .float 0f0.1177698745
  .float 0f0.0397883996
 .float 0f0.1573106846
  .float 0f0.0401030251
 .float 0f0.1970990842
  .float 0f0.0404883305
 .float 0f0.2372021093
  .float 0f0.0409489241
 .float 0f0.2776904398
  .float 0f0.0414905278
 .float 0f0.3186393639
  .float 0f0.0421201736
 .float 0f0.3601298917
  .float 0f0.0428464596
 .float 0f0.4022500653
  .float 0f0.0436798862
 .float 0f0.4450965249
  .float 0f0.0446332951
 .float 0f0.4887764111
  .float 0f0.0457224560
 .float 0f0.5334097062
  .float 0f0.0469668501
 .float 0f0.5791321622
  .float 0f0.0483907378
 .float 0f0.6260990123
  .float 0f0.0500246333
 .float 0f0.6744897501
  .float 0f0.0519073777
 .float 0f0.7245143834
  .float 0f0.0540891171
 .float 0f0.7764217611
  .float 0f0.0566356808
 .float 0f0.8305108782
  .float 0f0.0596351973
 .float 0f0.8871465590
  .float 0f0.0632084129
 .float 0f0.9467817563
  .float 0f0.0675253978
 .float 0f1.0099901692
  .float 0f0.0728338133
 .float 0f1.0775155670
  .float 0f0.0795093789
 .float 0f1.1503493803
  .float 0f0.0881521381
 .float 0f1.2298587592
  .float 0f0.0997862406
 .float 0f1.3180108973
  .float 0f0.1163234064
 .float 0f1.4177971379
  .float 0f0.1418191784
 .float 0f1.5341205443
  .float 0f0.1867921447
 .float 0f1.6759397227
  .float 0f0.2911428266
 .float 0f1.8627318674    // a[31]
  .float 0f0.0000000000
 .float 0f2.1538746940    //a[32]    - used start of tail
 ntlable:                // t = t(a-d)+1.0 ;  n=a/h-t (a+(r-t)*h)=(n+1.r)*h
 .float 0f0.0            // i=0  tail
 .float 0f0.0
 .float 0f-1.0007673828
 .float 0f1.0007673828
 .float 0f-0.0061455873
 .float 0f1.0023068704
 .float 0f0.9807616010
 .float 0f1.0038606184
 .float 0f1.9568020299
 .float 0f1.0054384541
 .float 0f2.9187552126
 .float 0f1.0070506988
 .float 0f3.8633097841
 .float 0f1.0087083958
 .float 0f4.7870398881
 .float 0f1.0104235698
 .float 0f5.6863784172
 .float 0f1.0122095319
 .float 0f6.5575872964
 .float 0f1.0140812473
 .float 0f7.3967234622
 .float 0f1.0160557881
 .float 0f8.1995995924
 .float 0f1.0181529007
 .float 0f8.9617377739
 .float 0f1.0203957318
 .float 0f9.6783151890
 .float 0f1.0228117673
 .float 0f10.3440978833
 .float 0f1.0254340733
 .float 0f10.9533608968
 .float 0f1.0283029560
 .float 0f11.4997891054
 .float 0f1.0314682249
 .float 0f11.9763530027
 .float 0f1.0349923344
 .float 0f12.3751508302
 .float 0f1.0389548297
 .float 0f12.6872037861
 .float 0f1.0434587838
 .float 0f12.9021858427
 .float 0f1.0486403492
 .float 0f13.0080591596
 .float 0f1.0546833385
 .float 0f12.9905698868
 .float 0f1.0618422239
 .float 0f12.8325295105
 .float 0f1.0704798276
 .float 0f12.5127550545
 .float 0f1.0811319498
 .float 0f12.0044367580
 .float 0f1.0946244354
 .float 0f11.2724877098
 .float 0f1.1123000789
 .float 0f10.2689326468
 .float 0f1.1364979994
 .float 0f8.9241092257
 .float 0f1.1716885601
 .float 0f7.1275068856
 .float 0f1.2276240548
 .float 0f4.6764191867
 .float 0f1.3304980279
 .float 0f1.0723666016
 .float 0f1.5847030938
//hs
hslable:
  .float 0f0.0          // tail i=0
 .float 0f1.0           // s[0] used tail
  .float 0f0.0392061716  // from table
 .float 0f1.0000000000
  .float 0f0.0393270496  // from table
 .float 0f1.0000000000
  .float 0f0.0395099949  // from table
 .float 0f1.0000000000
  .float 0f0.0397570269
 .float 0f1.0000000000
  .float 0f0.0400709276
 .float 0f1.0000000000
  .float 0f0.0404553261
 .float 0f1.0000000000
  .float 0f0.0409148089
 .float 0f1.0000000000
  .float 0f0.0414550711
 .float 0f1.0000000000
  .float 0f0.0420831105
 .float 0f1.0000000000
  .float 0f0.0428074814
 .float 0f1.0000000000
  .float 0f0.0436386273
 .float 0f1.0000000000
  .float 0f0.0445893180
 .float 0f1.0000000000
  .float 0f0.0456752278
 .float 0f1.0000000000
  .float 0f0.0469157137
 .float 0f1.0000000000
  .float 0f0.0483348698
 .float 0f1.0000000000
  .float 0f0.0499629842
 .float 0f1.0000000000
  .float 0f0.0518385865
 .float 0f1.0000000000
  .float 0f0.0540113819
 .float 0f1.0000000000
  .float 0f0.0565465619
 .float 0f1.0000000000
  .float 0f0.0595313042
 .float 0f1.0000000000
  .float 0f0.0630848897
 .float 0f1.0000000000
  .float 0f0.0673750349
 .float 0f1.0000000000
  .float 0f0.0726454356
 .float 0f1.0000000000
  .float 0f0.0792647141
 .float 0f1.0000000000
  .float 0f0.0878192233
 .float 0f1.0000000000
  .float 0f0.0993039833
 .float 0f1.0000000000
  .float 0f0.1155599414
 .float 0f1.0000000000
  .float 0f0.1404343835
 .float 0f1.0000000000
  .float 0f0.1836141833
 .float 0f1.0000000000
  .float 0f0.2790016348
 .float 0f1.0000000000
  .float 0f0.7010474248
 .float 0f1.0000000000
  .float 0f0.0          // tail i=32     h s -ve
 .float 0f-1.0           // s[32] used tail
  .float 0f-0.0392061716
 .float 0f-1.0000000000
  .float 0f-0.0393270496
 .float 0f-1.0000000000
  .float 0f-0.0395099949
 .float 0f-1.0000000000
  .float 0f-0.0397570269
 .float 0f-1.0000000000
  .float 0f-0.0400709276
 .float 0f-1.0000000000
  .float 0f-0.0404553261
 .float 0f-1.0000000000
  .float 0f-0.0409148089
 .float 0f-1.0000000000
  .float 0f-0.0414550711
 .float 0f-1.0000000000
  .float 0f-0.0420831105
 .float 0f-1.0000000000
  .float 0f-0.0428074814
 .float 0f-1.0000000000
  .float 0f-0.0436386273
 .float 0f-1.0000000000
  .float 0f-0.0445893180
 .float 0f-1.0000000000
  .float 0f-0.0456752278
 .float 0f-1.0000000000
  .float 0f-0.0469157137
 .float 0f-1.0000000000
  .float 0f-0.0483348698
 .float 0f-1.0000000000
  .float 0f-0.0499629842
 .float 0f-1.0000000000
  .float 0f-0.0518385865
 .float 0f-1.0000000000
  .float 0f-0.0540113819
 .float 0f-1.0000000000
  .float 0f-0.0565465619
 .float 0f-1.0000000000
  .float 0f-0.0595313042
 .float 0f-1.0000000000
  .float 0f-0.0630848897
 .float 0f-1.0000000000
  .float 0f-0.0673750349
 .float 0f-1.0000000000
  .float 0f-0.0726454356
 .float 0f-1.0000000000
  .float 0f-0.0792647141
 .float 0f-1.0000000000
  .float 0f-0.0878192233
 .float 0f-1.0000000000
  .float 0f-0.0993039833
 .float 0f-1.0000000000
  .float 0f-0.1155599414
 .float 0f-1.0000000000
  .float 0f-0.1404343835
 .float 0f-1.0000000000
  .float 0f-0.1836141833
 .float 0f-1.0000000000
  .float 0f-0.2790016348
 .float 0f-1.0000000000
  .float 0f-0.7010474248
 .float 0f-1.0000000000
  //   d:  5..32 of Ahrens Dieter
dlable:
      .float  0f0.2636843221   //   d6 = f(2^-6)-f(2^-5) so start tail
      .float  0f0.2425084523
      .float  0f0.2255674438
      .float  0f0.2116341657
      .float  0f0.1999242674
      .float  0f0.1899107584
      .float  0f0.1812251810
      .float  0f0.1736014003
      .float  0f0.1668419086
      .float  0f0.1607967291
      .float  0f0.1553497174
      .float  0f0.1504093838
      .float  0f0.1459025768
      .float  0f0.1417700327
      .float  0f0.1379631736
      .float  0f0.1344417615
      .float  0f0.1311721502
      .float  0f0.1281259651
      .float  0f0.1252790900
      .float  0f0.1226108828
      .float  0f0.1201035596
      .float  0f0.1177417070
      .float  0f0.1155118922
      .float  0f0.1134023487
      .float  0f0.1114027204
      .float  0f0.1095038520
      .float  0f0.1076976165