orka_b455160b/orka_tensors_cpu/src/orka-numerics-tensors-simd_cpu.adb

--  SPDX-License-Identifier: Apache-2.0
--
--  Copyright (c) 2021 onox <denkpadje@gmail.com>
--
--  Licensed under the Apache License, Version 2.0 (the "License");
--  you may not use this file except in compliance with the License.
--  You may obtain a copy of the License at
--
--      http://www.apache.org/licenses/LICENSE-2.0
--
--  Unless required by applicable law or agreed to in writing, software
--  distributed under the License is distributed on an "AS IS" BASIS,
--  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
--  See the License for the specific language governing permissions and
--  limitations under the License.

with Ada.Characters.Latin_1;
with Ada.Strings.Unbounded;
with Ada.Unchecked_Conversion;

with Orka.Numerics.Tensors.Operations;

package body Orka.Numerics.Tensors.SIMD_CPU is

   procedure Swap_Rows (Ab : in out CPU_Tensor; I, J : Index_Type) is
   begin
      if I /= J then
         declare
            Row_I : constant CPU_Tensor := Ab (I);
            Old_J : constant CPU_Tensor := Ab (J);
         begin
            Set (Ab, J, Row_I);
            Set (Ab, I, Old_J);
         end;
      end if;
   end Swap_Rows;

   procedure Scale_Row (Ab : in out CPU_Tensor; I : Index_Type; Scale : Element) is
      Row_I : constant CPU_Tensor := Ab (I);
   begin
      if Scale /= 1.0 then
         Set (Ab, I, Scale * Row_I);
      end if;
   end Scale_Row;

   procedure Replace_Row (Ab : in out CPU_Tensor; Scale : Element; I, J : Index_Type) is
      Row_I : constant CPU_Tensor := Ab (I);
      Row_J : constant CPU_Tensor := Ab (J);
   begin
      if Scale /= 0.0 then
         Set (Ab, J, Row_J - Scale * Row_I);
      end if;
   end Replace_Row;

   procedure Forward_Substitute (Ab : in out CPU_Tensor; Index, Pivot_Index : Index_Type) is
      Rows        : constant Natural := Ab.Rows;
      Pivot_Value : constant Element := Ab ((Index, Pivot_Index));
   begin
      --  Create zeros below the pivot position
      for Row_Index in Index + 1 .. Rows loop
         Replace_Row (Ab, Ab ((Row_Index, Pivot_Index)) / Pivot_Value, Index, Row_Index);
      end loop;
   end Forward_Substitute;

   procedure Back_Substitute (Ab : in out CPU_Tensor; Index, Pivot_Index : Index_Type) is
   begin
      Scale_Row (Ab, Index, 1.0 / Ab ((Index, Pivot_Index)));

      --  Create zeros above the pivot position
      for Row_Index in 1 .. Index - 1 loop
         Replace_Row (Ab, Ab ((Row_Index, Pivot_Index)), Index, Row_Index);
      end loop;
   end Back_Substitute;

   function Create_QR
     (Q, R         : CPU_Tensor;
      Determinancy : Matrix_Determinancy) return CPU_QR_Factorization
   is (Q_Size       => Q.Size,
       R_Size       => R.Size,
       Q            => Q,
       R            => R,
       Determinancy => Determinancy);

   procedure Make_Upper_Triangular (Object : in out CPU_Tensor; Offset : Integer := 0) is
      Rows    : constant Natural := Object.Rows;
      Columns : constant Natural := Object.Columns;
   begin
      if Offset >= -(Rows - 2) then
         --  Make matrix upper triangular by zeroing out the elements in the
         --  lower triangular part
         for Row_Index in Index_Type'First + 1 - Integer'Min (1, Offset) .. Rows loop
            for Column_Index in 1 .. Natural'Min (Row_Index - 1 + Offset, Columns) loop
               Object.Set ((Row_Index, Column_Index), 0.0);
            end loop;
         end loop;
      end if;
   end Make_Upper_Triangular;

   package Operations is new Orka.Numerics.Tensors.Operations
     (CPU_Tensor, Make_Upper_Triangular, Scale_Row, Swap_Rows, Forward_Substitute, Back_Substitute,
      Expression_Type, CPU_QR_Factorization, Create_QR, Q, R);

   ----------------------------------------------------------------------------

   function Convert is new Ada.Unchecked_Conversion (Vector_Type, Integer_Vector_Type);
   function Convert is new Ada.Unchecked_Conversion (Integer_Vector_Type, Vector_Type);

   function From_Last (Offset : Natural) return Vector_Index_Type is
     (Vector_Index_Type'Val (Vector_Index_Type'Pos (Vector_Index_Type'Last) - Offset));

   function Data_Vectors (Index : Natural) return Natural is
     ((Index + (Vector_Type'Length - 1)) / Vector_Type'Length);

   function Data_Vectors (Shape : Tensor_Shape) return Natural is
     (Data_Vectors (Elements (Shape)));

   function Data_Offset (Index : Positive) return Vector_Index_Type is
     (Vector_Index_Type'Val ((Index - 1) mod Vector_Type'Length
        + Vector_Index_Type'Pos (Vector_Index_Type'First)));

   function Data_Padding (Size, Count : Natural) return Natural is
     (Size * Vector_Type'Length - Count);

   function To_Index (Index : Tensor_Index; Columns : Positive) return Index_Type is
     (To_Index (Index, (Index (1), Columns)))
   with Pre => Index'Length = 2;
   --  Index (1) isn't used for 2-D Index except in Pre condition of To_Index in parent package

   function Reset_Padding
     (Object  : CPU_Tensor;
      Padding : Natural;
      Value   : Element_Type) return Vector_Type
   is
      Last_Vector : Vector_Type := Object.Data (Object.Data'Last);
   begin
      if Padding /= 0 then
         Last_Vector (From_Last (Padding - 1) .. Last_Vector'Last) := (others => Value);
      end if;

      return Last_Vector;
   end Reset_Padding;

   function Disable_Padding (Object : CPU_Tensor) return Vector_Type is
      Padding : constant Natural :=
        Data_Padding (Size => Object.Size, Count => Object.Elements);
   begin
      return Reset_Padding (Object, Padding, 0.0);
   end Disable_Padding;

   function Without_Data
     (Object : CPU_Tensor;
      Kind   : Data_Type := Float_Type) return CPU_Tensor;
   --  This silly definition is needed to avoid "length check failed" in GNAT FSF 11.1

   function Without_Data
     (Object : CPU_Tensor;
      Kind   : Data_Type := Float_Type) return CPU_Tensor
   is
     ((Axes  => Object.Axes,
       Size  => Object.Size,
       Kind  => Kind,
       Shape => Object.Shape,
       Data  => <>));

   function Without_Data
     (Shape : Tensor_Shape;
      Kind  : Data_Type := Float_Type) return CPU_Tensor
   is
     ((Axes  => Shape'Length,
       Size  => Data_Vectors (Shape),
       Kind  => Kind,
       Shape => Shape,
       Data  => <>));

   ----------------------------------------------------------------------------
   --                              Expressions                               --
   ----------------------------------------------------------------------------

   function Identity_Vector_Type (Value : Element) return Vector_Type is (others => Value);
   function Identity_Element     (Value : Element) return Element     is (Value);

   function Apply is new Generic_Apply (Vector_Type, Identity_Vector_Type);

   function Apply is new Generic_Apply (Element, Identity_Element,
     Min => Element'Min, Max => Element'Max, Sqrt => EF.Sqrt);

   ----------------------------------------------------------------------------

   overriding function Kind (Object : CPU_Tensor) return Data_Type is (Object.Kind);

   overriding
   function Get (Object : CPU_Tensor; Index : Index_Type) return Element renames Operations.Get;

   overriding
   function Get (Object : CPU_Tensor; Index : Index_Type) return Boolean renames Operations.Get;

   overriding function Get (Object : CPU_Tensor; Index : Index_Type) return CPU_Tensor is
      Count : constant Positive := Object.Columns;
      Shape : constant Tensor_Shape := (1 => Count);
   begin
      if Index > Object.Rows then
         raise Constraint_Error with
           "Stop index (" & Trim (Index) & ") out of bounds (1 .. " & Trim (Object.Rows) & ")";
      end if;

      --  Returning the row of a 2D tensor as a vector instead of a (1, n) 2D tensor
      return Result : CPU_Tensor := Without_Data (Shape, Object.Kind) do
         declare
            Result_Data : Element_Array (1 .. Count)
              with Import, Convention => Ada, Address => Result.Data'Address;

            Object_Data : Element_Array (1 .. Object.Elements)
              with Import, Convention => Ada, Address => Object.Data'Address;

            Base_Index : constant Natural := (Index - 1) * Count;
         begin
            Result_Data (1 .. Count) := Object_Data (Base_Index + 1 .. Base_Index + Count);
         end;
      end return;
   end Get;

   overriding procedure Set
     (Object : in out CPU_Tensor;
      Index  : Index_Type;
      Value  : CPU_Tensor) renames Operations.Set;

   overriding procedure Set
     (Object : in out CPU_Tensor;
      Index  : Range_Type;
      Value  : CPU_Tensor) renames Operations.Set;

   overriding
   procedure Set (Object : in out CPU_Tensor; Index : Tensor_Range; Value : CPU_Tensor) is
      Full_Index : constant Tensor_Range := Full_Range (Object.Shape, Index);
      Full_Value : constant Tensor_Shape := Full_Shape (Object.Axes, Value.Shape, Right);

      pragma Assert (Full_Value = Shape (Full_Index));
   begin
      --  If the value (and shape of index) has the full depth/height/width except
      --  for the first axis, then the memory to which the data will be written
      --  is contiguous, which means it has no gaps.
      --
      --  For example, if shape of Value is (2, 3) and you have the following
      --  object and index (in brackets):
      --
      --  1 [ 2  3  4]  5
      --  6 [ 7  8  9] 10
      --  11 12 13 14  15
      --
      --  then there is a gap (positions 5 and 6). Howerver, if the shape
      --  of Value is (2, 5) (with a matching Index) then there are no gaps.
      --
      --  Another case in which there are are no gaps is when all but the last
      --  axis have a shape equal to 1. For example if the index is
      --  ((2, 2), (7, 9)), which has the shape (1, 3).
      if Is_Equal (Object.Shape, Full_Value, 1)
        or else (for all D in Full_Value'First .. Full_Value'Last - 1 => Full_Value (D) = 1)
      then
         declare
            Start_Index : Tensor_Index (Full_Index'Range);
            Stop_Index  : Tensor_Index (Full_Index'Range);
         begin
            for Axis in Full_Index'Range loop
               Start_Index (Axis) := Full_Index (Axis).Start;
               Stop_Index (Axis)  := Full_Index (Axis).Stop;
            end loop;

            declare
               Start_Index_Flattened : constant Index_Type := To_Index (Start_Index, Object.Shape);
               Stop_Index_Flattened  : constant Index_Type := To_Index (Stop_Index, Object.Shape);

               Count : constant Natural := Value.Elements;
               pragma Assert (Stop_Index_Flattened - Start_Index_Flattened + 1 = Count);

               Row_Data : Element_Array (1 .. Count)
                 with Import, Convention => Ada, Address => Value.Data'Address;

               Object_Data : Element_Array (1 .. Object.Elements)
                 with Import, Convention => Ada, Address => Object.Data'Address;
            begin
               Object_Data (Start_Index_Flattened .. Stop_Index_Flattened) := Row_Data;
            end;
         end;
      else
         raise Program_Error with "Not implemented yet";  --  FIXME
      end if;
   end Set;

   function Flattened_Index (Object : CPU_Tensor; Index : Tensor_Index) return Index_Type is
   begin
      for Axis in Index'Range loop
         declare
            Index_Dim : constant Natural := Index (Axis);
            Shape_Dim : constant Natural := Object.Shape (Axis);
         begin
            if Index_Dim > Shape_Dim then
               raise Constraint_Error with
                 "Index (" & Trim (Index_Dim) & ") out of bounds (1 .. " & Trim (Shape_Dim) & ")";
            end if;
         end;
      end loop;

      return To_Index (Index, Object.Shape);
   end Flattened_Index;

   overriding procedure Set
     (Object : in out CPU_Tensor;
      Index  : Index_Type;
      Value  : Element) renames Operations.Set;

   overriding procedure Set
     (Object : in out CPU_Tensor;
      Index  : Index_Type;
      Value  : Boolean) renames Operations.Set;

   overriding procedure Set (Object : in out CPU_Tensor; Index : Tensor_Index; Value : Element) is
      Index_Flattened : constant Index_Type := Flattened_Index (Object, Index);
   begin
      Object.Data (Data_Vectors (Index_Flattened)) (Data_Offset (Index_Flattened)) := Value;
   end Set;

   overriding procedure Set (Object : in out CPU_Tensor; Index : Tensor_Index; Value : Boolean) is
      Index_Flattened : constant Index_Type := Flattened_Index (Object, Index);

      Zero_Vector : constant Vector_Type := (others => 0.0);

      Mask : constant Integer_Vector_Type := Convert (Zero_Vector = Zero_Vector);

      Booleans : Integer_Vector_Type := Convert (Object.Data (Data_Vectors (Index_Flattened)));
   begin
      Booleans (Data_Offset (Index_Flattened)) := (if Value then Mask (Mask'First) else 0);
      Object.Data (Data_Vectors (Index_Flattened)) := Convert (Booleans);
   end Set;

   overriding function Get (Object : CPU_Tensor; Index : Tensor_Index) return Element is
      Object_Data : Element_Array (1 .. Object.Elements)
        with Import, Convention => Ada, Address => Object.Data'Address;
   begin
      return Object_Data (Flattened_Index (Object, Index));
   end Get;

   overriding function Get (Object : CPU_Tensor; Index : Tensor_Index) return Boolean is
      Index_Flattened : constant Index_Type := Flattened_Index (Object, Index);

      One_Vector : constant Integer_Vector_Type := (others => 1);

      Mask       : constant Integer_Vector_Type :=
        Convert (Object.Data (Data_Vectors (Index_Flattened)));
      Ones_Zeros : constant Integer_Vector_Type := Mask and One_Vector;
   begin
      return Ones_Zeros (Data_Offset (Index_Flattened)) = 1;
   end Get;

   overriding
   function Get (Object : CPU_Tensor; Index : Range_Type) return CPU_Tensor renames Operations.Get;

   overriding function Get (Object : CPU_Tensor; Index : Tensor_Range) return CPU_Tensor is
      Rows : constant Natural := Object.Rows;

      Row_Start : constant Index_Type := Index (1).Start;
      Row_Stop  : constant Index_Type := Index (1).Stop;

      Result_Rows : constant Positive := Row_Stop - Row_Start + 1;
   begin
      case Object.Axes is
         when 1 =>
            declare
               Count : constant Positive := Result_Rows;
               Size  : constant Positive := Data_Vectors (Count);
               Shape : constant Tensor_Shape := (1 => Count);
            begin
               if Row_Stop > Rows then
                  raise Constraint_Error with
                    "Stop index (" & Trim (Row_Stop) & ") out of bounds (1 .. " &
                    Trim (Rows) & ")";
               end if;

               return Result : CPU_Tensor := Without_Data (Shape, Object.Kind) do
                  if Data_Offset (Row_Start) = Vector_Index_Type'First then
                     Result.Data (1 .. Size) :=
                       Object.Data (Data_Vectors (Row_Start) .. Data_Vectors (Row_Stop));
                  else
                     declare
                        Result_Data : Element_Array (1 .. Count)
                          with Import, Convention => Ada, Address => Result.Data'Address;

                        Object_Data : Element_Array (1 .. Object.Elements)
                          with Import, Convention => Ada, Address => Object.Data'Address;
                     begin
                        Result_Data (1 .. Count) := Object_Data (Row_Start .. Row_Stop);
                     end;
                  end if;
               end return;
            end;
         when 2 =>
            declare
               Columns : constant Natural :=
                 (if 2 in Object.Shape'Range then Object.Columns else 1);

               Index_Shape : constant Tensor_Shape := Shape (Index);
               Result_Columns : constant Positive :=
                 (if 2 in Index_Shape'Range then Index_Shape (2) else Columns);

               Shape : constant Tensor_Shape :=
                 (if Result_Rows = 1 then
                    (1 => Result_Columns)
                  else
                    (1 => Result_Rows, 2 => Result_Columns));

               Column_Start : constant Index_Type :=
                 (if 2 in Index'Range then Index (2).Start else 1);
               Column_Stop  : constant Index_Type :=
                 (if 2 in Index'Range then Index (2).Stop else Columns);
            begin
               if Row_Stop > Rows then
                  raise Constraint_Error with
                    "Stop index (" & Trim (Row_Stop) & ") out of bounds (1 .. " &
                    Trim (Rows) & ")";
               end if;

               if Column_Stop > Columns then
                  raise Constraint_Error with
                    "Stop index (" & Trim (Column_Stop) & ") out of bounds (1 .. " &
                    Trim (Columns) & ")";
               end if;

               return Result : CPU_Tensor := Without_Data (Shape, Object.Kind) do
                  declare
                     Result_Data : Element_Array (1 .. Result.Elements)
                       with Import, Convention => Ada, Address => Result.Data'Address;

                     Object_Data : Element_Array (1 .. Object.Elements)
                       with Import, Convention => Ada, Address => Object.Data'Address;
                  begin
                     for I in 1 .. Result_Rows loop
                        declare
                           Result_Index : constant Natural := (I - 1) * Result_Columns;

                           Current_Row : constant Natural := Row_Start - 1 + I;
                           Base_Index  : constant Natural := (Current_Row - 1) * Columns;
                        begin
                           Result_Data (Result_Index + 1 .. Result_Index + Result_Columns) :=
                             Object_Data (Base_Index + Column_Start .. Base_Index + Column_Stop);
                        end;
                     end loop;
                  end;
               end return;
            end;
         when others =>
            raise Not_Implemented_Yet;  --  FIXME
      end case;
   end Get;

   overriding function Get (Object : CPU_Tensor; Index : CPU_Tensor) return CPU_Tensor is
      type Integer_Vector_Array is array (Index_Type range <>) of Integer_Vector_Type;

      One_Vector : constant Integer_Vector_Type := (others => 1);

      Indices : Integer_Vector_Array (1 .. Index.Size);

      Offset : Integer_Vector_Type := (others => 0);
   begin
      for Index_Vector in Index.Data'Range loop
         declare
            --  1. Create Integer_Vector of 1s and 0s
            Mask : constant Integer_Vector_Type := Convert (Index.Data (Index_Vector));

            PS  : Integer_Vector_Type := Mask and One_Vector;
            Sum : Integer_Vector_Type := Offset + PS;
         begin
            --  2. Compute prefix sum
            for I in 2 .. Vector_Type'Length loop
               PS  := Shift_Elements_Left (PS);
               Sum := Sum + PS;
            end loop;

            --  Index is > 0 if valid according to mask
            Indices (Index_Vector) := Mask and Sum;

            --  The last number (sum of all the elements in the vector)
            --  is the offset for the next vector
            Offset := (others => Sum (Sum'Last));
         end;
      end loop;

      --  3. Last number of last vector is number of elements in Result tensor
      return Result : CPU_Tensor :=
        Without_Data (Tensor_Shape'(1 => Natural (Offset (Offset'Last))))
      do
         if Result.Elements > 0 then
            declare
               Result_Data : Element_Array (1 .. Result.Elements)
                 with Import, Convention => Ada, Address => Result.Data'Address;
            begin
               --  4. Iterate over the prefix sum and assign the value from Object.Data
               --  to Result.Data at the index found in the prefix sum
               for Index_Vector in Index.Data'Range loop
                  declare
                     Indices_Vector : Integer_Vector_Type renames Indices (Index_Vector);
                  begin
                     for I in Indices_Vector'Range loop
                        if Indices_Vector (I) > 0 then
                           Result_Data (Natural (Indices_Vector (I))) :=
                             Object.Data (Index_Vector) (I);
                        end if;
                     end loop;
                  end;
               end loop;
            end;
         end if;
      end return;
   end Get;

   ----------------------------------------------------------------------------

   overriding
   function Image (Object : CPU_Tensor) return String is
      package L1 renames Ada.Characters.Latin_1;
      package SU renames Ada.Strings.Unbounded;

      Row_Count : constant := 5;
      Count     : constant Natural := Object.Elements;

      Result : SU.Unbounded_String;
   begin
      SU.Append (Result, "tensor([");
      case Object.Axes is
         when 1 =>
            for I in 1 .. Count loop
               declare
                  First_Element_Of_Row : constant Boolean := (I - 1) mod Row_Count = 0;
                  Last_Element_Of_Row  : constant Boolean := (I - 0) mod Row_Count = 0;

                  Value : Element_Type renames Object.Data (Data_Vectors (I)) (Data_Offset (I));
               begin
                  if First_Element_Of_Row then
                     SU.Append (Result, (if I = 1 then "" else "        "));
                  end if;
                  case Object.Kind is
                     when Float_Type | Int_Type =>
                        SU.Append (Result, (if Value'Valid then Value'Image else "     invalid"));
                     when Bool_Type =>
                        SU.Append (Result, "       " &
                          (if not Value'Valid or else Value /= 0.0 then " True" else "False"));
                  end case;
                  if I < Count then
                     SU.Append (Result, ",");
                     if Last_Element_Of_Row then
                        SU.Append (Result, L1.LF);
                     end if;
                  end if;
               end;
            end loop;
         when 2 =>
            declare
               Rows    : constant Natural := Object.Rows;
               Columns : constant Natural := Object.Columns;

               Object_Data : Element_Array (1 .. Count)
                 with Import, Convention => Ada, Address => Object.Data'Address;
            begin
               for I in 1 .. Rows loop
                  SU.Append (Result, (if I = 1 then "" else "        "));
                  SU.Append (Result, "[");
                  for J in 1 .. Columns loop
                     declare
                        Value : Element_Type renames Object_Data ((I - 1) * Columns + J);
                     begin
                        case Object.Kind is
                           when Float_Type | Int_Type =>
                              SU.Append (Result,
                                (if Value'Valid then Value'Image else "     invalid"));
                           when Bool_Type =>
                              SU.Append (Result, "       " &
                                (if not Value'Valid or else Value /= 0.0 then
                                   " True"
                                 else
                                   "False"));
                        end case;
                     end;
                     if J < Columns then
                        SU.Append (Result, ",");
                     end if;
                  end loop;
                  SU.Append (Result, "]");
                  if I < Rows then
                     SU.Append (Result, ",");
                     SU.Append (Result, L1.LF);
                  end if;
               end loop;
            end;
         when others =>
            raise Not_Implemented_Yet;  --  FIXME
      end case;
      SU.Append (Result, "])");

      return SU.To_String (Result);
   end Image;

   overriding
   function Shape (Object : CPU_Tensor) return Tensor_Shape is (Object.Shape);

   overriding
   function Elements (Object : CPU_Tensor) return Natural is (Elements (Object.Shape));

   overriding
   function Axes (Object : CPU_Tensor) return Tensor_Axis is (Object.Axes);

   overriding
   function Empty (Shape : Tensor_Shape) return CPU_Tensor is (Without_Data (Shape));

   overriding
   function Fill (Shape : Tensor_Shape; Value : Element) return CPU_Tensor is
      Vector : constant Vector_Type := (others => Value);
   begin
      return Result : CPU_Tensor := Without_Data (Shape) do
         Result.Data := (others => Vector);
      end return;
   end Fill;

   overriding function Zeros (Shape : Tensor_Shape) return CPU_Tensor renames Operations.Zeros;
   overriding function Zeros (Elements : Positive)  return CPU_Tensor renames Operations.Zeros;

   overriding function Ones (Shape : Tensor_Shape) return CPU_Tensor renames Operations.Ones;
   overriding function Ones (Elements : Positive)  return CPU_Tensor renames Operations.Ones;

   overriding
   function To_Tensor (Elements : Element_Array; Shape : Tensor_Shape) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Shape) do
         declare
            Result_Data : Element_Array (1 .. Result.Elements)
              with Import, Convention => Ada, Address => Result.Data'Address;
         begin
            Result_Data := Elements;
         end;
      end return;
   end To_Tensor;

   overriding
   function To_Tensor (Elements : Element_Array) return CPU_Tensor renames Operations.To_Tensor;

   overriding
   function To_Boolean_Tensor
     (Booleans : Boolean_Array;
      Shape    : Tensor_Shape) return CPU_Tensor
   is
      Elements : Element_Array (Booleans'Range);
   begin
      for Index in Elements'Range loop
         Elements (Index) := (if Booleans (Index) then 1.0 else 0.0);
      end loop;

      return To_Tensor (Elements, Shape) > 0.0;
   end To_Boolean_Tensor;

   overriding
   function To_Boolean_Tensor (Booleans : Boolean_Array) return CPU_Tensor
     renames Operations.To_Boolean_Tensor;

   overriding
   function Linear_Space
     (Start, Stop : Element;
      Count       : Positive;
      Interval    : Interval_Kind := Closed) return CPU_Tensor
   is
      Shape : constant Tensor_Shape := (1 => Count);

      Step : constant Element :=
        (if Count > 1 then (Stop - Start) / Element (Count - (case Interval is
                                                                when Closed    => 1,
                                                                when Half_Open => 0))
                      else 0.0);
   begin
      return Result : CPU_Tensor := Without_Data (Shape) do
         declare
            Result_Data : Element_Array (1 .. Elements (Shape))
              with Import, Convention => Ada, Address => Result.Data'Address;
         begin
            for Index in 1 .. Count loop
               Result_Data (Index) := Start + Step * Element (Index - 1);
            end loop;
         end;
      end return;
   end Linear_Space;

   overriding
   function Log_Space
     (Start, Stop : Element;
      Count       : Positive;
      Interval    : Interval_Kind := Closed;
      Base        : Element := 10.0) return CPU_Tensor
   is
      Shape : constant Tensor_Shape := (1 => Count);

      Step : constant Element :=
        (if Count > 1 then (Stop - Start) / Element (Count - (case Interval is
                                                                when Closed    => 1,
                                                                when Half_Open => 0))
                      else 0.0);
   begin
      return Result : CPU_Tensor := Without_Data (Shape) do
         declare
            Result_Data : Element_Array (1 .. Elements (Shape))
              with Import, Convention => Ada, Address => Result.Data'Address;

            use EF;
         begin
            for Index in 1 .. Count loop
               Result_Data (Index) := Base ** (Start + Step * Element (Index - 1));
            end loop;
         end;
      end return;
   end Log_Space;

   overriding
   function Geometric_Space
     (Start, Stop : Element;
      Count       : Positive;
      Interval    : Interval_Kind := Closed;
      Base        : Element := 10.0) return CPU_Tensor renames Operations.Geometric_Space;

   overriding
   function Array_Range (Start, Stop : Element; Step : Element := 1.0) return CPU_Tensor
     renames Operations.Array_Range;

   overriding
   function Array_Range (Stop : Element) return CPU_Tensor renames Operations.Array_Range;

   overriding
   function Identity (Size : Positive; Offset : Integer := 0) return CPU_Tensor
     renames Operations.Identity;

   overriding
   function Identity (Rows, Columns : Positive; Offset : Integer := 0) return CPU_Tensor is
      Shape : constant Tensor_Shape := (1 => Rows, 2 => Columns);

      Max_Size : constant Positive := Positive'Max (Rows, Columns);

      Zero_Vector : constant Vector_Type := (others => 0.0);
   begin
      return Result : CPU_Tensor :=
        (Axes => Shape'Length,
         Size       => Data_Vectors (Shape),
         Kind       => Float_Type,
         Shape      => Shape,
         Data       => (others => Zero_Vector))
      do
         if Offset in -(Max_Size - 1) .. Max_Size - 1 then
            declare
               Result_Data : Element_Array (1 .. Elements (Shape))
                 with Import, Convention => Ada, Address => Result.Data'Address;

               Index : Integer := Offset + 1;
            begin
               for I in 1 .. Rows loop
                  if Index in Result_Data'Range then
                     Result_Data (Index) := 1.0;
                  end if;
                  Index := Index + Columns + 1;
               end loop;
            end;
         end if;
      end return;
   end Identity;

   overriding
   function Upper_Triangular (Object : CPU_Tensor; Offset : Integer := 0) return CPU_Tensor
     renames Operations.Upper_Triangular;

   overriding
   function Main_Diagonal (Object : CPU_Tensor; Offset : Integer := 0) return CPU_Tensor is
      Rows    : constant Positive := Object.Rows;
      Columns : constant Positive := Object.Columns;

      Shape : constant Tensor_Shape := (1 => Positive'Min (Rows, Columns));
   begin
      return Result : CPU_Tensor := Without_Data (Shape) do
         declare
            Result_Data : Element_Array (1 .. Elements (Shape))
              with Import, Convention => Ada, Address => Result.Data'Address;

            Object_Data : Element_Array (1 .. Object.Elements)
              with Import, Convention => Ada, Address => Object.Data'Address;

            Index : Integer := Offset + 1;
         begin
            for I in Result_Data'Range loop
               if Index in Object_Data'Range then
                  Result_Data (I) := Object_Data (Index);
               else
                  Result_Data (I) := 0.0;
               end if;
               Index := Index + Columns + 1;
            end loop;
         end;
      end return;
   end Main_Diagonal;

   overriding
   function Diagonal (Elements : Element_Array; Offset : Integer := 0) return CPU_Tensor is
      Size : constant Positive := Elements'Length;

      Shape : constant Tensor_Shape := (1 .. 2 => Size);

      Zero_Vector : constant Vector_Type := (others => 0.0);
   begin
      return Result : CPU_Tensor :=
        (Axes => Shape'Length,
         Size       => Data_Vectors (Shape),
         Kind       => Float_Type,
         Shape      => Shape,
         Data       => (others => Zero_Vector))
      do
         if Offset in -(Size - 1) .. Size - 1 then
            declare
               Result_Data : Element_Array (1 .. Size * Size)
                 with Import, Convention => Ada, Address => Result.Data'Address;

               Index : Integer := Offset + 1;
            begin
               for I in Elements'Range loop
                  if Index in Result_Data'Range then
                     Result_Data (Index) := Elements (I);
                  end if;
                  Index := Index + Size + 1;
               end loop;
            end;
         end if;
      end return;
   end Diagonal;

   overriding
   function Diagonal (Elements : CPU_Tensor; Offset : Integer := 0) return CPU_Tensor is
      Object_Data : Element_Array (1 .. Elements.Elements)
        with Import, Convention => Ada, Address => Elements.Data'Address;
   begin
      return Diagonal (Object_Data, Offset);
   end Diagonal;

   overriding
   function Trace (Object : CPU_Tensor; Offset : Integer := 0) return Element
     renames Operations.Trace;

   overriding
   function Reshape (Object : CPU_Tensor; Shape : Tensor_Shape) return CPU_Tensor is
     (Axes => Shape'Length,
      Size       => Object.Size,
      Kind       => Object.Kind,
      Shape      => Shape,
      Data       => Object.Data);

   overriding
   function Reshape (Object : CPU_Tensor; Elements : Positive) return CPU_Tensor
     renames Operations.Reshape;

   overriding
   function Flatten (Object : CPU_Tensor) return CPU_Tensor renames Operations.Flatten;

   overriding
   function Concatenate
     (Left, Right : CPU_Tensor;
      Axis   : Tensor_Axis) return CPU_Tensor
   is
      Shape : constant Tensor_Shape := Add (Left.Shape, Right.Shape, Axis);
      pragma Assert (Elements (Shape) = Left.Elements + Right.Elements);
   begin
      return Result : CPU_Tensor := Without_Data (Shape, Left.Kind) do
         declare
            Result_Data : Element_Array (1 .. Result.Elements)
              with Import, Convention => Ada, Address => Result.Data'Address;

            Left_Data : Element_Array (1 .. Left.Elements)
              with Import, Convention => Ada, Address => Left.Data'Address;

            Right_Data : Element_Array (1 .. Right.Elements)
              with Import, Convention => Ada, Address => Right.Data'Address;
         begin
            case Axis is
               when 1 =>
                  Result_Data (1 .. Left.Elements) := Left_Data;
                  Result_Data (Left.Elements + 1 .. Result_Data'Last) := Right_Data;
               when 2 =>
                  declare
                     Rows : constant Positive := Left.Rows;

                     Columns_Left  : constant Positive := Left.Columns;
                     Columns_Right : constant Positive := Right.Columns;
                  begin
                     for Index in 1 .. Rows loop
                        declare
                           Base_Left_Index : constant Natural :=
                             (Index - 1) * (Columns_Left + Columns_Right);
                           Base_Right_Index : constant Natural := Base_Left_Index + Columns_Left;
                           Left_Index  : constant Natural := (Index - 1) * Columns_Left;
                           Right_Index : constant Natural := (Index - 1) * Columns_Right;
                        begin
                           Result_Data (Base_Left_Index + 1 .. Base_Left_Index + Columns_Left)
                             := Left_Data (Left_Index + 1 .. Left_Index + Columns_Left);
                           Result_Data (Base_Right_Index + 1 .. Base_Right_Index + Columns_Right)
                             := Right_Data (Right_Index + 1 .. Right_Index + Columns_Right);
                        end;
                     end loop;
                  end;
               when others =>
                  raise Not_Implemented_Yet;  --  FIXME
            end case;
         end;
      end return;
   end Concatenate;

   overriding
   function "&" (Left, Right : CPU_Tensor) return CPU_Tensor renames Operations."&";

   ----------------------------------------------------------------------------
   --                            Matrix operations                           --
   ----------------------------------------------------------------------------

   procedure Multiply_Add
     (Result : in out Element_Array;
      Left   : Element;
      Right  : Element_Array)
   is
      Left_Vector : constant Vector_Type := (others => Left);

      subtype Vector_Elements is Element_Array (1 .. Vector_Type'Length);

      function Convert is new Ada.Unchecked_Conversion (Vector_Elements, Vector_Type);
      function Convert is new Ada.Unchecked_Conversion (Vector_Type, Vector_Elements);

      Last_Vector_Index : constant Positive := Data_Vectors (Right'Length);

      Padding : constant Natural :=
        Data_Padding (Last_Vector_Index, Right'Length);
   begin
      for Index in 1 .. Last_Vector_Index - (if Padding > 0 then 1 else 0) loop
         declare
            Vector_Index : constant Integer := (Index - 1) * Vector_Type'Length - 1;

            Result_Index : constant Natural := Vector_Index + Result'First;
            Right_Index  : constant Natural := Vector_Index + Right'First;

            Result_Vector : constant Vector_Type :=
              Convert (Result (Result_Index + 1 .. Result_Index + Vector_Type'Length));
            Right_Vector : constant Vector_Type :=
              Convert (Right (Right_Index + 1 .. Right_Index + Vector_Type'Length));
         begin
            Result (Result_Index + 1 .. Result_Index + Vector_Type'Length) :=
              Convert (Result_Vector + Left_Vector * Right_Vector);
         end;
      end loop;

      if Padding > 0 then
         declare
            Vector_Index : constant Integer := (Last_Vector_Index - 1) * Vector_Type'Length - 1;

            Result_Index : constant Natural := Vector_Index + Result'First;
            Right_Index  : constant Natural := Vector_Index + Right'First;

            Last_Count : constant Positive := Vector_Type'Length - Padding;
         begin
            for I in 1 .. Last_Count loop
               Result (Result_Index + I) :=
                 Result (Result_Index + I) + Left * Right (Right_Index + I);
            end loop;
         end;
      end if;
   end Multiply_Add;

   overriding
   function "*" (Left, Right : CPU_Tensor) return CPU_Tensor is
      --  m x n * n x p
      --      ^   ^
      --      |___|
      --     (Count)
      Left_Rows     : constant Natural := (if Left.Axes = 2 then Left.Rows else 1);
      Count         : constant Natural := Right.Rows;
      Right_Columns : constant Natural := (if Right.Axes = 2 then Right.Columns else 1);

      Shape : constant Tensor_Shape :=
         (case Right.Axes is
            when 1 => (1 => Left_Rows),
            when 2 => (1 => Left_Rows, 2 => Right_Columns),
            when others => raise Not_Implemented_Yet);  --  FIXME
   begin
      --  Matrix-matrix or matrix-vector or vector-matrix multiplication
      return Result : CPU_Tensor := Zeros (Shape) do
         declare
            Result_Data : Element_Array (1 .. Result.Elements)
              with Import, Convention => Ada, Address => Result.Data'Address;
            Left_Data : Element_Array (1 .. Left.Elements)
              with Import, Convention => Ada, Address => Left.Data'Address;
            Right_Data : Element_Array (1 .. Right.Elements)
              with Import, Convention => Ada, Address => Right.Data'Address;
         begin
            for Row_Index in 1 .. Left_Rows loop
               --  Result (Row_Index) := Left (Row_Index) * Right;
               declare
                  Result_Index : constant Natural := To_Index ((Row_Index, 1), Right_Columns) - 1;
               begin
                  for Column_Index in 1 .. Count loop
                     declare
                        Right_Index : constant Natural :=
                          To_Index ((Column_Index, 1), Right_Columns) - 1;
                     begin
                        --  Left_Value := Left (Row_Index, Column_Index)
                        --  Result (Row_Index) := @ + Left_Value * Right (Row_Index)
                        Multiply_Add
                          (Result_Data (Result_Index + 1 .. Result_Index + Right_Columns),
                           Left_Data (To_Index ((Row_Index, Column_Index), Count)),
                           Right_Data (Right_Index + 1 .. Right_Index + Right_Columns));
                     end;
                  end loop;
               end;
            end loop;
         end;
      end return;
   end "*";

   overriding
   function "*" (Left, Right : CPU_Tensor) return Element is
      Result : Element_Type := 0.0;

      Padding : constant Natural :=
        Data_Padding (Size => Left.Size, Count => Left.Elements);

      Last_Left  : constant Vector_Type := Reset_Padding (Left, Padding, 0.0);
      Last_Right : constant Vector_Type := Reset_Padding (Right, Padding, 0.0);

      type Sum_Index_Type is mod 8;

      Sums : array (Sum_Index_Type) of Vector_Type := (others => (others => 0.0));
      --  TODO Do pairwise summation recursively
   begin
      for Index in Left.Data'First .. Left.Data'Last - 1 loop
         declare
            Sum_Index : constant Sum_Index_Type := Sum_Index_Type (Index mod Sums'Length);
         begin
            Sums (Sum_Index) := Sums (Sum_Index) + Left.Data (Index) * Right.Data (Index);
         end;
      end loop;

      declare
         Sum_Index : constant Sum_Index_Type := Sum_Index_Type (Left.Data'Last mod Sums'Length);
      begin
         Sums (Sum_Index) := Sums (Sum_Index) + Last_Left * Last_Right;
      end;

      for Value of Sums loop
         Result := Result + Sum (Value);
      end loop;

      return Element (Result);
   end "*";

   overriding function "**" (Left : CPU_Tensor; Right : Integer) return CPU_Tensor
     renames Operations."**";

   overriding
   function Outer (Left, Right : CPU_Tensor) return CPU_Tensor is
      Shape : constant Tensor_Shape := (1 => Left.Elements, 2 => Right.Elements);
   begin
      return Result : CPU_Tensor := Without_Data (Shape) do
         declare
            Result_Data : Element_Array (1 .. Result.Elements)
              with Import, Convention => Ada, Address => Result.Data'Address;

            Left_Data : Element_Array (1 .. Left.Elements)
              with Import, Convention => Ada, Address => Left.Data'Address;

            Columns : constant Positive := Right.Elements;
         begin
            for Index in 1 .. Left.Elements loop
               declare
                  Row : constant CPU_Tensor := Left_Data (Index) * Right;

                  Row_Data : Element_Array (1 .. Columns)
                    with Import, Convention => Ada, Address => Row.Data'Address;

                  Base_Index : constant Natural := (Index - 1) * Columns;
               begin
                  Result_Data (Base_Index + 1 .. Base_Index + Columns) := Row_Data;
               end;
            end loop;
         end;
      end return;
   end Outer;

   overriding
   function Inverse (Object : CPU_Tensor) return CPU_Tensor renames Operations.Inverse;

   overriding
   function Transpose (Object : CPU_Tensor) return CPU_Tensor is
      Shape : constant Tensor_Shape :=
        (1 => Object.Columns,
         2 => Object.Rows);
   begin
      return Result : CPU_Tensor := Without_Data (Shape) do
         declare
            Result_Data : Element_Array (1 .. Result.Elements)
              with Import, Convention => Ada, Address => Result.Data'Address;

            Object_Data : Element_Array (1 .. Object.Elements)
              with Import, Convention => Ada, Address => Object.Data'Address;

            Rows    : constant Natural := Object.Rows;
            Columns : constant Natural := Object.Columns;

            Result_Columns : Natural renames Rows;
         begin
            for Row_Index in 1 .. Rows loop
               for Column_Index in 1 .. Columns loop
                  Result_Data (To_Index ((Column_Index, Row_Index), Result_Columns)) :=
                    Object_Data (To_Index ((Row_Index, Column_Index), Columns));
               end loop;
            end loop;
         end;
      end return;
   end Transpose;

   ----------------------------------------------------------------------------

   overriding
   function Solve (A, B : CPU_Tensor; Solution : out Solution_Kind) return CPU_Tensor
     renames Operations.Solve;

   overriding
   function Solve (A, B : CPU_Tensor; Form : Triangular_Form) return CPU_Tensor
     renames Operations.Solve;

   overriding
   function Divide_By (B, A : CPU_Tensor) return CPU_Tensor
     renames Operations.Divide_By;

   overriding
   function Divide_By (B, A : CPU_Tensor; Form : Triangular_Form) return CPU_Tensor
     renames Operations.Divide_By;

   overriding
   function QR (Object : CPU_Tensor) return CPU_Tensor
     renames Operations.QR;

   overriding
   function QR (Object : CPU_Tensor; Mode : QR_Mode := Reduced) return QR_Factorization'Class
     renames Operations.QR;

   overriding
   function QR_For_Least_Squares (Object : CPU_Tensor) return QR_Factorization'Class
     renames Operations.QR_For_Least_Squares;

   overriding
   function Least_Squares (Object : QR_Factorization'Class; B : CPU_Tensor) return CPU_Tensor
     renames Operations.Least_Squares;

   overriding
   function Least_Squares (A, B : CPU_Tensor) return CPU_Tensor
     renames Operations.Least_Squares;

   overriding
   function Constrained_Least_Squares (A, B, C, D : CPU_Tensor) return CPU_Tensor
     renames Operations.Constrained_Least_Squares;

   overriding
   function Cholesky (Object : CPU_Tensor; Form : Triangular_Form := Lower) return CPU_Tensor
     renames Operations.Cholesky;

   overriding
   function Cholesky_Update
     (R, V : CPU_Tensor;
      Mode : Update_Mode) return CPU_Tensor renames Operations.Cholesky_Update;

   ----------------------------------------------------------------------------
   --                            Vector operations                           --
   ----------------------------------------------------------------------------

   overriding
   function Norm (Object : CPU_Tensor) return Element renames Operations.Norm;

   overriding
   function Normalize (Object : CPU_Tensor) return CPU_Tensor renames Operations.Normalize;

   overriding
   function Standardize (Object : CPU_Tensor) return CPU_Tensor renames Operations.Standardize;

   overriding
   function Correlation_Coefficient (Left, Right : CPU_Tensor) return Correlation_Element
     renames Operations.Correlation_Coefficient;

   ----------------------------------------------------------------------------
   --                         Element-wise operations                        --
   ----------------------------------------------------------------------------

   overriding function "+" (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) + Right.Data (Index);
         end loop;
      end return;
   end "+";

   overriding function "-" (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) - Right.Data (Index);
         end loop;
      end return;
   end "-";

   overriding function "/" (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) / Right.Data (Index);
         end loop;
      end return;
   end "/";

   overriding function Divide_Or_Zero (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Divide_Or_Zero (Left.Data (Index), Right.Data (Index));
         end loop;
      end return;
   end Divide_Or_Zero;

   function Power (Left, Right : Element) return Element is
     (if Left = 0.0 and Right = 0.0 then 1.0 else EF."**" (Left, Right));
   --  0.0 ** X = 1.0 if X = 0.0 (EF."**" raises Argument_Error instead of returning 1.0)
   --  0.0 ** X = 0.0 if X > 0.0

   overriding function "**" (Left, Right : CPU_Tensor) return CPU_Tensor is
--     (Exp (Multiply (Right, Log (Left))));

      Padding : constant Natural :=
        Data_Padding (Size => Left.Size, Count => Elements (Left.Shape));

      Last_Left  : constant Vector_Type := Reset_Padding (Left, Padding, 0.0);
      Last_Right : constant Vector_Type := Reset_Padding (Right, Padding, 1.0);
   begin
      return Result : CPU_Tensor := Without_Data (Left) do
         for Index in Result.Data'First .. Result.Data'Last - 1 loop
            for Offset in Vector_Type'Range loop
               declare
                  Left_Element  : Element renames Left.Data (Index) (Offset);
                  Right_Element : Element renames Right.Data (Index) (Offset);
               begin
                  Result.Data (Index) (Offset) := Power (Left_Element, Right_Element);
               end;
            end loop;
         end loop;

         for Offset in Vector_Type'Range loop
            Result.Data (Result.Data'Last) (Offset) :=
              Power (Last_Left (Offset), Last_Right (Offset));
         end loop;
      end return;
   end "**";

   overriding function "**" (Left : CPU_Tensor; Right : Element) return CPU_Tensor is
--     (Exp (Right * Log (Left)));
      use EF;
   begin
      if Right = 0.0 then
         return Ones (Left.Shape);
      elsif Right = 1.0 then
         return Left;
      else
         return Result : CPU_Tensor := Without_Data (Left) do
            for Index in Result.Data'Range loop
               for Offset in Vector_Type'Range loop
                  Result.Data (Index) (Offset) := Left.Data (Index) (Offset) ** Right;
               end loop;
            end loop;
         end return;
      end if;
   end "**";

   overriding function "**" (Left : Element; Right : CPU_Tensor) return CPU_Tensor is
--     (Exp (Right * EF.Log (Left)));

      Padding : constant Natural :=
        Data_Padding (Size => Right.Size, Count => Elements (Right.Shape));

      --  EF."**" raises Constraint_Error if Left = 0.0 and element
      --  in the padding happens to be < 0.0
      Last_Right : constant Vector_Type := Reset_Padding (Right, Padding, 1.0);
   begin
      if Left = 1.0 then
         return Ones (Right.Shape);
      else
         return Result : CPU_Tensor := Without_Data (Right) do
            for Index in Result.Data'First .. Result.Data'Last - 1 loop
               for Offset in Vector_Type'Range loop
                  declare
                     Right_Element : Element renames Right.Data (Index) (Offset);
                  begin
                     Result.Data (Index) (Offset) := Power (Left, Right_Element);
                  end;
               end loop;
            end loop;

            for Offset in Vector_Type'Range loop
               Result.Data (Result.Data'Last) (Offset) := Power (Left, Last_Right (Offset));
            end loop;
         end return;
      end if;
   end "**";

   overriding function "*" (Left : CPU_Tensor; Right : Element) return CPU_Tensor is
      Right_Vector : constant Vector_Type := (others => Right);
   begin
      return Result : CPU_Tensor := Without_Data (Left) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) * Right_Vector;
         end loop;
      end return;
   end "*";

   overriding function "/" (Left : Element; Right : CPU_Tensor) return CPU_Tensor is
      Left_Vector : constant Vector_Type := (others => Left);
   begin
      return Result : CPU_Tensor := Without_Data (Right) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left_Vector / Right.Data (Index);
         end loop;
      end return;
   end "/";

   overriding function "/" (Left : CPU_Tensor; Right : Element) return CPU_Tensor is
      Right_Vector : constant Vector_Type := (others => Right);
   begin
      return Result : CPU_Tensor := Without_Data (Left) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) / Right_Vector;
         end loop;
      end return;
   end "/";

   overriding function "+" (Left : CPU_Tensor; Right : Element) return CPU_Tensor is
      Right_Vector : constant Vector_Type := (others => Right);
   begin
      return Result : CPU_Tensor := Without_Data (Left) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) + Right_Vector;
         end loop;
      end return;
   end "+";

   overriding function "-" (Left : Element; Right : CPU_Tensor) return CPU_Tensor is
      Left_Vector : constant Vector_Type := (others => Left);
   begin
      return Result : CPU_Tensor := Without_Data (Right) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left_Vector - Right.Data (Index);
         end loop;
      end return;
   end "-";

   overriding function "-" (Object : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Object) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := -Object.Data (Index);
         end loop;
      end return;
   end "-";

   overriding
   function "*" (Left : Element; Right : CPU_Tensor) return CPU_Tensor renames Operations."*";

   overriding
   function "+" (Left : Element; Right : CPU_Tensor) return CPU_Tensor renames Operations."+";

   overriding
   function "-" (Left : CPU_Tensor; Right : Element) return CPU_Tensor renames Operations."-";

   overriding
   function "mod" (Left, Right : CPU_Tensor) return CPU_Tensor renames Operations."mod";

   overriding
   function "rem" (Left, Right : CPU_Tensor) return CPU_Tensor renames Operations."rem";

   overriding
   function "mod" (Left : CPU_Tensor; Right : Element) return CPU_Tensor renames Operations."mod";

   overriding
   function "rem" (Left : CPU_Tensor; Right : Element) return CPU_Tensor renames Operations."rem";

   overriding function "abs" (Object : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Object) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := abs Object.Data (Index);
         end loop;
      end return;
   end "abs";

   overriding function Multiply (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) * Right.Data (Index);
         end loop;
      end return;
   end Multiply;

   overriding
   function Power (Left : CPU_Tensor; Right : Integer) return CPU_Tensor renames Operations.Power;

   overriding
   function Min (Left : Element; Right : CPU_Tensor) return CPU_Tensor renames Operations.Min;

   overriding
   function Max (Left : Element; Right : CPU_Tensor) return CPU_Tensor renames Operations.Max;

   overriding function Min (Left : CPU_Tensor; Right : Element) return CPU_Tensor is
      Right_Vector : constant Vector_Type := (others => Right);
   begin
      return Result : CPU_Tensor := Without_Data (Left) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Min (Left.Data (Index), Right_Vector);
         end loop;
      end return;
   end Min;

   overriding function Max (Left : CPU_Tensor; Right : Element) return CPU_Tensor is
      Right_Vector : constant Vector_Type := (others => Right);
   begin
      return Result : CPU_Tensor := Without_Data (Left) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Max (Left.Data (Index), Right_Vector);
         end loop;
      end return;
   end Max;

   overriding function Sqrt (Object : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Object) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Sqrt (Object.Data (Index));
         end loop;
      end return;
   end Sqrt;

   overriding function Ceil (Object : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Object) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Ceil (Object.Data (Index));
         end loop;
      end return;
   end Ceil;

   overriding function Floor (Object : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Object) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Floor (Object.Data (Index));
         end loop;
      end return;
   end Floor;

   overriding function Round (Object : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Object) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Round (Object.Data (Index));
         end loop;
      end return;
   end Round;

   overriding function Truncate (Object : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Object) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Truncate (Object.Data (Index));
         end loop;
      end return;
   end Truncate;

   overriding function Exp (Object : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Object) do
         for Index in Result.Data'Range loop
            for Offset in Vector_Type'Range loop
               Result.Data (Index) (Offset) := EF.Exp (Object.Data (Index) (Offset));
            end loop;
         end loop;
      end return;
   end Exp;

   overriding function Log (Object : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Object) do
         for Index in Result.Data'Range loop
            for Offset in Vector_Type'Range loop
               Result.Data (Index) (Offset) := EF.Log (Object.Data (Index) (Offset));
            end loop;
         end loop;
      end return;
   end Log;

   overriding function Log10 (Object : CPU_Tensor) return CPU_Tensor renames Operations.Log10;

   overriding function Log2 (Object : CPU_Tensor) return CPU_Tensor renames Operations.Log2;

   ----------------------------------------------------------------------------
   --                              Trigonometry                              --
   ----------------------------------------------------------------------------

   overriding function Sin (Object : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Object) do
         for Index in Result.Data'Range loop
            for Offset in Vector_Type'Range loop
               Result.Data (Index) (Offset) := EF.Sin (Object.Data (Index) (Offset));
            end loop;
         end loop;
      end return;
   end Sin;

   overriding function Cos (Object : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Object) do
         for Index in Result.Data'Range loop
            for Offset in Vector_Type'Range loop
               Result.Data (Index) (Offset) := EF.Cos (Object.Data (Index) (Offset));
            end loop;
         end loop;
      end return;
   end Cos;

   overriding function Tan (Object : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Object) do
         for Index in Result.Data'Range loop
            for Offset in Vector_Type'Range loop
               Result.Data (Index) (Offset) := EF.Tan (Object.Data (Index) (Offset));
            end loop;
         end loop;
      end return;
   end Tan;

   overriding function Arcsin (Object : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Object) do
         for Index in Result.Data'Range loop
            for Offset in Vector_Type'Range loop
               Result.Data (Index) (Offset) := EF.Arcsin (Object.Data (Index) (Offset));
            end loop;
         end loop;
      end return;
   end Arcsin;

   overriding function Arccos (Object : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Object) do
         for Index in Result.Data'Range loop
            for Offset in Vector_Type'Range loop
               Result.Data (Index) (Offset) := EF.Arccos (Object.Data (Index) (Offset));
            end loop;
         end loop;
      end return;
   end Arccos;

   overriding function Arctan (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left) do
         for Index in Result.Data'Range loop
            for Offset in Vector_Type'Range loop
               Result.Data (Index) (Offset) :=
                 EF.Arctan (Left.Data (Index) (Offset), Right.Data (Index) (Offset));
            end loop;
         end loop;
      end return;
   end Arctan;

   overriding function Degrees (Object : CPU_Tensor) return CPU_Tensor renames Operations.Degrees;

   overriding function Radians (Object : CPU_Tensor) return CPU_Tensor renames Operations.Radians;

   ----------------------------------------------------------------------------
   --                               Reductions                               --
   ----------------------------------------------------------------------------

   overriding
   function Reduce_Associative
     (Object    : CPU_Tensor;
      Subject   : Expression'Class;
      Initial   : Element) return Element
   is
      CPU_Subject : constant Expression_Type := Expression_Type (Subject);

      Max_Length_Sequential : constant Positive := 128 / Vector_Type'Length;

      function Pairwise (Lower, Upper : Positive) return Vector_Type is
         Length : constant Positive := Upper - Lower + 1;
      begin
         if Length <= Max_Length_Sequential then
            declare
               Result : Vector_Type := Object.Data (Lower);
            begin
               for Index in Lower + 1 .. Upper loop
                  Result := Apply (CPU_Subject, Result, Object.Data (Index));
               end loop;
               return Result;
            end;
         else
            declare
               Half_Index : constant Positive := Lower + Length / 2;
            begin
               return Apply
                 (CPU_Subject,
                  Pairwise (Lower, Half_Index - 1),
                  Pairwise (Half_Index, Upper));
            end;
         end if;
      end Pairwise;

      Padding : constant Natural :=
        Data_Padding (Size => Object.Size, Count => Object.Elements);

      Result : Element_Type := Initial;
   begin
      if Object.Data'Last = 0 then
         return Result;
      end if;

      if Object.Elements > Vector_Type'Length then
         for Element of Pairwise (Object.Data'First, Object.Data'Last - 1) loop
            Result := Apply (CPU_Subject, Result, Element);
         end loop;
      else
         --  Do not perform pairwise applying of expression because Object.Data has only 1 vector,
         --  which is added to Result below
         null;
      end if;

      for Index in Vector_Index_Type'First .. From_Last (Padding) loop
         Result := Apply (CPU_Subject, Result, Object.Data (Object.Data'Last) (Index));
      end loop;

      return Result;
   end Reduce_Associative;

   overriding
   function Reduce_Associative
     (Object  : CPU_Tensor;
      Subject : Expression'Class;
      Initial : Element;
      Axis    : Tensor_Axis) return CPU_Tensor is
   begin
      raise Program_Error;
      return Zeros ((1 => 1));  --  FIXME
   end Reduce_Associative;

   overriding
   function Reduce
     (Object    : CPU_Tensor;
      Subject   : Expression'Class;
      Initial   : Element) return Element
   is
      CPU_Subject : constant Expression_Type := Expression_Type (Subject);

      Data : Element_Array (1 .. Object.Elements)
        with Import, Convention => Ada, Address => Object.Data'Address;

      Result : Element_Type := Initial;
   begin
      for Value of Data loop
         Result := Apply (CPU_Subject, Result, Value);
      end loop;

      return Result;
   end Reduce;

   overriding
   function Reduce
     (Object  : CPU_Tensor;
      Subject : Expression'Class;
      Initial : Element;
      Axis    : Tensor_Axis) return CPU_Tensor is
   begin
      raise Program_Error;
      return Zeros ((1 => 1));  --  FIXME
   end Reduce;

   overriding function Sum (Object : CPU_Tensor) return Element renames Operations.Sum;

   overriding function Product (Object : CPU_Tensor) return Element renames Operations.Product;

   overriding
   function Sum (Object : CPU_Tensor; Axis : Tensor_Axis) return CPU_Tensor
     renames Operations.Sum;

   overriding
   function Product (Object : CPU_Tensor; Axis : Tensor_Axis) return CPU_Tensor
     renames Operations.Product;

   ----------------------------------------------------------------------------
   --                               Statistics                               --
   ----------------------------------------------------------------------------

   overriding function Min (Object : CPU_Tensor) return Element renames Operations.Min;

   overriding function Max (Object : CPU_Tensor) return Element renames Operations.Max;

   overriding function Quantile (Object : CPU_Tensor; P : Probability) return Element
     renames Operations.Quantile;

   overriding function Median (Object : CPU_Tensor) return Element
     renames Operations.Median;

   overriding function Mean (Object : CPU_Tensor) return Element
     renames Operations.Mean;

   overriding
   function Variance (Object : CPU_Tensor; Offset : Natural := 0) return Element
     renames Operations.Variance;

   overriding
   function Standard_Deviation (Object : CPU_Tensor; Offset : Natural := 0) return Element
     renames Operations.Standard_Deviation;

   overriding
   function Min (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Min (Left.Data (Index), Right.Data (Index));
         end loop;
      end return;
   end Min;

   overriding
   function Max (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Max (Left.Data (Index), Right.Data (Index));
         end loop;
      end return;
   end Max;

   overriding
   function Min (Object : CPU_Tensor; Axis : Tensor_Axis) return CPU_Tensor
     renames Operations.Min;

   overriding
   function Max (Object : CPU_Tensor; Axis : Tensor_Axis) return CPU_Tensor
     renames Operations.Max;

   overriding
   function Quantile
     (Object : CPU_Tensor;
      P      : Probability;
      Axis   : Tensor_Axis) return CPU_Tensor is
   begin
      raise Program_Error;
      return Zeros ((1 => 1));  --  FIXME
   end Quantile;

   overriding
   function Mean (Object : CPU_Tensor; Axis : Tensor_Axis) return CPU_Tensor
     renames Operations.Mean;

   overriding
   function Variance
     (Object : CPU_Tensor;
      Axis   : Tensor_Axis;
      Offset : Natural := 0) return CPU_Tensor
   renames Operations.Variance;

   overriding
   function Median (Object : CPU_Tensor; Axis : Tensor_Axis) return CPU_Tensor
     renames Operations.Median;

   overriding
   function Standard_Deviation
     (Object : CPU_Tensor;
      Axis   : Tensor_Axis;
      Offset : Natural := 0) return CPU_Tensor
   renames Operations.Standard_Deviation;

   ----------------------------------------------------------------------------
   --                                Logical                                 --
   ----------------------------------------------------------------------------

   overriding function And_Not (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left, Kind => Bool_Type) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := And_Not (Left.Data (Index), Right.Data (Index));
         end loop;
         Result.Data (Result.Data'Last) := Disable_Padding (Result);
      end return;
   end And_Not;

   overriding function "and" (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left, Kind => Left.Kind) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) and Right.Data (Index);
         end loop;
         Result.Data (Result.Data'Last) := Disable_Padding (Result);
      end return;
   end "and";

   overriding function "and" (Left : Element; Right : CPU_Tensor) return CPU_Tensor is
      Left_Vector : constant Vector_Type := (others => Left);
   begin
      return Result : CPU_Tensor := Without_Data (Right, Kind => Float_Type) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left_Vector and Right.Data (Index);
         end loop;
         Result.Data (Result.Data'Last) := Disable_Padding (Result);
      end return;
   end "and";

   overriding function "or"  (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left, Kind => Bool_Type) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) or Right.Data (Index);
         end loop;
         Result.Data (Result.Data'Last) := Disable_Padding (Result);
      end return;
   end "or";

   overriding function "xor" (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left, Kind => Bool_Type) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) xor Right.Data (Index);
         end loop;
         Result.Data (Result.Data'Last) := Disable_Padding (Result);
      end return;
   end "xor";

   overriding function "not"  (Object : CPU_Tensor) return CPU_Tensor is
      Zero_Vector : constant Vector_Type := (others => 0.0);

      Mask : constant Vector_Type := Zero_Vector = Zero_Vector;
   begin
      return Result : CPU_Tensor := Without_Data (Object, Kind => Bool_Type) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := And_Not (Object.Data (Index), Mask);
         end loop;
         Result.Data (Result.Data'Last) := Disable_Padding (Result);
      end return;
   end "not";

   ----------------------------------------------------------------------------
   --                              Comparisons                               --
   ----------------------------------------------------------------------------

   overriding
   function "="  (Left : CPU_Tensor; Right : Element) return CPU_Tensor renames Operations."=";

   overriding
   function "/=" (Left : CPU_Tensor; Right : Element) return CPU_Tensor renames Operations."/=";

   overriding function ">"  (Left : CPU_Tensor; Right : Element) return CPU_Tensor is
      Right_Vector : constant Vector_Type := (others => Right);
   begin
      return Result : CPU_Tensor := Without_Data (Left, Kind => Bool_Type) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) > Right_Vector;
         end loop;
         Result.Data (Result.Data'Last) := Disable_Padding (Result);
      end return;
   end ">";

   overriding function "<"  (Left : CPU_Tensor; Right : Element) return CPU_Tensor is
      Right_Vector : constant Vector_Type := (others => Right);
   begin
      return Result : CPU_Tensor := Without_Data (Left, Kind => Bool_Type) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) < Right_Vector;
         end loop;
         Result.Data (Result.Data'Last) := Disable_Padding (Result);
      end return;
   end "<";

   overriding function ">=" (Left : CPU_Tensor; Right : Element) return CPU_Tensor is
      Right_Vector : constant Vector_Type := (others => Right);
   begin
      return Result : CPU_Tensor := Without_Data (Left, Kind => Bool_Type) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) >= Right_Vector;
         end loop;
         Result.Data (Result.Data'Last) := Disable_Padding (Result);
      end return;
   end ">=";

   overriding function "<=" (Left : CPU_Tensor; Right : Element) return CPU_Tensor is
      Right_Vector : constant Vector_Type := (others => Right);
   begin
      return Result : CPU_Tensor := Without_Data (Left, Kind => Bool_Type) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) <= Right_Vector;
         end loop;
         Result.Data (Result.Data'Last) := Disable_Padding (Result);
      end return;
   end "<=";

   ----------------------------------------------------------------------------

   overriding
   function "="  (Left : Element; Right : CPU_Tensor) return CPU_Tensor renames Operations."=";

   overriding
   function "/=" (Left : Element; Right : CPU_Tensor) return CPU_Tensor renames Operations."/=";

   overriding
   function ">"  (Left : Element; Right : CPU_Tensor) return CPU_Tensor renames Operations.">";

   overriding
   function "<"  (Left : Element; Right : CPU_Tensor) return CPU_Tensor renames Operations."<";

   overriding
   function ">=" (Left : Element; Right : CPU_Tensor) return CPU_Tensor renames Operations.">=";

   overriding
   function "<=" (Left : Element; Right : CPU_Tensor) return CPU_Tensor renames Operations."<=";

   ----------------------------------------------------------------------------

   overriding function "=" (Left, Right : CPU_Tensor) return Boolean renames Operations."=";

   overriding function "="  (Left, Right : CPU_Tensor) return CPU_Tensor is
      Zero_Vector : constant Vector_Type := (others => 0.0);

      Mask : constant Vector_Type := Zero_Vector = Zero_Vector;
   begin
      case Left.Kind is
         when Float_Type =>
            return (abs (Left - Right) <= Element_Type'Model_Epsilon);
         when Int_Type | Bool_Type =>
            return Result : CPU_Tensor := Without_Data (Left, Kind => Bool_Type) do
               for Index in Result.Data'Range loop
                  Result.Data (Index) :=
                    And_Not (Left.Data (Index) xor Right.Data (Index), Mask);
               end loop;
               Result.Data (Result.Data'Last) := Disable_Padding (Result);
            end return;
      end case;
   end "=";

   overriding function "/=" (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      case Left.Kind is
         when Float_Type =>
            return (abs (Left - Right) > Element_Type'Model_Epsilon);
         when Int_Type | Bool_Type =>
            return Result : CPU_Tensor := Without_Data (Left, Kind => Bool_Type) do
               for Index in Result.Data'Range loop
                  Result.Data (Index) := Left.Data (Index) xor Right.Data (Index);
               end loop;
               Result.Data (Result.Data'Last) := Disable_Padding (Result);
            end return;
      end case;
   end "/=";

   overriding function ">"  (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left, Kind => Bool_Type) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) > Right.Data (Index);
         end loop;
         Result.Data (Result.Data'Last) := Disable_Padding (Result);
      end return;
   end ">";

   overriding function "<"  (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left, Kind => Bool_Type) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) < Right.Data (Index);
         end loop;
         Result.Data (Result.Data'Last) := Disable_Padding (Result);
      end return;
   end "<";

   overriding function ">=" (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left, Kind => Bool_Type) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) >= Right.Data (Index);
         end loop;
         Result.Data (Result.Data'Last) := Disable_Padding (Result);
      end return;
   end ">=";

   overriding function "<=" (Left, Right : CPU_Tensor) return CPU_Tensor is
   begin
      return Result : CPU_Tensor := Without_Data (Left, Kind => Bool_Type) do
         for Index in Result.Data'Range loop
            Result.Data (Index) := Left.Data (Index) <= Right.Data (Index);
         end loop;
         Result.Data (Result.Data'Last) := Disable_Padding (Result);
      end return;
   end "<=";

   ----------------------------------------------------------------------------

   overriding
   function All_Close
     (Left, Right        : CPU_Tensor;
      Relative_Tolerance : Element := 1.0e-05;
      Absolute_Tolerance : Element := Element_Type'Model_Epsilon) return Boolean
   renames Operations.All_Close;

   overriding
   function Any_True (Object : CPU_Tensor; Axis : Tensor_Axis) return CPU_Tensor is
   begin
      raise Program_Error;
      return Zeros ((1 => 1));  --  FIXME
   end Any_True;

   overriding
   function Any_True (Object : CPU_Tensor) return Boolean is
      Padding : constant Natural :=
        Data_Padding (Size => Object.Size, Count => Object.Elements);

      function All_Zeros (Vector : Vector_Type) return Boolean is
        (All_Zeros (Convert (Vector), Convert (Vector)));

      Zero_Vector : constant Vector_Type := (others => 0.0);
      Mask : Integer_Vector_Type := Convert (Zero_Vector = Zero_Vector);
   begin
      for Index in 1 .. Padding loop
         Mask := Shift_Elements_Right (Mask);
      end loop;

      return (for some Index in Object.Data'First .. Object.Data'Last - 1 =>
                not All_Zeros (Object.Data (Index)))
             or not All_Zeros (Object.Data (Object.Data'Last) and Convert (Mask));
   end Any_True;

   overriding
   function All_True (Object : CPU_Tensor; Axis : Tensor_Axis) return CPU_Tensor is
   begin
      raise Program_Error;
      return Zeros ((1 => 1));  --  FIXME
   end All_True;

   overriding
   function All_True (Object : CPU_Tensor) return Boolean is
      Padding : constant Natural :=
        Data_Padding (Size => Object.Size, Count => Object.Elements);

      function All_Ones (Vector : Vector_Type) return Boolean is
        (All_Ones (Convert (Vector), Convert (Vector = Vector)));

      Zero_Vector : constant Vector_Type := (others => 0.0);
      Mask : Integer_Vector_Type := Convert (Zero_Vector = Zero_Vector);
   begin
      for Index in 1 .. Vector_Type'Length - Padding loop
         Mask := Shift_Elements_Left (Mask);
      end loop;

      return (for all Index in Object.Data'First .. Object.Data'Last - 1 =>
                All_Ones (Object.Data (Index)))
             and All_Ones (Object.Data (Object.Data'Last) or Convert (Mask));
   end All_True;

   procedure Reset_Random (Seed : Duration) is
   begin
      Reset (Random_State, Seed);
   end Reset_Random;

   overriding function Random_Uniform (Shape : Tensor_Shape) return CPU_Tensor is
      Value : Vector_Type;
   begin
      return Result : CPU_Tensor := Without_Data (Shape) do
         for I in Result.Data'Range loop
            Next (Random_State, Value);
            Result.Data (I) := Value;
         end loop;
      end return;
   end Random_Uniform;

end Orka.Numerics.Tensors.SIMD_CPU;