From 452a2361a43b26089b1ba755f4224935b5b5e033 Mon Sep 17 00:00:00 2001
From: Jedidiah Barber <contact@jedbarber.id.au>
Date: Fri, 18 Nov 2022 04:09:42 +1300
Subject: Initial commit

---
 src/complex_fixed_points.adb | 529 +++++++++++++++++++++++++++++++++++++++++++
 src/complex_fixed_points.ads | 117 ++++++++++
 2 files changed, 646 insertions(+)
 create mode 100644 src/complex_fixed_points.adb
 create mode 100644 src/complex_fixed_points.ads

(limited to 'src')

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;
+
diff --git a/src/complex_fixed_points.ads b/src/complex_fixed_points.ads
new file mode 100644
index 0000000..97667b1
--- /dev/null
+++ b/src/complex_fixed_points.ads
@@ -0,0 +1,117 @@
+
+generic
+    type Real is delta <> digits <>;
+package Complex_Fixed_Points is
+
+    pragma Pure (Complex_Fixed_Points);
+
+    type Complex is private;
+    type Imaginary is private;
+    pragma Preelaborable_Initialization (Imaginary);
+
+    i : constant Imaginary;
+    j : constant Imaginary;
+
+    Zero : constant Complex;
+    One  : constant Complex;
+
+    function Re (X : in Complex) return Real'Base;
+    function Im (X : in Complex) return Real'Base;
+    function Im (X : in Imaginary) return Real'Base;
+
+    procedure Set_Re
+           (X  : in out Complex;
+            Re : in     Real'Base);
+    procedure Set_Im
+           (X  : in out Complex;
+            Im : in     Real'Base);
+    procedure Set_Im
+           (X  : in out Imaginary;
+            Im : in     Real'Base);
+
+    function Cartesian (Re, Im : in Real'Base) return Complex;
+    function Cartesian (Re     : in Real'Base) return Complex;
+    function Cartesian (Im     : in Imaginary) return Complex;
+
+    function Modulus (X     : in Complex) return Real'Base;
+    function "abs"   (Right : in Complex) return Real'Base renames Modulus;
+
+    --  Like hell am I writing fixed point trig functions right now.
+
+    --  function Argument (X : in Complex) return Real'Base;
+    --  function Argument (X : in Complex; Cycle : in Real'Base) return Real'Base;
+
+    --  function Polar (Modulus, Argument        : in Real'Base) return Complex;
+    --  function Polar (Modulus, Argument, Cycle : in Real'Base) return Complex;
+
+    function "+" (Right : in Complex) return Complex;
+    function "-" (Right : in Complex) return Complex;
+    function Conjugate (X : in Complex) return Complex;
+
+    function "+" (Left, Right : in Complex) return Complex;
+    function "-" (Left, Right : in Complex) return Complex;
+    function "*" (Left, Right : in Complex) return Complex;
+    function "/" (Left, Right : in Complex) return Complex;
+
+    function "**" (Left : in Complex; Right : in Integer) return Complex;
+
+    function "+"       (Right : in Imaginary) return Imaginary;
+    function "-"       (Right : in Imaginary) return Imaginary;
+    function Conjugate (X     : in Imaginary) return Imaginary renames "-";
+    function "abs"     (Right : in Imaginary) return Real'Base;
+
+    function "+" (Left, Right : in Imaginary) return Imaginary;
+    function "-" (Left, Right : in Imaginary) return Imaginary;
+    function "*" (Left, Right : in Imaginary) return Real'Base;
+    function "/" (Left, Right : in Imaginary) return Real'Base;
+
+    function "**" (Left : in Imaginary; Right : in Integer) return Complex;
+
+    function "<"  (Left, Right : in Imaginary) return Boolean;
+    function "<=" (Left, Right : in Imaginary) return Boolean;
+    function ">"  (Left, Right : in Imaginary) return Boolean;
+    function ">=" (Left, Right : in Imaginary) return Boolean;
+
+    function "+" (Left : in Complex;   Right : in Real'Base) return Complex;
+    function "+" (Left : in Real'Base; Right : in Complex)   return Complex;
+    function "-" (Left : in Complex;   Right : in Real'Base) return Complex;
+    function "-" (Left : in Real'Base; Right : in Complex)   return Complex;
+    function "*" (Left : in Complex;   Right : in Real'Base) return Complex;
+    function "*" (Left : in Real'Base; Right : in Complex)   return Complex;
+    function "/" (Left : in Complex;   Right : in Real'Base) return Complex;
+    function "/" (Left : in Real'Base; Right : in Complex)   return Complex;
+
+    function "+" (Left : in Complex;   Right : in Imaginary) return Complex;
+    function "+" (Left : in Imaginary; Right : in Complex)   return Complex;
+    function "-" (Left : in Complex;   Right : in Imaginary) return Complex;
+    function "-" (Left : in Imaginary; Right : in Complex)   return Complex;
+    function "*" (Left : in Complex;   Right : in Imaginary) return Complex;
+    function "*" (Left : in Imaginary; Right : in Complex)   return Complex;
+    function "/" (Left : in Complex;   Right : in Imaginary) return Complex;
+    function "/" (Left : in Imaginary; Right : in Complex)   return Complex;
+
+    function "+" (Left : in Imaginary; Right : in Real'Base) return Complex;
+    function "+" (Left : in Real'Base; Right : in Imaginary) return Complex;
+    function "-" (Left : in Imaginary; Right : in Real'Base) return Complex;
+    function "-" (Left : in Real'Base; Right : in Imaginary) return Complex;
+    function "*" (Left : in Imaginary; Right : in Real'Base) return Imaginary;
+    function "*" (Left : in Real'Base; Right : in Imaginary) return Imaginary;
+    function "/" (Left : in Imaginary; Right : in Real'Base) return Imaginary;
+    function "/" (Left : in Real'Base; Right : in Imaginary) return Imaginary;
+
+private
+
+    type Complex is record
+        Re, Im : Real'Base;
+    end record;
+
+    type Imaginary is new Real'Base;
+
+    i : constant Imaginary := 1.0;
+    j : constant Imaginary := 1.0;
+
+    Zero : constant Complex := (Re => 0.0, Im => 0.0);
+    One  : constant Complex := (Re => 1.0, Im => 0.0);
+
+end Complex_Fixed_Points;
+
-- 
cgit