1 files changed, 110 insertions, 0 deletions
diff --git a/2048/CH4/EX4.5/respol.sci b/2048/CH4/EX4.5/respol.sci
new file mode 100755
index 000000000..ea6818967
--- /dev/null
+++ b/2048/CH4/EX4.5/respol.sci
@@ -0,0 +1,110 @@
+// Computation of residues 
+// 4.5
+// Numerator and denominator coefficients
+// are passed in decreasing powers of z(say)
+
+function [res,pol,q] = respol(num,den)
+len = length(num);
+if num(len) == 0
+    num = num(1:len-1);
+end
+
+[resi,q] = pfe(num,den);
+res = resi(:,2);
+res = int(res) + (clean(res - int(res),1.d-04));
+pol = resi(:,1);
+pol = int(pol) + (clean(pol - int(pol),1.d-04));
+endfunction;
+
+///////////////////////////////////////////////////
+// Partial fraction expansion
+
+function [resid1,q] = pfe(num,den)
+x = poly(0,'x');
+s = %s;
+
+num2 = flip(num);
+den2 = flip(den);
+num = poly(num2,'s','coeff');
+den = poly(den2,'s','coeff');
+[fac,g] = factors(den);
+polf = polfact(den);
+n = 1;  
+
+r = clean(real(roots(den)),1.d-5);
+//disp(r);
+l = length(r);
+r = gsort(r,'g','i');
+r = [r; 0];
+j = 1;
+t1 = 1; q = [];
+rr = 0;
+rr1 = [0 0];
+m = 1;
+
+  for i = j:l
+     if abs(r(i)- r(i+1)) < 0.01 then
+       r(i);
+       r(i+1);
+       n = n+1;
+       m = n;
+       //pause
+       t1 = i;
+       //disp('Repeated roots')
+     else
+       m = n;
+       //pause
+       n = 1;
+    end
+    i;
+    if n == 1 then
+      rr1 = [rr1; r(i) m];
+      //pause
+    end;
+       j = t1 + 1;
+  end;
+rr2 = rr1(2:$,:);
+[r1,c1] = size(rr2);
+den1 = 1;
+
+for i = 1:r1
+  den1 = den1 * ((s-rr2(i,1))^(rr2(i,2)));
+end;
+[rem,quo] = pdiv(num,den);
+q = quo;
+if quo ~= 0
+  num = rem;
+end
+
+tf = num/den1;
+res1 = 0;
+res3 = [s 0];
+res5 = [0 0];
+for i = 1:r1
+  j = rr2(i,2) + 1;
+  tf1 = tf; //strictly proper
+  k = rr2(i,2);
+  tf2 = ((s-rr2(i,1))^k)*tf1;
+  rr2(i,1);
+  res1 = horner(tf2,rr2(i,1));
+  res2 = [(s - rr2(i,1))^(rr2(i,2)) res1];
+  res4 = [rr2(i,1) res1];
+  res3 = [res3; res2];
+  res5 = [res5; res4];
+  res = res1;
+    for m = 2:j-1
+    k;
+    rr2(i,1);
+      tf1 = derivat(tf2)/factorial(m-1); //ith derivative
+      res = horner(tf1,rr2(i,1));
+      res2 = [(s - rr2(i,1))^(j-m) res];
+      res4 = [rr2(i,1) res];
+      res5 = [res5; res4];
+      res3 = [res3; res2];
+      tf2 = tf1;
+    end;
+end;      
+resid = res3(2:$,:); //with s terms
+resid1 = res5(2:$,:); //only poles(in decreasing no. of repetitions)
+endfunction;
+////////////////////////////////////////////////////////////