text:
# Class: ComputationRegister
#
# Description:
#  This class implements a register of the computations (see
#  RigidMotionsParameterSpaceDecompostion) which is based on SQLite.
#  Thanks to the register the computations can be restored after a crash.
#
# Author:
#  Kacper Pluta - kacper.pluta@esiee.fr
#  Laboratoire d'Informatique Gaspard-Monge - LIGM, A3SI, France
#
# Date:
#  11/12/2016
#
# License:
#  Simplified BSD License
#
# Copyright (c) 2016, Kacper Pluta
# All rights reserved.

# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#   * Redistributions of source code must retain the above copyright
#     notice, this list of conditions and the following disclaimer.
#   * Redistributions in binary form must reproduce the above copyright
#     notice, this list of conditions and the following disclaimer in the
#     documentation and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL Kacper Pluta BE LIABLE FOR ANY
# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
module ComputationRegister()
  option object;

  (* Connection object returned by Database[SQLite]:-Open *)
  local connection;
  local version;

# Comments:
#  We should use Database[SQLite]:-Bind() with care. Values passed should be given as local or
#  global but never as values returned from other procedures.

# Method: ModuleCopy
#   Standard constructor / copy constructor
#
# Parameters:
#   self::ComputationRegister      - a new object to be constructed
#   proto::ComputationRegister     - a prototype object from which self is derived
#   dbPath::string                 - a path to a database
#
# Comment:
#   The database is open in such a way that journal_mode is set to WAL and
#   synchronous to NORMAL. These option should not be changed!
#   See SQLite documentation for more information. Note that second database
#   is created in memory to speed up inserts. In particular, small chunks of
#   data are first inserted into the memory database and then moved in one
#   transaction into the main database.
#
# Output:
#   An object of type ComputationRegister.
#
# Exceptions:
#  "There is no database" occurs when a given file (dbPath) does not exists.
#
  export ModuleCopy::static := proc( self::ComputationRegister,
                                     proto::ComputationRegister,
                                     dbPath::string, $ )
  local fileStatus := false, version, stmt, openStatus := true;
    if _passed = 2 then
      self:-connection := proto:-connection;
      self:-version := proto:-version;
    else
      fileStatus:=FileTools:-Exists(dbPath);
      # wait until DB is not busy
      while openStatus do
        try
          openStatus := false;
          self:-connection := Database[SQLite]:-Open(dbPath);
        catch "database is locked":
          openStatus := true;
        end try;
      od;
      Database[SQLite]:-Attach(self:-connection, ":memory:", "cacheDB");
      Database[SQLite]:-Execute(self:-connection, "PRAGMA synchronous = NORMAL;");
      Database[SQLite]:-Execute(self:-connection, "PRAGMA journal_mode = WAL;");
      Database[SQLite]:-Execute(self:-connection, "PRAGMA cacheDB.auto_vacuum = FULL;");

      #Create tables
      if not fileStatus then
        Database[SQLite]:-Execute(self:-connection,"CREATE TABLE Quadric (ID PRIMARY KEY NOT " ||
        "NULL UNIQUE, polynom TEXT NOT NULL UNIQUE check(length(polynom) > 0));");
        Database[SQLite]:-Execute(self:-connection,"CREATE TABLE RealAlgebraicNumber (ID " ||
        "INTEGER PRIMARY KEY UNIQUE, polynom TEXT NOT NULL check(length(polynom) > 0), " ||
        "IntervalL TEXT NOT NULL check(length(IntervalL) > 0), IntervalR " ||
        "TEXT NOT NULL check(length(IntervalR) > 0));");
        Database[SQLite]:-Execute(self:-connection,"CREATE TABLE Events (RANumID INTEGER " ||
        "REFERENCES RealAlgebraicNumber (ID) ON DELETE CASCADE ON UPDATE CASCADE NOT NULL, " ||
        "QuadID INTEGER REFERENCES Quadric (ID) ON DELETE CASCADE ON UPDATE CASCADE NOT NULL);");
        Database[SQLite]:-Execute(self:-connection,"CREATE TABLE SamplePoint (A TEXT NOT NULL " ||
        "check(length(A) > 0), B TEXT NOT NULL check(length(B) > 0), C TEXT NOT NULL " ||
        "check(length(C) > 0));");
        Database[SQLite]:-Execute(self:-connection, "CREATE TABLE ComputedNumbers (RANumID " ||
        "INTEGER REFERENCES RealAlgebraicNumber (ID) ON DELETE CASCADE ON UPDATE CASCADE NOT " ||
        "NULL);");
      fi;

      stmt := Database[SQLite]:-Prepare(self:-connection, "PRAGMA user_version;");
      while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
      self:-version := Database[SQLite]:-Fetch(stmt, 0);
      Database[SQLite]:-Finalize(stmt);

      if self:-version = 0 then
        Database[SQLite]:-Execute(self:-connection,"CREATE TABLE cacheDB.RealAlgebraicNumber (ID " ||
        "INTEGER PRIMARY KEY UNIQUE, polynom TEXT NOT NULL check(length(polynom) > 0), " ||
        "IntervalL TEXT NOT NULL check(length(IntervalL) > 0), IntervalR TEXT NOT NULL " ||
        "check(length(IntervalR) > 0));");
        Database[SQLite]:-Execute(self:-connection, "CREATE TABLE cacheDB.Quadric (ID PRIMARY KEY " ||
        "NOT NULL UNIQUE, polynom TEXT NOT NULL UNIQUE check(length(polynom) > 0));");
        Database[SQLite]:-Execute(self:-connection, "CREATE TABLE cacheDB.Events (RANumID INTEGER " ||
        "REFERENCES RealAlgebraicNumber (ID) ON DELETE CASCADE ON UPDATE CASCADE NOT NULL, " ||
        "QuadID INTEGER REFERENCES Quadric (ID) ON DELETE CASCADE ON UPDATE CASCADE NOT NULL);");
        Database[SQLite]:-Execute(self:-connection, "CREATE TABLE cacheDB.SamplePoint (" ||
        "A TEXT NOT NULL check(length(A) > 0), B TEXT NOT NULL check(length(B) > 0), C TEXT " ||
        "NOT NULL check(length(C) > 0));");
      elif self:-version = 1 then
        Database[SQLite]:-Execute(self:-connection,"CREATE TABLE cacheDB.SamplePointSignature " ||
        "(SP_ID NOT NULL, Signature TEXT NOT NULL UNIQUE check(length(Signature) > 0));");
      elif self:-version >= 2 then
        Database[SQLite]:-Execute(self:-connection, "CREATE TABLE cacheDB.NMM (NMM TEXT NOT NULL " ||
        "UNIQUE check(length(NMM) > 0), SP_ID NOT NULL, T1 TEXT NOT NULL " ||
        "check(length(T1) > 0), T2 TEXT NOT NULL check(length(T2) > 0), T3 TEXT NOT NULL " ||
        "check(length(T3) > 0));");
      fi;
    fi;
    return self;
  end proc;

  export Close::static := proc( self::ComputationRegister )
    Database[SQLite]:-Close(self:-connection);
  end proc;

  export DatabaseVersion::static := proc( self::ComputationRegister )
    return self:-version;
  end proc;

# Method: ModulePrint
#   Standard printout of an object of type ComputationRegister.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#
  export ModulePrint::static := proc( self::ComputationRegister )
    print(self:-connection);
  end proc;

# Method: InsertQuadric
#   Used to insert quadrics into a database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#   id::integer                    - an identifier of a given quadric
#   quadric::polynom               - a second degree polynomial
#
# Comments:
#   Each quadric is inserted into a cache, memory stored, database.
#   SynchronizeQuadrics has to be called to move inserted polynomials
#   into the register.
#
  export InsertQuadric::static := proc(self::ComputationRegister, id::integer, quadric::polynom)
    local stmt, poly := sprintf("%a", quadric);
    if self:-version > 0 then
      error "Adding new quadrics is blocked! Re-run computations with a new database.";
    fi;
    stmt := Database[SQLite]:-Prepare(self:-connection,"INSERT OR IGNORE INTO " ||
                                 "cacheDB.Quadric(ID, polynom) VALUES (?, ?);");
    Database[SQLite]:-Bind(stmt, 1, id);
    Database[SQLite]:-Bind(stmt, 2, poly);
    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    Database[SQLite]:-Finalize(stmt);
  end proc;

# Method: InsertEvent
#   Used to insert quadrics into a database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#   idNum::integer                 - an identifier of a given event
#   event::EventType               - a given event
# Comments:
#   Each event is inserted into a cache, memory stored, database.
#   SynchronizeEvents has to be called to move inserted events
#   into the register.
#
  export InsertEvent::static := proc(self::ComputationRegister, idNum::integer,
                                     event::EventType)
    local x::integer, num, quadrics, stmt, poly, interA, interB;
    if self:-version > 0 then
      error "Adding new events is blocked! Re-run computations with a new database.";
    fi;
    num := GetRealAlgebraicNumber(event); quadrics := GetQuadrics(event);
    poly := sprintf("%a", GetPolynomial(num));
    interA := sprintf("%a", GetInterval(num)[1]);
    interB := sprintf("%a", GetInterval(num)[2]);
    stmt := Database[SQLite]:-Prepare(self:-connection,"INSERT OR IGNORE INTO " ||
            "cacheDB.RealAlgebraicNumber(ID, polynom, IntervalL, IntervalR) " ||
                                             "VALUES (?, ?, ?, ?);");
    Database[SQLite]:-Bind(stmt, 1, idNum);
    Database[SQLite]:-Bind(stmt, 2, poly);
    Database[SQLite]:-Bind(stmt, 3, interA);
    Database[SQLite]:-Bind(stmt, 4, interB);

    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    Database[SQLite]:-Finalize(stmt);

    for x in quadrics do
      stmt := Database[SQLite]:-Prepare(self:-connection,"INSERT OR IGNORE INTO " ||
                             "cacheDB.Events(RANumID, QuadID) VALUES(?, ?);");
      Database[SQLite]:-Bind(stmt, 1, idNum);
      Database[SQLite]:-Bind(stmt, 2, x);
     while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
     Database[SQLite]:-Finalize(stmt);
    od;
  end proc;

# Method: InsertSignature
#   Used to insert unique signatures of rotational sample points into a database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#   id::integer                    - an identifier of a given rotational sample point
#   sig::string                    - a given signature -- ordered indices of the critical planes
#
  export InsertSignature::static := proc(self::ComputationRegister, id::integer, sig::string)
    local stmt;
    if self:-version <> 1 then
      return NULL;
    fi;
    stmt := Database[SQLite]:-Prepare(self:-connection,"INSERT OR IGNORE INTO " ||
            "cacheDB.SamplePointSignature(SP_ID, Signature) VALUES (?, ?);");
    Database[SQLite]:-Bind(stmt, 1, id);
    Database[SQLite]:-Bind(stmt, 2, sig);

    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    Database[SQLite]:-Finalize(stmt);
  end proc;

# Method: SynchronizeSamplePointsSignatures
#   Used to synchronize a cache database with the database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#
  export SynchronizeSamplePointsSignatures::static := proc(self::ComputationRegister)
    local stmt;
    if self:-version <> 1 then
      return NULL;
    fi;
    stmt := Database[SQLite]:-Prepare(self:-connection,"INSERT INTO SamplePointSignature " ||
            "(SP_ID, Signature) SELECT * FROM cacheDB.SamplePointSignature AS SC WHERE NOT " ||
            "EXISTS(SELECT 1 FROM SamplePointSignature AS S WHERE SC.Signature = S.Signature);");
    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    Database[SQLite]:-Finalize(stmt);

    #clean up cacheDB
    stmt := Database[SQLite]:-Execute(self:-connection,"DELETE FROM cacheDB.SamplePointSignature;");
end proc;

# Method: CloseSignaturesAddition
#   Used to mark that signatures were computed.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#
  export CloseSignaturesAddition::static := proc(self::ComputationRegister)
    if self:-version <> 1 then
      return NULL;
    fi;

    Database[SQLite]:-Execute(self:-connection, "PRAGMA user_version = 2;");
    self:-version := 2;
end proc;

# Method: InsertNMM
#   Used to insert a nieghborhood motion map into the database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#   ID::integer                    - an id of the rotational sample point
#   NMM::list                      - a list of lists which represent a neighborhood motion map
#   T::list                        - a list or rationals which represent a translational sample
#                                    point
#
  export InsertNMM::static := proc(self::ComputationRegister, ID::integer, NMM::list, T::list)
    local stmt, NMMString, transX, transY, transZ;
    if self:-version < 2 then
      return NULL;
    fi;
    stmt := Database[SQLite]:-Prepare(self:-connection,"INSERT OR IGNORE INTO " ||
                                             "cacheDB.NMM(SP_ID, NMM, T1, T2, T3) VALUES (?, ?, " ||
                                             "?, ?, ?);");
    NMMString := sprintf("%a", NMM); transX := sprintf("%a", T[1]); transY := sprintf("%a", T[2]);
    transZ := sprintf("%a", T[3]);
    Database[SQLite]:-Bind(stmt, 1, ID);
    Database[SQLite]:-Bind(stmt, 2, NMMString);
    Database[SQLite]:-Bind(stmt, 3, transX);
    Database[SQLite]:-Bind(stmt, 4, transY);
    Database[SQLite]:-Bind(stmt, 5, transZ);

    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    Database[SQLite]:-Finalize(stmt);
  end proc;

# Method: SynchronizeNMM
#   Used to synchronize a cache database with the database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#
  export SynchronizeNMM::static := proc(self::ComputationRegister)
    local stmt;
    if self:-version < 2 then
      return NULL;
    fi;
    stmt := Database[SQLite]:-Prepare(self:-connection,"INSERT INTO NMM (NMM, SP_ID, T1, T2, " ||
            "T3) SELECT * FROM cacheDB.NMM AS CNM WHERE NOT EXISTS(SELECT 1 FROM NMM AS NM " ||
                                                                     "WHERE NM.NMM = CNM.NMM);");
    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    Database[SQLite]:-Finalize(stmt);

    #clean up cacheDB
    stmt := Database[SQLite]:-Execute(self:-connection,"DELETE FROM cacheDB.NMM;");
  end proc;

# Method: SynchronizeQuadrics
#   Synchronize quadrics between memory, cache, database and a given database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#
#
  export SynchronizeQuadrics::static := proc(self::ComputationRegister)
    local stmt;
    if self:-version > 0 then
      return NULL;
    fi;
    stmt := Database[SQLite]:-Prepare(self:-connection,"INSERT INTO Quadric SELECT * FROM " ||
            "cacheDB.Quadric as qc WHERE NOT EXISTS(SELECT 1 FROM Quadric AS Q WHERE " ||
                                                           "q.POLYNOM = qc.POLYNOM);");
    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    Database[SQLite]:-Finalize(stmt);
  end proc;

