diff mbox series

[Ada] Primitive functions that require one formal and return an array

Message ID 20170906100211.GA9448@adacore.com
State New
Headers show
Series [Ada] Primitive functions that require one formal and return an array | expand

Commit Message

Arnaud Charlet Sept. 6, 2017, 10:02 a.m. UTC
Primitive functions whose first formal is a controlling parameter, whose
other formals have defaults and whose result is an array type can lead to
ambiguities when the result of such a call is the prefix of an indexed
component. The interpretation that analyzes Obj.F (X, Y) into F (Obj)(X, Y)
is only legal if the first parameter of F is a controlling parameter. This
additional guard was previously missing from the predicate, leading to
malformed trees and a compiler crash.

Compiling huckel.adb must yield:

huckel.adb:135:27: expected type "Real" defined at huckel.ads:9
huckel.adb:135:27: found type "Ada.Numerics.Generic_Real_Arrays.Real_Matrix"
                      from instance at huckel.ads:16

-- Huckel package
-- This is a translation from Fortran II code documented in the
-- book "Computing Methods for Quantum Organic Chemistry"
with Ada.Numerics.Generic_Real_Arrays;

package Huckel is
   type Real is digits 15;
   type Molecule (Atoms : Positive) is tagged private;
   function Input return Molecule;
   procedure Compute_Energies(Item : in out Molecule);
   procedure Output(Item : in Molecule);

private
   package Matrices is new Ada.Numerics.Generic_Real_Arrays(Real);
   use Matrices;
   type Molecule (Atoms : Positive) is tagged record
      Orbitals        : Positive;
      Atomic_Matrix   : Real_Matrix(1..Atoms, 1..Atoms);
      Atomic_Diagonal : Real_Vector(1..Atoms);
      Unit_Matrix     : Real_Matrix(1..Atoms, 1..Atoms);
      Bond_Orders     : Real_Matrix(1..Atoms, 1..Atoms);
      Free_Valences   : Real_vector(1..Atoms);
   end record;
end Huckel;
---
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Integer_Text_IO; use Ada.Integer_Text_IO;
with Ada.Text_IO;
with Ada.Numerics.Generic_Elementary_Functions;

