Решение на Polynomial от Михаил Младенов
Към профила на Михаил Младенов
Резултати
- 15 точки от тестове
- 2 бонус точки
- 17 точки общо
- 15 успешни тест(а)
- 0 неуспешни тест(а)
Код
Лог от изпълнението
Compiling solution v0.1.0 (file:///tmp/d20171121-6053-1gpavm8/solution)
Finished dev [unoptimized + debuginfo] target(s) in 5.56 secs
Running target/debug/deps/solution-200db9172ea1f728
running 0 tests
test result: ok. 0 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out
Running target/debug/deps/solution_test-e3c9eb714e09105e
running 15 tests
test solution_test::test_add_poly ... ok
test solution_test::test_add_poly_zero_one ... ok
test solution_test::test_arithmetic_properties ... ok
test solution_test::test_create_poly ... ok
test solution_test::test_div_poly_f64 ... ok
test solution_test::test_div_poly_f64_zero ... ok
test solution_test::test_fp_comparison ... ok
test solution_test::test_has_point ... ok
test solution_test::test_lagrange_poly_1 ... ok
test solution_test::test_lagrange_poly_2 ... ok
test solution_test::test_lagrange_poly_err_eq_x ... ok
test solution_test::test_mul_poly ... ok
test solution_test::test_mul_poly_f64 ... ok
test solution_test::test_mul_poly_f64_zero ... ok
test solution_test::test_mul_poly_zero_one ... ok
test result: ok. 15 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out
Doc-tests solution
running 0 tests
test result: ok. 0 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out
История (2 версии и 5 коментара)
Михаил качи решение на 13.11.2017 05:47 (преди почти 8 години)
use std::ops::Mul;
use std::ops::Div;
use std::ops::Add;
#[derive(Debug, Clone)]
pub struct Polynomial {
coefficients: Vec<f64>
}
const EPS: f64 = 1e-10;
impl Polynomial {
pub fn has(&self, point: &(f64, f64)) -> bool {
(self.evaluate_at(point.0) - point.1).abs() < EPS
}
pub fn interpolate(points: Vec<(f64, f64)>) -> Option<Self> {
if points.len() < 2 {
None
}
else
{
let mut result = Polynomial::default();
let mut l = Polynomial::from(vec![1.0]);
- result.coefficients.reserve_exact(points.len());
// Първо пресмятаме $l := \prod\limits_{k=0}^n (x - x_k)$, след това от
// него получаваме на всяка стъпка $l_j$ имайки в предвид равенството
// $l_j = \frac{\frac{l}{x - x_j}}{\frac{l}{x - x_j}(x_j)}$
for j in 0..points.len() {
l = l * Polynomial::from(vec![1.0, -points[j].0]);
}
for j in 0..points.len() {
let mut ll = l.clone();
let mut l_j = Polynomial::default();
// WTF?!? error: use of unstable library feature vec_resize_default
// l_j.coefficients.resize_default(l.coefficients.len() - 1);
// За съжаление трябва да мина на този гаден вариант :/
l_j.coefficients.resize(l.coefficients.len() - 1, 0.0);
let mut k = l.coefficients.len() - 1;
while 0 < k {
l_j.coefficients[k - 1] = ll.coefficients[k];
ll.coefficients[k - 1] += ll.coefficients[k] * points[j].0;
k -= 1;
}
let d = l_j.evaluate_at(points[j].0);
// Това значи че е имало повтарящи се точки
if d.abs() < EPS {
return None;
}
result = result + (l_j * points[j].1)/d;
}
result.trim();
Some(result)
}
}
}
impl From<Vec<f64>> for Polynomial {
fn from(mut coefs: Vec<f64>) -> Self {
if coefs.is_empty() {
Polynomial {
coefficients : vec![0.0]
}
}
else {
coefs.reverse();
let mut p = Polynomial {
coefficients : coefs
};
p.trim();
p
}
}
}
impl Default for Polynomial {
fn default() -> Self {
Polynomial {
coefficients : vec![0.0]
}
}
}
impl PartialEq for Polynomial {
fn eq(&self, rhs: &Self) -> bool {
if self.coefficients.len() != rhs.coefficients.len() {
false
}
else {
for i in 0..self.coefficients.len() {
if (self.coefficients[i] - rhs.coefficients[i]).abs() >= EPS {
return false;
}
}
true
}
}
}
impl Mul<f64> for Polynomial {
type Output = Polynomial;
// Защо self, a не mut self или &self?
fn mul(self, rhs: f64) -> Self::Output {
let mut result = Polynomial {
coefficients: Vec::with_capacity(self.coefficients.len())
};
for v in self.coefficients.iter() {
result.coefficients.push(v * rhs);
}
result.trim();
result
Вместо да инстанцираш Polynomial, да push-ваш коефициенти в него, и после да го trim-ваш, би могъл да си алокираш вектор, и после просто да викнеш Polynomial::from(...). Това по принцип ни беше идеята -- trim-ването на коефициентите да се направи при конструиране на полинома.
Не е голяма работа, просто нещо дребно, което може леко да съкрати кода и да ти даде малко повече reliability -- ако винаги конструираш вектор по един и същ начин, и винаги го нормализираш при конструиране, ще можеш смело да викаш from, без на всяко място, където изграждаш вектор, да мислиш дали си забравил нещо.
}
}
impl Div<f64> for Polynomial {
type Output = Polynomial;
fn div(self, rhs: f64) -> Self::Output {
let mut result = Polynomial {
coefficients: Vec::with_capacity(self.coefficients.len())
};
for v in self.coefficients.iter() {
result.coefficients.push(v / rhs);
}
result.trim();
result
}
}
const FFT_POLY_MUL_THRESHOLD: usize = 120;
impl Mul for Polynomial {
type Output = Polynomial;
fn mul(self, rhs: Polynomial) -> Self::Output {
let m = self.coefficients.len();
let n = rhs.coefficients.len();
if m == 1 {
rhs * self.coefficients[0]
}
else if n == 1 {
self * rhs.coefficients[0]
}
else if n < FFT_POLY_MUL_THRESHOLD || m < FFT_POLY_MUL_THRESHOLD {
// O(n^2) умножение
mult_polynomials_quadratic(&self, &rhs)
}
else {
// O(nlog(n)) умножение чрез бързата трансформации на Фуурие, това всъщност е
// оправдание за да мога да ползвам удобната стандартна функция
// fn next_power_of_two(self) -> usize :)
mult_polynomials_fft(&self, &rhs)
}
}
}
impl Add for Polynomial {
type Output = Polynomial;
// Тук интерфейса е self, нито mut self нито , &self което е странно.
На пръв поглед е странно, но идеята е + да взима ownership и да върне стойност. Има AddAssign (+=) което е 'in place' и взима &mut self.
Не ме притеснява, че взима ownership, a че след като го вземем няма как да го върнем (защото функцията приема два аргумента а връщаме само един полином) и освен това не правим нищо с тези обекти след като им вземем ownership (тоест ще ги оставим да умрат, защото ако се опитаме да пишем в тях компилатора се оплаква, че пипаме мютъбъл вериъбъл). Ако беше &self, &rhs старите обекти щяха да си останат. А ако беше mut self, mut rhs, поне след взимането на ownership, смъртта им нямаше да е напразна защото поне има как да ползваме алокираните за тях ресурси и да конструираме резултата в тях. Тоест нещо като sometype_t operator+(sometype_t&&, sometype_t&&) в C++ (Само като пример го давам там би било тъпо да се ползва такова нещо защото има много по мощен начин за избягване на копирания и временни обекти от такива аритметични изрази. Чрез expression templates можем в compile time да парснем дървото на израза и винаги да изберем оптималния начин за смятане избягващ максимален брой копирания).
Все пак не съм сигурен дали разбирам как точно работи това. Тоест има ли как да извлечем каквато и да е полза от това че тези обекти им взимаме ownership в тая функция, а те не са mutable?
Всъщност, можеш да дефинираш функцията като mut self, и това няма да гръмне. В случая, параметъра е еквивалентно на деклариране на променлива, така че просто self е все едно let self = ..., а mut self ще е let mut self = ..., но и двете работят със сигнатурата. Пробвай. Признавам, че може би трябваше да го напишем в условието така, или да го обясним на някоя лекция.
И да, това, че можеш да мутираш self, означава че, на теория, не е нужно да конструираш нов полином, а можеш да използваш десния. Но да, признавам, че и на мен ми е малко неинтуитивно.
fn add(self, rhs: Polynomial) -> Self::Output {
use ::std::cmp::max;
use ::std::cmp::min;
let mut result = Polynomial {
coefficients: Vec::with_capacity(max(self.coefficients.len(),rhs.coefficients.len()))
};
let m = min(self.coefficients.len(), rhs.coefficients.len());
for i in 0..m {
result.coefficients.push(self.coefficients[i] + rhs.coefficients[i]);
}
for i in m..self.coefficients.len() {
result.coefficients.push(self.coefficients[i]);
}
for i in m..rhs.coefficients.len() {
result.coefficients.push(rhs.coefficients[i]);
}
result.trim();
result
}
}
pub fn mult_polynomials_quadratic(a: &Polynomial, b: &Polynomial) -> Polynomial {
let n = a.coefficients.len();
let m = b.coefficients.len();
let d = n + m ;
let mut c = Polynomial { coefficients: Vec::with_capacity(d) };
for j in 1..d {
let mut c_j: f64 = 0.0;
for k in 0..j {
if k < n && j - 1 - k < m {
c_j += a.coefficients[k] * b.coefficients[j - 1 - k];
}
}
c.coefficients.push(c_j);
}
c.trim();
c
}
#[derive(Default, Clone, Copy, Debug, PartialEq)]
struct Complex {
real: f64,
imag: f64
}
impl From<(f64,f64)> for Complex {
fn from(p: (f64,f64)) -> Self{
Complex {
real: p.0,
imag: p.1
}
}
}
impl Mul for Complex {
type Output = Complex;
fn mul(self, rhs: Complex) -> Self::Output {
Complex {
real: self.real * rhs.real - self.imag * rhs.imag,
imag: self.imag * rhs.real + self.real * rhs.imag
}
}
}
impl Add for Complex {
type Output = Complex;
fn add(self, rhs: Complex) -> Self::Output {
Complex {
real: self.real + rhs.real,
imag: self.imag + rhs.imag
}
}
}
use std::ops::DivAssign;
impl DivAssign<f64> for Complex {
fn div_assign(&mut self, rhs: f64) {
self.real /= rhs;
self.imag /= rhs;
}
}
use std::ops::Sub;
impl Sub for Complex {
type Output = Complex;
fn sub(self, rhs: Complex) -> Self::Output {
Complex {
real: self.real - rhs.real,
imag: self.imag - rhs.imag
}
}
}
pub fn mult_polynomials_fft(a: &Polynomial, b: &Polynomial) -> Polynomial {
let mut d = a.coefficients.len() + b.coefficients.len();
d = d.next_power_of_two();
// Знаем, че коефициентите на произведението а полиноми са циклична конволюция на
// коефициентите на множителите. От друга страна дискретната трансформация на Фуурие
// от конволюцията е Адамаровото произведение на дискретните трансформации на Фуурие
// от множителите. Това означава, че ако вземем обратната трансформация на Фуурие от
// Адамаровото произведение на дискретните трансформации на Фуурие на двата множителя,
// ще получим точно произведението на полиномите.
let mut ac: Vec<Complex> = Vec::with_capacity(d);
let mut bc: Vec<Complex> = Vec::with_capacity(d);
for i in 0..d {
if i < a.coefficients.len() {
ac.push(Complex::from((a.coefficients[i], 0.0)));
}
else {
ac.push(Complex::from((0.0, 0.0)));;
}
}
for i in 0..d {
if i < b.coefficients.len() {
bc.push(Complex::from((b.coefficients[i], 0.0)));
}
else {
bc.push(Complex::from((0.0, 0.0)));
}
}
let mut ffta: Vec<Complex>;
let mut fftb: Vec<Complex>;
// Тъй като resize_default не работи ми се приисква да направя така:
unsafe {
let mut tmp1 = Vec::<Complex>::with_capacity(d);
let mut tmp2 = Vec::<Complex>::with_capacity(d);
ffta = Vec::from_raw_parts(tmp1.as_mut_ptr(), d, d);
fftb = Vec::from_raw_parts(tmp2.as_mut_ptr(), d, d);
::std::mem::forget(tmp1);
::std::mem::forget(tmp2);
}
// Надявам се да не ме биете много...
fft_power_of_two(&ac, &mut ffta, false);
fft_power_of_two(&bc, &mut fftb, false);
for i in 0..d {
ffta[i] = ffta[i] * fftb[i];
}
fft_power_of_two(&ffta, &mut fftb, true);
let mut result = Polynomial {
coefficients: Vec::with_capacity(d)
};
for i in 0..d {
result.coefficients.push(fftb[i].real);
}
result.trim();
result
}
fn fft_power_of_two(points: &Vec<Complex>, transforms: &mut Vec<Complex>, inverse: bool) {
let n = points.len();
let logn : usize = (n as f64).log(2.0).round() as usize;
let bits_in_usize = ::std::mem::size_of::<usize>() * 8;
// Разместваме точките, както биха били в дъното на рекурсията, ако използвахме
// рекурсивната версия на алгоритъма.
for i in 0..n {
transforms[i] = points[reverse_bits(i) >> (bits_in_usize - logn)];
}
// Сега започваме да комбинираме субтрансформациите чрез butterfly операции
use std::f64::consts::PI;
// Итерираме по височината на въображаемото дърво на рекурсията.
for height in 0..logn {
let subfft_size: usize = 1 << (height + 1);
let phase: f64 = (if inverse {-1.0} else {1.0}) * 2.0 * PI / subfft_size as f64;
let primitive_root_of_unity = Complex::from((phase.cos(), phase.sin()));
let mut subfft_offset: usize = 0;
while subfft_offset < n {
let mut twiddle_factor = Complex::from((1.0, 0.0));
let half_subfft = subfft_size/2;
// Сега крилатите гъсеници :)
for k in 0..half_subfft {
let e = transforms[subfft_offset + k];
let o = twiddle_factor * transforms[subfft_offset + k + half_subfft];
twiddle_factor = twiddle_factor * primitive_root_of_unity;
transforms[subfft_offset + k] = e + o;
transforms[subfft_offset + k + half_subfft] = e - o;
}
subfft_offset += subfft_size;
}
}
if inverse {
for i in 0..n {
transforms[i] /= n as f64;
}
}
}
fn reverse_bits(mut n: usize) -> usize {
let mut mask = 0b0101010101010101010101010101010101010101010101010101010101010101usize;
n = ((n & mask) << 1) | ((n >> 1) & mask);
mask = 0b0011001100110011001100110011001100110011001100110011001100110011usize;
n = ((n & mask) << 2) | ((n >> 2) & mask);
mask = 0x0f0f0f0f0f0f0f0fusize;
n = ((n & mask) << 4) | ((n >> 4) & mask);
mask = 0x00ff00ff00ff00ffusize;
n = ((n & mask) << 8) | ((n >> 8) & mask);
mask = 0x0000ffff0000ffffusize;
n = ((n & mask) << 16) | ((n >> 16) & mask);
n = (n << 32) | (n >> 32);
n
}
impl Polynomial {
fn trim(&mut self) {
let mut k = self.coefficients.len() - 1;
while 0 < k && self.coefficients[k].abs() < EPS {
k -= 1;
}
self.coefficients.truncate(k + 1);
}
fn evaluate_at(&self, x: f64) -> f64 {
let mut iter = self.coefficients.iter().rev();
// Всички полиноми имат поне един елемент по начина, по който ги създаваме
let mut value : f64 = *iter.next().unwrap();
for c in iter {
value = value * x + c;
}
value
}
}
Впечатляваща работа по задачата, и отвъд нея :). Оптимизациите ми се струват малко overkill за задачката, но се радвам, че ти е била достатъчно интересна да положиш усилия да ги приложиш. Получаваш 2 бонус точки от мен за работата. Не мисля, че много ти пука за точките, но въпроса е принципен :).
Вместо да инстанцираш
Polynomial, да push-ваш коефициенти в него, и после да гоtrim-ваш, би могъл да си алокираш вектор, и после просто да викнешPolynomial::from(...). Това по принцип ни беше идеята -- trim-ването на коефициентите да се направи при конструиране на полинома.Не е голяма работа, просто нещо дребно, което може леко да съкрати кода и да ти даде малко повече reliability -- ако винаги конструираш вектор по един и същ начин, и винаги го нормализираш при конструиране, ще можеш смело да викаш
from, без на всяко място, където изграждаш вектор, да мислиш дали си забравил нещо.На пръв поглед е странно, но идеята е
+да взима ownership и да върне стойност. ИмаAddAssign (+=)което е 'in place' и взима&mut self.Не ме притеснява, че взима ownership, a че след като го вземем няма как да го върнем (защото функцията приема два аргумента а връщаме само един полином) и освен това не правим нищо с тези обекти след като им вземем ownership (тоест ще ги оставим да умрат, защото ако се опитаме да пишем в тях компилатора се оплаква, че пипаме мютъбъл вериъбъл). Ако беше &self, &rhs старите обекти щяха да си останат. А ако беше mut self, mut rhs, поне след взимането на ownership, смъртта им нямаше да е напразна защото поне има как да ползваме алокираните за тях ресурси и да конструираме резултата в тях. Тоест нещо като sometype_t operator+(sometype_t&&, sometype_t&&) в C++ (Само като пример го давам там би било тъпо да се ползва такова нещо защото има много по мощен начин за избягване на копирания и временни обекти от такива аритметични изрази. Чрез expression templates можем в compile time да парснем дървото на израза и винаги да изберем оптималния начин за смятане избягващ максимален брой копирания).
Все пак не съм сигурен дали разбирам как точно работи това. Тоест има ли как да извлечем каквато и да е полза от това че тези обекти им взимаме ownership в тая функция, а те не са mutable?
Всъщност, можеш да дефинираш функцията като
mut self, и това няма да гръмне. В случая, параметъра е еквивалентно на деклариране на променлива, така че простоselfе все едноlet self = ..., аmut selfще еlet mut self = ..., но и двете работят със сигнатурата. Пробвай. Признавам, че може би трябваше да го напишем в условието така, или да го обясним на някоя лекция.И да, това, че можеш да мутираш self, означава че, на теория, не е нужно да конструираш нов полином, а можеш да използваш десния. Но да, признавам, че и на мен ми е малко неинтуитивно.