# Method: SynchronizeEvents
#   Synchronize events between memory, cache, database and a given database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#
# Comments:
#   Cache database is cleared. Therefore, the method should be called only when all events were
#   inserted.
  export SynchronizeEvents::static := proc(self::ComputationRegister)
    local stmt;
    if self:-version > 0 then
      return NULL;
    fi;
    stmt := Database[SQLite]:-Prepare(self:-connection,"INSERT INTO RealAlgebraicNumber " ||
                                "SELECT * FROM cacheDB.RealAlgebraicNumber AS rc WHERE NOT " ||
                                "EXISTS(SELECT 1 FROM RealAlgebraicNumber AS R WHERE " ||
                                "r.ID = rc.ID AND r.POLYNOM = rc.POLYNOM AND " ||
                                "r.INTERVALL = rc.INTERVALL AND " ||
                                "r.INTERVALR = rc.INTERVALR);");
    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    Database[SQLite]:-Finalize(stmt);

    stmt := Database[SQLite]:-Prepare(self:-connection,"INSERT INTO Events SELECT * FROM " ||
                           "cacheDB.Events AS ev WHERE NOT EXISTS( SELECT 1 FROM Events AS E " ||
                           "WHERE e.RANUMID = ev.RANUMID AND e.QUADID = ev.QUADID);");
    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    Database[SQLite]:-Finalize(stmt);

    #clean up cacheDB
    stmt := Database[SQLite]:-Execute(self:-connection,"DELETE FROM cacheDB.RealAlgebraicNumber;");
  end proc;

# Method: InsertComputedNumber
#    Inserts ids' of computed events into a given database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#   id::integer                    - an integer which stands for the identifier of a cluster
#
#
  export InsertComputedNumber::static := proc(self::ComputationRegister, id::integer)
    local stmt;
    if self:-version > 0 then
      error "Adding new computed events is blocked! Re-run computations with a new database.";
    fi;
    stmt := Database[SQLite]:-Prepare(self:-connection,"INSERT OR IGNORE INTO " ||
                                 "ComputedNumbers (RANumID) VALUES (?);");
    Database[SQLite]:-Bind(stmt, 1, id);
    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    Database[SQLite]:-Finalize(stmt);
  end proc;

# Method: InsertSamplePoint
#    Inserts a sample point into a given database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#   samp::list                     - a list of three rational numbers which represent a sample point
#
#
  export InsertSamplePoint::static := proc(self::ComputationRegister, samp::list)
    local stmt, a, b, c;
    if self:-version > 0 then
      error "Adding new computed sample points is blocked! Re-run computations with a new database.";
    fi;
    a := sprintf("%a", samp[1]); b := sprintf("%a", samp[2]); c := sprintf("%a", samp[3]);
    stmt := Database[SQLite]:-Prepare(self:-connection,"INSERT INTO cacheDB.SamplePoint " ||
                                                       "(A, B, C) VALUES (?, ?, ?);");
    Database[SQLite]:-Bind(stmt, 1, a);
    Database[SQLite]:-Bind(stmt, 2, b);
    Database[SQLite]:-Bind(stmt, 3, c);

    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    Database[SQLite]:-Finalize(stmt);
  end proc;

# Method: SynchronizeSamplePoints
#    Synchronize sample points between cache, memory, database and a given database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#
# Comments:
#   Cache database is cleared.
  export SynchronizeSamplePoints::static := proc(self::ComputationRegister)
    local stmt;
    if self:-version > 0 then
      return NULL;
    fi;
    stmt := Database[SQLite]:-Prepare(self:-connection, "INSERT INTO SamplePoint SELECT * " ||
                                                                 "FROM cacheDB.SamplePoint;");
    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    Database[SQLite]:-Finalize(stmt);

    #clean up cacheDB
    stmt := Database[SQLite]:-Execute(self:-connection,"DELETE FROM cacheDB.SamplePoint;");
  end proc;

# Method: FetchComputedNumbers
#    Fetch ids' of events which were processed already from the database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#
#
# Output:
#   A list of integers.
#
  export FetchComputedNumbers::static := proc(self::ComputationRegister)
    local stmt := Database[SQLite]:-Prepare(self:-connection, "SELECT * FROM ComputedNumbers;");
    local result := convert(Database[SQLite]:-FetchAll(stmt), list);
    Database[SQLite]:-Finalize(stmt);
    return result;
  end proc;

# Method: FetchQuadrics
#    Fetch quadrics from the database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#
#
# Output:
#   A list of second degree polynomials.
#
  export FetchQuadrics::static := proc(self::ComputationRegister)
    local stmt := Database[SQLite]:-Prepare(self:-connection, "SELECT polynom FROM Quadric " ||
                                            "ORDER BY ID;");
    local result := map(parse, convert(Database[SQLite]:-FetchAll(stmt), list));
    Database[SQLite]:-Finalize(stmt);
    return result;
  end proc;

export FetchLowerEventID := proc(self::ComputationRegister)
    local lowerID;
    local stmt := Database[SQLite]:-Prepare(self:-connection, "SELECT MIN(ID) FROM RealAlgebraicNumber;");
    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    lowerID := Database[SQLite]:-Fetch(stmt, 0);
    Database[SQLite]:-Finalize(stmt);
    return lowerID;
end proc;

# Method: FetchEvents
#    Fetch events from the database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#   first::integer                 - an id of the first event from a given range
#   last::integer                  - an id of the last event from a given range
#
#
# Output:
#   An Array of EventType.
#
  export FetchEvents::static := proc(self::ComputationRegister, first::integer, last::integer)
    local row, events:=Array([]);
    local stmt := Database[SQLite]:-Prepare(self:-connection, "SELECT ID, polynom, IntervalL, " ||
    "IntervalR, group_concat(QuadID) FROM RealAlgebraicNumber JOIN EVENTS " ||
    "ON ID = RANUMID WHERE ID BETWEEN ? AND ? GROUP BY ID ORDER BY ID;");
    Database[SQLite]:-Bind(stmt, 1, first);
    Database[SQLite]:-Bind(stmt, 2, last);

    #Slow but Fetching all can kill with memory consumption
    while Database[SQLite]:-Step(stmt) <> Database[SQLite]:-RESULT_DONE do
      row := Database[SQLite]:-FetchRow(stmt);
      events(row[1]) := EventType(RealAlgebraicNumber(parse(row[2]), parse(row[3]),
                    parse(row[4])), [parse(row[5])]);
    od;
    Database[SQLite]:-Finalize(stmt);
    return events;
  end proc;

# Method: FetchSamplePointsWithoutSignature
#   Used to fetch a given number of sample points from the database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#   first::integer                 - an id of a sample point which represent the beginning of the
#                                    range of ids to be fetched
#   last::integer                  - an id of a sample point which represent the end of the range of
#                                    ids to be fetched
# Comment:
#   Only sample points without signatures are fetch.
#
  export FetchSamplePointsWithoutSignature::static := proc(self::ComputationRegister,
                                                       first::integer, last::integer)
    local row, buffer:=Array([]), stmt;
     if self:-version < 1 then
      error "This version of FetchSamplePoints requires a database in version 1.";
    fi;
    stmt := Database[SQLite]:-Prepare(self:-connection, "SELECT ID, A, B, C FROM SamplePoint AS S " ||
            "WHERE S.ID BETWEEN ? AND ? AND NOT EXISTS( SELECT 1 FROM SamplePointSignature AS SS " ||
                                                                         "WHERE S.ID = SS.SP_ID);");
    Database[SQLite]:-Bind(stmt, 1, first);
    Database[SQLite]:-Bind(stmt, 2, last);

    #Slow but Fetching all can kill with memory consumption
    while Database[SQLite]:-Step(stmt) <> Database[SQLite]:-RESULT_DONE do
      row := Database[SQLite]:-FetchRow(stmt);
      ArrayTools:-Append(buffer, [row[1], op(map(parse, row[2..()]))], inplace=true);
    od;
    Database[SQLite]:-Finalize(stmt);
    return buffer;
  end proc;

# Method: FetchTopologicallyDistinctSamplePoints
#   Used to fetch a given number of topologically different sample points from the database. These
#   are the sample points which induce different order of the critical planes in the remainder range
#   -- they have different signatures.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#   first::integer                 - an id of a sample point which represent the beginning of the
#                                    range of ids to be fetched
#   last::integer                  - an id of a sample point which represent the end of the range of
#                                    ids to be fetched
#
  export FetchTopologicallyDistinctSamplePoints::static := proc(self::ComputationRegister,
                                                            first::integer, last::integer)
    local row, buffer:=Array([]), stmt;
     if self:-version < 1 then
      error "This version of FetchTopologicallyDistinctSamplePoints requires a database " ||
            "in version 1.";
    fi;
    stmt := Database[SQLite]:-Prepare(self:-connection, "SELECT S.ID, A, B, C FROM " ||
    "SamplePoint AS S, SamplePointSignature AS SSP WHERE S.ID = SSP.SP_ID AND SSP.ID " ||
    "BETWEEN ? AND ?;");
    Database[SQLite]:-Bind(stmt, 1, first);
    Database[SQLite]:-Bind(stmt, 2, last);

    #Slow but Fetching all can kill with memory consumption
    while Database[SQLite]:-Step(stmt) <> Database[SQLite]:-RESULT_DONE do
      row := Database[SQLite]:-FetchRow(stmt);
      ArrayTools:-Append(buffer, [row[1], op(map(parse, row[2..()]))], inplace=true);
    od;
    Database[SQLite]:-Finalize(stmt);
    return buffer;
  end proc;

# Method: NumberOfEvents
#   Returns number of events in the database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#
# Output:
#   Number of events in the database.
#
  export NumberOfEvents::static := proc(self::ComputationRegister)
    local stmt, num;
    stmt := Database[SQLite]:-Prepare(self:-connection, "SELECT COUNT(ID) " ||
                                            "FROM RealAlgebraicNumber;");
    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    num := Database[SQLite]:-Fetch(stmt, 0);
    Database[SQLite]:-Finalize(stmt);
    return num;
  end proc;

# Method: NumberOfEvents
#   Updates the database from version 0 into 1 by adding tables which stores neighborhood motion
#   maps, translational sample points and signatures of rotational sample points. After the updated
#   it is not possible to add new events or rotational sample points into the database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#
  export PrepareSamplePoints::static := proc(self::ComputationRegister)
    local stmt, toCompute;
    stmt := Database[SQLite]:-Prepare(self:-connection,"SELECT COUNT(ID) FROM RealAlgebraicNumber "
    || "WHERE ID NOT IN (SELECT RANUMID FROM ComputedNumbers) AND ID NOT IN (SELECT MAX(ID) " ||
    "FROM RealAlgebraicNumber);");
    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    toCompute := Database[SQLite]:-Fetch(stmt, 0);
    Database[SQLite]:-Finalize(stmt);

    if self:-version = 0 and toCompute = 0 then
      Database[SQLite]:-Execute(self:-connection, "ALTER TABLE SamplePoint RENAME TO " ||
      "sqlitestudio_temp_table; CREATE TABLE SamplePoint (A TEXT NOT NULL check(length(A) > 0), " ||
      "B TEXT NOT NULL check(length(B) > 0), C TEXT NOT NULL check(length(C) > 0), " ||
      "ID INTEGER PRIMARY KEY AUTOINCREMENT); INSERT INTO SamplePoint (A, B, C) SELECT A, B, C " ||
      "FROM sqlitestudio_temp_table; DROP TABLE sqlitestudio_temp_table;");
      Database[SQLite]:-Execute(self:-connection, "PRAGMA user_version = 1;");
      self:-version := 1;
      Database[SQLite]:-Execute(self:-connection, "CREATE TABLE SamplePointSignature (ID INTEGER " ||
      "PRIMARY KEY AUTOINCREMENT, SP_ID REFERENCES SamplePoint (ID) ON DELETE CASCADE ON UPDATE " ||
      "CASCADE NOT NULL, Signature TEXT NOT NULL UNIQUE check(length(Signature) > 0));");
      Database[SQLite]:-Execute(self:-connection, "CREATE TABLE cacheDB.SamplePointSignature " ||
      "(SP_ID NOT NULL, Signature TEXT NOT NULL UNIQUE check(length(Signature) > 0));");
      Database[SQLite]:-Execute(self:-connection,"CREATE TABLE NMM (ID INTEGER PRIMARY KEY " ||
      "AUTOINCREMENT, NMM TEXT NOT NULL UNIQUE check(length(NMM) > 0), SP_ID INTEGER " ||
      "REFERENCES SamplePoint (ID) ON DELETE CASCADE ON UPDATE CASCADE, T1 TEXT NOT NULL " ||
      "check(length(T1) > 0), T2 TEXT NOT NULL check(length(T2) > 0), T3 TEXT NOT NULL " ||
      "check(length(T3) > 0));");
      Database[SQLite]:-Execute(self:-connection,"CREATE TABLE cacheDB.NMM (NMM TEXT NOT NULL " ||
      "UNIQUE check(length(NMM) > 0), SP_ID NOT NULL, T1 TEXT NOT NULL check(length(T1) > 0), " ||
      "T2 TEXT NOT NULL check(length(T2) > 0), T3 TEXT NOT NULL check(length(T3) > 0));");
    elif toCompute <> 0 then
      Close(self);
      error "Before running computation of NMM it is necessary to compute all sample points! " ||
      "Please, run first RigidMotionsParameterSpaceDecompostion:-LaunchResumeComputations().";
    fi;
  end proc;

# Method: NumberOfSamplePoints
#   Returns number of the rotational sample points in the database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#
# Output:
#   Number of the rotational sample points in the database.
#
  export NumberOfSamplePoints::static := proc(self::ComputationRegister)
    local stmt, num::integer;
    stmt := Database[SQLite]:-Prepare(self:-connection, "SELECT COUNT(*) FROM SamplePoint;");
    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    num::integer := Database[SQLite]:-Fetch(stmt, 0);
    Database[SQLite]:-Finalize(stmt);
    return num;
  end proc;

# Method: NumberOfTopologicallyDistinctSamplePoints
#   Returns number of the topologically distinct rotational sample points in the database.
#
# Parameters:
#   self::ComputationRegister      - an instance of ComputationRegister
#
# Output:
#   Number of the topologically distinct sample points in the database.
#
  export NumberOfTopologicallyDistinctSamplePoints::static := proc(self::ComputationRegister)
    local stmt, num::integer;
    stmt := Database[SQLite]:-Prepare(self:-connection, "SELECT COUNT(*) FROM SamplePoint AS S, " ||
               "SamplePointSignature AS SSP WHERE S.ID = SSP.SP_ID;");
    while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
    num::integer := Database[SQLite]:-Fetch(stmt, 0);
    Database[SQLite]:-Finalize(stmt);
    return num;
  end proc;

  export DropRedundantSamplePoints::static := proc(self::ComputationRegister)
    local stmt;
    if self:-version = 2 then
      Database[SQLite]:-Execute(self:-connection, "CREATE TABLE SamplePoints (A TEXT NOT NULL " ||
      "CHECK (length(A) > 0), B TEXT NOT NULL CHECK (length(B) > 0), C TEXT NOT NULL CHECK " ||
      "(length(C) > 0), ID INTEGER PRIMARY KEY);");

      # Copy only topologically unique sample points
      stmt := Database[SQLite]:-Prepare(self:-connection, "INSERT INTO SamplePoints SELECT A, B," ||
      "C, S.ID FROM SamplePoint AS S, SamplePointSignature AS SSP WHERE S.ID = SSP.SP_ID;");
      while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
      Database[SQLite]:-Finalize(stmt);

      Database[SQLite]:-Execute(self:-connection, "CREATE TABLE SamplePointSignatures (ID INTEGER "
      || "PRIMARY KEY AUTOINCREMENT, SP_ID REFERENCES SamplePoints (ID) ON DELETE CASCADE ON " ||
      "UPDATE CASCADE NOT NULL, Signature TEXT NOT NULL UNIQUE CHECK (length(Signature) > 0));");

      # Copy data
      stmt := Database[SQLite]:-Prepare(self:-connection, "INSERT INTO SamplePointSignatures "||
      "SELECT * FROM SamplePointSignature;");
      while Database[SQLite]:-Step(stmt) = Database[SQLite]:-RESULT_BUSY do; od;
      Database[SQLite]:-Finalize(stmt);

      Database[SQLite]:-Execute(self:-connection, "DROP TABLE SamplePointSignature; " ||
      "DROP TABLE SamplePoint; VACUUM;");

     # Recover old names
     Database[SQLite]:-Execute(self:-connection, "ALTER TABLE SamplePoints RENAME TO " ||
     "SamplePoint;");

     Database[SQLite]:-Execute(self:-connection, "ALTER TABLE SamplePointSignatures RENAME TO " ||
     "SamplePointSignature;");

      Database[SQLite]:-Execute(self:-connection, "PRAGMA user_version = 3;");
      self:-version := 3;
    fi;
  end proc;