package body Huckel is
   package Real_IO is new Ada.Text_IO.Float_IO(Real);
   use Real_Io;

   -----------
   -- Input --
   -----------

   function Input return Molecule is
      Num_Atoms : Positive;
      Num_Orbs  : Positive;
   begin
      Get(Item => Num_Atoms);
      Get(Item => Num_Orbs);
      declare
         Temp : Molecule(Atoms => Num_Atoms);
      begin
         Temp.Orbitals := Num_Orbs;
         -- Read the atomic matrix into the upper semi-matrix of Atomic_Matrix
         for I in 1..Num_Atoms loop
            for J in 1..I loop
               Get(Item => Temp.Atomic_Matrix(J, I));
               -- Print the input matrix in lower semi-matrix format
               Put(Item =>
                 Temp.Atomic_Matrix(J,I), Aft => 0, Fore => 2, Exp => 0);
               -- Make all bonding terms negative
               Temp.Atomic_Matrix(I, J) := -Temp.Atomic_Matrix(I,J);
            end loop;
            New_Line;
         end loop;
         return Temp;
      end;
   end Input;

   ------------
   -- Modify --
   ------------
   procedure Modify(Item : in out Molecule) is
      Num_Mods     : natural;
      I, J         : Positive;
      Modification : Real;
   begin
      Get(Item => Num_Mods);
      if Num_Mods > 0 then
         New_Line(3);
         Put_Line("Modifications");
         for Num in 1..Num_Mods loop
            Get(Item => I);
            Get(Item => J);
            Get(Item => Modification);
            Put(Item => I, Width => 3);
            Put(Item => J, Width => 6);
            Put(Item => Modification, Aft => 3, Fore => 7, Exp => 0);
            New_Line;
            if I = J then
               Item.Atomic_Diagonal(J) := Modification;
            elsif I < J then
               Item.Atomic_Matrix(I, J) := Modification;
            else
               Item.Atomic_Matrix(J, I) := Modification;
            end if;
         end loop;
      end if;
   end Modify;

   ----------
   -- Pahy --
   ----------
   procedure Pahy(Item : in out Molecule) is

   begin
      for J in 1..Item.Atoms loop
         for I in 1..J loop
            Item.Atomic_Matrix(I, J) := Item.Atomic_Matrix(J, I);
            Item.Atomic_Diagonal(J) := Item.Atomic_Matrix(J,J);
         end loop;
      end loop;
   end Pahy;

   ------------
   -- Scofi1 --
   ------------

   procedure Scofi1(Item : in out Molecule) is
    package elem_funcs is new Ada.Numerics.Generic_Elementary_Functions(real);
      use elem_funcs;
      Max : Real := 0.0;
      J_up : Natural;
      Aii  : Real;
      Ajj  : Real;
      Aod  : Real;
      Asq  : Real;
      Eps  : constant Real := 1.0e-16;
      diffr : Real;
      sign  : Real;
      tden  : Real;
      Tank  : Real;
      C    : Real;
      S      : Real;
      xj     : Real;
   begin
      -- initialize unit matrix
      Item.Unit_Matrix := (Others => (Others => 0.0));
      for I in 1..Item.Atoms loop
         Item.Unit_Matrix(I, I) := 1.0;
      end loop;
      for I in 2..Item.Atoms loop
         J_Up := I - 1;
         for J in 1..J_Up loop
            Aii := Item.Atomic_Diagonal(I);
            Ajj := Item.Atomic_Diagonal(J);
            Aod := Item.Atomic_Matrix(J, I);
            Asq := Aod * Aod;
            if Asq > Max then
               Max := Asq;
            end if;
            if Asq > EPS then
               diffr := Aii - Ajj;
               if Diffr < 0.0 then
                  Sign := -2.0;
                  diffr := -Diffr;
               else
                  Sign := 2.0;
               end if;
               Tden := Diffr + sqrt(Diffr * Diffr + 4.0 * Asq);
               tank := Sign * Aod / Tden;
               C := 1.0 / sqrt(1.0 + Tank * tank);
               S := C * Tank;
               for K in 1..Item.Atoms loop
                  xj := c * Unit_Matrix(J,K) - S * Item.Unit_Matrix(I,K);
                  Item.Unit_Matrix(I, K) := S * Item.Unit_Matrix(J,K) +
                    C * Item.Unit_Matrix(I, K);
                  Item.Unit_Matrix(J, K) := Xj;
                  if K < J then
                     Xj := C * Item.Atomic_Matrix(k, j) -
                             S * Item.Atomic_Matrix(K, i);
                     Item.Atomic_Matrix(K, I) := S * Item.Atomic_Matrix(K,J) +
                       C * Item.Atomic_Matrix(k,I);
                  elsif K > J then
                     if K < I then
                        xj := C * Item.Atomic_Matrix(J,K) -
                          S * Item.Atomic_Matrix(K, I);
                        Item.Atomic_Matrix(K, I) := S * Item.Atomic_Matrix(J,K)
                          + C * Item.Atomic_Matrix(k, I);
                        Item.Atomic_Matrix(J,K) := xj;
                     elsif K > I then
                        Xj := C * Item.Atomic_Matrix(J,K) -
                          S * Item.Atomic_Matrix(I,K);
                        Item.Atomic_Matrix(I,k) := S * Item.Atomic_Matrix(J,K)
                          + C * Item.Atomic_Matrix(I,K);
                        Item.Atomic_Matrix(J,K) := Xj;
                     end if;
                  end if;
               end loop;
               Item.Atomic_Diagonal(I) := C * C * Aii + S * S * Ajj +
                 2.0 * S * C * Aod;
               Item.Atomic_Diagonal(J) := C * C * Ajj + S * S * Aii -
                 2.0 * S * C * Aod;
            end if;
         end loop;
      end loop;

   end Scofi1;

   -----------
   -- Order --
   -----------

   procedure Order(Item : in out Molecule) is
      atest  : Real;
      jtest  : Positive;
      umtest : Real;
   begin
      for K in 1..Item.Atoms loop
         atest := Item.Atomic_Diagonal(K);
         Jtest := K;
         -- find the minimum diagonal value in the range k..Item.Atoms
         for J in K..Item.Atoms loop
            if Item.Atomic_Diagonal(J) < atest then
               atest := Item.Atomic_Diagonal(J);
               Jtest := J;
            end if;
            Item.Atomic_Diagonal(J) := Item.Atomic_Diagonal(K);
            Item.Atomic_Diagonal(K) := atest;
         end loop;
         for I in Item.Unit_Matrix'Range(1) loop
            umtest := Item.Unit_Matrix(Jtest, I);
            Item.Unit_Matrix(Jtest, I) := Item.Unit_Matrix(K, I);
            Item.Unit_Matrix(K,I) := umtest;
         end loop;
      end loop;
   end Order;

   ----------
   -- Pmat --
   ----------

   procedure Pmat(Item : in out Molecule) is
      Sum : Real;
   begin
      for R in 1..Item.Atoms loop
         for S in 1..R loop
            Sum := 0.0;
            for J in 1..Item.Atoms loop
               Sum := Sum + Item.Bond_Orders(J, R) * Item.Bond_Orders(J, S);
            end loop;
            Item.Bond_Orders(S,R) := 2.0 * Sum;
         end loop;
      end loop;
   end Pmat;

   ----------
   -- Fval --
   ----------

   procedure Fval(Item : in out Molecule) is
      xnr : Real;
   begin
      for J in 1..Item.Atoms loop
         xnr := 0.0;
         for I in 1..Item.Atoms loop
            if I < J and then Item.Atomic_Matrix(I,J) + 0.1 <= 0.0 then
               xnr := xnr + Item.Bond_Orders(I, J);
            elsif Item.Atomic_Matrix(I,J) + 0.1 < 0.0 then
               xnr := xnr + Item.Bond_Orders(J,I);
            end if;
         end loop;
         Item.Free_Valences(J) := 1.732 - xnr;
      end loop;
   end Fval;

   ----------------------
   -- Compute_Energies --
   ----------------------

   procedure Compute_Energies (Item : in out Molecule) is
   begin
      Pahy(Item);
      Modify(Item);
      Scofi1(Item);
      Order(Item);
      Pmat(Item);
      Fval(Item);
   end Compute_Energies;


   ------------
   -- Output --
   ------------

   procedure Output (Item : in Molecule) is
      Ind  : Boolean := False;
      L    : Natural := 0;
      Nupp : Natural;
      Klow : Natural;
      Kupp : Natural;
      Sum  : Real;
      type Colcount is mod 10;
      Cols : Colcount := 0;
   begin
      while not Ind loop
         L := L + 1;
         if 10 * L >= Item.Atoms then
            Nupp := Item.Atoms - 10 * (L - 1);
            Ind := True;
         else
            Nupp := 10;
         end if;
         New_Line;
         Put_Line("Energy Levels");
         Klow := 10 * (L - 1) + 1;
         Kupp := Nupp + 10 * (L - 1);
         Put("    J=   ");
         for K in Klow..Kupp loop
            put(Item => K, Width => 10);
         end loop;
         New_Line(2);
         for K in Klow..kupp loop
            Put(Item => Item.Atomic_Diagonal(K),
              aft => 4, Fore => 6, Exp => 0);
         end loop;
         New_Line(2);
         Put_Line("Huckel Orbitals");
         Put("    J=   ");
         for K in Klow..Kupp loop
            put(Item => K, Width => 10);
         end loop;
         New_Line(2);
         for I in 1..Item.Atoms loop
            Put(Item => I, Width => 4);
            for K in Klow..kupp loop
               Put(Item => Item.Unit_Matrix(K,I),
                 Aft => 6, fore => 4, Exp => 0);
            end loop;
            New_Line;
         end loop;
      end loop;
      Sum := 0.0;
      for J in 1..Item.Orbitals loop
         Sum := Sum + Item.Atomic_Diagonal(J);
      end loop;
      Sum := 2.0 * Sum;
      New_Line;
      Put("Total PI-Electron Energy =");
      Put(Item => Sum, Aft => 4, Fore => 6, Exp => 0);
      New_Line(2);
      Put_Line("Charge Densities");
      for I in 1..Item.Atoms loop
         if Cols = 0 then
            New_Line;
         end if;
         Cols := Cols + 1;
         Put(Item => Item.Bond_Orders(I,I), Aft => 4, Fore => 6, Exp => 0);
      end loop;
      New_Line(2);
      Put_Line("Free Valences");
      Cols := 0;
      for I in 1..Item.Atoms loop
         if Cols = 0 then
            New_Line;
         end if;
         Cols := Cols + 1;
         Put(Item => Item.Free_Valences(I), Aft => 4, Fore => 6, Exp => 0);
      end loop;
      New_Line(2);
      Put_Line("Bond-Order Matrix");
      Cols := 0;
      for I in 1..Item.Atoms loop
         for J in 1..I loop
            if Cols = 0 then
               New_Line;
            end if;
            Cols := Cols + 1;
            Put(Item => Item.Bond_Orders(J, I), Aft => 4, Fore => 6, Exp => 0);
         end loop;
      end loop;

   end Output;

