diff options
author | Jedidiah Barber <contact@jedbarber.id.au> | 2022-11-16 01:06:13 +1300 |
---|---|---|
committer | Jedidiah Barber <contact@jedbarber.id.au> | 2022-11-16 01:06:13 +1300 |
commit | 0d8f4e413ba7a1f64753d73b757d7cc45e363369 (patch) | |
tree | 0ae7c956f21015b507f5090afa0f852fa63a84df /src/complex_fixed_points.adb |
Initial attempt
Diffstat (limited to 'src/complex_fixed_points.adb')
-rw-r--r-- | src/complex_fixed_points.adb | 529 |
1 files changed, 529 insertions, 0 deletions
diff --git a/src/complex_fixed_points.adb b/src/complex_fixed_points.adb new file mode 100644 index 0000000..a0c7be1 --- /dev/null +++ b/src/complex_fixed_points.adb @@ -0,0 +1,529 @@ + +package body Complex_Fixed_Points is + + -- I have no idea what visibility conflict is requiring these two to be defined + + function "*" + (Left, Right : in Real'Base) + return Real'Base is + begin + return Standard."*" (Left, Right); + end "*"; + + function "/" + (Left, Right : in Real'Base) + return Real'Base is + begin + return Standard."/" (Left, Right); + end "/"; + + + + function Re + (X : in Complex) + return Real'Base is + begin + return X.Re; + end Re; + + function Im + (X : in Complex) + return Real'Base is + begin + return X.Im; + end Im; + + function Im + (X : in Imaginary) + return Real'Base is + begin + return Real'Base (X); + end Im; + + + + procedure Set_Re + (X : in out Complex; + Re : in Real'Base) is + begin + X.Re := Re; + end Set_Re; + + procedure Set_Im + (X : in out Complex; + Im : in Real'Base) is + begin + X.Im := Im; + end Set_Im; + + procedure Set_Im + (X : in out Imaginary; + Im : in Real'Base) is + begin + X := Imaginary (Im); + end Set_Im; + + + + function Cartesian + (Re, Im : in Real'Base) + return Complex is + begin + return (Re => Re, Im => Im); + end Cartesian; + + function Cartesian + (Re : in Real'Base) + return Complex is + begin + return (Re => Re, Im => 0.0); + end Cartesian; + + function Cartesian + (Im : in Imaginary) + return Complex is + begin + return (Re => 0.0, Im => Real'Base (Im)); + end Cartesian; + + + + function Sqrt + (X : in Real'Base) + return Real'Base + is + Old_Result : Real'Base; + Result : Real'Base := 1.0; + Current_Square : Real'Base := 1.0; + Next_L : Real'Base := 3.0; + begin + -- If we don't check for this an overflow will happen due to taking things to the limit + if X = 0.0 then + return 0.0; + end if; + + -- Find the smallest integer greater than Sqrt (X) + while Current_Square < X loop + Current_Square := Current_Square + Next_L; + Next_L := Next_L + 2.0; + Result := Result + 1.0; + end loop; + + -- Lucky guess + if Current_Square = X then + return Result; + end if; + + -- Babylonian aka Newton's Method + -- To be honest I'm kinda suspicious this may not necessarily terminate + loop + Old_Result := Result; + Result := (Old_Result + X / Old_Result) / 2.0; + exit when Real (Result) = Real (Old_Result); + end loop; + return Result; + end Sqrt; + + function Modulus + (X : in Complex) + return Real'Base is + begin + return Sqrt (X.Re * X.Re + X.Im * X.Im); + end Modulus; + + + + function "+" + (Right : in Complex) + return Complex is + begin + return Right; + end "+"; + + function "-" + (Right : in Complex) + return Complex is + begin + return (Re => -Right.Re, Im => -Right.Im); + end "-"; + + function Conjugate + (X : in Complex) + return Complex is + begin + return (Re => X.Re, Im => -X.Im); + end Conjugate; + + + + function "+" + (Left, Right : in Complex) + return Complex is + begin + return (Re => Left.Re + Right.Re, + Im => Left.Im + Right.Im); + end "+"; + + function "-" + (Left, Right : in Complex) + return Complex is + begin + return (Re => Left.Re - Right.Re, + Im => Left.Im - Right.Im); + end "-"; + + function "*" + (Left, Right : in Complex) + return Complex is + begin + return (Re => Left.Re * Right.Re - Left.Im * Right.Im, + Im => Left.Re * Right.Im + Left.Im * Right.Re); + end "*"; + + function "/" + (Left, Right : in Complex) + return Complex + is + Denominator : Real'Base := Right.Re * Right.Re + Right.Im * Right.Im; + begin + return (Re => (Left.Re * Right.Re + Left.Im * Right.Im) / Denominator, + Im => (Left.Im * Right.Re - Left.Re * Right.Im) / Denominator); + end "/"; + + + + function "**" + (Left : in Complex; + Right : in Integer) + return Complex + is + Result : Complex; + Exponent : Integer := Right; + begin + if Exponent >= 1 then + Result := Left; + while Exponent > 1 loop + Result := Result * Left; + Exponent := Exponent - 1; + end loop; + else + Result := One; + while Exponent < 0 loop + Result := Result / Left; + Exponent := Exponent + 1; + end loop; + end if; + return Result; + end "**"; + + + + function "+" + (Right : in Imaginary) + return Imaginary is + begin + return Right; + end "+"; + + function "-" + (Right : in Imaginary) + return Imaginary is + begin + return Imaginary (-Real'Base (Right)); + end "-"; + + function "abs" + (Right : in Imaginary) + return Real'Base is + begin + return abs (Real'Base (Right)); + end "abs"; + + + + function "+" + (Left, Right : in Imaginary) + return Imaginary is + begin + return Imaginary (Real'Base (Left) + Real'Base (Right)); + end "+"; + + function "-" + (Left, Right : in Imaginary) + return Imaginary is + begin + return Imaginary (Real'Base (Left) - Real'Base (Right)); + end "-"; + + function "*" + (Left, Right : in Imaginary) + return Real'Base is + begin + return -(Real'Base (Left) * Real'Base (Right)); + end "*"; + + function "/" + (Left, Right : in Imaginary) + return Real'Base is + begin + return Real'Base (Left) / Real'Base (Right); + end "/"; + + + + function "**" + (Left : in Imaginary; + Right : in Integer) + return Complex + is + Result : Complex; + Exponent : Integer := Right; + begin + if Exponent >= 1 then + Result := Cartesian (Left); + while Exponent > 1 loop + Result := Result * Left; + Exponent := Exponent - 1; + end loop; + else + Result := One; + while Exponent < 0 loop + Result := Result / Left; + Exponent := Exponent + 1; + end loop; + end if; + return Result; + end "**"; + + + + function "<" + (Left, Right : in Imaginary) + return Boolean is + begin + return Real'Base (Left) < Real'Base (Right); + end "<"; + + function "<=" + (Left, Right : in Imaginary) + return Boolean is + begin + return Real'Base (Left) <= Real'Base (Right); + end "<="; + + function ">" + (Left, Right : in Imaginary) + return Boolean is + begin + return Real'Base (Left) > Real'Base (Right); + end ">"; + + function ">=" + (Left, Right : in Imaginary) + return Boolean is + begin + return Real'Base (Left) >= Real'Base (Right); + end ">="; + + + + function "+" + (Left : in Complex; + Right : in Real'Base) + return Complex is + begin + return (Re => Left.Re + Right, Im => Left.Im); + end "+"; + + function "+" + (Left : in Real'Base; + Right : in Complex) + return Complex is + begin + return Right + Left; + end "+"; + + function "-" + (Left : in Complex; + Right : in Real'Base) + return Complex is + begin + return (Re => Left.Re - Right, Im => Left.Im); + end "-"; + + function "-" + (Left : in Real'Base; + Right : in Complex) + return Complex is + begin + return (Re => Left - Right.Re, Im => -Right.Im); + end "-"; + + function "*" + (Left : in Complex; + Right : in Real'Base) + return Complex is + begin + return (Re => Left.Re * Right, Im => Left.Im * Right); + end "*"; + + function "*" + (Left : in Real'Base; + Right : in Complex) + return Complex is + begin + return Right * Left; + end "*"; + + function "/" + (Left : in Complex; + Right : in Real'Base) + return Complex is + begin + return (Re => Left.Re / Right, Im => Left.Im / Right); + end "/"; + + function "/" + (Left : in Real'Base; + Right : in Complex) + return Complex is + begin + return Cartesian (Left) / Right; + end "/"; + + + + function "+" + (Left : in Complex; + Right : in Imaginary) + return Complex is + begin + return (Re => Left.Re, Im => Left.Im + Real'Base (Right)); + end "+"; + + function "+" + (Left : in Imaginary; + Right : in Complex) + return Complex is + begin + return Right + Left; + end "+"; + + function "-" + (Left : in Complex; + Right : in Imaginary) + return Complex is + begin + return (Re => Left.Re, Im => Left.Im - Real'Base (Right)); + end "-"; + + function "-" + (Left : in Imaginary; + Right : in Complex) + return Complex is + begin + return (Re => -Right.Re, Im => Real'Base (Left) - Right.Im); + end "-"; + + function "*" + (Left : in Complex; + Right : in Imaginary) + return Complex is + begin + return (Re => -Left.Im * Real'Base (Right), + Im => Left.Re * Real'Base (Right)); + end "*"; + + function "*" + (Left : in Imaginary; + Right : in Complex) + return Complex is + begin + return Right * Left; + end "*"; + + function "/" + (Left : in Complex; + Right : in Imaginary) + return Complex is + begin + return (Re => Left.Im / Real'Base (Right), + Im => -Left.Re / Real'Base (Right)); + end "/"; + + function "/" + (Left : in Imaginary; + Right : in Complex) + return Complex is + begin + return Cartesian (Left) / Right; + end "/"; + + + + function "+" + (Left : in Imaginary; + Right : in Real'Base) + return Complex is + begin + return (Re => Right, Im => Real'Base (Left)); + end "+"; + + function "+" + (Left : in Real'Base; + Right : in Imaginary) + return Complex is + begin + return (Re => Left, Im => Real'Base (Right)); + end "+"; + + function "-" + (Left : in Imaginary; + Right : in Real'Base) + return Complex is + begin + return (Re => -Right, Im => Real'Base (Left)); + end "-"; + + function "-" + (Left : in Real'Base; + Right : in Imaginary) + return Complex is + begin + return (Re => Left, Im => -Real'Base (Right)); + end "-"; + + function "*" + (Left : in Imaginary; + Right : in Real'Base) + return Imaginary is + begin + return Imaginary (Standard."*" (Right, Real'Base (Left))); + end "*"; + + function "*" + (Left : in Real'Base; + Right : in Imaginary) + return Imaginary is + begin + return Imaginary (Standard."*" (Left, Real'Base (Right))); + end "*"; + + function "/" + (Left : in Imaginary; + Right : in Real'Base) + return Imaginary is + begin + return Imaginary (Standard."/" (Real'Base (Left), Right)); + end "/"; + + function "/" + (Left : in Real'Base; + Right : in Imaginary) + return Imaginary is + begin + return -Imaginary (Standard."/" (Left, Real'Base (Right))); + end "/"; + +end Complex_Fixed_Points; + |