text:
# Class: RealAlgebraicNumber
#
# Description:
#  Implementation of real algebraic numbers together with their comparison.
#  This implementation is inspired by the implementation in CGAL 4.7.
#
# Author:
#  Kacper Pluta - kacper.pluta@esiee.fr
#  Laboratoire d'Informatique Gaspard-Monge - LIGM, A3SI, France
#
# Date:
#  11/12/2015
#
# License:
#  Simplified BSD License
#
# Copyright (c) 2015, 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 RealAlgebraicNumber()
  option object;
  (* Univariete polynomaial *)
  local poly::polynom;
  (* Lower bound of the range in which exists only one real root of poly. *)
  local a::rational;
  (* Upper bound of the range in which exists only one real root of poly. *)
  local b::rational;
  (* Real algebraic number is rational when a = b and sign of poly at a/b is 0. *)
  local isRational_;
  (* Note that isolating interval has to be open iff a real algebraic number is not rational and
  closed, a = b, otherwise.*)

# Method: ModuleCopy
#   Standard constructor / copy constructor
#
# Parameters:
#   self::RealAlgebraicNumber     - a new object to be constructed
#   proto::RealAlgebraicNumber    - a prototype object from which self is derived
#   poly::polynom                 - a univariate polynomial
#   a::rational                   - a lower bound of the range in which exists only one real root of poly
#   b::rational                   - an upper bound of the range in which exists only one real root of poly
#
# Output:
#   An object of type RealAlgebraicNumber.
#
# Exceptions:
#  "Invalid range. A range is valid when: a <= b."
#  "Degree of %1 is invalid."
#
  export ModuleCopy::static := proc( self::RealAlgebraicNumber,
                                     proto::RealAlgebraicNumber,
                                     poly::polynom,
                                     a::rational,
                                     b::rational, $ )
    local signAtA, signAtB;
    if _passed = 2 then
      self:-poly := proto:-poly;
      self:-a := proto:-a;
      self:-b := proto:-b;
      self:-isRational_ := proto:-isRational_;
    else
      if upperbound( indets( poly ) ) > 1 then
        error "%1 is not univariate!", poly;
      end if;
      if a > b then
        error "Invalid range. A range is valid when: a <= b.";
      end if;
      if gcd( poly, diff( poly, op( indets( poly ) ) ) ) <> 1 then
        error "Polynomial: %1 is not square-free.", poly;
      end if;
      self:-poly := poly;
      if degree(poly) >= 1 then
        signAtA := signum( eval( poly, indets( poly )[1] = a ) );
        signAtB := signum( eval( poly, indets( poly )[1] = b ) );
        self:-a := a;
        self:-b := b;
        self:-isRational_ := evalb( self:-a = self:-b and signAtA = 0 );
        if signAtA = 0 and signAtB <> 0 then
          WARNING("Incorrect interval. Sign of univariate polynomial on one side of the interval"
          " evaluated to zero but not on the another. Interval fixed.");
          self:-b := self:-a;
          self:-isRational_ := true:
        elif signAtB = 0 and signAtA <> 0 then
          WARNING("Incorrect interval. Sign of univariate polynomial on one side of the interval"
          " evaluated to zero but not on the another. Interval fixed.");
          self:-a := self:-b;
          self:-isRational_ := true:
        elif signAtA = signAtB and self:-a <> self:-b then
          error "Interval incorrect! No root in the interval: (%1, %2), for %3 .", self:-a, self:-b,
          self:-poly;
        fi:
      elif degree( poly ) = 0 then
        self:- denom( poly ) * 'a'  - numer( poly );
        self:-a := poly;
        self:-b := poly;
        self:-isRational_ := true;
      else
        error "Degree of %1 is invalid.", poly;
      end if;
    end if;
    return self;
  end proc:

# Method: ModulePrint
#   Standard printout of an object of type RealAlgebraicNumber.
#
# Parameters:
#   self::RealAlgebraicNumber      - a real algebraic number
#
  export ModulePrint::static := proc( self::RealAlgebraicNumber )
    if(self:-a = self:-b) then
        nprintf( "( %a, [%a, %a] )", self:-poly, self:-a, self:-b );
    else
        nprintf( "( %a, ]%a, %a[ )", self:-poly, self:-a, self:-b );
    end if;
  end proc;

# Method: ModuleApply
#   Define standard constructor.
#
  export ModuleApply::static := proc()
   Object(RealAlgebraicNumber, args)
  end proc;

# Method: ModuleDeconstruct
#   Provides information how to recreate an object after being serialized.
#
# Parameters:
#   self::RealAlgebraicNumber        - a real algebraic number
#
  export ModuleDeconstruct := proc( self::RealAlgebraicNumber )
    ('RealAlgebraicNumber')(self:-poly, self:-a, self:-b)
  end proc;

# Method: GetPolynomial
#   A getter method to access the univariate polynomial of RealAlgebraicNumber.
#
# Parameters:
#   self::RealAlgebraicNumber      - a real algebraic number
#
# Output:
#   Univariate polynomial stored in self.
#
  export GetPolynomial::static := proc( self::RealAlgebraicNumber )
    return self:-poly;
  end proc:

# Method: GetInterval
#   A getter method to access the range isolating a root of univariate polynomial.
#
# Parameters:
#   self::RealAlgebraicNumber      - a real algebraic number
#
# Output:
#   The range isolating a root of univariate polynomial -- self:-poly.
#   Upper and lower bounds of a range are rationals. When lower = upper
#   then a real algebraic number is rational.
#
  export GetInterval::static := proc( self::RealAlgebraicNumber )
    return [ self:-a, self:-b ];
  end proc:

# Method: IsRational
#   A method to check if a real algebraic number is rational.
#
# Parameters:
#   self::RealAlgebraicNumber      - a real algebraic number
#
# Output:
#   True when a real algebraic number is rational, false
#   otherwise. A real algebraic number is meant as rational when a = b and
#   when sign of poly at a/b is zero.
#
  export IsRational::static := proc( self::RealAlgebraicNumber )
    return self:-isRational_;
  end proc:

# Method: CompareRational
#   A method used to compare a real algebraic number with a
#   rational number.
#
# Parameters:
#   self::RealAlgebraicNumber      - a real algebraic number
#   m::rational                    - a rational number
#
# Output:
#   -1 when a real algebraic number is smaller than a rational
#    number, 0 when they are equal and 1 when a real algebraic
#    number is bigger than a rational.
#
  local CompareRational::static := proc( self::RealAlgebraicNumber, m::rational )
    local refined:
    refined := StrongRefineAt(self,m);
    if evalb( refined:-a < m ) then
      return -1;
    elif evalb( refined:-a > m ) then
      return 1;
    elif evalb( signum( eval( refined:-poly, op( indets( refined:-poly ) ) = m ) ) = 0 ) then
      return 0;
    end if;
  end proc:

# Method: RefineAt
#   A method used to refine a real algebraic number using a rational
#   number for adaptation of a range isolating a root of poly.
#
# Parameters:
#   self::RealAlgebraicNumber      - a real algebraic number
#   m::rational                    - a rational number
#
# Output:
#   A RealAlgebraicNumber obtaind from self refined at m.
#
# Comment:
#   Not that the type can change to rational.
#
  local RefineAt::static := proc( self::RealAlgebraicNumber, m::rational )
    local signAtM, f::polynom, g::polynom;
    local var := op( indets( self:-poly ) );
    if self:-isRational_ or m <= self:-a or self:-b <= m then
      return self;
    end if;
    signAtM := signum( eval( self:-poly, var = m ) );
    if evalb( signAtM = 0 ) then
     g := denom( m ) * var  - numer( m );
     return Object( RealAlgebraicNumber, g, m, m );
    elif evalb( signum( eval( self:-poly, var = self:-a ) ) = signAtM ) then
      return Object( RealAlgebraicNumber, self:-poly, m, self:-b );
    elif evalb( signum( eval( self:-poly, var = self:-b ) ) = signAtM ) then
      return Object( RealAlgebraicNumber, self:-poly, self:-a, m );
    else
      return self;
    end if;
  end proc:

# Method: BisectRange
#   A method used to compare a real algebraic number with a
#   rational number.
#
# Parameters:
#   self::RealAlgebraicNumber      - a real algebraic number
#
# Output:
#   Refine an isolating range at ( self:-a + self:-b ) / 2
#
  local BisectRange::static := proc( self::RealAlgebraicNumber )
    return RefineAt( self, ( self:-a + self:-b ) / 2 );
  end proc:

# Method: StrongRefineAt
#   A method used to refine a real algebraic number using a rational
#   number for adaptation of a range isolating a root of poly.
#
# Parameters:
#   self::RealAlgebraicNumber      - a real algebraic number
#   m::rational                    - a rational number
#
  local StrongRefineAt::static := proc( self::RealAlgebraicNumber, m::rational )
    local refined:
    if self:-isRational_ or signum( eval( self:-poly, indets( self:-poly )[1] = m ) ) = 0 then
      return self;
    fi:
    refined := self;
    while refined:-a <= m and m <= refined:-b do
      refined := BisectRange(refined);
    od:
    return refined:
  end proc:

  export DisjointRanges::static := proc(a::RealAlgebraicNumber, b::RealAlgebraicNumber)
    local ll := a, rr := b, i;

    (* No intersection.*)
    if evalb( ll:-b < rr:-a ) or evalb( ll:-a > rr:-b ) then
      return [a,b];
    fi:

    for i from 1 while 1 = 1 do
      ll := BisectRange( ll ):
      rr := BisectRange( rr ):

      (* No intersection.*)
      if evalb( ll:-b < rr:-a ) or evalb( ll:-a > rr:-b ) then
        return [ll,rr];
      fi:
    od:
  end proc:

# Method: Compare
#   A method used to compare two real algebraic numbers.
#
# Parameters:
#   l::RealAlgebraicNumber      - a real algebraic number
#   r::RealAlgebraicNumber      - a real algebraic number
#
# Output:
#   -1 when l is smaller than r, 0 when they are equal and 1 when l is bigger than r.
#
  export Compare::static := proc( l::RealAlgebraicNumber, r::RealAlgebraicNumber, $ )
    local i::integer, a::rational, b::rational, F1::polynom, F2::polynom, G::polynom;
    local ll::RealAlgebraicNumber, rr::RealAlgebraicNumber;

    if indets(l:-poly) <> indets(r:-poly) then
      error "Univariate polynomials have different variables: %1 and %2.", indets(l:-poly),
      indets(r:-poly);
    fi;

    if evalb( l:-poly = r:-poly and l:-a = r:-a and l:-b = r:-b ) then
        return 0;
    end if;

    (* When rationals *)
    if r:-isRational_ then
      return CompareRational( l, r:-a );
    elif l:-isRational_ then
      return -CompareRational( r, l:-a );
    end if;

    (* Check if there is no intersection of the ranges *)
    if evalb( l:-b < r:-a ) then
      return -1;
    elif evalb( l:-a > r:-b ) then
      return 1;
    end if:

    (* Get the intersecting interval *)
    if evalb( l:-a > r:-a ) then
      a := l:-a;
    else
      a := r:-a;
    end if;
    if evalb( l:-b < r:-b ) then
      b := l:-b;
    else
      b := r:-b;
    end if;

    (* refine at the intersecting interval *)
    ll := RefineAt( l, a ):
    ll := RefineAt( ll, b ):
    rr := RefineAt( r, a ):
    rr := RefineAt( rr, b ):

    (* Refiment can change type to rational. *)
    if rr:-isRational_ then
      return CompareRational( ll, rr:-a );
    elif ll:-isRational_ then
      return -CompareRational( rr, ll:-a );
    end if;

    (* Check if there is no intersection after refiment. *)
    if evalb( ll:-b < rr:-a ) then
      return -1;
    elif evalb( ll:-a > rr:-b ) then
      return 1;
    end if;

    (* The number of roots of the GCD of two polynomials is equal to the number of common roots.
       use this to simplify the problem in the intersecting range.*)
    G := gcd( ll:-poly, rr:-poly );
    F1 := simplify( ll:-poly / G );
    F2 := simplify( rr:-poly / G );

    if evalb( signum( eval( G, op( indets( G ) ) = ll:-a ) ) <> signum( eval( G,
      op( indets( G ) ) = ll:-b ) ) ) then
      ll := Object( ll, G, ll:-a, ll:-b ):
    else
      ll := Object( ll, F1, ll:-a, ll:-b ):
    end if:

    if evalb( signum( eval( G, op( indets( G ) ) = rr:-a ) ) <> signum( eval( G,
      op( indets( G ) ) = rr:-b ) ) ) then
      rr := Object( rr, G, rr:-a, rr:-b ):
    else
      rr := Object( rr, F2, rr:-a, rr:-b ):
    end if:

    (* Use of GCD can change type to rational. *)
    if rr:-isRational_ then
      return CompareRational( ll, rr:-a );
    elif ll:-isRational_ then
      return -CompareRational( rr, ll:-a );
    end if;

    (* Check for equality. *)
    if evalb( signum( eval( G, op( indets( G ) ) = a ) ) <> signum( eval( G,
      op( indets( G ) ) = b ) ) ) then
      return 0;
    end if;

    (* Refiment until disjoitness. *)
    for i from 1 while 1 = 1 do
      ll := BisectRange( ll ):
      rr := BisectRange( rr ):

      (* Rationals after refiment. *)
      if rr:-isRational_ then
        return CompareRational( ll, rr:-a );
      elif ll:-isRational_ then
        return -CompareRational( rr, ll:-a );
      end if;

      (* No intersection after refiment. *)
      if evalb( ll:-b < rr:-a ) then
        return -1;
      elif evalb( ll:-a > rr:-b ) then
        return 1;
      end if:
    end do:
  end proc:

# Method: < operator
#   A method used to compare two real algebraic numbers.
#
# Parameters:
#   l::RealAlgebraicNumber      - a real algebraic number
#   r::RealAlgebraicNumber      - a real algebraic number
#
# Output:
#   true when l is smaller than r and false otherwise.
#
  export `<`::static := proc( l, r, $ )
   if ( _npassed <> 2 or not l::RealAlgebraicNumber or not r::RealAlgebraicNumber ) then
    return false;
  end if;
  if Compare( l, r ) = -1 then
    return true;
  else
    return false;
  end if;
  end proc:

# Method: <= operator
#   A method used to compare two real algebraic numbers.
#
# Parameters:
#   l::RealAlgebraicNumber      - a real algebraic number
#   r::RealAlgebraicNumber      - a real algebraic number
#
# Output:
#   true when l is smaller or equal to r and false otherwise.
#
  export `<=`::static := proc( l, r, $ )
   if ( _npassed <> 2 or not l::RealAlgebraicNumber or not r::RealAlgebraicNumber ) then
    return false;
  end if;
  if Compare( l, r ) <= 0 then
    return true;
  else
    return false;
  end if;
  end proc:

# Method: = operator
#   A method used to compare two real algebraic numbers.
#
# Parameters:
#   l::RealAlgebraicNumber      - a real algebraic number
#   r::RealAlgebraicNumber      - a real algebraic number
#
# Output:
#   true when l is equal to r and false otherwise.
#
  export `=`::static := proc( l, r, $ )
    if Compare( l, r ) = 0 then
      return true;
    else
      return false;
    end if;
  end proc:

end module:

max_stars_repo_name: copyme/MapleTools

max_issues_repo_name: copyme/MapleTools

max_forks_repo_name: copyme/MapleTools


text:
func $foo ( var %i i32 ) i32
func $fact ( var %n i32) i32
func $fact ( var %n i32) i32
func $foo ( var %i i32 ) i32

func $fact (
 var %n i32) i32 {
 if (ne i32 i32 (dread i32 %n, constval i32 1)) {
    call &foo(
      sub i32 (dread i32 %n, constval i32 1))
    return (regread i32 %%retval)}
   else {
     return (constval i32 1) } }

func $foo ( var %i i32 ) i32 {
  return (dread i32 %i) }
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