end Huckel;

Tested on x86_64-pc-linux-gnu, committed on trunk

2017-09-06  Ed Schonberg  <schonberg@adacore.com>

	* sem_util.adb (Needs_One_Formal): The first formal of such a
	function must be a controlling formal, so that Obj.F (X, Y)
	can have the interpretation F(Obj)(X, Y).
	* sem_util.ads: Clarify documentation.
diff mbox series

Patch

Index: sem_util.adb
===================================================================
--- sem_util.adb	(revision 251762)
+++ sem_util.adb	(working copy)
@@ -17050,6 +17050,7 @@ 
       if Ada_Version >= Ada_2005
         and then Present (First_Formal (E))
         and then No (Default_Value (First_Formal (E)))
+        and then Is_Controlling_Formal (First_Formal (E))
       then
          Formal := Next_Formal (First_Formal (E));
          while Present (Formal) loop
Index: sem_util.ads
===================================================================
--- sem_util.ads	(revision 251760)
+++ sem_util.ads	(working copy)
@@ -2012,9 +2012,10 @@ 
    --  entity E. If no such instance exits, return Empty.
 
    function Needs_One_Actual (E : Entity_Id) return Boolean;
-   --  Returns True if a function has defaults for all but its first
-   --  formal. Used in Ada 2005 mode to solve the syntactic ambiguity that
-   --  results from an indexing of a function call written in prefix form.
+   --  Returns True if a function has defaults for all but its first formal,
+   --  which is a controlling formal. Used in Ada 2005 mode to solve the
+   --  syntactic ambiguity that results from an indexing of a function call
+   --  that returns an array, so that Obj.F (X, Y) may mean F (Ob) (X, Y).
 
    function New_Copy_List_Tree (List : List_Id) return List_Id;
    --  Copy recursively an analyzed list of nodes. Uses New_Copy_Tree defined