Polynomial long division: Difference between revisions
Undo revision 361341 by Peak (talk)
(add FreeBASIC) |
|||
| (24 intermediate revisions by 13 users not shown) | |||
Line 28:
* 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 87 ⟶ 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:
<
procedure Long_Division is
Line 199 ⟶ 245:
Put ("Q: "); Output (Test_Q);
Put ("R: "); Output (Test_R);
end Long_Division;</
output:
Line 208 ⟶ 254:
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}}==
<
{
q r d←⍵
Line 218 ⟶ 381:
} ⍬ ⍺ ⍵
}
</syntaxhighlight>
{{out}}
<pre> N←¯42 0 ¯12 1
Line 228 ⟶ 391:
=={{header|BBC BASIC}}==
<
DIM D%(3) : D%() = -3, 1, 0, 0
DIM q%(3), r%(3)
Line 271 ⟶ 434:
IF n%<0 THEN = " - " + STR$(-n%)
IF n%=1 THEN = " + "
= " + " + STR$(n%)</
'''Output:'''
<pre>
Line 282 ⟶ 445:
{{libheader|GNU Scientific Library}}
<
#include <stdlib.h>
#include <stdarg.h>
Line 388 ⟶ 551:
return r;
}</
<
{
int i;
Line 412 ⟶ 575:
return 0;
}</
===Another version===
Without outside libs, for clarity. Note that polys are stored and show with zero-degree term first:<
#include <stdlib.h>
#include <stdarg.h>
Line 518 ⟶ 681:
return 0;
}</
=={{header|C sharp|C#}}==
{{trans|Java}}
<
namespace PolynomialLongDivision {
Line 644 ⟶ 807:
}
}
}</
{{out}}
<pre>Numerator : x^3 - 12.0x^2 - 42.0
Line 653 ⟶ 816:
=={{header|C++}}==
<
#include <iostream>
#include <iterator>
Line 752 ⟶ 915:
}
</syntaxhighlight>
=={{header|Clojure}}==
Line 760 ⟶ 923:
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.
<
(let [grade1 (reduce +' term1)
grade2 (reduce +' term2)
Line 854 ⟶ 1,017:
(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} {}))))</
=={{header|Common Lisp}}==
Line 860 ⟶ 1,023:
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.
<
(do ((sum '())) ((and (endp p1) (endp p2)) (nreverse sum))
(let ((pd1 (if (endp p1) -1 (caar p1)))
Line 900 ⟶ 1,063:
(if (endp quotient) (return (values sum remainder))
(setf dividend remainder
sum (add quotient sum)))))))</
The [[wp:Polynomial_long_division#Example|wikipedia example]]:
<
'((1 . 1) (0 . -3))) ; x - 3
((2 . 1) (1 . -9) (0 . -27)) ; x^2 - 9x - 27
((0 . -123)) ; -123</
=={{header|D}}==
<
Tuple!(double[], double[]) polyDiv(in double[] inN, in double[] inD)
Line 951 ⟶ 1,114:
immutable D = [-3.0, 1.0, 0.0, 0.0];
writefln("%s / %s = %s remainder %s", N, D, polyDiv(N, D)[]);
}</
{{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}}==
Line 1,056 ⟶ 1,466:
}</pre>
<
def d := makePolynomial([-3, 1])
println("Numerator: ", n)
Line 1,062 ⟶ 1,472:
def [q, r] := n.quotRem(d, stdout)
println("Quotient: ", q)
println("Remainder: ", r)</
Output:
Line 1,082 ⟶ 1,492:
=={{header|Elixir}}==
{{trans|Ruby}}
<
def division(_, []), do: raise ArgumentError, "denominator is zero"
def division(_, [0]), do: raise ArgumentError, "denominator is zero"
Line 1,112 ⟶ 1,522:
{q, r} = Polynomial.division(f, g)
IO.puts "#{inspect f} / #{inspect g} => #{inspect q} remainder #{inspect r}"
end)</
{{out}}
Line 1,122 ⟶ 1,532:
</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}}==
<
{ -42 0 -12 1 } { -3 1 } p/mod ptrim [ . ] bi@</
{{out}}
<pre>
Line 1,135 ⟶ 1,595:
{{works with|Fortran|95 and later}}
<
implicit none
Line 1,214 ⟶ 1,674:
end subroutine poly_print
end module Polynom</
<
use Polynom
implicit none
Line 1,233 ⟶ 1,693:
deallocate(q, r)
end program PolyDivTest</
=={{Header|FreeBASIC}}==
<
type polyterm
Line 1,251 ⟶ 1,711:
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)))
Line 1,335 ⟶ 1,795:
poly_print Q() 'quotient
poly_print R() 'remainder</
{{out}}
<pre>x^2 - 9x - 27
Line 1,342 ⟶ 1,802:
=={{header|GAP}}==
GAP has built-in functions for computations with polynomials.
<
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 ]</
=={{header|Go}}==
By the convention and pseudocode given in the task:
<
import "fmt"
Line 1,395 ⟶ 1,855:
}
return q, nn, true
}</
Output:
<pre>
Line 1,408 ⟶ 1,868:
Translated from the OCaml code elsewhere on the page.
{{works with|GHC|6.10}}
<
shift n l = l ++ replicate n 0
Line 1,429 ⟶ 1,889:
ks = map (* k) $ shift ddif s
q' = zipWith' (+) q $ shift ddif [k]
f' = norm $ tail $ zipWith' (-) f ks</
And this is the also-translated pretty printing function.
<
where term v 0 = show v
term 1 1 = "x"
Line 1,443 ⟶ 1,903:
terms [] = []
terms (0:t) = terms t
terms (h:t) = (term h (length t)) : (terms t)</
=={{header|J}}==
Line 1,449 ⟶ 1,909:
From http://www.jsoftware.com/jwiki/Phrases/Polynomials
<
Wikipedia example:
<
This produces the result:
┌────────┬────┐
Line 1,466 ⟶ 1,926:
<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;
Line 2,010 ⟶ 2,470:
}
}
</syntaxhighlight>
{{out}}
<pre>
Line 2,024 ⟶ 2,484:
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
Line 2,037 ⟶ 2,550:
println(p, " divided by ", q, " is ", d, " with remainder ", r, ".")
</syntaxhighlight>
{{out}}
Line 2,045 ⟶ 2,558:
=={{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) {
Line 2,057 ⟶ 2,571:
return Int.MIN_VALUE
}
fun polyShiftRight(p: DoubleArray, places: Int): DoubleArray {
if (places <= 0) return p
Line 2,071 ⟶ 2,585:
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) {
Line 2,103 ⟶ 2,617:
return Solution(q, n2)
}
fun polyShow(p: DoubleArray) {
val pd = polyDegree(p)
Line 2,120 ⟶ 2,634:
println()
}
fun main(args: Array<String>) {
val n = doubleArrayOf(-42.0, 0.0, -12.0, 1.0)
Line 2,134 ⟶ 2,648:
print("Remainder : ")
polyShow(r)
}</
{{out}}
<pre>
Output:
Numerator : x^3 - 12.0x^2 - 42.0
Denominator : x - 3.0
Line 2,143 ⟶ 2,659:
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
Line 2,167 ⟶ 2,730:
> expand( (x^2+2)*q + r - p ); # check
0
</syntaxhighlight>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<
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).
<
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
Line 2,184 ⟶ 2,843:
let zip op p q =
let d = (List.length p) - (List.length q) in
List.map2 op (pad (-d) p) (pad d q)</
Then the main polynomial division function
<
let rec aux f s q =
let ddif = (deg f) - (deg s) in
Line 2,195 ⟶ 2,854:
and f' = norm (List.tl (zip (-.) f ks)) in
aux f' s q' in
aux (norm f) (norm g) []</
For output we need a pretty-printing function
<
let term v p = match (v, p) with
| ( _, 0) -> string_of_float v
Line 2,208 ⟶ 2,867:
| h :: t ->
if h = 0.0 then (terms t) else (term h (List.length t)) :: (terms t) in
String.concat " + " (terms l)</
and then the example
<
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"
(str_poly f) (str_poly g) (str_poly q) (str_poly r)</
gives the output:
<pre>
Line 2,227 ⟶ 2,886:
Octave has already facilities to divide two polynomials (<code>deconv(n,d)</code>); and the reason to adopt the convention of keeping the highest power coefficient first, is to make the code compatible with builtin functions: we can use <tt>polyout</tt> to output the result.
<
gd = length(d);
pv = zeros(1, length(n));
Line 2,263 ⟶ 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.
<
my(rem=a%b);
[(a - rem)/b, rem]
};
poldiv(x^9+1, x^3+x-3)</
Alternately, use the built-in function <code>divrem</code>:
<
=={{header|Perl}}==
Line 2,279 ⟶ 2,938:
{{trans|Octave}}
<
use List::Util qw(min);
Line 2,300 ⟶ 2,959:
return ( [0], $rn );
}
}</
<
{
my @c = @_;
Line 2,311 ⟶ 2,970:
}
print "\n";
}</
<
($q, $r) = poly_long_div([1, -12, 0, -42], [1, -3]);
Line 2,330 ⟶ 2,989:
($q, $r) = poly_long_div([1,-4,6,5,3], [1,2,1]);
poly_print(@$q);
poly_print(@$r);</
=={{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
Line 2,340 ⟶ 3,003:
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
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
rem[mi] -= d2[mi]*qk
end for
dn = degree(
end while
return {
end function
function poly(sequence si)
-- display helper
string r = ""
for t=length(si) to 1 by -1 do
Line 2,383 ⟶ 3,048:
return r
end function
constant tests = {{{-42,0,-12,1},{-3,1}},
{{-3,1},{-42,0,-12,1}},
Line 2,399 ⟶ 3,057:
{{-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,
{
printf(1,fmt,
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
Line 2,419 ⟶ 3,078:
=={{header|PicoLisp}}==
<
(let I NIL
(for (N . C) P
Line 2,436 ⟶ 3,095:
(set (nth Q (inc Diff)) F)
(setq N (mapcar '((N E) (- N (* E F))) N E)) ) ) )
(list Q N) ) ) )</
Output:
<pre>: (divPoly (-42 0 -12 1) (-3 1 0 0))
Line 2,443 ⟶ 3,102:
=={{header|Python}}==
{{works with|Python 2.x}}
<
from itertools import izip
def degree(poly):
Line 2,463 ⟶ 3,121:
mult = q[dN - dD] = N[-1] / float(d[-1])
d = [coeff*mult for coeff in d]
N = [
dN = degree(N)
r = N
Line 2,476 ⟶ 3,134:
D = [-3, 1, 0, 0]
print " %s / %s =" % (N,D),
print " %s remainder %s" % poly_div(N, D)</
Sample output:
Line 2,484 ⟶ 3,142:
=={{header|R}}==
{{trans|Octave}}
<
gd <- length(d)
pv <- vector("numeric", length(n))
Line 2,517 ⟶ 3,175:
r <- polylongdiv(c(1,-12,0,-42), c(1,-3))
print.polynomial(r$q)
print.polynomial(r$r)</
=={{header|Racket}}==
<
#lang racket
(define (deg p)
Line 2,556 ⟶ 3,214:
'#(-27 -9 1)
'#(-123 0)
</syntaxhighlight>
=={{header|Raku}}==
Line 2,562 ⟶ 3,220:
{{Works with|rakudo|2018.10}}
{{trans|Perl}} for the core algorithm; original code for LaTeX pretty-printing.
<syntaxhighlight lang="raku"
return [0], |@n if +@n < +@d;
Line 2,585 ⟶ 3,243:
printf Q"%s , & %s \\\\\n", poly_long_div( @a, @b ).map: { poly_print($_) };
}
say '\end{array}</math>';</
Output:
Line 2,597 ⟶ 3,255:
=={{header|REXX}}==
<
z='1 -12 0 -42' /* Numerator */
n='1 -3' /* Denominator */
Line 2,640 ⟶ 3,298:
d=d-1
End
Return strip(space(res,0),'L','+')</
{{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:
<
dd = degree(denominator)
raise ArgumentError, "denominator is zero" if dd < 0
Line 2,692 ⟶ 3,359:
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]</
Implementing the algorithms on the [[wp:Polynomial long division|wikipedia page]] -- uglier code but nicer user interface
<
if g.length == 0 or (g.length == 1 and g[0] == 0)
raise ArgumentError, "denominator is zero"
Line 2,762 ⟶ 3,429:
q, r = polynomial_division(f, g)
puts "#{f} / #{g} => #{q} remainder #{r}"
# => [1, -12, 0, -42] / [1, 1, -3] => [1, -13] remainder [16, -81]</
Best of both worlds: {{trans|Tcl}}
<
if g.length == 0 or (g.length == 1 and g[0] == 0)
raise ArgumentError, "denominator is zero"
Line 2,793 ⟶ 3,460:
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]</
=={{header|Sidef}}==
{{trans|Perl}}
<
var n = rn.map{_}
Line 2,816 ⟶ 3,483:
return([0], rn)
}</
Example:
<
var l = c.len
c.each_kv {|i, n|
Line 2,841 ⟶ 3,508:
poly_print(r)
print "\n"
}</
{{out}}
Line 2,860 ⟶ 3,527:
=={{header|Slate}}==
<
p@(Polynomial traits) new &capacity: n
Line 2,921 ⟶ 3,588:
{q. n}]
ifFalse: [{p newFrom: #(0). p copy}]
].</
=={{header|Smalltalk}}==
{{works with|GNU Smalltalk}}
<
|coeffs|
Polynomial class >> new [ ^ super basicNew init ]
Line 2,993 ⟶ 3,660:
]
displayNl [ self display. Character nl display ]
].</
<
res := OrderedCollection new.
Line 3,005 ⟶ 3,672:
res do: [ :o |
(o at: 1) display. ' with rest: ' display. (o at: 2) displayNl
]</
=={{header|SPAD}}==
Line 3,011 ⟶ 3,678:
{{works with|OpenAxiom}}
{{works with|Axiom}}
<
2
(1) [quotient = x - 9x - 27,remainder = - 123]
Type: Record(quotient: Polynomial(Integer),remainder: Polynomial(Integer))</
Domain:[http://fricas.github.io/api/PolynomialCategory.html#l-polynomial-category-monic-divide]
Line 3,024 ⟶ 3,691:
{{trans|Kotlin}}
<
static func / (lhs: Self, rhs: Self) -> Self
}
Line 3,153 ⟶ 3,820:
polyPrint(sol.quotient)
print("Remainder: ", terminator: "")
polyPrint(sol.remainder)</
{{out}}
Line 3,166 ⟶ 3,833:
{{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 3,205 ⟶ 3,872:
lassign [poldiv {-42. 0. -12. 1.} {-3. 1. 0. 0.}] Q R
puts [list Q = $Q]
puts [list R = $R]</
=={{header|Ursala}}==
Line 3,213 ⟶ 3,880:
of the algorithm in terms of list operations (fold, zip, map, distribute, etc.) instead
of array indexing, hence not unnecessarily verbose.
<
#import flo
Line 3,221 ⟶ 3,888:
@lrrPX ==!| zipp0.; @x not zeroid+ ==@h->hr ~&t,
(^lryPX/~&lrrl2C minus^*p/~&rrr times*lrlPD)^/div@bzPrrPlXO ~&,
@r ^|\~& ~&i&& :/0.)</
test program:
<
example = polydiv(<-42.,0.,-12.,1.>,<-3.,1.,0.,0.>)</
output:
<pre>(
Line 3,232 ⟶ 3,899:
=={{header|VBA}}==
{{trans|Phix}}<
Function degree(p As Variant)
For i = UBound(p) To 1 Step -1
Line 3,327 ⟶ 3,994:
Debug.Print polyn(Array(num, denom, quot, rmdr))
Next i
End Sub</
<pre> x3 - 12x2 - 42 / x - 3 = x2 - 9x - 27 rem -123
x - 3 / x3 - 12x2 - 42 = 0 rem x - 3
Line 3,335 ⟶ 4,002:
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}}==
<
_assert_(degree(b)>=0,"degree(%s) < 0".fmt(b));
q:=List.createLong(a.len(),0.0);
Line 3,358 ⟶ 4,179:
if(str[0]=="+") str[1,*]; // leave leading space
else String("-",str[2,*]);
}</
<
println("Quotient = ",polyString(q));
println("Remainder = ",polyString(r));</
{{out}}
<pre>
| |||