ext: 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:
(*
 Copyright (C) 2020 Susi Lehtola

 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_pgslb_params *params;

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

# Equation (4) and (8)
pgslb_f0 := (s, q) -> 5/3*s^2 + exp(-params_a_pgslb_mu * s^2) + params_a_pgslb_beta*q^2:
pgslb_f := (x, u) -> pgslb_f0(X2S*x, X2S^2*u):

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


text:
restart:
Digits := 100:
with (numapprox):with(orthopoly):
read "common-procedures.mpl";
mkdir("TEMPATAN");

#First, quick phase
#Approach for the first function:
#atan(x) = atan(b(i)) + atan ( (x-bi) / (1+x*b(i) ) )
#We try to have 64 correct bits

#Constants

#if (x >= xmax) : return arctan(x)/2)=neareast(Pi/2)

log[2](evalf(min(solve(arctan(x)=nearest(Pi/2),x)))):
xmax:=2^54:

#halfPi=Pi/2 (remark that maxvalue < Pi/2 so that atan(x) < Pi/2

halfPi:=nearest(Pi/2):
halfPi_to_plus_infinity:=nearest(nearest(Pi/2)*(1+2^(-53))):

# invPi, for the atanpi function

invpih,invpil := hi_lo(1/Pi);

# if (xmin<xmin) : return arctan(x)=x

xmin=2^(-27);

ieeehexa(halfPi); ieeehexa( halfPi_to_plus_infinity ); ieeehexa(2^(-27));ieeehexa(2^(54));ieeehexa(2^(-10));

#Interval parameters
#Some constants :
e:=2^(-6.3):
marge := 2^(-30):

#The intervals are (a[i],a[i+1]). the "middle" is b(i) :
maxx := bi -> simplify( solve( (x-bi) / (1+x*bi) = e ,x) ):
minx := bi -> simplify(solve( (x-bi) /(1+x*bi)=-e ,x)):

nextbi := proc (x) evalf( max(solve( minx(bi) = x ,bi) )*(1-marge)); end:

allbi := proc (n)
  local xi,nbi,x,i,j;
  global b,a, nb_of_ai, nb_of_bi, value_for_dicho;
  x := e;
  nbi := 0;
  i := 0;
  while(i<n and nbi < 1/e and nbi >= 0) do
     nbi := nearest ( nextbi(x) );
     b[i] := evalf( nbi );
     a[i] := x;
     x   := evalf(maxx(nbi));
     i:=i+1;
  od;
  j:=0;
  while ( 2^j < i ) do j:=j+1 od:
  nb_of_ai := i;
  nb_of_bi := i;
  b[i-1] := nearest(1/e+4):
  value_for_dicho := 2^j;
  return i,b[0],b[i-1];
end:

allbi(100);

#The polynome and its size ( 4 because 4 terms : x^2/3, ...)

Q9:= poly_exact(x-1/3*x^3+1/5*x^5 -1/7*x^7+1/9*x^9):
deg_poly := 4:
coef_poly[0]:=nearest (-1/3):
coef_poly[1]:=nearest (1/5):
coef_poly[2]:=nearest (-1/7):
coef_poly[3]:=nearest (1/9):

Q := x^2 * coef_poly[0] + x^4*coef_poly[1] + x^6*coef_poly[2] + x^8 * coef_poly[3]:

Qprime := poly_exact( x * coef_poly[0] + x^2*coef_poly[1] + x^3*coef_poly[2] + x^4 * coef_poly[3]):

log[2] (infnorm( (arctan(x)-x*(1+Q) )/x,x=2^(-53)..e));
# we can note that arctan(x) = sum((-1)^i*x^(2*i+1)/(2*i+1), i = (0 .. infinity))
# so ( arctan(x) - Q9(x) ) / x = sum((-1)^i*x^(2*i)/(2*i+1), i = (5 .. infinity))
# so abs ( (arctan(x) -Q9(x) )/x )  <= x^10 * sum((-1)^(i+1)*x^(2*i)/(2*i+1+10), i = (0 .. infinity))
#                                   <= x^10 * 1/11
#                                   <= 2^(-66) if x<2^(-6.3)

#Conputation of the RN constant

#Epsilon = relative error
#Delta = absolute error
#1st : Error about reduction:

EpsilonXminusBi := 2^(-105):

# 1 + x*b[i] :
Epsilon1A22xbi :=  2^(-104):

# Xred = (x-b[i]) / ( 1 + x*b[i] )

EpsilonXred := EpsilonXminusBi+Epsilon1A22xbi + 2^(-105):

DeltaXred := EpsilonXred * e:

log2(EpsilonXred);

#Polynomial evaluation :
#atan ~= x - x^3/3 + x^5/5 - x^7/7 + x^9/9
#        ~= x * ( 1 + Q(x^2) )
#    = (xhi+xlo) * ( 1 + Q)
#        ~= xhi + (xhi*Q + xlo)

#  Warning : the approx about x !!

#calc of Q:

EpsilonXred2 := (1+2^(-53) + 2^(-105) )^2 * (1+2^(-53)) -1 :

log2(EpsilonXred2);

errlist:=errlist_quickphase_horner( degree( x^4 ),0,0,EpsilonXred2, 2**(-53)):

rounding_error1:= compute_horner_rounding_error(Qprime,x,e,errlist,true):

EpsilonQ := rounding_error1[1]:
log2(EpsilonQ);

# Since q will be multiplied in x(1+q), its absolute error term is :

Qmax:=infnorm(Qprime,x=-e^2..e^2):

DeltaQ := e*Qmax*EpsilonQ;
log2(DeltaQ);

# calc of xi*Q + xlo + atan(b[i])lo

#First the truncation of xredlo*q : at most half an ulp of xredhi*q
DeltaTruncXredloQ :=  2^(-53) * e * Qmax:
log2(DeltaTruncXredloQ);
#then the sum of atanbilo and xredlo : the largest term is atanbilo
DeltaAdd1  := 0.5*ulp(2^(-53)*(Pi/2 + e)):
log2(DeltaAdd1);

#then compute xredhi*q
DeltaXredhiQ := 0.5*ulp(e*Qmax):
log2(DeltaXredhiQ);

# and the second addition
DeltaAdd2 :=  0.5*ulp(e*Qmax +  2^(-53)*(Pi/2 + e) ):
log2(DeltaAdd2);

DeltaAtanlolo := DeltaTruncXredloQ + DeltaAdd1 + DeltaXredhiQ + DeltaAdd2:
log2(DeltaAtanlolo);

# In the second reconstruction there is one more FP add which is
# aligned to the previous:
DeltaReconst := DeltaAtanlolo + DeltaAdd2 :
log2(DeltaReconst);

# error due to the approx of atan by P : < e^11/11
Deltaapprox := infnorm( (arctan(x)-x*(1+Q)), x=0..e):

log[2](Deltaapprox);

DeltaTotal := Deltaapprox + DeltaQ + DeltaReconst :
log2(DeltaTotal);
EpsilonTotal := evalf( DeltaTotal / arctan(xx) ):

EpsilonFinal := infnorm( EpsilonTotal,xx=2^(-6.3)..1,'xmaxmax'):
EpsilonFinal_i_10 := infnorm( EpsilonTotal,xx=a[10]..2^60,'xmaxmax'):
log2( EpsilonFinal );
log2( EpsilonFinal_i_10 );

# Computation of e and e_i_30 (which is e when i>30 (ie x>1)
E := evalf(compute_rn_constant(EpsilonFinal));
E_i_10 := evalf(compute_rn_constant(EpsilonFinal_i_10));

#An other Constant : when there is no reduction.
#We compute the polynom and then atan(x) = x + x.Q

Epsilonx2 := 2^(-53):
errlist:=errlist_quickphase_horner( degree( Qprime ),0,0,Epsilonx2, 2**(-53)):
rounding_error:= compute_horner_rounding_error(Qprime,x,e,errlist,true):

epsilon_Q := rounding_error[1]:
log[2](epsilon_Q);

epsilon_x_Q := 2^(-52):

delta_final := 2^(-105)*x + epsilon_Q*x^3 + epsilon_x_Q*x^3 - arctan(x) + x*(1+Q):
epsilon_final := infnorm(delta_final/x,x=2^(-10)..e):
epsilon_final_m10 := infnorm(delta_final/x,x=2^(-24)..2^(-10)):

log[2](epsilon_final);
E_no_reduction := evalf(compute_rn_constant(epsilon_final));
E_no_reduction_m10 := evalf(compute_rn_constant(epsilon_final_m10));

# Output

#test
filename:="TEMPATAN/atan_fast.h":

fd:=fopen(filename, WRITE, TEXT):
fprintf(fd, "\n/*File generated by maple/atan.mpl */\n"):

fprintf(fd, "#ifndef _ATAN_FAST_H\n#define _ATAN_FAST_H\n\n");

fprintf(fd, "#include \"crlibm.h\"\n#include \"crlibm_private.h\"\n"):

fprintf(fd,"#ifdef WORDS_BIGENDIAN\n"):
fprintf(fd, "static const db_number HALFPI = {{0x3FF921FB,0x54442D18}};"):
fprintf(fd, "static const db_number HALFPI_TO_PLUS_INFINITY = {{0x3FF921FB,0x54442D19}};"):
fprintf(fd,"\n#else\n"):
fprintf(fd, "static const db_number HALFPI = {{0x54442D18,0x3FF921FB}};"):
fprintf(fd, "static const db_number HALFPI_TO_PLUS_INFINITY = {{0x54442D19,0x3FF921FB}};"):
fprintf(fd,"\n"):
fprintf(fd,"#endif\n"):

fprintf(fd, "#define MIN_REDUCTION_NEEDED %1.50f\n",e):

fprintf(fd, "#define INVPIH %1.50f\n", invpih):
fprintf(fd, "#define INVPIL %1.50f\n", invpil):

fprintf(fd,"#define nb_of_ai %d\n",nb_of_ai):
fprintf(fd,"#define nb_of_bi %d\n",nb_of_bi):

# we multiply by 2^(-20) in order to prevent from rounding
fprintf(fd, "static const double rncst[4] ={\n"):
fprintf(fd, "%1.50f , /* i<10 */ \n",E*(1+2^(-20))) :
fprintf(fd, "%1.50f , /* i>10 */ \n",E_i_10*(1+2^(-20))) :
fprintf(fd, "%1.50f , /* e > 2^-10 */ \n",E_no_reduction*(1+2^(-20))) :
fprintf(fd, "%1.50f , /* e < 2^-10 */ \n };\n",E_no_reduction_m10*(1+2^(-20))) :

fprintf(fd, "static const double epsilon[4] ={\n"):
fprintf(fd, "%1.50e ,\n",EpsilonFinal *(1+2^(-20))) :
fprintf(fd, "%1.50e ,\n",EpsilonFinal_i_10*(1+2^(-20))) :
fprintf(fd, "%1.50e ,\n",epsilon_final*(1+2^(-20))) :
fprintf(fd, "%1.50e ,\n };\n",epsilon_final_m10*(1+2^(-20))) :

fprintf(fd, "#define DEGREE %d\n", deg_poly):
fprintf(fd, "static double const coef_poly[%d] = \n{\n",deg_poly):

for i from deg_poly-1 to 0 by -1 do
    fprintf(fd, "/* coef for degree %d */ %1.50f, \n" , 2*i+3, coef_poly[i]):
od:
fprintf(fd," }; \n");

fprintf(fd,"#define A 0\n"):
fprintf(fd,"#define B 1\n"):
fprintf(fd,"#define ATAN_BHI 2\n"):
fprintf(fd,"#define ATAN_BLO 3\n"):

fprintf(fd,"\n#ifdef WORDS_BIGENDIAN\n "):

for isbig from 1 to 0 by -1 do

    if(isbig=0) then
        fprintf(fd,"#else\n");
    fi:

    if(not (nb_of_ai = nb_of_bi)) then
        printf("Warning : nb_of_ai != nb_of_bi, this should not work");
    fi:

    fprintf(fd,"\n/* limits of the intervals [a[i],b[i]] */\n");
    fprintf(fd, "static db_number const arctan_table[%d][4] = \n{\n" , nb_of_ai );

   for i from 0 to nb_of_ai - 1 do
        fprintf(fd, "{\n/*a[%d]     */ ",i);
        printendian(fd,a[i],isbig);
        fprintf(fd," ,\n/*b[%d] :   */ ",i):
        printendian(fd,b[i],isbig):
        fprintf(fd," ,\n/*atan_b[%d]*/ ",i):
        printendian(fd, arctan(b[i]) ,isbig):
        fprintf(fd,",");
        printendian(fd, arctan(b[i])-nearest(arctan(b[i])) ,isbig);
        fprintf(fd," ,\n}\n,");
    od;
    fprintf(fd,"\n};\n");
od:

fprintf(fd,"\n#endif\n\n"):

fprintf(fd, "#endif /* def _ATAN_FAST_H */\n");

fclose(fd);

# Scs phase

#Now we want a precision of 130 bits in order to have correct rounding in all cases.
evalf(log[2](e^21/21));
#Test about Rounding :
#error on the polynom :
P_scs := convert(series(arctan(x),x=0,21),polynom);
evalf(log[2](e^21/21));
#Intervals parameters
#We choose the same e and the same intervals as in the quick phase.

#We have the same b[i]s and a[i]s but we save atan[b[i]] on 4 doubles

# Output
filename:="TEMPATAN/atan_accurate.h":
fd1:=fopen(filename, WRITE, TEXT):

fprintf(fd1, "#include \"crlibm.h\"\n#include \"crlibm_private.h\" \n#include \"atan_fast.h\"\n"):
fprintf(fd1, "\n/*File generated by maple/atan.mpl */\n"):

Write_SCS_poly(fd1, "constant_poly", P_scs);

fprintf(fd1,"#define constant_poly_ptr (scs_ptr)&constant_poly\n");

# The 1/Pi SCS constant
  fprintf(fd, "static const scs InvPiSCS=\n"):
  WriteSCS(fd, evalf(1/Pi)):
  fprintf(fd, ";\n#define InvPiSCS_ptr  (scs_ptr)(& InvPiSCS)\n\n"):

fprintf(fd1,"#ifdef WORDS_BIGENDIAN\n\n"):
for isbig from 1 to 0 by -1 do
   if isbig = 0 then
      fprintf(fd1, "\n#else\n"):
   fi:

   fprintf(fd1,"static const db_number atan_blolo[%d] = {\n",nb_of_ai);
   for i from 0 to nb_of_ai-1 do
      fprintf(fd1,"/* %d */  ",i):
      temp0 := nearest(arctan(b[i]));
      temp := nearest( arctan(b[i]) - nearest(arctan(b[i])));
      temp1 := nearest( arctan(b[i]) - temp -temp0);
      printendian(fd1,temp1,isbig):
      fprintf(fd1,", \n");
   od:

   fprintf(fd1,"};\n"):

od:
fprintf(fd1,"\n#endif  /* WORDS_BIGENDIAN */ \n");

fclose(fd1);

max_stars_repo_name: squarePenguin/parvsl

max_issues_repo_name: squarePenguin/parvsl

max_forks_repo_name: squarePenguin/parvsl


text:
# Various utilities used by Hakaru which don't fit logically elsewhere
Utilities := module ()
  option package;

  # Returns true for a sequence for which all the elements are equal up to the
  # specified equality, or false otherwise. A 0-length sequence produces an
  # error.
  export the := proc(as, eq:=`=`, $)
    if nops(as)=0 then error "sequence %1 must have at least one operand", as; end if;
    andmap(x -> eq(x,op(1,as)), [op(2..-1,as)]);
  end proc;

  # zip together an arbitrary number of lists;
  # - if any lists are shorter, the function recieves all available arguments
  #    e.g. zip_k(F, [a,b], [c]) -> [F(a, c), F(b)]
  # - if any lists contain subslists, those are also zipped together, so
  #   zip_k(F, as) and zip_k(F, [as]) are equivalent. To avoid this, use a 0-th
  #   operand other than `[]`.
  #    e.g. zip_k(F, [a,b], [ [c,d], [e,f] ]) -> [F(a, c, d), F(b, e, f)]
  export zip_k := proc(f)
    map(f@op@ListTools[Flatten], foldl((a,b)->zip(`[]`,a,b,[]), _rest));
  end proc;

  # Merges two records, `rspec' and `rdef'. If `rspec' contains any fields of
  # `rdef', the fields of `rspec' are taken.
  export merge_record_default := proc(rdef::record, rspec::record, $)
    local s;
    s   := {exports(rdef)} intersect {exports(rspec)};
    Record[ Record[rdef](op(map(`-`,s))) , rspec ]();
  end proc;

  # Creates a new record with `field' updated to the given function applied to
  # `field' of the given record.
  export upd_field := proc(r::record, field::name, k, $)
    Record[r](-field,field=k(r[field]))
  end proc;

  # Simplistic negation of relations. Used by Hakaru:-flatten_piecewise, KB, and
  # Domain.
  export negated_relation:= table([`<`, `<=`, `=`, `<>`] =~ [`>=`, `>`, `<>`, `=`]);

  # Takes the bool type (true/false) to mean universal and empty relations
  # respectively.  i.e. negate R, where R is an 'atomic' relation of a KB.
  # This is used by KB and Partition
  export negate_rel:= proc(R::Relation, $)::Relation;
    if R :: truefalse then
      not R
    elif R :: relation then
      negated_relation[op(0,R)](op(R));
    else
      # This is to appease 'piecewise', which won't be happy with Not
      # However, KB doesn't really care - it's already expecting {Not,not}
      # 'Technically' this is a KB 'constructor'!
      not(R);
    end if;
  end proc;

  # Differs from negate_rel in that
  #  [ negate_rel (x) ] = "Not x"  but
  #  [ negate_rels(x) ] = "x"
  export negate_rels := proc(e, $)
    subsindets(e, { specfunc(relation, `Not`), `not`(relation) }, negate_rel@op );
  end proc;

  # Creates an N-ary operator `Plus` such:
  #  Plus(Plus(a,b),c)=Plus(a,Plus(a,b))
  #  Plus(a,Iden)=a
  #  Plus(a,Zero)=Zero
  #  Plus(a)=a
  #  Plus(x,x)=x
  export Mk_Plus := proc(plus,iden,zero,$) proc()
    local as := {args};
    as := map(a->if op(0,a)=plus then op(a) else a end if, as);
    if zero in as then return zero end if;
    as := remove(`=`,as,iden);
    if nops(as)=0 then iden
    elif nops(as)=1 then op(1,as)
    else plus(op(as))
    end if;
  end proc; end proc;

  # Replacements for `and`, `or`, `not` which do not
  # evaluate "x = y" to "false" when "x,y" are unbound vars.
  export bool_And := Mk_Plus(And,true ,false);
  export bool_Or  := Mk_Plus(Or ,false,true);
  export bool_Not := proc(a,$)
    if a :: Relation and not (a :: `::`) then
      subsindets(negate_rel(a), `not`, Not@op);
    else
      Not(a)
    end if;
  end proc;

  # boolean_if should be equivalent to `if`, but it assumes
  # all its arguments are boolean conditions, so it basically
  # simplifies "cond and th or not cond and el"
  export boolean_if := proc(cond, th, el, $)
    use
      a = ((x,y)-> `if`(x=true,y, `if`(x=false,x,
                   `if`(y=true,x, `if`(y=false,y, And(x,y)))))),
      o = ((x,y)-> `if`(x=false,y, `if`(x=true,x,
                   `if`(y=false,x, `if`(y=true,y, Or (x,y)))))),
      n = (x    -> `if`(x=false,true,
                   `if`(x=true,false,             Not(x))))
    in
      o(a(cond,th), a(n(cond),el))
    end use
  end proc;

  # In order to hopefully produce a simplification,
  #   assert_deny will on some occasions repeatedly apply
  #   `simplify@ln` in the hopes of producing an 'improved'
  #   constraint. This metric gives the stopping condition
  # - when the simplification ceases to improve the constraint
  # - which is when the metric is made no less by the `simplify@ln`.
  # Since this is a `length`, this is strictly decreasing,
  #  so such a loop will 'provably' always terminate
  export log_metric := proc(e, x, $)
    local m, L;
    m := select(depends, indets(e, 'exp(anything)'), indets(x,name));
    length(subsindets(map2(op, 1, m), name, _->L));
  end proc;

  # Used by both KB and Domain
  export try_improve_exp := proc(b0, x, ctx, $)
    local b := b0, log_b;
    do
      try log_b := map(simplify@ln, b) assuming op(ctx); catch: break; end try;

      if log_metric(log_b, x) < log_metric(b, x)
         and (andmap(e->is(e,real)=true, log_b) assuming op(ctx)) then
        b := log_b;
      else
        break;
      end if;
    end do;
    b;
  end proc;

  # Given an expression containing products and sums, i.e. polynomials
  # and a function , applies this expression to each factor and summand
  export for_poly := proc(e, f, $)
    if e :: '{`+`,`*`}' then map(for_poly, e, f)
    elif e :: 'specfunc({product,Product,sum,Sum})' then
      applyop(for_poly, 1, e, f)
    else f(e)
    end if
  end proc;

  # Apply a function to each summand at the top level of the expression
  export distrib_over_sum := proc(f,e,$) `+`(op(map(f,convert(e, 'list',`+`)))) end proc;

  # Like convert(e, 'list', `*`) but tries to keep the elements positive
  export list_of_mul := proc(e, $)
    local rest, should_negate, can_negate, fsn;
    rest := convert(e, 'list', `*`);
    rest := map((f -> [f, signum(f),
                       `if`(f::'{specfunc({Sum,sum}),anything^odd}',
                            applyop(`-`,1,f),
                       `if`(f::'specfunc(ln)',
                            applyop(`/`,1,f),
                            -f))]),
                rest);
    should_negate, rest := selectremove(type, rest,
      '[anything, -1, Not(`*`)]');
    if nops(should_negate) :: even then
      [seq(op(3,fsn), fsn=should_negate),
       seq(op(1,fsn), fsn=rest)]
    else
      can_negate, rest := selectremove(type, rest,
        '[{`+`, specfunc({Sum,sum,ln}), `^`}, Not(1), Not(`*`)]');
      if nops(can_negate) > 0 then
        [seq(op(3,fsn), fsn=should_negate),
         op([1,3], can_negate),
         seq(op(1,fsn), fsn=subsop(1=NULL, can_negate)),
         seq(op(1,fsn), fsn=rest)]
      else
        [seq(op(3,fsn), fsn=subsop(-1=NULL, should_negate)),
         op([-1,1], should_negate),
         seq(op(1,fsn), fsn=rest)]
      end if
    end if
  end proc;

  export mysolve := proc(constraints)
    # This wrapper around "solve" works around the problem that Maple sometimes
    # thinks there is no solution to a set of constraints because it doesn't
    # recognize the solution to each constraint is the same.  For example--
    # This fails     : solve({c*2^(-1/2-alpha) = sqrt(2)/2, c*4^(-alpha) = 2^(-alpha)}, {c}) assuming alpha>0;
    # This also fails: solve(simplify({c*2^(-1/2-alpha) = sqrt(2)/2, c*4^(-alpha) = 2^(-alpha)}), {c}) assuming alpha>0;
    # But this works : map(solve, {c*2^(-1/2-alpha) = sqrt(2)/2, c*4^(-alpha) = 2^(-alpha)}, {c}) assuming alpha>0;
    # And the difference of the two solutions returned simplifies to zero.

    local result;
    if nops(constraints) = 0 then return NULL end if;
    result := solve(constraints, _rest);
    if result <> NULL or not (constraints :: {set,list}) then
      return result
    end if;
    result := mysolve(subsop(1=NULL,constraints), _rest);
    if result <> NULL
       and op(1,constraints) :: 'anything=anything'
       and simplify(eval(op([1,1],constraints) - op([1,2],constraints),
                         result)) <> 0 then
      return NULL
    end if;
    result
  end proc;

  # A debugging utility that's like `and` except it calls `userinfo` if there is
  # disagreement
  export and_info := proc(e :: {list,set})
    local s, r;
    s, r := selectremove(evalb, e);
    if nops(r) = 0 then return true end if;
    if nops(s) > 0 then userinfo(op([_rest, s, r])) end if;
    return false;
  end proc;

  # Simplifies relations using `Logic:-Normalize' and places the result in a
  # canonical form. The input can be {and,not,or} but the output is always
  # {And,Or,Not}.
  #  - false, true are expressed as And and Or
  #  - When `norty' is DNF/CNF then the outer constructor is always Or/And
  #
  # TODO:
  # A better interface for `simpl_relation` - something that allows it to treat
  #   arbitrary constructor forms as And, Or, Not. Typically, we do a `subs` right
  #   before the call to `simpl_relation` to put things in the right form, then
  #   another `subs` immediately inside `simpl_relation` to turn `And,Or` into `&and,&or`.
  export simpl_relation :=
  proc( expr_ :: {relation, boolean, specfunc({`And`,`Not`,`Or`}), `and`, `not`, `or`}
      , { norty := 'DNF' }
      , $) :: { specfunc(specfunc({Relation, specfunc(relation, Not)}, `Or`), `And`)
              , specfunc(specfunc({Relation, specfunc(relation, Not)}, `And`), `Or`)
              };
      local expr := expr_, outty, outmk, inty, inmk, ty_ord ;

      # Necessary because we convert `And' (etc.) to `Logic:-&and' (etc.), and
      # e.g. `Logic:-&and(a<b,c<d)' is not a valid piecewise condition.
      expr := chillFns('piecewise', expr);

      expr := foldr( proc(v,e) subsindets(e, op(v)) end proc
                   , expr
                   , [ { specfunc(relation, `Not`), `not`(relation) }
                     , x-> negate_rel(op(1,x)) ]
                   , [ { specfunc(`Not`), `not` }
                     , x->Logic:-`&not`(op(1,x)) ]
                   , [ { specfunc(`Or`), `or` }
                     , x->Logic:-`&or`(op(x)) ]
                   , [ { specfunc(`And`), `and` }
                     , x->Logic:-`&and`(op(x)) ] );
      expr := Logic:-Normalize(expr, form=norty);
      expr := foldr( proc(v,e) subsindets(e, op(v)) end proc
                   , expr
                   , [ specfunc(Logic:-`&and`), x->`And`(op(x)) ]
                   , [ specfunc(Logic:-`&or`) , x->`Or`(op(x)) ]
                   , [ specfunc(Logic:-`&not`), x->negate_rel(op(1,x))  ] );

      if expr :: identical(false) then
        return `if`(norty='DNF', `Or`(), `And`(`Or`()));
      elif expr :: identical(true) then
        return `if`(norty='DNF', `Or`(`And`()), `And`());
      end if;

      ty_ord := `if`(norty='DNF', [1,2], [2,1]);
      outty, inty := [ 'specfunc(Or)', 'specfunc(And)' ][ty_ord][];
      outmk, inmk := [ `Or`, `And` ][ty_ord][];

      if not expr :: outty then expr := outmk(expr) end if;
      expr := map(x -> if not x :: inty then inmk(x) else x end if, expr);
      expr := warmFns('piecewise', expr);
  end proc;

  # Classify a relation by matching either the LHS or RHS against
  # a given predicate, and placing the relation in a canonical form
  # where the LHS (3rd component of the result) is the side matching
  # the predicate, or FAIL of neither side matches.
  export classify_relation := proc(r0::relation, vars0, $)
    ::{identical(FAIL), [identical(B_LO,B_HI,B_EQ,B_NEQ),
                         satisfies(q->q in {indices(flip_relation,nolist)}),
                         anything, anything]};
    local r_k, r_s, in_vars, vars := vars0, r := r0;
    if vars :: ({set,list})({name,list(name)}) then
      vars := map(x->`if`(x::list,op(1,x),x),vars);
      vars := {op(vars)}; in_vars := x->x in vars;
    elif vars :: procedure then
      in_vars := vars;
    elif vars :: type then
      in_vars := x->type(x,vars);
    else error "unknown argument: %1 "
               "(expecting variables, variable membership check)"
               , vars;
    end if;
    if r :: {`<=`,`<`, `=`, `<>`} then
      if in_vars(rhs(r)) then
        r_k := flip_relation[op(0,r)]; r_s := rhs(r), lhs(r);
      elif in_vars(lhs(r)) then
        r_k := op(0,r); r_s := op(r);
      else return FAIL end if;
      [ classify_relop[r_k], r_k, r_s ];
    else FAIL
    end if;
  end proc;
  local flip_relation := table([`=`=`=`,`<>`=`<>`,`<=`=`>=`,
                                `<`=`>`,`>=`=`<=`,`>`=`<`]);
  local classify_relop := table([`=`=B_EQ,`<>`=B_NEQ,`<=`=B_HI,
                                 `<`=B_HI,`>=`=B_LO,`>`=B_LO]);

  # Functionalized, curried `if' statement
  export Case :=
  proc(ty,f,g)
    proc(x)
      if x::ty then f(x) else g(x) end if
    end proc
  end proc;

  # A replacement for `coulditbe .. assuming ..' which uses `is' internally
  # (since `is' is actually stronger than `coulditbe'); tries to handle `Or`s
  # correctly (which don't do well with `assuming'); and catches some exceptions
  # which we are reasonably sure have a valid interpretation.
  local do_rel_coulditbe := proc(a,as_::{set,list,Relation},$)
    option remember, system;
    local os, rs, as := as_;
    if as::{set,list} then
      as := And(op(as));
    end if;
    as := simpl_relation(as,norty='DNF');

    if nops(as)=1 then
      try
        not(is(bool_Not(a)) assuming op(1,as));
      catch "when calling '%1'. Received: 'contradictory assumptions'" :
        # technically this means the KB was already contradictory, we
        # just didn't know?
        return false;
      catch "when calling '%1'. Received: "
            "'side relations must be polynomials in"
            " (name or function) variables'":
        # This is seemingly a Maple bug - the condition could still be, but we
        # don't know, so conservatively return true.
        WARNING( sprintf( "siderels bug:\n\t'%s'\n"
                          "when calling coulditbe(%%1) assuming (%%2)"
                          , StringTools[FormatMessage](lastexception[2..-1])),
                 a, as_ );
        return true;
      catch "when calling '%3'. Received: "
            "'when calling '%2'. Received: "
            "'expression independent of, %0''":
        error expr_indp_errMsg(), a, as_;
      catch "when calling '%2'. Received: 'expression independent of, %0'":
        error expr_indp_errMsg(), a, as_;
      end try;
    else
      ormap(o1->rel_coulditbe(a,o1)=true, as);
    end if;
  end proc;

  local expr_indp_errMsg := proc($)
    sprintf("Something strange happened(%s)\n"
            "\tin coulditbe(%%1) assuming %%2"
            ,StringTools[FormatMessage](lastexception[2..-1]))
  end proc;

  # Given a symmetric relation `x0', applies `test' the operands of the relation
  # in a symmetric manner. If it succeeds, returns `x0' with the side satisifying
  # `test' on the left-hand. If neither side satisfies `test', returns FAIL.
  #  e.g.:
  #   is_lhs(x=3, (a,_)->type(a,integer))  =  3=x
  #   is_lhs(x=3, (a,_)->type(a,name))     =  x=3
  #   is_lhs(x=3, (a,_)->type(a,complex))  =  FAIL
  export is_lhs := proc(test,x)
    if test(lhs(x), rhs(x),_rest) then return x end if;
    if test(rhs(x), lhs(x),_rest) then return op(0,x)(rhs(x),lhs(x)) end if;
    FAIL;
  end proc;

  # `curry(chillFns,foo)' is the function which converts all occurences of
  # `foo(as)' to `f[as]' in the given expression. `warmFns' is the inverse.
  #
  # For some reason making these curried also requires `KB:-chill' and
  # `KB:-warm' to be eta-expanded (i.e. `chill := x->chillFns(chilled)(x)')
  export chillFns := (fns, e) -> subsindets(e, 'specfunc'(fns), c->op(0,c)[op(c)]);
  export warmFns  := (fns, e) -> subsindets(e, 'specindex'(fns), c->map(curry(warmFns,fns), op(0,c)(op(c))));

  # Print profiling information for a function when `infolevel' for that
  # function is set to at least 3, and assertlevel>0; e.g.
  #  > do_func := proc() .. func .. end: func := ProfileFn(do_func):
  #  > infolevel[func] := 3:
  #  > kernelopts(assertlevel) := 1:
  #  > func(a,b,c);
  #
  #  - profiling information is only printed for runs of the function
  #    which use a minimum amount of time, bytesalloc, or bytesused
  #  - setting infolevel >= 5 prints profiling information for all runs of the
  #    function.
  #  - less than a minimum amount of a resource consumed omits the actual usage
  #    with both infolevel >= 3 < 5 and infolevel >= 5 (just "less than <min>"
  #    is printed instead)
  #  - profiling information is not printed for recursives calls of a function
  #    (in other words, if a recursive function is profiled, only non-recursive
  #     instances will have profiling attached)
  export ProfileFn := module()
    local get_prof := proc()
      [ time[real](), kernelopts(bytesalloc)/2^20, kernelopts(bytesused)/2^20 ];
    end proc;
    local ppr_prof1 := proc(m,a,nm,fmt0,$)
      local v,p,fmt;
      v,p,fmt := `if`(a>=m, ["",a,fmt0], ["less than ",m,"%a"])[];
      cat(v,sprintf(fmt,p)," ",nm);
    end proc;
    local ppr_prof := proc(min_t, t, nms, $)
      local prefix, prefixl, strs;
      prefix := "took "; prefixl := length(prefix);
      strs :=
      zip_k(ppr_prof1, min_t, t, nms);
      cat(op(
        ListTools:-Interleave(
          ["\n\t"$nops(strs)],
          [prefix,cat(" "$prefixl)$(nops(strs)-1)],
          strs)));
    end proc;
    export ModuleApply :=
    proc(fn,min_t:=0.1,min_ba:=100,min_bu:=100)
      proc()
        local t, min_prof, profs, res, ctx, fncall;
        if kernelopts(assertlevel) > 0 and
           not (assigned(_Env_ProfileFn_inside[fn])) and
           assigned(infolevel['procname']) and
           infolevel['procname'] >= 3
        then
          _Env_ProfileFn_inside[fn] := true;
          t[0] := get_prof();
          res  := fn(args);
          t[1] := get_prof();
          t[2] := zip(`-`,t[1],t[0]);
          ctx  := map(op@getassumptions,indets([fn,args],Name));
          fncall := sprintf("%a(%q)",'procname',args);
          ctx    := sprintf("assuming (%q)", op(ctx));
          min_prof := [min_t,min_ba,min_bu];
          profs  := ppr_prof(min_prof, t[2],
                             [["seconds","%.3f"],
                              ["MiB alloc", "%.3f"],
                              ["MiB total used", "%.3f"]]);
          userinfo(`if`(`or`(zip(`>`,t[2], min_prof)[]),3,5),
                   'procname',
                   printf("Evaluating\n\t%s\n\t%s%s\n",fncall,ctx,profs));
          res;
        else
          fn(args)
        end if;
      end proc;
    end proc;
  end module;

  export rel_coulditbe := ProfileFn(do_rel_coulditbe, 1);

  # The closure of `f:set->set' with respect to `x:set'. That is, all values
  # which are outputs of `f' applied to `x' or an output of `f'.
  export cl := proc(f,x,$)
    local rs := f(x) minus x;
    if rs={} then x else cl(f,rs) union x end if;
  end proc;

  # The first `n' elements of a list, or the entire list if it has fewer than
  # `n' elements.
  export take := proc(xs,n,$)
    op(1,[ListTools[LengthSplit](xs,n)]);
  end proc;

  # A copy of `:-profile' with
  #  - better pretty printing of the result
  #  - bug fix for the fact that multiple names of sub-modules with the same
  #    'root' name cannot be profiled (e.g. M1:-x and M2:-x, commonly an issue
  #    with ModuleApply)
  #  - automatic collection of the names present in the Hakaru library
  # To profile `f(as)', call `Profile(f,_args=[as])'. Profiling information
  #  if printed to stdout
  # Extra arguments to `Profile' are passed to `GetProf' and `PPrProf', which
  #  determine which stats to profile and how to pretty print profiling
  #  information.
  # Environment variables affecting this function:
  #  _Env_Profile_count_ppr  = number of profiled names to print
  #  _Env_Profile_remove_ppr = selector function (: name->bool) to omit
  #    certain names from the profiling output
  export Profile := module()
    option package;
    export ModuleApply, GetProf, PPrProf, modules_to_profile, names_to_profile;
    local ModuleLoad, profile_flag_to_ord, name_to_string;

    modules_to_profile := proc()
      kernelopts(opaquemodules=false):
      { 'BindingTools'
      # , 'Hakaru'
      , 'KB'
      , 'Partition'
      , 'Partition:-Simpl'
      , 'Loop'
      , 'Domain'
      # , 'NewSLO'
      , 'NewSLO:-fromLO'
      , 'NewSLO:-toLO'
      , 'NewSLO:-reduce'
      , 'NewSLO:-reduce_Integrals'
      , 'NewSLO:-reduce_Partition'
      , 'NewSLO:-simplify_factor_assuming'
      , 'Summary'
      , entries(Domain:-Improve:-Simplifiers, nolist)
      };
    end proc;

    names_to_profile := proc()
      option remember, system;
      local ns;
      kernelopts(opaquemodules=false):
      ns := cl(curry(map,x->CodeTools:-Profiling:-getMemberFuncs(x,true)),
               modules_to_profile());
      map(proc(n)
            if n::`module`('ModuleApply') then 'n[ModuleApply]' elif n::procedure then n else NULL end if
          end proc, ns);
    end proc;

    name_to_string := (x->convert(x,'string'));

    profile_flag_to_ord := table(
       ['alpha' = (proc (a, b) lexorder(a[1],b[1]) end proc)
       ,'ralpha' = (proc (a, b) not lexorder(a[1],b[1]) end proc)
       ,'time' = (proc (a, b) evalb(a[4] < b[4]) end proc)
       ,'rtime' = (proc (a, b) evalb(b[4] < a[4]) end proc)
       ,'bytes' = (proc (a, b) evalb(a[6] < b[6]) end proc)
       ,'rbytes' = (proc (a, b) evalb(b[6] < a[6]) end proc)
       ,'load' = (proc (a, b) evalb(a[6]*a[6]*a[4] < b[6]*b[6]*b[4]) end proc)
       ,'rload' = (proc (a, b) evalb(b[6]*b[6]*b[4] < a[6]*a[6]*a[4]) end proc)
       ]);

    GetProf := proc(fns::({list,set})(satisfies(x->member(x,:-profile_profiled)))
                   ,{_flag:='rload'}
                   )
      local i, totaltime, totalbytes, totaldepth, totalcalls, timepc, bytespc, numargs, displist, totaltimepc, totalbytespc
           , ix, get_timepc, get_bytespc, get_nm, nm, flag ;
      global profile_time, profile_bytes, profile_maxdepth, profile_calls, profile_profiled;

      flag := _flag;
      totaltime, totalbytes, totalcalls, totaldepth := 0$4;
      numargs := nops(fns);

      for i to nops(profile_profiled) do
        ix := name_to_string(profile_profiled[i]);
        totaltime := totaltime+profile_time[ix];
        totalbytes := totalbytes+profile_bytes[ix];
        totalcalls := totalcalls+profile_calls[ix];
        totaldepth := totaldepth+profile_maxdepth[ix];
      end do;

      if totaltime = 0 then
        get_timepc := i->0;
        totaltimepc := 0;
      else
        get_timepc := i->100*profile_time[name_to_string(profile_profiled[i])]/totaltime;
        totaltimepc := 100;
      end if;
      for i to nops(profile_profiled) do
        timepc[name_to_string(profile_profiled[i])] := get_timepc(i);
      end do;

      if totalbytes = 0 then
        get_bytespc := i->0;
        totalbytespc := 0;
      else
        get_bytespc := i->100*profile_bytes[name_to_string(profile_profiled[i])]/totalbytes;
        totalbytespc := 100;
      end if;
      for i to nops(profile_profiled) do
        bytespc[name_to_string(profile_profiled[i])] := get_bytespc(i);
      end do;

      displist := [];
      if 0 < numargs then
        get_nm := i->op(i,fns);
      else
        numargs := nops(profile_profiled);
        get_nm := i->profile_profiled[i];
      end if;
      for i to numargs do
        nm := get_nm(i); ix := name_to_string(nm);
        displist := [op(displist),
                     [nm, map(q->q[ix],[profile_maxdepth, profile_calls,
                                        profile_time, timepc, profile_bytes, bytespc])[]  ]];
      end do;

      displist := sort(displist, profile_flag_to_ord[flag]);
      displist, [totaldepth,totalcalls,totaltime,totaltimepc,totalbytes,totalbytespc];
    end proc;

    PPrProf := proc(dat,tot,{_max_name_len:=55,_prof_name:=`function`})
      local i, k, pr_nm;
      k := _max_name_len; pr_nm := _prof_name;

      printf(cat(pr_nm,(" "$k-StringTools[Length](convert(pr_nm,string))),
                 `depth    calls     time    time%%     `));
      printf(`    bytes   bytes%%\n`);
      printf(cat("--"$k+1,`\n`));
      for i from 1 to nops(dat) do
        printf(cat(`%-`,convert(k,string),`a%7d%9d%9.3f%9.2f%14d%9.2f\n`),
               op(1..7,dat[i]));
      end do;
      printf(cat("--"$k+1,`\n`));
      printf(cat(`%-`,convert(k,string),`a%7d%9d%9.3f%9.2f%14d%9.2f\n\n`),
             `total:`,op(1..6,tot));
        return NULL
    end proc;

    ModuleApply := proc(f,{_args:=[]})
      local names_prof, res, nppr, rm_ppr, prf, tot, as; as := _args;
      if assigned(_Env_Profile_count_ppr) then
        nppr := _Env_Profile_count_ppr;
      else nppr := 25; end if;
      if assigned(_Env_Profile_remove_ppr) then
        rm_ppr := _Env_Profile_remove_ppr;
      else rm_ppr := (x->andmap(q->op(q,x)<0.001,[3,4,6])); end if;

      unprofile();
      names_prof := names_to_profile();
      profile(op(names_prof));
      if   f::string then res := NULL; read(f);
      else                res := f(op(as)); end if;
      prf, tot := GetProf(names_prof,_rest);
      unprofile();
      prf := remove(rm_ppr, prf);
      PPrProf(take(prf,nppr),tot,_rest);
      res;
    end proc;

    ModuleLoad := proc($)
      unprotect(:-convert);
      :-convert := overload([proc(x::name,_::identical(string),$) option overload(callseq_only); sprintf("%a",x) end proc
                            ,:-convert]);
      protect(:-convert);
      NULL;
    end proc;
    ModuleLoad();
  end module;

  local ModuleLoad := proc($)
    BindingTools:-load_types_from_table(
      table(
        [(Relation=
          ''{`::`, boolean, `in`, specfunc(anything,{Or,Not,And})}'')
        ,(Name = 'And(name, Not({constant,undefined}))')
        ]));
  end proc;
  ModuleLoad();
end module;

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


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 *)

