Polynomial long division: Difference between revisions
Undo revision 361341 by Peak (talk)
m (→{{header|Octave}}: specify the builtin to make the poly long div) |
|||
| (122 intermediate revisions by 64 users not shown) | |||
Line 2:
:<cite>In algebra, [[wp:Polynomial long division|polynomial long division]] is an algorithm for dividing a polynomial by another polynomial of the same or lower degree.</cite>
Let us suppose a polynomial is represented by a vector, <math>x</math> (i.e., an ordered collection of [[wp:Coefficient|coefficients]]) so that the <math>i</math><sup>th</sup> element keeps the coefficient of <math>x^i</math>, and the multiplication by a monomial is a ''shift'' of the vector's elements "towards right" (injecting
Then a pseudocode for the polynomial long division using the conventions described above could be:
Line 13:
<span class="co1">// '''N''', '''D''', '''q''', '''r''' are vectors</span>
'''if''' degree('''D''') < 0 '''then''' ''error''
'''
'''while'''
'''
<span class="co1">//
'''d''' ←
'''r''' ← '''
'''return''' ('''q''', '''r''')
Line 32 ⟶ 27:
* Error handling (for allocations or for wrong inputs) is not mandatory.
* Conventions can be different; in particular, note that if the first coefficient in the vector is the highest power of x for the polynomial represented by the vector, then the algorithm becomes simpler
'''Example for clarification'''
Line 92 ⟶ 88:
q: -27 -9 1 → x<sup>2</sup> - 9x - 27
r: -123 0 0 → -123
;Related task:
:* [[Polynomial derivative]]
=={{header|11l}}==
{{trans|Python}}
<syntaxhighlight lang="11l">F degree(&poly)
L !poly.empty & poly.last == 0
poly.pop()
R poly.len - 1
F poly_div(&n, &D)
V dD = degree(&D)
V dN = degree(&n)
I dD < 0
exit(1)
[Float] q
I dN >= dD
q = [0.0] * dN
L dN >= dD
V d = [0.0] * (dN - dD) [+] D
V mult = n.last / Float(d.last)
q[dN - dD] = mult
d = d.map(coeff -> coeff * @mult)
n = zip(n, d).map((coeffN, coeffd) -> coeffN - coeffd)
dN = degree(&n)
E
q = [0.0]
R (q, n)
print(‘POLYNOMIAL LONG DIVISION’)
V n = [-42.0, 0.0, -12.0, 1.0]
V D = [-3.0, 1.0, 0.0, 0.0]
print(‘ #. / #. =’.format(n, D), end' ‘ ’)
V (q, r) = poly_div(&n, &D)
print(‘ #. remainder #.’.format(q, r))</syntaxhighlight>
{{out}}
<pre>
POLYNOMIAL LONG DIVISION
[-42, 0, -12, 1] / [-3, 1, 0, 0] = [-27, -9, 1] remainder [-123]
</pre>
=={{header|Ada}}==
long_division.adb:
<syntaxhighlight lang="ada">with Ada.Text_IO; use Ada.Text_IO;
procedure Long_Division is
package Int_IO is new Ada.Text_IO.Integer_IO (Integer);
use Int_IO;
type Degrees is range -1 .. Integer'Last;
subtype Valid_Degrees is Degrees range 0 .. Degrees'Last;
type Polynom is array (Valid_Degrees range <>) of Integer;
function Degree (P : Polynom) return Degrees is
begin
for I in reverse P'Range loop
if P (I) /= 0 then
return I;
end if;
end loop;
return -1;
end Degree;
function Shift_Right (P : Polynom; D : Valid_Degrees) return Polynom is
Result : Polynom (0 .. P'Last + D) := (others => 0);
begin
Result (Result'Last - P'Length + 1 .. Result'Last) := P;
return Result;
end Shift_Right;
function "*" (Left : Polynom; Right : Integer) return Polynom is
Result : Polynom (Left'Range);
begin
for I in Result'Range loop
Result (I) := Left (I) * Right;
end loop;
return Result;
end "*";
function "-" (Left, Right : Polynom) return Polynom is
Result : Polynom (Left'Range);
begin
for I in Result'Range loop
if I in Right'Range then
Result (I) := Left (I) - Right (I);
else
Result (I) := Left (I);
end if;
end loop;
return Result;
end "-";
procedure Poly_Long_Division (Num, Denom : Polynom; Q, R : out Polynom) is
N : Polynom := Num;
D : Polynom := Denom;
begin
if Degree (D) < 0 then
raise Constraint_Error;
end if;
Q := (others => 0);
while Degree (N) >= Degree (D) loop
declare
T : Polynom := Shift_Right (D, Degree (N) - Degree (D));
begin
Q (Degree (N) - Degree (D)) := N (Degree (N)) / T (Degree (T));
T := T * Q (Degree (N) - Degree (D));
N := N - T;
end;
end loop;
R := N;
end Poly_Long_Division;
procedure Output (P : Polynom) is
First : Boolean := True;
begin
for I in reverse P'Range loop
if P (I) /= 0 then
if First then
First := False;
else
Put (" + ");
end if;
if I > 0 then
if P (I) /= 1 then
Put (P (I), 0);
Put ("*");
end if;
Put ("x");
if I > 1 then
Put ("^");
Put (Integer (I), 0);
end if;
elsif P (I) /= 0 then
Put (P (I), 0);
end if;
end if;
end loop;
New_Line;
end Output;
Test_N : constant Polynom := (0 => -42, 1 => 0, 2 => -12, 3 => 1);
Test_D : constant Polynom := (0 => -3, 1 => 1);
Test_Q : Polynom (Test_N'Range);
Test_R : Polynom (Test_N'Range);
begin
Poly_Long_Division (Test_N, Test_D, Test_Q, Test_R);
Put_Line ("Dividing Polynoms:");
Put ("N: "); Output (Test_N);
Put ("D: "); Output (Test_D);
Put_Line ("-------------------------");
Put ("Q: "); Output (Test_Q);
Put ("R: "); Output (Test_R);
end Long_Division;</syntaxhighlight>
output:
<pre>Dividing Polynoms:
N: x^3 + -12*x^2 + -42
D: x + -3
-------------------------
Q: x^2 + -9*x + -27
R: -123</pre>
=={{header|ALGOL 68}}==
<syntaxhighlight lang="algol68">
BEGIN # polynomial division #
# in this polynomials are represented by []INT items where #
# the coefficients are in order of increasing powers, i.e., #
# element 0 = coefficient of x^0, element 1 = coefficient of #
# x^1, etc. #
# returns the degree of the polynomial p, the highest index of #
# p where the element is non-zero or - max int if all #
# elements of p are 0 #
OP DEGREE = ( []INT p )INT:
BEGIN
INT result := - max int;
FOR i FROM LWB p TO UPB p DO
IF p[ i ] /= 0 THEN result := i FI
OD;
result
END # DEGREE # ;
MODE POLYNOMIALDIVISIONRESULT = STRUCT( FLEX[ 1 : 0 ]INT q, r );
# in-place multiplication of the elements of a by b returns a #
OP *:= = ( REF[]INT a, INT b )REF[]INT:
BEGIN
FOR i FROM LWB a TO UPB a DO
a[ i ] *:= b
OD;
a
END # *:= # ;
# subtracts the corresponding elements of b from those of a, #
# a and b must have the same bounds - returns a #
OP -:= = ( REF[]INT a, []INT b )REF[]INT:
BEGIN
FOR i FROM LWB a TO UPB a DO
a[ i ] -:= b[ i ]
OD;
a
END # -:= # ;
# returns the polynomial a right-shifted by shift, the bounds #
# are unchanged, so high order elements are lost #
OP SHR = ( []INT a, INT shift )[]INT:
BEGIN
INT da = DEGREE a;
[ LWB a : UPB a ]INT result;
FOR i FROM LWB result TO shift - ( LWB result + 1 ) DO result[ i ] := 0 OD;
FOR i FROM shift - LWB result TO UPB result DO result[ i ] := a[ i - shift ] OD;
result
END # SHR # ;
# polynomial long disivion of n in by d in, returns q and r #
OP / = ( []INT n in, d in )POLYNOMIALDIVISIONRESULT:
IF DEGREE d < 0 THEN
print( ( "polynomial division by polynomial with negative degree", newline ) );
stop
ELSE
[ LWB d in : UPB d in ]INT d := d in;
[ LWB n in : UPB n in ]INT n := n in;
[ LWB n in : UPB n in ]INT q; FOR i FROM LWB q TO UPB q DO q[ i ] := 0 OD;
INT dd in = DEGREE d in;
WHILE DEGREE n >= dd in DO
d := d in SHR ( DEGREE n - dd in );
q[ DEGREE n - dd in ] := n[ DEGREE n ] OVER d[ DEGREE d ];
# DEGREE d is now DEGREE n #
d *:= q[ DEGREE n - dd in ];
n -:= d
OD;
( q, n )
FI # / # ;
# displays the polynomial p #
OP SHOWPOLYNOMIAL = ( []INT p )VOID:
BEGIN
BOOL first := TRUE;
FOR i FROM UPB p BY - 1 TO LWB p DO
IF INT e = p[ i ];
e /= 0
THEN
print( ( IF e < 0 AND first THEN "-"
ELIF e < 0 THEN " - "
ELIF first THEN ""
ELSE " + "
FI
, IF ABS e = 1 THEN "" ELSE whole( ABS e, 0 ) FI
)
);
IF i > 0 THEN
print( ( "x" ) );
IF i > 1 THEN print( ( "^", whole( i, 0 ) ) ) FI
FI;
first := FALSE
FI
OD;
IF first THEN
# degree is negative #
print( ( "(negative degree)" ) )
FI
END # SHOWPOLYNOMIAL # ;
[]INT n = ( []INT( -42, 0, -12, 1 ) )[ AT 0 ];
[]INT d = ( []INT( -3, 1, 0, 0 ) )[ AT 0 ];
POLYNOMIALDIVISIONRESULT qr = n / d;
SHOWPOLYNOMIAL n; print( ( " divided by " ) ); SHOWPOLYNOMIAL d;
print( ( newline, " -> Q: " ) ); SHOWPOLYNOMIAL q OF qr;
print( ( newline, " R: " ) ); SHOWPOLYNOMIAL r OF qr
END
</syntaxhighlight>
{{out}}
<pre>
x^3 - 12x^2 - 42 divided by x - 3
-> Q: x^2 - 9x - 27
R: -123
</pre>
=={{header|APL}}==
<syntaxhighlight lang="apl">div←{
{
q r d←⍵
(≢d) > n←≢r : q r
c ← (⊃⌽r) ÷ ⊃⌽d
∇ (c,q) ((¯1↓r) - c × ¯1↓(-n)↑d) d
} ⍬ ⍺ ⍵
}
</syntaxhighlight>
{{out}}
<pre> N←¯42 0 ¯12 1
D←¯3 1
⍪N div D
¯27 ¯9 1
¯123
</pre>
=={{header|BBC BASIC}}==
<syntaxhighlight lang="bbcbasic"> DIM N%(3) : N%() = -42, 0, -12, 1
DIM D%(3) : D%() = -3, 1, 0, 0
DIM q%(3), r%(3)
PROC_poly_long_div(N%(), D%(), q%(), r%())
PRINT "Quotient = "; FNcoeff(q%(2)) "x^2" FNcoeff(q%(1)) "x" FNcoeff(q%(0))
PRINT "Remainder = " ; r%(0)
END
DEF PROC_poly_long_div(N%(), D%(), q%(), r%())
LOCAL d%(), i%, s%
DIM d%(DIM(N%(),1))
s% = FNdegree(N%()) - FNdegree(D%())
IF s% >= 0 THEN
q%() = 0
WHILE s% >= 0
FOR i% = 0 TO DIM(d%(),1) - s%
d%(i%+s%) = D%(i%)
NEXT
q%(s%) = N%(FNdegree(N%())) DIV d%(FNdegree(d%()))
d%() = d%() * q%(s%)
N%() -= d%()
s% = FNdegree(N%()) - FNdegree(D%())
ENDWHILE
r%() = N%()
ELSE
q%() = 0
r%() = N%()
ENDIF
ENDPROC
DEF FNdegree(a%())
LOCAL i%
i% = DIM(a%(),1)
WHILE a%(i%)=0
i% -= 1
IF i%<0 EXIT WHILE
ENDWHILE
= i%
DEF FNcoeff(n%)
IF n%=0 THEN = ""
IF n%<0 THEN = " - " + STR$(-n%)
IF n%=1 THEN = " + "
= " + " + STR$(n%)</syntaxhighlight>
'''Output:'''
<pre>
Quotient = + x^2 - 9x - 27
Remainder = -123
</pre>
=={{header|C}}==
{{trans|Fortran}}
{{libheader|
<
#include <stdlib.h>
#include <stdarg.h>
Line 203 ⟶ 551:
return r;
}</
<
{
int i;
Line 227 ⟶ 575:
return 0;
}</
===Another version===
Without outside libs, for clarity. Note that polys are stored and show with zero-degree term first:<syntaxhighlight lang="c">#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
typedef struct {
int power;
double * coef;
} poly_t, *poly;
#define E(x, i) (x)->coef[i]
/* passing in negative power to have a zeroed poly */
poly p_new(int power, ...)
{
int i, zeroed = 0;
va_list ap;
if (power < 0) {
power = -power;
zeroed = 1;
}
poly p = malloc(sizeof(poly_t));
p->power = power;
p->coef = malloc(sizeof(double) * ++power);
if (zeroed)
for (i = 0; i < power; i++) p->coef[i] = 0;
else {
va_start(ap, power);
for (i = 0; i < power; i++)
E(p, i) = va_arg(ap, double);
va_end(ap);
}
return p;
}
void p_del(poly p)
{
free(p->coef);
free(p);
}
void p_print(poly p)
{
int i;
for (i = 0; i <= p->power; i++)
printf("%g ", E(p, i));
printf("\n");
}
poly p_copy(poly p)
{
poly q = p_new(-p->power);
memcpy(q->coef, p->coef, sizeof(double) * (1 + p->power));
return q;
}
/* p: poly; d: divisor; r: remainder; returns quotient */
poly p_div(poly p, poly d, poly* r)
{
poly q;
int i, j;
int power = p->power - d->power;
double ratio;
if (power < 0) return 0;
q = p_new(-power);
*r= p_copy(p);
for (i = p->power; i >= d->power; i--) {
E(q, i - d->power) = ratio = E(*r, i) / E(d, d->power);
E(*r ,i) = 0;
for (j = 0; j < d->power; j++)
E(*r, i - d->power + j) -= E(d, j) * ratio;
}
while (! E(*r, --(*r)->power));
return q;
}
int main()
{
poly p = p_new(3, 1., 2., 3., 4.);
poly d = p_new(2, 1., 2., 1.);
poly r;
poly q = p_div(p, d, &r);
printf("poly: "); p_print(p);
printf("div: "); p_print(d);
printf("quot: "); p_print(q);
printf("rem: "); p_print(r);
p_del(p);
p_del(q);
p_del(r);
p_del(d);
return 0;
}</syntaxhighlight>
=={{header|C sharp|C#}}==
{{trans|Java}}
<syntaxhighlight lang="csharp">using System;
namespace PolynomialLongDivision {
class Solution {
public Solution(double[] q, double[] r) {
Quotient = q;
Remainder = r;
}
public double[] Quotient { get; }
public double[] Remainder { get; }
}
class Program {
static int PolyDegree(double[] p) {
for (int i = p.Length - 1; i >= 0; --i) {
if (p[i] != 0.0) return i;
}
return int.MinValue;
}
static double[] PolyShiftRight(double[] p, int places) {
if (places <= 0) return p;
int pd = PolyDegree(p);
if (pd + places >= p.Length) {
throw new ArgumentOutOfRangeException("The number of places to be shifted is too large");
}
double[] d = new double[p.Length];
p.CopyTo(d, 0);
for (int i = pd; i >= 0; --i) {
d[i + places] = d[i];
d[i] = 0.0;
}
return d;
}
static void PolyMultiply(double[] p, double m) {
for (int i = 0; i < p.Length; ++i) {
p[i] *= m;
}
}
static void PolySubtract(double[] p, double[] s) {
for (int i = 0; i < p.Length; ++i) {
p[i] -= s[i];
}
}
static Solution PolyLongDiv(double[] n, double[] d) {
if (n.Length != d.Length) {
throw new ArgumentException("Numerator and denominator vectors must have the same size");
}
int nd = PolyDegree(n);
int dd = PolyDegree(d);
if (dd < 0) {
throw new ArgumentException("Divisor must have at least one one-zero coefficient");
}
if (nd < dd) {
throw new ArgumentException("The degree of the divisor cannot exceed that of the numerator");
}
double[] n2 = new double[n.Length];
n.CopyTo(n2, 0);
double[] q = new double[n.Length];
while (nd >= dd) {
double[] d2 = PolyShiftRight(d, nd - dd);
q[nd - dd] = n2[nd] / d2[nd];
PolyMultiply(d2, q[nd - dd]);
PolySubtract(n2, d2);
nd = PolyDegree(n2);
}
return new Solution(q, n2);
}
static void PolyShow(double[] p) {
int pd = PolyDegree(p);
for (int i = pd; i >= 0; --i) {
double coeff = p[i];
if (coeff == 0.0) continue;
if (coeff == 1.0) {
if (i < pd) {
Console.Write(" + ");
}
} else if (coeff == -1.0) {
if (i < pd) {
Console.Write(" - ");
} else {
Console.Write("-");
}
} else if (coeff < 0.0) {
if (i < pd) {
Console.Write(" - {0:F1}", -coeff);
} else {
Console.Write("{0:F1}", coeff);
}
} else {
if (i < pd) {
Console.Write(" + {0:F1}", coeff);
} else {
Console.Write("{0:F1}", coeff);
}
}
if (i > 1) Console.Write("x^{0}", i);
else if (i == 1) Console.Write("x");
}
Console.WriteLine();
}
static void Main(string[] args) {
double[] n = { -42.0, 0.0, -12.0, 1.0 };
double[] d = { -3.0, 1.0, 0.0, 0.0 };
Console.Write("Numerator : ");
PolyShow(n);
Console.Write("Denominator : ");
PolyShow(d);
Console.WriteLine("-------------------------------------");
Solution sol = PolyLongDiv(n, d);
Console.Write("Quotient : ");
PolyShow(sol.Quotient);
Console.Write("Remainder : ");
PolyShow(sol.Remainder);
}
}
}</syntaxhighlight>
{{out}}
<pre>Numerator : x^3 - 12.0x^2 - 42.0
Denominator : x - 3.0
-------------------------------------
Quotient : x^2 - 9.0x - 27.0
Remainder : -123.0</pre>
=={{header|C++}}==
<syntaxhighlight lang="cpp">
#include <iostream>
#include <iterator>
#include <vector>
using namespace std;
typedef vector<double> Poly;
// does: prints all members of vector
// input: c - ASCII char with the name of the vector
// A - reference to polynomial (vector)
void Print(char name, const Poly &A) {
cout << name << "(" << A.size()-1 << ") = [ ";
copy(A.begin(), A.end(), ostream_iterator<decltype(A[0])>(cout, " "));
cout << "]\n";
}
int main() {
Poly N, D, d, q, r; // vectors - N / D == q && N % D == r
size_t dN, dD, dd, dq, dr; // degrees of vectors
size_t i; // loop counter
// setting the degrees of vectors
cout << "Enter the degree of N: ";
cin >> dN;
cout << "Enter the degree of D: ";
cin >> dD;
dq = dN-dD;
dr = dN-dD;
if( dD < 1 || dN < 1 ) {
cerr << "Error: degree of D and N must be positive.\n";
return 1;
}
// allocation and initialization of vectors
N.resize(dN+1);
cout << "Enter the coefficients of N:"<<endl;
for ( i = 0; i <= dN; i++ ) {
cout << "N[" << i << "]= ";
cin >> N[i];
}
D.resize(dN+1);
cout << "Enter the coefficients of D:"<<endl;
for ( i = 0; i <= dD; i++ ) {
cout << "D[" << i << "]= ";
cin >> D[i];
}
d.resize(dN+1);
q.resize(dq+1);
r.resize(dr+1);
cout << "-- Procedure --" << endl << endl;
if( dN >= dD ) {
while(dN >= dD) {
// d equals D shifted right
d.assign(d.size(), 0);
for( i = 0 ; i <= dD ; i++ )
d[i+dN-dD] = D[i];
dd = dN;
Print( 'd', d );
// calculating one element of q
q[dN-dD] = N[dN]/d[dd];
Print( 'q', q );
// d equals d * q[dN-dD]
for( i = 0 ; i < dq + 1 ; i++ )
d[i] = d[i] * q[dN-dD];
Print( 'd', d );
// N equals N - d
for( i = 0 ; i < dN + 1 ; i++ )
N[i] = N[i] - d[i];
dN--;
Print( 'N', N );
cout << "-----------------------" << endl << endl;
}
}
// r equals N
for( i = 0 ; i <= dN ; i++ )
r[i] = N[i];
cout << "=========================" << endl << endl;
cout << "-- Result --" << endl << endl;
Print( 'q', q );
Print( 'r', r );
}
</syntaxhighlight>
=={{header|Clojure}}==
This example performs ''multivariate'' polynomial division using [https://en.wikipedia.org/wiki/Buchberger%27s_algorithm Buchberger's algorithm] to decompose a polynomial into its [https://en.wikipedia.org/wiki/Gr%C3%B6bner_basis Gröbner bases]. Polynomials are represented as hash-maps of monomials with tuples of exponents as keys and their corresponding coefficients as values: e.g. 2xy + 3x + 5y + 7 is represented as {[1 1] 2, [1 0] 3, [0 1] 5, [0 0] 7}.
Since this algorithm is much more efficient when the input is in [https://en.wikipedia.org/wiki/Monomial_order#Graded_reverse_lexicographic_order graded reverse lexicographic (grevlex) order] a comparator is included to be used with Clojure's sorted-map—<code>(into (sorted-map-by grevlex) ...)</code>—as well as necessary functions to compute polynomial multiplication, monomial complements, and S-polynomials.
<syntaxhighlight lang="clojure">(defn grevlex [term1 term2]
(let [grade1 (reduce +' term1)
grade2 (reduce +' term2)
comp (- grade2 grade1)] ;; total degree
(if (not= 0 comp)
comp
(loop [term1 term1
term2 term2]
(if (empty? term1)
0
(let [grade1 (last term1)
grade2 (last term2)
comp (- grade1 grade2)] ;; differs from grlex because terms are flipped from above
(if (not= 0 comp)
comp
(recur (pop term1)
(pop term2)))))))))
(defn mul
;; transducer
([poly1] ;; completion
(fn
([] poly1)
([poly2] (mul poly1 poly2))
([poly2 & more] (mul poly1 poly2 more))))
([poly1 poly2]
(let [product (atom (transient (sorted-map-by grevlex)))]
(doall ;; `for` is lazy so must to be forced for side-effects
(for [term1 poly1
term2 poly2
:let [vars (mapv +' (key term1) (key term2))
coeff (* (val term1) (val term2))]]
(if (contains? @product vars)
(swap! product assoc! vars (+ (get @product vars) coeff))
(swap! product assoc! vars coeff))))
(->> product
(deref)
(persistent!)
(denull))))
([poly1 poly2 & more]
(reduce mul (mul poly1 poly2) more)))
(defn compl [term1 term2]
(map (fn [x y]
(cond
(and (zero? x) (not= 0 y)) nil
(< x y) nil
(>= x y) (- x y)))
term1
term2))
(defn s-poly [f g]
(let [f-vars (first f)
g-vars (first g)
lcm (compl f-vars g-vars)]
(if (not-any? nil? lcm)
{(vec lcm)
(/ (second f) (second g))})))
(defn divide [f g]
(loop [f f
g g
result (transient {})
remainder {}]
(if (empty? f)
(list (persistent! result)
(->> remainder
(filter #(not (nil? %)))
(into (sorted-map-by grevlex))))
(let [term1 (first f)
term2 (first g)
s-term (s-poly term1 term2)]
(if (nil? s-term)
(recur (dissoc f (first term1))
(dissoc g (first term2))
result
(conj remainder term1))
(recur (sub f (mul g s-term))
g
(conj! result s-term)
remainder))))))
(deftest divide-tests
(is (= (divide {[1 1] 2, [1 0] 3, [0 1] 5, [0 0] 7}
{[1 1] 2, [1 0] 3, [0 1] 5, [0 0] 7})
'({[0 0] 1} {})))
(is (= (divide {[1 1] 2, [1 0] 3, [0 1] 5, [0 0] 7}
{[0 0] 1})
'({[1 1] 2, [1 0] 3, [0 1] 5, [0 0] 7} {})))
(is (= (divide {[1 1] 2, [1 0] 10, [0 1] 3, [0 0] 15}
{[0 1] 1, [0 0] 5})
'({[1 0] 2, [0 0] 3} {})))
(is (= (divide {[1 1] 2, [1 0] 10, [0 1] 3, [0 0] 15}
{[1 0] 2, [0 0] 3})
'({[0 1] 1, [0 0] 5} {}))))</syntaxhighlight>
=={{header|Common Lisp}}==
Polynomials are represented as lists of degree/coefficient pairs ordered by degree (highest degree first), and pairs with zero coefficients can be omitted. <code>Multiply</code> and <code>divide</code> perform long multiplication and long division, respectively. <code>multiply</code> returns one value, the product, and <code>divide</code> returns two, the quotient and the remainder.
<syntaxhighlight lang="lisp">(defun add (p1 p2)
(do ((sum '())) ((and (endp p1) (endp p2)) (nreverse sum))
(let ((pd1 (if (endp p1) -1 (caar p1)))
(pd2 (if (endp p2) -1 (caar p2))))
(multiple-value-bind (c1 c2)
(cond
((> pd1 pd2) (values (cdr (pop p1)) 0))
((< pd1 pd2) (values 0 (cdr (pop p2))))
(t (values (cdr (pop p1)) (cdr (pop p2)))))
(let ((csum (+ c1 c2)))
(unless (zerop csum)
(setf sum (acons (max pd1 pd2) csum sum))))))))
(defun multiply (p1 p2)
(flet ((*p2 (p)
(destructuring-bind (d . c) p
(loop for (pd . pc) in p2
collecting (cons (+ d pd) (* c pc))))))
(reduce 'add (mapcar #'*p2 p1) :initial-value '())))
(defun subtract (p1 p2)
(add p1 (multiply '((0 . -1)) p2)))
(defun divide (dividend divisor &aux (sum '()))
(assert (not (endp divisor)) (divisor)
'division-by-zero
:operation 'divide
:operands (list dividend divisor))
(flet ((floor1 (dividend divisor)
(if (endp dividend) (values '() ())
(destructuring-bind (d1 . c1) (first dividend)
(destructuring-bind (d2 . c2) (first divisor)
(if (> d2 d1) (values '() dividend)
(let* ((quot (list (cons (- d1 d2) (/ c1 c2))))
(rem (subtract dividend (multiply divisor quot))))
(values quot rem))))))))
(loop (multiple-value-bind (quotient remainder)
(floor1 dividend divisor)
(if (endp quotient) (return (values sum remainder))
(setf dividend remainder
sum (add quotient sum)))))))</syntaxhighlight>
The [[wp:Polynomial_long_division#Example|wikipedia example]]:
<syntaxhighlight lang="lisp">> (divide '((3 . 1) (2 . -12) (0 . -42)) ; x^3 - 12x^2 - 42
'((1 . 1) (0 . -3))) ; x - 3
((2 . 1) (1 . -9) (0 . -27)) ; x^2 - 9x - 27
((0 . -123)) ; -123</syntaxhighlight>
=={{header|D}}==
<syntaxhighlight lang="d">import std.stdio, std.range, std.algorithm, std.typecons, std.conv;
Tuple!(double[], double[]) polyDiv(in double[] inN, in double[] inD)
nothrow pure @safe {
// Code smell: a function that does two things.
static int trimAndDegree(T)(ref T[] poly) nothrow pure @safe @nogc {
poly = poly.retro.find!q{ a != b }(0.0).retro;
return poly.length.signed - 1;
}
auto N = inN.dup;
const(double)[] D = inD;
const dD = trimAndDegree(D);
auto dN = trimAndDegree(N);
double[] q;
if (dD < 0)
throw new Error("ZeroDivisionError");
if (dN >= dD) {
q = [0.0].replicate(dN);
while (dN >= dD) {
auto d = [0.0].replicate(dN - dD) ~ D;
immutable mult = q[dN - dD] = N[$ - 1] / d[$ - 1];
d[] *= mult;
N[] -= d[];
dN = trimAndDegree(N);
}
} else
q = [0.0];
return tuple(q, N);
}
int trimAndDegree1(T)(ref T[] poly) nothrow pure @safe @nogc {
poly.length -= poly.retro.countUntil!q{ a != 0 };
return poly.length.signed - 1;
}
void main() {
immutable N = [-42.0, 0.0, -12.0, 1.0];
immutable D = [-3.0, 1.0, 0.0, 0.0];
writefln("%s / %s = %s remainder %s", N, D, polyDiv(N, D)[]);
}</syntaxhighlight>
{{out}}
<pre>[-42, 0, -12, 1] / [-3, 1, 0, 0] = [-27, -9, 1] remainder [-123]</pre>
=={{header|Delphi}}==
{{libheader| System.SysUtils}}
{{Trans|C#}}
<syntaxhighlight lang="delphi">
program Polynomial_long_division;
{$APPTYPE CONSOLE}
uses
System.SysUtils;
type
PPolySolution = ^TPolySolution;
TPolynomio = record
private
class function Degree(p: TPolynomio): Integer; static;
class function ShiftRight(p: TPolynomio; places: Integer): TPolynomio; static;
class function PolyMultiply(p: TPolynomio; m: double): TPolynomio; static;
class function PolySubtract(p, s: TPolynomio): TPolynomio; static;
class function PolyLongDiv(n, d: TPolynomio): PPolySolution; static;
function GetSize: Integer;
public
value: TArray<Double>;
class operator RightShift(p: TPolynomio; b: Integer): TPolynomio;
class operator Multiply(p: TPolynomio; m: double): TPolynomio;
class operator Subtract(p, s: TPolynomio): TPolynomio;
class operator Divide(p, s: TPolynomio): PPolySolution;
class operator Implicit(a: TArray<Double>): TPolynomio;
class operator Implicit(a: TPolynomio): string;
procedure Assign(other: TPolynomio); overload;
procedure Assign(other: TArray<Double>); overload;
property Size: Integer read GetSize;
function ToString: string;
end;
TPolySolution = record
Quotient, Remainder: TPolynomio;
constructor Create(q, r: TPolynomio);
end;
{ TPolynomio }
procedure TPolynomio.Assign(other: TPolynomio);
begin
Assign(other.value);
end;
procedure TPolynomio.Assign(other: TArray<Double>);
begin
SetLength(value, length(other));
for var i := 0 to High(other) do
value[i] := other[i];
end;
class function TPolynomio.Degree(p: TPolynomio): Integer;
begin
var len := high(p.value);
for var i := len downto 0 do
begin
if p.value[i] <> 0.0 then
exit(i);
end;
Result := -1;
end;
class operator TPolynomio.Divide(p, s: TPolynomio): PPolySolution;
begin
Result := PolyLongDiv(p, s);
end;
function TPolynomio.GetSize: Integer;
begin
Result := Length(value);
end;
class operator TPolynomio.Implicit(a: TPolynomio): string;
begin
Result := a.toString;
end;
class operator TPolynomio.Implicit(a: TArray<Double>): TPolynomio;
begin
Result.Assign(a);
end;
class operator TPolynomio.Multiply(p: TPolynomio; m: double): TPolynomio;
begin
Result := TPolynomio.PolyMultiply(p, m);
end;
class function TPolynomio.PolyLongDiv(n, d: TPolynomio): PPolySolution;
var
Solution: TPolySolution;
begin
if length(n.value) <> Length(d.value) then
raise Exception.Create('Numerator and denominator vectors must have the same size');
var nd := Degree(n);
var dd := Degree(d);
if dd < 0 then
raise Exception.Create('Divisor must have at least one one-zero coefficient');
if nd < dd then
raise Exception.Create('The degree of the divisor cannot exceed that of the numerator');
var n2, q: TPolynomio;
n2.Assign(n);
SetLength(q.value, length(n.value));
while nd >= dd do
begin
var d2 := d shr (nd - dd);
q.value[nd - dd] := n2.value[nd] / d2.value[nd];
d2 := d2 * q.value[nd - dd];
n2 := n2 - d2;
nd := Degree(n2);
end;
new(Result);
Result^.Create(q, n2);
end;
class function TPolynomio.PolyMultiply(p: TPolynomio; m: double): TPolynomio;
begin
Result.Assign(p);
for var i := 0 to High(p.value) do
Result.value[i] := p.value[i] * m;
end;
class operator TPolynomio.RightShift(p: TPolynomio; b: Integer): TPolynomio;
begin
Result := TPolynomio.ShiftRight(p, b);
end;
class function TPolynomio.ShiftRight(p: TPolynomio; places: Integer): TPolynomio;
begin
Result.Assign(p);
if places <= 0 then
exit;
var pd := Degree(p);
Result.Assign(p);
for var i := pd downto 0 do
begin
Result.value[i + places] := Result.value[i];
Result.value[i] := 0.0;
end;
end;
class operator TPolynomio.Subtract(p, s: TPolynomio): TPolynomio;
begin
Result := TPolynomio.PolySubtract(p, s);
end;
class function TPolynomio.PolySubtract(p, s: TPolynomio): TPolynomio;
begin
Result.Assign(p);
for var i := 0 to High(p.value) do
Result.value[i] := p.value[i] - s.value[i];
end;
function TPolynomio.ToString: string;
begin
Result := '';
var pd := Degree(self);
for var i := pd downto 0 do
begin
var coeff := value[i];
if coeff = 0.0 then
Continue;
if coeff = 1.0 then
begin
if i < pd then
Result := Result + ' + ';
end
else
begin
if coeff = -1 then
begin
if i < pd then
Result := Result + ' - '
else
Result := Result + '-';
end
else
begin
if coeff < 0.0 then
begin
if i < pd then
Result := Result + format(' - %.1f', [-coeff])
else
Result := Result + format('%.1f', [coeff]);
end
else
begin
if i < pd then
Result := Result + format(' + %.1f', [coeff])
else
Result := Result + format('%.1f', [coeff]);
end;
end;
end;
if i > 1 then
Result := Result + 'x^' + i.tostring
else if i = 1 then
Result := Result + 'x';
end;
end;
{ TPolySolution }
constructor TPolySolution.Create(q, r: TPolynomio);
begin
Quotient.Assign(q);
Remainder.Assign(r);
end;
// Just for force implicitty string conversion
procedure Writeln(s: string);
begin
System.Writeln(s);
end;
var
n, d: TPolynomio;
Solution: PPolySolution;
begin
n := [-42.0, 0.0, -12.0, 1.0];
d := [-3.0, 1.0, 0.0, 0.0];
Write('Numerator : ');
Writeln(n);
Write('Denominator : ');
Writeln(d);
Writeln('-------------------------------------');
Solution := n / d;
Write('Quotient : ');
Writeln(Solution^.Quotient);
Write('Remainder : ');
Writeln(Solution^.Remainder);
FreeMem(Solution, sizeof(TPolySolution));
Readln;
end.</syntaxhighlight>
=={{header|E}}==
{{lines too long|E}}
This program has some unnecessary features contributing to its length:
* It creates polynomial objects rather than performing its operations directly on arrays.
* It includes code for printing polynomials nicely.
* It prints the intermediate steps of the division.
<pre>pragma.syntax("0.9")
pragma.enable("accumulator")
def superscript(x, out) {
if (x >= 10) { superscript(x // 10) }
out.print("⁰¹²³⁴⁵⁶⁷⁸⁹"[x %% 10])
}
def makePolynomial(initCoeffs :List) {
def degree := {
var i := initCoeffs.size() - 1
while (i >= 0 && initCoeffs[i] <=> 0) { i -= 1 }
if (i < 0) { -Infinity } else { i }
}
def coeffs := initCoeffs(0, if (degree == -Infinity) { [] } else { degree + 1 })
def polynomial {
/** Print the polynomial (not necessary for the task) */
to __printOn(out) {
out.print("(λx.")
var first := true
for i in (0..!(coeffs.size())).descending() {
def coeff := coeffs[i]
if (coeff <=> 0) { continue }
out.print(" ")
if (coeff <=> 1 && !(i <=> 0)) {
# no coefficient written if it's 1 and not the constant term
} else if (first) { out.print(coeff)
} else if (coeff > 0) { out.print("+ ", coeff)
} else { out.print("- ", -coeff)
}
if (i <=> 0) { # no x if it's the constant term
} else if (i <=> 1) { out.print("x")
} else { out.print("x"); superscript(i, out)
}
first := false
}
out.print(")")
}
/** Evaluate the polynomial (not necessary for the task) */
to run(x) {
return accum 0 for i => c in coeffs { _ + c * x**i }
}
to degree() { return degree }
to coeffs() { return coeffs }
to highestCoeff() { return coeffs[degree] }
/** Could support another polynomial, but not part of this task.
Computes this * x**power. */
to timesXToThe(power) {
return makePolynomial([0] * power + coeffs)
}
/** Multiply (by a scalar only). */
to multiply(scalar) {
return makePolynomial(accum [] for x in coeffs { _.with(x * scalar) })
}
/** Subtract (by another polynomial only). */
to subtract(other) {
def oc := other.coeffs() :List
return makePolynomial(accum [] for i in 0..(coeffs.size().max(oc.size())) { _.with(coeffs.fetch(i, fn{0}) - oc.fetch(i, fn{0})) })
}
/** Polynomial long division. */
to quotRem(denominator, trace) {
var numerator := polynomial
require(denominator.degree() >= 0)
if (numerator.degree() < denominator.degree()) {
return [makePolynomial([]), denominator]
} else {
var quotientCoeffs := [0] * (numerator.degree() - denominator.degree())
while (numerator.degree() >= denominator.degree()) {
trace.print(" ", numerator, "\n")
def qCoeff := numerator.highestCoeff() / denominator.highestCoeff()
def qPower := numerator.degree() - denominator.degree()
quotientCoeffs with= (qPower, qCoeff)
def d := denominator.timesXToThe(qPower) * qCoeff
trace.print("- ", d, " (= ", denominator, " * ", qCoeff, "x"); superscript(qPower, trace); trace.print(")\n")
numerator -= d
trace.print(" -------------------------- (Quotient so far: ", makePolynomial(quotientCoeffs), ")\n")
}
return [makePolynomial(quotientCoeffs), numerator]
}
}
}
return polynomial
}</pre>
<syntaxhighlight lang="e">def n := makePolynomial([-42, 0, -12, 1])
def d := makePolynomial([-3, 1])
println("Numerator: ", n)
println("Denominator: ", d)
def [q, r] := n.quotRem(d, stdout)
println("Quotient: ", q)
println("Remainder: ", r)</syntaxhighlight>
Output:
Numerator: (λx. x³ - 12x² - 42)
Denominator: (λx. x - 3)
(λx. x³ - 12x² - 42)
- (λx. x³ - 3.0x²) (= (λx. x - 3) * 1.0x²)
-------------------------- (Quotient so far: (λx. x²))
(λx. -9.0x² - 42.0)
- (λx. -9.0x² + 27.0x) (= (λx. x - 3) * -9.0x¹)
-------------------------- (Quotient so far: (λx. x² - 9.0x))
(λx. -27.0x - 42.0)
- (λx. -27.0x + 81.0) (= (λx. x - 3) * -27.0x⁰)
-------------------------- (Quotient so far: (λx. x² - 9.0x - 27.0))
Quotient: (λx. x² - 9.0x - 27.0)
Remainder: (λx. -123.0)
=={{header|Elixir}}==
{{trans|Ruby}}
<syntaxhighlight lang="elixir">defmodule Polynomial do
def division(_, []), do: raise ArgumentError, "denominator is zero"
def division(_, [0]), do: raise ArgumentError, "denominator is zero"
def division(f, g) when length(f) < length(g), do: {[0], f}
def division(f, g) do
{q, r} = division(g, [], f)
if q==[], do: q = [0]
if r==[], do: r = [0]
{q, r}
end
defp division(g, q, r) when length(r) < length(g), do: {q, r}
defp division(g, q, r) do
p = hd(r) / hd(g)
r2 = Enum.zip(r, g)
|> Enum.with_index
|> Enum.reduce(r, fn {{pn,pg},i},acc ->
List.replace_at(acc, i, pn - p * pg)
end)
division(g, q++[p], tl(r2))
end
end
[ { [1, -12, 0, -42], [1, -3] },
{ [1, -12, 0, -42], [1, 1, -3] },
{ [1, 3, 2], [1, 1] },
{ [1, -4, 6, 5, 3], [1, 2, 1] } ]
|> Enum.each(fn {f,g} ->
{q, r} = Polynomial.division(f, g)
IO.puts "#{inspect f} / #{inspect g} => #{inspect q} remainder #{inspect r}"
end)</syntaxhighlight>
{{out}}
<pre>
[1, -12, 0, -42] / [1, -3] => [1.0, -9.0, -27.0] remainder [-123.0]
[1, -12, 0, -42] / [1, 1, -3] => [1.0, -13.0] remainder [16.0, -81.0]
[1, 3, 2] / [1, 1] => [1.0, 2.0] remainder [0.0]
[1, -4, 6, 5, 3] / [1, 2, 1] => [1.0, -6.0, 17.0] remainder [-23.0, -14.0]
</pre>
=={{header|F_Sharp|F#}}==
{{trans|Ocaml}}
<syntaxhighlight lang="fsharp">
let rec shift n l = if n <= 0 then l else shift (n-1) (l @ [0.0])
let rec pad n l = if n <= 0 then l else pad (n-1) (0.0 :: l)
let rec norm = function | 0.0 :: tl -> norm tl | x -> x
let deg l = List.length (norm l) - 1
let zip op p q =
let d = (List.length p) - (List.length q) in
List.map2 op (pad (-d) p) (pad d q)
let polydiv f g =
let rec aux f s q =
let ddif = (deg f) - (deg s) in
if ddif < 0 then (q, f) else
let k = (List.head f) / (List.head s) in
let ks = List.map ((*) k) (shift ddif s) in
let q' = zip (+) q (shift ddif [k])
let f' = norm (List.tail (zip (-) f ks)) in
aux f' s q' in
aux (norm f) (norm g) []
let str_poly l =
let term v p = match (v, p) with
| ( _, 0) -> string v
| (1.0, 1) -> "x"
| ( _, 1) -> (string v) + "*x"
| (1.0, _) -> "x^" + (string p)
| _ -> (string v) + "*x^" + (string p) in
let rec terms = function
| [] -> []
| h :: t ->
if h = 0.0 then (terms t) else (term h (List.length t)) :: (terms t) in
String.concat " + " (terms l)
let _ =
let f,g = [1.0; -4.0; 6.0; 5.0; 3.0], [1.0; 2.0; 1.0] in
let q, r = polydiv f g in
Printf.printf
" (%s) div (%s)\ngives\nquotient:\t(%s)\nremainder:\t(%s)\n"
(str_poly f) (str_poly g) (str_poly q) (str_poly r)
</syntaxhighlight>
{{out}}
<pre>
(x^4 + -4*x^3 + 6*x^2 + 5*x + 3) div (x^2 + 2*x + 1)
gives
quotient: (x^2 + -6*x + 17)
remainder: (-23*x + -14)
</pre>
=={{header|Factor}}==
<syntaxhighlight lang="factor">USE: math.polynomials
{ -42 0 -12 1 } { -3 1 } p/mod ptrim [ . ] bi@</syntaxhighlight>
{{out}}
<pre>
V{ -27 -9 1 }
V{ -123 }
</pre>
=={{header|Fortran}}==
{{works with|Fortran|95 and later}}
<
implicit none
Line 311 ⟶ 1,674:
end subroutine poly_print
end module Polynom</
<
use Polynom
implicit none
Line 330 ⟶ 1,693:
deallocate(q, r)
end program PolyDivTest</
=={{Header|FreeBASIC}}==
<syntaxhighlight lang="freebasic">#define EPS 1.0e-20
type polyterm
degree as uinteger
coeff as double
end type
sub poly_print( P() as double )
dim as string outstr = "", sri
for i as integer = ubound(P) to 0 step -1
if outstr<>"" then
if P(i)>0 then outstr = outstr + " + "
if P(i)<0 then outstr = outstr + " - "
end if
if P(i)=0 then continue for
if abs(P(i))<>1 or i=0 then
if outstr="" then
outstr = outstr + str((P(i)))
else
outstr = outstr + str(abs(P(i)))
end if
end if
if i>0 then outstr=outstr+"x"
sri= str(i)
if i>1 then outstr=outstr + "^" + sri
next i
print outstr
end sub
function lc_deg( B() as double ) as polyterm
'gets the coefficent and degree of the leading term in a polynomial
dim as polyterm ret
for i as uinteger = ubound(B) to 0 step -1
if B(i)<>0 then
ret.degree = i
ret.coeff = B(i)
return ret
end if
next i
return ret
end function
sub poly_multiply( byval k as polyterm, P() as double )
'in-place multiplication of polynomial by a polynomial term
dim i as integer
for i = ubound(P) to k.degree step -1
P(i) = k.coeff*P(i-k.degree)
next i
for i = k.degree-1 to 0 step -1
P(i)=0
next i
end sub
sub poly_subtract( P() as double, Q() as double )
'in place subtraction of one polynomial from another
dim as uinteger deg = ubound(P)
for i as uinteger = 0 to deg
P(i) -= Q(i)
if abs(P(i))<EPS then P(i)=0 'stupid floating point subtraction, grumble grumble
next i
end sub
sub poly_add( P() as double, byval t as polyterm )
'in-place addition of a polynomial term to a polynomial
P(t.degree) += t.coeff
end sub
sub poly_copy( source() as double, target() as double )
for i as uinteger = 0 to ubound(source)
target(i) = source(i)
next i
end sub
sub polydiv( A() as double, B() as double, Q() as double, R() as double )
dim as polyterm s
dim as double sB(0 to ubound(B))
poly_copy A(), R()
dim as uinteger d = ubound(B), degr = lc_deg(R()).degree
dim as double c = lc_deg(B()).coeff
while degr >= d
s.coeff = lc_deg(R()).coeff/c
s.degree = degr - d
poly_add Q(), s
poly_copy B(), sB()
redim preserve sB(0 to s.degree+ubound(sB)) as double
poly_multiply s, sB()
poly_subtract R(), sB()
degr = lc_deg(R()).degree
redim sB(0 to ubound(B))
wend
end sub
dim as double N(0 to 4) = {-42, 0, -12, 1} 'x^3 - 12x^2 - 42
dim as double D(0 to 2) = {-3, 1} ' x - 3
dim as double Q(0 to ubound(N)), R(0 to ubound(N))
polydiv( N(), D(), Q(), R() )
poly_print Q() 'quotient
poly_print R() 'remainder</syntaxhighlight>
{{out}}
<pre>x^2 - 9x - 27
-123</pre>
=={{header|GAP}}==
GAP has built-in functions for computations with polynomials.
<syntaxhighlight lang="gap">x := Indeterminate(Rationals, "x");
p := x^11 + 3*x^8 + 7*x^2 + 3;
q := x^7 + 5*x^3 + 1;
QuotientRemainder(p, q);
# [ x^4+3*x-5, -16*x^4+25*x^3+7*x^2-3*x+8 ]</syntaxhighlight>
=={{header|Go}}==
By the convention and pseudocode given in the task:
<syntaxhighlight lang="go">package main
import "fmt"
func main() {
n := []float64{-42, 0, -12, 1}
d := []float64{-3, 1}
fmt.Println("N:", n)
fmt.Println("D:", d)
q, r, ok := pld(n, d)
if ok {
fmt.Println("Q:", q)
fmt.Println("R:", r)
} else {
fmt.Println("error")
}
}
func degree(p []float64) int {
for d := len(p) - 1; d >= 0; d-- {
if p[d] != 0 {
return d
}
}
return -1
}
func pld(nn, dd []float64) (q, r []float64, ok bool) {
if degree(dd) < 0 {
return
}
nn = append(r, nn...)
if degree(nn) >= degree(dd) {
q = make([]float64, degree(nn)-degree(dd)+1)
for degree(nn) >= degree(dd) {
d := make([]float64, degree(nn)+1)
copy(d[degree(nn)-degree(dd):], dd)
q[degree(nn)-degree(dd)] = nn[degree(nn)] / d[degree(d)]
for i := range d {
d[i] *= q[degree(nn)-degree(dd)]
nn[i] -= d[i]
}
}
}
return q, nn, true
}</syntaxhighlight>
Output:
<pre>
N: [-42 0 -12 1]
D: [-3 1]
Q: [-27 -9 1]
R: [-123 0 0 0]
</pre>
=={{header|Haskell}}==
Translated from the OCaml code elsewhere on the page.
{{works with|GHC|6.10}}
<syntaxhighlight lang="haskell">import Data.List
shift n l = l ++ replicate n 0
pad n l = replicate n 0 ++ l
norm :: Fractional a => [a] -> [a]
norm = dropWhile (== 0)
deg l = length (norm l) - 1
zipWith' op p q = zipWith op (pad (-d) p) (pad d q)
where d = (length p) - (length q)
polydiv f g = aux (norm f) (norm g) []
where aux f s q | ddif < 0 = (q, f)
| otherwise = aux f' s q'
where ddif = (deg f) - (deg s)
k = (head f) / (head s)
ks = map (* k) $ shift ddif s
q' = zipWith' (+) q $ shift ddif [k]
f' = norm $ tail $ zipWith' (-) f ks</syntaxhighlight>
And this is the also-translated pretty printing function.
<syntaxhighlight lang="haskell">str_poly l = intercalate " + " $ terms l
where term v 0 = show v
term 1 1 = "x"
term v 1 = (show v) ++ "x"
term 1 p = "x^" ++ (show p)
term v p = (show v) ++ "x^" ++ (show p)
terms :: Fractional a => [a] -> [String]
terms [] = []
terms (0:t) = terms t
terms (h:t) = (term h (length t)) : (terms t)</syntaxhighlight>
=={{header|J}}==
From http://www.jsoftware.com/jwiki/Phrases/Polynomials
<syntaxhighlight lang="j">divmod=:[: (}: ; {:) ([ (] -/@,:&}. (* {:)) ] , %&{.~)^:(>:@-~&#)&.|.~</syntaxhighlight>
Wikipedia example:
<syntaxhighlight lang="j">_42 0 _12 1 divmod _3 1</syntaxhighlight>
This produces the result:
┌────────┬────┐
│_27 _9 1│_123│
└────────┴────┘
This means that <math>-42-12 x^2+x^3</math> divided by <math>-3+x</math> produces <math>-27-9 x+x^2</math> with a remainder of <math>-123</math>.
=={{header|Java}}==
Replace existing translation.
<br><br>
When generalized, the coefficients of polynomial division are fractions. This implementation supports integer and fraction coefficients.
<br><br>
To test and validate the results, polynomial multiplication and addition are also implemented.
<syntaxhighlight lang="java">
import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.List;
public class PolynomialLongDivision {
public static void main(String[] args) {
RunDivideTest(new Polynomial(1, 3, -12, 2, -42, 0), new Polynomial(1, 1, -3, 0));
RunDivideTest(new Polynomial(5, 2, 4, 1, 1, 0), new Polynomial(2, 1, 3, 0));
RunDivideTest(new Polynomial(5, 10, 4, 7, 1, 0), new Polynomial(2, 4, 2, 2, 3, 0));
RunDivideTest(new Polynomial(2,7,-24,6,2,5,-108,4,3,3,-120,2,-126,0), new Polynomial(2, 4, 2, 2, 3, 0));
}
private static void RunDivideTest(Polynomial p1, Polynomial p2) {
Polynomial[] result = p1.divide(p2);
System.out.printf("Compute: (%s) / (%s) = %s reminder %s%n", p1, p2, result[0], result[1]);
System.out.printf("Test: (%s) * (%s) + (%s) = %s%n%n", result[0], p2, result[1], result[0].multiply(p2).add(result[1]));
}
private static final class Polynomial {
private List<Term> polynomialTerms;
// Format - coeff, exp, coeff, exp, (repeating in pairs) . . .
public Polynomial(long ... values) {
if ( values.length % 2 != 0 ) {
throw new IllegalArgumentException("ERROR 102: Polynomial constructor. Length must be even. Length = " + values.length);
}
polynomialTerms = new ArrayList<>();
for ( int i = 0 ; i < values.length ; i += 2 ) {
polynomialTerms.add(new Term(BigInteger.valueOf(values[i]), values[i+1]));
}
Collections.sort(polynomialTerms, new TermSorter());
}
public Polynomial() {
// zero
polynomialTerms = new ArrayList<>();
polynomialTerms.add(new Term(BigInteger.ZERO, 0));
}
private Polynomial(List<Term> termList) {
if ( termList.size() != 0 ) {
// Remove zero terms if needed
for ( int i = 0 ; i < termList.size() ; i++ ) {
if ( termList.get(i).coefficient.compareTo(Integer.ZERO_INT) == 0 ) {
termList.remove(i);
}
}
}
if ( termList.size() == 0 ) {
// zero
termList.add(new Term(BigInteger.ZERO,0));
}
polynomialTerms = termList;
Collections.sort(polynomialTerms, new TermSorter());
}
public Polynomial[] divide(Polynomial v) {
Polynomial q = new Polynomial();
Polynomial r = this;
Number lcv = v.leadingCoefficient();
long dv = v.degree();
while ( r.degree() >= dv ) {
Number lcr = r.leadingCoefficient();
Number s = lcr.divide(lcv);
Term term = new Term(s, r.degree() - dv);
q = q.add(term);
r = r.add(v.multiply(term.negate()));
}
return new Polynomial[] {q, r};
}
public Polynomial add(Polynomial polynomial) {
List<Term> termList = new ArrayList<>();
int thisCount = polynomialTerms.size();
int polyCount = polynomial.polynomialTerms.size();
while ( thisCount > 0 || polyCount > 0 ) {
Term thisTerm = thisCount == 0 ? null : polynomialTerms.get(thisCount-1);
Term polyTerm = polyCount == 0 ? null : polynomial.polynomialTerms.get(polyCount-1);
if ( thisTerm == null ) {
termList.add(polyTerm);
polyCount--;
}
else if (polyTerm == null ) {
termList.add(thisTerm);
thisCount--;
}
else if ( thisTerm.degree() == polyTerm.degree() ) {
Term t = thisTerm.add(polyTerm);
if ( t.coefficient.compareTo(Integer.ZERO_INT) != 0 ) {
termList.add(t);
}
thisCount--;
polyCount--;
}
else if ( thisTerm.degree() < polyTerm.degree() ) {
termList.add(thisTerm);
thisCount--;
}
else {
termList.add(polyTerm);
polyCount--;
}
}
return new Polynomial(termList);
}
public Polynomial add(Term term) {
List<Term> termList = new ArrayList<>();
boolean added = false;
for ( int index = 0 ; index < polynomialTerms.size() ; index++ ) {
Term currentTerm = polynomialTerms.get(index);
if ( currentTerm.exponent == term.exponent ) {
added = true;
if ( currentTerm.coefficient.add(term.coefficient).compareTo(Integer.ZERO_INT) != 0 ) {
termList.add(currentTerm.add(term));
}
}
else {
termList.add(currentTerm);
}
}
if ( ! added ) {
termList.add(term);
}
return new Polynomial(termList);
}
public Polynomial multiply(Polynomial polynomial) {
List<Term> termList = new ArrayList<>();
for ( int i = 0 ; i < polynomialTerms.size() ; i++ ) {
Term ci = polynomialTerms.get(i);
for ( int j = 0 ; j < polynomial.polynomialTerms.size() ; j++ ) {
Term cj = polynomial.polynomialTerms.get(j);
Term currentTerm = ci.multiply(cj);
boolean added = false;
for ( int k = 0 ; k < termList.size() ; k++ ) {
if ( currentTerm.exponent == termList.get(k).exponent ) {
added = true;
Term t = termList.remove(k).add(currentTerm);
if ( t.coefficient.compareTo(Integer.ZERO_INT) != 0 ) {
termList.add(t);
}
break;
}
}
if ( ! added ) {
termList.add(currentTerm);
}
}
}
return new Polynomial(termList);
}
public Polynomial multiply(Term term) {
List<Term> termList = new ArrayList<>();
for ( int index = 0 ; index < polynomialTerms.size() ; index++ ) {
Term currentTerm = polynomialTerms.get(index);
termList.add(currentTerm.multiply(term));
}
return new Polynomial(termList);
}
public Number leadingCoefficient() {
return polynomialTerms.get(0).coefficient;
}
public long degree() {
return polynomialTerms.get(0).exponent;
}
@Override
public String toString() {
StringBuilder sb = new StringBuilder();
boolean first = true;
for ( Term term : polynomialTerms ) {
if ( first ) {
sb.append(term);
first = false;
}
else {
sb.append(" ");
if ( term.coefficient.compareTo(Integer.ZERO_INT) > 0 ) {
sb.append("+ ");
sb.append(term);
}
else {
sb.append("- ");
sb.append(term.negate());
}
}
}
return sb.toString();
}
}
private static final class TermSorter implements Comparator<Term> {
@Override
public int compare(Term o1, Term o2) {
return (int) (o2.exponent - o1.exponent);
}
}
private static final class Term {
Number coefficient;
long exponent;
public Term(BigInteger c, long e) {
coefficient = new Integer(c);
exponent = e;
}
public Term(Number c, long e) {
coefficient = c;
exponent = e;
}
public Term multiply(Term term) {
return new Term(coefficient.multiply(term.coefficient), exponent + term.exponent);
}
public Term add(Term term) {
if ( exponent != term.exponent ) {
throw new RuntimeException("ERROR 102: Exponents not equal.");
}
return new Term(coefficient.add(term.coefficient), exponent);
}
public Term negate() {
return new Term(coefficient.negate(), exponent);
}
public long degree() {
return exponent;
}
@Override
public String toString() {
if ( coefficient.compareTo(Integer.ZERO_INT) == 0 ) {
return "0";
}
if ( exponent == 0 ) {
return "" + coefficient;
}
if ( coefficient.compareTo(Integer.ONE_INT) == 0 ) {
if ( exponent == 1 ) {
return "x";
}
else {
return "x^" + exponent;
}
}
if ( exponent == 1 ) {
return coefficient + "x";
}
return coefficient + "x^" + exponent;
}
}
private static abstract class Number {
public abstract int compareTo(Number in);
public abstract Number negate();
public abstract Number add(Number in);
public abstract Number multiply(Number in);
public abstract Number inverse();
public abstract boolean isInteger();
public abstract boolean isFraction();
public Number subtract(Number in) {
return add(in.negate());
}
public Number divide(Number in) {
return multiply(in.inverse());
}
}
public static class Fraction extends Number {
private final Integer numerator;
private final Integer denominator;
public Fraction(Integer n, Integer d) {
numerator = n;
denominator = d;
}
@Override
public int compareTo(Number in) {
if ( in.isInteger() ) {
Integer result = ((Integer) in).multiply(denominator);
return numerator.compareTo(result);
}
else if ( in.isFraction() ) {
Fraction inFrac = (Fraction) in;
Integer left = numerator.multiply(inFrac.denominator);
Integer right = denominator.multiply(inFrac.numerator);
return left.compareTo(right);
}
throw new RuntimeException("ERROR: Unknown number type in Fraction.compareTo");
}
@Override
public Number negate() {
if ( denominator.integer.signum() < 0 ) {
return new Fraction(numerator, (Integer) denominator.negate());
}
return new Fraction((Integer) numerator.negate(), denominator);
}
@Override
public Number add(Number in) {
if ( in.isInteger() ) {
//x/y+z = (x+yz)/y
return new Fraction((Integer) ((Integer) in).multiply(denominator).add(numerator), denominator);
}
else if ( in.isFraction() ) {
Fraction inFrac = (Fraction) in;
// compute a/b + x/y
// Let q = gcd(b,y)
// Result = ( (a*y + x*b)/q ) / ( b*y/q )
Integer x = inFrac.numerator;
Integer y = inFrac.denominator;
Integer q = y.gcd(denominator);
Integer temp1 = numerator.multiply(y);
Integer temp2 = denominator.multiply(x);
Integer newDenom = denominator.multiply(y).divide(q);
if ( newDenom.compareTo(Integer.ONE_INT) == 0 ) {
return temp1.add(temp2);
}
Integer newNum = (Integer) temp1.add(temp2).divide(q);
Integer gcd2 = newDenom.gcd(newNum);
if ( gcd2.compareTo(Integer.ONE_INT) == 0 ) {
return new Fraction(newNum, newDenom);
}
newNum = newNum.divide(gcd2);
newDenom = newDenom.divide(gcd2);
if ( newDenom.compareTo(Integer.ONE_INT) == 0 ) {
return newNum;
}
else if ( newDenom.compareTo(Integer.MINUS_ONE_INT) == 0 ) {
return newNum.negate();
}
return new Fraction(newNum, newDenom);
}
throw new RuntimeException("ERROR: Unknown number type in Fraction.compareTo");
}
@Override
public Number multiply(Number in) {
if ( in.isInteger() ) {
//x/y*z = x*z/y
Integer temp = numerator.multiply((Integer) in);
Integer gcd = temp.gcd(denominator);
if ( gcd.compareTo(Integer.ONE_INT) == 0 || gcd.compareTo(Integer.MINUS_ONE_INT) == 0 ) {
return new Fraction(temp, denominator);
}
Integer newTop = temp.divide(gcd);
Integer newBot = denominator.divide(gcd);
if ( newBot.compareTo(Integer.ONE_INT) == 0 ) {
return newTop;
}
if ( newBot.compareTo(Integer.MINUS_ONE_INT) == 0 ) {
return newTop.negate();
}
return new Fraction(newTop, newBot);
}
else if ( in.isFraction() ) {
Fraction inFrac = (Fraction) in;
// compute a/b * x/y
Integer tempTop = numerator.multiply(inFrac.numerator);
Integer tempBot = denominator.multiply(inFrac.denominator);
Integer gcd = tempTop.gcd(tempBot);
if ( gcd.compareTo(Integer.ONE_INT) == 0 || gcd.compareTo(Integer.MINUS_ONE_INT) == 0 ) {
return new Fraction(tempTop, tempBot);
}
Integer newTop = tempTop.divide(gcd);
Integer newBot = tempBot.divide(gcd);
if ( newBot.compareTo(Integer.ONE_INT) == 0 ) {
return newTop;
}
if ( newBot.compareTo(Integer.MINUS_ONE_INT) == 0 ) {
return newTop.negate();
}
return new Fraction(newTop, newBot);
}
throw new RuntimeException("ERROR: Unknown number type in Fraction.compareTo");
}
@Override
public boolean isInteger() {
return false;
}
@Override
public boolean isFraction() {
return true;
}
@Override
public String toString() {
return numerator.toString() + "/" + denominator.toString();
}
@Override
public Number inverse() {
if ( numerator.equals(Integer.ONE_INT) ) {
return denominator;
}
else if ( numerator.equals(Integer.MINUS_ONE_INT) ) {
return denominator.negate();
}
else if ( numerator.integer.signum() < 0 ) {
return new Fraction((Integer) denominator.negate(), (Integer) numerator.negate());
}
return new Fraction(denominator, numerator);
}
}
public static class Integer extends Number {
private BigInteger integer;
public static final Integer MINUS_ONE_INT = new Integer(new BigInteger("-1"));
public static final Integer ONE_INT = new Integer(new BigInteger("1"));
public static final Integer ZERO_INT = new Integer(new BigInteger("0"));
public Integer(BigInteger number) {
this.integer = number;
}
public int compareTo(Integer val) {
return integer.compareTo(val.integer);
}
@Override
public int compareTo(Number in) {
if ( in.isInteger() ) {
return compareTo((Integer) in);
}
else if ( in.isFraction() ) {
Fraction frac = (Fraction) in;
BigInteger result = integer.multiply(frac.denominator.integer);
return result.compareTo(frac.numerator.integer);
}
throw new RuntimeException("ERROR: Unknown number type in Integer.compareTo");
}
@Override
public Number negate() {
return new Integer(integer.negate());
}
public Integer add(Integer in) {
return new Integer(integer.add(in.integer));
}
@Override
public Number add(Number in) {
if ( in.isInteger() ) {
return add((Integer) in);
}
else if ( in.isFraction() ) {
Fraction f = (Fraction) in;
Integer top = f.numerator;
Integer bot = f.denominator;
return new Fraction((Integer) multiply(bot).add(top), bot);
}
throw new RuntimeException("ERROR: Unknown number type in Integer.add");
}
@Override
public Number multiply(Number in) {
if ( in.isInteger() ) {
return multiply((Integer) in);
}
else if ( in.isFraction() ) {
// a * x/y = ax/y
Integer x = ((Fraction) in).numerator;
Integer y = ((Fraction) in).denominator;
Integer temp = (Integer) multiply(x);
Integer gcd = temp.gcd(y);
if ( gcd.compareTo(Integer.ONE_INT) == 0 || gcd.compareTo(Integer.MINUS_ONE_INT) == 0 ) {
return new Fraction(temp, y);
}
Integer newTop = temp.divide(gcd);
Integer newBot = y.divide(gcd);
if ( newBot.compareTo(Integer.ONE_INT) == 0 ) {
return newTop;
}
if ( newBot.compareTo(Integer.MINUS_ONE_INT) == 0 ) {
return newTop.negate();
}
return new Fraction(newTop, newBot);
}
throw new RuntimeException("ERROR: Unknown number type in Integer.add");
}
public Integer gcd(Integer in) {
return new Integer(integer.gcd(in.integer));
}
public Integer divide(Integer in) {
return new Integer(integer.divide(in.integer));
}
public Integer multiply(Integer in) {
return new Integer(integer.multiply(in.integer));
}
@Override
public boolean isInteger() {
return true;
}
@Override
public boolean isFraction() {
return false;
}
@Override
public String toString() {
return integer.toString();
}
@Override
public Number inverse() {
if ( equals(ZERO_INT) ) {
throw new RuntimeException("Attempting to take the inverse of zero in IntegerExpression");
}
else if ( this.compareTo(ONE_INT) == 0 ) {
return ONE_INT;
}
else if ( this.compareTo(MINUS_ONE_INT) == 0 ) {
return MINUS_ONE_INT;
}
return new Fraction(ONE_INT, this);
}
}
}
</syntaxhighlight>
{{out}}
<pre>
Compute: (x^3 - 12x^2 - 42) / (x - 3) = x^2 - 9x - 27 reminder -123
Test: (x^2 - 9x - 27) * (x - 3) + (-123) = x^3 - 12x^2 - 42
Compute: (5x^2 + 4x + 1) / (2x + 3) = 5/2x - 7/4 reminder 25/4
Test: (5/2x - 7/4) * (2x + 3) + (25/4) = 5x^2 + 4x + 1
Compute: (5x^10 + 4x^7 + 1) / (2x^4 + 2x^2 + 3) = 5/2x^6 - 5/2x^4 + 2x^3 - 5/4x^2 - 2x + 5 reminder -2x^3 - 25/4x^2 + 6x - 14
Test: (5/2x^6 - 5/2x^4 + 2x^3 - 5/4x^2 - 2x + 5) * (2x^4 + 2x^2 + 3) + (-2x^3 - 25/4x^2 + 6x - 14) = 5x^10 + 4x^7 + 1
Compute: (2x^7 - 24x^6 + 2x^5 - 108x^4 + 3x^3 - 120x^2 - 126) / (2x^4 + 2x^2 + 3) = x^3 - 12x^2 - 42 reminder 0
Test: (x^3 - 12x^2 - 42) * (2x^4 + 2x^2 + 3) + (0) = 2x^7 - 24x^6 + 2x^5 - 108x^4 + 3x^3 - 120x^2 - 126
</pre>
=={{header|jq}}==
{{works with|jq}}
'''Also works with gojq, the Go implementation of jq, and with fq'''
'''Adapted from the second version in the [[#Wren|Wren]] entry.'''
In this entry, polynomials are represented by JSON arrays exactly as
in the task description; that is, using jq notation, `.[$i]` corresponds to
the coefficient of {\displaystyle x^i}.
<syntaxhighlight lang=jq>
# Emit the canonical form of the polynomical represented by the input array
def canonical:
if length == 0 then .
elif .[-1] == 0 then .[:-1]|canonical
else .
end;
# string representation
def poly2s: "Polynomial(\(join(",")))";
# Polynomial division
# Output [ quotient, remainder]
def divrem($divisor):
($divisor|canonical) as $divisor
| { curr: canonical}
| .base = ((.curr|length) - ($divisor|length))
| until( .base < 0;
(.curr[-1] / $divisor[-1]) as $res
| .result += [$res]
| .curr |= .[0:-1]
| reduce range (0;$divisor|length-1) as $i (.;
.curr[.base + $i] += (- $res * $divisor[$i]) )
| .base += -1
)
| [(.result | reverse), (.curr | canonical)];
def demo($num; $den):
{$num, $den,
res: ($num | divrem($den)) }
| .quot = .res[0]
| .rem = .res[1]
| del(.res)
| map_values(poly2s)
| "\(.num) / \(.den) = \(.quot) remainder \(.rem)";
demo( [-42, 0, -12, 1]; [-3, 1, 0, 0])
</syntaxhighlight>
{{output}}
<pre>
Polynomial(-42,0,-12,1) / Polynomial(-3,1,0,0) = Polynomial(-27,-9,1)
remainder Polynomial(-123)
</pre>
=={{header|Julia}}==
This task is straightforward with the help of Julia's [https://github.com/Keno/Polynomials.jl Polynomials] package.
<syntaxhighlight lang="julia">
using Polynomials
p = Poly([-42,0,-12,1])
q = Poly([-3,1])
d, r = divrem(p,q)
println(p, " divided by ", q, " is ", d, " with remainder ", r, ".")
</syntaxhighlight>
{{out}}
<pre>
-42 - 12x^2 + x^3 divided by -3 + x is -27.0 - 9.0x + x^2 with remainder -123.0.
</pre>
=={{header|Kotlin}}==
===Version 1===
<syntaxhighlight lang="scala">// version 1.1.51
typealias IAE = IllegalArgumentException
data class Solution(val quotient: DoubleArray, val remainder: DoubleArray)
fun polyDegree(p: DoubleArray): Int {
for (i in p.size - 1 downTo 0) {
if (p[i] != 0.0) return i
}
return Int.MIN_VALUE
}
fun polyShiftRight(p: DoubleArray, places: Int): DoubleArray {
if (places <= 0) return p
val pd = polyDegree(p)
if (pd + places >= p.size) {
throw IAE("The number of places to be shifted is too large")
}
val d = p.copyOf()
for (i in pd downTo 0) {
d[i + places] = d[i]
d[i] = 0.0
}
return d
}
fun polyMultiply(p: DoubleArray, m: Double) {
for (i in 0 until p.size) p[i] *= m
}
fun polySubtract(p: DoubleArray, s: DoubleArray) {
for (i in 0 until p.size) p[i] -= s[i]
}
fun polyLongDiv(n: DoubleArray, d: DoubleArray): Solution {
if (n.size != d.size) {
throw IAE("Numerator and denominator vectors must have the same size")
}
var nd = polyDegree(n)
val dd = polyDegree(d)
if (dd < 0) {
throw IAE("Divisor must have at least one one-zero coefficient")
}
if (nd < dd) {
throw IAE("The degree of the divisor cannot exceed that of the numerator")
}
val n2 = n.copyOf()
val q = DoubleArray(n.size) // all elements zero by default
while (nd >= dd) {
val d2 = polyShiftRight(d, nd - dd)
q[nd - dd] = n2[nd] / d2[nd]
polyMultiply(d2, q[nd - dd])
polySubtract(n2, d2)
nd = polyDegree(n2)
}
return Solution(q, n2)
}
fun polyShow(p: DoubleArray) {
val pd = polyDegree(p)
for (i in pd downTo 0) {
val coeff = p[i]
if (coeff == 0.0) continue
print (when {
coeff == 1.0 -> if (i < pd) " + " else ""
coeff == -1.0 -> if (i < pd) " - " else "-"
coeff < 0.0 -> if (i < pd) " - ${-coeff}" else "$coeff"
else -> if (i < pd) " + $coeff" else "$coeff"
})
if (i > 1) print("x^$i")
else if (i == 1) print("x")
}
println()
}
fun main(args: Array<String>) {
val n = doubleArrayOf(-42.0, 0.0, -12.0, 1.0)
val d = doubleArrayOf( -3.0, 1.0, 0.0, 0.0)
print("Numerator : ")
polyShow(n)
print("Denominator : ")
polyShow(d)
println("-------------------------------------")
val (q, r) = polyLongDiv(n, d)
print("Quotient : ")
polyShow(q)
print("Remainder : ")
polyShow(r)
}</syntaxhighlight>
{{out}}
<pre>
Output:
Numerator : x^3 - 12.0x^2 - 42.0
Denominator : x - 3.0
-------------------------------------
Quotient : x^2 - 9.0x - 27.0
Remainder : -123.0
</pre>
<br>
===Version 2===
More succinct version that provides an easy-to-use API.
<syntaxhighlight lang="scala">class Polynom(private vararg val factors: Double) {
operator fun div(divisor: Polynom): Pair<Polynom, Polynom> {
var curr = canonical().factors
val right = divisor.canonical().factors
val result = mutableListOf<Double>()
for (base in curr.size - right.size downTo 0) {
val res = curr.last() / right.last()
result += res
curr = curr.copyOfRange(0, curr.size - 1)
for (i in 0 until right.size - 1)
curr[base + i] -= res * right[i]
}
val quot = Polynom(*result.asReversed().toDoubleArray())
val rem = Polynom(*curr).canonical()
return Pair(quot, rem)
}
private fun canonical(): Polynom {
if (factors.last() != 0.0) return this
for (newLen in factors.size downTo 1)
if (factors[newLen - 1] != 0.0)
return Polynom(*factors.copyOfRange(0, newLen))
return Polynom(factors[0])
}
override fun toString() = "Polynom(${factors.joinToString(" ")})"
}
fun main() {
val num = Polynom(-42.0, 0.0, -12.0, 1.0)
val den = Polynom(-3.0, 1.0, 0.0, 0.0)
val (quot, rem) = num / den
print("$num / $den = $quot remainder $rem")
}</syntaxhighlight>
{{out}}
<pre>
Polynom(-42.0 0.0 -12.0 1.0) / Polynom(-3.0 1.0 0.0 0.0) = Polynom(-27.0 -9.0 1.0) remainder Polynom(-123.0)
</pre>
=={{header|Maple}}==
As Maple is a symbolic computation system, polynomial arithmetic is, of course, provided by the language runtime. The remainder (rem) and quotient (quo) operations each allow for the other to be computed simultaneously by passing an unassigned name as an optional fourth argument. Since rem and quo deal also with multivariate polynomials, the indeterminate is passed as the third argument.
<syntaxhighlight lang="maple">
> p := randpoly( x ); # pick a random polynomial in x
5 4 3 2
p := -56 - 7 x + 22 x - 55 x - 94 x + 87 x
> rem( p, x^2 + 2, x, 'q' ); # remainder
220 + 169 x
> q; # quotient
3 2
-7 x + 22 x - 41 x - 138
> quo( p, x^2 + 2, x, 'r' ); # quotient
3 2
-7 x + 22 x - 41 x - 138
> r; # remainder
220 + 169 x
> expand( (x^2+2)*q + r - p ); # check
0
</syntaxhighlight>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">PolynomialQuotientRemainder[x^3-12 x^2-42,x-3,x]</syntaxhighlight>
output:
<pre>{-27 - 9 x + x^2, -123}</pre>
=={{header|Nim}}==
<syntaxhighlight lang="nim">const MinusInfinity = -1
type
Polynomial = seq[int]
Term = tuple[coeff, exp: int]
func degree(p: Polynomial): int =
## Return the degree of a polynomial.
## "p" is supposed to be normalized.
result = if p.len > 0: p.len - 1 else: MinusInfinity
func normalize(p: var Polynomial) =
## Normalize a polynomial, removing useless zeroes.
while p[^1] == 0: discard p.pop()
func `shr`(p: Polynomial; n: int): Polynomial =
## Shift a polynomial of "n" positions to the right.
result.setLen(p.len + n)
result[n..^1] = p
func `*=`(p: var Polynomial; n: int) =
## Multiply in place a polynomial by an integer.
for item in p.mitems: item *= n
p.normalize()
func `-=`(a: var Polynomial; b: Polynomial) =
## Substract in place a polynomial from another polynomial.
for i, val in b: a[i] -= val
a.normalize()
func longdiv(a, b: Polynomial): tuple[q, r: Polynomial] =
## Compute the long division of a polynomial by another.
## Return the quotient and the remainder as polynomials.
result.r = a
if b.degree < 0: raise newException(DivByZeroDefect, "divisor cannot be zero.")
result.q.setLen(a.len)
while (let k = result.r.degree - b.degree; k >= 0):
var d = b shr k
result.q[k] = result.r[^1] div d[^1]
d *= result.q[k]
result.r -= d
result.q.normalize()
const Superscripts: array['0'..'9', string] = ["⁰", "¹", "²", "³", "⁴", "⁵", "⁶", "⁷", "⁸", "⁹"]
func superscript(n: Natural): string =
## Return the Unicode string to use to represent an exponent.
if n == 1:
return ""
for d in $n:
result.add(Superscripts[d])
func `$`(term: Term): string =
## Return the string representation of a term.
if term.coeff == 0: "0"
elif term.exp == 0: $term.coeff
else:
let base = 'x' & superscript(term.exp)
if term.coeff == 1: base
elif term.coeff == -1: '-' & base
else: $term.coeff & base
func `$`(poly: Polynomial): string =
## return the string representation of a polynomial.
for idx in countdown(poly.high, 0):
let coeff = poly[idx]
var term: Term = (coeff: coeff, exp: idx)
if result.len == 0:
result.add $term
else:
if coeff > 0:
result.add '+'
result.add $term
elif coeff < 0:
term.coeff = -term.coeff
result.add '-'
result.add $term
const
N = @[-42, 0, -12, 1]
D = @[-3, 1]
let (q, r) = longdiv(N, D)
echo "N = ", N
echo "D = ", D
echo "q = ", q
echo "r = ", r</syntaxhighlight>
{{out}}
<pre>N = x³-12x²-42
D = x-3
q = x²-9x-27
r = -123</pre>
=={{header|OCaml}}==
First define some utility operations on polynomials as lists (with highest power coefficient first).
<syntaxhighlight lang="ocaml">let rec shift n l = if n <= 0 then l else shift (pred n) (l @ [0.0])
let rec pad n l = if n <= 0 then l else pad (pred n) (0.0 :: l)
let rec norm = function | 0.0 :: tl -> norm tl | x -> x
let deg l = List.length (norm l) - 1
let zip op p q =
let
Then the main polynomial division function
<syntaxhighlight lang="ocaml">let polydiv f g =
let rec aux f s q =
let
if
let k = (List.hd f) /. (List.hd s) in
let
let
and f' = norm (List.tl (zip (-.) f ks)) in
aux f' s q' in
aux (norm f) (norm g) []</syntaxhighlight>
For output we need a pretty-printing function
<syntaxhighlight lang="ocaml">let str_poly l =
let term v p = match (v, p) with
| ( _, 0) -> string_of_float v
Line 376 ⟶ 2,867:
| h :: t ->
if h = 0.0 then (terms t) else (term h (List.length t)) :: (terms t) in
String.concat " + " (terms l)</syntaxhighlight>
and then the example
<syntaxhighlight lang="ocaml">let _ =
let f = [1.0; -4.0; 6.0; 5.0; 3.0] and g = [1.0; 2.0; 1.0] in
let q, r = polydiv f g in
Printf.printf
" (%s) div (%s)\ngives\nquotient:\t(%s)\nremainder:\t(%s)\n"
gives the output:
<pre>
gives
quotient: (x^2 + -6.*x + 17.)
</pre>
=={{header|Octave}}==
Octave has already facilities to divide two polynomials (<code>deconv(n,d)</code>); and the reason to adopt the convention of keeping the
<
gd = length(d);
pv = zeros(1, length(n));
Line 436 ⟶ 2,922:
[q, r] = poly_long_div([1,3], [1,-12,0,-42]);
polyout(q, 'x');
polyout(r, 'x');</
=={{header|PARI/GP}}==
This uses the built-in PARI polynomials.
<syntaxhighlight lang="parigp">poldiv(a,b)={
my(rem=a%b);
[(a - rem)/b, rem]
};
poldiv(x^9+1, x^3+x-3)</syntaxhighlight>
Alternately, use the built-in function <code>divrem</code>:
<syntaxhighlight lang="parigp">divrem(x^9+1, x^3+x-3)~</syntaxhighlight>
=={{header|Perl}}==
This solution keeps the highest power coefficient first, like [[Polynomial long division#OCaml|OCaml solution]] and [[Polynomial long division#Octave|Octave solution]].
{{trans|Octave}}
<syntaxhighlight lang="perl">use strict;
use List::Util qw(min);
sub poly_long_div
{
my ($rn, $rd) = @_;
my @n = @$rn;
my $gd = scalar(@$rd);
if ( scalar(@n) >= $gd ) {
my @q = ();
while ( scalar(@n) >= $gd ) {
my $piv = $n[0]/$rd->[0];
push @q, $piv;
$n[$_] -= $rd->[$_] * $piv foreach ( 0 .. min(scalar(@n), $gd)-1 );
shift @n;
}
return ( \@q, \@n );
} else {
return ( [0], $rn );
}
}</syntaxhighlight>
<syntaxhighlight lang="perl">sub poly_print
{
my @c = @_;
my $l = scalar(@c);
for(my $i=0; $i < $l; $i++) {
print $c[$i];
print "x^" . ($l-$i-1) . " + " if ($i < ($l-1));
}
print "\n";
}</syntaxhighlight>
<syntaxhighlight lang="perl">my ($q, $r);
($q, $r) = poly_long_div([1, -12, 0, -42], [1, -3]);
poly_print(@$q);
poly_print(@$r);
print "\n";
($q, $r) = poly_long_div([1,-12,0,-42], [1,1,-3]);
poly_print(@$q);
poly_print(@$r);
print "\n";
($q, $r) = poly_long_div([1,3,2], [1,1]);
poly_print(@$q);
poly_print(@$r);
print "\n";
# the example from the OCaml solution
($q, $r) = poly_long_div([1,-4,6,5,3], [1,2,1]);
poly_print(@$q);
poly_print(@$r);</syntaxhighlight>
=={{header|Phix}}==
<!--(phixonline)-->
<syntaxhighlight lang="phix">
-- demo\rosetta\Polynomial_long_division.exw
with javascript_semantics
function degree(sequence p)
for i=length(p) to 1 by -1 do
if p[i]!=0 then return i end if
end for
return -1
end function
function poly_div(sequence n, d)
d = deep_copy(d)
while length(d)<length(n) do d &=0 end while
integer dn = degree(n),
dd = degree(d)
if dd<0 then throw("divide by zero") end if
sequence quo = repeat(0,dn),
rem = deep_copy(n)
while dn>=dd do
integer k = dn-dd, qk = rem[dn]/d[dd]
sequence d2 = d[1..length(d)-k]
quo[k+1] = qk
for i=1 to length(d2) do
integer mi = -i
rem[mi] -= d2[mi]*qk
end for
dn = degree(rem)
end while
return {quo,rem}
end function
function poly(sequence si)
-- display helper
string r = ""
for t=length(si) to 1 by -1 do
integer sit = si[t]
if sit!=0 then
if sit=1 and t>1 then
r &= iff(r=""? "":" + ")
elsif sit=-1 and t>1 then
r &= iff(r=""?"-":" - ")
else
if r!="" then
r &= iff(sit<0?" - ":" + ")
sit = abs(sit)
end if
r &= sprintf("%d",sit)
end if
r &= iff(t>1?"x"&iff(t>2?sprintf("^%d",t-1):""):"")
end if
end for
if r="" then r="0" end if
return r
end function
constant tests = {{{-42,0,-12,1},{-3,1}},
{{-3,1},{-42,0,-12,1}},
{{-42,0,-12,1},{-3,1,1}},
{{2,3,1},{1,1}},
{{3,5,6,-4,1},{1,2,1}},
{{3,0,7,0,0,0,0,0,3,0,0,1},{1,0,0,5,0,0,0,1}},
{{-56,87,-94,-55,22,-7},{2,0,1}},
}
constant fmt = "%40s / %-16s = %25s rem %s\n"
for i=1 to length(tests) do
sequence {num,den} = tests[i],
{quo,rem} = poly_div(num,den)
printf(1,fmt,apply({num,den,quo,rem},poly))
end for
</syntaxhighlight>
{{out}}
<pre style="font-size: 12px">
x^3 - 12x^2 - 42 / x - 3 = x^2 - 9x - 27 rem -123
x - 3 / x^3 - 12x^2 - 42 = 0 rem x - 3
x^3 - 12x^2 - 42 / x^2 + x - 3 = x - 13 rem 16x - 81
x^2 + 3x + 2 / x + 1 = x + 2 rem 0
x^4 - 4x^3 + 6x^2 + 5x + 3 / x^2 + 2x + 1 = x^2 - 6x + 17 rem -23x - 14
x^11 + 3x^8 + 7x^2 + 3 / x^7 + 5x^3 + 1 = x^4 + 3x - 5 rem -16x^4 + 25x^3 + 7x^2 - 3x + 8
-7x^5 + 22x^4 - 55x^3 - 94x^2 + 87x - 56 / x^2 + 2 = -7x^3 + 22x^2 - 41x - 138 rem 169x + 220
</pre>
=={{header|PicoLisp}}==
<syntaxhighlight lang="picolisp">(de degree (P)
(let I NIL
(for (N . C) P
(or (=0 C) (setq I N)) )
(dec I) ) )
(de divPoly (N D)
(if (lt0 (degree D))
(quit "Div/0" D)
(let (Q NIL Diff)
(while (ge0 (setq Diff (- (degree N) (degree D))))
(setq Q (need (- -1 Diff) Q 0))
(let E D
(do Diff (push 'E 0))
(let F (/ (get N (inc (degree N))) (get E (inc (degree E))))
(set (nth Q (inc Diff)) F)
(setq N (mapcar '((N E) (- N (* E F))) N E)) ) ) )
(list Q N) ) ) )</syntaxhighlight>
Output:
<pre>: (divPoly (-42 0 -12 1) (-3 1 0 0))
-> ((-27 -9 1) (-123 0 0 0))</pre>
=={{header|Python}}==
{{works with|Python 2.x}}
<
from itertools import izip
def degree(poly):
return
def poly_div(N, D):
Line 474 ⟶ 3,134:
D = [-3, 1, 0, 0]
print " %s / %s =" % (N,D),
print " %s remainder %s" % poly_div(N, D)</
Sample output:
<pre>POLYNOMIAL LONG DIVISION
[-42, 0, -12, 1] / [-3, 1, 0, 0] = [-27.0, -9.0, 1.0] remainder [-123.0]</pre>
=={{header|R}}==
{{trans|Octave}}
<syntaxhighlight lang="r">polylongdiv <- function(n,d) {
gd <- length(d)
pv <- vector("numeric", length(n))
pv[1:gd] <- d
if ( length(n) >= gd ) {
q <- c()
while ( length(n) >= gd ) {
q <- c(q, n[1]/pv[1])
n <- n - pv * (n[1]/pv[1])
n <- n[2:length(n)]
pv <- pv[1:(length(pv)-1)]
}
list(q=q, r=n)
} else {
list(q=c(0), r=n)
}
}
# an utility function to print polynomial
print.polynomial <- function(p) {
i <- length(p)-1
for(a in p) {
if ( i == 0 ) {
cat(a, "\n")
} else {
cat(a, "x^", i, " + ", sep="")
}
i <- i - 1
}
}
r <- polylongdiv(c(1,-12,0,-42), c(1,-3))
print.polynomial(r$q)
print.polynomial(r$r)</syntaxhighlight>
=={{header|Racket}}==
<syntaxhighlight lang="racket">
#lang racket
(define (deg p)
(for/fold ([d -inf.0]) ([(pi i) (in-indexed p)])
(if (zero? pi) d i)))
(define (lead p) (vector-ref p (deg p)))
(define (mono c d) (build-vector (+ d 1) (λ(i) (if (= i d) c 0))))
(define (poly*cx^n c n p) (vector-append (make-vector n 0) (for/vector ([pi p]) (* c pi))))
(define (poly+ p q) (poly/lin 1 p 1 q))
(define (poly- p q) (poly/lin 1 p -1 q))
(define (poly/lin a p b q)
(cond [(< (deg p) 0) q]
[(< (deg q) 0) p]
[(< (deg p) (deg q)) (poly/lin b q a p)]
[else (define ap+bq (for/vector #:length (+ (deg p) 1) #:fill 0
([pi p] [qi q]) (+ (* a pi) (* b qi))))
(for ([i (in-range (+ (deg q) 1) (+ (deg p) 1))])
(vector-set! ap+bq i (* a (vector-ref p i))))
ap+bq]))
(define (poly/ n d)
(define N (deg n))
(define D (deg d))
(cond
[(< N 0) (error 'poly/ "can't divide by zero")]
[(< N D) (values 0 n)]
[else (define c (/ (lead n) (lead d)))
(define q (mono c (- N D)))
(define r (poly- n (poly*cx^n c (- N D) d)))
(define-values (q1 r1) (poly/ r d))
(values (poly+ q q1) r1)]))
; Example:
(poly/ #(-42 0 -12 1) #(-3 1))
; Output:
'#(-27 -9 1)
'#(-123 0)
</syntaxhighlight>
=={{header|Raku}}==
(formerly Perl 6)
{{Works with|rakudo|2018.10}}
{{trans|Perl}} for the core algorithm; original code for LaTeX pretty-printing.
<syntaxhighlight lang="raku" line>sub poly_long_div ( @n is copy, @d ) {
return [0], |@n if +@n < +@d;
my @q = gather while +@n >= +@d {
@n = @n Z- flat ( ( @d X* take ( @n[0] / @d[0] ) ), 0 xx * );
@n.shift;
}
return @q, @n;
}
sub xP ( $power ) { $power>1 ?? "x^$power" !! $power==1 ?? 'x' !! '' }
sub poly_print ( @c ) { join ' + ', @c.kv.map: { $^v ~ xP( @c.end - $^k ) } }
my @polys = [ [ 1, -12, 0, -42 ], [ 1, -3 ] ],
[ [ 1, -12, 0, -42 ], [ 1, 1, -3 ] ],
[ [ 1, 3, 2 ], [ 1, 1 ] ],
[ [ 1, -4, 6, 5, 3 ], [ 1, 2, 1 ] ];
say '<math>\begin{array}{rr}';
for @polys -> [ @a, @b ] {
printf Q"%s , & %s \\\\\n", poly_long_div( @a, @b ).map: { poly_print($_) };
}
say '\end{array}</math>';</syntaxhighlight>
Output:
<math>\begin{array}{rr}
1x^2 + -9x + -27 , & -123 \\
1x + -13 , & 16x + -81 \\
1x + 2 , & 0 \\
1x^2 + -6x + 17 , & -23x + -14 \\
\end{array}</math>
=={{header|REXX}}==
<syntaxhighlight lang="rexx">/* REXX needed by some... */
z='1 -12 0 -42' /* Numerator */
n='1 -3' /* Denominator */
zx=z
nx=n copies('0 ',words(z)-words(n))
qx='' /* Quotient */
Do Until words(zx)<words(n)
Parse Value div(zx,nx) With q zx
qx=qx q
nx=subword(nx,1,words(nx)-1)
End
Say '('show(z)')/('show(n)')=('show(qx)')'
Say 'Remainder:' show(zx)
Exit
div: Procedure
Parse Arg z,n
q=word(z,1)/word(n,1)
zz=''
Do i=1 To words(z)
zz=zz word(z,i)-q*word(n,i)
End
Return q subword(zz,2)
show: Procedure
Parse Arg poly
d=words(poly)-1
res=''
Do i=1 To words(poly)
Select
When d>1 Then fact='*x**'d
When d=1 Then fact='*x'
Otherwise fact=''
End
Select
When word(poly,i)=0 Then p=''
When word(poly,i)=1 Then p='+'substr(fact,2)
When word(poly,i)=-1 Then p='-'substr(fact,2)
When word(poly,i)<0 Then p=word(poly,i)||fact
Otherwise p='+'word(poly,i)||fact
End
res=res p
d=d-1
End
Return strip(space(res,0),'L','+')</syntaxhighlight>
{{out}}
<pre>(x**3-12*x**2-42)/(x-3)=(x**2-9*x-27)
Remainder: -123 </pre>
=={{header|RPL}}==
{{works with|HP|49}}
'-42-12*X^2+X^3' 'X-3' DIV2
{{out}}
<pre>
2: 'X^2-9*X-27'
1: -123
</pre>
=={{header|Ruby}}==
Implementing the algorithm given in the task description:
<syntaxhighlight lang="ruby">def polynomial_long_division(numerator, denominator)
dd = degree(denominator)
raise ArgumentError, "denominator is zero" if dd < 0
if dd == 0
return [multiply(numerator, 1.0/denominator[0]), [0]*numerator.length]
end
q = [0] * numerator.length
while (dn = degree(numerator)) >= dd
d = shift_right(denominator, dn - dd)
q[dn-dd] = numerator[dn] / d[degree(d)]
d = multiply(d, q[dn-dd])
numerator = subtract(numerator, d)
end
[q, numerator]
end
def degree(ary)
idx = ary.rindex(&:nonzero?)
idx ? idx : -1
end
def shift_right(ary, n)
[0]*n + ary[0, ary.length - n]
end
def subtract(a1, a2)
a1.zip(a2).collect {|v1,v2| v1 - v2}
end
def multiply(ary, num)
ary.collect {|x| x * num}
end
f = [-42, 0, -12, 1]
g = [-3, 1, 0, 0]
q, r = polynomial_long_division(f, g)
puts "#{f} / #{g} => #{q} remainder #{r}"
# => [-42, 0, -12, 1] / [-3, 1, 0, 0] => [-27, -9, 1, 0] remainder [-123, 0, 0, 0]
g = [-3, 1, 1, 0]
q, r = polynomial_long_division(f, g)
puts "#{f} / #{g} => #{q} remainder #{r}"
# => [-42, 0, -12, 1] / [-3, 1, 1, 0] => [-13, 1, 0, 0] remainder [-81, 16, 0, 0]</syntaxhighlight>
Implementing the algorithms on the [[wp:Polynomial long division|wikipedia page]] -- uglier code but nicer user interface
<syntaxhighlight lang="ruby">def polynomial_division(f, g)
if g.length == 0 or (g.length == 1 and g[0] == 0)
raise ArgumentError, "denominator is zero"
elsif g.length == 1
[f.collect {|x| Float(x)/g[0]}, [0]]
elsif g.length == 2
synthetic_division(f, g)
else
higher_degree_synthetic_division(f, g)
end
end
def synthetic_division(f, g)
board = [f] << Array.new(f.length) << Array.new(f.length)
board[2][0] = board[0][0]
1.upto(f.length - 1).each do |i|
board[1][i] = board[2][i-1] * -g[1]
board[2][i] = board[0][i] + board[1][i]
end
[board[2][0..-2], [board[2][-1]]]
end
# an ugly mess of array index arithmetic
# http://en.wikipedia.org/wiki/Polynomial_long_division#Higher_degree_synthetic_division
def higher_degree_synthetic_division(f, g)
# [use] the negative coefficients of the denominator following the leading term
lhs = g[1..-1].collect {|x| -x}
board = [f]
q = []
1.upto(f.length - lhs.length).each do |i|
n = 2*i - 1
# underline the leading coefficient of the right-hand side, multiply it by
# the left-hand coefficients and write the products beneath the next columns
# on the right.
q << board[n-1][i-1]
board << Array.new(f.length).fill(0, i) # row n
(lhs.length).times do |j|
board[n][i+j] = q[-1]*lhs[j]
end
# perform an addition
board << Array.new(f.length).fill(0, i) # row n+1
(lhs.length + 1).times do |j|
board[n+1][i+j] = board[n-1][i+j] + board[n][i+j] if i+j < f.length
end
end
# the remaining numbers in the bottom row correspond to the coefficients of the remainder
r = board[-1].compact
q = [0] if q.empty?
[q, r]
end
f = [1, -12, 0, -42]
g = [1, -3]
q, r = polynomial_division(f, g)
puts "#{f} / #{g} => #{q} remainder #{r}"
# => [1, -12, 0, -42] / [1, -3] => [1, -9, -27] remainder [-123]
g = [1, 1, -3]
q, r = polynomial_division(f, g)
puts "#{f} / #{g} => #{q} remainder #{r}"
# => [1, -12, 0, -42] / [1, 1, -3] => [1, -13] remainder [16, -81]</syntaxhighlight>
Best of both worlds: {{trans|Tcl}}
<syntaxhighlight lang="ruby">def polynomial_division(f, g)
if g.length == 0 or (g.length == 1 and g[0] == 0)
raise ArgumentError, "denominator is zero"
end
return [[0], f] if f.length < g.length
q, n = [], f.dup
while n.length >= g.length
q << Float(n[0]) / g[0]
n[0, g.length].zip(g).each_with_index do |pair, i|
n[i] = pair[0] - q[-1] * pair[1]
end
n.shift
end
q = [0] if q.empty?
n = [0] if n.empty?
[q, n]
end
f = [1, -12, 0, -42]
g = [1, -3]
q, r = polynomial_division(f, g)
puts "#{f} / #{g} => #{q} remainder #{r}"
# => [1, -12, 0, -42] / [1, -3] => [1.0, -9.0, -27.0] remainder [-123.0]
g = [1, 1, -3]
q, r = polynomial_division(f, g)
puts "#{f} / #{g} => #{q} remainder #{r}"
# => [1, -12, 0, -42] / [1, 1, -3] => [1.0, -13.0] remainder [16.0, -81.0]</syntaxhighlight>
=={{header|Sidef}}==
{{trans|Perl}}
<syntaxhighlight lang="ruby">func poly_long_div(rn, rd) {
var n = rn.map{_}
var gd = rd.len
if (n.len >= gd) {
return(gather {
while (n.len >= gd) {
var piv = n[0]/rd[0]
take(piv)
{ |i|
n[i] -= (rd[i] * piv)
} << ^(n.len `min` gd)
n.shift
}
}, n)
}
return([0], rn)
}</syntaxhighlight>
Example:
<syntaxhighlight lang="ruby">func poly_print(c) {
var l = c.len
c.each_kv {|i, n|
print n
print("x^", (l - i - 1), " + ") if (i < l-1)
}
print "\n";
}
var poly = [
Pair([1,-12,0,-42], [1, -3]),
Pair([1,-12,0,-42], [1,1,-3]),
Pair( [1,3,2], [1,1]),
Pair( [1,-4,6,5,3], [1,2,1]),
]
poly.each { |pair|
var (q, r) = poly_long_div(pair.first, pair.second)
poly_print(q)
poly_print(r)
print "\n"
}</syntaxhighlight>
{{out}}
<pre>
1x^2 + -9x^1 + -27
-123
1x^1 + -13
16x^1 + -81
1x^1 + 2
0
1x^2 + -6x^1 + 17
-23x^1 + -14
</pre>
=={{header|Slate}}==
<syntaxhighlight lang="slate">define: #Polynomial &parents: {Comparable} &slots: {#coefficients -> ExtensibleArray new}.
p@(Polynomial traits) new &capacity: n
[
p cloneSettingSlots: #(coefficients) to: {p coefficients new &capacity: n}
].
p@(Polynomial traits) newFrom: seq@(Sequence traits)
[
p clone `>> [coefficients: (seq as: p coefficients). normalize. ]
].
p@(Polynomial traits) copy
[
p cloneSettingSlots: #(coefficients) to: {p coefficients copy}
].
p1@(Polynomial traits) >= p2@(Polynomial traits)
[p1 degree >= p2 degree].
p@(Polynomial traits) degree
[p coefficients indexOfLastSatisfying: [| :n | n isZero not]].
p@(Polynomial traits) normalize
[
[p degree isPositive /\ [p coefficients last isZero]]
whileTrue: [p coefficients removeLast]
].
p@(Polynomial traits) * n@(Number traits)
[
p newFrom: (p coefficients collect: [| :x | x * n])
].
p@(Polynomial traits) / n@(Number traits)
[
p newFrom: (p coefficients collect: [| :x | x / n])
].
p1@(Polynomial traits) minusCoefficients: p2@(Polynomial traits)
[
p1 newFrom: (p1 coefficients with: p2 coefficients collect: #- `er)
].
p@(Polynomial traits) / denom@(Polynomial traits)
[
p >= denom
ifTrue:
[| n q |
n: p copy.
q: p new.
[n >= denom]
whileTrue:
[| piv |
piv: p coefficients last / denom coefficients last.
q coefficients add: piv.
n: (n minusCoefficients: denom * piv).
n normalize].
n coefficients isEmpty ifTrue: [n coefficients add: 0].
{q. n}]
ifFalse: [{p newFrom: #(0). p copy}]
].</syntaxhighlight>
=={{header|Smalltalk}}==
{{works with|GNU Smalltalk}}
<syntaxhighlight lang="smalltalk">Object subclass: Polynomial [
|coeffs|
Polynomial class >> new [ ^ super basicNew init ]
init [ coeffs := OrderedCollection new. ^ self ]
Polynomial class >> newWithCoefficients: coefficients [
|r|
r := super basicNew.
^ r initWithCoefficients: coefficients
]
initWithCoefficients: coefficients [
coeffs := coefficients asOrderedCollection.
^ self
]
/ denominator [ |n q|
n := self deepCopy.
self >= denominator
ifTrue: [
q := Polynomial new.
[ n >= denominator ]
whileTrue: [ |piv|
piv := (n coeff: 0) / (denominator coeff: 0).
q addCoefficient: piv.
n := n - (denominator * piv).
n clean
].
^ { q . (n degree) > 0 ifTrue: [ n ] ifFalse: [ n addCoefficient: 0. n ] }
]
ifFalse: [
^ { Polynomial newWithCoefficients: #( 0 ) . self deepCopy }
]
]
* constant [ |r| r := self deepCopy.
1 to: (coeffs size) do: [ :i |
r at: i put: ((r at: i) * constant)
].
^ r
]
at: index [ ^ coeffs at: index ]
at: index put: obj [ ^ coeffs at: index put: obj ]
>= anotherPoly [
^ (self degree) >= (anotherPoly degree)
]
degree [ ^ coeffs size ]
- anotherPoly [ "This is not a real subtraction between Polynomial: it is an
internal method ..."
|a|
a := self deepCopy.
1 to: ( (coeffs size) min: (anotherPoly degree) ) do: [ :i |
a at: i put: ( (a at: i) - (anotherPoly at: i) )
].
^ a
]
coeff: index [ ^ coeffs at: (index + 1) ]
addCoefficient: coeff [ coeffs add: coeff ]
clean [
[ (coeffs size) > 0
ifTrue: [ (coeffs at: 1) = 0 ] ifFalse: [ false ] ]
whileTrue: [ coeffs removeFirst ].
]
display [
1 to: (coeffs size) do: [ :i |
(coeffs at: i) display.
i < (coeffs size)
ifTrue: [ ('x^%1 + ' % {(coeffs size) - i} ) display ]
]
]
displayNl [ self display. Character nl display ]
].</syntaxhighlight>
<syntaxhighlight lang="smalltalk">|res|
res := OrderedCollection new.
res add: ((Polynomial newWithCoefficients: #( 1 -12 0 -42) ) /
(Polynomial newWithCoefficients: #( 1 -3 ) )) ;
add: ((Polynomial newWithCoefficients: #( 1 -12 0 -42) ) /
(Polynomial newWithCoefficients: #( 1 1 -3 ) )).
res do: [ :o |
(o at: 1) display. ' with rest: ' display. (o at: 2) displayNl
]</syntaxhighlight>
=={{header|SPAD}}==
{{works with|FriCAS}}
{{works with|OpenAxiom}}
{{works with|Axiom}}
<syntaxhighlight lang="spad">(1) -> monicDivide(x^3-12*x^2-42,x-3,'x)
2
(1) [quotient = x - 9x - 27,remainder = - 123]
Type: Record(quotient: Polynomial(Integer),remainder: Polynomial(Integer))</syntaxhighlight>
Domain:[http://fricas.github.io/api/PolynomialCategory.html#l-polynomial-category-monic-divide]
=={{header|Swift}}==
{{trans|Kotlin}}
<syntaxhighlight lang="swift">protocol Dividable {
static func / (lhs: Self, rhs: Self) -> Self
}
extension Int: Dividable { }
struct Solution<T> {
var quotient: [T]
var remainder: [T]
}
func polyDegree<T: SignedNumeric>(_ p: [T]) -> Int {
for i in stride(from: p.count - 1, through: 0, by: -1) where p[i] != 0 {
return i
}
return Int.min
}
func polyShiftRight<T: SignedNumeric>(p: [T], places: Int) -> [T] {
guard places > 0 else {
return p
}
let deg = polyDegree(p)
assert(deg + places < p.count, "Number of places to shift too large")
var res = p
for i in stride(from: deg, through: 0, by: -1) {
res[i + places] = res[i]
res[i] = 0
}
return res
}
func polyMul<T: SignedNumeric>(_ p: inout [T], by: T) {
for i in 0..<p.count {
p[i] *= by
}
}
func polySub<T: SignedNumeric>(_ p: inout [T], by: [T]) {
for i in 0..<p.count {
p[i] -= by[i]
}
}
func polyLongDiv<T: SignedNumeric & Dividable>(numerator n: [T], denominator d: [T]) -> Solution<T>? {
guard n.count == d.count else {
return nil
}
var nDeg = polyDegree(n)
let dDeg = polyDegree(d)
guard dDeg >= 0, nDeg >= dDeg else {
return nil
}
var n2 = n
var quo = [T](repeating: 0, count: n.count)
while nDeg >= dDeg {
let i = nDeg - dDeg
var d2 = polyShiftRight(p: d, places: i)
quo[i] = n2[nDeg] / d2[nDeg]
polyMul(&d2, by: quo[i])
polySub(&n2, by: d2)
nDeg = polyDegree(n2)
}
return Solution(quotient: quo, remainder: n2)
}
func polyPrint<T: SignedNumeric & Comparable>(_ p: [T]) {
let deg = polyDegree(p)
for i in stride(from: deg, through: 0, by: -1) where p[i] != 0 {
let coeff = p[i]
switch coeff {
case 1 where i < deg:
print(" + ", terminator: "")
case 1:
print("", terminator: "")
case -1 where i < deg:
print(" - ", terminator: "")
case -1:
print("-", terminator: "")
case _ where coeff < 0 && i < deg:
print(" - \(-coeff)", terminator: "")
case _ where i < deg:
print(" + \(coeff)", terminator: "")
case _:
print("\(coeff)", terminator: "")
}
if i > 1 {
print("x^\(i)", terminator: "")
} else if i == 1 {
print("x", terminator: "")
}
}
print()
}
let n = [-42, 0, -12, 1]
let d = [-3, 1, 0, 0]
print("Numerator: ", terminator: "")
polyPrint(n)
print("Denominator: ", terminator: "")
polyPrint(d)
guard let sol = polyLongDiv(numerator: n, denominator: d) else {
fatalError()
}
print("----------")
print("Quotient: ", terminator: "")
polyPrint(sol.quotient)
print("Remainder: ", terminator: "")
polyPrint(sol.remainder)</syntaxhighlight>
{{out}}
<pre>Numerator: x^3 - 12x^2 - 42
Denominator: x - 3
----------
Quotient: x^2 - 9x - 27
Remainder: -123</pre>
=={{header|Tcl}}==
{{works with|Tcl|8.5 and later}}
<
# Result is a list of two polynomials, q and r, where n = qd + r
# and the degree of r is less than the degree of b.
Line 522 ⟶ 3,872:
lassign [poldiv {-42. 0. -12. 1.} {-3. 1. 0. 0.}] Q R
puts [list Q = $Q]
puts [list R = $R]</
=={{header|Ursala}}==
The input is a pair of lists of coefficients in order of increasing degree.
Trailing zeros can be omitted. The output is a pair of lists (q,r), the quotient
and remainder polynomial coefficients. This is a straightforward implementation
of the algorithm in terms of list operations (fold, zip, map, distribute, etc.) instead
of array indexing, hence not unnecessarily verbose.
<syntaxhighlight lang="ursala">#import std
#import flo
polydiv =
zeroid~-l~~; leql?rlX\~&NlX ^H\(@rNrNSPXlHDlS |\ :/0.) @NlX //=> ?(
@lrrPX ==!| zipp0.; @x not zeroid+ ==@h->hr ~&t,
(^lryPX/~&lrrl2C minus^*p/~&rrr times*lrlPD)^/div@bzPrrPlXO ~&,
@r ^|\~& ~&i&& :/0.)</syntaxhighlight>
test program:
<syntaxhighlight lang="ursala">#cast %eLW
example = polydiv(<-42.,0.,-12.,1.>,<-3.,1.,0.,0.>)</syntaxhighlight>
output:
<pre>(
<-2.700000e+01,-9.000000e+00,1.000000e+00>,
<-1.230000e+02>)</pre>
=={{header|VBA}}==
{{trans|Phix}}<syntaxhighlight lang="vb">Option Base 1
Function degree(p As Variant)
For i = UBound(p) To 1 Step -1
If p(i) <> 0 Then
degree = i
Exit Function
End If
Next i
degree = -1
End Function
Function poly_div(ByVal n As Variant, ByVal d As Variant) As Variant
If UBound(d) < UBound(n) Then
ReDim Preserve d(UBound(n))
End If
Dim dn As Integer: dn = degree(n)
Dim dd As Integer: dd = degree(d)
If dd < 0 Then
poly_div = CVErr(xlErrDiv0)
Exit Function
End If
Dim quot() As Integer
ReDim quot(dn)
Do While dn >= dd
Dim k As Integer: k = dn - dd
Dim qk As Integer: qk = n(dn) / d(dd)
quot(k + 1) = qk
Dim d2() As Variant
d2 = d
ReDim Preserve d2(UBound(d) - k)
For i = 1 To UBound(d2)
n(UBound(n) + 1 - i) = n(UBound(n) + 1 - i) - d2(UBound(d2) + 1 - i) * qk
Next i
dn = degree(n)
Loop
poly_div = Array(quot, n) '-- (n is now the remainder)
End Function
Function poly(si As Variant) As String
'-- display helper
Dim r As String
For t = UBound(si) To 1 Step -1
Dim sit As Integer: sit = si(t)
If sit <> 0 Then
If sit = 1 And t > 1 Then
r = r & IIf(r = "", "", " + ")
Else
If sit = -1 And t > 1 Then
r = r & IIf(r = "", "-", " - ")
Else
If r <> "" Then
r = r & IIf(sit < 0, " - ", " + ")
sit = Abs(sit)
End If
r = r & CStr(sit)
End If
End If
r = r & IIf(t > 1, "x" & IIf(t > 2, t - 1, ""), "")
End If
Next t
If r = "" Then r = "0"
poly = r
End Function
Function polyn(s As Variant) As String
Dim t() As String
ReDim t(2 * UBound(s))
For i = 1 To 2 * UBound(s) Step 2
t(i) = poly(s((i + 1) / 2))
Next i
t(1) = String$(45 - Len(t(1)) - Len(t(3)), " ") & t(1)
t(2) = "/"
t(4) = "="
t(6) = "rem"
polyn = Join(t, " ")
End Function
Public Sub main()
Dim tests(7) As Variant
tests(1) = Array(Array(-42, 0, -12, 1), Array(-3, 1))
tests(2) = Array(Array(-3, 1), Array(-42, 0, -12, 1))
tests(3) = Array(Array(-42, 0, -12, 1), Array(-3, 1, 1))
tests(4) = Array(Array(2, 3, 1), Array(1, 1))
tests(5) = Array(Array(3, 5, 6, -4, 1), Array(1, 2, 1))
tests(6) = Array(Array(3, 0, 7, 0, 0, 0, 0, 0, 3, 0, 0, 1), Array(1, 0, 0, 5, 0, 0, 0, 1))
tests(7) = Array(Array(-56, 87, -94, -55, 22, -7), Array(2, 0, 1))
Dim num As Variant, denom As Variant, quot As Variant, rmdr As Variant
For i = 1 To 7
num = tests(i)(1)
denom = tests(i)(2)
tmp = poly_div(num, denom)
quot = tmp(1)
rmdr = tmp(2)
Debug.Print polyn(Array(num, denom, quot, rmdr))
Next i
End Sub</syntaxhighlight>{{out}}
<pre> x3 - 12x2 - 42 / x - 3 = x2 - 9x - 27 rem -123
x - 3 / x3 - 12x2 - 42 = 0 rem x - 3
x3 - 12x2 - 42 / x2 + x - 3 = x - 13 rem 16x - 81
x2 + 3x + 2 / x + 1 = x + 2 rem 0
x4 - 4x3 + 6x2 + 5x + 3 / x2 + 2x + 1 = x2 - 6x + 17 rem -23x - 14
x11 + 3x8 + 7x2 + 3 / x7 + 5x3 + 1 = x4 + 3x - 5 rem -16x4 + 25x3 + 7x2 - 3x + 8
-7x5 + 22x4 - 55x3 - 94x2 + 87x - 56 / x2 + 2 = -7x3 + 22x2 - 41x - 138 rem 169x + 220 </pre>
=={{header|Wren}}==
{{trans|Kotlin}}
===Version 1===
{{libheader|Wren-dynamic}}
<syntaxhighlight lang="wren">import "./dynamic" for Tuple
var Solution = Tuple.create("Solution", ["quotient", "remainder"])
var polyDegree = Fn.new { |p|
for (i in p.count-1..0) if (p[i] != 0) return i
return -2.pow(31)
}
var polyShiftRight = Fn.new { |p, places|
if (places <= 0) return p
var pd = polyDegree.call(p)
if (pd + places >= p.count) {
Fiber.abort("The number of places to be shifted is too large.")
}
var d = p.toList
for (i in pd..0) {
d[i + places] = d[i]
d[i] = 0
}
return d
}
var polyMultiply = Fn.new { |p, m|
for (i in 0...p.count) p[i] = p[i] * m
}
var polySubtract = Fn.new { |p, s|
for (i in 0...p.count) p[i] = p[i] - s[i]
}
var polyLongDiv = Fn.new { |n, d|
if (n.count != d.count) {
Fiber.abort("Numerator and denominator vectors must have the same size")
}
var nd = polyDegree.call(n)
var dd = polyDegree.call(d)
if (dd < 0) {
Fiber.abort("Divisor must have at least one one-zero coefficient")
}
if (nd < dd) {
Fiber.abort("The degree of the divisor cannot exceed that of the numerator")
}
var n2 = n.toList
var q = List.filled(n.count, 0)
while (nd >= dd) {
var d2 = polyShiftRight.call(d, nd - dd)
q[nd - dd] = n2[nd] / d2[nd]
polyMultiply.call(d2, q[nd - dd])
polySubtract.call(n2, d2)
nd = polyDegree.call(n2)
}
return Solution.new(q, n2)
}
var polyShow = Fn.new { |p|
var pd = polyDegree.call(p)
for (i in pd..0) {
var coeff = p[i]
if (coeff != 0) {
System.write(
(coeff == 1) ? ((i < pd) ? " + " : "") :
(coeff == -1) ? ((i < pd) ? " - " : "-") :
(coeff < 0) ? ((i < pd) ? " - %(-coeff)" : "%(coeff)") :
((i < pd) ? " + %( coeff)" : "%(coeff)")
)
if (i > 1) {
System.write("x^%(i)")
} else if (i == 1) {
System.write("x")
}
}
}
System.print()
}
var n = [-42, 0, -12, 1]
var d = [ -3, 1, 0, 0]
System.write("Numerator : ")
polyShow.call(n)
System.write("Denominator : ")
polyShow.call(d)
System.print("-------------------------------------")
var sol = polyLongDiv.call(n, d)
System.write("Quotient : ")
polyShow.call(sol.quotient)
System.write("Remainder : ")
polyShow.call(sol.remainder)</syntaxhighlight>
{{out}}
<pre>
Numerator : x^3 - 12x^2 - 42
Denominator : x - 3
-------------------------------------
Quotient : x^2 - 9x - 27
Remainder : -123
</pre>
===Version 2===
<syntaxhighlight lang="wren">class Polynom {
construct new(factors) {
_factors = factors.toList
}
factors { _factors.toList }
/(divisor) {
var curr = canonical().factors
var right = divisor.canonical().factors
var result = []
var base = curr.count - right.count
while (base >= 0) {
var res = curr[-1] / right[-1]
result.add(res)
curr = curr[0...-1]
for (i in 0...right.count-1) {
curr[base + i] = curr[base + i] - res * right[i]
}
base = base - 1
}
var quot = Polynom.new(result[-1..0])
var rem = Polynom.new(curr).canonical()
return [quot, rem]
}
canonical() {
if (_factors[-1] != 0) return this
var newLen = factors.count
while (newLen > 0) {
if (_factors[newLen-1] != 0) return Polynom.new(_factors[0...newLen])
newLen = newLen - 1
}
return Polynom.new(_factors[0..0])
}
toString { "Polynomial(%(_factors.join(", ")))" }
}
var num = Polynom.new([-42, 0, -12, 1])
var den = Polynom.new([-3, 1, 0, 0])
var res = num / den
var quot = res[0]
var rem = res[1]
System.print("%(num) / %(den) = %(quot) remainder %(rem)")</syntaxhighlight>
{{out}}
<pre>
Polynomial(-42, 0, -12, 1) / Polynomial(-3, 1, 0, 0) = Polynomial(-27, -9, 1) remainder Polynomial(-123)
</pre>
=={{header|zkl}}==
<syntaxhighlight lang="zkl">fcn polyLongDivision(a,b){ // (a0 + a1x + a2x^2 + a3x^3 ...)
_assert_(degree(b)>=0,"degree(%s) < 0".fmt(b));
q:=List.createLong(a.len(),0.0);
while((ad:=degree(a)) >= (bd:=degree(b))){
z,d,m := ad-bd, List.createLong(z,0.0).extend(b), a[ad]/b[bd];;
q[z]=m;
d,a = d.apply('*(m)), a.zipWith('-,d);
}
return(q,a); // may have trailing zero elements
}
fcn degree(v){ // -1,0,..len(v)-1, -1 if v==0
v.len() - v.copy().reverse().filter1n('!=(0)) - 1;
}
fcn polyString(terms){ // (a0,a1,a2...)-->"a0 + a1x + a2x^2 ..."
str:=[0..].zipWith('wrap(n,a){ if(a) "+ %sx^%s ".fmt(a,n) else "" },terms)
.pump(String)
.replace("x^0 "," ").replace(" 1x"," x").replace("x^1 ","x ")
.replace("+ -", "- ");
if(not str) return(" "); // all zeros
if(str[0]=="+") str[1,*]; // leave leading space
else String("-",str[2,*]);
}</syntaxhighlight>
<syntaxhighlight lang="zkl">q,r:=polyLongDivision(T(-42.0, 0.0, -12.0, 1.0),T(-3.0, 1.0));
println("Quotient = ",polyString(q));
println("Remainder = ",polyString(r));</syntaxhighlight>
{{out}}
<pre>
Quotient = -27 - 9x + x^2
Remainder = -123
</pre>
| |||