1 use crate::error::{FendError, Interrupt}; 2 use crate::num::real::{self, Real}; 3 use crate::num::Exact; 4 use crate::num::{Base, FormattingStyle}; 5 use crate::result::FResult; 6 use std::cmp::Ordering; 7 use std::ops::Neg; 8 use std::{fmt, io}; 9 10 #[derive(Clone, Hash)] 11 pub(crate) struct Complex { 12 real: Real, 13 imag: Real, 14 } 15 16 impl fmt::Debug for Complex { fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result17 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { 18 write!(f, "{:?}", self.real)?; 19 if !self.imag.is_definitely_zero() { 20 write!(f, " + {:?}i", self.imag)?; 21 } 22 Ok(()) 23 } 24 } 25 26 #[derive(Copy, Clone, Eq, PartialEq)] 27 pub(crate) enum UseParentheses { 28 No, 29 IfComplex, 30 IfComplexOrFraction, 31 } 32 33 impl Complex { compare<I: Interrupt>(&self, other: &Self, int: &I) -> FResult<Option<Ordering>>34 pub(crate) fn compare<I: Interrupt>(&self, other: &Self, int: &I) -> FResult<Option<Ordering>> { 35 if self.imag().is_zero() && other.imag().is_zero() { 36 Ok(Some(self.real().compare(&other.real(), int)?)) 37 } else { 38 Ok(None) 39 } 40 } 41 serialize(&self, write: &mut impl io::Write) -> FResult<()>42 pub(crate) fn serialize(&self, write: &mut impl io::Write) -> FResult<()> { 43 self.real.serialize(write)?; 44 self.imag.serialize(write)?; 45 Ok(()) 46 } 47 deserialize(read: &mut impl io::Read) -> FResult<Self>48 pub(crate) fn deserialize(read: &mut impl io::Read) -> FResult<Self> { 49 Ok(Self { 50 real: Real::deserialize(read)?, 51 imag: Real::deserialize(read)?, 52 }) 53 } 54 try_as_usize<I: Interrupt>(self, int: &I) -> FResult<usize>55 pub(crate) fn try_as_usize<I: Interrupt>(self, int: &I) -> FResult<usize> { 56 if !self.imag.is_zero() { 57 return Err(FendError::ComplexToInteger); 58 } 59 self.real.try_as_usize(int) 60 } 61 try_as_i64<I: Interrupt>(self, int: &I) -> FResult<i64>62 pub(crate) fn try_as_i64<I: Interrupt>(self, int: &I) -> FResult<i64> { 63 if !self.imag.is_zero() { 64 return Err(FendError::ComplexToInteger); 65 } 66 self.real.try_as_i64(int) 67 } 68 69 #[inline] real(&self) -> Real70 pub(crate) fn real(&self) -> Real { 71 self.real.clone() 72 } 73 74 #[inline] imag(&self) -> Real75 pub(crate) fn imag(&self) -> Real { 76 self.imag.clone() 77 } 78 conjugate(self) -> Self79 pub(crate) fn conjugate(self) -> Self { 80 Self { 81 real: self.real, 82 imag: -self.imag, 83 } 84 } 85 factorial<I: Interrupt>(self, int: &I) -> FResult<Self>86 pub(crate) fn factorial<I: Interrupt>(self, int: &I) -> FResult<Self> { 87 if !self.imag.is_zero() { 88 return Err(FendError::FactorialComplex); 89 } 90 Ok(Self { 91 real: self.real.factorial(int)?, 92 imag: self.imag, 93 }) 94 } 95 exp<I: Interrupt>(self, int: &I) -> FResult<Exact<Self>>96 pub(crate) fn exp<I: Interrupt>(self, int: &I) -> FResult<Exact<Self>> { 97 // e^(a + bi) = e^a * e^(bi) = e^a * (cos(b) + i * sin(b)) 98 let r = self.real.exp(int)?; 99 let exact = r.exact; 100 101 Ok(Exact::new( 102 Self { 103 real: r.clone().mul(self.imag.clone().cos(int)?.re(), int)?.value, 104 imag: r.mul(self.imag.sin(int)?.re(), int)?.value, 105 }, 106 exact, 107 )) 108 } 109 pow_n<I: Interrupt>(self, n: usize, int: &I) -> FResult<Exact<Self>>110 fn pow_n<I: Interrupt>(self, n: usize, int: &I) -> FResult<Exact<Self>> { 111 if n == 0 { 112 return Ok(Exact::new(Self::from(1), true)); 113 } 114 115 let mut result = Exact::new(Self::from(1), true); 116 let mut base = Exact::new(self, true); 117 let mut exp = n; 118 while exp > 0 { 119 if exp & 1 == 1 { 120 result = result.mul(&base, int)?; 121 } 122 base = base.clone().mul(&base, int)?; 123 exp >>= 1; 124 } 125 126 Ok(result) 127 } 128 pow<I: Interrupt>(self, rhs: Self, int: &I) -> FResult<Exact<Self>>129 pub(crate) fn pow<I: Interrupt>(self, rhs: Self, int: &I) -> FResult<Exact<Self>> { 130 if !rhs.real.is_integer() || !rhs.imag.is_integer() { 131 return self.frac_pow(rhs, int); 132 } 133 134 if self.imag.is_zero() && rhs.imag.is_zero() { 135 let real = self.real.pow(rhs.real, int)?; 136 Ok(Exact::new( 137 Self { 138 real: real.value, 139 imag: 0.into(), 140 }, 141 real.exact, 142 )) 143 } else { 144 let rem = rhs.clone().real.modulo(4.into(), int); 145 // Reduced case: (ix)^y = x^y * i^y 146 if self.real.is_zero() && rhs.imag.is_zero() { 147 if let Ok(n) = rhs.real.try_as_usize(int) { 148 return self.pow_n(n, int); 149 } 150 151 let mut result = Exact::new( 152 match rem.and_then(|rem| rem.try_as_usize(int)) { 153 Ok(0) => 1.into(), 154 Ok(1) => Self { 155 real: 0.into(), 156 imag: 1.into(), 157 }, 158 Ok(2) => Self { 159 real: Real::from(1).neg(), 160 imag: 0.into(), 161 }, 162 Ok(3) => Self { 163 real: 0.into(), 164 imag: Real::from(1).neg(), 165 }, 166 Ok(_) | Err(_) => unreachable!("modulo 4 should always be 0, 1, 2, or 3"), 167 }, 168 true, 169 ); 170 171 if !self.imag.is_definitely_one() { 172 result = self 173 .imag 174 .pow(2.into(), int)? 175 .apply(Self::from) 176 .mul(&result, int)?; 177 } 178 179 return Ok(result); 180 } 181 182 // z^w = e^(w * ln(z)) 183 let exp = self.ln(int)?.mul(&Exact::new(rhs, true), int)?; 184 Ok(exp.value.exp(int)?.combine(exp.exact)) 185 } 186 } 187 i() -> Self188 pub(crate) fn i() -> Self { 189 Self { 190 real: 0.into(), 191 imag: 1.into(), 192 } 193 } 194 pi() -> Self195 pub(crate) fn pi() -> Self { 196 Self { 197 real: Real::pi(), 198 imag: 0.into(), 199 } 200 } 201 abs<I: Interrupt>(self, int: &I) -> FResult<Exact<Real>>202 pub(crate) fn abs<I: Interrupt>(self, int: &I) -> FResult<Exact<Real>> { 203 Ok(if self.imag.is_zero() { 204 if self.real.is_neg() { 205 Exact::new(-self.real, true) 206 } else { 207 Exact::new(self.real, true) 208 } 209 } else if self.real.is_zero() { 210 if self.imag.is_neg() { 211 Exact::new(-self.imag, true) 212 } else { 213 Exact::new(self.imag, true) 214 } 215 } else { 216 let power = self.real.pow(2.into(), int)?; 217 let power2 = self.imag.pow(2.into(), int)?; 218 let real = power.add(power2, int)?; 219 let result = real.value.root_n(&Real::from(2), int)?; 220 result.combine(real.exact) 221 }) 222 } 223 floor<I: Interrupt>(self, int: &I) -> FResult<Exact<Real>>224 pub(crate) fn floor<I: Interrupt>(self, int: &I) -> FResult<Exact<Real>> { 225 Ok(Exact::new(self.expect_real()?.floor(int)?, true)) 226 } 227 ceil<I: Interrupt>(self, int: &I) -> FResult<Exact<Real>>228 pub(crate) fn ceil<I: Interrupt>(self, int: &I) -> FResult<Exact<Real>> { 229 Ok(Exact::new(self.expect_real()?.ceil(int)?, true)) 230 } 231 round<I: Interrupt>(self, int: &I) -> FResult<Exact<Real>>232 pub(crate) fn round<I: Interrupt>(self, int: &I) -> FResult<Exact<Real>> { 233 Ok(Exact::new(self.expect_real()?.round(int)?, true)) 234 } 235 arg<I: Interrupt>(self, int: &I) -> FResult<Exact<Real>>236 pub(crate) fn arg<I: Interrupt>(self, int: &I) -> FResult<Exact<Real>> { 237 Ok(Exact::new(self.imag.atan2(self.real, int)?, false)) 238 } 239 format<I: Interrupt>( &self, exact: bool, style: FormattingStyle, base: Base, use_parentheses: UseParentheses, int: &I, ) -> FResult<Exact<Formatted>>240 pub(crate) fn format<I: Interrupt>( 241 &self, 242 exact: bool, 243 style: FormattingStyle, 244 base: Base, 245 use_parentheses: UseParentheses, 246 int: &I, 247 ) -> FResult<Exact<Formatted>> { 248 let style = if !exact && style == FormattingStyle::Auto { 249 FormattingStyle::DecimalPlaces(10) 250 } else if !self.imag.is_zero() && style == FormattingStyle::Auto { 251 FormattingStyle::Exact 252 } else { 253 style 254 }; 255 256 if self.imag.is_zero() { 257 let use_parens = use_parentheses == UseParentheses::IfComplexOrFraction; 258 let x = self.real.format(base, style, false, use_parens, int)?; 259 return Ok(Exact::new( 260 Formatted { 261 first_component: x.value, 262 separator: "", 263 second_component: None, 264 use_parentheses: false, 265 }, 266 exact && x.exact, 267 )); 268 } 269 270 Ok(if self.real.is_zero() { 271 let use_parens = use_parentheses == UseParentheses::IfComplexOrFraction; 272 let x = self.imag.format(base, style, true, use_parens, int)?; 273 Exact::new( 274 Formatted { 275 first_component: x.value, 276 separator: "", 277 second_component: None, 278 use_parentheses: false, 279 }, 280 exact && x.exact, 281 ) 282 } else { 283 let mut exact = exact; 284 let real_part = self.real.format(base, style, false, false, int)?; 285 exact = exact && real_part.exact; 286 let (positive, imag_part) = if self.imag.is_pos() { 287 (true, self.imag.format(base, style, true, false, int)?) 288 } else { 289 ( 290 false, 291 (-self.imag.clone()).format(base, style, true, false, int)?, 292 ) 293 }; 294 exact = exact && imag_part.exact; 295 let separator = if positive { " + " } else { " - " }; 296 Exact::new( 297 Formatted { 298 first_component: real_part.value, 299 separator, 300 second_component: Some(imag_part.value), 301 use_parentheses: use_parentheses == UseParentheses::IfComplex 302 || use_parentheses == UseParentheses::IfComplexOrFraction, 303 }, 304 exact, 305 ) 306 }) 307 } frac_pow<I: Interrupt>(self, n: Self, int: &I) -> FResult<Exact<Self>>308 pub(crate) fn frac_pow<I: Interrupt>(self, n: Self, int: &I) -> FResult<Exact<Self>> { 309 if self.imag.is_zero() && n.imag.is_zero() && !self.real.is_neg() { 310 Ok(self.real.pow(n.real, int)?.apply(Self::from)) 311 } else { 312 let exponent = self.ln(int)?.mul(&Exact::new(n, true), int)?; 313 Ok(exponent.value.exp(int)?.combine(exponent.exact)) 314 } 315 } 316 expect_real(self) -> FResult<Real>317 fn expect_real(self) -> FResult<Real> { 318 if self.imag.is_zero() { 319 Ok(self.real) 320 } else { 321 Err(FendError::ExpectedARealNumber) 322 } 323 } 324 sin<I: Interrupt>(self, int: &I) -> FResult<Exact<Self>>325 pub(crate) fn sin<I: Interrupt>(self, int: &I) -> FResult<Exact<Self>> { 326 // sin(a + bi) = sin(a) * cosh(b) + i * cos(a) * sinh(b) 327 if self.imag.is_zero() { 328 Ok(self.real.sin(int)?.apply(Self::from)) 329 } else { 330 let cosh = Exact::new(self.imag.clone().cosh(int)?, false); 331 let sinh = Exact::new(self.imag.sinh(int)?, false); 332 333 let real = self.real.clone().sin(int)?.mul(cosh.re(), int)?; 334 let imag = self.real.cos(int)?.mul(sinh.re(), int)?; 335 336 Ok(Exact::new( 337 Self { 338 real: real.value, 339 imag: imag.value, 340 }, 341 real.exact && imag.exact, 342 )) 343 } 344 } 345 cos<I: Interrupt>(self, int: &I) -> FResult<Exact<Self>>346 pub(crate) fn cos<I: Interrupt>(self, int: &I) -> FResult<Exact<Self>> { 347 // cos(a + bi) = cos(a) * cosh(b) - i * sin(a) * sinh(b) 348 if self.imag.is_zero() { 349 Ok(self.real.cos(int)?.apply(Self::from)) 350 } else { 351 let cosh = Exact::new(self.imag.clone().cosh(int)?, false); 352 let sinh = Exact::new(self.imag.sinh(int)?, false); 353 let exact_real = Exact::new(self.real, true); 354 355 let real = exact_real.value.clone().cos(int)?.mul(cosh.re(), int)?; 356 let imag = exact_real.value.sin(int)?.mul(sinh.re(), int)?.neg(); 357 Ok(Exact::new( 358 Self { 359 real: real.value, 360 imag: imag.value, 361 }, 362 real.exact && imag.exact, 363 )) 364 } 365 } 366 tan<I: Interrupt>(self, int: &I) -> FResult<Exact<Self>>367 pub(crate) fn tan<I: Interrupt>(self, int: &I) -> FResult<Exact<Self>> { 368 let num = self.clone().sin(int)?; 369 let den = self.cos(int)?; 370 num.div(den, int) 371 } 372 373 /// Calculates ln(i * z + sqrt(1 - z^2)) 374 /// This is used to implement asin and acos for all complex numbers asin_ln<I: Interrupt>(self, int: &I) -> FResult<Exact<Self>>375 fn asin_ln<I: Interrupt>(self, int: &I) -> FResult<Exact<Self>> { 376 let half = Exact::new(Self::from(1), true).div(Exact::new(Self::from(2), true), int)?; 377 let exact = Exact::new(self, true); 378 let i = Exact::new(Self::i(), true); 379 380 let sqrt = Exact::new(Self::from(1), true) 381 .add(exact.clone().mul(&exact, int)?.neg(), int)? 382 .try_and_then(|x| x.frac_pow(half.value, int))?; 383 i.mul(&exact, int)? 384 .add(sqrt, int)? 385 .try_and_then(|x| x.ln(int)) 386 } 387 asin<I: Interrupt>(self, int: &I) -> FResult<Self>388 pub(crate) fn asin<I: Interrupt>(self, int: &I) -> FResult<Self> { 389 // Real asin is defined for -1 <= x <= 1 390 if self.imag.is_zero() && self.real.between_plus_minus_one_incl(int)? { 391 Ok(Self::from(self.real.asin(int)?)) 392 } else { 393 // asin(z) = -i * ln(i * z + sqrt(1 - z^2)) 394 Ok(self 395 .asin_ln(int)? 396 .mul(&Exact::new(Self::i(), true), int)? 397 .neg() 398 .value) 399 } 400 } 401 acos<I: Interrupt>(self, int: &I) -> FResult<Self>402 pub(crate) fn acos<I: Interrupt>(self, int: &I) -> FResult<Self> { 403 // Real acos is defined for -1 <= x <= 1 404 if self.imag.is_zero() && self.real.between_plus_minus_one_incl(int)? { 405 Ok(Self::from(self.real.acos(int)?)) 406 } else { 407 // acos(z) = pi/2 + i * ln(i * z + sqrt(1 - z^2)) 408 let half_pi = Exact::new(Self::pi(), true).div(Exact::new(Self::from(2), true), int)?; 409 Ok(half_pi 410 .add( 411 self.asin_ln(int)?.mul(&Exact::new(Self::i(), true), int)?, 412 int, 413 )? 414 .value) 415 } 416 } 417 atan<I: Interrupt>(self, int: &I) -> FResult<Self>418 pub(crate) fn atan<I: Interrupt>(self, int: &I) -> FResult<Self> { 419 // Real atan is defined for all real numbers 420 if self.imag.is_zero() { 421 Ok(Self::from(self.real.atan(int)?)) 422 } else { 423 // i/2 * (ln(-iz+1) - ln(iz+1)) 424 let half_i = Exact::new(Self::i(), true).div(Exact::new(Self::from(2), true), int)?; 425 let z1 = Exact::new(self.clone(), true) 426 .mul(&Exact::new(Self::i().neg(), true), int)? 427 .add(Exact::new(Self::from(1), true), int)?; 428 let z2 = Exact::new(self, true) 429 .mul(&Exact::new(Self::i(), true), int)? 430 .add(Exact::new(Self::from(1), true), int)?; 431 432 Ok(half_i 433 .mul( 434 &z1.try_and_then(|z| z.ln(int))? 435 .add(z2.try_and_then(|z| z.ln(int))?.neg(), int)?, 436 int, 437 )? 438 .value) 439 } 440 } 441 sinh<I: Interrupt>(self, int: &I) -> FResult<Self>442 pub(crate) fn sinh<I: Interrupt>(self, int: &I) -> FResult<Self> { 443 if self.imag.is_zero() { 444 Ok(Self::from(self.real.sinh(int)?)) 445 } else { 446 // sinh(a+bi)=sinh(a)cos(b)+icosh(a)sin(b) 447 let sinh = Exact::new(self.real.clone().sinh(int)?, false); 448 let cos = self.imag.clone().cos(int)?; 449 let cosh = Exact::new(self.real.cosh(int)?, false); 450 let sin = self.imag.sin(int)?; 451 452 Ok(Self { 453 real: sinh.mul(cos.re(), int)?.value, 454 imag: cosh.mul(sin.re(), int)?.value, 455 }) 456 } 457 } 458 cosh<I: Interrupt>(self, int: &I) -> FResult<Self>459 pub(crate) fn cosh<I: Interrupt>(self, int: &I) -> FResult<Self> { 460 if self.imag.is_zero() { 461 Ok(Self::from(self.real.cosh(int)?)) 462 } else { 463 // cosh(a+bi)=cosh(a)cos(b)+isinh(a)sin(b) 464 let cosh = Exact::new(self.real.clone().cosh(int)?, false); 465 let cos = self.imag.clone().cos(int)?; 466 let sinh = Exact::new(self.real.sinh(int)?, false); 467 let sin = self.imag.sin(int)?; 468 469 Ok(Self { 470 real: cosh.mul(cos.re(), int)?.value, 471 imag: sinh.mul(sin.re(), int)?.value, 472 }) 473 } 474 } 475 tanh<I: Interrupt>(self, int: &I) -> FResult<Self>476 pub(crate) fn tanh<I: Interrupt>(self, int: &I) -> FResult<Self> { 477 if self.imag.is_zero() { 478 Ok(Self::from(self.real.tanh(int)?)) 479 } else { 480 // tanh(a+bi)=sinh(a+bi)/cosh(a+bi) 481 Exact::new(self.clone().sinh(int)?, false) 482 .div(Exact::new(self.cosh(int)?, false), int) 483 .map(|x| x.value) 484 } 485 } 486 asinh<I: Interrupt>(self, int: &I) -> FResult<Self>487 pub(crate) fn asinh<I: Interrupt>(self, int: &I) -> FResult<Self> { 488 // Real asinh is defined for all real numbers 489 if self.imag.is_zero() { 490 Ok(Self::from(self.real.asinh(int)?)) 491 } else { 492 // asinh(z)=ln(z+sqrt(z^2+1)) 493 let exact = Exact::new(self, true); 494 let half = Exact::new(Self::from(1), true).div(Exact::new(Self::from(2), true), int)?; 495 let sqrt = exact 496 .clone() 497 .mul(&exact, int)? 498 .add(Exact::new(Self::from(1), true), int)? 499 .try_and_then(|x| x.frac_pow(half.value, int))?; 500 sqrt.add(exact, int)? 501 .try_and_then(|x| x.ln(int)) 502 .map(|x| x.value) 503 } 504 } 505 acosh<I: Interrupt>(self, int: &I) -> FResult<Self>506 pub(crate) fn acosh<I: Interrupt>(self, int: &I) -> FResult<Self> { 507 // Real acosh is defined for x >= 1 508 if self.imag.is_zero() && self.real.compare(&1.into(), int)? != Ordering::Less { 509 Ok(Self::from(self.real.acosh(int)?)) 510 } else { 511 // acosh(z)=ln(z+sqrt(z^2-1)) 512 let exact = Exact::new(self, true); 513 let half = Exact::new(Self::from(1), true).div(Exact::new(Self::from(2), true), int)?; 514 let sqrt = exact 515 .clone() 516 .mul(&exact, int)? 517 .add(Exact::new(Self::from(1), true).neg(), int)? 518 .try_and_then(|x| x.frac_pow(half.value, int))?; 519 sqrt.add(exact, int)? 520 .try_and_then(|x| x.ln(int)) 521 .map(|x| x.value) 522 } 523 } 524 atanh<I: Interrupt>(self, int: &I) -> FResult<Self>525 pub(crate) fn atanh<I: Interrupt>(self, int: &I) -> FResult<Self> { 526 // Real atanh is defined for -1 < x < 1 527 // Undefined for x = 1, -1 528 if self.imag.is_zero() && self.real.between_plus_minus_one_excl(int)? { 529 Ok(Self::from(self.real.atanh(int)?)) 530 } else { 531 // atanh(z)=ln(sqrt(-(z-1)/(z-1))) 532 let exact = Exact::new(self, true); 533 let one = Exact::new(Self::from(1), true); 534 let half = Exact::new(Self::from(1), true).div(Exact::new(Self::from(2), true), int)?; 535 536 exact 537 .clone() 538 .add(one.clone(), int)? 539 .neg() 540 .div(exact.add(one.neg(), int)?, int)? 541 .try_and_then(|x| x.frac_pow(half.value, int))? 542 .try_and_then(|z| z.ln(int)) 543 .map(|x| x.value) 544 } 545 } 546 ln<I: Interrupt>(self, int: &I) -> FResult<Exact<Self>>547 pub(crate) fn ln<I: Interrupt>(self, int: &I) -> FResult<Exact<Self>> { 548 if self.imag.is_zero() && self.real.is_pos() { 549 Ok(self.real.ln(int)?.apply(Self::from)) 550 } else { 551 // ln(z) = ln(|z|) + i * arg(z) 552 let abs = self.clone().abs(int)?; 553 let real = abs.value.ln(int)?.combine(abs.exact); 554 let arg = self.arg(int)?; 555 556 Ok(Exact::new( 557 Self { 558 real: real.value, 559 imag: arg.value, 560 }, 561 real.exact && arg.exact, 562 )) 563 } 564 } 565 log<I: Interrupt>(self, base: Self, int: &I) -> FResult<Self>566 pub(crate) fn log<I: Interrupt>(self, base: Self, int: &I) -> FResult<Self> { 567 // log_n(z) = ln(z) / ln(n) 568 let ln = self.ln(int)?; 569 let ln2 = base.ln(int)?; 570 Ok(ln.div(ln2, int)?.value) 571 } 572 log2<I: Interrupt>(self, int: &I) -> FResult<Self>573 pub(crate) fn log2<I: Interrupt>(self, int: &I) -> FResult<Self> { 574 if self.imag.is_zero() && self.real.is_pos() { 575 Ok(Self::from(self.real.log2(int)?)) 576 } else { 577 self.log(Self::from(2), int) 578 } 579 } log10<I: Interrupt>(self, int: &I) -> FResult<Self>580 pub(crate) fn log10<I: Interrupt>(self, int: &I) -> FResult<Self> { 581 if self.imag.is_zero() && self.real.is_pos() { 582 Ok(Self::from(self.real.log10(int)?)) 583 } else { 584 self.log(Self::from(10), int) 585 } 586 } 587 is_definitely_one(&self) -> bool588 pub(crate) fn is_definitely_one(&self) -> bool { 589 self.real.is_definitely_one() && self.imag.is_definitely_zero() 590 } 591 modulo<I: Interrupt>(self, rhs: Self, int: &I) -> FResult<Self>592 pub(crate) fn modulo<I: Interrupt>(self, rhs: Self, int: &I) -> FResult<Self> { 593 Ok(Self::from( 594 self.expect_real()?.modulo(rhs.expect_real()?, int)?, 595 )) 596 } 597 bitwise<I: Interrupt>( self, rhs: Self, op: crate::ast::BitwiseBop, int: &I, ) -> FResult<Self>598 pub(crate) fn bitwise<I: Interrupt>( 599 self, 600 rhs: Self, 601 op: crate::ast::BitwiseBop, 602 int: &I, 603 ) -> FResult<Self> { 604 Ok(Self::from(self.expect_real()?.bitwise( 605 rhs.expect_real()?, 606 op, 607 int, 608 )?)) 609 } 610 combination<I: Interrupt>(self, rhs: Self, int: &I) -> FResult<Self>611 pub(crate) fn combination<I: Interrupt>(self, rhs: Self, int: &I) -> FResult<Self> { 612 Ok(Self::from( 613 self.expect_real()?.combination(rhs.expect_real()?, int)?, 614 )) 615 } 616 permutation<I: Interrupt>(self, rhs: Self, int: &I) -> FResult<Self>617 pub(crate) fn permutation<I: Interrupt>(self, rhs: Self, int: &I) -> FResult<Self> { 618 Ok(Self::from( 619 self.expect_real()?.permutation(rhs.expect_real()?, int)?, 620 )) 621 } 622 } 623 624 impl Exact<Complex> { add<I: Interrupt>(self, rhs: Self, int: &I) -> FResult<Self>625 pub(crate) fn add<I: Interrupt>(self, rhs: Self, int: &I) -> FResult<Self> { 626 let (self_real, self_imag) = self.apply(|x| (x.real, x.imag)).pair(); 627 let (rhs_real, rhs_imag) = rhs.apply(|x| (x.real, x.imag)).pair(); 628 let real = self_real.add(rhs_real, int)?; 629 let imag = self_imag.add(rhs_imag, int)?; 630 Ok(Self::new( 631 Complex { 632 real: real.value, 633 imag: imag.value, 634 }, 635 real.exact && imag.exact, 636 )) 637 } 638 mul<I: Interrupt>(self, rhs: &Self, int: &I) -> FResult<Self>639 pub(crate) fn mul<I: Interrupt>(self, rhs: &Self, int: &I) -> FResult<Self> { 640 // (a + bi) * (c + di) 641 // => ac + bci + adi - bd 642 // => (ac - bd) + (bc + ad)i 643 let (self_real, self_imag) = self.apply(|x| (x.real, x.imag)).pair(); 644 let (rhs_real, rhs_imag) = rhs.clone().apply(|x| (x.real, x.imag)).pair(); 645 646 let prod1 = self_real.clone().mul(rhs_real.re(), int)?; 647 let prod2 = self_imag.clone().mul(rhs_imag.re(), int)?; 648 let real_part = prod1.add(-prod2, int)?; 649 let prod3 = self_real.mul(rhs_imag.re(), int)?; 650 let prod4 = self_imag.mul(rhs_real.re(), int)?; 651 let imag_part = prod3.add(prod4, int)?; 652 Ok(Self::new( 653 Complex { 654 real: real_part.value, 655 imag: imag_part.value, 656 }, 657 real_part.exact && imag_part.exact, 658 )) 659 } 660 div<I: Interrupt>(self, rhs: Self, int: &I) -> FResult<Self>661 pub(crate) fn div<I: Interrupt>(self, rhs: Self, int: &I) -> FResult<Self> { 662 // (u + vi) / (x + yi) = (1/(x^2 + y^2)) * ((ux + vy) + (vx - uy)i) 663 let (u, v) = self.apply(|x| (x.real, x.imag)).pair(); 664 let (x, y) = rhs.apply(|x| (x.real, x.imag)).pair(); 665 // if both numbers are real, use this simplified algorithm 666 if v.exact && v.value.is_zero() && y.exact && y.value.is_zero() { 667 return Ok(u.div(&x, int)?.apply(|real| Complex { 668 real, 669 imag: 0.into(), 670 })); 671 } 672 let prod1 = x.clone().mul(x.re(), int)?; 673 let prod2 = y.clone().mul(y.re(), int)?; 674 let sum = prod1.add(prod2, int)?; 675 let real_part = Exact::new(Real::from(1), true).div(&sum, int)?; 676 let prod3 = u.clone().mul(x.re(), int)?; 677 let prod4 = v.clone().mul(y.re(), int)?; 678 let real2 = prod3.add(prod4, int)?; 679 let prod5 = v.mul(x.re(), int)?; 680 let prod6 = u.mul(y.re(), int)?; 681 let imag2 = prod5.add(-prod6, int)?; 682 let multiplicand = Self::new( 683 Complex { 684 real: real2.value, 685 imag: imag2.value, 686 }, 687 real2.exact && imag2.exact, 688 ); 689 let result = Self::new( 690 Complex { 691 real: real_part.value, 692 imag: 0.into(), 693 }, 694 real_part.exact, 695 ) 696 .mul(&multiplicand, int)?; 697 Ok(result) 698 } 699 } 700 701 impl Neg for Complex { 702 type Output = Self; 703 neg(self) -> Self704 fn neg(self) -> Self { 705 Self { 706 real: -self.real, 707 imag: -self.imag, 708 } 709 } 710 } 711 712 impl Neg for &Complex { 713 type Output = Complex; 714 neg(self) -> Complex715 fn neg(self) -> Complex { 716 -self.clone() 717 } 718 } 719 720 impl From<u64> for Complex { from(i: u64) -> Self721 fn from(i: u64) -> Self { 722 Self { 723 real: i.into(), 724 imag: 0.into(), 725 } 726 } 727 } 728 729 impl From<Real> for Complex { from(i: Real) -> Self730 fn from(i: Real) -> Self { 731 Self { 732 real: i, 733 imag: 0.into(), 734 } 735 } 736 } 737 738 #[derive(Debug)] 739 pub(crate) struct Formatted { 740 first_component: real::Formatted, 741 separator: &'static str, 742 second_component: Option<real::Formatted>, 743 use_parentheses: bool, 744 } 745 746 impl fmt::Display for Formatted { fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result747 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { 748 if self.use_parentheses { 749 write!(f, "(")?; 750 } 751 write!(f, "{}{}", self.first_component, self.separator)?; 752 if let Some(second_component) = &self.second_component { 753 write!(f, "{second_component}")?; 754 } 755 if self.use_parentheses { 756 write!(f, ")")?; 757 } 758 Ok(()) 759 } 760 } 761