theta0 := 1.0008:
theta1 := 0.1926:
theta2 := 1.8962:

f0 := s -> s^2/(1 + s)^2:
f  := x -> theta0 + f0(X2S*x)* (theta1 + f0(X2S*x) * theta2):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
 func $foo (
  var %i i32
#var %i1 i32, var %j1 i32, var %k1 i32
  ) i32 {
   return (
     extractbits i32 1 23(
       constval i32 -998))}

 func $foo1 (
  var %i i32, var %j i32, var %k i32,
  var %i1 i32, var %j1 i32, var %k1 i32
  ) i32 {
   return (
     extractbits i32 0 5(dread i32 %i))}

 func $foo11 (
  var %i i32, var %j i32, var %k i32,
  var %i1 i32, var %j1 i32, var %k1 i32
  ) u32 {
   return (
     extractbits u32 0 5(dread i32 %i))}

 func $foo2 (
  var %i i32, var %j i32, var %k i32
  ) i32 {
   return (
     extractbits i32 3 20(constval i32 0x111111111))}

 func $foo3 (
  var %i i64, var %j i64, var %k i32
  ) i64 {
   return (
     extractbits i64 5 50(dread i64 %i))}

 func $foo4 (
  var %i i64, var %j i64, var %k i32
  )i64  {
   return (
     extractbits i64 6 44 (dread i64 %i))}

 func $foo5 (
  var %i i64, var %j i64, var %k i32
  ) i64 {
   return (
     extractbits i64 1 10(constval i64 0x11111))}

 func $foo9 (
  var %i i8, var %j u8, var %k i32) u32 {
   return (
     extractbits u32 1 31(dread u32 %j))}

 func $foo51 (
  var %i u64
  ) u64 {
   return (
     extractbits u64 0 15(constval i64 0x11111))}

 func $foo52 (
  var %i i32
  ) i32 {
   return (
     sext i32 11 (constval i32 99))}

 func $foo53 (
  var %i u64
  ) u64 {
   return (
     zext u64 11 (constval u64 88))}

 func $foo54 (
  var %i i32
  ) i32 {
   return (
     sext i32 11 (dread i32 %i))}

 func $foo55 (
  var %i u64
  ) u64 {
   return (
     zext u64 41 (dread u64 %i))}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


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 *)

a1 := 199.81:
a2 := 4.3476:
c1 := 0.8524:
c2 := 1.2264:

# This is Eq. (40) of the paper.
f0 := s -> c1*(1 - exp(-a1*s^2)) + c2*(1 - exp(-a2*s^4)):

f  := x -> f0(X2S*x):


text:
with(Algebraic):

findElementaryIntegral := proc(maxint, maxdeg, raddeg, radexpt, x, ngoal)

	local ngiveup, nsincelast, totaltime, nfound, ntrials, nrepeats, num, den, rad, integrand, integral, radsqf, fac, r;

	Seed := randomize();

	ngiveup := 2500;

	randinrange := rand(-maxint..maxint);
	negrandinrange := rand(-maxint..-1);
	posrandinrange := rand(1..maxint);
	randf := rand(0.0..1.0);

	interface(prettyprint=0);

	found := [];
	nfound := 0;
	ntrials := 0;
	nrepeats := 0;
	nsincelast := 0;
	while (nfound < ngoal) do

		num := [seq(randinrange()*x^k, k=0..maxdeg)];
		num := convert(num, `+`);
		den := [seq(randinrange()*x^k, k=0..maxdeg)];
		den := convert(den, `+`);
		rad := [seq(randinrange()*x^k, k=1..raddeg-1)];
		rad := convert(rad, `+`);
		if randf() < 0.5 then
			rad := x^raddeg + negrandinrange() + rad;
		else
			rad := x^raddeg + posrandinrange() + rad;
		end if;

		if num*den*rad = 0 then next end if;

		radsqf := Squarefree(rad)[2];
		m := 1/radexpt;

		if nops(radsqf) = 1 then
			if degree(radsqf[1][1],x) = 1 then next end if;
		end if;

		bad := false;
		for fac in radsqf do
			r := fac[2] mod m;
			if r = 0 then
				bad := true;
			end if;
		end do;

		if bad then next end if;

		integrand := normal(num/den)/(rad)^radexpt;

		if evalb(integrand in found) then
			nrepeats += 1;
			if nrepeats > 4 then
				break
			end if;
			next
		end if;

		ntrials += 1;
#		print ("integrand = ", integrand);
		t0 := Now(ProcessClock);
		integral := int(convert(integrand, RootOf), x);
		t1 := Now(ProcessClock);
		totaltime := totaltime + t1 - t0;
#		print("integral = ", integral);

		if not has(integral,[int,EllipticPi,EllipticE,EllipticF])  then
			nfound += 1;
			nsincelast := 0;
			found := [op(found), integrand];
			print("integrand = ", integrand);
			print("integral  = ", convert(integral, radical));
		else
			nsincelast += 1;
		end if;

		if nsincelast > ngiveup then break end if;

		if ntrials mod 100 = 0 then
		   print("ntrials = ", ntrials, ", nfound = ", nfound, ", % elementary = ", evalf(100.0*nfound/ntrials));
		end if;
	end;

	print("FINISHED! params = ", maxint, maxdeg, raddeg, radexpt, x, ngoal, " ntrials = ", ntrials, ", nfound = ", nfound, ", % elementary = ", evalf(100.0*nfound/ntrials), ", mean time = ", evalf(totaltime/ntrials));
end proc;


text:
  #Disintegration *
  #---------------*
  #Abbreviations:
  #     wrt = "with respect to".
  #     wrt var = "a variable wrt which disintegration is performed". There may be
  #          more than one of these, in which case they are passed (in the 2nd arg)
  #          as Pairs, possibly nested.
  #     wrt-var type: Each wrt var may be continuous (Lebesgue), discrete (Counting),
  #          point evaluation (Dirac), and there may be other types added later.
  disint:= module()
  export ModuleApply, `do`;
  local
    #Dispatch table for wrt-var types. For each, there is a "cond constructor",
    #a "disintegrator", and a "disintegrator_arg_extractor". The cond constructor
    #builds the associated relation in the master piecewise; the disintegrator
    #does the differentiation or whatever operator is analogous to differentiation
    #for that type after the measure has been `improve`d; and the
    #disintegrator_arg_extractor builds the 2nd arg to the disintegrator from
    #the t_disint_var passed in the 2nd arg of disint.
    #
    #APIs:
    #     Indices:
    #          Must be the op(0, ...) of the expression passed with the
    #          wrt var.
    #     cond_constructor: (algebraic,name)-> boolean (most likely, a relation)
    #     disintegrator: (algebraic, {name, name=anything})-> algebraic
    #     disintegrator_arg_extractor:
    #          (`&M`(name, t_wrt_var_type))-> {name, name= anything}
    Wrt_var_types::static:= table([
         Lebesgue= Record(
              cond_constructor= `<=`,
              disintegrator= diff,
              disintegrator_arg_extractor= (A-> op(1,A)),
              disintegrator_type_extractor= [ (ll-> kb-> KB:-genLebesgue(op(1,ll), op(op(2,ll)), kb)), AlmostEveryReal ]
         ),
         Counting= Record(
              cond_constructor= `<=`,
              disintegrator= LREtools[delta],
              disintegrator_arg_extractor= (A-> op(1,A)),
              disintegrator_type_extractor= [ (ll-> kb-> KB:-genSummation(op(1,ll), op(op(2,ll)), kb)), HInt ]
         ),
         #Ret is aka Dirac.
         Ret= Record(
              cond_constructor= `=`,
              disintegrator='eval',
              disintegrator_arg_extractor= (A-> op(1,A)= op([2,1], A)),
              #E.g., (x &M Ret(3)) --> (x= 3).

              disintegrator_type_extractor= [ (ll -> kb -> (op(1,ll), kb)) ]
              # disintegrator_type_extractor= (ll-> kb-> KB:-genLet(op(1,ll), op([2,1],ll), kb))
         )
    ]),

    #types for disint wrt vars (2nd arg to disint)
    t_wrt_var_type,
    t_disint_var, #Curiosity: giving this a type `type` causes a kernel crash
                  #during update-archive.
    t_disint_var_pair,
    ModuleLoad::static:= proc($) #Needed to declare types.
         TypeTools:-AddType(
              t_wrt_var_type,
              satisfies(t-> assigned(Wrt_var_types[op(0,t)]))
         );
         TypeTools:-AddType(t_disint_var, {name, name &M t_wrt_var_type});
         TypeTools:-AddType(     #Caution: recursive type: Make sure base cases
           t_disint_var_pair, #are on left (a la McCarthy rule).
              'Pair'(Or(t_disint_var, t_disint_var_pair) $ 2)
         )
    end proc,
    #end of types for disint

    DV::table,  #wrt vars, with their types and conditions
    p::symbol,  #"pair"--abstract representation
    #`path` is layers of fst(...) and snd(...) built by traversing tree
    #(Weird Maple syntax note: Module prefixes seem to be required for
    #assertion type checking. Failure to include them causes kernel crash
    #during execution.)
    path::{specfunc({Hakaru:-fst, Hakaru:-snd}), symbol},

    #Parses the 2nd arg---the wrt vars.
    # works by side-effect: accumulates "paths" to variables in T
    # via the module variable DV.
    traverse_var_tree::static:= proc(
         T::{t_disint_var, t_disint_var_pair}, $
    )::identical(NULL);
    local
         v::name, #the wrt var
         M::NewSLO:-disint:-t_wrt_var_type,
         pp, #iterator over [fst, snd]---the deconstructors of Pair
         vM, vK, mks, mk_ty, mk_ctx
    ;
         if T::t_disint_var then
              #Add a default wrt-var type if none appears.
              (v,M):= op(`if`(T::name, T &M 'Lebesgue'((-1,1)*~infinity), T));
              vK   := op(0,M);
              vM   := v &M M;
              mks := Wrt_var_types[vK]:-disintegrator_type_extractor;
              mk_ty := op(1,mks);
              if nops(mks)=2 then mk_ctx := ((e)->lam(v,op(2,mks)(),e));
              else                mk_ctx := ((e)->e);
              end if;
              DV[v]:= Record(
                   'wrt_var_type'= M,
                   'path'= path,
                   'disintegrator_arg'=
                        Wrt_var_types[vK]:-disintegrator_arg_extractor(vM),
                   'disintegrator_mktype'=mk_ty(vM),
                   'disintegrator_mkctx' =mk_ctx
              );
              path:= op(path) #Peel off outer layer: fst or snd.
         else #T::Pair(..., ...)---deconstruct recursively.
              for pp in [fst, snd] do path:= pp(path); thisproc(pp(T)) end do
         end if;
         NULL
    end proc, #disint:-traverse_var_tree

    subs_disint_data := proc(sub_vars, expr, $)
        local V:= [indices(DV, 'nolist')], v
            , vars := [indices(sub_vars, 'nolist')]
            , eqns := [seq(v=sub_vars[v], v=vars)]
            , sb := (x->subs(eqns,x));

        for v in V do
            DV[v]:-disintegrator_arg := sb(DV[v]:-disintegrator_arg);
        end do;

        sb(expr);
    end proc;

  ;
    ModuleApply := proc()::t_Hakaru;
      local todo,expr,t;
      expr, todo := disint:-`do`(args)[];
      for t in todo do expr := t(expr); end do; expr;
    end proc;

   `do` := proc(
    m::t_Hakaru,
    #var &M wrt-var type, or Pairs thereof
    A::{t_disint_var, t_disint_var_pair},
    ctx::t_kb_atoms:= [] #context: parameter assumptions, "knowledge"
   ):: [t_Hakaru, list(appliable)];
   local
    mc,  #final integral to be passed to improve @ toLO; then result
         #of each disintegration step
    kb, var_rn := table(), mc_prts,
    V, #wrt vars
    v::name, #iterator over V
    improve_opts := [],
    todo := NULL,
    atodo := proc(z) todo := eval(todo), z; end proc,
    di, da
   ;
    if not {_rest} in {{'do_lam'},{}} then
      error "bad extra args: %1", {_rest};
    end if;

    #Init module variables.
    DV:= table();
    p:= gensym('p');
    path:= fst(p);

    traverse_var_tree(A); # collect information about A in DV
    V:= [indices(DV, 'nolist')];

    mc:= Bind(
         m, p,
         #The piecewise condition is a conjunction of conditions, each
         #condition formed from a (var,path) pair from DV.
         piecewise(
          And(seq(
               Wrt_var_types[op(0, DV[v]:-wrt_var_type)]:-
                cond_constructor(DV[v]:-path, v),
               v= V
          )),
          Ret(snd(p)),
          Msum()
         ));

    kb := empty;
    for v in ListTools[Reverse](V) do
        var_rn[v], kb := DV[v]:-disintegrator_mktype(kb);
    end do;
    kb := build_kb(ctx,"disint",kb);
    mc := subs_disint_data(var_rn, mc);

    atodo(x->toLO(x));
    atodo(x->improve(x, _ctx=kb,improve_opts));
    #Does the order of application of the disintegrators matter?
    #Theoretically, I think not, it's just like differentiation. As far
    #as Maple's ability to do the computation, maybe it matters.
    for v in V do
      di   := Wrt_var_types[op(0, DV[v]:-wrt_var_type)]:-disintegrator;
      da   := DV[v]:-disintegrator_arg;
      atodo(x->applyop(di, 2, x, da));
      atodo(x->improve(x, _ctx=kb,improve_opts));
    end do;

    atodo(x->fromLO(x, _ctx= kb));
    ## simplify integrals which do not mention the LO var (i.e. integrals in
    ## weights). this is a hack, we should do this inside of `improve'* in
    ## the correct place.  It is required to see the correct output for
    ## d7_normalFB1.
    #
    ## * Or, some other simplifier, which does only specific things - note
    ## the simplification we hope to see here (in d7) can only be done after
    ## the application of the `diff'. So just chucking it into `improve'
    ## might not be the correct thing to do. Calling `improve' at all might
    ## be wrong since we roundtrip through pulling off domains and replacing
    ## them (this is sort of expensive..), although we may also want such
    ## simplifications here - the application of the `diff' might give a new
    ## domain problem that can be improved significanly.
    # 'Simplify' partitions in weights by converting them to piecewise and
    # letting piecewise simplify do the work; and evaluate integrals in weights
    atodo(x->
          subsindets( x, 'Weight(anything, anything)'
                    , mw -> applyop(x->simplify_assuming(
                        subsindets( subsindets(x,Partition,PartitionToPW)
                                  , And(specfunc(`Int`)
                                       ,satisfies(e->indets(e,specfunc(`exp`))<>{})
                                       )
                                  , value
                                  )
                        , kb ), 1, mw)));

    atodo(x->subs_disint_data( table([seq(var_rn[v]=v, v=V)]), x));

    if 'do_lam' in {_rest} then
      for v in ListTools[Reverse](V) do
        atodo(x->DV[v]:-disintegrator_mkctx(x));
      end do;
    end if;

    [mc, [todo]];
   end proc; #disint:-ModuleApply
   ModuleLoad();
  end module; #disint
  ###################### end of Carl's code ######################

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/hakaru


