use crate::core::integer::FullInt;
use crate::{
Circle, CircleConstants, ExponentConstants, FractionConstants, Integer, Scalar, ScalarConstants,
};
use i256::I256;
use num_traits::{AsPrimitive, WrappingAdd, WrappingMul, WrappingNeg, WrappingSub};
use core::ops::*;
trait RandomFraction {
fn random() -> Self;
}
impl RandomFraction for i8 {
fn random() -> Self {
rand::random()
}
}
impl RandomFraction for i16 {
fn random() -> Self {
rand::random()
}
}
impl RandomFraction for i32 {
fn random() -> Self {
rand::random()
}
}
impl RandomFraction for i64 {
fn random() -> Self {
rand::random()
}
}
impl RandomFraction for i128 {
fn random() -> Self {
rand::random()
}
}
#[allow(private_bounds)]
impl<
F: Integer
+ FractionConstants
+ FullInt
+ RandomFraction
+ Shl<isize, Output = F>
+ Shr<isize, Output = F>
+ Shl<F, Output = F>
+ Shr<F, Output = F>
+ Shl<E, Output = F>
+ Shr<E, Output = F>
+ WrappingNeg
+ WrappingAdd
+ WrappingMul
+ WrappingSub,
E: Integer
+ ExponentConstants
+ FullInt
+ Shl<isize, Output = E>
+ Shr<isize, Output = E>
+ Shl<E, Output = E>
+ Shr<E, Output = E>
+ Shl<F, Output = E>
+ Shr<F, Output = E>
+ WrappingNeg
+ WrappingAdd
+ WrappingMul
+ WrappingSub,
> Scalar<F, E>
where
Scalar<F, E>: ScalarConstants,
u8: AsPrimitive<F>,
u16: AsPrimitive<F>,
u32: AsPrimitive<F>,
u64: AsPrimitive<F>,
u128: AsPrimitive<F>,
usize: AsPrimitive<F>,
i8: AsPrimitive<F>,
i16: AsPrimitive<F>,
i32: AsPrimitive<F>,
i64: AsPrimitive<F>,
i128: AsPrimitive<F>,
isize: AsPrimitive<F>,
I256: From<F>,
u8: AsPrimitive<E>,
u16: AsPrimitive<E>,
u32: AsPrimitive<E>,
u64: AsPrimitive<E>,
u128: AsPrimitive<E>,
usize: AsPrimitive<E>,
i8: AsPrimitive<E>,
i16: AsPrimitive<E>,
i32: AsPrimitive<E>,
i64: AsPrimitive<E>,
i128: AsPrimitive<E>,
isize: AsPrimitive<E>,
I256: From<E>,
{
#[inline]
pub fn random() -> Self {
let mut result = Self {
fraction: F::random(),
exponent: 0.as_(),
};
loop {
let leading: E = result
.fraction
.leading_ones()
.max(result.fraction.leading_zeros())
.as_();
if leading > E::ONE {
let new_exponent = result.exponent.wrapping_sub(&leading);
let shift: isize = leading.as_();
if new_exponent.is_negative() {
result.exponent = new_exponent.wrapping_add(&E::ONE);
result.fraction = result.fraction << shift.wrapping_sub(1);
let mask: F = (F::ONE << shift.wrapping_sub(1)).wrapping_sub(&F::ONE);
let random_fill: F = F::random() & mask;
result.fraction = result.fraction | random_fill;
} else {
result.exponent = E::AMBIGUOUS_EXPONENT;
result.fraction = F::random();
let mut leading: E = result
.fraction
.leading_ones()
.max(result.fraction.leading_zeros())
.as_();
while leading != E::TWO {
result.fraction = F::random();
leading = result
.fraction
.leading_ones()
.max(result.fraction.leading_zeros())
.as_();
}
break;
}
} else {
break;
}
}
result
}
#[inline]
pub fn random_gauss() -> Self {
let mut s: Self;
let mut u: Self;
let mut v: Self;
loop {
u = Self::random();
v = Self::random();
s = u.square() + v.square();
if !s.exponent.is_positive() {
break;
}
}
if !s.is_normal() {
let mut normal = Self::random();
normal.normalize_vanished();
return normal;
}
let multiplier: Self = -2 * s.ln() / s;
if !multiplier.is_normal() {
let mut normal = Self::random();
normal.normalize_exploded();
return normal;
}
multiplier.sqrt() * u
}
}
#[allow(private_bounds)]
impl<
F: Integer
+ FractionConstants
+ FullInt
+ RandomFraction
+ Shl<isize, Output = F>
+ Shr<isize, Output = F>
+ Shl<F, Output = F>
+ Shr<F, Output = F>
+ Shl<E, Output = F>
+ Shr<E, Output = F>
+ WrappingNeg
+ WrappingAdd
+ WrappingMul
+ WrappingSub,
E: Integer
+ ExponentConstants
+ FullInt
+ Shl<isize, Output = E>
+ Shr<isize, Output = E>
+ Shl<E, Output = E>
+ Shr<E, Output = E>
+ Shl<F, Output = E>
+ Shr<F, Output = E>
+ WrappingNeg
+ WrappingAdd
+ WrappingMul
+ WrappingSub,
> Circle<F, E>
where
Circle<F, E>: CircleConstants,
Scalar<F, E>: ScalarConstants,
u8: AsPrimitive<F>,
u16: AsPrimitive<F>,
u32: AsPrimitive<F>,
u64: AsPrimitive<F>,
u128: AsPrimitive<F>,
usize: AsPrimitive<F>,
i8: AsPrimitive<F>,
i16: AsPrimitive<F>,
i32: AsPrimitive<F>,
i64: AsPrimitive<F>,
i128: AsPrimitive<F>,
isize: AsPrimitive<F>,
I256: From<F>,
u8: AsPrimitive<E>,
u16: AsPrimitive<E>,
u32: AsPrimitive<E>,
u64: AsPrimitive<E>,
u128: AsPrimitive<E>,
usize: AsPrimitive<E>,
i8: AsPrimitive<E>,
i16: AsPrimitive<E>,
i32: AsPrimitive<E>,
i64: AsPrimitive<E>,
i128: AsPrimitive<E>,
isize: AsPrimitive<E>,
I256: From<E>,
{
#[inline]
pub fn random() -> Self {
loop {
let mut result = Self {
real: F::random(),
imaginary: F::random(),
exponent: 0.as_(),
};
let leading_r: E = result
.real
.leading_ones()
.max(result.real.leading_zeros())
.as_();
let leading_i: E = result
.imaginary
.leading_ones()
.max(result.imaginary.leading_zeros())
.as_();
let leading = leading_r.min(leading_i);
if leading > E::ONE {
let new_exponent = result.exponent.wrapping_sub(&leading);
let shift: isize = leading.as_();
if new_exponent.is_negative() {
result.exponent = new_exponent.wrapping_add(&E::ONE);
result.real = result.real << shift.wrapping_sub(1);
result.imaginary = result.imaginary << shift.wrapping_sub(1);
let mask: F =
(F::POS_ONE_FRACTION << shift.wrapping_sub(1)).wrapping_sub(&F::ONE);
let random_fill_r: F = F::random() & mask;
let random_fill_i: F = F::random() & mask;
result.real = result.real | random_fill_r;
result.imaginary = result.imaginary | random_fill_i;
} else {
result.exponent = E::AMBIGUOUS_EXPONENT;
loop {
result.real = F::random();
result.imaginary = F::random();
let leading_r: E = result
.real
.leading_ones()
.max(result.real.leading_zeros())
.as_();
let leading_i: E = result
.imaginary
.leading_ones()
.max(result.imaginary.leading_zeros())
.as_();
if leading_r.min(leading_i) == E::TWO {
break;
}
}
}
}
if !result.magnitude_squared().exponent.is_positive() {
return result;
}
}
}
pub fn random_gauss() -> Self {
let mut s: Scalar<F, E>;
let mut u: Scalar<F, E>;
let mut v: Scalar<F, E>;
loop {
u = Scalar::<F, E>::random();
v = Scalar::<F, E>::random();
s = u.square() + v.square();
if !s.exponent.is_positive() {
break;
}
}
if !s.is_normal() {
let mut normal = Self::random();
normal.normalize_vanished();
return normal;
}
let multiplier: Scalar<F, E> = -2 * s.ln() / s;
if !multiplier.is_normal() {
let mut normal = Self::random();
normal.normalize_exploded();
return normal;
}
let factor = multiplier.sqrt();
Circle::<F, E>::from((factor * u, factor * v))
}
}