Infinite lazy polynomials

gistfile1.js

/* 
 * Some powerseries fun based on:
 * http://jliszka.github.io/2013/10/31/infinite-lazy-polynomials.html 
 * */

if (!_.sum) {
    _.sum = function(list) { 
        return _.reduce(list, function(memo, val) { 
            return memo + val;
        }, 0); 
    };
};
if (!_.product) {
    _.product = function(list) { 
        return _.reduce(list, function(memo, val) { 
            return memo * val;
        }, 1); 
    };
};

var Poly = function(cf) {
    var poly = {};

    poly.getc = _.memoize(cf);

    poly.repr = function() { return _.map(_.range(10), poly.getc); };

    // poly + poly
    poly.add = function(other) {
        return Poly(function(n) { return poly.getc(n) + other.getc(n); });
    };
    poly.plus = poly.add;

    // poly - poly
    poly.subtract = function(other) {
        return Poly(function(n) { return poly.getc(n) - other.getc(n); });
    };
    poly.take = poly.subtract;
    poly.minus = poly.subtract;

    // -poly
    poly.unary_minus = function() {
        return Poly(function(n) { return -poly.getc(n); });
    };

    // poly * scalar
    poly.scalar_multiply = function(scalar) {
        return Poly(function(n) { return poly.getc(n) * scalar; });
    };
    poly.scalar_times = poly.scalar_multiply;

    // poly / scalar
    poly.scalar_divide = function(scalar) {
        return Poly(function(n) { return poly.getc(n) / scalar; });
    };
    
    // poly * poly
    poly.multiply = function(other) {
        return Poly(function(n) {
            return _.sum(_.map(_.range(n+1), function(i) {
                return poly.getc(i) * other.getc(n-i);
            }));
        });
    };
    poly.times = poly.multiply;

    // poly ^ scalar
    poly.power = _.memoize(function(scalar) {
        if (scalar == 0) { 
            return Poly.one 
        } else if (scalar % 1 == 0) {
            var s2 = poly.power(Math.floor(scalar / 2));    
            if (scalar % 2 == 0) {
                return s2.times(s2);
            } else {
                return poly.times(s2).times(s2);
            };
        } else {
            var a = poly.getc(0);
            var ar = Math.pow(a, scalar);
            var q = poly.scalar_divide(a).minus(Poly.one);
            var coeff = function(n) {
                return _.product(_.map(_.range(n), function(i) {
                            return scalar - i; 
                       })) / _.product(_.range(1, n+1));
            };
            return Poly(function(n) {
                return _.sum(_.map(_.range(n+1), function(i) {
                    return coeff(i) * q.power(i).getc(n);
                })) * ar;
            });
        };
    });
    poly.pow = poly.power;

    // 1 / poly
    poly.inverse = function() {
        var a = poly.getc(0);
        var q = Poly.one.minus(poly.scalar_divide(a));
        return Poly(function(n) {
            return _.sum(_.map(_.range(n+1), function(i) {
                return q.pow(i).getc(n);
            })) / a;
        });
    };
    poly.inv = poly.inverse;

    // poly / poly
    poly.divide = function(other) {
        return poly.times(other.inverse());
    };

    // e ^ poly
    poly.exp = function() {
        var a = poly.getc(0);
        var q = poly.minus(Poly.constant(a));
        var ea = Math.exp(a);
        var fact = function(n) { return _.product(_.range(1,n+1))};
        return Poly(function(n) {
            return _.sum(_.map(_.range(n+1), function(i) {
                return q.pow(i).getc(n) / fact(i);
            })) * ea;
        });
    };

    // ln(poly)
    poly.log = function() {
        var a = poly.getc(0);
        var logA = Poly.constant(Math.log(a));
        var q = poly.scalar_divide(a).minus(Poly.one);
        var alternatingHarmonic = function(n) {
            return Poly.constant(n % 2 == 0 ? -1.0 / n : 1);
        };
        return logA.plus(Poly(function(n) {
            return _.sum(_.map(_.range(1, n+1), function(i) {
                return alternatingHarmonic(i).times(q.pow(i)).getc(n)
            }));
        }));
    };
    poly.ln = poly.log;

    return poly;
};

Poly.one = Poly(function(n) {
    return n==0 ? 1 : 0;
});
Poly.x = Poly(function(n) {
    return n==1 ? 1 : 0;
});
Poly.constant = function(n) {
    return Poly.one.scalar_multiply(n)
};