text:
# File: RigidMotionsParameterSpaceDecompostionRecursive.mpl
#
# Description:
#  This file contains functions used to obtain an arrangement 6 dimensional parameter space of 3D
#  digitized rigid motions.
#  This code has been written for research propose and its aim is to calculate a particular
#  arrangement of quadrics. Therefore, it can or it cannot be useful in study of generic
#  arrangements. The final output are sample points of full dimensional open cells.
#
#  The code was written in relation with the paper: Kacper Pluta, Guillaume Moroz, Yukiko
#  Kenmochi, Pascal Romon, Quadric arrangement in classifying rigid motions of a 3D digital image,
#  2016, https://hal.archives-ouvertes.fr/hal-01334257 referred late on as [Quadrics:2016].
#
# Author:
#  Kacper Pluta - kacper.pluta@esiee.fr
#  Laboratoire d'Informatique Gaspard-Monge - LIGM, A3SI, France
#
# Date:
#  12/24/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.
#
#
RigidMotionsParameterSpaceDecompostionRecursive := module()
  option package;

  (*Threshold which controls when to synchronize databases.*)
  local RECORDS_TO_SYNCH := 1000;

  uses   RigidMotionsParameterSpaceCommon;
  local  ComputeEventsAType2D, ComputeEventsBType2D, ComputeEventsAType1D,
         ComputeAsymptoticAAEvents2D, ComputeEventsAlgebraicNumbers2D, ComputeSamplePoints2D;

  export IsAsymptotic2D, LaunchComputeSamplePoints2D;

# Procedure: ComputeEventsAType2D
#   Compute events such that a sweep line is tangent to a conic with respect to the first provided
#   variable.
#
# Parameters:
#   Q2D2       - a set of conics
#   grid       - if true then computations are performed in the parallel computation framework
#                called Grid (see Maple documentation)
#   vars2D       - a list of the variables
#
# Output:
#   list of couples: univariate polynomial and index of the generating conic: [poly, [index]]
#
# Comment/Limitations:
#  - Univariate polynomials are expressed in the first variable.
ComputeEventsAType2D := proc(Q2D2, grid::boolean, vars2D::list)
  uses ArrayTools;
  local result := Array([]), s;
  if nops(vars2D) < 2 then
    error "Only systems in at least two variables are supported. The variables are: %1.", vars2D;
  fi:
  s := proc(i::integer, vars2D::list)
    local sys, univ, sol:
    local q := Q2D2[i];
    sys := [q, diff( q, vars2D[2] )];
    univ := UnivariatePolynomial(sys, vars2D):
    return SerializeEvents(GenerateEvents(univ, [i]));
  end proc;
  if grid then
    map[inplace](proc(x) Extend(result, x, inplace=true) end proc, [Grid:-Seq( s( i, vars2D ),
                                                                   i=1..nops( Q2D2 ) )]);
  else
    map[inplace](proc(x) Extend(result, x, inplace=true) end proc, [seq( s( i, vars2D ),
                                                             i=1..nops( Q2D2 ) )]);
  fi;
  return ReconstructEvents(result);
end proc:

# Procedure: ComputeEventsBType2D
#   Compute events such that two conics intersects in a point.
#
# Parameters:
#   Q2D2       - a set of conics
#   grid       - if true then computations are performed in the parallel computation framework
#                called Grid (see Maple documentation)
#   vars2D       - a list of the variables
#
# Output:
#   list of couples: univariate polynomial and indices of the generating conics: [poly,
#   [index1,index2]]
#
# Comment/Limitations:
#  - Univariate polynomials are expressed in the first variable.
ComputeEventsBType2D := proc(Q2D2, grid::boolean, vars2D::list)
  uses ArrayTools;
  local result := Array([]), s:
  s := proc (i, j, vars2D::list)
    local p, sol, univ, sys;
    sys := [Q2D2[i], Q2D2[j]];
    univ := UnivariatePolynomial(sys, vars2D);
    return SerializeEvents(GenerateEvents(univ, [i, j]));
  end proc;
  if grid then
    map[inplace](proc(x) Extend(result, x, inplace=true) end proc, [Grid:-Seq( seq( s(i, j, vars2D),
                                                       j=i+1..nops( Q2D2 ) ), i=1..nops( Q2D2 ) )]);
  else
    map[inplace](proc(x) Extend(result, x, inplace=true) end proc, [seq( seq( s(i, j, vars2D),
                                                 j=i+1..nops( Q2D2 ) ), i=1..nops( Q2D2 ) )]);
  fi;
  return ReconstructEvents(result);
end proc:

# Procedure: IsAsymptotic2D
#   Checks if a curve is an asymptote.
#
# Parameters:
#   p          - a curve given as a polynomial
#   var        - variable to be reduced
#
# Output:
#   Monomial such that there exists a line tangent to the curve at infinity represented in a
#   variable different than 'var' or a constant.
IsAsymptotic2D := proc(p::polynom, var)
  return lcoeff(p, var);
end proc:

# Procedure: ComputeAsymptoticAAEvents2D
#   Compute real algebraic numbers which corresponds to asymptotic cases given by one curve.
#
# Parameters:
#   Q2D2          - a set of conics
#   vars2D       - a list of the variables
#
# Output: Real
#   A list of Events (see EventType).
ComputeAsymptoticAAEvents2D := proc(Q2D2, vars2D::list)
  uses ArrayTools;
  local result := Array([]), s;
  s:=proc(i::integer, vars2D::list)
    local asy;
    asy := RigidMotionsParameterSpaceDecompostionRecursive:-IsAsymptotic2D(Q2D2[i], vars2D[-1]);
    return GenerateEvents(asy, [i]);
  end proc:
  map[inplace](proc(x) Extend(result, x, inplace=true) end proc, [seq(s(i, vars2D),
                                                                i=1..nops(Q2D2))]);
  return result;
end:

# Procedure: ComputeEventsAlgebraicNumbers2D
#   Compute and sort events
#
# Parameters:
#   Q2D2     - set of conics
#   grid       - if true then computations are performed in the parallel computation framework
#                called Grid (see Maple documentation)
#   vars2D       - a list of the variables
#
# Output:
#   Sorted Array of real algebraic numbers
ComputeEventsAlgebraicNumbers2D := proc(Q2D2, grid::boolean, vars2D::list)
  uses ArrayTools;
  local events2D := Array([]);
  Extend(events2D, ComputeEventsAType2D(Q2D2, grid, vars2D), inplace=true);
  Extend(events2D, ComputeEventsBType2D(Q2D2, grid, vars2D), inplace=true);
  Extend(events2D, ComputeAsymptoticAAEvents2D(Q2D2, vars2D), inplace=true);
  return AlgebraicSort(events2D);
end proc:

# Procedure: ComputeEventsAType1D
#   Compute and sort events as algebraic numbers
#
# Parameters:
#   Q2D2     - a list of conics
#
# Output:
#   Sorted list of real algebraic numbers without repetitions.
ComputeEventsAType1D := proc(Q2D2::list)
  local q, factored, sqrFree, rootsF, rf, numbers := Array([]);
  for q in Q2D2 do
    if RootFinding:-HasRealRoots(q) then
      factored := factors( q )[2,..,1];
      for sqrFree in factored do
        rootsF := RootFinding:-Isolate(sqrFree, output='interval');
        for rf in rootsF do
          ArrayTools:-Append(numbers, RealAlgebraicNumber(sqrFree, op(rf)[2][1],
                             op(rf)[2][2]));
        od;
      od;
   fi;
  od;
  numbers := AlgebraicSort(numbers);
  return ListTools:-MakeUnique(convert(numbers, list), 1, proc(a,b) evalb(Compare(a,b) = 0) end proc);
end proc:

# Procedure: ComputeSamplePoints2D
#   Computes sample points for rotational part of rigid motions
#
# Parameters:
#   Q2D         - a list of conics
#   events2D    - each element contains a list of algebraically unique events
#   first       - integer value which indicates the first event to proceed.
#   last        - integer value which indicates the last event to proceed.
#   vars2D      - list of variables in which conics are expressed
#   db          - an instance of the class ComputationRegister
#
# Output:
#   It populates a database, given by databasePath, with sample points.
ComputeSamplePoints2D := proc(Q2D, events2D::Array, first::integer, last::integer,
                              vars2D::list, amid::rational, db::ComputationRegister)
  local i::integer, j::integer, x::list, midpoint::rational, sys::list, records := 0;
  local disjointEvent::list, oneD::list, oneDNeg::list, ranumI, ranumJ;
  if first < 0 or last < 0 or last < first or upperbound(events2D) <= last then
    error "Bounds of the array range are incorrect.":
  end if:
  for i from first to last do
    sys := Q2D[GetQuadrics(events2D[i])];
    ranumI := GetRealAlgebraicNumber(events2D[i]);
    ranumJ := GetRealAlgebraicNumber(events2D[i+1]);
    disjointEvent:=DisjointRanges(ranumI, ranumJ);
    midpoint := (GetInterval(disjointEvent[1])[2] + GetInterval(disjointEvent[2])[1])/2:
    # never call eval with sets!
    sys := eval(sys, vars2D[1] = midpoint);
    oneD := ComputeEventsAType1D(sys);
    oneD := select(proc(x) return evalb(GetInterval(x)[2] >= 0); end proc, oneD);

    if upperbound(oneD) > 0 then
      for j from 1 to upperbound(oneD) - 1 do
        disjointEvent := DisjointRanges(oneD[j],oneD[j+1]);
        InsertSamplePoint(db, [amid, midpoint, (GetInterval(disjointEvent[1])[2] +
                                            GetInterval(disjointEvent[2])[1])/2]);
        records := records + 1;
        if records mod RECORDS_TO_SYNCH = 0 then
          SynchronizeSamplePoints(db);
        fi;
      od:
      InsertSamplePoint(db, [amid, midpoint, GetInterval(oneD[-1])[2] + 1/2]);
      records := records + 1;
      if records mod RECORDS_TO_SYNCH = 0 then
        SynchronizeSamplePoints(db);
      fi;
    fi;
 od:
end proc:

# Procedure: ComputeSamplePoints2D
#   Computes sample points for rotational part of rigid motions using the grid framework
#
#
# Parameters:
#   s         - a list of conics
#   midpoint  - the first dimensional midpoint obtained from the 3D decomposition
#   nodes     - number of nodes used in the parallel computations
#   grid      - a control variable for parallel computations. If true and additional conditions on
#               the size of the problem are fulfilled the problem is solved in the grid framework
#   variables - list of variables in which the problem is expressed
#   db        - an instance of ComputationRegister which provides interface to the database
# Output:
#   It populates a database with sample points.
LaunchComputeSamplePoints2D := proc (s::list, midpoint::rational, nodes::integer,
                                           grid::boolean, variables::list, db::ComputationRegister)
  local events2D, Q2D := ListTools:-MakeUnique([op(variables), op(s)]);
  if grid and nops(s) > 20 then
     events2D := ComputeEventsAlgebraicNumbers2D(Q2D, true, variables);
  else
     events2D := ComputeEventsAlgebraicNumbers2D(Q2D, false, variables);
  fi;
  events2D := remove[flatten](proc(x) evalb(GetInterval(GetRealAlgebraicNumber(x))[2] < 0) end proc,
                            events2D);
  if upperbound(events2D) = 0 then
    return NULL;
  fi;
  events2D := ReduceEvents(events2D);
  AdjustEvents(events2D, upperbound(Q2D), variables);
  ComputeSamplePoints2D(Q2D, events2D, 1, upperbound(events2D) - 1, variables, midpoint, db);
end proc:

end module;


text:
type $RR <class {@e1 i32, @e2 f32, @e3 f64}>
type $SS <class <$RR> { &method1(agg)agg,
                        &method2(void)void}>
func $foo (
	var %x <$SS>) i32 {
  dassign %x 2 ( constval i32 32 )
  virtualcall &method2(addrof ptr %x)
  superclasscall &method2(addrof ptr %x)
  interfacecall &method2(addrof ptr %x)
  return ( dread i32 %x 2 ) }
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
#@ Not autoload

with(LinearAlgebra):

p := 2;
adams_n_max := 5;
adams_s_max := 20;

unprotect('t');
unassign('t');

xi[0] := 1;

protect('t');

HH_vars := plex(seq(xi[i],i=1..adams_n_max));

# The cobar complex involves variables xi[n,s] corresponding to the copy of
# xi[n] in the s'th tensor factor.  We use the ordering with
# xi[1,1] >> xi[2,1] >> xi[3,1] >> ... >> xi[1,2] >> xi[2,2] >> xi[3,2] >> ...
# so the variables from each tensor factor dominate those from the next tensor
# factor.  It is not clear whether this is the best choice.
H_cobar_vars :=
 plex(seq(seq(xi[i,j],i=1..adams_n_max),j=1..adams_s_max));

HH_cmp := (a,b) -> TestOrder(a,b,HH_vars):
H_cobar_cmp := (a,b) -> TestOrder(a,b,H_cobar_vars):

H_degree_rule := {
  seq(xi[n] = e^(2^n-1) * xi[n],n=1..adams_n_max),
  seq(seq(xi[n,i] = e^(2^n-1) * xi[n,i],n=1..adams_n_max),i=1..adams_s_max)
}:

H_degree := (u) -> degree(subs(H_degree_rule,u),e);

H_bidegree := proc(u)
 local s,t;
 t := H_degree(u);
 s := adams_s_max;
 while s > 0 and not(has(u,{seq(xi[i,s],i=1..adams_n_max)})) do
  s := s-1;
 od;
 return [s,t];
end:

# Basis for Z[xi[i] : i >= n] in degree d
T_basis := proc(d::integer,n::posint := 1)
 option remember;
 local m,r,i;
 m := 2^n-1;
 r := floor(d/m);
 if d < 0 then
  return [];
 elif d = 0 then
  return [1];
 elif r = 0 then
  return [];
 else
  map(op,[seq(xi[n]^i *~ T_basis(d-m*i,n+1),i=0..r)]);
 fi;
end:

# Basis for Z[v[i] : i >= n] in degree <= d
T_lower_basis := proc(d::integer,n::posint := 1)
 option remember;
 local m,r,i;
 m := 2^n-1;
 r := floor(d/m);
 if d < 0 then
  return [];
 elif r = 0 then
  return [1];
 else
  map(op,[seq(xi[n]^i *~ T_lower_basis(d-m*i,n+1),i=0..r)]);
 fi;
end:

# Basis for the s-fold tensor power of Z[t_1,t_2,...] in degree d
# The copy of xi[i] in the j'th tensor factor is represented by xi[i,j]
T_power_basis := proc(s::nonnegint,d::integer)
 local B,m,R1,R2,u,v;
 if s = 0 then return `if`(d = 0,[1],[]); fi;

 R1 := {seq(xi[n] = xi[n,1],n=1..adams_n_max)};
 R2 := {seq(seq(xi[n,i] = xi[n,i+1],i=1..s-1),n=1..adams_n_max)};
 B := NULL;

 for u in subs(R1,T_lower_basis(d)) do
  m := d - H_degree(u);
  for v in subs(R2,T_power_basis(s-1,m)) do
   B := B,(u*v);
  od;
 od;
 return [B];
end:

# Basis for the s-fold tensor power of the augmentation ideal in
# Z[t_1,t_2,...] in degree d
T_reduced_power_basis := proc(s::nonnegint,d::integer)
 local B,m,R1,R2,u,v;
 if s = 0 then return `if`(d = 0,[1],[]); fi;

 R1 := {seq(xi[n] = xi[n,1],n=1..adams_n_max)};
 R2 := {seq(seq(xi[n,i] = xi[n,i+1],i=1..s-1),n=1..adams_n_max)};
 B := NULL;

 for u in subs(R1,T_lower_basis(d)) do
  if u <> 1 then
   m := d - H_degree(u);
   for v in subs(R2,T_reduced_power_basis(s-1,m)) do
    B := B,(u*v);
   od;
  fi;
 od;
 return [B];
end:

H_cobar_basis := (s,d) -> T_reduced_power_basis(s,d);

# Hopf algebroid coproduct on the generators xi[n]
psi_xi := proc(n)
 local a,b;
 option remember;
 if n = 0 then return 1; fi;

 return xi[n,1] + xi[n,2] + add(xi[n-i,1]^(2^i)*xi[i,2],i=1..n-1);
end:

d_H_cobar_rule := proc(s,i)
 local R0,R1,R2,R3;
 if i = 0 then
  R0 := {seq(xi[j]=xi[j,1],j=1..adams_n_max)};
  R1 := {seq(xi[j]=xi[j,2],j=1..adams_n_max)};
  R2 := {seq(seq(xi[j,k]=xi[j,k+1],k=1..s),j=1..adams_n_max)};
  return {op(R1),op(R2)};
 else
  R0 := {seq(xi[j,1]=xi[j,i],j=1..adams_n_max),
         seq(xi[j,2]=xi[j,i+1],j=1..adams_n_max)};
  R1 := {seq(seq(xi[j,k]=xi[j,k+1],k=i+1..s),j=1..adams_n_max)};
  R2 := {seq(xi[j,i] = subs(R0,psi_xi(j)),j=1..adams_n_max)};
  if i = 1 then
   R3 := {seq(xi[j] = modp(expand(subs(R0,psi_xi(j))),2),j=1..adams_n_max)};
  else
   R3 := {seq(xi[j] = xi[j,1],j=1..adams_n_max)};
  fi;
  return {op(R1),op(R2),op(R3)};
 fi;
end:

d_H_cobar := (s) -> (u) ->
 modp(expand(add(subs(d_H_cobar_rule(s,i),u),i=0..s+1)),2);

d_H_cobar_matrix := proc(s,d)
 local B1,B2,cf;
 B1 := H_cobar_basis(s,d);
 B2 := H_cobar_basis(s+1,d);
 cf := proc(u)
  local sol;
  sol := solve({coeffs(u - add(c[i]*B2[i],i=1..nops(B2)),indets(B2))});
  subs(sol,[seq(c[i],i=1..nops(B2))]);
 end;
 Transpose(Matrix(map(cf,map(d_H_cobar(s),B1))));
end:

mu_H_cobar := (s1,s2) -> proc(a,b)
 local i,b0;
 b0 := b;
 for i from 0 to s1-1 do
  b0 := modp(expand(subs(d_H_cobar_rule(s2+i,0),b0)),2);
 od:
 return modp(expand(a * b0),2);
end:

analyse_H_cobar := proc(s::nonnegint,d::integer)
 local i,R,R0,B1,B2,n1,n2,M,L0,P0,x,x1,L,P,Q,nx,LB2,y,HB,HE,T,U,L1,K,Ki;
 global H_cobar_data;

 if d + s >= 2*(p^(adams_n_max + 1) - 1) then
  error("adams_n_max is too small");
 fi;

 if s > adams_s_max then
  error("adams_s_max is too small");
 fi;

 R := table():

 if s = 0 then
  if d = 0 then
   R["chain_basis"]     := [1];
   R["chain_rank"]      := 1;
   R["cycle_basis"]     := [1];
   R["cycle_rank"]      := 1;
   R["boundary_basis"]  := [];
   R["boundary_rank"]   := 0;
   R["pivot_data"]      := [[1,infinity]];
   R["non_cycle_basis"] := [];
   R["homology_basis"]  := [1];
   R["homology_exponents"] := [infinity];
  else
   R["chain_basis"]     := H_cobar_basis(0,d);
   R["chain_rank"]      := nops(R["chain_basis"]);
   R["cycle_basis"]     := [];
   R["cycle_rank"]      := 0;
   R["boundary_basis"]  := [];
   R["boundary_rank"]   := 0;
   R["pivot_data"]      := [];
   R["non_cycle_basis"] := R["chain_basis"];
   R["homology_basis"]  := [];
   R["homology_exponents"] := [];
  fi;

  H_cobar_data[s,d] := eval(R);
  return eval(R);
 fi;

 if type(H_cobar_data[s-1,d+1],table) then
  R0 := H_cobar_data[s-1,d+1];
  B1 := R0["non_cycle_basis"];
 else
  B1 := H_cobar_basis(s-1,s+d):
 fi;

 B2 := H_cobar_basis(s,s+d):
 n1 := nops(B1);
 n2 := nops(B2);
 R["chain_basis"] := B2;
 R["chain_rank"]  := n2;

 if n2 = 0 then
  R["cycle_basis"]     := [];
  R["cycle_rank"]      := 0;
  R["boundary_basis"]  := [];
  R["boundary_rank"]   := 0;
  R["pivot_data"]      := [];
  R["non_cycle_basis"] := [];
  R["homology_basis"]  := [];
  R["homology_exponents"] := [];
  H_cobar_data[s,d] := eval(R);
  return eval(R);
 fi;

 M := Transpose(Matrix(map(coeff_list,map(d_H_cobar(s-1),B1),B2)));
 L0,P0,x := op(Zpl_reduce(Transpose(M),p)):
 L := Transpose(L0):
 P := Transpose(P0):
 nx := nops(x):
 R["cycle_rank"] := nx;

 LB2 := convert(L0.Vector(B2),list):
 R["boundary_basis"] := [seq(LB2[i],i=1..nx)];
 R["cycle_basis"]    := [seq(LB2[i]/p^xi[i][2],i=1..nx)];
 R["pivot_data"]     := x;

 y := sort([op({seq(i,i=1..n2)} minus {seq(xi[i][1],i=1..nx)})]);
 R["non_cycle_basis"] := [seq(B2[i],i in y)];

 HB := [];
 HE := [];

 for i from 1 to nx do
  if xi[i][2] > 0 then
   HB := [op(HB),R["cycle_basis"][i]];
   HE := [op(HE),p^xi[i][2]];
  fi;
 od;

 R["homology_basis"] := HB;
 R["homology_exponents"] := HE;

 T := <IdentityMatrix(nx)|Matrix(nx,n1-nx)>;
 U := Matrix(nx,n2):
 for i from 1 to nx do U[i,xi[i][1]] := 1; od:
 L1 := U.L.Transpose(T);
 Q := P.Transpose(T).(1/L1).U;
 K := <SubMatrix(L,1..n2,1..nx)|Matrix(n2,n2-nx)>;
 y := sort([op({seq(i,i=1..n2)} minus {seq(xi[i][1],i=1..nx)})]);
 for i from 1 to n2 - nx do
  K[y[i],nx+i] := 1;
 od:
 Ki := 1/K;

 H_cobar_data[s,d] := eval(R);
 return eval(R);
end:

save_H2_data := proc()
 local file;
 file := sprintf("%s/H_2.m",data_dir);
 save(adams_n_max,adams_s_max,psi_xi,H_basis,H_lower_basis,file);
end:

load_H2_data := proc(p)
 local file;
 file := sprintf("%s/H_2.m",data_dir,p);
 load(file);
end:

save_H2_cobar_data := proc()
 local file;
 file := sprintf("%s/H_cobar_data_2.m",data_dir);
 save(H_cobar_data,file);
end:

load_H2_cobar_data := proc()
 local file;
 file := sprintf("%s/H_cobar_data_2.m",data_dir);
 load(file);
end:

adams_representative := proc(u)
 local a,b,c,sa,sb,v,i;
 if type(u,integer) then
  return modp(u,2);
 elif type(u,`+`) then
  return modp(map(adams_representative,u),2);
 elif type(u,`*`) then
  c,v := selectremove(type,u,integer);
  if modp(c,2) = 0 then
   return 0;
  else
   v := sort([op(v)]);
   a := adams_representative(v[1]);
   b := adams_representative(mul(v[i],i=2..nops(v)));
   sa := H_bidegree(a)[1];
   sb := H_bidegree(b)[1];
   return modp(mu_H_cobar(sa,sb)(a,b),2)
  fi;
 elif type(u,`^`) and type(op(2,u),posint) then
  a := adams_representative(op(1,u));
  b := a;
  sa := H_bidegree(a)[1];
  for i from 2 to op(2,u) do
   b := modp(mu_H_cobar(sa,(i-1)*sa)(a,b),2);
  od;
  return b;
 else
  return procname(args);
 fi;
