home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / PLYDIV.C < prev    next >
Encoding:
Text File  |  1984-06-23  |  1.9 KB  |  84 lines
  1.    plydiv(p,s,ndp,nds,q,r,ndq,ndr)
  2.  
  3.       /*this function divides two polynomials and retains the   */
  4.       /*remainder. Divisor, dividend, quotient, and remainder   */
  5.       /*are of the form a(1)*x**nda + a(2)*x**(nda-1) + ...     */
  6.       /*+ a(nda)*x + a(nda+1).                                  */
  7.  
  8.       int ndp,nds,*ndq,*ndr;
  9.       float p[],s[],q[],r[];
  10.  
  11.     {
  12.       int nv,ntq,i,ip,j,js,k,kq,it,n,nq;
  13.       float xx,zz;
  14.       double fabs();
  15.  
  16.       nv = 0;
  17.  /* if q[0]=999999. the function calculates only the remainder */
  18.       if (q[0] == 999999.) nv = 1;
  19.       *ndq = ndp - nds;
  20.       ntq = *ndq ;
  21.       j = ndp;
  22.  
  23.       for(ip = 0; ip<= j; ip++)
  24.         r[ip] = p[ip];
  25.       if(*ndq <  0) goto r1;
  26.       ip = 0;
  27.       js = 0;
  28.       kq = 0;
  29.       it = 0;
  30.       q[kq] = p[ip]/s[js];
  31.  
  32.       r[ip] = p[ip] - s[js]*q[kq];
  33.       while(js != (nds))
  34.        {
  35.         js = js + 1;
  36.         ip = ip + 1;
  37.         r[ip] = p[ip] - s[js]*q[kq];
  38.        }
  39.  
  40. a5:   if(ntq <= kq) goto r2;
  41.       kq = kq + 1;
  42.       it = it + 1;
  43.       ip = it;
  44.       js = 0;
  45.       if (r[ip] == 0.) goto a11;
  46.       if (s[0] == 0.) goto a9;
  47.       xx = r[ip]/s[0];
  48. a8:   if(fabs(xx) <= 0.000005) goto a11;
  49. a9:   nq = 0;
  50.       if(nv != 1) nq = kq;
  51.       q[nq] = r[ip]/s[js];
  52. a13:  r[ip] = r[ip] - s[js]*q[nq];
  53. a14:  if(js == (nds)) goto a5;
  54.       js = js + 1;
  55.       ip = ip + 1;
  56.       goto a13;
  57. a11:  if(nv != 1) q[kq] = 0.0;
  58.       goto a5;
  59. r1:   if(nv != 1) q[0] = 0.0;
  60.       *ndq = 0;
  61.       *ndr = ndp;
  62.       return;
  63. r2:   j = ndp;
  64.       i = 0;
  65. a17:  if (r[i] == 0.) goto a19;
  66.       if (s[i] == 0.) goto a16;
  67.       zz = r[i]/s[i];
  68. a18:  if (fabs(zz) <= 0.000005) goto a19;
  69. a16:  *ndr = j - i;
  70.       n = *ndr;
  71.       for(k = 0; k <= n; k++)
  72.       {
  73.        r[k] = r[i];
  74.        i = i + 1;
  75.       }
  76.       return;
  77. a19:  if(j <= i) goto r22;
  78.       i = i + 1;
  79.       goto a17;
  80. r22:  *ndr = 0;
  81.       r[0] = 0.0;
  82.       return;
  83.    }
  84.