end module;

max_stars_repo_name: copyme/MapleTools

max_issues_repo_name: copyme/MapleTools

max_forks_repo_name: copyme/MapleTools


text:
Foo:=module()
    export fun1,fun2;
    `fun1`:=proc()
        fun11:=proc()

        end proc:
        return fun11;
    end proc:
    fun2:=proc()
        return module()
        end module:
    end proc:
end module:

A[B]:=proc()

end proc:

A[1]:=proc()

end proc:

A["fadfa"]:=proc()

end proc:

:-C[D]:=proc()

end proc:

:-U['V']:=proc()

end proc:

A:-B:-C:=proc()

end proc:

U[V][W]:=proc()
    fadfa:=proc()

    end proc:

end proc:

jkladfjakdflkja:=proc()

end proc:

(*
    OptMethods.mpl
*)

# Maple 符号线性规划
Opt[SLP]:=proc(A,b)
    uses LinearAlgebra;
    local E,x,FF,AA,BB,m,n;
    n:=RowDimension(A);
    m:=ColumnDimension(A);
    E:=IdentityMatrix(n);
    FF:=Vector([[0$m][],[1$n][]]);
    AA:=<<A|-E>,<-A|-E>>;
    BB:=convert(<b,-b>,Vector);
    x:=Opt[Simplex](FF,AA,BB);
    x:=x[1..m];
    return x;
end proc:

# Maple 单纯形法矩阵接口
# min f*x
# A*x<=b
Opt[Simplex]:=proc(f,A,b)
    uses LinearAlgebra;
    local m,n,v,V,x,AA,CC;
    m:=RowDimension(A);
    n:=ColumnDimension(A);
    V:=Vector([seq(v[k],k=1..n)]);
    AA:=A.V;
    CC:=map(k->(AA[k]<=b[k]),{seq(1..m)});
    x:=simplex[minimize](f^%T.V,CC,UNRESTRICTED);
    x:=map(k->rhs(k),[x[]]);
    return Vector(x);
end proc:

# Maple 进行整数线性规划
Opt[NILP]:=proc(A,b)
    uses LinearAlgebra,Optimization;
    local E,x,FF,AA,BB,m,n;
    n:=RowDimension(A);
    m:=ColumnDimension(A);
    E:=IdentityMatrix(n);
    FF:=Vector([[0$m][],[1$n][]]);
    AA:=<<A|-E>,<-A|-E>>;
    # 根据帮助文档的说法，sparse 矩阵更快
    # 实际上并不是 -_-||
    # AA:=Matrix(AA,storage=sparse);
    BB:=convert(<b,-b>,Vector);
    x:=LPSolve(FF,[AA,BB],integervariables=[seq(1..m)]);
    x:=map(round,x[2][1..m]);
    return x;
end proc:

# Matlab 整数线性规划
# 调用代价约为 0.02 秒
Opt[MatlabILP]:=proc(A,b)
    uses Matlab;
    local x;
    setvar("A",A);
    setvar("b",b);
    # 直接运行具体代码，避免产生路径问题
    evalM("[n,m]=size(A);E=eye(n);AA=[A,-E;-A,-E];BB=[b;-b];f=[zeros(m,1);ones(n,1)];x=intlinprog(f,1:m,AA,BB,[],[],[],[],optimoptions('intlinprog','Display','off'));x=x(1:m);");
    x:=getvar("x");
    if type(x,float) then
        x:=Vector([round(x)]);
    else
        x:=map(round,x);
    end if;
    return x;
end proc:

# Mathematica 符号线性规划
# 调用代价约为 0.7 秒
Opt[mmaSLP]:=proc(A,b)
    uses LinearAlgebra;
    local E,x,FF,AA,BB,m,n;
    n:=RowDimension(A);
    m:=ColumnDimension(A);
    E:=IdentityMatrix(n);
    FF:=Vector([[0$m][],[1$n][]]);
    AA:=<<A|-E>,<-A|-E>>;
    BB:=convert(<b,-b>,Vector);
    x:=Opt[mmaLP](FF,-AA,-BB);
    x:=x[1..m];
    return x;
end proc:

Opt[mmaLP]:=proc(f,A,b)
    uses FileTools;
    local cmd,fd,r;
    cmd:=sprintf("Print[LinearProgramming[%s,%s,%s,-Infinity,Method->\"Simplex\"]];",
        mmaVec(f),mmaMat(A),mmaVec(b));
    fd:=Text[Open]("LP.wl",create,overwrite);
    Text[WriteString](fd,cmd);
    Text[Close](fd);
    r:=ssystem("wolframscript -file LP.wl")[2];
    try
        r:=parse(cat("[",r[2..-2],"]"));
    catch :
        error r;
    finally
        Remove("LP.wl");
    end try;
    return Vector(r);# end proc
end proc:# proc

(*
mmaVec:=proc(x)
    uses StringTools;
    local k,n,buffer,fmt;
    n:=numelems(x);
    buffer:=StringBuffer();
    buffer:-clear();
    buffer:-append("{");
    for k from 1 to n do
        buffer:-append(sprintf("%d",x[k]));
        if k<n then
            buffer:-append(",");
        end if;
    end do;
    buffer:-append("}");
    return buffer:-value('clear');
end proc:
*)

# mmaMat:=proc(A)
#     uses LinearAlgebra,StringTools;
#     local k,m,buffer;
#     m:=RowDimension(A);
#     buffer:=StringBuffer();
#     buffer:-clear();
#     buffer:-append("{");
#     for k from 1 to m do
#         buffer:-append(mmaVec(A[k]));
#         if k<m then
#             buffer:-append(",");
#         end if;
#     end do;
#     buffer:-append("}");
#     return buffer:-value('clear');
# end proc:

printf("\d \d \" fadfa hahah:=proc()end proc \" \"\"\n");
printf("hahah:=proc()end proc\n");

fadsfasdf:=proc()
    return proc()

    end proc
end proc:

fadf:=proc()
    fadsfafadfasd:=module()
    end:
end:

gjhdg:=module()

end:

ext: mpl


text:
flavor 1
srclang 3
id 65535
numfuncs 2
import "triangle.mplt"
import "/home/bravewtz/Desktop/openarkcompiler/libjava-core/libjava-core.mplt"
entryfunc &Ltriangle_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V
fileinfo {
  @INFO_filename "triangle.jar"}
srcfileinfo {
  1 "triangle.java"}