end:

adams_assassin := proc(err)
 local s,t,B0,B1,M0,V0,MV,sol;
 s,t := op(H_bidegree(err));
 B0 := H_cobar_basis(s-1,t);
 B1 := H_cobar_basis(s,t);
 M0 := d_H_cobar_matrix(s-1,t);
 V0 := Vector(coeff_list(err,B1));
 MV := LinearAlgebra[Modular][Mod](2,<M0|V0>,integer[]):
 sol := LinearAlgebra[Modular][LinearSolve](2,MV,1,inplace=false):
 return ([sol][1] . Vector(B0));
end:

adams_representative(h[0]) := xi[1,1];
adams_representative(h[1]) := xi[1,1]^2;
adams_representative(h[2]) := xi[1,1]^4;
adams_representative(h[3]) := xi[1,1]^8;
adams_representative(c[0]) := xi[1,1]^2*xi[1,2]^7*xi[1,3]^2+xi[1,1]^2*xi[1,2]^4*xi[1,3]^2*xi[2,2]-xi[1,2]^4*xi[1,3]^4*xi[2,1]+xi[1,1]^2*xi[1,2]^3*xi[2,3]^2+xi[1,1]^2*xi[1,2]*xi[1,3]^2*xi[2,2]^2+xi[1,1]^2*xi[1,2]*xi[1,3]^2*xi[2,3]^2+xi[1,1]^2*xi[1,3]^2*xi[3,2]+xi[1,1]^2*xi[2,2]*xi[2,3]^2;

adams_relations := [
 [h[0]*h[1],xi[1,1]^3+xi[2,1]],
 [h[0]^2*h[2]+h[1]^3,xi[1,1]^2*xi[1,2]*xi[2,1]+xi[1,1]^2*xi[1,2]*xi[2,2]+xi[2,1]*xi[2,2]],
 [h[1]*h[2],xi[1,1]^6+xi[2,1]^2],
 [h[0]*h[2]^2,xi[1,1]^5*xi[1,2]^4+xi[1,1]^2*xi[1,2]^7+xi[1,1]^2*xi[1,2]^4*xi[2,1]+xi[1,1]^2*xi[1,2]^4*xi[2,2]+xi[1,2]^6*xi[2,1]+xi[1,1]^2*xi[1,2]*xi[2,2]^2+xi[1,1]^2*xi[3,2]+xi[2,1]*xi[2,2]^2],
 [h[0]^4*h[3],xi[1,1]^2*xi[1,2]^5*xi[1,3]^3*xi[1,4]^2+xi[1,1]^2*xi[1,2]^3*xi[1,3]^2*xi[1,4]^5+xi[1,1]^2*xi[1,2]^2*xi[1,3]^5*xi[1,4]^3+xi[1,1]^2*xi[1,2]^2*xi[1,3]^2*xi[1,4]^6+xi[1,1]^6*xi[1,2]*xi[1,3]*xi[1,4]*xi[2,1]+xi[1,1]^4*xi[1,2]^3*xi[1,3]^2*xi[2,4]+xi[1,1]^4*xi[1,2]^3*xi[1,3]*xi[1,4]*xi[2,2]+xi[1,1]^4*xi[1,2]^3*xi[1,3]*xi[1,4]*xi[2,3]+xi[1,1]^2*xi[1,2]^5*xi[1,3]*xi[1,4]*xi[2,2]+xi[1,1]^2*xi[1,2]^4*xi[1,3]^2*xi[1,4]*xi[2,3]+xi[1,1]^2*xi[1,2]^4*xi[1,3]*xi[1,4]^2*xi[2,3]+xi[1,1]^2*xi[1,2]^3*xi[1,3]^3*xi[1,4]*xi[2,3]+xi[1,1]^2*xi[1,2]^3*xi[1,3]^3*xi[1,4]*xi[2,4]+xi[1,1]^2*xi[1,2]^3*xi[1,3]^2*xi[1,4]^2*xi[2,2]+xi[1,1]^2*xi[1,2]^3*xi[1,4]^4*xi[2,3]+xi[1,1]^2*xi[1,2]^2*xi[1,3]^4*xi[1,4]*xi[2,3]+xi[1,1]^2*xi[1,2]^2*xi[1,3]^4*xi[1,4]*xi[2,4]+xi[1,1]^2*xi[1,2]^2*xi[1,3]^3*xi[1,4]^2*xi[2,2]+xi[1,1]^2*xi[1,2]^2*xi[1,3]^3*xi[1,4]^2*xi[2,3]+xi[1,1]^2*xi[1,2]^2*xi[1,3]^2*xi[1,4]^3*xi[2,3]+xi[1,1]^2*xi[1,2]^2*xi[1,3]^2*xi[1,4]^3*xi[2,4]+xi[1,1]^2*xi[1,3]^4*xi[1,4]^3*xi[2,2]+xi[1,1]^2*xi[1,3]^2*xi[1,4]^5*xi[2,2]+xi[1,2]^4*xi[1,3]^3*xi[1,4]^2*xi[2,1]+xi[1,2]^2*xi[1,3]^2*xi[1,4]^5*xi[2,1]+xi[1,1]^4*xi[1,2]*xi[1,4]*xi[2,2]*xi[2,3]+xi[1,1]^4*xi[1,3]*xi[1,4]*xi[2,2]^2+xi[1,1]^2*xi[1,2]^3*xi[1,3]*xi[2,3]*xi[2,4]+xi[1,1]^2*xi[1,2]^3*xi[1,4]*xi[2,3]^2+xi[1,1]^2*xi[1,2]^3*xi[1,4]*xi[2,3]*xi[2,4]+xi[1,1]^2*xi[1,2]^2*xi[1,3]^2*xi[2,1]*xi[2,4]+xi[1,1]^2*xi[1,2]^2*xi[1,3]^2*xi[2,3]*xi[2,4]+xi[1,1]^2*xi[1,2]^2*xi[1,3]^2*xi[2,4]^2+xi[1,1]^2*xi[1,2]^2*xi[1,3]*xi[1,4]*xi[2,1]*xi[2,2]+xi[1,1]^2*xi[1,2]^2*xi[1,3]*xi[1,4]*xi[2,1]*xi[2,3]+xi[1,1]^2*xi[1,2]^2*xi[1,3]*xi[1,4]*xi[2,2]^2+xi[1,1]^2*xi[1,2]^2*xi[1,3]*xi[1,4]*xi[2,3]^2+xi[1,1]^2*xi[1,2]^2*xi[1,4]^2*xi[2,2]*xi[2,3]+xi[1,1]^2*xi[1,2]^2*xi[1,4]^2*xi[2,3]^2+xi[1,1]^2*xi[1,2]^2*xi[1,4]^2*xi[2,3]*xi[2,4]+xi[1,1]^2*xi[1,2]*xi[1,3]^2*xi[1,4]*xi[2,2]*xi[2,3]+xi[1,1]^2*xi[1,2]*xi[1,3]^2*xi[1,4]*xi[2,2]*xi[2,4]+xi[1,1]^2*xi[1,2]*xi[1,3]*xi[1,4]^2*xi[2,2]^2+xi[1,1]^2*xi[1,3]^3*xi[1,4]*xi[2,2]^2+xi[1,1]^2*xi[1,3]^3*xi[1,4]*xi[2,2]*xi[2,3]+xi[1,1]^2*xi[1,3]^3*xi[1,4]*xi[2,2]*xi[2,4]+xi[1,1]^2*xi[1,3]^2*xi[1,4]^2*xi[2,2]*xi[2,3]+xi[1,1]^2*xi[1,4]^4*xi[2,2]*xi[2,3]+xi[1,2]^4*xi[1,3]*xi[1,4]*xi[2,1]^2+xi[1,2]^4*xi[1,3]*xi[1,4]*xi[2,1]*xi[2,2]+xi[1,2]^2*xi[1,3]^3*xi[1,4]*xi[2,1]*xi[2,3]+xi[1,2]^2*xi[1,3]^3*xi[1,4]*xi[2,1]*xi[2,4]+xi[1,2]^2*xi[1,3]^2*xi[1,4]^2*xi[2,1]*xi[2,2]+xi[1,2]^2*xi[1,4]^4*xi[2,1]*xi[2,3]+xi[1,1]^2*xi[1,2]*xi[1,3]*xi[1,4]*xi[3,1]+xi[1,1]^2*xi[1,2]*xi[1,3]*xi[1,4]*xi[3,2]+xi[1,1]^2*xi[1,2]*xi[2,2]*xi[2,3]*xi[2,4]+xi[1,1]^2*xi[1,3]^2*xi[1,4]*xi[3,2]+xi[1,1]^2*xi[1,3]*xi[1,4]^2*xi[3,2]+xi[1,1]^2*xi[1,3]*xi[2,2]*xi[2,3]*xi[2,4]+xi[1,1]^2*xi[1,4]*xi[2,1]*xi[2,2]*xi[2,3]+xi[1,1]^2*xi[1,4]*xi[2,2]^2*xi[2,3]+xi[1,1]^2*xi[1,4]*xi[2,2]*xi[2,3]*xi[2,4]+xi[1,2]^3*xi[1,3]*xi[1,4]*xi[3,1]+xi[1,2]^2*xi[1,3]*xi[2,1]*xi[2,3]*xi[2,4]+xi[1,2]^2*xi[1,4]*xi[2,1]^2*xi[2,3]+xi[1,2]^2*xi[1,4]*xi[2,1]*xi[2,3]^2+xi[1,2]^2*xi[1,4]*xi[2,1]*xi[2,3]*xi[2,4]+xi[1,2]*xi[1,3]^2*xi[2,1]^2*xi[2,4]+xi[1,2]*xi[1,3]*xi[1,4]*xi[2,1]^3+xi[1,2]*xi[1,3]*xi[1,4]*xi[2,1]^2*xi[2,2]+xi[1,2]*xi[1,3]*xi[1,4]*xi[2,1]^2*xi[2,3]+xi[1,3]^2*xi[1,4]*xi[2,1]*xi[2,2]*xi[2,3]+xi[1,3]^2*xi[1,4]*xi[2,1]*xi[2,2]*xi[2,4]+xi[1,3]*xi[1,4]^2*xi[2,1]*xi[2,2]^2+xi[1,2]*xi[1,4]*xi[2,3]*xi[3,1]+xi[1,3]*xi[1,4]*xi[2,1]*xi[3,2]+xi[2,1]*xi[2,2]*xi[2,3]*xi[2,4]],
[h[1]^2*h[3]+h[2]^3,xi[1,1]^4*xi[1,2]^2*xi[2,1]^2+xi[1,1]^4*xi[1,2]^2*xi[2,2]^2+xi[2,1]^2*xi[2,2]^2],
[h[0]*c[0],xi[1,1]^2*xi[1,2]^5*xi[1,3]^2*xi[2,3]+xi[1,1]^2*xi[1,2]^4*xi[1,3]^3*xi[2,3]+xi[1,1]^2*xi[1,2]^4*xi[2,3]^2+xi[1,1]^2*xi[1,2]^3*xi[1,3]*xi[2,2]^2+xi[1,1]^2*xi[1,2]*xi[1,3]^3*xi[2,2]^2+xi[1,1]^2*xi[1,3]^4*xi[2,2]^2+xi[1,2]^4*xi[1,3]^2*xi[2,1]*xi[2,3]+xi[1,1]^2*xi[1,2]*xi[1,3]^2*xi[3,2]+xi[1,1]^2*xi[1,2]*xi[2,2]^2*xi[2,3]+xi[1,1]^2*xi[1,3]^3*xi[3,2]+xi[1,1]^2*xi[1,3]*xi[2,2]^3+xi[1,1]^2*xi[1,3]*xi[2,2]^2*xi[2,3]+xi[1,3]^3*xi[2,1]*xi[2,2]^2+xi[1,3]^2*xi[2,1]*xi[3,2]+xi[2,1]*xi[2,2]^2*xi[2,3]]
]:

`is_admissible/Steenrod2` := proc(u)
 local v;

 if type(u,`+`) or type(u,list) or type(u,set) then
  return `and`(op(map(`is_admissible/Steenrod2`,[op(u)])));
 elif type(u,`*`) then
  v := select(type,[op(u)],specfunc(nonnegint,Sq));
  if nops(v) = 0 then
   return true;
  elif nops(v) = 1 then
   return `is_admissible/Steenrod2`(v[1]);
  else
   return FAIL;
  fi;
 elif type(u,specfunc(nonnegint,Sq)) then
  return `and`(true,seq(evalb(op(j,u) >= 2*op(j+1,u)),j=1..nops(u)-1));
 else
  return true;
 fi;
end:

adem_relation := proc(k::nonnegint,j::nonnegint)
 Sq(k,j) + add(modp(binomial(j-m-1,k-2*m),2) * Sq(j+k-m,m),m=0..floor(k/2));
end:

`reduce_Sq` := proc()
 option remember;
 local a,r,m,n,u,v,j,k;

 a := select(i -> i > 0,[args]);
 r := nops(a);

 if `is_admissible/Steenrod2`(Sq(op(a))) then
  return Sq(op(a));
 fi;

 n := 1;
 while n < r and a[n] >= 2*a[n+1] do
  n := n+1;
 od:

 u := seq(a[i],i=1..n-1);
 k := a[n];
 j := a[n+1];
 v := seq(a[i],i=n+2..r);
 return
  modp(add(modp(binomial(j-m-1,k-2*m),2) * reduce_Sq(u,j+k-m,m,v),m=0..floor(k/2)),2);
end:

`reduce/Steenrod2` := (u) -> modp(expand(eval(subs(Sq=reduce_Sq,u))),2);

`mu0/Steenrod2` := proc(u,v)
 if type(u,specfunc(nonnegint,Sq)) and type(v,specfunc(nonnegint,Sq)) then
  return reduce_Sq(op(u),op(v));
 else
  return FAIL;
 fi;
end:

`mu/Steenrod2` := apply_linear_assoc_mod(`mu0/Steenrod2`,Sq(),2);