quotient.cpp

   1 // Divide polynomials
   2
   3 #include "stdafx.h"
   4 #include "defs.h"
   5
   6 void
   7 eval_quotient(void)
   8 {
   9         push(cadr(p1));                 // 1st arg, p(x)
  10         eval();
  11
  12         push(caddr(p1));                // 2nd arg, q(x)
  13         eval();
  14
  15         push(cadddr(p1));               // 3rd arg, x
  16         eval();
  17
  18         p1 = pop();                     // default x
  19         if (p1 == symbol(NIL))
  20                 p1 = symbol(SYMBOL_X);
  21         push(p1);
  22
  23         divpoly();
  24 }
  25
  26 //-----------------------------------------------------------------------------
  27 //
  28 //      Divide polynomials
  29 //
  30 //      Input:          tos-3           Dividend
  31 //
  32 //                      tos-2           Divisor
  33 //
  34 //                      tos-1           x
  35 //
  36 //      Output:         tos-1           Quotient
  37 //
  38 //-----------------------------------------------------------------------------
  39
  40 #define DIVIDEND p1
  41 #define DIVISOR p2
  42 #define X p3
  43 #define Q p4
  44 #define QUOTIENT p5
  45
  46 void
  47 divpoly(void)
  48 {
  49         int h, i, m, n, x;
  50         U **dividend, **divisor;
  51
  52         save();
  53
  54         X = pop();
  55         DIVISOR = pop();
  56         DIVIDEND = pop();
  57
  58         h = tos;
  59
  60         dividend = stack + tos;
  61
  62         push(DIVIDEND);
  63         push(X);
  64         m = coeff() - 1;        // m is dividend's power
  65
  66         divisor = stack + tos;
  67
  68         push(DIVISOR);
  69         push(X);
  70         n = coeff() - 1;        // n is divisor's power
  71
  72         x = m - n;
  73
  74         push_integer(0);
  75         QUOTIENT = pop();
  76
  77         while (x >= 0) {
  78
  79                 push(dividend[m]);
  80                 push(divisor[n]);
  81                 divide();
  82                 Q = pop();
  83
  84                 for (i = 0; i <= n; i++) {
  85                         push(dividend[x + i]);
  86                         push(divisor[i]);
  87                         push(Q);
  88                         multiply();
  89                         subtract();
  90                         dividend[x + i] = pop();
  91                 }
  92
  93                 push(QUOTIENT);
  94                 push(Q);
  95                 push(X);
  96                 push_integer(x);
  97                 power();
  98                 multiply();
  99                 add();
 100                 QUOTIENT = pop();
 101
 102                 m--;
 103                 x--;
 104         }
 105
 106         tos = h;
 107
 108         push(QUOTIENT);
 109
 110         restore();
 111 }
 112
 113 #if SELFTEST
 114
 115 static char *s[] = {
 116
 117         "quotient(x^2+1,x+1)-x+1",
 118         "0",
 119
 120         "quotient(a*x^2+b*x+c,d*x+e)-(-a*e/(d^2)+a*x/d+b/d)",
 121         "0",
 122 };
 123
 124 void
 125 test_quotient(void)
 126 {
 127         test(__FILE__, s, sizeof s / sizeof (char *));
 128 }
 129
 130 #endif