javaclass $Ltriangle_3B <$Ltriangle_3B> public
func &Ltriangle_3B_7C_3Cinit_3E_7C_28_29V public constructor (var %_this <* <$Ltriangle_3B>>) void
func &Ltriangle_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V public static (var %Reg10_R743 <* <[] <* <$Ljava_2Flang_2FString_3B>>>>) void
var $_C_STR_2000 fstatic <[3] u64> readonly = [0, 0x300000000, 0x2000000020]
var $Ljava_2Flang_2FSystem_3B_7Cout extern <* <$Ljava_2Fio_2FPrintStream_3B>> final public static
var $__cinf_Ljava_2Flang_2FString_3B <$__class_meta__>
func &MCC_GetOrInsertLiteral () <* <$Ljava_2Flang_2FString_3B>>
func &Ltriangle_3B_7C_3Cinit_3E_7C_28_29V public constructor (var %_this <* <$Ltriangle_3B>>) void {
  funcid 48153
  var %Reg1_R43694 <* <$Ltriangle_3B>>
  var %Reg1_R57 <* <$Ljava_2Flang_2FObject_3B>>

  dassign %Reg1_R43694 0 (dread ref %_this)
  #INSTIDX : 0||0000:  aload_0
  #INSTIDX : 1||0001:  invokespecial
  dassign %Reg1_R57 0 (retype ref <* <$Ljava_2Flang_2FObject_3B>> (dread ref %Reg1_R43694))
  callassigned &Ljava_2Flang_2FObject_3B_7C_3Cinit_3E_7C_28_29V (dread ref %Reg1_R57) {}
  #INSTIDX : 4||0004:  return
  return ()
}
func &Ltriangle_3B_7Cmain_7C_28ALjava_2Flang_2FString_3B_29V public static (var %Reg10_R743 <* <[] <* <$Ljava_2Flang_2FString_3B>>>>) void {
  funcid 48154
  var %Reg0_I i32
  var %Reg6_I i32
  var %Reg0_R9160 <* <[] <* <[] i32>>>>
  var %Reg7_R12 ref
  var %Reg8_I i32
  var %Reg9_I i32
  var %Reg7_R9160 <* <[] <* <[] i32>>>>
  var %Reg0_R12 ref
  var %Reg1_I i32
  var %Reg1_R12 ref
  var %Reg2_I i32
  var %Reg2_R12 ref
  var %Reg0_R562 <* <$Ljava_2Fio_2FPrintStream_3B>>
  var %Reg1_R43 <* <$Ljava_2Flang_2FString_3B>>
  var %L_STR_15536 <* <$Ljava_2Flang_2FString_3B>>
  var %Reg1_R1112 <* <$Ljava_2Flang_2FStringBuilder_3B>>
  var %Reg3_I i32
  var %Reg2_R43 <* <$Ljava_2Flang_2FString_3B>>

  intrinsiccallwithtype <$Ltriangle_3B> JAVA_CLINIT_CHECK ()
  #INSTIDX : 0||0000:  bipush
  dassign %Reg0_I 0 (constval i32 10)
  #INSTIDX : 2||0002:  istore_1
  dassign %Reg6_I 0 (dread i32 %Reg0_I)
  #INSTIDX : 3||0003:  iload_1
  #INSTIDX : 4||0004:  iload_1
  #INSTIDX : 5||0005:  multianewarray
  dassign %Reg6_I 0 (cvt i32 i32 (dread i32 %Reg6_I))
  dassign %Reg6_I 0 (cvt i32 i32 (dread i32 %Reg6_I))
  #INSTIDX : 9||0009:  astore_2
  dassign %Reg7_R12 0 (dread ref %Reg0_R9160)
  #INSTIDX : 10||000a:  iconst_0
  dassign %Reg0_I 0 (constval i32 0)
  #INSTIDX : 11||000b:  istore_3
  dassign %Reg8_I 0 (dread i32 %Reg0_I)
  dassign %Reg8_I 0 (cvt i32 i32 (dread i32 %Reg8_I))
@label0   #INSTIDX : 12||000c:  iload_3
  #INSTIDX : 13||000d:  iload_1
  #INSTIDX : 14||000e:  if_icmpge
  brtrue @label6 (ge i32 i32 (dread i32 %Reg8_I, dread i32 %Reg6_I))
  #INSTIDX : 17||0011:  iconst_0
  dassign %Reg0_I 0 (constval i32 0)
  #INSTIDX : 18||0012:  istore 4
  dassign %Reg9_I 0 (dread i32 %Reg0_I)
  dassign %Reg9_I 0 (cvt i32 i32 (dread i32 %Reg9_I))
@label1   #INSTIDX : 20||0014:  iload 4
  #INSTIDX : 22||0016:  iload_3
  #INSTIDX : 23||0017:  if_icmpgt
  brtrue @label5 (gt i32 i32 (dread i32 %Reg9_I, dread i32 %Reg8_I))
  #INSTIDX : 26||001a:  iload 4
  #INSTIDX : 28||001c:  ifeq
  brtrue @label2 (eq i32 i32 (dread i32 %Reg9_I, constval i32 0))
  #INSTIDX : 31||001f:  iload_3
  #INSTIDX : 32||0020:  iload 4
  #INSTIDX : 34||0022:  if_icmpne
  brtrue @label3 (ne i32 i32 (dread i32 %Reg8_I, dread i32 %Reg9_I))
@label2   #INSTIDX : 37||0025:  aload_2
  #INSTIDX : 38||0026:  iload_3
  #INSTIDX : 39||0027:  aaload
  dassign %Reg7_R9160 0 (retype ref <* <[] <* <[] i32>>>> (dread ref %Reg7_R12))
  dassign %Reg0_R12 0 (iread ref <* ref> 0 (array 1 ptr <* <[] ref>> (dread ref %Reg7_R9160, dread i32 %Reg8_I)))
  #INSTIDX : 40||0028:  iload 4
  #INSTIDX : 42||002a:  iconst_1
  dassign %Reg1_I 0 (constval i32 1)
  #INSTIDX : 43||002b:  iastore
  iassign <* i32> 0 (
    array 1 ptr <* <[] i32>> (dread ref %Reg0_R12, dread i32 %Reg9_I),
    dread i32 %Reg1_I)
  #INSTIDX : 44||002c:  goto
  goto @label4
@label3   #INSTIDX : 47||002f:  aload_2
  #INSTIDX : 48||0030:  iload_3
  #INSTIDX : 49||0031:  aaload
  dassign %Reg7_R9160 0 (retype ref <* <[] <* <[] i32>>>> (dread ref %Reg7_R12))
  dassign %Reg0_R12 0 (iread ref <* ref> 0 (array 1 ptr <* <[] ref>> (dread ref %Reg7_R9160, dread i32 %Reg8_I)))
  #INSTIDX : 50||0032:  iload 4
  #INSTIDX : 52||0034:  aload_2
  #INSTIDX : 53||0035:  iload_3
  #INSTIDX : 54||0036:  iconst_1
  dassign %Reg1_I 0 (constval i32 1)
  #INSTIDX : 55||0037:  isub
  dassign %Reg1_I 0 (sub i32 (dread i32 %Reg8_I, dread i32 %Reg1_I))
  #INSTIDX : 56||0038:  aaload
  dassign %Reg7_R9160 0 (retype ref <* <[] <* <[] i32>>>> (dread ref %Reg7_R12))
  dassign %Reg1_R12 0 (iread ref <* ref> 0 (array 1 ptr <* <[] ref>> (dread ref %Reg7_R9160, dread i32 %Reg1_I)))
  #INSTIDX : 57||0039:  iload 4
  #INSTIDX : 59||003b:  iconst_1
  dassign %Reg2_I 0 (constval i32 1)
  #INSTIDX : 60||003c:  isub
  dassign %Reg2_I 0 (sub i32 (dread i32 %Reg9_I, dread i32 %Reg2_I))
  #INSTIDX : 61||003d:  iaload
  dassign %Reg1_I 0 (iread i32 <* i32> 0 (array 1 ptr <* <[] i32>> (dread ref %Reg1_R12, dread i32 %Reg2_I)))
  #INSTIDX : 62||003e:  aload_2
  #INSTIDX : 63||003f:  iload_3
  #INSTIDX : 64||0040:  iconst_1
  dassign %Reg2_I 0 (constval i32 1)
  #INSTIDX : 65||0041:  isub
  dassign %Reg2_I 0 (sub i32 (dread i32 %Reg8_I, dread i32 %Reg2_I))
  #INSTIDX : 66||0042:  aaload
  dassign %Reg7_R9160 0 (retype ref <* <[] <* <[] i32>>>> (dread ref %Reg7_R12))
  dassign %Reg2_R12 0 (iread ref <* ref> 0 (array 1 ptr <* <[] ref>> (dread ref %Reg7_R9160, dread i32 %Reg2_I)))
  #INSTIDX : 67||0043:  iload 4
  #INSTIDX : 69||0045:  iaload
  dassign %Reg2_I 0 (iread i32 <* i32> 0 (array 1 ptr <* <[] i32>> (dread ref %Reg2_R12, dread i32 %Reg9_I)))
  #INSTIDX : 70||0046:  iadd
  dassign %Reg1_I 0 (add i32 (dread i32 %Reg1_I, dread i32 %Reg2_I))
  #INSTIDX : 71||0047:  iastore
  iassign <* i32> 0 (
    array 1 ptr <* <[] i32>> (dread ref %Reg0_R12, dread i32 %Reg9_I),
    dread i32 %Reg1_I)
@label4   #INSTIDX : 72||0048:  iinc
  dassign %Reg9_I 0 (add i32 (dread i32 %Reg9_I, constval i32 1))
  dassign %Reg9_I 0 (cvt i32 i32 (dread i32 %Reg9_I))
  #INSTIDX : 75||004b:  goto
  goto @label1
@label5   #INSTIDX : 78||004e:  iinc
  dassign %Reg8_I 0 (add i32 (dread i32 %Reg8_I, constval i32 1))
  dassign %Reg8_I 0 (cvt i32 i32 (dread i32 %Reg8_I))
  #INSTIDX : 81||0051:  goto
  goto @label0
@label6   #INSTIDX : 84||0054:  iconst_0
  dassign %Reg0_I 0 (constval i32 0)
  #INSTIDX : 85||0055:  istore_3
  dassign %Reg8_I 0 (dread i32 %Reg0_I)
  dassign %Reg8_I 0 (cvt i32 i32 (dread i32 %Reg8_I))
@label7   #INSTIDX : 86||0056:  iload_3
  #INSTIDX : 87||0057:  iload_1
  #INSTIDX : 88||0058:  if_icmpge
  brtrue @label12 (ge i32 i32 (dread i32 %Reg8_I, dread i32 %Reg6_I))
  #INSTIDX : 91||005b:  iload_3
  #INSTIDX : 92||005c:  istore 4
  dassign %Reg9_I 0 (dread i32 %Reg8_I)
@label8   #INSTIDX : 94||005e:  iload 4
  #INSTIDX : 96||0060:  iload_1
  #INSTIDX : 97||0061:  if_icmpge
  brtrue @label9 (ge i32 i32 (dread i32 %Reg9_I, dread i32 %Reg6_I))
  #INSTIDX : 100||0064:  getstatic
  intrinsiccallwithtype <$Ljava_2Flang_2FSystem_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg0_R562 0 (dread ref $Ljava_2Flang_2FSystem_3B_7Cout)
  #INSTIDX : 103||0067:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_2000) { dassign %L_STR_15536 0 }
  dassign %Reg1_R43 0 (dread ptr %L_STR_15536)
  #INSTIDX : 105||0069:  invokevirtual
  virtualcallassigned &Ljava_2Fio_2FPrintStream_3B_7Cprint_7C_28Ljava_2Flang_2FString_3B_29V (dread ref %Reg0_R562, dread ref %Reg1_R43) {}
  #INSTIDX : 108||006c:  iinc
  dassign %Reg9_I 0 (add i32 (dread i32 %Reg9_I, constval i32 1))
  #INSTIDX : 111||006f:  goto
  goto @label8
@label9   #INSTIDX : 114||0072:  iload_1
  #INSTIDX : 115||0073:  iload_3
  #INSTIDX : 116||0074:  isub
  dassign %Reg6_I 0 (cvt i32 i32 (dread i32 %Reg6_I))
  dassign %Reg0_I 0 (sub i32 (dread i32 %Reg6_I, dread i32 %Reg8_I))
  #INSTIDX : 117||0075:  istore 4
  dassign %Reg9_I 0 (dread i32 %Reg0_I)
  dassign %Reg9_I 0 (cvt i32 i32 (dread i32 %Reg9_I))
@label10   #INSTIDX : 119||0077:  iload 4
  #INSTIDX : 121||0079:  iload_1
  #INSTIDX : 122||007a:  if_icmpgt
  brtrue @label11 (gt i32 i32 (dread i32 %Reg9_I, dread i32 %Reg6_I))
  #INSTIDX : 125||007d:  getstatic
  intrinsiccallwithtype <$Ljava_2Flang_2FSystem_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg0_R562 0 (dread ref $Ljava_2Flang_2FSystem_3B_7Cout)
  #INSTIDX : 128||0080:  new
  intrinsiccallwithtype <$Ljava_2Flang_2FStringBuilder_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg1_R1112 0 (gcmalloc ref <$Ljava_2Flang_2FStringBuilder_3B>)
  #INSTIDX : 131||0083:  dup
  #INSTIDX : 132||0084:  invokespecial
  callassigned &Ljava_2Flang_2FStringBuilder_3B_7C_3Cinit_3E_7C_28_29V (dread ref %Reg1_R1112) {}
  #INSTIDX : 135||0087:  aload_2
  #INSTIDX : 136||0088:  iload_3
  #INSTIDX : 137||0089:  aaload
  dassign %Reg7_R9160 0 (retype ref <* <[] <* <[] i32>>>> (dread ref %Reg7_R12))
  dassign %Reg2_R12 0 (iread ref <* ref> 0 (array 1 ptr <* <[] ref>> (dread ref %Reg7_R9160, dread i32 %Reg8_I)))
  #INSTIDX : 138||008a:  iload_1
  #INSTIDX : 139||008b:  iload 4
  #INSTIDX : 141||008d:  isub
  dassign %Reg6_I 0 (cvt i32 i32 (dread i32 %Reg6_I))
  dassign %Reg3_I 0 (sub i32 (dread i32 %Reg6_I, dread i32 %Reg9_I))
  #INSTIDX : 142||008e:  iaload
  dassign %Reg2_I 0 (iread i32 <* i32> 0 (array 1 ptr <* <[] i32>> (dread ref %Reg2_R12, dread i32 %Reg3_I)))
  #INSTIDX : 143||008f:  invokevirtual
  virtualcallassigned &Ljava_2Flang_2FStringBuilder_3B_7Cappend_7C_28I_29Ljava_2Flang_2FStringBuilder_3B (dread ref %Reg1_R1112, dread i32 %Reg2_I) { dassign %Reg1_R1112 0 }
  #INSTIDX : 146||0092:  ldc
  callassigned &MCC_GetOrInsertLiteral (addrof ptr $_C_STR_2000) { dassign %L_STR_15536 0 }
  dassign %Reg2_R43 0 (dread ptr %L_STR_15536)
  #INSTIDX : 148||0094:  invokevirtual
  virtualcallassigned &Ljava_2Flang_2FStringBuilder_3B_7Cappend_7C_28Ljava_2Flang_2FString_3B_29Ljava_2Flang_2FStringBuilder_3B (dread ref %Reg1_R1112, dread ref %Reg2_R43) { dassign %Reg1_R1112 0 }
  #INSTIDX : 151||0097:  invokevirtual
  virtualcallassigned &Ljava_2Flang_2FStringBuilder_3B_7CtoString_7C_28_29Ljava_2Flang_2FString_3B (dread ref %Reg1_R1112) { dassign %Reg1_R43 0 }
  #INSTIDX : 154||009a:  invokevirtual
  virtualcallassigned &Ljava_2Fio_2FPrintStream_3B_7Cprint_7C_28Ljava_2Flang_2FString_3B_29V (dread ref %Reg0_R562, dread ref %Reg1_R43) {}
  #INSTIDX : 157||009d:  iinc
  dassign %Reg9_I 0 (add i32 (dread i32 %Reg9_I, constval i32 1))
  dassign %Reg9_I 0 (cvt i32 i32 (dread i32 %Reg9_I))
  #INSTIDX : 160||00a0:  goto
  goto @label10
@label11   #INSTIDX : 163||00a3:  getstatic
  intrinsiccallwithtype <$Ljava_2Flang_2FSystem_3B> JAVA_CLINIT_CHECK ()
  dassign %Reg0_R562 0 (dread ref $Ljava_2Flang_2FSystem_3B_7Cout)
  #INSTIDX : 166||00a6:  invokevirtual
  virtualcallassigned &Ljava_2Fio_2FPrintStream_3B_7Cprintln_7C_28_29V (dread ref %Reg0_R562) {}
  #INSTIDX : 169||00a9:  iinc
  dassign %Reg8_I 0 (add i32 (dread i32 %Reg8_I, constval i32 1))
  dassign %Reg8_I 0 (cvt i32 i32 (dread i32 %Reg8_I))
  #INSTIDX : 172||00ac:  goto
  goto @label7
@label12   #INSTIDX : 175||00af:  return
  return ()
}

max_stars_repo_name: Bravewtz/OpenArkCompiler

max_issues_repo_name: Bravewtz/OpenArkCompiler

max_forks_repo_name: Bravewtz/OpenArkCompiler


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)

$include "gga_c_zvpbeint.mpl"
$include "gga_c_pbeloc.mpl"

params_a_alpha := 0.5:
params_a_omega := 2:

(* redefine nu of zbpbeint *)
zvpbeint_nu := (rs, z, t) ->
  2*(4/(3*Pi^2))^(1/18) * rs^(1/3):

(* Note that f_pbe here is, in fact, pbeloc *)
f  := (rs, z, xt, xs0, xs1) ->
  zvpbeint_ff(rs, z, 0) * f_pbe(rs, z, xt, xs0, xs1):


text:
impulse process ACV_confirmBlockOption {
  ACV_loadChunks()
  /clone ${origin + (1 1 -6)} ${origin + (7 2 -6)} 0 1 1
  /fill 0 1 2 6 2 2 air
  /tellraw @a [{"text":"[","color":"gold"},{"text":"Info","color":"gold","italic":"true"},{"text":"] ","color":"gold"},{"text":"Block-Option confirmed","color":"yellow"}]
}


text:
# Example from Section 5.3
read "../ComputeIdentifiableFunctions.mpl";

#------------------------------------------------------------------------------
# Auxiliary function for counting nonconstant coefficients

CountNonNumeric := proc(p, var)
    local result, c:
    result := 0:
    for c in coeffs(p, var) do
        if not type(c, numeric) then
            result := result + 1:
        end if:
    end do:
    return result:
end proc:

#------------------------------------------------------------------------------

model_original := [
    diff(x1(t), t) = (a * x1(t) - b * x1(t) * x2(t)),
    diff(x2(t), t) = (-c * x2(t) + d * x1(t) * x2(t)),
    diff(f(t), t) = 0,
    y1(t) = (e * x1(t) + f(t) * x2(t)),
    y2(t) = f(t)
]:

printf("We are interested in computing ME-identifiable functions and the bound for ME-identifiablility\n"):
ME_identifiable_functions := MultiExperimentIdentifiableFunctions(model_original, no_bound=true, simplified_generators=true)[3]:
printf("The field of multi-experiment identifiable functions is generated by %a\n", ME_identifiable_functions):

printf("To find a good bound for the number of experiments we consider the following model (as described in Section 5.3)\n"):
model := [
    diff(x1(t), t) = (a * x1(t) - b * x1(t) * x2(t)),
    diff(x2(t), t) = (-c * x2(t) + d * x1(t) * x2(t)),
    y(t) = (e * x1(t) + f * x2(t))
]:
print(model):

printf("First we compute the bound and IO-equations considering f just as a parameter\n"):
result := MultiExperimentIdentifiableFunctions(model, simplified_generators=true):
printf("If we consider a different system with f being an unknown parameter with the same value for different experiments, the bound from the theorem will be %a\n", result[1]):
# In the notation of Section 5.3
io_eq_coeffs := [expand(denom(result[2][1])), op(map(p -> expand(simplify(p * denom(result[2][1]))), result[2]))]:
print(io_eq_coeffs):
ns := sort(map(p -> CountNonNumeric(p, f), io_eq_coeffs)):
n := ns[1]:
printf("The number of non-numeric coefficients with respect to f (in the notation of Section 5.3: n, n1, n2, ...) in the coefficients of the only input-output equation are: %a\n", ns):
printf("Section 5.3 implies that the number of experiments for the case when f is a known parameter different in different experiments is %a\n", max(n + min(ns[2...]), max(ns[2..]))):

max_stars_repo_name: iliailmer/AllIdentifiableFunctions

max_issues_repo_name: iliailmer/AllIdentifiableFunctions

max_forks_repo_name: iliailmer/AllIdentifiableFunctions


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: mgga_exc *)

b97mv_par_n := 6:

b97mv_gamma_x := 0.004:
b97mv_par_x := [
    [  0.85,  0, 0],
    [  1.007, 0, 1],
    [  0.259, 1, 0],
    [  0,   0, 0],
    [  0,   0, 0],
    [  0,   0, 0]
]:

b97mv_gamma_ss := 0.2:
b97mv_par_ss := [
    [  0.443,  0, 0],
    [ -1.437,  0, 4],
    [ -4.535,  1, 0],
    [ -3.39,   2, 0],
    [  4.278,  4, 3],
    [  0,    0, 0]
]:

b97mv_gamma_os := 0.006:
b97mv_par_os := [
    [  1.000,  0, 0],
    [  1.358,  1, 0],
    [  2.924,  2, 0],
    [ -8.812,  2, 1],
    [ -1.39,   6, 0],
    [  9.142,  6, 1]
]:

$include "lda_x_erf.mpl"
$include "b97mv.mpl"

wb97mv_f := (rs, z, xs0, xs1, ts0, ts1) ->
   my_piecewise3(screen_dens_zeta(rs,  z), 0, (1 + z)/2 * lda_x_erf_spin(rs*(2/(1 + z))^(1/3),  1)
    * b97mv_g(b97mv_gamma_x,  b97mv_wx_ss, b97mv_par_x,  b97mv_par_n, xs0, ts0, 0))
+  my_piecewise3(screen_dens_zeta(rs, -z), 0, (1 - z)/2 * lda_x_erf_spin(rs*(2/(1 - z))^(1/3),  1)
    * b97mv_g(b97mv_gamma_x,  b97mv_wx_ss, b97mv_par_x,  b97mv_par_n, xs1, ts1, 0)):

f :=  (rs, z, xt, xs0, xs1, us0, us1, ts0, ts1) ->
  wb97mv_f(rs, z, xs0, xs1, ts0, ts1) +
  b97mv_f(rs, z, xs0, xs1, ts0, ts1):


