Polynomial long division: Difference between revisions

Content added Content deleted
(Updated D entry)
Line 511: Line 511:
<lang d>import std.stdio, std.range, std.algorithm, std.typecons;
<lang d>import std.stdio, std.range, std.algorithm, std.typecons;


Tuple!(double[],double[]) polyDiv(double[] N, in double[] D) {
Tuple!(double[],double[]) polyDiv(in double[] inN, in double[] inD)
/*pure nothrow*/ {
static nothrow pure int degree(T)(auto ref T[] poly) {
// code smell: a function that does two things
while (poly.length && poly[$-1] == 0)
poly = poly[0 .. $-1];
static int trimAndDegree(T)(ref T[] poly) /*nothrow pure*/ {
return cast(int)poly.length - 1;
poly.length -= poly.retro().countUntil!q{a != 0}();
return (cast(int)poly.length) - 1;
}
}


const dD = degree(D);
double[] N = inN.dup;
auto dN = degree(N);
const(double)[] D = inD;
const dD = trimAndDegree(D);
auto dN = trimAndDegree(N);
double[] q, r;
double[] q, r;
if (dD < 0)
if (dD < 0)
throw new Exception("ZeroDivisionError");
throw new Exception("ZeroDivisionError");
if (dN >= dD) {
if (dN >= dD) {
q = array(take(repeat(0.0), dN));
q = repeat(0.0).take(dN).array();
while (dN >= dD) {
while (dN >= dD) {
auto d = array(take(repeat(0.0), dN - dD)) ~ D;
auto d = repeat(0.0).take(dN - dD).array() ~ D;
auto mult = q[dN - dD] = N[$-1] / d[$-1];
const mult = q[dN - dD] = N[$ - 1] / d[$ - 1];
d[] *= mult;
d[] *= mult;
N[] -= d[];
N[] -= d[];
dN = degree(N);
dN = trimAndDegree(N);
}
}
r = N;
} else {
} else {
q = [0.0];
q = [0.0];
r = N;
}
}
r = N;
return tuple(q, r);
return tuple(q, r);
}
}


void main() {
void main() {
auto N = [-42.0, 0.0, -12.0, 1.0];
immutable N = [-42.0, 0.0, -12.0, 1.0];
auto D = [-3.0, 1.0, 0.0, 0.0];
immutable D = [-3.0, 1.0, 0.0, 0.0];
const qr = polyDiv(N, D);
writefln("%s / %s = %s remainder %s", N, D, polyDiv(N,D).tupleof);
writefln("%s / %s = %s remainder %s", N, D, qr[0], qr[1]);
}</lang>
}</lang>
{{out}}
Output:
<pre>[-42, 0, -12, 1] / [-3, 1, 0, 0] = [-27, -9, 1] remainder [-123]</pre>
<pre>[-42, 0, -12, 1] / [-3, 1, 0, 0] = [-27, -9, 1] remainder [-123]</pre>