summaryrefslogtreecommitdiff
path: root/src/complex_fixed_points.adb
diff options
context:
space:
mode:
Diffstat (limited to 'src/complex_fixed_points.adb')
-rw-r--r--src/complex_fixed_points.adb529
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;
+
generated by cgit (git 2.30.0) at 2024-05-15 12:53:52 +0000