Comments (4)
Thanks! I'll have a look when I get a chance.
from polynomials.jl.
Regarding my point two above, here's a proof-of-concept implementation for composing a Polynomial
with a degree one ImmutablePolynomial
:
function weighted_horner(x, p, f)
deg = length(p) - 1
w = f(Val(:init), deg)
out = (p[end] * w) * one(x)
for i ∈ (deg - 1):-1:0
w *= f(Val(:mul), i)
w = div(w, f(Val(:div), i), RoundToZero)
out = muladd(out, x, p[begin + i] * w)
end
out
end
struct PolynomialCompositionHelper
k::Int
end
bin_coef_bottom(h::PolynomialCompositionHelper) = h.k
function (h::PolynomialCompositionHelper)(::Val{:init}, deg::Int)
k = bin_coef_bottom(h)
binomial(deg + k, k)
end
(::PolynomialCompositionHelper)(::Val{:mul}, deg::Int) = deg + 1
(h::PolynomialCompositionHelper)(::Val{:div}, deg::Int) =
deg + 1 + bin_coef_bottom(h)
function polynomial_composed_with_shifted_identity(
p::Vector{T},
shift::T,
) where {T}
len = length(p)
out = Vector{T}(undef, len)
for k ∈ 0:(len - 1)
h = PolynomialCompositionHelper(k)
out[begin + k] = weighted_horner(shift, (@view p[(begin + k):end]), h)
end
out
end
function compose_polynomial_with_linear!(p::Vector{T}, scale::T) where {T}
s = scale
for i ∈ 1:(length(p) - 1)
p[begin + i] *= s
s *= scale
end
p
end
# Composition `p ∘ r`. `p` and `r` are the coefficients of polynomials
# in the standard basis.
composed_polynomials(p::Vector{T}, r::NTuple{2,T}) where {T} =
# Why this works:
#
# 1. `r` can be decomposed like so: `r = s ∘ l`, where
# `r(x) = a + b*x`, `s(x) = a + x`, `l(x) = b*x`
#
# 2. `p ∘ r = p ∘ (s ∘ l) = (p ∘ s) ∘ l`
compose_polynomial_with_linear!(
polynomial_composed_with_shifted_identity(p, first(r)),
last(r),
)
using Polynomials
(p::Polynomial{T, v})(r::ImmutablePolynomial{T, v, 2}) where {T, v} =
Polynomial{T, v}(composed_polynomials(coeffs(p), coeffs(r)))
I didn't compare accuracy yet, however the results appear to be approximately correct:
julia> include("polynomial_composition.jl") # includes the above code
julia> p = Polynomial(rand(6));
julia> q = ImmutablePolynomial((rand(), rand()))
ImmutablePolynomial(0.7923048129450314 + 0.4495999592340951*x)
julia> p(q) # new implementation
Polynomial(2.2449893995110872 + 2.4472444477738295*x + 1.902292181529416*x^2 + 0.8098154846551124*x^3 + 0.18248093200464596*x^4 + 0.017903450489708504*x^5)
julia> p(Polynomial(q)) # old implementation
Polynomial(2.244989399511087 + 2.447244447773829*x + 1.9022921815294158*x^2 + 0.8098154846551123*x^3 + 0.18248093200464593*x^4 + 0.017903450489708504*x^5)
The above implementation is faster than what's currently available in Polynomials.jl, however I'll wait until your upcoming PR is merged before making concrete comparisons.
from polynomials.jl.
Thanks! I'll have to look this over to see how to phase it in. It seems like there should be some big performance optimizations when the immutable polynomial type is involved.
from polynomials.jl.
After playing some more with my implementation above, I realized that the approach is possibly problematic: in weighted_horner
, the weight w
is mathematically an integer, however it's huge for large-degree polynomials, and f(Val(:init), deg)
overflows Int
starting when deg
is somewhere between 60 and 70.
So I guess this should be adapted slightly and restricted to something like Polynomial{<:Union{AbstractFloat,BigInt}}
.
from polynomials.jl.
Related Issues (20)
- docs: advertise composition of polynomials HOT 1
- `ImmutablePolynomial` constructors seem to prevent good type inference HOT 4
- bad type inference for composition of `ImmutablePolynomial`s HOT 2
- Should `AbstractPolynomial` subtype `Function`? HOT 2
- Question: is multipliciation of polynomials implemented with FFT HOT 6
- type stability: `(::Polynomial)(::ImmutablePolynomial)` composition yields `ImmutablePolynomial` HOT 2
- polynomials should have `big` methods? HOT 2
- Should `coeffs(::ImmutablePolynomial)` be type-stable? HOT 1
- Feature request: decompose a polynomial HOT 1
- Polynomials.fit fails with InexactError for Rational HOT 2
- Question about using HOT 2
- Warnings and Test failure with Julia 1.10 HOT 2
- Polynomial/Polynomial = Rational or? HOT 1
- (Low priority) Axis offset warning when using OffsetArrays HOT 3
- Conversions to and from RationalFunction broken HOT 3
- MakieCore doesn't need to be in `deps` HOT 12
- Issues with differentiating with respect to polynomial coefficients
- Polynomials.fit changes behavior somewhere between (inclusive) Julia 1.7.3 and Julia 1.10.4, the latter being self-inconsistent (problematic)
- Polynomials requiring MakieCore HOT 5
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from polynomials.jl.