text:
# Kinematik-Berechnung 3. Arm KAS5
# Beschreibung
# Berechnung der Kinematik der parallelen Struktur: Stellung des Kompensationsmechanismus in Abhängigkeit der gemessenen Winkel der Haupt-Armkinematik
# Zu Geometrie und Winkeln, siehe Bild (Unterlagen)
#
# Modell m3: Kurbel N7 ist starr mit Körper K2 gekoppelt (über die Achse)
#
# Berechne alle notwendigen Winkel in kintmp in Abhängigkeit konstanter Größen und der verallgemeinerten Koordinaten
#
# Berechne zuerst nur Sinus und Cosinus jedes Winkels. Diese werden für die Rotationsmatrizen benötigt.
# Die Berechnung der Winkel selbst (zur Kontrolle) benötigt die arctan-Funktion, die längere Ausdrücke nur sehr langsam auswertet.
# Daher werden zunächst alle Winkel mit %arctan (inert-Funktion berechnet) und am Ende erst ausgewertet.
# Die Zeitableitung der Winkel (in späteren Skripten) wird aus der inert-Funktion berechnet.
#
# Hinweise:
# Der Befehl arctan hat die Funktionalität von atan2
#
# Eingabe:
# "../codeexport/KAS5m3_definitions"
# * enthält die Definitionen aller Konstanten und verallgemeinerter Koordinaten
#
# Ausgabe:
# "../codeexport/%s/tmp/kinematic_constraints_maple_inert.m"
# * Erzeugt die Ausdrücke für Winkel der Parallelstruktur (kintmp), Ersetzungsausdrücke für sin,cos (kintmp_subsexp) und die Federlänge (lpar) in Abhängigkeit der verallgemeinerten Koordinaten.
# Quellen
# [1] Wang, Jiegao and Gosselin, ClémentM.: A New Approach for the Dynamic Analysis of Parallel Manipulators (1998)
# Autor
# Moritz Schappler, schappler@irt.uni-hannover.de, 2015-11
# Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
# Initialisierung
restart:
kin_constraints_exist := true: # Für Speicherung
;
with(StringTools): # Für Zeitausgabe
with(LinearAlgebra):
with(codegen):
with(CodeGeneration):
codegen_act := true:
codegen_opt := 1: # Geringerer Optimierungsgrad. Sonst zu lange.
codegen_debug := 0: # Zur Code-Generierung auch für Nicht-Inert-Ausdrücke
;
read "../helper/proc_MatlabExport":
read "../transformation/proc_rotx":
read "../transformation/proc_roty":
read "../transformation/proc_rotz":
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
read "../helper/proc_intersect_circle":
with(RealDomain): # Schränkt alle Funktionen auf den reellen Bereich ein. Muss nach Definition von MatlabExport kommen. Sonst geht dieses nicht.
;
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
# Variable mit Winkeln der Nebenstruktur nur in Abhängigkeit der verallgemeinerten Koordinaten
kintmp_qs := Matrix(RowDimension(kintmp_s),1):
# Konstante Winkel bereits hineinschreiben
for i from 1 to RowDimension(kintmp_s) do
  if diff(kintmp_s(i,1), t) = 0 then
    kintmp_qs(i,1) := kintmp_s(i,1):
  end if:
end do:
# Variablen definieren für die Hilfswinkel
gamma5_s := kintmp_s(5,1):
delta_s := kintmp_s(6..23,1):
delta_qs := kintmp_qs(6..23,1):
# Ersetzungsausdrücke definieren.
# Speichere Sinus und Cosinus der Winkel direkt ab, da diese in den Rotationsmatrizen direkt auftreten.
# Spalte 1: Zu suchender Ausdruck (sin oder cos eines Winkels)
# Spalte 2: Einzusetzender Ausdruck.
# Dadurch werden arctan-Ausdrücke in der direkten Kinematik reduziert.
# Ähnliches Vorgehen wie in [1].
kintmp_subsexp := Matrix(2*RowDimension(kintmp_s),2):
for i from 1 to RowDimension(kintmp_s) do
  kintmp_subsexp(2*i-1, 1) := sin(kintmp_s(i,1)):
  kintmp_subsexp(2*i,   1) := cos(kintmp_s(i,1)):
  # Initialisierung der rechten Spalte mit gleichen Werten. Später nur Ersetzung, wenn Vorteilhaft.
  kintmp_subsexp(2*i-1, 2) := kintmp_subsexp(2*i-1, 1):
  kintmp_subsexp(2*i,   2) := kintmp_subsexp(2*i,   1):
end do:
# Annahmen über Konstanten: Längen sind alle Positiv.
# (Wird anscheinend für arctan gebraucht). Ansonsten steht csgn in exportiertem Code.
assume(l1>0):
assume(l2>0):
assume(l3>0):
assume(l4>0):
assume(l5>0):
assume(l6>0):
assume(l7>0):
assume(l8>0):
assume(l9>0):
assume(l10>0):
assume(l11>0):
assume(l12>0):
assume(l13>0):
assume(l14>0):
assume(l15>0):
assume(l16>0):
assume(l17>0):
assume(l18>0):
assume(l19>0):
assume(l20>0):
assume(l21>0):
assume(l22>0):
# Schalter mit Optionen für die Code-Generierung.
# Wähle, ob trigonometrische arctan-Ausdrücke ausgewertet werden, oder ob Sinus- und Cosinussatz angewendet werden sollen.
# Dies hat großen Einfluss auf die Rechendauer.
# Die Optimierungen sind so gewählt, dass
# * Kreisschnittpunktsfunktionen nicht mit den komplizierten Eingabedaten aufgerufen, sondern mit allgemeinen Platzhaltern, die dann ersetzt werden,
# * Keine Zwischenoptimierungen durchgeführt werden. Durch simplify und combine werden die symbolischen Ausdrücke größer, was später die Rechenzeit enorm erhöht
# * Sinus und Cosinus der Parallelwinkel ersetzt werden und kein sin(arctan(...)) ausgerechnet wird.
Schalter_Opt := <1;1;0;0;0;0;0;1>:
printf("Beginn der Berechnungen. %s\n", FormatTime("%Y-%m-%d %H:%M:%S")):
st := time():
# Reihenfolge, in der die Winkel der Parallelstruktur berechnet werden. Hilft beim exportieren von Code, wenn beim Debuggen in der gleichen Reihenfolge ausgegeben wird.
Reihenfolge_kintmp := <12; 5; 21; 4; 9; 16; 3; 8; 24; 11; 2; 1; 6; 7; 10; 13; 14; 15; 17; 18; 19; 20; 22; 23>:
# Berechne delta7
# Gl. (30)
delta_qs(7) := kintmp_qs(22) + convert_t_s(theta(4)):
kintmp_qs(12) := delta_qs(7):
# Berechne gamma5
# Gl.(2)
gamma5_qs := - kintmp_qs(12) + Pi + convert_t_s(theta(5)):
kintmp_qs(5) := gamma5_qs:
kintmp_subsexp(9,2) := sin(gamma5_qs):
kintmp_subsexp(10,2) := cos(gamma5_qs):
# Berechne delta16
# Kreisschnittpunkt berechnen
r_4_G4_G5 := <l11;0;0>:
r_12_G5_B := <l14;0;0>:
R_4_12 := rotz(Pi+delta_qs(7)):
r_4_G5_B := R_4_12 . r_12_G5_B:
save r_4_G5_B, "../codeexport/KAS5m3_kinematic_constraints_r_4_G5_B_debug_maple":
R_3_4 := rotz(convert_t_s(theta(4))):
R_4_3 := Transpose(R_3_4):
# Vektor D-G4 (ausgedrückt in KS 3)
r_3_D_G4 := <l5-l6;0;0>:
r_4_G4_D := R_4_3 . (-r_3_D_G4):
# Vektor G4-G5 (ausgedrückt in KS 4)
r_4_G4_G5 := <l11;0;0>:
r_4_G5_D := -r_4_G4_G5 + r_4_G4_D:
save r_4_G5_D, "../codeexport/KAS5m3_kinematic_constraints_r_4_G5_D_noopt_debug_maple":
if Schalter_Opt(7) = 1 then
  r_4_G5_D := simplify(r_4_G5_D, trig):
  r_4_G5_D := simplify(r_4_G5_D):
  r_4_G5_D := combine(r_4_G5_D, trig):
end if:
save r_4_G5_D, "../codeexport/KAS5m3_kinematic_constraints_r_4_G5_D_opt_debug_maple":
if Schalter_Opt(8) = 1 then
  # Rufe die Kreisschnittpunkt-Funktion mit allgemeinen Werten auf und ersetze das Ergebnis.
  subs_ic := <r_4_G5_B_Allg_x; r_4_G5_B_Allg_y; r_4_G5_D_Allg_x; r_4_G5_D_Allg_y>:
  subs_ic_KAS := <r_4_G5_B(1..2,1); r_4_G5_D(1..2,1)>:
  ReturnValue_Allg := intersect_circle(<r_4_G5_B_Allg_x; r_4_G5_B_Allg_y>, <r_4_G5_D_Allg_x; r_4_G5_D_Allg_y>, l13, l22):
  ReturnValue := ReturnValue_Allg:
  for i from 1 to 4 do
    ReturnValue := subs({subs_ic(i) = subs_ic_KAS(i)}, ReturnValue):
  end do:
else
  ReturnValue := intersect_circle(r_4_G5_B, r_4_G5_D, l13, l22):
