From 535a9a04ee008e7bdb510669c7efcb87475e7703 Mon Sep 17 00:00:00 2001 From: James Stine Date: Wed, 15 Jun 2022 11:44:09 -0500 Subject: [PATCH] Add r=4 C code --- pipelined/srt/stine/Makefile | 18 ++ pipelined/srt/stine/disp.c | 60 +++++++ pipelined/srt/stine/disp.h | 18 ++ pipelined/srt/stine/srt4div | Bin 0 -> 27024 bytes pipelined/srt/stine/srt4div.c | 325 ++++++++++++++++++++++++++++++++++ 5 files changed, 421 insertions(+) create mode 100644 pipelined/srt/stine/Makefile create mode 100755 pipelined/srt/stine/disp.c create mode 100755 pipelined/srt/stine/disp.h create mode 100755 pipelined/srt/stine/srt4div create mode 100755 pipelined/srt/stine/srt4div.c diff --git a/pipelined/srt/stine/Makefile b/pipelined/srt/stine/Makefile new file mode 100644 index 00000000..a019b388 --- /dev/null +++ b/pipelined/srt/stine/Makefile @@ -0,0 +1,18 @@ + +CC = gcc +CFLAGS = -lm +LIBS = +OBJS = disp.o srt4div.o + +srt4div: $(OBJS) + $(CC) -g -O3 -o srt4div $(OBJS) $(CFLAGS) + +disp.o: disp.h disp.c + $(CC) -g -c -o disp.o disp.c $(CFLAGS) + +srt4div.o: srt4div.c + $(CC) -g -c -o srt4div.o srt4div.c $(CFLAGS) + +clean: + rm -f *.o *~ + rm -f core diff --git a/pipelined/srt/stine/disp.c b/pipelined/srt/stine/disp.c new file mode 100755 index 00000000..66652cf0 --- /dev/null +++ b/pipelined/srt/stine/disp.c @@ -0,0 +1,60 @@ +#include "disp.h" + +double rnd_zero(double x, double bits) { + if (x < 0) + return ceiling(x, bits); + else + return flr(x, bits); +} + +double rne(double x, double precision) { + double scale, x_round; + scale = pow(2.0, precision); + x_round = rint(x * scale) / scale; + return x_round; +} + +double flr(double x, double precision) { + double scale, x_round; + scale = pow(2.0, precision); + x_round = floor(x * scale) / scale; + return x_round; +} + +double ceiling(double x, double precision) { + double scale, x_round; + scale = pow(2.0, precision); + x_round = ceil(x * scale) / scale; + return x_round; +} + +void disp_bin(double x, int bits_to_left, int bits_to_right, FILE *out_file) { + + double diff; + int i; + if (fabs(x) < pow(2.0, -bits_to_right)) { + for (i = -bits_to_left + 1; i <= bits_to_right; i++) { + fprintf(out_file,"0"); + } + return; + } + if (x < 0.0) { + // fprintf(out_file, "-"); + // x = - x; + x = pow(2.0, ((double) bits_to_left)) + x; + } + for (i = -bits_to_left + 1; i <= bits_to_right; i++) { + diff = pow(2.0, -i); + if (x < diff) { + fprintf(out_file, "0"); + } + else { + fprintf(out_file, "1"); + x -= diff; + } + if (i == 0) { + fprintf(out_file, "."); + } + } +} + diff --git a/pipelined/srt/stine/disp.h b/pipelined/srt/stine/disp.h new file mode 100755 index 00000000..bf2c7887 --- /dev/null +++ b/pipelined/srt/stine/disp.h @@ -0,0 +1,18 @@ +#include +#include +#include + +#ifndef DISP +#define DISP + +double rnd_zero(double x, double bits); + +double rne(double x, double precision); + +double flr(double x, double precision); + +double ceiling(double x, double precision); + +void disp_bin(double x, int bits_to_left, int bits_to_right, FILE *out_file); + +#endif diff --git a/pipelined/srt/stine/srt4div b/pipelined/srt/stine/srt4div new file mode 100755 index 0000000000000000000000000000000000000000..ff76638fae90ef9bb09e93e26e1f5f7e79619555 GIT binary patch literal 27024 zcmeHQ3wT^roj-SGZZdgJo^8?xofcY3X_5v*X@Su6F=Z-ENlSw*6mBOolVoHv6Xub& ztkOaSiKTp0KNYpR6c<+%)WV`$9t-IQLSfY{%3}On-H22ZD-zRIiUPC0|9RZGbCam+ z=YG52Zce}4|NDRb=YP&U_uR*v?)I%;>(Df&ij!Tzh#DvtoE9l)+b$7+7PgF);5eV1 z$MQhCI8Vr11b|eOK}RlhOTH73>}m*^0L~I@rl5+DAlXe;?kx}!1zDM($xbI!(suY} zwa`;g$<6wf2syLC@5o>jwAxW_+M0!3TWc;Hka=OU{mPEmWVc<~ZI^ZmGG!;B$|r?E zf6JslJG6z(>~`JK&JNeWClv*$4l382NQM5A`r4%3*^7i9I}8anQ&82n5q6Z9KewoE z-X!artRD``JSwQhSwl3^v9zfn8mfy%V(FcAJC`l3TiR5gjMt06DK0PmXiQzXp`E&A zfOx7q9eV_={E@6m%1M9WJG*~$^lzViGkML&?{C=6SnAsI-K0ZtlMU%mA%CjN6-}om z>F|7Lw}B{g5U-LHP!`Uq zvl-0WCM?wXcpb=jf|YapG$_(6YV$y6wwPBHM3gohIe!h_v`1Pkwsq+lOO#)He38_8ra5bFeuu!1)k!S0)k z&H(C!`yM*KmJ${e*p=(oty*ndRDV(Z#mvN)o9HDb)?{KIXhNnOIBzr@T{u4-{GW!i zB_Dqd@aahZsF2pnfOOt1!C!ht1Xqwe_EX90KkylxW-j(C&{QDIX`N*x{W!OGvI0%Wg?|V44;7}a3#$JQ$-%Nk&ijhLUGm*+WW~(+ThW783qH*W z*7;({!dG(>X%1WXg~F5{M=gAhg@4w=Rk_Q;=hts3xYxqx^PuG4XW>&FsvNTLRZdCvkcE$qH_KrQzsTe=cErLjw(tim zd<+P)9JBCCO)g`{EqpbWEaEA2Dg#p)n99IZ2BtFb-Mn$7F?_EKyS;KgHK=;__u#LXKMz`DORvYcJ(akpcVjEp&qvzY`**3b` zN;kf_uK$%A{r#`|_r5XOw%OM>(D;&ne{&tyA(`q^c#v?qv)(i3V@Luzcz<&(Z2kN7 z3IfZ|rmE5PuR)2rz92K?nX`xVo|6#zBP85^@#_RGc;DYY>OcCM75<}Voqp{F|10mO zs^H*D(t$fO)XDRz>}h=LY5obsEPY|Se{b{cg!=nWq>BCfn(qWL9C<&J84jUyUeLb@ zTDt+ZX8mM8vI9aw+F^&mS~lh>_npO|b+ps-)%hrb6G-98ID{7d9JkV0&nMo?WVo%r-#1Ds zP~Per&AJ<|M60sonCC{)<6)N}L<$Or*ubBvHkwmy>`zwvw9??-|Ke)T%(S-l?RWoF zYr}s*wM%%C`#<_de|pH@-#*CGc*=8Y$z8ni%|}Uoe4lUV?Lq41G0&m_^uzTzk=>vo zqvjsRNN7La*Kma9qYHg+^!rYq#OU@7`uksab_saHH$r`{@3gkhclt8l(5*io8Kf~F z$Ix)!Vp)&x4OBEX!DIU@Q803lD7AQp?C4a08z?fj*XF|2XVI zM+YW>It1zOf04IkcnO?9IIM)Ttinne2hQ)eUPA0gJ-6S9n$Z+}5ATu&-%#HU-iJ7p zLS2i;y-xkczl>WG2K8FjtnVC)t3xL^!X;w~b!F+T2#F;VqwqV{<`s0EL zTZ8!rQscqA3d|nNG3CS-%tFPtWya*&g1JlSj|*lHD_Ju$8%d1^b3QP8FjpxjwqU+6 zZQ{(Bd|NO@%Fno9uCfO6Ai3qid=;2Im}gX_wqQyX-aCf^o}ru>WxrqLSAW8{_xGYHHc%y(3!wqPz-d|PG~ zE50q5yA^+2Fb~3mnVEJ{{Bgm20-9!KW|G?9fYj##6Tw7DIt*NIR!(fed>YS| zC(6tnif;?%J;k>NbCTsOmwNj@`iQ?jeJJGL=bFC|PH5TIII#DlkGOU>V2yOhzpo5y z#P_kP^zSdp;0dn3|9iBMqWXsO$1<5`DLuo@z_|L}0o-Ckw|xNmE<7o>>c9P4Qwv^6 z>ln0tYV-O^R;xy69fsD^Hm!%VT3ASuug9VF8JpG}(AwvFJiD&;4PPvYQG6N?@sYTV#G+=`jFbTMqy36cR*t>BAo(#r^vG z9?FO`en=SKW*Q&MF@C}{UI9hXgmNFg`muoab16BVBT7#5oieBFd1v8%EH_%1A7s$bVB5xO-&Sd(1Pa{@`HDNp6 z;@NJ)tpKuZR=7k~NM3SN*hd4D=l|W`@cjQX-c6YKUntyZ!X2&maqrdoH2nGy{8qZT zw)tCN`g=t9OIiI55^rDiqjhR99M9>`5j0c;6NjQ1eF80{ZbbdxNB!+oKc#o;%eWH~ zohyg052CGn{ojf_irX_vsPDLcU+dt$_NVvxj;OpI?SBE;`vo42P;+TE8D5TBg%X#Y zmNJOg@NChXx2(JZ(VW-u%2$l-!JOvIpc8Pv{+!k1wKsVNzG0f=X5g3NB>1c_8Qux) zy}l#3OFg{Ycf@n|J#=#|Z&HtY?tBq}QR~peZPaQp_Ni<5Kojjhju6_UzlHE8)8WId zw0t@Og$*)Bo`#Dk#{K$cYux9ExZl2T;-)r;xCikf%p5!?f6bfvylIl#)ORtqL}%X) zMKQTh_tF_6&r+PBNBO`xE3H$Z5)ud3E{Ctth$!LzrOkwgbc;<+l3_o&`(&0PK z+W}Y0q$_@Pc$S3w^(L$Phq*!Z$4vLR`FmOn$f$JR3T-|u`S-P-4f!vu;uouwKD-%% zm!1>PH~I5j^NB8v>!~u8fvF5kWnd};QyG}bz*GjNGVuR31N8ml1f@ms=s5|>ZbgT@ z-nzOgyz`@7g)Ffni?=5OUExc;^ON3Kx+k0nq~Zy0C>)FT;QLX+!8IcIJ~fbv#A8Ws zZz3ElWbne~hu8*hvzM)*1LHW@i@ij1WJ80ADkmfgc?%0!)1nVuAG?7pZ(hV5fV{Ag zGi1OKuDjTKvxFN7Hr6*Ti*_#deke;6vRqBJ*h{Jaq{Qf6nU})V@dRQ|MPgmv4dE^{ zr!VZ?6s8Yi@y#on+&%Y(@4OS~aJJUP-e9~3pUH+pR$igfLhTkg+vID_s(A|~ts)~h zijOp?99QsJL61tvS{|tT<)&@zCstf_EWh~+?I#p}KXi)$d6kq` z4Lp8g+b16W>Wmf%rMyPUy$VYHEJ>fIpyV%*^a4rGlk`FbCBI(M^^#s9=_Qh0Bn?k z#82Um+B5*D3e|Nt{+>X69mMpNmfh)CU6SY6iikukkAVFL$amz(W2XE&;BP`MCX@df z@~`|6c`@Y#x$7aQ>UV}Ts(|8K41Oze*+Ad(nDPY$O1>WQZNJK7zK=I|S@}zOJV3Cf z%2WoXGBA~asSHeIU@8Ms8Tg;bfcpKG`n?t{b*Q{7fccxJS0!C6q544?UTN?$G+#h` zdBltQofy7o;{{7LUhvHpFY5PVw2q>p>{Ddm9_^%gZHojmy)V zZoEKp*kU*q85RBY75S6Erb*H$^gn@5xrF{tP+LxrLUrIXV=N5OjMJ5 z9)Pcm6t}T&|CE%Jpqf<+lPKbTzuqftmJsa6XW~5=W-t zT{v4>@F^VK+Fdwa!tpY#4~ok;UZXt^d^yLxS`i#w%JGF-82DuzU!q+O+siq=OnU%$ zGsj!BCg`+q+^;M{q*a>zfXGxDKaI{!TV$8dc#qu*DSyoHKhRr#8EWh-wg2 zo9k(KnX`t%b{z2Q+aa43p!Pa>8q^HZ0TvbQAe=gDWW!va2Fp25P}><( zGtoihv7FCDuaXgpeGaRKM&3LKGYg>4>J0zs7X5 z7Jk%pM5p>h(-EEcYC57RMNLN(!+&KuibMZ-mM5npnq8*J=_sG3BNt7#E@Br^H&Haj zVmfj)0M93hi>AhWnvPsFF>=M-O0kl;sP{^UQTjd(_j0Td=hp$kJaiHkKL$770_7lt zAJ}x7X7bDOerU2+l8hh7ka_rdMh2s%3Q?@)7ohY0Nh|&dAnBK*GKYC^GuQoKmJp*> z$HCp}D5)nE(q(18gVl{1;WVE)=nCYRMry3=ZQvm;GDt)o+-$~t0PU1c7R7VO1!Mg2 z9-wyc3Y_r;Gof*+loXr?4(=`>yHhhr;yiHgZoQ42DdO!u(4uc(XG#fo9yq*PU&GFn z3*yLbeIq+lCA9qdCJ@sFu|?kmwV8t0K9d5~&dedQkmkNK^M#~MznT5cEh2bicj238 z$lH=Uv|Ep`w=0C~q1}Z)fb6W4J-C~>J_R?kyyVx#+o&cD{)iJY?O*}zkmRb60nH9U zlSkko8}!bgrfe zHAfAuludR{p3X=$AFA$4=UPQ*dq719+qBSi*#t(x=IZuMRQ@6mR)UXBX8%|T4R}p#JE*!DO1yP zb#o)HnM_FCIKFR0ik;ama+zI39g>^(ylne;$Edu!5fm#J*Q15wdbCJ&P_gQul1q3W zmkL$hOQq{7SJYB>dSt5R=2PeKJ`$Z*CfkT!)tApL&StTEe799(yRA}mTh#(yWwpKc zre#|Rf{4^kDn0(%=xji+Lx7i%UvnJEUZ3RCAq$Z63 zwn(Ep4>zORSxwXaP+CF&OuBOiG6xfnctgY2!@wQd9ET^rj@bP06(AZvou35eJW7XS zI5_l$tiYjGA60sud#hgjHTTQ9=Kh+#uum^?U#FKi&aDJ;OfNj4YxhsWT|A-ghCV&t z{Y8DIdy8JNMbC50tI&&|()CC5eB7rv8}y|FGWxjy+6Q(04n6NWz2F*`i!bm@ORw z{ys}(>!*|n%2)wkNAbf%89NC**D138Jz@Gh}IUh}`wM`-EQT zZqQ2(=mjS%M|Z;~_2RUicSnw*RYT%AYHx6NQ(7+VPiCTw^8bHXOS@j5=Panl5X=8n zmt%gl?m48pZr45Tm&OA(&Srk*xvvt?b(Oo>{loj|pywUZ3l8XcTXYvrERUnPN}qO2 zcct`VWYBq?UPu5_LN$k5Bwe9nvnrRo=~N#K49=OQ}?^%ar)kY3K!DsVGL z{=NJ>aDGyduLG6u(T%d4$5lNigZ|JjuU=GQ5&Ulhel;RUQCY`(= zp1{;WE#<8(IHc!0F0Mplw_4aI+)r_7(JyrEI>$838HN$>>@-pgn}X4XW$Ugs65&9| z=uIUUwg#hZ#*)~53<^fFJDx~&q&rz67BX%LC*oWi45K5E3^VL{*5&1gbSxR^is6$A zdIA)Rbat|iNGfTh;zl&wnX)L792=ZP@tI1t=5Q>;^&yD{x{_uk#5aQJgb|EK(>*cH z55_~`U>_=pr}=hk@o0#elMJWWP9qUdOIK8NBGT23+R(n9K;kBRDPnXc!q~MfX;72+ zrfVeJk-$eXy!FJwM<$Fn+nduSQcMQ=!lL0~?=~YE58fnWzK9t;hzxB{_8J{XVb(Km z1e~x2`jg+!YIUC*rYqI$4*F?ZCG* zb=+eTLcUF0U2h_u!e=+}L|rlkMY0ltU^I~IR$_cx1VOOVLNof*D3XK(Th5_#mNv1*2&U z6|N|W@rrc*{23|@sa?HczWZM&+*#kvI0=XzO9xYYpn>WMq`GmkN0TH;q?$sJJ{75o z2C0a|A*KP$?L%k-kUvsrh2~ns011)Ru zf*`)P!pg-nUddPQdz5q?5Doh^`r8DD;-k`zKaEXe>U|Ao;_wmtk^e(B|95fyYIC1J zrT=9h6UG1aN%&YZPpE$cd<|a>9cuz&D--#D0XiDq;@ID7{42MuWNJTG zQv9<(Y)$1N`bKCX{srKBm9Yi-z@KOz&6l}8?sa7TmB|+HC-NWR{4%E2=N-}@0e%^) zWo0%$yC(5-m-J(QDtb5Zu}iCZuczYqZVNz_Z|&J8GGFPozXTo?etwMwbW9CHZPCuO5WV<$UEe3+aje zME$iM{4$(hZMJ^-9&cBfxqtNn>5AS3kiNYiZ{YsR#yxzZC*c&o{V_EC`g5Xqc22_Y z2cPn0ulLj7d$Y|nv2RS`X8`8I0_mNMK1_{<@#Sj$l8H5F z>)419NF)NgjBqTK*u^>%fgUU!(mg%9U}B*RXr^*h`OF`*TptW$?WUEjKEt{(2zRkwwm8;hK47C?+i@A#|c9b2r zrL8H&CbEBQUtO`+t*K4B)6%}TIgIQoHOv=Fcqf#VFUNw5mVjtCqR8C{*VJmtF5vi9 zyQYBZL#JQ~Cbr{M%L$u8?%uyfC>b}p1F?`?X45q`X3!dkuRdw1$=4dRo3G7L&c?l_ z>=K4<#hTdYNG4?*e4Rp1mie00)TDT6iIc-H3~W`bx|8-Wo^WGgS_Rw8th*BDif9zc zy@zuItUkG`2kS||R6@XRg`u_KL@%q4VR>1PMNvH#n!VvfYL|u8k&a*=@JL9qR<2ry z)oB;w>fKmox_|_bpP) zbrnwVT>R0eB+7o7v{z7Grbv|+Jao<}OF>!f?fc=Mm!nNn zA;^D+Z&FnJ>bo5URs8(*Aaw2V2f)JN^_+7fMYPztU6iR>TiUDsSL;#*&zb_^z+V4%ZT9MZM!}G*R{?wd zKepL#ryE2R1*6nBOHt!Tu}%QN`#gRQ5>I{ir~dy6yS}0o{4H!e7JGG{xuu!BWyh_u} z-_jyXR6WXpvQsqOKN3&X?`@NIEdnkPz&v_QOd#d5e!E^-i;yp{=~Jv!6r3UL1vh(K z+9L4nSr*5Xonq5{uAN7+&$I}QDte${DR$k9C$ax3U6@chRbm7=uCCRS*q5~mgI;+( zkmlp;ua@@q{O#E!EEY{tziP*Ilh}Wi9-N>!m7*Zc`W2g??G)#lhp!a|Ei(RdLDDl3 oHGb6mPwBOV$PGuY=YXBIIal=xDcy5TZ2!iMLgT|W0~^Wy6DEJmH~;_u literal 0 HcmV?d00001 diff --git a/pipelined/srt/stine/srt4div.c b/pipelined/srt/stine/srt4div.c new file mode 100755 index 00000000..cbf0d822 --- /dev/null +++ b/pipelined/srt/stine/srt4div.c @@ -0,0 +1,325 @@ +#include "disp.h" +#include + +int qslc (double prem, double d) { + + int q; + + 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) + q = 2; + else if (prem>=2.0) + q = 1; + else if (prem>=-2.0) + q = 0; + else if (prem >= -5) + q = -1; + else + q = -2; + return q; + } + + if ((d>=4.0)&&(d<5.0)) { + if (prem>=6.66667) + q = 2; + else if (prem>=2.0) + q = 1; + else if (prem>=-2.0) + q = 0; + else if (prem >= -6.66667) + q = -1; + else + q = -2; + return q; + } + + if ((d>=5.0)&&(d<6.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>=6.0)&&(d<7.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>=7.0)&&(d<8.0)) { + if (prem>=11.0) + q = 2; + else if (prem>=4.0) + q = 1; + else if (prem>=-4.0) + q = 0; + else if (prem >= -11.0) + q = -1; + else + q = -2; + return q; + } + + if ((d>=8.0)&&(d<9.0)) { + if (prem>=12.0) + q = 2; + else if (prem>=4.0) + q = 1; + else if (prem>=-4.0) + q = 0; + else if (prem >= -12.0) + q = -1; + else + q = -2; + 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, + and the number of iterations. It assumes that 0.5 <= D < 1. + +*/ + +int main(int argc, char* argv[]) { + + double P, N, D, Q, RQ, RD, RREM, scale; + int q; + int num_iter, i; + int prec; + if (argc < 5) { + fprintf(stderr, + "Usage: %s numerator denominator num_iterations prec\n", + argv[0]); + exit(1); + } + sscanf(argv[1],"%lg", &N); + sscanf(argv[2],"%lg", &D); + sscanf(argv[3],"%d", &num_iter); + sscanf(argv[4],"%d", &prec); + // Round to precision + N = rne(N, prec); + D = rne(D, prec); + printf("N = "); + disp_bin(N, 3, prec, stdout); + printf("\n"); + printf("D = "); + disp_bin(D, 3, prec, stdout); + printf("\n"); + + Q = 0; + P = N*0.25; + 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; + printf("4*W[n] = "); + disp_bin(4*P,3,prec,stdout); + printf("\n"); + printf("q*D = "); + disp_bin(q*D,3,prec,stdout); + printf("\n"); + printf("W[n+1] = "); + disp_bin(P ,3,prec,stdout); + printf("\n"); + // Recurrence + P = 4*P - q*D; + // OTFC + 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 = "); + disp_bin(Q, 3, prec, stdout); + printf(", W = "); + disp_bin(P, 3, prec, stdout); + printf("\n\n"); + } + if (P < 0) { + Q = Q - scale; + P = P + D; + printf("\nCorrecting Negative Remainder\n"); + printf("Q = %1.18lf, W = %1.18lf\n", Q, P); + printf("Q = "); + disp_bin(Q, 3, prec, stdout); + printf(", W = "); + disp_bin(P, 3, prec, stdout); + printf("\n"); + } + RQ = flr(N/D, (double) prec); + RD = Q*4; + printf("true = %1.18lf, computed = %1.18lf, \n", RQ, RD); + printf("true = "); + disp_bin(RQ, 3, prec, stdout); + printf(", computed = "); + disp_bin(RD, 3, prec, stdout); + printf("\n\n"); + printf("REM = %1.18lf \n", P); + printf("REM = "); + disp_bin(P, 3, prec, stdout); + printf("\n\n"); + + return 0; + +}