From 1108268557b784623588e537b4c97da1434345f9 Mon Sep 17 00:00:00 2001 From: James Stine Date: Mon, 20 Jun 2022 11:32:40 -0500 Subject: [PATCH] Update C program for r=4 division by recurrence to match Table in EL --- pipelined/srt/stine/notes | 13 ++ pipelined/srt/stine/srt4div | Bin 27024 -> 26896 bytes pipelined/srt/stine/srt4div.c | 231 ++++++++++------------------------ 3 files changed, 79 insertions(+), 165 deletions(-) create mode 100644 pipelined/srt/stine/notes diff --git a/pipelined/srt/stine/notes b/pipelined/srt/stine/notes new file mode 100644 index 000000000..32cb772c6 --- /dev/null +++ b/pipelined/srt/stine/notes @@ -0,0 +1,13 @@ +Dividend x --(0.10101111), divisord --(0.11000101)(i -- 16(0.1100)2- 12) + +X = 175 (xAF) +D = 197 (xC5) + +X = 175/256 = 0.68359375 +D = 197/256 = 0.76953125 + +Note: Add lg(r) extra iterations due to shifting of computed q + q_{computed} = q / radix + +./srt4div 0.68359375 0.76953125 4 10 + diff --git a/pipelined/srt/stine/srt4div b/pipelined/srt/stine/srt4div index ff76638fae90ef9bb09e93e26e1f5f7e79619555..b44c6dc0d18016698d3b7391786f51596957c3b8 100755 GIT binary patch delta 6234 zcmZ`-4Rlo1wLa(0oXnlcOlBs3cP0tRBoZ2eW|DxWK|{jN)JY{0Tm}ZA0wILNhd_9s z_9;#%h$KVWgx$K($E%U$>ytj}0z;oEicV@P!KIk0mEd}6tPeZbq9j@ct(~{eIrk2e zEPJwY_TRVn*=L`hd*}APNqgU*y`JpG6$h zBMySqv^leW&h*buK6lri=K|G_Wmn(w{;}@XU%)<7iAdi_j?kTiP?y4|lnJkyW8#iyLgshZ4IGLis&r`H8ViaN`WF00-9Zrfp{~mmG9wpt`FP(WxFBZor7myJX~%WNk$5$*t)MvQ>&WOUNF*Rd z0JyG81W-tt?qDpb%Du5g^_|hEFQTFEgu|9u@K=l7F}6e^`NlBAMT_Q5l?u^cO0v3O znyPs-lB`B48S2uZLDd|YN!BgUftvgD2o|&`utU}0h-!E2P12^ouTcBHs6CLRRpDLK z+O2lSQ%PDIIIJ4ngxbL*tpjRMyAid+Nm?h|ZPloIpwGH6WxuOPom)%DM7;P$D0RJ( z*=EOlCN31lZY_CK6|D7wuWA9EITa^OokMf_=jzy7x{ zD?N`s1i|$DnI}<2w$Ew`?y70#Y_RM6HhVBq(~3fD#XsWlSPM?%W%e`pN%~{-B`D3f zhc1REu>T`?G9Soq%O#}MdHc|W(;XOz1jfD`ZFPPcop6>$F_|0Tnano&2IOYd(*3Y9 z>q+X^pPi@t-3$5I*+qV~37xk-IN@}$jUtRcFm}&o6a%qVr=P9Kh50$zN@Fez>1JH)_98bm zg~c*RiwJVN46=v^x%Npsl8GQJ5hV9pg1iWIIhu0TBc^f6BZ}~>M@-|BX=w!qJVpkO zAU8FGO5yx{={#3DmrCdR^V578-1&aKdhUet#>DFVUK8~{sqvxduV4%2fx>uxark$g zT{{j2ySDeX1S96+60z$7mn)*D0)zbdvTMgdbJ!n@1bTxJ7w)1Z*YjQOXw2ff^GvYo z)F6UL?7lDJ@fZ2!6nh#u0pAr>gAv3wMG@aL1>a_%!n>jvmcxb-=AQyvV8C`oJrZ>U zQI097ysIdmL>)!cpYf7PPVtMs#p98{(W@sS5N&mqfWv1Aw|p6ow+5XLMuS~}=w!?= zVusPSDH@xXgu-Om!(srv0Lw6XFp}k*LrA3IC>9{vb>>AA7DW1PxQz9&@%mi+9DzM- zNs@Co=9S)W%r?A-CwZ49ypN){oPj&Q5BC>%xVI1T+%L(2!W_V-@|p8N_FF6zS0i#c z{+GOj3s}Z=i8*y5e8S~)HQ+6=Lccefj=+b$nc@BT2;l3>ClPDnINV%6KyQA#jGYP` zM?JaFk7o$;Kx63+k;i43$6KEz<)P_$jA0)5U-CQ5V;2_SY94IufBYNLDU3BcmNJ5kijmXGR9{g3~&D@U31lL8fitr@}7jw#xyttEA=U;I9$+aDy5No?XIoLEFERywhg%RsCma5@- z%px!hV};o(aue=PNcZ|=_t&^ClC?*=NqD>-%x0w-q*hl&E><3^2#x?aY4 zzJPC;e4hM;Emnapk*4P(f#bsVs&3=!8q1QooJ!VFnahcCW0#Ij$|ap|;vV+z$u8kD z<0Aiu6F!ske;@Z1`fQdyu`P(+6*x{Kfg{y{mor($Z)CR{E}a&bzR z$DKH}mabE=egVg`$1g>UW!Gbg+kb!?dzhLW|8W^Vl!!mM3`gZotK;!8{QOU0s8~~i zpANiJJi|87Z+*-3O!*la_SNxqAf%{{_k?)r1R=zrQ7L2!upeY(9wmrHf( zGN^8!uG4F!K1|AVK)Fs!bO7ljliDSG|Kv@rqq(LVTae>#+Wqgw;|Gv-z88;2k@g}T zMY`twc)Sy7{a8GH3F&_$bzv!sK8VN5kv@vF5ox;XPx1H`WOjUrud+zz#9-O%QaTDd zXKRdL^fCC4+1JoP7@DnBA4PctqZus^nSJ=b6h~|3%0QZHkFwB_qO7o6@O_rskKjK* z#hhIBDsmqqS5P<~v$+?gtN2X`yd9G`03CA*s?KAM&+?ZABdY0XGl%i0ixKs}NAdVF z{=*<)Z?K?!V7&N>I@9>+(uuboZ-^J|b)6dWIrw2sf!kgr)%bQLS`ZG`6wpREUE{}6 zTtVIo-g$m^Xo<8dI_|r$VqSr}7pKY_&xN;s=$f~L`rzEW0y+U#=53%MSX+B9y##-# z-A+%#vAPU+PO(F4oy+zOKErr(*x{Ix4Ug32(9;9Q>Rwi8FWk5AUlngATwK^h55dAk z%jhh4Y*7dUe75L1`U(^Vyh-i@ivk%+Sp5WJnX52A^KX!|^S!YQCC6N3Q0Yfr_gz17|BA!*4gvFI&Gr`!+TG+aE(|TgDiP*(BU?Mgl-oKzz z%Cpq|C#3T+{`a^&(blVeajzc(>xV}cOJ1OTMOPq~eUu#wzyBy=z1F30gF0B>T$A*Ba@gj+Lp1|s&Rae%3BUpO_7i#gn4>gx=% zM$(-sO4C=zZYAAmoHz9y4YNh0+bxB`Fk4ExGbJj>ZVRKNNh!o`LtL&D*5noPWV#Eu z$j=t(E|!)?)*7;k{QZtm zYo^WOE8QktDB=s47yzkPg zkik6(AKa1&X^Xd384)}paGDWVhV>?ttU7g5+ST!LpntJOb0O`PJXfB`RwTqjPeY{# zE07QI#T|ImIvQTH&PJILl{CZ6w^mXc{M)T1MZ=h}@)dSwRz$3dQ7Mg?)E}tNu@v=r z*tw(u>KZ??*~-*1rhd#U_rP^a>(*rVF}r$#*_6_BmJ(v-KIT+vGFkR8GjC%)qL@v0 zuxw5cWapFvZJc?BV9VFW#+R3 z>u&o>ksrENHO*vO*S56ZH?X^DA%*A~GrYQNvDcN3MM-Kt5T8E?Z+D~2+eae9mqkAcAK)nG1%HL->{an#xfL%_r(~N2^Y73q zrlXZ->(yf#nOv*ujs0hw)Czb4>#%_L;Z=xn;it>SDUsjP^E5sY&yZ=7wL0E41%J0b z*3_$TtFE6lrAqe*J%Rta!%g>!d>)0hD>e3?sQ2O1E`z?vwQL0CA|d|;hE{5HeBj(l zCw20#sdYCmu3xZlb!Ay)^}wH3SsA?y<7@ryQU1vgn?b!$(d>6m@Q;w#3|a~u&3>y- ze=7UnaC4QlUzhsfLUUE(!EHRizY7)X!gLIdt@DGg#Ra>quKg`Ow+r9cc{At;RJ8b~ z8D%qQsY{|n4?C?H(A(m&8C%3&1ix(g9`66~7QglC=kJ;8t8B4K=@_eoyVl=eH9o@P zF<4&({p-uY)#|bu2O$JATdNZ9AtC*x@Tq~0)>?|^;n;>M>nHUx-u*DXp~}{-f4g08 zhRU`ZtY?-=Py(Xe7oUM45JILq8-@9)f9Ust+4{3)pHEbokR0u%=!Oj`KGISV(+R^!m*LX2aC>PTp zl@{Lfu3o5`Bo?7l6bX0*%7h-Mo+!hT&{c1O!JvD@o_{^`N#8F|Ui`_AfBup;M98)k(-hdR%XhIJYRM703l)pRY&LjLHxn%D z^$4ye#OX$P-e(>du0R-?Q5{c(QZ1!wd>1i$QrBH2j~a%500qOH&=JLJ{K zIu&{(8+{J?nqcT3b?BG7sKJ@L+Bd(P2sw(GwZ(?2qEuCE_J& z!p6pep-!V@gIt8yc6T$*O#(Tm+FCoG$3fSVuTW=&V*M;ibyA3)fd5A@+M*H+iE!F`;IN13qR&Q zjFjup$L(zMUYH-V`!ze7$%WfuTzAHEMWe9KWzMX`<;eg;cuO94v(9+?w(dktPJ{>2 zZPJ(& z!-D31L6a9xfxNpfQ(^KU`o6(qlgB`$~bYbVd-66|!j3PQLiW0@RxYMH!^waJ;h z;pCYt|m)l0`9+>9-&$<0R7S4mLydb!c%rq#F)4L9-g;S$^fvQaOvY z$*J^10ORrnNj@jpCrAc$ajDeANT%I_(}QDS`$CquH6xdITemWY4r%^c}Z|zDGpwO5;hOV9``WKjBRF9L3XGtkXfPO>X<$ zG~V{zB6^mH-k-+XHoB5MwXx{+w}>tm(XZKg+b`L9+fh44=X!?hTD$*(waM)swexnT zipW$Exs^pa3!cVp99s(+jOKAv<(wGJCg$batW3~6XX6PxW#b9lupx~Dy4@+EuQ#_o zam3%cy{p0Brq9hsVt=deT5m@&K{MrG!eXtdS$CJu)bHyn}2F7dGp6!Ni*wOFLK$XF(h&!cPs z^3nb7Z+vwBkV2fX;f>@AJ}UU2I_viOyja8$5wR`h=6pm3nXN6Yr$C#`feF*`C3ZzL z_Xq-G`EA8Muvvi9S-#M5t>=BMZ)5-dg6}G9y4hJBnTJtrwQsihj&t7?)rT$;nIh7; z6mK(#bUvS@r877o9XHzn?ebgkA!|Om$h+`bEM{UC{wh}CKZ%&g4n#lcJ5Jkt-Sd3M z&HJCZA(k0?9u1iHp2Ezqn{36^`-ULJlMxE2k@*^@Z6) z;c~sdk1s>a^|V-pY__=xcYlO|*z%QG2jOc^R_GZu1KyDumO;DR6`xg^$SOV7k;sJU zS;AIdySzO<#3(|(n4$H3WI{-W2>FjB?0ZA&smPHFtlz!?o-cCUsdDjW?n&lsm*>Ru zzsdvB_Qd!nj?n<9WKAOKe6l=wZob{EO1JNkQP;5QF*SYI7?Xg?W*4%J~mwc^s9&666+A~%Cc1^1GJXOzA^{J|!uSwNjsOp8PUZmQMRD$`ExWng3>k=jaqQ}o<>=Ol71GA22k!m zxtSu>@aNHJ2O8{O!?}_fQz|e%593H>B0dEPxgD=LcVGpqk98vOv@jBoniNkD&Z&eT?Ay7qH7sQx-dJTIqJsp)57 zOIfycC|8*Ag_@TDyj+$|>!7d9OT%y(>vpg%@LGdQg`aIjRy=e%P|o z0=pJ0#ZH`AkWJ6Rz9Px#tOxs zst$O*!c7nL_Efwr(ao@9(NCo7Hn_N`7QbjLs$4;n;Kj;1v*IU06Y*7wd`4qhG^nzrn%KH)4$gLmJRV^tWL>nHf^P z$6E)l8Sy;tWHocsrme(u6ls2on#3ynMnaByWK6Uqd>pgi!<_uaB8hqH!;QqqB1U2m zrV-~~P*35j2E#dS5>tuBsY_*slxY4$oL`ps;{ZegCYFTf;^cX zMy%(()Yb~$S?Hq&41*Y>PgOIT(tt>wg4D^{P3S#e^$iW#h%V)B+#ZS(V+=BCcjeIe zN%NAbY!aG`b0%AoD6O(l8xCs@sj_+T*_2$067ktFiN z{-e=}LofyEs<+_IIbHppfz6XG;FdTPRywA#}?Q+6X8W1)w^+4s)yuWw#RxUNsL_;A`CVmC0 z0~LNlr<_nH>lK@9mdfmMdcUmSCY!KX-D=s$gtAIbg}odR)X&3`Mu|sFQkyYqNaF+eT|JP9;{3FG18M7XT<7k)gFYGS9tNM zdS->gr9GcDdjcbaP z$AlpJnH@hvYbSTTdLZI&H?Lc2Q6Xu&|<_`h?@i^F~6_osr?s$pX``t zjf%HS!r!O1HTe;_PPHdZ%8mV~!@&i|YAfBr^J#@Ms~z%lX!r0B?XNMb97bJyc>fb=}pPVvQ~ctfe-aFsAN!IllBMoYEOHp_5y!(2nlGNI0Z%Nt7RVVJ(L+|a#3 e*eYPl#!} +// QSLC is for division by recuerrence for +// r=4 using a CPA - See Table 5.9 EL int qslc (double prem, double d) { int q; + // For Debugging printf("d --> %lg\n", d); printf("rw --> %lg\n", prem); - if ((d>=0.0)&&(d<1.0)) { - if (prem>=1.0) - q = 2; - else if (prem>=0.25) - q = 1; - else if (prem>=-0.25) - q = 0; - else if (prem >= -1) - q = -1; - else - q = -2; - return q; - } - - if ((d>=1.0)&&(d<2.0)) { - if (prem>=2.0) - q = 2; - else if (prem>=0.66667) - q = 1; - else if (prem>=-0.6667) - q = 0; - else if (prem >= -2) - q = -1; - else - q = -2; - return q; - } - - if ((d>=2.0)&&(d<3.0)) { - if (prem>=4.0) - q = 2; - else if (prem>=1.25) - q = 1; - else if (prem>=-1.25) - q = 0; - else if (prem >= -4) - q = -1; - else - q = -2; - return q; - } - - if ((d>=3.0)&&(d<4.0)) { - if (prem>=5.0) + + if ((d>=8.0)&&(d<9.0)) { + if (prem>=6.0) q = 2; else if (prem>=2.0) q = 1; else if (prem>=-2.0) q = 0; - else if (prem >= -5) + else if (prem >= -6) q = -1; else q = -2; return q; } - if ((d>=4.0)&&(d<5.0)) { - if (prem>=6.66667) + if ((d>=9.0)&&(d<10.0)) { + if (prem>=7) q = 2; else if (prem>=2.0) q = 1; else if (prem>=-2.0) q = 0; - else if (prem >= -6.66667) + else if (prem >= 7.0) q = -1; else q = -2; return q; } - if ((d>=5.0)&&(d<6.0)) { + if ((d>=10.0)&&(d<11.0)) { if (prem>=8.0) q = 2; else if (prem>=2.0) @@ -91,7 +53,21 @@ int qslc (double prem, double d) { return q; } - if ((d>=6.0)&&(d<7.0)) { + if ((d>=11.0)&&(d<12.0)) { + if (prem>=8.0) + q = 2; + else if (prem>=2.0) + q = 1; + else if (prem>=-2.0) + q = 0; + else if (prem >= -8.0) + q = -1; + else + q = -2; + return q; + } + + if ((d>=12.0)&&(d<13.0)) { if (prem>=10.0) q = 2; else if (prem>=4.0) @@ -105,21 +81,35 @@ int qslc (double prem, double d) { return q; } - if ((d>=7.0)&&(d<8.0)) { - if (prem>=11.0) + if ((d>=13.0)&&(d<14.0)) { + if (prem>=10.0) q = 2; else if (prem>=4.0) q = 1; else if (prem>=-4.0) q = 0; - else if (prem >= -11.0) + else if (prem >= -10.0) q = -1; else q = -2; return q; } - if ((d>=8.0)&&(d<9.0)) { + if ((d>=14.0)&&(d<15.0)) { + if (prem>=10.0) + q = 2; + else if (prem>=4.0) + q = 1; + else if (prem>=-4.0) + q = 0; + else if (prem >= -10.0) + q = -1; + else + q = -2; + return q; + } + + if ((d>=15.0)&&(d<16.0)) { if (prem>=12.0) q = 2; else if (prem>=4.0) @@ -133,106 +123,9 @@ int qslc (double prem, double d) { return q; } - if ((d>=9.0)&&(d<10.0)) { - if (prem>=15.0) - q = 2; - else if (prem>=4.0) - q = 1; - else if (prem>=-4.0) - q = 0; - else if (prem >= -15.0) - q = -1; - else - q = -2; - return q; - } - - if ((d>=10.0)&&(d<11.0)) { - if (prem>=15.0) - q = 2; - else if (prem>=4.0) - q = 1; - else if (prem>=-4.0) - q = 0; - else if (prem >= -15.0) - q = -1; - else - q = -2; - return q; - } - - if ((d>=11.0)&&(d<12.0)) { - if (prem>=16.0) - q = 2; - else if (prem>=4.0) - q = 1; - else if (prem>=-4.0) - q = 0; - else if (prem >= -16.0) - q = -1; - else - q = -2; - return q; - } - - if ((d>=12.0)&&(d<13.0)) { - if (prem>=20.0) - q = 2; - else if (prem>=8.0) - q = 1; - else if (prem>=-8.0) - q = 0; - else if (prem >= -20.0) - q = -1; - else - q = -2; - return q; - } - - if ((d>=13.0)&&(d<14.0)) { - if (prem>=20.0) - q = 2; - else if (prem>=8.0) - q = 1; - else if (prem>=-8.0) - q = 0; - else if (prem >= -20.0) - q = -1; - else - q = -2; - return q; - } - - if ((d>=14.0)&&(d<15.0)) { - if (prem>=20.0) - q = 2; - else if (prem>=8.0) - q = 1; - else if (prem>=-8.0) - q = 0; - else if (prem >= -20.0) - q = -1; - else - q = -2; - return q; - } - - if ((d>=15.0)&&(d<16.0)) { - if (prem>=24.0) - q = 2; - else if (prem>=8.0) - q = 1; - else if (prem>=-8.0) - q = 0; - else if (prem >= -24.0) - q = -1; - else - q = -2; - return q; - } - } + /* This routine performs a radix-4 SRT division algorithm. The user inputs the numerator, the denominator, @@ -246,6 +139,8 @@ int main(int argc, char* argv[]) { int q; int num_iter, i; int prec; + int radix = 4; + if (argc < 5) { fprintf(stderr, "Usage: %s numerator denominator num_iterations prec\n", @@ -267,27 +162,29 @@ int main(int argc, char* argv[]) { printf("\n"); Q = 0; - P = N*0.25; + P = N * pow(2.0, -log2(radix)); printf("N = %lg, D = %lg, N/D = %lg, num_iter = %d \n\n", N, D, N/D, num_iter); for (scale = 1, i = 0; i < num_iter; i++) { // Shift by r - scale = scale*0.25; - q = qslc(flr((4*P)*16,3), D*16); - //q = -q; + scale = scale * pow(2.0, -log2(radix)); + // (4*P)*8 because of footnote in Table 5.9, page 296 EL + // i.e., real value = shown value / 8 + // D*16 since we use 4 bits of D (1 bit known) + q = qslc(flr((radix * P) * 8, 3), D*16); printf("4*W[n] = "); - disp_bin(4*P,3,prec,stdout); + disp_bin(radix*P, 3, prec, stdout); printf("\n"); printf("q*D = "); - disp_bin(q*D,3,prec,stdout); + disp_bin(q*D, 3, prec, stdout); printf("\n"); printf("W[n+1] = "); - disp_bin(P ,3,prec,stdout); + disp_bin(P ,3, prec, stdout); printf("\n"); // Recurrence - P = 4*P - q*D; + P = radix * P - q * D; // OTFC - Q = Q + q*scale; + Q = Q + q * scale; printf("i = %d, q = %d, Q = %1.18lf, W = %1.18lf\n", i, q, Q, P); printf("i = %d, q = %d", i, q); printf(", Q = "); @@ -296,8 +193,9 @@ int main(int argc, char* argv[]) { disp_bin(P, 3, prec, stdout); printf("\n\n"); } + // Is shifted partial remainder negative? if (P < 0) { - Q = Q - scale; + Q = Q - pow(2.0, -prec); P = P + D; printf("\nCorrecting Negative Remainder\n"); printf("Q = %1.18lf, W = %1.18lf\n", Q, P); @@ -306,9 +204,12 @@ int main(int argc, char* argv[]) { printf(", W = "); disp_bin(P, 3, prec, stdout); printf("\n"); - } - RQ = flr(N/D, (double) prec); - RD = Q*4; + } + + // Output Results + RQ = N/D; + // Since q_{computed} = q / radix, multiply by radix + RD = Q * radix; printf("true = %1.18lf, computed = %1.18lf, \n", RQ, RD); printf("true = "); disp_bin(RQ, 3, prec, stdout);