end if:
printf("Erster Kreisschnittpunkt berechnet. %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
r_4_G5_C := <ReturnValue(1..2, 2); 0>: # Alternative 2
;
save r_4_G5_C, "../codeexport/KAS5m3_kinematic_constraints_r_4_G5_C_debug_noopt_maple":
if Schalter_Opt(4) = 1 then
  r_4_G5_C := simplify(r_4_G5_C, trig):
  r_4_G5_C := simplify(r_4_G5_C):
  r_4_G5_C := combine(r_4_G5_C, trig):
end if:
printf("CPU-Zeit bis hier: %1.2fs\n", time()-st):
save r_4_G5_C, "../codeexport/KAS5m3_kinematic_constraints_r_4_G5_C_debug_opt_maple":
r_4_D_C := - r_4_G5_D + r_4_G5_C:
r_3_D_C := R_3_4 . r_4_D_C:
save r_3_D_C, "../codeexport/KAS5m3_kinematic_constraints_r_3_D_C_debug_noopt_maple":
if Schalter_Opt(6) = 1 then
  r_3_D_C := simplify(r_3_D_C, trig):
  r_3_D_C := simplify(r_3_D_C):
  r_3_D_C := combine(r_3_D_C, trig):
end if:
save r_3_D_C, "../codeexport/KAS5m3_kinematic_constraints_r_3_D_C_debug_opt_maple":
# Gl. (m3.5)
delta_qs(16) := %arctan(-r_3_D_C(2), r_3_D_C(1)):
sin_delta16_qs := -r_3_D_C(2)/l22:
cos_delta16_qs := r_3_D_C(1)/l22:
kintmp_qs(21) := delta_qs(16):
kintmp_subsexp(41,2) := sin_delta16_qs:
kintmp_subsexp(42,2) := cos_delta16_qs:
# Berechne gamma3
# Berechnungen für Hebel E-F-G3
# Zweiter Kreismittelpunkt für Schnittpunkt-Algorithmus
r_3_G3_D := <l6;0;0>:
r_3_D_G3 := -r_3_G3_D:
r_10_D_E := rotz(delta_qs(9)) . <l21;0;0>:
# Stelle Rotationsmatrix mit Zeitabhängigen Winkel delta16 auf (ohne Ersetzung mit delta16(q), da dort arctan vorkommt.
# Ersetze sin(delta16), cos(delta16) mit Ausdrücken von oben. Bei arctan(sin(...)) entstehen kompliziertere Ausdrücke
if Schalter_Opt(1) = 1 then
  R_3_10_s := rotz(- delta_s(16)):
  R_3_10 := subs_kintmp_exp(R_3_10_s):
else
  R_3_10 := rotz(- delta_qs(16)):
end if:
r_3_D_E := R_3_10 . r_10_D_E:
save r_3_D_E, sprintf("../codeexport/%s/tmp/kinematic_constraints_r_3_D_E_noopt_debug_maple", robot_name):
if Schalter_Opt(5) = 1 then
  r_3_D_E := simplify(r_3_D_E, trig):
  r_3_D_E := simplify(r_3_D_E):
  r_3_D_E := combine(r_3_D_E, trig):
end if:
save r_3_D_E, sprintf("../codeexport/%s/tmp/kinematic_constraints_r_3_D_E_opt_debug_maple", robot_name):
if Schalter_Opt(2) = 1 then
  # Rufe die Kreisschnittpunkt-Funktion mit allgemeinen Werten auf und ersetze das Ergebnis.
  # dies führt zu einer wesentlichen Verringerungen der Rechenzeit, da das Einsetzen und Weiterrechnen der komplizierten Ausdrücke in der Unterfunktion entfällt.
  subs_ic := <r_3_D_E_Allg_x; r_3_D_E_Allg_y; r_3_D_G3_Allg_x; r_3_D_G3_Allg_y>:
  subs_ic_KAS := <r_3_D_E(1..2,1); r_3_D_G3(1..2,1)>:
  ReturnValue_Allg := intersect_circle(<r_3_D_E_Allg_x; r_3_D_E_Allg_y>, <r_3_D_G3_Allg_x; r_3_D_G3_Allg_y>, l20, l4):
  ReturnValue := ReturnValue_Allg:
  for i from 1 to 4 do
    ReturnValue := subs({subs_ic(i) = subs_ic_KAS(i)}, ReturnValue):
  end do:
else
  ReturnValue := intersect_circle(r_3_D_E, r_3_D_G3, l20, l4):
end if:
printf("Zweiter Kreisschnittpunkt berechnet. %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
r_3_D_F := <ReturnValue(1..2, 2); 0>: # Alternative 2
;
save r_3_D_F, sprintf("../codeexport/%s/tmp/kinematic_constraints_r_3_D_F_noopt_debug_maple", robot_name):
if Schalter_Opt(3) = 1 then
  r_3_D_F := simplify(r_3_D_F, trig):
  r_3_D_F := simplify(r_3_D_F):
  r_3_D_F := combine(r_3_D_F, trig):
end if:
r_3_F_G3 := -r_3_D_F + r_3_D_G3:
printf("CPU-Zeit bis hier: %1.2fs\n", time()-st):
save r_3_D_F, sprintf("../codeexport/%s/tmp/kinematic_constraints_r_3_D_F_opt_debug_maple", robot_name):
save ReturnValue, r_3_D_F, r_3_F_G3, sprintf("../codeexport/%s/tmp/kinematic_constraints_gamma3_qs_debug_maple", robot_name):
# Gl. (m3.6)
gamma3_qs := %arctan(r_3_F_G3(2,1), -r_3_F_G3(1,1)):
sin_gamma3_qs := r_3_F_G3(2,1) / l4:
cos_gamma3_qs := (-r_3_F_G3(1,1)) / l4:
kintmp_qs(4) := gamma3_qs:
kintmp_subsexp(7,2) := sin_gamma3_qs:
kintmp_subsexp(8,2) := cos_gamma3_qs:
save ReturnValue, r_3_D_F, r_3_F_G3, gamma3_qs, sprintf("../codeexport/%s/tmp/kinematic_constraints_gamma3_qs_debug_maple", robot_name):
save kintmp_qs, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_debug_nach_gamma3.m", robot_name):
printf("Winkel gamma3 berechnet. %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
# Berechne delta4
# Gl. (22)
delta_qs(4) := Pi/2 - delta_qs(16) - delta_qs(12):
kintmp_qs(9) := delta_qs(4):
# Berechne delta11
# Gl. (3)
delta_qs(11) := delta_qs(10) - delta_qs(4):
kintmp_qs(16) := delta_qs(11):
# Berechne beta2
# Gl. (m3.7)
# Siehe auch: KAS5m3_sym_codegen_kinematic_constraints_hilfsberechnungen.mw
d7b2 := %arctan((sin(delta_qs(7)) * l14 + r_4_G5_C(2)), (cos(delta_qs(7)) * l14 + r_4_G5_C(1))):
beta2_qs := d7b2 - delta_qs(7):
kintmp_qs(3) := beta2_qs:
# Nutze Additionstheoreme zur Bestimmung von Sinus und Cosinus.
cos_d7b2 := (cos(delta_qs(7)) * l14 + r_4_G5_C(1)) / l13:
sin_d7b2 := (sin(delta_qs(7)) * l14 + r_4_G5_C(2)) / l13:
sin_beta2_qs := sin_d7b2*cos(delta_qs(7)) - cos_d7b2*sin(delta_qs(7)):
cos_beta2_qs := cos_d7b2*cos(delta_qs(7)) + sin_d7b2*sin(delta_qs(7)):
kintmp_subsexp(5,2) := sin_beta2_qs:
kintmp_subsexp(6,2) := cos_beta2_qs:
save kintmp_qs, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_debug_nach_beta2.m", robot_name):
# Berechne delta3
# Gl. (m3.9)
# Siehe auch: KAS5m3_sym_codegen_kinematic_constraints_hilfsberechnungen.mw
r_4_D_B := -r_4_G5_D + r_4_G5_B:
t7 := convert_t_s(theta(4)) + delta_qs(16):
# Nutze Additionstheoreme zur Bestimmung von Sinus und Cosinus.
# Dadurch kann Sinus/Cosinus von delta16 direkt genutzt werden und es tritt kein arctan auf.
cos_t7 := cos(convert_t_s(theta(4)))*cos_delta16_qs - sin(convert_t_s(theta(4)))*sin_delta16_qs:
sin_t7 := sin(convert_t_s(theta(4)))*cos_delta16_qs + cos(convert_t_s(theta(4)))*sin_delta16_qs:
# Hilfswinkel t8
t8 := %arctan((sin_t7 * l22 + r_4_D_B(2)), (cos_t7 * l22 - r_4_D_B(1))):
sin_t8 := (sin_t7 * l22 + r_4_D_B(2))/l13:
cos_t8 := (cos_t7 * l22 - r_4_D_B(1))/l13:
delta_qs(3) := t8 - t7:
cos_delta3_qs := (cos_t8)*(cos_t7) + (sin_t8)*(sin_t7):
sin_delta3_qs := (sin_t8)*(cos_t7) - (cos_t8)*(sin_t7):
kintmp_qs(8) := delta_qs(3):
kintmp_subsexp(15,2) := sin_delta3_qs:
kintmp_subsexp(16,2) := cos_delta3_qs:
save kintmp_qs, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_debug_nach_delta3.m", robot_name):
# Berechne rho3
# Gl. (11)
rho3_qs := gamma3_qs + delta_qs(8) - Pi/2:
# Nutze Additionstheoreme zur Bestimmung von Sinus und Cosinus.
# Dadurch kann Sinus/Cosinus von gamma3 direkt genutzt werden und es tritt kein arctan auf.
sin_rho3_qs := (sin_gamma3_qs)*(cos(delta_qs(8) - Pi/2)) + (cos_gamma3_qs)*(sin(delta_qs(8) - Pi/2)):
cos_rho3_qs := (cos_gamma3_qs)*(cos(delta_qs(8) - Pi/2)) - (sin_gamma3_qs)*(sin(delta_qs(8) - Pi/2)):
kintmp_qs(24) := rho3_qs:
kintmp_subsexp(47,2) := sin_rho3_qs:
kintmp_subsexp(48,2) := cos_rho3_qs:
save kintmp_qs, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_debug_nach_rho3.m", robot_name):
printf("Winkel rho3 berechnet. %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
# Berechne delta6
# Gl. (m3.8)
# Siehe auch: KAS5m3_sym_codegen_kinematic_constraints_hilfsberechnungen.mw
r_3_G3_E := r_3_G3_D + r_3_D_E:
t1 := -rho3_qs + delta_qs(8):
t2 := %arctan((-sin_t1 * l4 + r_3_G3_E(1)), (cos_t1 * l4 + r_3_G3_E(2))):
# Nutze Additionstheoreme zur Bestimmung von Sinus und Cosinus.
cos_t1 :=  cos_rho3_qs*cos(delta_qs(8)) + sin_rho3_qs*sin(delta_qs(8)):
sin_t1 := -sin_rho3_qs*cos(delta_qs(8)) + cos_rho3_qs*sin(delta_qs(8)):
delta_qs(6) := t1 + t2:
kintmp_qs(11) := delta_qs(6):
# Nutze Additionstheoreme zur Bestimmung von Sinus und Cosinus.
cos_t2 := ( cos_t1 * l4 + r_3_G3_E(2)) / l20:
sin_t2 := (-sin_t1 * l4 + r_3_G3_E(1)) / l20:
cos_delta6_qs := cos_t2*cos_t1 - sin_t2*sin_t1:
sin_delta6_qs := sin_t2*cos_t1 + cos_t2*sin_t1:
kintmp_subsexp(21,2) := sin_delta6_qs:
kintmp_subsexp(22,2) := cos_delta6_qs:
save kintmp_qs, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_debug_nach_delta6.m", robot_name):
# Kompensationsfeder
R_5_6:=rotz(convert_t_s(theta(6))):
r_6_A_H := <l17; 0; 0>:
r_5_A_H:=R_5_6 . r_6_A_H:
r_5_G5_G6:=<l12;0;0>:
r_6_G6_H:=<0; -l18; 0>:
r_5_G6_H:=R_5_6 . r_6_G6_H:
r_5_G5_H:= r_5_G5_G6 + r_5_G6_H: # Weg G5-G6-H
;
R_12_5_s:=rotz(gamma5_s): # Rotationsmatrix mit allgemeinem zeitabhängigem Ausdruck für den Winkel gamma5
;
R_12_5 := subs_kintmp_exp(R_12_5_s): # sinus und cosinus von gamma5 direkt ersetzen, ohne den Winkel einzusetzen.
;
R_5_12:=Transpose(R_12_5):
r_12_B_G5:=<-l14;0;0>:
r_5_G5_B:=R_5_12 . (-r_12_B_G5):
r_5_AB_2 := r_5_A_H - r_5_G5_H + r_5_G5_B: # Weg AB aus A-H-G5-B
;
l16_c := sqrt(r_5_AB_2[1,1]^2 + r_5_AB_2[2,1]^2): # Länge AB entspricht Federlänge
;
lpar_qs := l16_c:
printf("Kompensationsfederlänge berechnet. %s.  CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
# Berechne beta1
# Gl. (29)
# Siehe auch: KAS5_m4_sym_codegen_kinematic_constraints_hilfsberechnungen.mw
t2 := cos(convert_t_s(theta(6))):
t1 := sin(convert_t_s(theta(6))):
t3 := %arctan(-(t2 * l18 - sin(gamma5_qs) * l14 + t1 * l17), (t2 * l17 + cos(gamma5_qs) * l14 - t1 * l18 - l12)):
beta1_qs := t3 - gamma5_qs:
cos_t3 := (t2 * l17 + cos(gamma5_qs) * l14 - t1 * l18 - l12) / l16_c:
sin_t3 := (-(t2 * l18 - sin(gamma5_qs) * l14 + t1 * l17)) / l16_c:
kintmp_qs(2) := beta1_qs:
# Nutze Additionstheoreme zur Bestimmung von Sinus und Cosinus.
# Zu erledigen: In den Sinus/Cosinus-Ausdrücken kommt l16 vor. Muss noch generiert werden.
cos_beta1_qs := cos_t3*cos(gamma5_qs) + sin_t3*sin(gamma5_qs):
sin_beta1_qs := sin_t3*cos(gamma5_qs) - cos_t3*sin(gamma5_qs):
kintmp_subsexp(3,2) := sin_beta1_qs:
kintmp_subsexp(4,2) := cos_beta1_qs:
# Exportiere Code für folgende Skripte
# Entferne Annahmen aus Variablen
l1 := 'l1':
l2 := 'l2':
l3 := 'l3':
l4 := 'l4':
l5 := 'l5':
l6 := 'l6':
l7 := 'l7':
l8 := 'l8':
l9 := 'l9':
l10 := 'l10':
l11 := 'l11':
l12 := 'l12':
l13 := 'l13':
l14 := 'l14':
l15 := 'l15':
l16 := 'l16':
l17 := 'l17':
l18 := 'l18':
l19 := 'l19':
l20 := 'l20':
l21 := 'l21':
l22 := 'l22':
# Die Annahmen sind im Ausdruck bereits in den Variablen gespeichert. Ersetze das "~"-Zeichen.
# Quelle: http://www.mapleprimes.com/questions/207601-Remove-Assumptions
nms:=convert(indets(kintmp_qs,name),list): # Liste der Symbole
;
nmsS:=convert~(nms,string): # Liste der Symbole als String
;
L:=StringTools:-Substitute~(nmsS,"~",""): #Removing "~"
;
L1:=parse~(L):
S:=nms=~L1:
kintmp_qs:=subs(S,kintmp_qs):
lpar_qs:=subs(S,lpar_qs):
kintmp_subsexp:=subs(S,kintmp_subsexp):
kintmp_qt := convert_s_t(kintmp_qs):
lpar_qt := convert_s_t(lpar_qs):
# Speichere Maple-Ausdruck (Eingabe-Format und internes Format)
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple", robot_name):
save kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_subsexp_maple.m", robot_name):
printf("Ausdrücke für kintmp_subsexp gespeichert (Maple). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
for i from 1 to RowDimension(kintmp_s) do
  tmp := kintmp_qs(i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d", robot_name, i):
  save tmp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert_kintmpq_%d.m", robot_name, i):
end do:
save kin_constraints_exist, kintmp_qs, kintmp_qt, lpar_qs, lpar_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert", robot_name):
save kin_constraints_exist, kintmp_qs, kintmp_qt, lpar_qs, lpar_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
save kintmp_qs, sprintf("../codeexport/%s/tmp/kinematic_constraints_kintmp_qs_maple_inert", robot_name):
# Exportieren des vollständigen Ausdruckes
if codegen_act then
  MatlabExport(kintmp_subsexp, sprintf("../codeexport/KAS5_kinematik_parallel_kintmp_subsexp_matlab_noopt.m"), 0):
  MatlabExport(kintmp_subsexp, sprintf("../codeexport/KAS5_kinematik_parallel_kintmp_subsexp_matlab_opt.m"), 1):
end if:
printf("Ausdrücke mit Inert-Arctan exportiert (Matlab). %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( kintmp_subsexp(..,2) )):
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
MatlabExport(Transpose(kc_symbols), sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_matlab.m", robot_name), 2):
printf("Fertig. %s. CPU-Zeit bis hier: %1.2fs.\n", FormatTime("%Y-%m-%d %H:%M:%S"), time()-st):


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: mgga_exc *)
(* prefix:
  mgga_k_pc07_params *params;

  assert(p->params != NULL);
  params = (mgga_k_pc07_params * ) (p->params);
*)

pc07_p := x -> X2S^2*x^2:
pc07_q := u -> X2S^2*u:

(* Equation (15) redefined with decaying exponentials to avoid inf/inf situations *)
pc07_fab0 := z -> exp(-params_a_a*params_a_b/z) * (1+exp(-params_a_a/(params_a_a-z)))^params_a_b/(exp(-params_a_a/z) + exp(-params_a_a/(params_a_a-z)))^params_a_b:
(* The function is ill-behaving at z=0 and z=a.
   However, it also goes to 0 and 1 very quickly near these points.
   pc07_fab0(params_a_a/40) is ~ 1.6e-17, which is smaller than machine epsilon.
*)
pc07_thr := 1/40:
pc07_zlo := pc07_thr*params_a_a:
pc07_zhi := (1-pc07_thr)*params_a_a:
pc07_fab := z -> my_piecewise5(z<=pc07_zlo, 0, z>=pc07_zhi, 1, pc07_fab0( m_min(pc07_zhi, m_max(pc07_zlo, z)) ) ):

(* Equation (7) *)
pc07_Delta := (x, u) ->
  8*pc07_q(u)^2/81 - pc07_p(x)*pc07_q(u)/9 + 8*pc07_p(x)^2/243:

pc07_f_W    := x -> 5*pc07_p(x)/3:

(* Equation (8) *)
pc07_GE4  := (x, u) ->
  1 + 5*pc07_p(x)/27 + 20*pc07_q(u)/9 + pc07_Delta(x, u):

(* Equation (11) *)
pc07_GE4_M := (x, u) ->
  pc07_GE4(x, u)/sqrt(1 + pc07_Delta(x, u)^2/(1 + pc07_f_W(x))^2):

(* Equation (17) *)
pc07_f := (x, u) ->
  pc07_f_W(x) + (pc07_GE4_M(x, u) - pc07_f_W(x))*pc07_fab(pc07_GE4_M(x, u) - pc07_f_W(x)):

f := (rs, z, xt, xs0, xs1, u0, u1, t0, t1) ->
  mgga_kinetic(pc07_f, rs, z, xs0, xs1, u0, u1):


text:
######################################################################

`is_element/rel` := (A::set,B::set) -> proc(R)
 local AB,a,b;

 global reason;

 AB := [seq(seq([a,b],b in B),a in A)];

 if not type(R,set) then
  reason := [convert(procname,string),"R is not a set",R];
  return false;
 fi;

 if R minus {op(AB)} <> {} then
  reason := [convert(procname,string),"R is not a subset of AxB",R,A,B];
  return false;
 fi;

 return true;
end;

`is_equal/rel` := (A::set,B::set) -> proc(R,S)
 global reason;

 if R <> S then
  reason := [convert(procname,string),"R <> S",R,S];
  return false;
 fi;

 return true;
end;

`is_leq/rel` := (A::set,B::set) -> proc(R,S)
 global reason;

 if R minus S <> {} then
  reason := [convert(procname,string),"R is not a subset of S",R,S];
  return false;
 fi;

 return true;
end;

`bot/rel` := (A::set,B::set) -> {};
`top/rel` := proc(A::set,B::set)
 local a,b;
 return {seq(seq([a,b],b in B),a in A)};
end:

`is_a_function/rel` := (A::set,B::set) -> proc(R)
 local r,N;

 r := `hash/rel`(A,B)(R);
 N := map(a -> nops(r[a]),A) minus {1};
 return evalb(N = {});
end;

`is_total/rel` := (A::set,B::set) -> (R) -> evalb(R = `top/rel`(A,B));

`op/rel` := (A::set,B::set) -> (R) -> map(ab -> [ab[2],ab[1]],R);

`hash/rel` := (A::set,B::set) -> proc(R)
 local r,a,b,ab;

 r := table();
 for a in A do r[a] := {}; od;
 for ab in R do
  a,b := op(ab);
  r[a] := {op(r[a]),b};
 od;

 return eval(r);
end:

`unhash/rel` := (A::set,B::set) -> proc(r)
 local a,b;
 {seq(seq([a,b],b in r[a]),a in A)};
end;

`id/rel` := proc(A::set)
 local a;
 return {seq([a,a],a in A)};
end:

`o/rel` := (A::set,B::set,C::set) -> proc(S,R)
 local r,s,sr,a;
 r := `hash/rel`(A,B)(R);
 s := `hash/rel`(B,C)(S);
 sr := table();
 for a in A do sr[a] := map(op,map(b -> s[b],r[a])); od;
 return `unhash/rel`(A,C)(sr);
end;

`list_elements/rel` := proc(A::set,B::set)
 local AB,a,b;
 AB := {seq(seq([a,b],b in B),a in A)};
 `list_elements/subsets`(AB);
end:

`count_elements/rel` := (A::set,B::set) -> 2^(nops(A)*nops(B));

`random_element/rel` := (A::set,B::set) -> proc(p_)
 local AB,p,R,r,ab,a,b;
 p := `if`(nargs > 0,p_,0.5);
 r := rand(0..10^6-1);
 AB := {seq(seq([a,b],b in B),a in A)};
 R := NULL;
 for ab in AB do
  if r() < 10^6 * p then R := R,ab; fi;
 od:
 R := {R};
 return R;
end:


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: mgga_exc *)

$include "mgga_c_scan.mpl"

(* we need to redefine two functions *)
scan_e0_g := (rs, z, t) -> (1 + 8*A(rs, z, t)*t^2)^(-1/4)/2
       + (1 + 80*A(rs, z, t)^2*t^4)^(-1/8)/2:
scan_g_infty := s -> (1 + 4*scan_chi_infty*s^2)^(-1/4)/2
       + (1 + 80*scan_chi_infty^2*s^4)^(1/8)/2:

(* and the new parameters *)
scan_b1c :=  0.030197:
scan_b2c := -0.06623:
scan_b3c :=  0.16672:

(* set parameters of f_alpha *)
params_a_c1 := 1.131:
params_a_c2 := 1.7:
params_a_d  := 1.37:


text:
TITLE BetterBreadBakery;

MAX

    Profit = 0.05 Sun + 0.08 Moon - 4000/30 ;

SUBJECT TO

    1/10 Sun + 1/5 Moon <= 10 * 140 ;

    Sun  <=  8000;
    Moon <=  5000;
END

max_stars_repo_name: BjarniMax/mpl-models

max_issues_repo_name: BjarniMax/mpl-models

max_forks_repo_name: BjarniMax/mpl-models


text:
# This file is about Singer systems of endomorphisms of the additive formal
# group.  Such a system has groups G(i,j) (for 0 <= i <= j <= n say),
# all of which are equal to the standard additive formal group.  There
# are morphisms G(j,k) -> G(i,l) whenever i <= j <= k <= l, with the
# evident functoriality condition.  The horizontal maps are linear
# and the vertical maps are short isogenies.  Every equal-length
# L-shaped subdiagram is like G -> G/a <- G//a for some a.

aa := proc(i::nonnegint,k::nonnegint)
 if i > k then return FAIL; fi;
 return v[k] * mul(v[t]^(2^(k-t-1)),t=i..k-1);
end:

bb := proc(i::nonnegint,k::nonnegint)
 if i > k then return FAIL; fi;
 return v[i-1]^(2^(k-i));
end:

# G(i,k) to G(i-1,k)
pp := (i,k) -> (x) -> bb(i,k) * x:

# G(i,k) to G(i,k+1)
qq := (i,k) -> (x) -> aa(i,k) * x + x^2:

# G(j,k) to G(i,k)
ppp := proc(i,j,k)
 local x,y,z;

 if i > j then
  return FAIL;
 elif i = j then
  return (x -> x);
 else
  y := ppp(i+1,j,k)(x);
  z := collect(modp(expand(pp(i+1,k)(y)),2),x);
  return unapply(z,x);
 fi;
end:

# G(j,k) to G(i,k)
ppp1 := proc(i,j,k) local x; unapply(mul(v[m]^(2^(k-m-1)),m=i..j-1) * x,x); end:

# G(i,k) to G(i,l)
qqq := proc(i,k,l)
 local x,y,z;

 if k > l then
  return FAIL;
 elif l = k then
  return (x -> x);
 else
  y := qqq(i,k,l-1)(x);
  z := collect(modp(expand(qq(i,l-1)(y)),2),x);
  return unapply(z,x);
 fi;
end:

W0V := (k) -> mul(V[i],i=0..k-1):

WW := proc(k,l,m)
  local J,P;
  J := [seq(j,j=k..l-1)];
  P := combinat[choose](J,l-k-m);
  return add(mul(V[u[t]]^(2^(k+m+t-1-u[t])),t=1..l-k-m),u in P);
end:

WWW := (j,k,l,m) -> WW(k,l,m) * W0V(j)^(2^(k+m-j) - 2^(l-j)):

WWWW := (i,j,k,l,m) -> WW(k,l,m) * W0V(j)^(2^(k+m-j)) / W0V(i)^(2^(l-i)):

# G(i,k) to G(i,l)
qqq1 := proc(j,k,l)
 local x;
 return unapply(sort(collect(expand(add(WWW(j,k,l,i) * x^(2^i),i=0..l-k)),x)),x):
end:

# G(j,k) to G(i,l)
ff := proc(i,j,k,l)
 local x;
 return unapply(sort(collect(expand(add(WWWW(i,j,k,l,m) * x^(2^m),m=0..l-k)),x)),x):
end:

vV_rule := {seq(v[i]=V[i]/W0V(i),i=0..100)}:


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: lda_exc *)

a :=  0.0311:
b := -0.048:
c :=  0.009:
d := -0.017:

f := (rs, zeta) -> a*log(rs) + b + c*rs*log(rs) + d*rs:


text:
InvClassify:=module()
    option  package;
    export  doClassify,                     # 进行分类
            InvSol,                         # 解对象
            RepSol;                         # 代表元对象
    local   ModuleLoad;

$include "headers.mpl"

$include "InvSol.mpl"
$include "RepSol.mpl"

$include "Basic.mpl"
$include "Classifyer.mpl"
$include "Combine.mpl"
$include "Condition.mpl"
$include "Fetch.mpl"
$include "InvOrder.mpl"
$include "Interaction.mpl"
$include "InvSimplify.mpl"
$include "Logout.mpl"
$include "Utils.mpl"

    ModuleLoad:=proc()
        PDETools:-declare(quiet):
    end proc:
    ModuleLoad();

    # 进行分类
    doClassify:=proc(vv::list)
        local As,A,eqs,sols,reps;
        As,A,eqs:=getTransMatAndPDE(vv);
        classify(A,As,eqs);
        sols:=getSols();
        reps:=buildReps(sols);
        return reps;
    end proc:

end module:


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: work_gga_x *)

mbeta := 0.018641:
mcsi  := 8: (* for harmonic potentials *)

f := x -> 1 + mbeta/X_FACTOR_2D_C*x^2/(1 + mcsi*mbeta*x*arcsinh(x)):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
# If X ~ Exp(α) then e^−X ~ Beta(1/α, 1).

# Proof:

PDF_exponential := (alpha, x) -> exp(-x/alpha)/alpha;

PDF_beta := (beta, gamma, x) -> GAMMA(beta+gamma) * x^(beta-1) * (x - 1)^(gamma-1) / GAMMA(beta) / GAMMA(gamma);

simplify(PDF_exponential(1/lambda, -log(x))*diff(log(x), x));

simplify(PDF_beta(lambda, 1, x));

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


text:
Launch:=proc (
    NumNodes::posint := "ALL_NODES",
    code := "",
    Printer::procedure := NULL,
    CheckAbort::procedure := NULL,
    Exports::{list, set} := NULL,
    Imports::{list, set} := NULL,
    {
        allexternal::truefalse := true,
        checkabort::procedure := CheckAbort,
        clear::truefalse := true,
        exports::{list, set} := Exports,
        imports::{list, set} := Imports,
        numnodes::{"ALL_NODES", posint} := NumNodes,
        printer::procedure := Printer
    })
    localGridInterface:-Launch(numnodes, code, gridcmdargs = _rest, _options);
end proc:


text:
`is_element/EP/Fbar` := (N::posint) -> (A::set) -> (TT) -> proc(x)
 global reason;
 local TT1,C,P,Q,T;

 TT1 := `big_sets/trees`(TT);
 if not is_table_on(TT1)(x) then
  reason := [convert(procname,string),"x is not a table on TT1",x,TT1];
  return false;
 fi;

 C := children_map(TT);
 for T in TT1 do
  Q := C[T];
  P := select(U -> nops(U) > 1,Q);
  if not `is_element/D/Fbar`(N)(P,Q)(x[T]) then
   reason := [convert(procname,string),"x[T] is not in D/Fbar(N)(P,Q)",x[T],P,Q,reason];
   return false;
  fi;
 od;

 return true;
end:

######################################################################

`is_equal/EP/Fbar` := (N::posint) -> (A::set) -> (TT) -> proc(x,y)
 local TT1,C,P,Q,T,U;

 TT1 := `big_sets/trees`(TT);

 C := children_map(TT);
 for T in TT1 do
  Q := C[T];
  P := select(U -> nops(U) > 1,Q);
  if not `is_equal/E/Fbar`(N)(P,Q)(x[T],y[T]) then
   return false;
  fi;
 od;

 return true;
end:

######################################################################

`theta/EP/tree_Fbar_alt` := (N::posint) -> (A::set) -> (TT) -> proc(x)
 local TT1,C,P,Q,T,U,y;

 TT1 := `big_sets/trees`(TT);

 y := table();
 C := children_map(TT);
 for T in TT1 do
  Q := C[T];
  P := select(U -> nops(U) > 1,Q);
  y[T] := `theta/E/Fbar`(N)(P,Q)(x[T],y[T]);
 od;

 return eval(y);
end:

######################################################################

`phi/EP/tree_Fbar` := (N::posint) -> (A::set) -> (TT) -> proc(tpqx)
 local s,t,p,q,x,y,z,t0,x0,C,P,TT1,UU,VV,T,U,V,u,y0,n0,r;
 t := table();
 p := table();
 q := table();
 x := table();

 C := children_map(TT);
 P := parent_map(TT);
 TT1 := `big_sets/trees`(TT);

 for T in TT1 do
  t0[T],p[T],q[T],x0[T] := op(tpqx[T]);
  for U in C[T] do
   x[U] := x0[T][U];
   if nops(U) > 1 then t[U] := t0[T][U]; fi;
  od;
 od;

 s := table();

 for T in TT1 do
  UU := select(U -> (U minus T = {}),TT1);
  for U in UU do
   VV := select(V -> (U minus V = {}),UU);
   VV := sort([op(VV)],(V1,V2) -> evalb(nops(V1) < nops(V2)));
   r := nops(VV) - 1;
   s[T,U] := q[U] * mul(p[VV[i]],i=2..r+1) * mul(t[V[i]],i=1..r);
  od;
 od;

 y := table();
 for T in TT do
  y[T] := table();
  for u in A do y[T][u] := 0; od;
  for U in C[T] do
   for u in U do
    y[T][u] := x[T][U];
   od;
  od;
  y0 := `sum/vector_function`(N)(T)(y[T]) /~ nops(T);
  for u in T do y[T][u] := y[T][u] -~ y0; od;
  n0 := `norm/vector_functions`(y[T]);
  for u in T do y[T][u] := y[T][u] /~ n0; od;
 od;

 z := table();
 for T in TT1 do
  z[T] := table();
  for u in T do z[T][u] := [0$N]; od:

  UU := select(U -> (U minus T = {}),TT1);
  for U in UU do
   for u in U do
    z[T][u] := z[T][u] +~ s[T,U] *~ y[U][u];
   od;
  od;
 od;

 return eval(z);
end:


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)

params_a_n := 19:

params_a_a := [
   13/12, 7/6,  8/6,  9/6, 10/6, 17/12, 9/6, 10/6,
   11/6, 10/6, 11/6, 12/6, 10/6, 11/6, 12/6,  7/6,
    8/6,  9/6, 10/6.0
]:

params_a_b := [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1]:
params_a_c := [0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0]:
params_a_d := [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0]:

params_a_omega := [
   +0.678831e+00, -0.175821e+01, +0.127676e+01, -0.160789e+01, +0.365610e+00, -0.181327e+00,
   +0.146973e+00, +0.147141e+00, -0.716917e-01, -0.407167e-01, +0.214625e-01, -0.768156e-03,
   +0.310377e-01, -0.720326e-01, +0.446562e-01, -0.266802e+00, +0.150822e+01, -0.194515e+01,
   +0.679078e+00
]:

$include "th.mpl"

f := (rs, z, xt, xs0, xs1) -> f_th(rs, z, xt, xs0, xs1):


text:
impulse process MMH_loadChunks {
/summon area_effect_cloud ~ ~ ~ {CustomName:"MMH_Chunkloader",Duration:1}
/spreadplayers 23 23 0 1 false @e[type=area_effect_cloud,name=MMH_Chunkloader]
/kill @e[type=area_effect_cloud,name=MMH_Chunkloader]
}


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: work_mgga_x *)

params_a_mu := MU_GE:
params_a_b  := 0.4:
params_a_c  := 2.14951:
params_a_e  := 1.987:

mgamma      := 0.000023:
ff          := 3:

mkappa := (x, t) ->
  1/(mgamma/params_a_mu^2 + mgamma/params_a_mu + 1):

$include "tpss_x.mpl"

(* Equation (8) *)

f   := (rs, x, t, u) -> 1 + fx(x, t)*exp(-mgamma*fx(x, t)/params_a_mu)/(1 + fx(x, t)) \
    + (1 - exp(-mgamma*fx(x, t)^2/params_a_mu^2))*(params_a_mu/fx(x, t) - 1):


text:
function pascal(n)
    if n < 7 then
        field_width = 4
    else
        field_width = floor((lgamma(n+1) - 2*lgamma((n+1)/2))/log(10) + 1)
        if field_width mod 2 == 1 then
            field_width = field_width + 1
        end
    end
    line_width = (n+1)*field_width
    p = list(n+1)
    for i in 0:n
        p[i] = 1
	for j in i-1:1:-1
	    p[j] = p[j-1] + p[j]
	end
        s = ""
	for j in 0:i
	    s = s + center(ffmt(p[j], 0), field_width)
	end
	center(s, line_width)
    end
end

pascal(13)

max_stars_repo_name: yugyfoog/MPL

max_issues_repo_name: yugyfoog/MPL

max_forks_repo_name: yugyfoog/MPL


text:
 # a function returning (i + j) * -998
 func $foo (
# var %i xxx,
  var %i i32,
   var %j i32) i32 {
   return (
     mul i32 (
       add i32 ( dread i32 %i, dread i32 %j),
       constval i32 -998))}

 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

max_stars_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_issues_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_forks_repo_name: harmonyos-mirror/OpenArkCompiler-test


text:
# This is supposed to be a toy 2D model for issues that arise in
# computational geometry for computer aided design.
# We work everywhere in the real projective plane, and consider
# only straight lines and conic curves and segments thereof.
# This is supposed to make everything nicely computable.

dp3 := proc(x,y) local i; add(x[i]*y[i],i=1..3); end;
nm3 := proc(x) local i; sqrt(add(x[i]^2,i=1..3)); end;
xp3 := (u,v) -> [u[2]*v[3]-u[3]*v[2],
                 u[3]*v[1]-u[1]*v[3],
                 u[1]*v[2]-u[2]*v[1]];
tp3 := (u,v,w) -> Determinant(Matrix([u,v,w]));
eq3 := (u,v) -> evalb(xp3(u,v) = [0$3]);
sp3 := (u) -> [u[1]/u[3],u[2]/u[3]];

`plane_stereo/RP` := (x) -> [x[1]/x[3],x[2]/x[3]];
`plane_unstereo/RP` := (u) -> [u[1],u[2],1];

`disc_stereo/RP` := (x) -> [x[1],x[2]] *~ (signum(x[3])/nm3(x));

`disc_unstereo/RP` := (u) ->
  [u[1],u[2],sqrt(1-u[1]^2-u[2]^2)];

######################################################################

`is_element/RP_points` := proc(x)
 type(x,[numeric,numeric,numeric]);
end:

`is_leq/RP_points` := NULL:
`list_elements/RP_points` := NULL:
`count_elements/RP_points` := NULL:

`random_element/RP_points` := () -> `random_element/R`(3)();

`dist/RP_points` := proc(x,y)
 local i;
 1 - add(x[i]*y[i],i=1..3)^2/add(x[i]^2,i=1..3)
end:

`is_equal/RP_points` := (x,y) -> evalb(xp3(x,y) = [0$3]);

`in_general_position/RP_points` :=
 (u,v,w) -> evalb(Determinant(Matrix([u,v,w])) = 0);

`disc_plot/point` := (x) -> point(`disc_stereo/RP`(x),args[2..-1]);

######################################################################

`is_element/RP_lines` := proc(x)
 type(x,[numeric,numeric,numeric]);
end:

`is_leq/RP_lines` := NULL:
`list_elements/RP_lines` := NULL:
`count_elements/RP_lines` := NULL:

`random_element/RP_lines` := () -> `random_element/R`(3)();

`dist/RP_lines` := proc(x,y)
 local i;
 1 - add(x[i]*y[i],i=1..3)^2/add(x[i]^2,i=1..3)
end:

`is_equal/RP_lines` := (x,y) -> evalb(xp3(x,y) = [0$3]);

`in_general_position/RP_lines` :=
 (u,v,w) -> evalb(Determinant(Matrix([u,v,w])) = 0);

`track/RP_lines` := proc(x,t)
 local u,v;
 u,v := op(NullSpace(Matrix(x)));
 return (1-t) *~ u +~ t *~ v;
end:

`disc_plot/line` := proc(x)
 local u,v,w,p,t,t0,R;
 u,v := op(NullSpace(Matrix(x)));
 w := cos(Pi*t) *~ u +~ sin(Pi*t) *~ v;
 t0 := fsolve(w[3]);
 if evalf(subs(t = t0,diff(w[3],t))) > 0 then
  R := (t0 + 0.0001) .. (t0 + 0.9999);
 else
  R := (t0 - 0.9999) .. (t0 - 0.0001);
 fi;
 p := simplify(`disc_stereo/RP`(w));
 plot([op(p),t=R],args[2..-1]);
end:

######################################################################

`is_incident/RP` := proc(x,y) local i; evalb(add(x[i]*y[i],i=1..3) = 0); end;

`point_join/RP` := (x,y) -> xp3(x,y);

`line_meet/RP`  := (x,y) -> xp3(x,y);

######################################################################

`is_element/RP_conics` := proc(q)
 local c;
 if not(type(q,[numeric$6])) then
  return false;
 fi;

 # We now need to check whether the associated quadratic form
 # has one negative and two positive eigenvalues.  This works
 # out as follows.

 c := [-q[1]-q[2]-q[3],
       (q[1]*q[2]+q[2]*q[3]+q[3]*q[1]) - (q[4]^2+q[5]^2+q[6]^2)/4,
       (q[1]*q[4]^2+q[2]*q[5]^2+q[3]*q[6]^2-q[4]*q[5]*q[6])/4
         - q[1]*q[2]*q[3]
      ];

 if c[3] > 0 and (c[1] <= 0 or c[2] <= 0) then
  return true;
 fi;

 return false;
end:

`is_leq/RP_conics` := NULL:
`list_elements/RP_conics` := NULL:
`count_elements/RP_conics` := NULL:

`random_element/RP_conics` := proc()
 local q;

 q := `random_element/R`(6)();
 while not(`is_element/RP_conics`(q)) do
  q := `random_element/R`(6)();
 od:

 return q;
end:

`dist/RP_conics` := proc(q,r)
 local i,j;
 sqrt(add(add((q[i]*r[j]-q[j]*r[i])^2,j=i+1..6),i=1..6)/
      (add(q[i]^2,i=1..6) * add(r[i]^2,i=1..6)));
end:

`is_equal/RP_conics` := proc(q,r)
 evalb(`dist/RP_conics`(q,r) = 0);
end:

`conic_matrix/RP` := (q) ->
  <<q[1]|q[6]/2|q[5]/2>,<q[6]/2|q[2]|q[4]/2>,<q[5]/2|q[4]/2|q[3]>>;

`conic_unmatrix/RP` := (Q) ->
 [Q[1,1],Q[2,2],Q[3,3],2*Q[2,3],2*Q[3,1],2*Q[1,2]];

`conic_eval/RP` := (q) -> (x) ->
 q[1]*x[1]^2 + q[2]*x[2]^2 + q[3]*x[3]^2 +
 q[4]*x[2]*x[3] + q[5]*x[3]*x[1] + q[6]*x[1]*x[2];

`conic_bilin/RP` := (q) -> (x,y) ->
 (`conic_eval/RP`(q)(x +~ y) - `conic_eval/RP`(q)(x -~ y))/4;

`conic_coeffs/RP` := (qx,x) -> [
 coeff(qx,x[1],2),
 coeff(qx,x[2],2),
 coeff(qx,x[3],2),
 coeff(coeff(qx,x[2],1),x[3],1),
 coeff(coeff(qx,x[3],1),x[1],1),
 coeff(coeff(qx,x[1],1),x[2],1)
];

`disc_implicitplot/conic` := proc(q)
 local u,v,qq;

 qq := numer(simplify(`conic_eval/RP`(q)(`disc_unstereo/RP`([u,v]))));
 implicitplot(qq,u=-1.1..1.1,v=-1.1..1.1,args[2..-1]);
end:

######################################################################

`is_element/RP_arcs` := proc(xxc)
 if not(type(xxc,[[numeric$3]$3])) then
  return false;
 fi;

 if Determinant(Matrix(xxc)) = 0 then return false; fi;

 return true;
end:

`random_element/RP_arcs` := proc()
 local xxc,i,j;
 xxc := [[0$3]$3];

 while Determinant(Matrix(xxc)) = 0 do
  xxc := [seq([seq(rand(-10..10)(),i=1..3)],j=1..3)];
 od;

 return xxc;
end:

`ratio/RP_arcs` := proc(xxc1,xxc2)
 local u,v,w,i;

 u := [seq(dp3(xxc1[i],xxc2[i]),i=1..3)];
 v := [seq(dp3(xxc1[i],xxc1[i]),i=1..3)];
 w := [seq(dp3(xxc2[i],xxc2[i]),i=1..3)];

 if (u *~ u = v *~ w) then
  return u /~ w;
 else
  return NULL;
 fi;
end:

`is_equal/RP_arcs` := proc(xxc1,xxc2)
 local r;
 r := `ratio/RP_arcs`(xxc1,xxc2);

 if r = NULL then return false; fi;
 if r[1] <> r[2] or r[2] <> r[3] then
  return false;
 fi;

 return true;
end:

`is_similar/RP_arcs` := proc(xxc1,xxc2)
 local r;
 r := `ratio/RP_arcs`(xxc1,xxc2);

 if r = NULL then return false; fi;
 if r[1]*r[2] - r[3]^2 <> 0 then
  return false;
 fi;

 return true;
end:

`is_leq/RP_arcs` := NULL:
`list_elements/RP_arcs` := NULL:
`count_elements/RP_arcs` := NULL:

`arc_eval/RP` := (xxc) -> (t) ->
 (1-t)^2 *~ xxc[1] +~ t^2 *~ xxc[2] +~ (t*(1-t)) *~ xxc[3];

`conic_arc_eval/RP` := (q,xxc) -> (t) ->
 `conic_eval/RP`(q)(`arc_eval/RP`(xxc)(t));

`conic_contains_point/RP` := proc(q,x)
 local err;
 err := abs(`conic_eval/RP`(q)(x));
 return evalb(err = 0);
end:

`conic_contains_arc/RP` := proc(q,xxc)
 local t,err;
 err := max(map(abs,[coeffs(expand(`conic_arc_eval/RP`(q,xxc)(t)),t)]));
 return evalb(err = 0);
end:

`normalise_similar/RP_arcs` := proc(xxc)
 local x,c,n;
 x[1] := xxc[1];
 x[2] := xxc[2];
 c := xxc[3];
 n[1] := nm3(x[1]);
 n[2] := nm3(x[2]);
 return [x[1]/~n[1],x[2]/~n[2],c/~sqrt(n[1]*n[2])];
end:

`arc_coeffs/RP` := proc(st,t)
 local su,u;
 su := factor((1+u)^2 *~ subs(t=u/(1+u),st));
 [map(coeff,su,u,0),
  map(coeff,su,u,2),
  map(coeff,su,u,1)];
end:

`arc_from_tangents/RP` := proc(x1,x2,u1,u2)
 local m1,m2,n,c;
 m1 := Determinant(Matrix([x1,x2,u2]));
 m2 := Determinant(Matrix([x1,x2,u1]));
 n := - Determinant(Matrix([x1,u1,u2])) * m2/m1;
 c := n *~ x2 +~ m2 *~ u2;
 return [m1 *~ x1,m2 *~ x2,n *~ c];
end:

`arc_inv/RP` := (xxc) -> proc(y)
 local u1,u2,u3;
 u1 := Determinant(Matrix(x1,x2,y));
 u2 := Determinant(Matrix(x2,c ,y));
 u3 := Determinant(Matrix(c ,x1,y));
 return u2/(u1+u2); # or u1/(u1+u3);
end:

`conic_from_arc/RP` := proc(xxc)
 local x1,x2,c,M;
 x1,x2,c := op(xxc);

 M := Matrix(
  [[x1[1]^2, x1[2]^2, x1[3]^2,
     x1[2]*x1[3], x1[1]*x1[3], x1[1]*x1[2]],
   [x2[1]^2, x2[2]^2, x2[3]^2,
     x2[2]*x2[3], x2[1]*x2[3], x2[1]*x2[2]],
   [2*c[1]*x1[1], 2*c[2]*x1[2], 2*c[3]*x1[3], c[2]*x1[3]+c[3]*x1[2],
     c[1]*x1[3]+c[3]*x1[1], c[1]*x1[2]+c[2]*x1[1]],
   [2*c[1]*x2[1], 2*c[2]*x2[2], 2*c[3]*x2[3], c[2]*x2[3]+c[3]*x2[2],
     c[1]*x2[3]+c[3]*x2[1], c[1]*x2[2]+c[2]*x2[1]],
   [c[1]^2+2*x1[1]*x2[1], c[2]^2+2*x1[2]*x2[2], c[3]^2+2*x1[3]*x2[3],
     c[2]*c[3]+x1[2]*x2[3]+x1[3]*x2[2],
      c[1]*c[3]+x1[1]*x2[3]+x1[3]*x2[1],
       c[1]*c[2]+x1[1]*x2[2]+x1[2]*x2[1]]]
 );
 return NullSpace(M)[1];
end:

`arc_from_rational_conic/RP` := proc(q)
 local x,y,z,f,sol,u0,u1,u2,u,t,Z;
 if not(type(q,[rational$6])) then
  return NULL;
 fi;

 f := `conic_eval/RP`(q)([x,y,z]);
 sol := isolve(f);
 if sol = NULL then return NULL; fi;

 u0 := subs(sol,[x,y,z]);
 u1 := eval(subs(igcd = (() -> 1),u0));
 u2 := u1;
 for Z in indets(u1) do
  if degree(u1[1],Z) = 1 then u2 := subs(Z = 1,u2); fi;
 od:
 u := subs({indets(u2)[1] = t,indets(u2)[2] = 1-t},u2);
 return `arc_coeffs/RP`(u,t);
end:

`arc_from_conic_eigenvalues/RP` := proc(q)
 local Q,E,V,i,v,xxc;
 Q := evalf(Matrix(`conic_matrix/RP`(q),shape=symmetric));
 E,V := Eigenvectors(Q);
 for i from 1 to 3 do
  v[i] := convert(Column(V,i)/sqrt(abs(E[i])),list);
 od;
 xxc := [(v[1] +~ v[2])/~2,(v[1] -~ v[2])/~2,v[3]];
 return xxc;
end:

`cut_arc/RP` := proc(xxc,t1,t2)
 local y1,y2,d;
 y1 := `arc_eval/RP`(xxc)(t1);
 y2 := `arc_eval/RP`(xxc)(t2);
 d := (2*(1-t1)*(1-t2)) *~ xxc[1] +~
      (2*t1*t2) *~ xxc[2] +~
      (t1 + t2 - 2*t1*t2) *~ xxc[3];
 return [y1,y2,d];
end:

`shift_arc/RP` := proc(xxc,u)
 [exp(u),exp(-u),1] *~ xxc;
end:

`bulge_arc/RP` := proc(xxc,u)
 local x1,x2,c;
 x1,x2,c := op(xxc);
 if u >= 0 then
  return [x1,x2,u*~c];
 else
  return [x1,x2,u*~c];
 fi;
end:

`disc_plot/arc` := proc(xxc)
 local u,v,tt,m,i,P,e,t0,t1;
 u := `arc_eval/RP`(xxc)(t);
 tt := sort(evalf([op({0,1,fsolve(u[3],t=0..1)})]));
 m := nops(tt);
 P := NULL;
 e := 10.^(-4);
 for i from 1 to m-1 do
  t0 := (1-e) *~ tt[i] +~ e *~ tt[i+1];
  t1 := e *~ tt[i] +~ (1-e) *~ tt[i+1];
  v := simplify(`disc_stereo/RP`(u)) assuming t > t0 and t < t1;
  P := P,plot([op(v),t=t0..t1],args[2..-1]);
 od:
 display(P);
end:

`disc_point_plot/arc` := proc(xxc)
 local i,N,P,u,v,p;
 N := 12;
 P := NULL;
 for i from 0 to N do
  u := evalf(`arc_eval/RP`(xxc)(i/N));
  v := `disc_stereo/RP`(u);
  p := point(v,args[2..-1]);
  P := P,p;
 od:
 display(P);
end: