text:
# Kinematik-Berechnung 3. Arm KAS5
# Beschreibung
#
# Modell m5:
# * Kurbel N7 ist starr mit Körper K2 gekoppelt (über die Achse) (wie m3),
# * Zahnradkopplung zwischen Körpern Z2 und Z3 (Unterschied zu m3)
#
# Setze die Zwangsbedingung (Zahnradkopplung) in die bereits berechneten Zwangsbedingungen für das KAS5m3 ein
#
# Berechnung der Kopplung der Zahnräder am Ellenbogen
# Autor
# Moritz Schappler, schappler@irt.uni-hannover.de, 2016-02
# Institut fuer Regelungstechnik, Leibniz Universitaet Hannover
# Initialisierung
restart:
kin_constraints_exist := true: # Für Speicherung
;
with(LinearAlgebra):
with(StringTools):
read "../helper/proc_convert_s_t":
read "../helper/proc_convert_t_s":
read "../helper/proc_MatlabExport":
read "../robot_codegen_constraints/proc_subs_kintmp_exp":
read "../robot_codegen_definitions/robot_env":
read sprintf("../codeexport/%s/tmp/tree_floatb_definitions", robot_name):
# Lese Ergebnisse der Kinematik von KAS5m3 aus.
# Siehe KAS5m3_sym_codegen_kinematic_constraints.mw
# kintmp_qs, kintmp_qt, lpar_qs, lpar_qt, kintmp_subsexp
read  sprintf("../codeexport/KAS5m3/tmp/kinematic_constraints_maple_inert.m", robot_name):
kintmp_qs_KAS5m3 := kintmp_qs:
kintmp_qt_KAS5m3 := kintmp_qt:
lpar_qs_KAS5m3 := lpar_qs:
lpar_qt_KAS5m3 := lpar_qt:
kintmp_subsexp_KAS5m3 := kintmp_subsexp:
# Wende Zwangsbedingungen aus Zahnradkontakt an
# Wandle die verallgemeinerten Koordinaten um. Benutze die Reihenfolge für KAS5m5.
# Funktion zur Ersetzung der Gelenkwinkel (allgemein -> m3)
subs_q_J_m3 := proc(A)
  local A_s, j;
  A_s := A;
  for j from 1 to RowDimension(qJ_s) do
    A_s := subs({qJ_s(j,1) = qJs_KAS5m3(j,1)}, A_s):
  end do:
  return A_s:
end proc:
qJs_KAS5m3:=Matrix(7, 1, [qJs1_m3, qJs2_m3, qJs3_m3, qJs4_m3, qJs5_m3, qJs6_m3, qJs7_m3]):
kintmp_subsexp_KAS5m3 := subs_q_J_m3(kintmp_subsexp_KAS5m3):
kintmp_qs_KAS5m3 := subs_q_J_m3(kintmp_qs_KAS5m3):
lpar_qs_KAS5m3 := subs_q_J_m3(lpar_qs_KAS5m3):
# Funktion zur Ersetzung der verallgemeinerten Koordinaten (m3 -> m5)
subs_q_m3_m5 := proc(A)
  local A_s, j;
  A_s := A;
  for j from 1 to RowDimension(qJ_s) do
    A_s := subs({qJs_KAS5m3(j,1) = qJs_KAS5m3m5subs(j,1)}, A_s):
  end do:
  return A_s:
end proc:
qJs_KAS5m3m5subs:=Matrix(7, 1, [qJs1_m5, qJs2_m5, qJs3_m5, qJs4_m5, qJs4_m5-qJs3_m5+delta18s, qJs5_m5, 0]):
kintmp_subsexp_KAS5m5 := subs_q_m3_m5(kintmp_subsexp_KAS5m3):
kintmp_qs_KAS5m5 := subs_q_m3_m5(kintmp_qs_KAS5m3):
lpar_qs_KAS5m5 := subs_q_m3_m5(lpar_qs_KAS5m3):
# Konvertiere Ergebnisse
# Rückersetzung der markierten Koordinaten gegen die normale Form (qJ_KAS5m5 -> qJ)
subs_q_m5_J := proc(A)
  local A_s, j;
  A_s := A;
  for j from 1 to RowDimension(qJs_KAS5m5) do
    A_s := subs({qJs_KAS5m5(j,1) = qJ_s(j,1)}, A_s):
  end do:
  return A_s:
end proc:
qJs_KAS5m5:=Matrix(5, 1, [qJs1_m5, qJs2_m5, qJs3_m5, qJs4_m5, qJs5_m5]):
kintmp_subsexp := subs_q_m5_J(kintmp_subsexp_KAS5m5):
kintmp_qs := subs_q_m5_J(kintmp_qs_KAS5m5):
lpar_qs := subs_q_m5_J(lpar_qs_KAS5m5):
kintmp_qt := convert_s_t(kintmp_qs):
lpar_qt := convert_s_t(lpar_qs):
# Speichern der Ergebnisse
save kin_constraints_exist, kintmp_qs, kintmp_qt, lpar_qs, lpar_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert", robot_name):
save kin_constraints_exist, kintmp_qs, kintmp_qt, lpar_qs, lpar_qt, kintmp_subsexp, sprintf("../codeexport/%s/tmp/kinematic_constraints_maple_inert.m", robot_name):
# Liste der Parameter speichern
# Liste mit abhängigen konstanten Kinematikparametern erstellen (wichtig für Matlab-Funktionsgenerierung)
read "../helper/proc_list_constant_expressions";
kc_symbols := Matrix(list_constant_expressions( kintmp_subsexp(..,2) )):
save kc_symbols, sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_maple", robot_name):
MatlabExport(Transpose(kc_symbols), sprintf("../codeexport/%s/tmp/kinematic_constraints_symbols_list_matlab.m", robot_name), 2):
printf("Zwangsbedingungen der Parallelstruktur von KAS5m3 nach KAS5m5 angepasst. %s\n", FormatTime("%Y-%m-%d %H:%M:%S")):


text:
with(simplex);

Constraints := [
# conservation of units at each node
e13+e14 = 300, # CHI
e23+e24 = 300, # DC
e13+e23 = 300, # HOU
e14+e24 = 300, # BOS
# maximum flow on individual edges
0 <= e13, e13 <= 200,
0 <= e14, e14 <= 200,
0 <= e23, e23 <= 280,
0 <= e24, e24 <= 350
];

Cost := 7*e13 + 6*e14 + 4*e23 + 6*e24;
# Invoke linear programming to solve problem
minimize (Cost, Constraints, NONNEGATIVE);

ext: mpl


text:
func $fact (
 var %n i32) i32 {
 var %foo <func(i32) i32>
 var %r1 i32
 var %r2 i32
 var %r3 i32
 if (ne i32 i32 (dread i32 %n, constval i32 1)) {
    icallassigned ( dread a64 %foo,
      sub i32 (dread i32 %n, constval i32 1)) {
      dassign %r1 0
    }
    # Not correct, but for parse and dump test.
    icallassigned ( dread a64 %foo,
      sub i32 (dread i32 %n, constval i32 1)) {
      dassign %r1 1
      dassign %r2 2
      dassign %r3 3
    }
    return (dread i32 %n)
  }
}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl

max_stars_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_issues_repo_name: harmonyos-mirror/OpenArkCompiler-test

max_forks_repo_name: harmonyos-mirror/OpenArkCompiler-test


text:
$ifndef _TCON_REFINE_
$define _TCON_REFINE_

$include "Condition.mpl"
$include "Utils.mpl"

# 简化变换方程的约束条件
tconRefine:=proc(s::TeqSol)
    s:-tcons:=[map[2](singleRefine,s:-rsol,s:-tcons[1]),
               map[2](singleRefine,s:-rsol,s:-tcons[2])];
    # s:-tcons:=map(x->epDeal~(x),s:-tcons);
    return s;
end proc:

# 对于 v[k] 删去 a[k]>0 和 a[k]<0 约束
singleRefine:=proc(rsol,con)
    local ind,n,rep,vs;
    n:=numelems(rsol);
    rep:=add(rsol[i]*v[i],i=1..n);
    vs:=indets(rep,name);
    if numelems(vs)<>1 then
        return con;
    else
        vs:=op([1,1],vs);
        return con minus {a[vs]>0,a[vs]<0};
    end if;
end proc:

# 处理每个可能情况的约束条件
epDeal:=proc(s::set)
    local ca,ce,r;
    ce,ca:=selectremove(has,s,epsilon);
    ce:=ceDeal~(ce);
    r:=ce union ca;
    r:=classifySolve(r);
    return r;
end proc:

# 单个关于epsilon的条件处理
ceDeal:=proc(e)
    local rs;
    # 有的时候solve会抽，对于这种情况不做处理
    try
        rs:=conSolve(e);
    catch :
        return e;
    end try;
    if numelems(rs)=0 then
        # 无解则原样返回
        return e;
    elif numelems(rs)=1 then
        # 只有一个解，则不论有多少条件都可以
        return rs[][];
    else
       return e;
    end if;
end proc:

# 求解不等式，
# 在解中删除自由变量
# 删除和epsilon有关的约束
conSolve:=proc(c)
    local r;
    r:=[RealDomain:-solve(c)];
    r:=map(x->remove(t->evalb(lhs(t)=rhs(t)),x),r);
    r:=map(x->remove(has,x,epsilon),r);
    r:=[{r[]}[]];# 去重
    return r;
end proc:

$endif


text:
Mozilla Public License, version 2.0 1. Definitions 1.1. “Contributor” means each individual or legal entity that creates, contributes to the creation of, or owns Covered Software.

1.2. “Contributor Version” means the combination of the Contributions of others (if any) used by a Contributor and that particular Contributor’s Contribution.

1.3. “Contribution” means Covered Software of a particular Contributor.

1.4. “Covered Software” means Source Code Form to which the initial Contributor has attached the notice in Exhibit A, the Executable Form of such Source Code Form, and Modifications of such Source Code Form, in each case including portions thereof.

1.5. “Incompatible With Secondary Licenses” means

that the initial Contributor has attached the notice described in Exhibit B to the Covered Software; or

that the Covered Software was made available under the terms of version 1.1 or earlier of the License, but not also under the terms of a Secondary License.

1.6. “Executable Form” means any form of the work other than Source Code Form.

1.7. “Larger Work” means a work that combines Covered Software with other material, in a separate file or files, that is not Covered Software.

1.8. “License” means this document.

1.9. “Licensable” means having the right to grant, to the maximum extent possible, whether at the time of the initial grant or subsequently, any and all of the rights conveyed by this License.

1.10. “Modifications” means any of the following:

any file in Source Code Form that results from an addition to, deletion from, or modification of the contents of Covered Software; or

any new file in Source Code Form that contains any Covered Software.

1.11. “Patent Claims” of a Contributor means any patent claim(s), including without limitation, method, process, and apparatus claims, in any patent Licensable by such Contributor that would be infringed, but for the grant of the License, by the making, using, selling, offering for sale, having made, import, or transfer of either its Contributions or its Contributor Version.

1.12. “Secondary License” means either the GNU General Public License, Version 2.0, the GNU Lesser General Public License, Version 2.1, the GNU Affero General Public License, Version 3.0, or any later versions of those licenses.

1.13. “Source Code Form” means the form of the work preferred for making modifications.

1.14. “You” (or “Your”) means an individual or a legal entity exercising rights under this License. For legal entities, “You” includes any entity that controls, is controlled by, or is under common control with You. For purposes of this definition, “control” means (a) the power, direct or indirect, to cause the direction or management of such entity, whether by contract or otherwise, or (b) ownership of more than fifty percent (50%) of the outstanding shares or beneficial ownership of such entity.

2. License Grants and Conditions 2.1. Grants Each Contributor hereby grants You a world-wide, royalty-free, non-exclusive license:

under intellectual property rights (other than patent or trademark) Licensable by such Contributor to use, reproduce, make available, modify, display, perform, distribute, and otherwise exploit its Contributions, either on an unmodified basis, with Modifications, or as part of a Larger Work; and

under Patent Claims of such Contributor to make, use, sell, offer for sale, have made, import, and otherwise transfer either its Contributions or its Contributor Version.

2.2. Effective Date The licenses granted in Section 2.1 with respect to any Contribution become effective for each Contribution on the date the Contributor first distributes such Contribution.

2.3. Limitations on Grant Scope The licenses granted in this Section 2 are the only rights granted under this License. No additional rights or licenses will be implied from the distribution or licensing of Covered Software under this License. Notwithstanding Section 2.1(b) above, no patent license is granted by a Contributor:

for any code that a Contributor has removed from Covered Software; or

for infringements caused by: (i) Your and any other third party’s modifications of Covered Software, or (ii) the combination of its Contributions with other software (except as part of its Contributor Version); or

under Patent Claims infringed by Covered Software in the absence of its Contributions.

This License does not grant any rights in the trademarks, service marks, or logos of any Contributor (except as may be necessary to comply with the notice requirements in Section 3.4).

2.4. Subsequent Licenses No Contributor makes additional grants as a result of Your choice to distribute the Covered Software under a subsequent version of this License (see Section 10.2) or under the terms of a Secondary License (if permitted under the terms of Section 3.3).

2.5. Representation Each Contributor represents that the Contributor believes its Contributions are its original creation(s) or it has sufficient rights to grant the rights to its Contributions conveyed by this License.

2.6. Fair Use This License is not intended to limit any rights You have under applicable copyright doctrines of fair use, fair dealing, or other equivalents.

2.7. Conditions Sections 3.1, 3.2, 3.3, and 3.4 are conditions of the licenses granted in Section 2.1.

3. Responsibilities 3.1. Distribution of Source Form All distribution of Covered Software in Source Code Form, including any Modifications that You create or to which You contribute, must be under the terms of this License. You must inform recipients that the Source Code Form of the Covered Software is governed by the terms of this License, and how they can obtain a copy of this License. You may not attempt to alter or restrict the recipients’ rights in the Source Code Form.

3.2. Distribution of Executable Form If You distribute Covered Software in Executable Form then:

such Covered Software must also be made available in Source Code Form, as described in Section 3.1, and You must inform recipients of the Executable Form how they can obtain a copy of such Source Code Form by reasonable means in a timely manner, at a charge no more than the cost of distribution to the recipient; and

You may distribute such Executable Form under the terms of this License, or sublicense it under different terms, provided that the license for the Executable Form does not attempt to limit or alter the recipients’ rights in the Source Code Form under this License.

3.3. Distribution of a Larger Work You may create and distribute a Larger Work under terms of Your choice, provided that You also comply with the requirements of this License for the Covered Software. If the Larger Work is a combination of Covered Software with a work governed by one or more Secondary Licenses, and the Covered Software is not Incompatible With Secondary Licenses, this License permits You to additionally distribute such Covered Software under the terms of such Secondary License(s), so that the recipient of the Larger Work may, at their option, further distribute the Covered Software under the terms of either this License or such Secondary License(s).

3.4. Notices You may not remove or alter the substance of any license notices (including copyright notices, patent notices, disclaimers of warranty, or limitations of liability) contained within the Source Code Form of the Covered Software, except that You may alter any license notices to the extent required to remedy known factual inaccuracies.

3.5. Application of Additional Terms You may choose to offer, and to charge a fee for, warranty, support, indemnity or liability obligations to one or more recipients of Covered Software. However, You may do so only on Your own behalf, and not on behalf of any Contributor. You must make it absolutely clear that any such warranty, support, indemnity, or liability obligation is offered by You alone, and You hereby agree to indemnify every Contributor for any liability incurred by such Contributor as a result of warranty, support, indemnity or liability terms You offer. You may include additional disclaimers of warranty and limitations of liability specific to any jurisdiction.

4. Inability to Comply Due to Statute or Regulation If it is impossible for You to comply with any of the terms of this License with respect to some or all of the Covered Software due to statute, judicial order, or regulation then You must: (a) comply with the terms of this License to the maximum extent possible; and (b) describe the limitations and the code they affect. Such description must be placed in a text file included with all distributions of the Covered Software under this License. Except to the extent prohibited by statute or regulation, such description must be sufficiently detailed for a recipient of ordinary skill to be able to understand it.

5. Termination 5.1. The rights granted under this License will terminate automatically if You fail to comply with any of its terms. However, if You become compliant, then the rights granted under this License from a particular Contributor are reinstated (a) provisionally, unless and until such Contributor explicitly and finally terminates Your grants, and (b) on an ongoing basis, if such Contributor fails to notify You of the non-compliance by some reasonable means prior to 60 days after You have come back into compliance. Moreover, Your grants from a particular Contributor are reinstated on an ongoing basis if such Contributor notifies You of the non-compliance by some reasonable means, this is the first time You have received notice of non-compliance with this License from such Contributor, and You become compliant prior to 30 days after Your receipt of the notice.

5.2. If You initiate litigation against any entity by asserting a patent infringement claim (excluding declaratory judgment actions, counter-claims, and cross-claims) alleging that a Contributor Version directly or indirectly infringes any patent, then the rights granted to You by any and all Contributors for the Covered Software under Section 2.1 of this License shall terminate.

5.3. In the event of termination under Sections 5.1 or 5.2 above, all end user license agreements (excluding distributors and resellers) which have been validly granted by You or Your distributors under this License prior to termination shall survive termination.

6. Disclaimer of Warranty Covered Software is provided under this License on an “as is” basis, without warranty of any kind, either expressed, implied, or statutory, including, without limitation, warranties that the Covered Software is free of defects, merchantable, fit for a particular purpose or non-infringing. The entire risk as to the quality and performance of the Covered Software is with You. Should any Covered Software prove defective in any respect, You (not any Contributor) assume the cost of any necessary servicing, repair, or correction. This disclaimer of warranty constitutes an essential part of this License. No use of any Covered Software is authorized under this License except under this disclaimer.

7. Limitation of Liability Under no circumstances and under no legal theory, whether tort (including negligence), contract, or otherwise, shall any Contributor, or anyone who distributes Covered Software as permitted above, be liable to You for any direct, indirect, special, incidental, or consequential damages of any character including, without limitation, damages for lost profits, loss of goodwill, work stoppage, computer failure or malfunction, or any and all other commercial damages or losses, even if such party shall have been informed of the possibility of such damages. This limitation of liability shall not apply to liability for death or personal injury resulting from such party’s negligence to the extent applicable law prohibits such limitation. Some jurisdictions do not allow the exclusion or limitation of incidental or consequential damages, so this exclusion and limitation may not apply to You.

8. Litigation Any litigation relating to this License may be brought only in the courts of a jurisdiction where the defendant maintains its principal place of business and such litigation shall be governed by laws of that jurisdiction, without reference to its conflict-of-law provisions. Nothing in this Section shall prevent a party’s ability to bring cross-claims or counter-claims.

9. Miscellaneous This License represents the complete agreement concerning the subject matter hereof. If any provision of this License is held to be unenforceable, such provision shall be reformed only to the extent necessary to make it enforceable. Any law or regulation which provides that the language of a contract shall be construed against the drafter shall not be used to construe this License against a Contributor.

10. Versions of the License 10.1. New Versions Mozilla Foundation is the license steward. Except as provided in Section 10.3, no one other than the license steward has the right to modify or publish new versions of this License. Each version will be given a distinguishing version number.

10.2. Effect of New Versions You may distribute the Covered Software under the terms of the version of the License under which You originally received the Covered Software, or under the terms of any subsequent version published by the license steward.

10.3. Modified Versions If you create software not governed by this License, and you want to create a new license for such software, you may create and use a modified version of this License if you rename the license and remove any references to the name of the license steward (except to note that such modified license differs from this License).

10.4. Distributing Source Code Form that is Incompatible With Secondary Licenses If You choose to distribute Source Code Form that is Incompatible With Secondary Licenses under the terms of this version of the License, the notice described in Exhibit B of this License must be attached.

Exhibit A - Source Code Form License Notice 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/.

If it is not possible or desirable to put the notice in a particular file, then You may include the notice in a location (such as a LICENSE file in a relevant directory) where a recipient would be likely to look for such a notice.

You may add additional accurate notices of copyright ownership.

Exhibit B - “Incompatible With Secondary Licenses” Notice This Source Code Form is “Incompatible With Secondary Licenses”, as defined by the Mozilla Public License, v. 2.0.

max_stars_repo_name: marrow/blueprint

max_issues_repo_name: marrow/blueprint

max_forks_repo_name: marrow/blueprint


text:
#######################################################################
# This file is part of the crlibm library, and is distributed under
# the  LGPL.
# To use:
# restart; read "exp-td.mpl";
Digits := 120:

interface(quiet=true):

read "common-procedures.mpl":
read "triple-double.mpl":

mkdir("TEMPPOW"):

# Some constants for special cases tests

two52 := 2^(52):
two53 := 2^(53):
twoM53 := 2^(-53):
twoM54 := 2^(-54):
twoM1021 := 2^(-1021):
two1021 := 2^(1021):
twoM1000 := 2^(-1000):
two1000 := 2^(1000):
two999 := 2^(999):
two11 := 2^(11):
two64 := 2^(64):
twoM64 := 2^(-64):

scale := 2^(12):
rescale := 2^(-12):
shiftConst := 2^(52) + 2^(51):

largest := 2^(1023) * ((2^(53) - 1) / 2^(52)):
smallest := 2^(-1023) * 1 * 2^(-51):

# Logarithm log2_12 for power

L := 7: # number of bits used to address the table

MAXINDEX    := round(2^L * (sqrt(2)-1)):

for i from 0 to MAXINDEX-1 do
    center[i] := 1 + i*2^(-L): # center[i] in [1, 2[
    t :=  evalf(1/center[i]):
    r[i] := round(t * 2**(floor(-log[2](abs(t))) + 23)) / 2**(floor(-log[2](abs(t))) + 23):
    (logih[i], logil[i], logill[i]) := hi_mi_lo(evalf(-log[2](r[i]))):
od:
for i from MAXINDEX to 2^L do
    # y has been divided by two, center[i] in [0.5, 1[
    center[i]:=(1 + i*2^(-L)) / 2:
    t :=  evalf(1/center[i]):
    r[i] := round(t * 2**(floor(-log[2](abs(t))) + 23)) / 2**(floor(-log[2](abs(t))) + 23):
    (logih[i], logil[i], logill[i]) := hi_mi_lo(evalf(-log[2](r[i]))):
od:

#Computation of ZMax.
for i from 0 to MAXINDEX-1 do
    __x := center[i] + 2^(-L-1) :
    zmax[i] := (__x*r[i]-1) :
    __x := center[i] - 2^(-L-1) :
    zmin[i] := (__x*r[i]-1) :
od:
for i from MAXINDEX to 2^L do
    __x := center[i] + 2^(-L-2) :
    zmax[i] := (__x*r[i]-1) :
    __x := center[i] - 2^(-L-2) :
    zmin[i] := (__x*r[i]-1) :
od:

zmaxmax:=0:
zminmin:=0:
for i from 0 to 2^L do
    if zmax[i] > zmaxmax then zmaxmax := zmax[i]: fi:
    if zmin[i] < zminmin then zminmin := zmin[i]: fi:
od:
printf("zminmin = -2^(%2f)   zmaxmax = 2^(%2f)\n", log2(-zminmin), log2(zmaxmax) ):

pLogExact := unapply(-6497523352912697/36028797018963968*X^8+464105351149111/2251799813685248*X^7-4331547231022885/18014398509481984*X^6+324866042375467/1125899906842624*X^5-6497320848556805/18014398509481984*X^4+8460053188225/17592186044416*X^3-58522663504933901518104329421789/81129638414606681695789005144064*X^2+29261331752466950759073917813481/20282409603651670423947251286016*X,X):

pLog := poly_exact2(pLogExact(x),3):

epsLog := numapprox[infnorm](pLog/log[2](1+x)-1,x=zminmin..zmaxmax):

printf("Relative error of pLog w.r.t log2(1+x) is 2^(%f)\n",log[2](abs(epsLog))):

# Exponential exp2_12 for power

for i from 0 to 2^6 - 1 do
	twoPowerIndex1hi[i], twoPowerIndex1lo[i], twoPowerIndex1lolo[i] := hi_mi_lo(evalf(2^(i/(2^6)))):
	twoPowerIndex2hi[i], twoPowerIndex2lo[i], twoPowerIndex2lolo[i] := hi_mi_lo(evalf(2^(i/(2^(12))))):
od:

pExpExact := unapply(2772236920359585/288230376151711744*X^4+1999746284130149/36028797018963968*X^3+8655072057804175/36028797018963968*X^2+6243314768165359/9007199254740992*X+1,X):

pExp := poly_exact(pExpExact(x)):

epsExp := numapprox[infnorm]((pExp/2^x)-1,x=-2^(-12)..2^(-12)):

printf("Relative error of pExp w.r.t 2^x is 2^(%f)\n",log[2](abs(epsExp))):

log2 := nearest(log(2)):

# Exponential exp2_33 for power

# Polynomial for approximating 2^x - 1 in x=-2^(-42)..2^(-42)

pExpXM := unapply(poly_exact2(x * ((6243314768165359 * 2^(-53) + 29397410857115 * 2^(-100)) + (x * 8655072057804175 * 2^(-55))),2),x):

epsExpXM := numapprox[infnorm](pExpXM(x)/(2^x-1)-1,x=-2^(-42)..2^(-42)):

printf("The relative error of pExpXM w.r.t. 2^x - 1 is 2^(%f)\n",log[2](abs(epsExpXM))):

# Polynomial for approximating 2^x - 1 in x=-2^(-95)..2^(-95)

pExpXL := unapply(poly_exact(x * 6243314768165359 * 2^(-53)), x):

epsExpXL := numapprox[infnorm](pExpXL(x)/(2^x-1)-1,x=-2^(-95)..2^(-95)):

printf("The relative error of pExpXL w.r.t. 2^x - 1 is 2^(%f)\n",log[2](abs(epsExpXL))):

# Polynomial for approximating 2^x in x=-2^(-12)..2^(-12)

pExpXH := unapply(poly_exact32(x * (506517869649829535567849302923399275789356375957 * 2^(-159) + (x * (38978979294391673005692521213079 * 2^(-107) + (x * (36024226132016099441525232746301 * 2^(-109) + (x * (3121261346907607936312652866425 * 2^(-108) + (x * (55385433661433492776134419224183 * 2^(-115) + (x * (5682899659966205 * 2^(-65) + (x * (4501812434047971 * 2^(-68) + (x * 6240991224781291 * 2^(-72))))))))))))))),2,4),x):

epsExpXH := numapprox[infnorm](pExpXH(x)/(2^x-1)-1,x=-2^(-12)..2^(-12)):

printf("The relative error of pExpXH w.r.t. 2^x - 1 is 2^(%f)\n",log[2](abs(epsExpXH))):

# Logarithm log2_13 for power

pLog13 := unapply(poly_exact32(x * (65890661388387311068680317907364672336343222485 * 2^(-155) + (x * ((-263562645553549244274721271629458689351564598207 * 2^(-158)) + (x * (39015109003289267678766993386435 * 2^(-106) + (x * ((-29261331752466950759075245039823 * 2^(-106)) + (x * (11704532700986780303630098000775 * 2^(-105) + (x * ((-78030218006578535357533995772145 * 2^(-108)) + (x * (66883044005638744592219245028355 * 2^(-108) + (x * ((-14630665876233475378839909031845 * 2^(-106)) + (x * (721924538728533 * 2^(-52) + (x * ((-324866042427843 * 2^(-51)) + (x * (4725324253587977 * 2^(-55) + (x * ((-8663094456247059 * 2^(-56)) + (x * (7996693900181249 * 2^(-56) + (x * ((-7425880934468497 * 2^(-56)) + (x * 7091529758988931 * 2^(-56))))))))))))))))))))))))))))),3,6),x):

epsLog13 := numapprox[infnorm](pLog13(x)/(log[2](1+x))-1,x=-2^(-8)..2^(-8)):

printf("The relative error of pLog13 w.r.t. log2(1+x) is 2^(%f)\n",log[2](abs(epsLog13))):

# exp2_30bits for exactness test

coeff_0 :=  9.99999999947486895024439945700578391551971435546875000000000000000000000000000000e-01:
coeff_1 :=  6.93147180274189311788290979166049510240554809570312500000000000000000000000000000e-01:
coeff_2 :=  2.40226513201275415632096610352164134383201599121093750000000000000000000000000000e-01:
coeff_3 :=  5.55041194035996443556513213479775004088878631591796875000000000000000000000000000e-02:
coeff_4 :=  9.61801251055323207228564541537707555107772350311279296875000000000000000000000000e-03:
coeff_5 :=  1.33325640280455024258565721595459763193503022193908691406250000000000000000000000e-03:
coeff_6 :=  1.54736006782907911617439000728779774362919852137565612792968750000000000000000000e-04:
coeff_7 :=  1.55294506644329091183710789270122631933190859854221343994140625000000000000000000e-05:

pExp2Exact := unapply(coeff_0 + X * (coeff_1 + X * (coeff_2 + X * (coeff_3 + X * (coeff_4 + X * (coeff_5 + X * (coeff_6 + X * coeff_7)))))),X):

pExp2 := poly_exact(pExp2Exact(x)):

epsExp2 := numapprox[infnorm]((pExp2/2^x)-1,x=-0.5..0.5):

printf("Relative error of pExp2 w.r.t 2^x is 2^(%f)\n",log[2](abs(epsExp2))):

# Overall accuracy estimate

# C'est pifometrique
epsOverall := 2^(-62):

bi := floor(-log[2](abs(epsOverall))):
approxBoundFactor := 2^(-(bi - 54)):

epsOverallAccurate := 2^(-120):

bi2 := floor(-log[2](abs(epsOverallAccurate))) - 1:
approxBoundFactorAccurate := 2^(-(bi2 - 54)):

# Print out of the .h file

filename:="TEMPPOW/pow.h":
fd:=fopen(filename, WRITE, TEXT):

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

fprintf(fd, "\n/*File generated by maple/pow.mpl*/\n"):

fprintf(fd, "\#define APPROXBOUNDFACTOR %1.50e\n", approxBoundFactor):
fprintf(fd, "\#define APPROXBOUNDFACTORACCURATE %1.50e\n", approxBoundFactorAccurate):
fprintf(fd, "\#define TWO52 %1.50e\n", two52):
fprintf(fd, "\#define TWO53 %1.50e\n", two53):
fprintf(fd, "\#define TWO11 %1.50e\n", two11):
fprintf(fd, "\#define TWOM53 %1.50e\n", twoM53):
fprintf(fd, "\#define TWOM54 %1.50e\n", twoM54):
fprintf(fd, "\#define TWOM1021 %1.50e\n", twoM1021):
fprintf(fd, "\#define TWO1021 %1.50e\n", two1021):
fprintf(fd, "\#define TWOM1000 %1.50e\n", twoM1000):
fprintf(fd, "\#define TWO999 %1.50e\n", two999):
fprintf(fd, "\#define TWO1000 %1.50e\n\n", two1000):
fprintf(fd, "\#define TWO64 %1.50e\n", two64):
fprintf(fd, "\#define TWOM64 %1.50e\n", twoM64):
fprintf(fd, "\#define SCALE %1.50e\n", scale):
fprintf(fd, "\#define RESCALE %1.50e\n", rescale):
fprintf(fd, "\#define SHIFTCONSTANT %1.50e\n", shiftConst):
fprintf(fd, "\#define LARGEST %1.50e\n",largest):
fprintf(fd, "\#define SMALLEST %1.50e\n\n",smallest):

(log2coeff1dh,log2coeff1dl) := hi_lo(coeff(pLog,x,1)):
(log2coeff2dh,log2coeff2dl) := hi_lo(coeff(pLog,x,2)):
fprintf(fd, "\#define log2coeff1h %1.50e\n",log2coeff1dh):
fprintf(fd, "\#define log2coeff1l %1.50e\n",log2coeff1dl):
fprintf(fd, "\#define log2coeff2h %1.50e\n",log2coeff2dh):
fprintf(fd, "\#define log2coeff2l %1.50e\n",log2coeff2dl):
for i from 3 to 8 do
	fprintf(fd, "\#define log2coeff%d %1.50e\n",i,coeff(pLog,x,i)):
od:
fprintf(fd,"\n"):

for i from 1 to 4 do
	fprintf(fd, "\#define exp2coeff%d %1.50e\n",i,coeff(pExp,x,i)):
od:
fprintf(fd,"\n"):

for i from 0 to 7 do
	fprintf(fd, "\#define exp2InaccuCoeff%d %1.50e\n",i,coeff(pExp2,x,i)):
od:
fprintf(fd,"\n"):

(exp2XMcoeff1dh, exp2XMcoeff1dl) := hi_lo(coeff(pExpXM(x),x,1)):
exp2XMcoeff2dh := nearest(coeff(pExpXM(x),x,2)):
fprintf(fd, "\#define exp2XMcoeff1h %1.50e\n",exp2XMcoeff1dh):
fprintf(fd, "\#define exp2XMcoeff1l %1.50e\n",exp2XMcoeff1dl):
fprintf(fd, "\#define exp2XMcoeff2h %1.50e\n\n",exp2XMcoeff2dh):

exp2XLcoeff1dh := nearest(coeff(pExpXL(x),x,1)):
fprintf(fd, "\#define exp2XLcoeff1h %1.50e\n\n",exp2XLcoeff1dh):

(exp2XHcoeff1dh,exp2XHcoeff1dm,exp2XHcoeff1dl) := hi_mi_lo(coeff(pExpXH(x),x,1)):
fprintf(fd, "\#define exp2XHcoeff1h %1.50e\n",exp2XHcoeff1dh):
fprintf(fd, "\#define exp2XHcoeff1m %1.50e\n",exp2XHcoeff1dm):
fprintf(fd, "\#define exp2XHcoeff1l %1.50e\n",exp2XHcoeff1dl):

for i from 2 to 5 do
	(exp2XHcoeffidh,exp2XHcoeffidl) := hi_lo(coeff(pExpXH(x),x,i)):
	fprintf(fd, "\#define exp2XHcoeff%dh %1.50e\n",i,exp2XHcoeffidh):
	fprintf(fd, "\#define exp2XHcoeff%dm %1.50e\n",i,exp2XHcoeffidl):
od:

for i from 6 to 8 do
	fprintf(fd, "\#define exp2XHcoeff%dh %1.50e\n",i,nearest(coeff(pExpXH(x),x,i))):
od:

fprintf(fd,"\n"):

(log213coeff1dh,log213coeff1dm,log213coeff1dl) := hi_mi_lo(coeff(pLog13(x),x,1)):
fprintf(fd, "\#define log213coeff1h %1.50e\n",log213coeff1dh):
fprintf(fd, "\#define log213coeff1m %1.50e\n",log213coeff1dm):
fprintf(fd, "\#define log213coeff1l %1.50e\n",log213coeff1dl):

(log213coeff2dh,log213coeff2dm,log213coeff2dl) := hi_mi_lo(coeff(pLog13(x),x,2)):
fprintf(fd, "\#define log213coeff2h %1.50e\n",log213coeff2dh):
fprintf(fd, "\#define log213coeff2m %1.50e\n",log213coeff2dm):
fprintf(fd, "\#define log213coeff2l %1.50e\n",log213coeff2dl):

for i from 3 to 8 do
	(log213coeffidh, log213coeffidl) := hi_lo(coeff(pLog13(x),x,i)):
	fprintf(fd, "\#define log213coeff%dh %1.50e\n",i,log213coeffidh):
	fprintf(fd, "\#define log213coeff%dm %1.50e\n",i,log213coeffidl):
od:

for i from 9 to 15 do
	fprintf(fd, "\#define log213coeff%dh %1.50e\n",i,nearest(coeff(pLog13(x),x,i))):
od:

fprintf(fd, "\n\#define LOG2 %1.50e\n\n", log2):

fprintf(fd, "typedef struct rri_tag {float ri; double logih; double logil; double logill;} rri;  \n"):
fprintf(fd, "static const rri argredtable[%d] = {\n", 2^L):
for i from 0 to 2^L-1 do
      fprintf(fd, "  { \n"):
      fprintf(fd, "    %1.50e,   /* r[%d] */ \n", r[i], i):
      fprintf(fd, "    %1.50e, /* logih[%d] */ \n", logih[i], i):
      fprintf(fd, "    %1.50e, /* logil[%d] */ \n", logil[i], i):
      fprintf(fd, "    %1.50e, /* logill[%d] */ \n", logill[i], i):
      fprintf(fd, "  } "):
      if(i<2^L-1) then  fprintf(fd, ", \n"): fi
od:
fprintf(fd, "}; \n \n"):

fprintf(fd, "typedef struct tPi_t_tag {double hi; double lo; double lolo;} tPi_t;  \n"):
fprintf(fd, "static const tPi_t twoPowerIndex1[%d] = {\n", 2^(6)):
for i from 0 to 2^(6)-1 do
      fprintf(fd, "  { \n"):
      fprintf(fd, "    %1.50e, /* twoPowerIndex1hi[%d] */ \n", twoPowerIndex1hi[i], i):
      fprintf(fd, "    %1.50e, /* twoPowerIndex1lo[%d] */ \n", twoPowerIndex1lo[i], i):
      fprintf(fd, "    %1.50e, /* twoPowerIndex1lolo[%d] */ \n", twoPowerIndex1lolo[i], i):
      fprintf(fd, "  } "):
      if(i<2^(6)-1) then  fprintf(fd, ", \n"): fi
od:
fprintf(fd, "}; \n \n"):
fprintf(fd, "static const tPi_t twoPowerIndex2[%d] = {\n", 2^(6)):
for i from 0 to 2^(6)-1 do
      fprintf(fd, "  { \n"):
      fprintf(fd, "    %1.50e, /* twoPowerIndex2hi[%d] */ \n", twoPowerIndex2hi[i], i):
      fprintf(fd, "    %1.50e, /* twoPowerIndex2lo[%d] */ \n", twoPowerIndex2lo[i], i):
      fprintf(fd, "    %1.50e, /* twoPowerIndex2lolo[%d] */ \n", twoPowerIndex2lolo[i], i):
      fprintf(fd, "  } "):
      if(i<2^(6)-1) then  fprintf(fd, ", \n"): fi
od:
fprintf(fd, "}; \n \n"):

fclose(fd):

printf("--------- DONE -----------\n");

max_stars_repo_name: squarePenguin/parvsl

max_issues_repo_name: squarePenguin/parvsl

max_forks_repo_name: squarePenguin/parvsl


text:
# Extract a domain from an expression
export Extract := module ()
    export ModuleApply := proc(e, kb, $) :: [ HDomain, anything ];
        local b, eb, s, es;
        b, eb := op(Bound(e));
        s, es := op(Shape(eb, kb ));
        [ DOMAIN(b, s), es ];
    end proc;

    # Extract a domain bound from an expression
    # This pops off the integration constructors recursively, keeping
    # track of the bounds in a KB (which literally become the DBound).
    export Bound := module ()
        export ModuleApply := proc(e, $) :: [ DomBound, anything ];
            local arg, vars, kb;
            arg, vars := do_extract(e)[];
            vars := Domain:-Bound:-withVarsIxs(DBound(vars));
            [ vars , arg ];
        end proc;

        local do_extract_arg := proc(kind, arg_)
            local x0, x, vars, arg := arg_, rng;
            x0  := ExtBound[kind]:-ExtractVar(_rest);   # variable name
            rng := ExtBound[kind]:-ExtractRange(_rest); # variable range
            x   := DInto(x0, rng, kind);                # the variable spec
            arg, vars := do_extract(arg)[];
            [ arg, [ op(vars), x ] ];
        end proc;

        local do_extract := proc(arg, $)
            local sub, prod, svars;
            if arg :: `*` then
                sub := map(do_extract, [op(arg)]);
                prod, svars := selectremove(x->op(2,x)=[],sub);
                if nops(svars) = 1 then
                    [ `*`(op([1,1],svars),op(map2(op,1,prod)))
                    , op([1,2], svars) ];
                else
                    [ arg, [] ];
                end if;
            elif Domain:-Has:-Bound(arg) then
                do_extract_arg(op(0,arg), op(arg));
            else
                [ arg, [] ]
            end if;
        end proc;
    end module;

    # Extract a domain shape from an expression
    # This extracts the domain shape from individual multiplicands of
    # expressions, and multiplies the subexpressions back together.
    # essentially this assumes a distributive law (of domain shapes over
    # products)
    export Shape := module ()
        export ModuleApply := proc(e) :: [ anything, anything ];
            local ixs, w, e1;
            ixs := [indices(ExtShape, 'nolist')];
            w, e1 := do_gets(ixs, true, e) [];
            if not ('no_simpl' in {_rest}) then
                w := simpl_shape(w);
            end if;
            [ w, e1 ];
        end proc;

        local do_get := proc(f, f_ty, e, $)
            local sub, inds, rest;
            if e::`*` then
              sub := map(x->do_get(f, f_ty,x), [op(e)]);
              [ `And`(op(map2(op,1,sub))), `*`(op(map2(op,2,sub))) ]
            elif e::`^` then
              inds, rest := do_get(f, f_ty, op(1,e)) [] ;
              [ inds, subsop(1=rest, e) ]
            elif e:: f_ty then
              f(e,z->ModuleApply(z,'no_simpl'))
            else
              [ true, e ]
            end if
        end proc;

        # apply a list of extractors, in order, until all fail to produce
        # any output .
        local do_gets := proc(todo::list, w, e, $)
            local t0, ts, w1, e1;
            if nops(todo) = 0 then
                [ w, e ]
            else
                t0 := op(1, todo);
                ts := op(subsop(1=NULL, todo));
                w1, e1 := do_get(ExtShape[t0]:-MakeCtx
                                ,ExtShape[t0]:-MapleType
                                ,e) [] ;
                ts := `if`(is(w1), [ts], [ts, t0]);
                do_gets( ts, bool_And(w1, w), e1 );
            end if;
        end proc;

        # todo: simplify the shape
        local simpl_shape := proc(e0,$)
            local e := Domain:-simpl_relation(e0);
            e := subsindets(e, specfunc(`Or`) , x->DSum(op(x)));
            e := subsindets(e, specfunc(`And`), x->DConstrain(op(x)));
            e;
        end proc;
    end module;
end module;

max_stars_repo_name: zaxtax/hakaru

max_issues_repo_name: zaxtax/hakaru

max_forks_repo_name: zaxtax/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: gga_exc *)
(* prefix:
  gga_x_fd_lb94_params *params;

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

(* replace: "fd_int0\(" -> "xc_integrate(func0, NULL, 0.0, " *)
(* replace: "fd_int1\(" -> "xc_integrate(func1, NULL, 0.0, " *)

`diff/fd_int0` := proc(g, x) diff(g, x) * fd_f_inter(0, g) end proc:
`diff/fd_int1` := proc(g, x) diff(g, x) * fd_f_inter(1, g) end proc:

fd_beta := params_a_beta:
fd_csi  := 2^(1/3):

fd_f_inter := (n, t) -> -3/4 * fd_beta*fd_csi*log(t)^n /
  (1 + 3*fd_beta*fd_csi*t*log(fd_csi*t + sqrt((fd_csi*t)^2 + 1))):

fd_f0 := s -> 1 - s*(fd_int0(s)*log(s) - fd_int1(s)):
fd_f  := x -> fd_f0(X2S*x):

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


text:
//impulse: /give @p spawn_egg 1 0 {display:{Name:"Light Bridge Spawner"},EntityTag:{id:"minecraft:chicken",CustomName:"MMH_LightBridgeSpawner",Silent:1,NoGravity:1}}
//MMH_loadChunks()
///clone ~1 ~1 ~1 ~3 ~3 ~3 20 1 16

#MMH
repeat process MMH_lightBridgeSpawner {
if: /testfor @e[type=Chicken,name=MMH_LightBridgeSpawner]
then {
  MMH_loadChunks()
  /execute @e[type=Chicken,name=MMH_LightBridgeSpawner] ~ ~ ~ scoreboard players operation Rotation MMH_Rotation = @p MMH_Rotation
  /scoreboard players test Rotation MMH_Rotation 0 0
  conditional: /execute @e[type=Chicken,name=MMH_LightBridgeSpawner] ~ ~ ~ clone 20 1 16 20 3 17 ~-2 ~-1 ~-1 masked
  /scoreboard players test Rotation MMH_Rotation 1 1
  conditional: /execute @e[type=Chicken,name=MMH_LightBridgeSpawner] ~ ~ ~ clone 21 1 16 22 3 16 ~ ~-1 ~-2 masked
  /scoreboard players test Rotation MMH_Rotation 2 2
  conditional: /execute @e[type=Chicken,name=MMH_LightBridgeSpawner] ~ ~ ~ clone 22 1 17 22 3 18 ~2 ~-1 ~ masked
  /scoreboard players test Rotation MMH_Rotation 3 3
  conditional: /execute @e[type=Chicken,name=MMH_LightBridgeSpawner] ~ ~ ~ clone 20 1 18 21 3 18 ~-1 ~-1 ~2 masked
  /tp @e[type=Chicken,name=MMH_LightBridgeSpawner] ~ -100 ~
}
}

//impulse: /summon area_effect_cloud ~2 ~1 ~ {CustomName:"ACV_LightBridge",Rotation:[-90.0f,0.0f],Duration:2147483647}
///execute @e[name=ACV_lightBridges] ~ ~ ~ setblock ~ ~ ~ redstone_block
//impulse: /summon area_effect_cloud ~2 ~-1 ~ {CustomName:"ACV_AntiBridge",Rotation:[-90.0f,0.0f],Duration:2147483647}
///execute @e[name=ACV_antiBridges] ~ ~ ~ setblock ~ ~ ~ redstone_block
//
//impulse: /summon area_effect_cloud ~ ~1 ~2 {CustomName:"ACV_LightBridge",Rotation:[0.0f,0.0f],Duration:2147483647}
///execute @e[name=ACV_lightBridges] ~ ~ ~ setblock ~ ~ ~ redstone_block
//impulse: /summon area_effect_cloud ~ ~-1 ~2 {CustomName:"ACV_AntiBridge",Rotation:[0.0f,0.0f],Duration:2147483647}
///execute @e[name=ACV_antiBridges] ~ ~ ~ setblock ~ ~ ~ redstone_block
//
//impulse: /summon area_effect_cloud ~-2 ~1 ~ {CustomName:"ACV_LightBridge",Rotation:[90.0f,0.0f],Duration:2147483647}
///execute @e[name=ACV_lightBridges] ~ ~ ~ setblock ~ ~ ~ redstone_block
//impulse: /summon area_effect_cloud ~-2 ~-1 ~ {CustomName:"ACV_AntiBridge",Rotation:[90.0f,0.0f],Duration:2147483647}
///execute @e[name=ACV_antiBridges] ~ ~ ~ setblock ~ ~ ~ redstone_block
//
//impulse: /summon area_effect_cloud ~ ~1 ~-2 {CustomName:"ACV_LightBridge",Rotation:[180.0f,0.0f],Duration:2147483647}
///execute @e[name=ACV_lightBridges] ~ ~ ~ setblock ~ ~ ~ redstone_block
//impulse: /summon area_effect_cloud ~ ~-1 ~-2 {CustomName:"ACV_AntiBridge",Rotation:[180.0f,0.0f],Duration:2147483647}
///execute @e[name=ACV_antiBridges] ~ ~ ~ setblock ~ ~ ~ redstone_block


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

`eta/partitions` := (A::set) -> `if`(nops(A)=1,{A},FAIL);

`gamma/partitions` := (A::set,B::set) -> (p) -> proc(pi,omega)
 local F,rho,u,v,b;

 F := fibres(A,B)(p);

 rho := NULL;
 for u in pi do
  if nops(u) = 1 then
   b := op(u);
   rho := rho,op(omega[b]);
  else
   rho := rho,`union`(seq(F[b],b in u));
  fi;
 od;

 rho := {rho};
 return rho;
end;


text:
function mandel_calc(a)
    z = 0
    for i in 1:100
        z = z^2 + a
    end
    return z
end

function mandelbrot()
    "starting mandelbrot"
    for i in 0:30
    	y = 0.066666*i - 1
	s = ""
	for j in 0:78
	    x = 0.0315*j - 2
	    if abs(mandel_calc((x, y))) < 2 then
	        s = s + "*"
	    else
	        s = s + " "
	    end
	end
	s
    end
end

mandelbrot()

max_stars_repo_name: yugyfoog/MPL

max_issues_repo_name: yugyfoog/MPL

max_forks_repo_name: yugyfoog/MPL


text:
 func $selecti32 (
  var %i i32,
   var %j i32,
   var %k i32) i32 {
   return (
     select i32 (
      dread i32 %i,
      dread i32 %j,
      dread i32 %k))}

 func $selecti64 (
  var %i i32,
   var %j i64,
   var %k i64) i64 {
   return (
     select i64 (
      dread i32 %i,
      dread i64 %j,
      dread i64 %k))}

 func $selectf32 (
  var %i i32,
   var %j f32,
   var %k f32) f32 {
   return (
     select f32 (
      dread i32 %i,
      dread f32 %j,
      dread f32 %k))}

 func $selectf64 (
  var %i i32,
   var %j f64,
   var %k f64) f64 {
   return (
     select f64 (
      dread i32 %i,
      dread f64 %j,
      dread f64 %k))}
 # 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: mgga_exc *)

$define xc_dimensions_2d

js17_lambda := 0.74:
js17_beta   := 30.0:

(*
Notes:
  rho -> 2 rho_sigma, from the spin sum rule
  tsu2d_unif_sigma/n_sigma^2 = 4 Pi, tau is missing 1/2 in paper
*)

(* Eq. (18) *)
js17_ff := s ->
  (1 + 90*(2*js17_lambda - 1)^2*s^2 + js17_beta*(2*js17_lambda - 1)^4*s^4)^(1/15):

(* Eq. (23) *)
js17_R := (s, t) ->
  1 + 128/21*(2*js17_lambda - 1)^2*s^2 +
  (3*(js17_lambda^2 - js17_lambda + 1/2)*(t - 4*Pi) - t)/(4*Pi):

js17_f := (x, u, t) -> 1/js17_ff(X2S_2D*x)
  + 2/5 * js17_R(X2S_2D*x, t)/js17_ff(X2S_2D*x)^3:

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


text:
//impulse: /give @p spawn_egg 1 0 {display:{Name:"Cube Spawner"},EntityTag:{id:"minecraft:chicken",CustomName:"MMH_CubeSpawner",Silent:1,NoGravity:1}}
//MMH_loadChunks()
///clone ~1 ~1 ~1 ~3 ~4 ~3 16 1 20

#MMH
repeat process MMH_weightedStorageCubeSpawner {
if: /testfor @e[type=Chicken,name=MMH_CubeSpawner]
then {
  MMH_loadChunks()
  /execute @e[type=Chicken,name=MMH_CubeSpawner] ~ ~ ~ clone 16 1 20 18 4 22 ~-1 ~ ~-1 masked
  /execute @e[type=Chicken,name=MMH_CubeSpawner] ~ ~ ~ fill ~ ~ ~ ~ ~ ~ air 0 replace pumpkin
  /tp @e[type=Chicken,name=MMH_CubeSpawner] ~ -100 ~
}
}

//repeat: /testfor @e[type=armor_stand,name=ACV_Cube,tag=1]
///blockdata ${-1} {SuccessCount:1}
//conditional: /summon armor_stand ~1 ~-2 ~-1 {CustomName:"ACV_Cube",DisabledSlots:2035728,ArmorItems:[{},{},{},{Count:1,Damage:6,id:"stone"}],NoGravity:1b,Invisible:1b,Invulnerable:1b,Tags:["1"]}


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

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

(* type: gga_exc *)
(* prefix:
  gga_c_wi_params *params;

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

f_num := xt -> params_a_a + params_a_b*xt^2*exp(-params_a_k*xt^2):
f_den := (rs, xt) -> params_a_c + rs*(1 + params_a_d*(4*Pi/3)^(1/3)*xt^(7/2)):

f := (rs, zeta, xt, xs0, xs1) -> f_num(xt)/f_den(rs, xt):


text:
if not(type(lib_dir,string)) then
 lib_dir := currentdir();
fi;

currentdir(lib_dir);
currentdir("../lib_checks");
lib_checks_dir := currentdir();
currentdir(lib_dir);
currentdir("../data");
data_dir := currentdir();
currentdir(lib_dir);

lib_read := proc(f)
 global lib_dir;

 if read_quietly <> true then
  printf("Reading %s\n",f);
 fi;
 read(cat(lib_dir,"/",f)):
end:

with(LinearAlgebra):
with(plots):
with(plottools):

lib_read("cartesian_product.mpl"):
lib_read("realisation.mpl"):
lib_read("apply.mpl"):
lib_read("op_props.mpl"):

lib_read("util.mpl"):
lib_read("colours.mpl"):

lib_read("BW1.mpl"):
lib_read("C.mpl"):
lib_read("CP.mpl"):
lib_read("R.mpl"):
lib_read("sphere.mpl"):
lib_read("RZ.mpl"):
lib_read("RZ_set.mpl"):
lib_read("SW.mpl"):
lib_read("SWW.mpl"):
lib_read("W.mpl"):
lib_read("W1.mpl"):
lib_read("Zpl.mpl"):
lib_read("a_infinity.mpl"):
lib_read("abelian_group_homology.mpl"):
lib_read("autorel.mpl"):
lib_read("big_subsets.mpl"):
lib_read("binary_standard_stasheff_trees.mpl"):
lib_read("binary_trees.mpl"):
lib_read("binops.mpl"):
lib_read("chromatic/Lambda2.mpl"):
lib_read("cubes.mpl"):
lib_read("cycord.mpl"):
lib_read("elliptic/edwards.mpl"):
lib_read("epi.mpl"):
lib_read("equiv.mpl"):
lib_read("euler_quartic.mpl"):
lib_read("extended_stasheff_trees.mpl"):
lib_read("free_algebra.mpl"):
lib_read("full_flags.mpl"):
lib_read("full_trees.mpl"):
lib_read("geometry/geometry.mpl"):
lib_read("geometry/boys.mpl"):
lib_read("geometry/boys_data.mpl"):
lib_read("graphs/digraphs.mpl"):
lib_read("graphs/graphs.mpl"):
lib_read("groebner.mpl"):
lib_read("group_homology.mpl"):
lib_read("groups/alternating_five.mpl"):
lib_read("groups/Out_SL2.mpl"):
lib_read("groups/SL2_11.mpl"):
lib_read("groups/holder_graph.mpl"):
lib_read("hasse_diagrams.mpl"):
lib_read("height_functions.mpl"):
lib_read("homology.mpl"):
lib_read("hurewicz.mpl"):
lib_read("interval_partitions.mpl"):
#lib_read("itloc.mpl"):
lib_read("koszul.mpl"):
lib_read("lie_algebra.mpl"):
lib_read("linord.mpl"):
lib_read("maps.mpl"):
lib_read("morava/fgl.mpl"):
lib_read("morse.mpl"):
lib_read("narayana.mpl"):
lib_read("nat_partitions.mpl"):
lib_read("nonempty_subsets.mpl"):
lib_read("numerical_poly.mpl"):
lib_read("one_cubes_prime.mpl"):
lib_read("operads/general.mpl"):
lib_read("operads/BW.mpl"):
lib_read("operads/F.mpl"):
lib_read("operads/F3.mpl"):
lib_read("operads/Phi.mpl"):
lib_read("operads/WK.mpl"):
lib_read("operads/Xi.mpl"):
lib_read("operads/ass.mpl"):
lib_read("operads/cacti/cactus_lobes.mpl"):
lib_read("operads/cacti/cactus_trees.mpl"):
lib_read("operads/cacti/cactus_planar_trees.mpl"):
lib_read("operads/cacti/cacti.mpl"):
lib_read("operads/chains/ACP.mpl"):
lib_read("operads/chains/E.mpl"):
lib_read("operads/chains/FFbar.mpl"):
lib_read("operads/chains/ICP.mpl"):
lib_read("operads/chains/ICP_alt.mpl"):
lib_read("operads/chains/P.mpl"):
lib_read("operads/chains/PQQ.mpl"):
lib_read("operads/chains/PQ_alt.mpl"):
lib_read("operads/chains/SCP.mpl"):
lib_read("operads/chains/SCP2.mpl"):
lib_read("operads/chains/SEM.mpl"):
lib_read("operads/chains/tree_FFbar.mpl"):
lib_read("operads/comm.mpl"):
lib_read("operads/cubes.mpl"):
lib_read("operads/full_trees.mpl"):
lib_read("operads/fulton/D.mpl"):
lib_read("operads/fulton/E.mpl"):
lib_read("operads/fulton/F1bar.mpl"):
lib_read("operads/fulton/Fbar.mpl"):
lib_read("operads/fulton/P2.mpl"):
lib_read("operads/fulton/Q.mpl"):
lib_read("operads/fulton/WFbar.mpl"):
lib_read("operads/fulton/tree_Fbar.mpl"):
lib_read("operads/fulton/tree_Fbar_alt.mpl"):
lib_read("operads/height_functions.mpl"):
lib_read("operads/lie.mpl"):
lib_read("operads/nonempty_subsets.mpl"):
lib_read("operads/one_cubes_prime.mpl"):
lib_read("operads/ord.mpl"):
lib_read("operads/ord_simplex_interior.mpl"):
lib_read("operads/parallel/PP.mpl"):
lib_read("operads/parallel/YC.mpl"):
lib_read("operads/partitions.mpl"):
lib_read("operads/prime_simplex.mpl"):
lib_read("operads/simplex.mpl"):
lib_read("operads/simplex_interior.mpl"):
lib_read("operads/stasheff_star.mpl"):
lib_read("operads/stasheff_trees.mpl"):
lib_read("operads/trees.mpl"):
lib_read("ord.mpl"):  # Duplicates linord.mpl
lib_read("ord_simplex_interior.mpl"):
lib_read("pair_subsets.mpl"):
lib_read("pairs.mpl"):
lib_read("partitions.mpl"):
lib_read("partord.mpl"):
lib_read("permutations.mpl"):
lib_read("plane_cad.mpl"):
lib_read("plots.mpl"):
lib_read("poisson.mpl"):
lib_read("posint_intervals.mpl"):
lib_read("preord.mpl"):
lib_read("prime_W.mpl"):
lib_read("prime_simplex.mpl"):
lib_read("prime_simplex_boundary.mpl"):
lib_read("proper_nonempty_subsets.mpl"):
lib_read("quaternions.mpl"):
lib_read("real_functions.mpl"):
lib_read("rel.mpl"):
lib_read("semioperads/general.mpl"):
lib_read("semioperads/circ.mpl"):
lib_read("semioperads/coloured.mpl"):
lib_read("semioperads/permutations.mpl"):
lib_read("shelling.mpl"):
lib_read("shuffles.mpl"):
lib_read("simplex.mpl"):
lib_read("simplex_interior.mpl"):
lib_read("simplices.mpl"):
lib_read("simplicial.mpl"):
lib_read("simplicial_alt.mpl"):
lib_read("simplicial_complexes/simplicial_complex.mpl"):
lib_read("simplicial_complexes/CP2_complex.mpl"):
lib_read("simplicial_complexes/RP2_complex.mpl"):
lib_read("simplicial_complexes/icosahedron_complex.mpl"):
lib_read("simplicial_complexes/octahedron_complex.mpl"):
lib_read("simplicial_complexes/mobius_complex.mpl"):
lib_read("simplicial_complexes/moore_complex.mpl"):
lib_read("single_cubes.mpl"):
lib_read("space_curves.mpl"):
lib_read("spaces/spaces.mpl"):
lib_read("spaces/complex_projective_space.mpl"):
lib_read("standard_stasheff_star.mpl"):
lib_read("standard_stasheff_trees.mpl"):
lib_read("stasheff_star.mpl"):
lib_read("stasheff_trees.mpl"):
lib_read("steiner_systems.mpl"):
lib_read("striped_preord.mpl"):
lib_read("subsets.mpl"):
lib_read("subspaces.mpl"):
lib_read("symfun.mpl"):
lib_read("tanabe/tanabe.mpl"):
lib_read("tanabe/twisted_product.mpl"):
lib_read("total_preord.mpl"):
lib_read("tree_height_functions.mpl"):
lib_read("trees.mpl"):
lib_read("triples.mpl"):
lib_read("vector_functions.mpl"):
lib_read("witt.mpl"):
lib_read("young.mpl"):


text:
# A Groebner basis is by definition a set of generators for an ideal
# I in a polynomial ring over a field k.  Given such a set, one can
# deduce a basis for the quotient ring P/I as a vector space over k.
# We call that a cobasis for I.  The function below calculates the
# cobasis, assuming that it is finite.  (It will just give an infinite
# loop if the cobasis is infinite.)

cobasis := proc(basis,vars)
 local lm,nf,L,L0,x,m,i;

 lm := map(LeadingMonomial,basis,vars);
 nf := (u) -> NormalForm(u,lm,vars);
 L := [1];
 for x in [op(vars)] do
  L0 := L;
  L := NULL;
  for m in L0 do
   i := 0;
   while nf(m*x^i) <> 0 do
    L := L,m*x^i;
    i := i+1;
   od;
  od:
  L := [L];
 od:

 return L;
end:


text:
make_upper_unitriangular := proc(n,p)
  local G,i,j,k,m,g,L,f;

  m := n * (n - 1) / 2;

  G := table():
  G["p"] := p;
  G["id"] := IdentityMatrix(n);
  G["o"] := (g,h) -> map(mods,g . h,G["p"]);
  G["inv"] := g -> map(mods,1/g,G["p"]);
  G["order"] := p ^ m;

  G["elementary"] := table():
  for i from 1 to n do
   for j from i+1 to n do
    g := Matrix(n,n);
    for k from 1 to n do g[k,k] := 1; od:
    g[i,j] := 1;
    G["elementary"][i,j] := g;
   od:
  od:

  L := [[]];
  for i from 1 to m do
   L := [seq(seq([op(u),j],j=0..p-1),u in L)];
  od:

  f := proc(u)
   local i,j,k,g;
   g := Matrix(n,n);
   for k from 1 to n do g[k,k] := 1; od:
   k := 1;
   for i from 1 to n do
    for j from i+1 to n do
     g[i,j] := mods(u[k],p); k := k + 1;
    od:
   od:
   return g;
  end:

  G["elements"] := map(f,L);

  return eval(G);
end:


text:
with(Threads:-Mutex):
with(Threads:-Task):
m:=Create();
c:=1;
fun:=proc()
    global c,m;
    Lock(m);
    print(c);
    c:=c+1;
    Unlock(m);
end proc:
Start(null,Tasks=[fun,seq([],i=1..10)]);
Destroy(m);


text:
# This is the problem 11914 the American Mathematical Monthly (May 2016)
# The following sum is zero.

S := Sum(Sum(x^n*y^m*Sum((-4)^(-k)*Binomial(n-k,k-1)*Sum((-2)^(-j)*Binomial(n+1-2*k,j-1)*Binomial(m-k,3*m-j),j=1..3*m), k=1..n), n=1..infinity), m=1..infinity);

R, ord := BinomSums[sumtores](S, u);
#  -> 0, []

max_stars_repo_name: lairez/binomsum

max_issues_repo_name: lairez/binomsum

max_forks_repo_name: lairez/binomsum


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

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

(* type: gga_exc *)

bayesian_theta0 := 1.0008:
bayesian_theta1 := 0.1926:
bayesian_theta2 := 1.8962:

bayesian_f0 := s -> s^2/(1 + s)^2:
bayesian_f  := x -> bayesian_theta0 + bayesian_f0(X2S*x)* (bayesian_theta1 + bayesian_f0(X2S*x) * bayesian_theta2):

f := (rs, zeta, xt, xs0, xs1) -> gga_exchange(bayesian_f, rs, zeta, xs0, xs1):


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

# make gensym global, so that it can be shared with other 'global' routines
gensym := module ()
  export ModuleApply;
  local gs_counter, utf8, blocks, radix, unicode;
  gs_counter := -1;
  utf8 := proc(n :: integer, $)
    local m;
    if n<128 then n
    elif n<2048 then 192+iquo(n,64,'m'), 128+m
    elif n<65536 then 224+iquo(n,4096,'m'), 128+iquo(m,64,'m'), 128+m
    elif n<2097152 then 240+iquo(n,262144,'m'), 128+iquo(m,4096,'m'), 128+iquo(m,64,'m'), 128+m
    elif n<67108864 then 248+iquo(n,16777216,'m'), 128+iquo(m,262144,'m'), 128+iquo(m,4096,'m'), 128+iquo(m,64,'m'), 128+m
    elif n<2147483648 then 248+iquo(n,1073741824,'m'), 128+iquo(m,16777216,'m'), 128+iquo(m,262144,'m'), 128+iquo(m,4096,'m'), 128+iquo(m,64,'m'), 128+m
    end if
  end proc;
  blocks := map((b -> block(convert(op(0,b), decimal, hex), op(1,b))),
                ["4e00"(20950)]);
  radix := `+`(op(map2(op, 2, blocks))) / 2;
  unicode := proc(nn, $)
    local n, b;
    n := nn;
    for b in blocks do
      if n < op(2,b) then return n + op(1,b) else n := n - op(2,b) end if
    end do
  end proc;
  ModuleApply := proc(x::name, $)
    gs_counter := gs_counter + 1;
    cat(x, op(map(StringTools:-Char, map(utf8 @ unicode, applyop(`+`, 1, map(`*`, convert(gs_counter, 'base', radix), 2), 1)))))
  end proc;
end module: # gensym

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

BindingTools := module ()
  option package;
  export generic_evalat, generic_evalatstar;

  generic_evalat := proc(vv::{name,list(name)}, mm, eqs, $)
    local v, m, eqsRemain, eq, rename, funs;
    funs := map2(op, 0, indets(mm, 'function'));
    eqsRemain := select(
      (eq -> op(1,eq) <> op(2,eq)
             and not has(op(1,eq), vv)
             and (not (op(1,eq) :: 'name')
                  or depends(mm, op(1,eq))
                  or member(op(1,eq), funs))),
      eqs);
    if nops(eqsRemain) = 0 then return vv, mm end if;
    m := mm;
    rename := proc(v::name, $)
      local vRename;
      if depends(eqsRemain, v) then
        vRename := gensym(v);
        m := subs(v=vRename, m);
        vRename
      else
        v
      end if
    end proc;
    if vv :: name then
      v := rename(vv)
    else
      v := map(rename, vv);
    end if;
    v, eval(m, eqsRemain);
  end proc:

  generic_evalatstar := proc(body, bound::list, eqs, $)
    local indefinite, e, n, b, j;
    e := foldl(((x,i) -> `if`(i::`=`,
                              lam(lhs(i), rhs(i), x),
                              lam(i, indefinite, x))),
               body, op(bound));
    e := eval(e, eqs); # Piggyback on `eval/lam`
    n := nops(bound);
    b := table();
    for j from n to 1 by -1 do
      b[j] := `if`(op(2,e)=indefinite, op(1,e), `=`(op(1..2,e)));
      e := op(3,e);
    end do;
    e, [entries(b, 'nolist', 'indexorder')]
  end proc:
end module: # BindingTools


text:
######################################################################
# Commutative operad

`is_element/comm` := (A::set) -> (z) -> type(z,identical(0));

`is_equal/comm` := (A::set) -> (x,y) -> true;

`is_leq/comm` := (A::set) -> (x,y) -> true;

`random_element/comm` := (A::set) -> () -> 0;

`list_elements/comm` := (A::set) -> [0];

`count_elements/comm` := (A::set) -> 1;

`eta/comm` := (A::set) -> `if`(nops(A)=1,0,FAIL);

`gamma/comm` := (A::set,B::set) -> (p) -> (z,zz) -> 0;


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

c     :=  0.7168:
alpha :=  2.804:
d     := 28.23705740248932030511071641312341561894: (* POW(CBRT(4/3) * 2*M_PI/3, 4) *)

csi  := s -> (3/2 * LambertW(s^(3/2) / (2*sqrt(6))))^(2/3):
fb   := s -> Pi/3 * s/(csi(s) * (d + csi(s)^2)^(1/4)):
flaa := s -> (1 + c*s^2)/(1 + c*s^2/fb(s)):
XX   := s -> 1 - alpha*s^2/(1 + alpha*s^2):

f    := x ->  XX(X2S*x) + (1 - XX(X2S*x))*flaa(X2S*x):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
interface(screenwidth=9999):
if not (NewSLO :: `module`) then
  WARNING("loading NewSLO failed");
  `quit`(3);
end if;

with(Hakaru):
with(NewSLO):

#####################################################################
#
# Dirichlet conjugacy tests
#
#####################################################################

(dice_index, dice_index_t) := Concrete("examples/dice_index.hk"):
dice_index := value(eval(dice_index)):

(gmm_gibbs, gmm_gibbs_t) := Concrete("examples/gmm_gibbs.hk"):
gmm_gibbs := value(eval(gmm_gibbs)):

(naive_bayes_gibbs, naive_bayes_gibbs_t) := Concrete("examples/naive_bayes/naive_bayes_gibbs.hk"):
naive_bayes_gibbs := value(eval(naive_bayes_gibbs)):

#####################################################################
TestEfficient( dice_index, dice_index_t, KB:-empty, label="dice index" ):
TestEfficient( gmm_gibbs, gmm_gibbs_t, KB:-empty, label="gmm gibbs" ):
TestEfficient( naive_bayes_gibbs, naive_bayes_gibbs_t, KB:-empty, label="naive bayes gibbs" ):


text:
#######################################################################
# This file is part of the crlibm library, and is distributed under
# the  LGPL.
# To use:
# restart; read "exp-td.mpl";
Digits := 120:

interface(quiet=true):

read "common-procedures.mpl":
read "triple-double.mpl":
mkdir("TEMPEXP"):

L := 12:

printf("   memory requirement for L = %d and two triple-double tables: %d bytes\n",L,48*2^(ceil(L/2)));

rmax := log(2) / (2^(L+1)):

printf("   maximal absolute value for rmax = 2^(%f)\n",log[2](rmax)):

MsLog2Div2L := evalf(-log(2)/(2^L)):

msLog2Div2Lh, msLog2Div2Lm, msLog2Div2Ll := hi_mi_lo(MsLog2Div2L):

epsMsLog2Div2L := evalf(abs(((msLog2Div2Lh + msLog2Div2Lm + msLog2Div2Ll) - MsLog2Div2L)/MsLog2Div2L)):
epsDDMsLog2Div2L := evalf(abs(((msLog2Div2Lh + msLog2Div2Lm) - MsLog2Div2L)/MsLog2Div2L)):

printf("   error made by storing MsLog2Div2L as a double-double: 2^(%f)\n",log[2](epsDDMsLog2Div2L)):
printf("   error made by storing MsLog2Div2L as a triple-double: 2^(%f)\n",log[2](epsMsLog2Div2L)):

gap := -floor(-log[2](abs(msLog2Div2Lm/msLog2Div2Lh))):

printf("   |msLog2Div2Lm| <= 2^(%f) * |msLog2Div2Lh|\n",gap):

log2InvMult2L := nearest(2^L / (log(2))):

shiftConst := 2^(52) + 2^(51):

indexmask1 := 2^(L/2) - 1:
indexmask2 := indexmask1 * 2^(L/2):

largest := 2^(1023) * ((2^(53) - 1) / 2^(52)):
smallest := 2^(-1023) * 1 * 2^(-51):

overflowbound := nearest(log(largest)):

overflowboundHex := ieeehexa(overflowbound):
overflowSimplebound := convert(overflowboundHex[1],decimal,hex):

underflowbound := nearest(log(2^(-1075))):

denormbound := nearest(log(2^(-1022) * 1)):

overUnderflowboundHex := ieeehexa(min(abs(underflowbound),min(abs(overflowbound),abs(denormbound)))):
overUnderflowSimplebound := convert(overUnderflowboundHex[1],decimal,hex):

twoPowerM1000 := 2^(-1000):
twoPower1000 := 2^(1000):

twoM52 := 2^(-52):
mTwoM53 := - 2^(-53):

for i from 0 to 2^(L/2) - 1 do
	twoPowerIndex1hi[i], twoPowerIndex1mi[i], twoPowerIndex1lo[i] := hi_mi_lo(evalf(2^(i/(2^L)))):
	twoPowerIndex2hi[i], twoPowerIndex2mi[i], twoPowerIndex2lo[i] := hi_mi_lo(evalf(2^(i/(2^(L/2))))):
od:

PolyDegreeQuick:=4:
printf("   degree of the polynomial used in the quick phase is %d\n",PolyDegreeQuick);

polyQuick:= poly_exact(1 + x + 0.5*x^2 + x^3 * (numapprox[minimax](((exp(x) - (1 + x + 0.5*x^2))/x^3),
				x=-rmax..rmax, [PolyDegreeQuick-3,0], 1 ,  'deltaApprox'))):

epsilonApproxQuick := numapprox[infnorm]( ((polyQuick-1)/(exp(x)-1))-1, x=-rmax..rmax):
printf("   approximation rel error for the quick phase is 2^(%2f)\n", log2(epsilonApproxQuick) ) :
deltaApproxQuick := numapprox[infnorm]( polyQuick-exp(x), x=-rmax..rmax):
printf("   approximation abs error for the quick phase is 2^(%2f)\n", log2(deltaApproxQuick) ) :

PolyDegreeAccurate:=7:

printf("   degree of the polynomial used in the accurate phase is %d\n",PolyDegreeAccurate):

DDNumberAccu:=5:

printf("   number of double doubles used for the coefficients is %d\n",DDNumberAccu):

polyAccurate:= poly_exact2(1 + x + 0.5*x^2 + x^3 * (numapprox[minimax](((exp(x) - (1 + x + 0.5*x^2))/x^3),
				x=-rmax..rmax, [PolyDegreeAccurate-3,0], 1 ,  'deltaApprox')),
				DDNumberAccu):

epsilonApproxAccurate := numapprox[infnorm]( ((polyAccurate-1)/(exp(x)-1))-1, x=-rmax..rmax):
printf("   approximation rel error for the accurate phase is 2^(%2f)\n", log2(epsilonApproxAccurate) ) :
deltaApproxAccurate := numapprox[infnorm]( polyAccurate-exp(x), x=-rmax..rmax):
printf("   approximation abs error for the accurate phase is 2^(%2f)\n", log2(deltaApproxAccurate) ) :

epsilonApproxRmAccurate := numapprox[infnorm]( (x/(exp(x)-1))-1, x=-rmax*2^(-52)..rmax*2^(-52)):
epsilonApproxRlAccurate := numapprox[infnorm]( (x/(exp(x)-1))-1, x=-rmax*2^(-105)..rmax*2^(-105)):

printf("   approximation rel error for approximating exp(rm) - 1 by rm is 2^(%2f)\n", log2(abs(epsilonApproxRmAccurate))):
printf("   approximation rel error for approximating exp(rl) - 1 by rl is 2^(%2f)\n", log2(abs(epsilonApproxRlAccurate))):

epsilonApproxAccurateSpecial := numapprox[infnorm]( ((polyAccurate-1)/(exp(x)-1))-1, x=-2^(-30)..2^(-30)):
printf("   approximation rel error for the accurate phase in the special interval (|r| \\leq 2^(-30)) is 2^(%2f)\n",
log2(epsilonApproxAccurateSpecial) ) :

epsilonApproxAccurateSpecial2 := numapprox[infnorm]( ((polyAccurate-1)/(exp(x)-1))-1, x=-2^(-18)..2^(-18)):
printf("   approximation rel error for the accurate phase in the special interval (|r| \\leq 2^(-18)) is 2^(%2f)\n",
log2(epsilonApproxAccurateSpecial2) ) :

epsilon_quick := 2^(-64): # The Gappa proof will show this bound

#-------------------------------------------------------------------
# Output

filename:="TEMPEXP/exp-td.h":
fd:=fopen(filename, WRITE, TEXT):

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

fprintf(fd, "\n/*File generated by maple/exp-td.mpl*/\n"):

fprintf(fd, "\#define L %d\n",L):
fprintf(fd, "\#define LHALF %d\n",L/2):
fprintf(fd, "\#define log2InvMult2L %1.50e\n",log2InvMult2L):
fprintf(fd, "\#define msLog2Div2Lh %1.50e\n",msLog2Div2Lh):
fprintf(fd, "\#define msLog2Div2Lm %1.50e\n",msLog2Div2Lm):
fprintf(fd, "\#define msLog2Div2Ll %1.50e\n",msLog2Div2Ll):
fprintf(fd, "\#define shiftConst %1.50e\n",shiftConst):
fprintf(fd, "\#define INDEXMASK1 0x%08x\n",indexmask1):
fprintf(fd, "\#define INDEXMASK2 0x%08x\n",indexmask2):
fprintf(fd, "\#define OVRUDRFLWSMPLBOUND 0x%08x\n",overUnderflowSimplebound):
fprintf(fd, "\#define OVRFLWBOUND %1.50e\n",overflowbound):
fprintf(fd, "\#define LARGEST %1.50e\n",largest):
fprintf(fd, "\#define SMALLEST %1.50e\n",smallest):
fprintf(fd, "\#define DENORMBOUND %1.50e\n",denormbound):
fprintf(fd, "\#define UNDERFLWBOUND %1.50e\n",underflowbound):
fprintf(fd, "\#define twoPowerM1000 %1.50e\n",twoPowerM1000):
fprintf(fd, "\#define twoPower1000 %1.50e\n",twoPower1000):
fprintf(fd, "\#define ROUNDCST %1.50e\n", compute_rn_constant(epsilon_quick)):
fprintf(fd, "\#define RDROUNDCST %1.50e\n", epsilon_quick):
fprintf(fd, "\#define twoM52 %1.50e\n", twoM52):
fprintf(fd, "\#define mTwoM53 %1.50e\n", mTwoM53):

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

for i from 3 to PolyDegreeQuick do
	fprintf(fd, "\#define c%d %1.50e\n",i,coeff(polyQuick,x,i)):
od:

for i from 3 to DDNumberAccu-1 do
	(hi,lo) := hi_lo(coeff(polyAccurate,x,i)):
	fprintf(fd, "\#define accPolyC%dh %1.50e\n",i,hi):
	fprintf(fd, "\#define accPolyC%dl %1.50e\n",i,lo):
od:

for i from DDNumberAccu to PolyDegreeAccurate do
	fprintf(fd, "\#define accPolyC%d %1.50e\n",i,coeff(polyAccurate,x,i)):
od:

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

# Print the tables
fprintf(fd, "typedef struct tPi_t_tag {double hi; double mi; double lo;} tPi_t;  \n"):
fprintf(fd, "static const tPi_t twoPowerIndex1[%d] = {\n", 2^(L/2)):
for i from 0 to 2^(L/2)-1 do
      fprintf(fd, "  { \n"):
      fprintf(fd, "    %1.50e, /* twoPowerIndex1hi[%d] */ \n", twoPowerIndex1hi[i], i):
      fprintf(fd, "    %1.50e, /* twoPowerIndex1mi[%d] */ \n", twoPowerIndex1mi[i], i):
      fprintf(fd, "    %1.50e, /* twoPowerIndex1lo[%d] */ \n", twoPowerIndex1lo[i], i):
      fprintf(fd, "  } "):
      if(i<2^(L/2)-1) then  fprintf(fd, ", \n"): fi
od:
fprintf(fd, "}; \n \n"):
fprintf(fd, "static const tPi_t twoPowerIndex2[%d] = {\n", 2^(L/2)):
for i from 0 to 2^(L/2)-1 do
      fprintf(fd, "  { \n"):
      fprintf(fd, "    %1.50e, /* twoPowerIndex2hi[%d] */ \n", twoPowerIndex2hi[i], i):
      fprintf(fd, "    %1.50e, /* twoPowerIndex2mi[%d] */ \n", twoPowerIndex2mi[i], i):
      fprintf(fd, "    %1.50e, /* twoPowerIndex2lo[%d] */ \n", twoPowerIndex2lo[i], i):
      fprintf(fd, "  } "):
      if(i<2^(L/2)-1) then  fprintf(fd, ", \n"): fi
od:
fprintf(fd, "}; \n \n"):

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

fclose(fd):

filename:="TEMPEXP/exp-td-accurate.sed":
fd:=fopen(filename, WRITE, TEXT):
fprintf(fd, " s/_rhmax/%1.50e/g\n", rmax):
fprintf(fd, " s/_rmmax/%1.50e/g\n", rmax*2^(-52)):
fprintf(fd, " s/_rlmax/%1.50e/g\n", rmax*2^(-105)):
fprintf(fd, " s/_epsilonApproxAccurate/%1.50e/g\n", epsilonApproxAccurate):
fprintf(fd, " s/_epsilonApproxRmAccurate/%1.50e/g\n", epsilonApproxRmAccurate):
fprintf(fd, " s/_epsilonApproxRlAccurate/%1.50e/g\n", epsilonApproxRlAccurate):
fprintf(fd, " s/_epsilonApproxSpecial/%1.50e/g\n", epsilonApproxAccurateSpecial):

for i from 3 to DDNumberAccu-1 do
	(hi,lo) := hi_lo(coeff(polyAccurate,x,i)):
	fprintf(fd, "s/_accPolyC%dh/%1.50e/g\n",i,hi):
	fprintf(fd, "s/_accPolyC%dl/%1.50e/g\n",i,lo):
od:

for i from DDNumberAccu to PolyDegreeAccurate do
	fprintf(fd, "s/_accPolyC%d/%1.50e/g\n",i,coeff(polyAccurate,x,i)):
od:

fclose(fd):

printf("----DONE---\n") :


text:
input := ImportMatrix("AoC-2021-1-input.txt" ):

scan :=  (v,d)->ifelse(d[2]<0,[v,0],
                        [v, ifelse(v<=d[1], d[2], d[2]+1)]);

rtable_scanblock(input, [], noindex, scan, [-1,-1] )[2]; # part 1 answer

windows := Array([seq(input[i]+input[i+1]+input[i+2],
   i=1..numelems(input)-2)]);

rtable_scanblock(windows, [], noindex, scan, [-1,-1])[2]; # part 2 answer


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

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

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

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

scan_p     := x -> X2S^2*x^2:
scan_alpha := (x, t) -> (t - x^2/8)/K_FACTOR_C:

(* The interpolating functions are nasty for a -> 1, so we need to
   truncate them. The natural choice is to cut off the functions to
   zero when the exponential term reaches machine epsilon.

   The left cutoff is |log epsilon|/(|log epsilon| + c1) < 1
   and the right one is (|log epsilon| + c2)/|log epsilon| > 1,
   so we don't even really need the step function.
*)
scan_f_alpha_left0 := a -> exp(-params_a_c1*a/(1 - a)):
scan_f_alpha_left_cutoff := -log(DBL_EPSILON)/(-log(DBL_EPSILON) + params_a_c1):
scan_f_alpha_left := a -> my_piecewise3(a > scan_f_alpha_left_cutoff, 0, scan_f_alpha_left0(m_min(scan_f_alpha_left_cutoff, a))):

scan_f_alpha_right0 := a -> -params_a_d*exp(params_a_c2/(1 - a)):
scan_f_alpha_right_cutoff := (-log(DBL_EPSILON/abs(params_a_d)) + params_a_c2)/(-log(DBL_EPSILON/abs(params_a_d))):
scan_f_alpha_right := a -> my_piecewise3(a < scan_f_alpha_right_cutoff, 0, scan_f_alpha_right0(m_max(scan_f_alpha_right_cutoff, a))):
scan_f_alpha := a -> my_piecewise3(
  a <= 1, scan_f_alpha_left(a), scan_f_alpha_right(a)
  ):

scan_h1x := x -> 1 + params_a_k1*(1 - params_a_k1/(params_a_k1 + x)):

scan_b2 := sqrt(5913/405000):
scan_b1 := (511/13500)/(2*scan_b2):
scan_b3 := 1/2:
scan_b4 := MU_GE^2/params_a_k1 - 1606/18225 - scan_b1^2:
scan_y  := (x, a) -> MU_GE*scan_p(x) + scan_b4*scan_p(x)^2*exp(-scan_b4*scan_p(x)/MU_GE)
  + (scan_b1*scan_p(x) + scan_b2*(1 - a)*exp(-scan_b3*(1 - a)^2))^2:

scan_a1 := 4.9479:
scan_gx := x -> 1 - exp(-scan_a1/sqrt(X2S*x)):

scan_h0x := 1.174:
scan_f   := (x, u, t) -> (scan_h1x(scan_y(x, scan_alpha(x, t)))*(1 - scan_f_alpha(scan_alpha(x, t)))
  + scan_h0x*scan_f_alpha(scan_alpha(x, t)))*scan_gx(x):

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


text:
#######################################################################
# This file is part of the crlibm library, and is distributed under
# the  LGPL.
# To use:
# restart; read "exp-td.mpl";
Digits := 145:

interface(quiet=true):

read "common-procedures.mpl":
read "triple-double.mpl":
read "gal.mpl":
mkdir("TEMPASIN"):

with(orthopoly,T):

truncPoly := proc(p,k)
local i, q:
q := 0:
convert(q,polynom):
for i from 0 to min(degree(p,x),k) do
    q := q + coeff(p,x,i) * x^i:
od:
return (q):
end:

intervals := 10:

for i from 1 to intervals do
	epsAccurate[i] := 0.5:
	epsQuick[i] := 0.5:
	polyAccurate[i] := 1:
	polyQuick[i] := 1:
od:
epsAccurateSpecial := 0.5:
epsQuickExtra := 0.5:

lowestIntervalMax := 0.185:
highestIntervalMin := 0.78:

polyAccurateTDCoeffsLowest := 13:
polyAccurateDDCoeffsLowest := 12:
polyAccurateDCoeffsLowest := 16:

polyAccurateDegreeLowest := polyAccurateTDCoeffsLowest + polyAccurateDDCoeffsLowest + polyAccurateDCoeffsLowest:

polyQuickDDCoeffsLowest := 12:
polyQuickDCoeffsLowest := 10:

polyQuickDegreeLowest := polyQuickDDCoeffsLowest + polyQuickDCoeffsLowest:

polyAccurateTDCoeffsMiddle := 7:
polyAccurateDDCoeffsMiddle := 9:
polyAccurateDCoeffsMiddle := 19:

polyAccurateDegreeMiddle := polyAccurateTDCoeffsMiddle + polyAccurateDDCoeffsMiddle + polyAccurateDCoeffsMiddle:

polyQuickDDCoeffsMiddle := 7:
polyQuickDCoeffsMiddle := 7:

polyQuickDegreeMiddle := polyQuickDDCoeffsMiddle + polyQuickDCoeffsMiddle:

polyAccurateTDCoeffsHighest := 9:
polyAccurateDDCoeffsHighest := 9:
polyAccurateDCoeffsHighest := 11:

polyAccurateDegreeHighest := polyAccurateTDCoeffsHighest + polyAccurateDDCoeffsHighest + polyAccurateDCoeffsHighest:

polyQuickDDCoeffsHighest := 9:
polyQuickDCoeffsHighest := 9:

polyQuickDegreeHighest := polyQuickDDCoeffsHighest + polyQuickDCoeffsHighest:

extrabound := hexa2ieee(["3F500000","00000000"]):
polyQuickDegreeExtra := 5:

bound[0] := 0:
b := nearest(lowestIntervalMax):
he := ieeehexa(b):
bound[1] := hexa2ieee([he[1],"00000000"]):
b := nearest(highestIntervalMin):
he := ieeehexa(b):
bound[intervals - 1] := hexa2ieee([he[1],"00000000"]):
bound[intervals] := 1:
linearWidth := (highestIntervalMin - lowestIntervalMax) / (intervals - 2):
scaleSum := 0:
scale[2] := 1.5:
scale[3] := 1.35:
scale[4] := 1.18:
scale[5] := 1.00:
scale[6] := 0.915:
scale[7] := 0.74:
scale[8] := 0.595:
scale[9] := 0.50:
for i from 2 to (intervals-1) do
	scaleSum := scaleSum + scale[i]:
od:
for i from 2 to (intervals-1) do
	scale[i] := scale[i] / scaleSum;
od:
current := lowestIntervalMax:
for i from 2 to (intervals-2) do
	b := nearest(current + (highestIntervalMin - lowestIntervalMax) * scale[i]):
	current := b:
	he := ieeehexa(b):
	bound[i] := hexa2ieee([he[1],"00000000"]):
od:
printf("Using %d intervals with bounds:\n",intervals):
for i from 0 to intervals do
	printf("bound[%d] = %1.30e\n",i,bound[i]):
od:
printf("Using an extra bound for truncating the quick poly to deg. %d in interval #0 for small args up to %f (2^(%f))\n",
	polyQuickDegreeExtra,extrabound,log[2](abs(extrabound)));

printf("Computing Gal's accurate table values for interval midpoints\n"):
for i from 2 to (intervals-1) do
	m := nearest((bound[i] + bound[i-1]) / 2):
	mhe := ieeehexa(m):
	printf("Interval %d: accurate table research start value %f (%s%s)\n",i,m,mhe[1],mhe[2]):
	midpointFloat[i] := galDoubleToDoubleDouble(nearest(m),arcsin,2^(-121),2^(15)):
od:

printf("Using the following floating point midpoints for intervals 2-%d:\n",intervals-2):
for i from 2 to (intervals-1) do
	mhe := ieeehexa(midpointFloat[i]):
	printf("midpointFloat[%d] = %f (%s%s)\n",i,midpointFloat[i],mhe[1],mhe[2]):
od:

printf("The reduced argument z is therefore bounded by:\n"):
printf("Interval 1: |z| < 2^(%f)\n",
	log[2](abs(bound[1]))):
for i from 2 to (intervals-1) do
	printf("Interval %d: |z| < 2^(%f)\n",i,
	log[2](max(abs(midpointFloat[i] - bound[i-1]),abs(midpointFloat[i] - bound[i])))):
od:
printf("Interval %d: |z| < 2^(%f)\n",intervals,
	log[2](abs(1 - bound[intervals-1]))):

printf("Using a %d degree polynomial for lowest interval (1) (accurate phase)\n",
	polyAccurateDegreeLowest):
printf("with %d triple-double, %d double-double and %d double coefficients\n",
	polyAccurateTDCoeffsLowest, polyAccurateDDCoeffsLowest, polyAccurateDCoeffsLowest):
printf("Using a %d degree polynomial for lowest interval (1) (quick phase)\n",
	polyQuickDegreeLowest):
printf("with %d double-double and %d double coefficients\n",
	polyQuickDDCoeffsLowest, polyQuickDCoeffsLowest):
printf("Using a %d degree polynomial for middle intervals (2-%d) (accurate phase)\n",
	polyAccurateDegreeMiddle,intervals-2):
printf("with %d triple-double, %d double-double and %d double coefficients\n",
	polyAccurateTDCoeffsMiddle, polyAccurateDDCoeffsMiddle, polyAccurateDCoeffsMiddle):
printf("Using a %d degree polynomial for middle intervals (2-%d) (quick phase)\n",
	polyQuickDegreeMiddle,intervals-2):
printf("with %d double-double and %d double coefficients\n",
	polyQuickDDCoeffsMiddle, polyQuickDCoeffsMiddle):
printf("Using a %d degree polynomial for highest interval (%d) (accurate phase)\n",
	polyAccurateDegreeHighest,intervals):
printf("with %d triple-double, %d double-double and %d double coefficients\n",
	polyAccurateTDCoeffsHighest, polyAccurateDDCoeffsHighest, polyAccurateDCoeffsHighest):
printf("Using a %d degree polynomial for highest interval (%d) (quick phase)\n",
	polyQuickDegreeHighest,intervals):
printf("with %d double-double and %d double coefficients\n",
	polyQuickDDCoeffsHighest, polyQuickDCoeffsHighest):

if true then

printf("Computing polynomials for interval 1 ([%f;%f])\n",bound[0],bound[1]):

f := unapply(convert(series((arcsin(sqrt(x))/sqrt(x) - 1)/x,x=0,300),polynom),x):

p := unapply(eval(numapprox[chebyshev](f(x),x=(bound[0])^2..(bound[1])^2,2^(-127))),x):

polyAccuExact[1] := truncPoly(subs(X=x^2,p(X))*x^3+x,polyAccurateDegreeLowest):

polyAccurate[1] := poly_exact32(polyAccuExact[1],polyAccurateTDCoeffsLowest,polyAccurateDDCoeffsLowest):

epsAccurate[1] := numapprox[infnorm]((polyAccurate[1]/arcsin(x))-1, x=bound[0]..bound[1]):
epsAccurateSpecial := numapprox[infnorm]((polyAccurate[1]/arcsin(x))-1, x=bound[0]..evalf(sin(2^(-18)))):

polyQuickExact[1] := truncPoly(polyAccuExact[1],polyQuickDegreeLowest):
polyQuick[1] := poly_exact2(polyQuickExact[1],polyQuickDDCoeffsLowest):

epsQuick[1] := numapprox[infnorm]((polyQuick[1]/arcsin(x))-1, x=bound[0]..bound[1]):

polyQuickExtraExact := truncPoly(polyAccuExact[1],polyQuickDegreeExtra):
polyQuickExtra := poly_exact2(polyQuickExtraExact,polyQuickDegreeExtra):
epsQuickExtra := numapprox[infnorm]((polyQuickExtra/arcsin(x))-1, x=bound[0]..extrabound):

end if:

for i from 2 to (intervals-1) do

printf("Computing polynomials for interval %d ([%f;%f])\n",i,bound[i-1],bound[i]):
printf("Reduced argument z will be in interval [%1.8e;%1.8e] ([-2^%f,2^%f]),\n",
	bound[i-1]-midpointFloat[i],bound[i]-midpointFloat[i],
	log[2](abs(bound[i-1]-midpointFloat[i])),log[2](abs(bound[i]-midpointFloat[i]))):

f := unapply(arcsin(x+midpointFloat[i]),x):

fhelp := unapply(convert(series((f(x)-arcsin(midpointFloat[i]))/x,x=0,polyAccurateDegreeMiddle*3),polynom),x):

polyAccuExact[i] := numapprox[minimax](fhelp(x),
				x=bound[i-1]-midpointFloat[i]..bound[i]-midpointFloat[i],
				[polyAccurateDegreeMiddle,0],1,'err')*x + nearestDD(evalf(arcsin(midpointFloat[i]))):

polyAccurate[i] := poly_exact32(polyAccuExact[i],polyAccurateTDCoeffsMiddle,polyAccurateDDCoeffsMiddle):

epsAccurate[i] := numapprox[infnorm]((polyAccurate[i]/f(x))-1,
					x=bound[i-1]-midpointFloat[i]..bound[i]-midpointFloat[i]):

polyQuickExact[i] := truncPoly(polyAccuExact[i],polyQuickDegreeMiddle):
polyQuick[i] := poly_exact2(polyQuickExact[i],polyQuickDDCoeffsMiddle):

epsQuick[i] := numapprox[infnorm]((polyQuick[i]/f(x))-1,
			 	   x=bound[i-1]-midpointFloat[i]..bound[i]-midpointFloat[i]):

od:

if true then

printf("Computing polynomials for interval %d ([%f;%f])\n",intervals,bound[intervals-1],bound[intervals]):

g := unapply(((arcsin(1 - x) - Pi/2)/sqrt(2*x)),x):
f := unapply(convert(series(((g(x)+1)/x),x=0,polyAccurateDegreeHighest*4),polynom),x):

polyAccuExact[intervals] := numapprox[minimax](f(x),x=(1-bound[intervals]+2^(-53))..(1-bound[intervals-1]),
						[polyAccurateDegreeHighest-1,0],1,'err')*x-1:

polyAccurate[intervals] := poly_exact32(polyAccuExact[intervals],polyAccurateTDCoeffsHighest,polyAccurateDDCoeffsHighest):

epsAccurate[intervals] := numapprox[infnorm](((unapply(polyAccurate[intervals],x)(x)/g(x))-1),
					x=(1-bound[intervals]+2^(-53))..(1-bound[intervals-1])):

polyQuickExact[intervals] := truncPoly(polyAccuExact[intervals],polyQuickDegreeHighest):
polyQuick[intervals] := poly_exact2(polyQuickExact[intervals],polyQuickDDCoeffsHighest):

epsQuick[intervals] := numapprox[infnorm](((unapply(polyQuick[intervals],x)(x)/g(x))-1),
					x=(1-bound[intervals]+2^(-53))..(1-bound[intervals-1])):

printf("Checking if the polynomial for interval %d is exactly -1 in z = %f...\n",intervals,1-bound[intervals]);

if (unapply(polyAccurate[intervals],x)(1-bound[intervals]) = -1) then
	printf("  Check passed!\n"):
else
	printf("  Check failed!\n"):
end if:

end if:

for i from 1 to intervals do
	printf("Relative error for accurate phase polynomial in interval %d ([%f;%f]) is 2^(%f)\n",
       		i,bound[i-1],bound[i],log[2](abs(epsAccurate[i]))):
	printf("Relative error for quick phase polynomial in interval %d ([%f;%f]) is 2^(%f)\n",
       		i,bound[i-1],bound[i],log[2](abs(epsQuick[i]))):
od:
printf("Relative error for accurate phase polynomial #1 in special interval [0;sin(2^(-18))]) is 2^(%f)\n",
       log[2](abs(epsAccurateSpecial))):

printf("Relative error for quick phase extra case truncated polynomial #1 in special interval [0;%f]) is 2^(%f)\n",
       extrabound,log[2](abs(epsQuickExtra))):

(PiHalfH, PiHalfM, PiHalfL) := hi_mi_lo(evalf(Pi/2)):

epsPiDD := evalf(((PiHalfH + PiHalfM) - Pi/2)/(Pi/2)):
epsPiTD := evalf(((PiHalfH + PiHalfM + PiHalfL) - Pi/2) / (Pi/2)):

printf("Relative error for storing Pi/2 as a double-double is 2^(%f)\n",evalf(log[2](abs(epsPiDD)))):
printf("Relative error for storing Pi/2 as a triple-double is 2^(%f)\n",evalf(log[2](abs(epsPiTD)))):

# Ce qui suit est pifometrique et doit etre prouve en Gappa ensuite

arithmeticalErrorQuick[1] := 2^(-61):
arithmeticalErrorQuick[2] := 2^(-61):
arithmeticalErrorQuick[3] := 2^(-61):
arithmeticalErrorQuick[4] := 2^(-61):
arithmeticalErrorQuick[5] := 2^(-61):
arithmeticalErrorQuick[6] := 2^(-61):
arithmeticalErrorQuick[7] := 2^(-61):
arithmeticalErrorQuick[8] := 2^(-61):
arithmeticalErrorQuick[9] := 2^(-61):
arithmeticalErrorQuick[10] := 2^(-68):

arithmeticalErrorQuickExtra := 2^(-80):

for i from 1 to intervals do
	estimatedOverallEpsQuick[i] := abs(epsQuick[i]) +
				       abs(arithmeticalErrorQuick[i]) +
                                       abs(epsQuick[i]) * abs(arithmeticalErrorQuick[i]):
	printf("Relative quick phase overall error bound to show in Gappa for interval %d ([%f;%f]) is 2^(%f)\n",
		i,bound[i-1],bound[i],log[2](abs(estimatedOverallEpsQuick[i]))):
od:

estimatedOverallEpsQuickExtra := abs(epsQuickExtra) +
				 abs(arithmeticalErrorQuickExtra) +
                                 abs(epsQuickExtra) * abs(arithmeticalErrorQuickExtra):
printf("Relative quick phase overall error bound to show for extra truncted poly in interval [0;%f]) is 2^(%f)\n",
	extrabound,log[2](abs(estimatedOverallEpsQuickExtra))):

printf("Write tables...\n"):

filename:="TEMPASIN/asin-td.h":
fd:=fopen(filename, WRITE, TEXT):

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

fprintf(fd, "\n/* File generated by maple/asin-td.mpl */\n\n"):

printf("Write table with bounds\n"):

fprintf(fd, "/* High order words of interval bounds (low order word is 0) */\n"):
for i from 1 to (intervals-1) do
	heb := ieeehexa(bound[i]):
	fprintf(fd, "\#define BOUND%d 0x%s\n",i,heb[1]):
od:

heb := ieeehexa(extrabound):
fprintf(fd, "\#define EXTRABOUND 0x%s\n",heb[1]):

printf("Write additional constants\n"):

fprintf(fd, "\n\n/* Pi/2 as a triple-double*/\n"):
fprintf(fd, "\#define PIHALFH %1.50e\n",PiHalfH):
fprintf(fd, "\#define PIHALFM %1.50e\n",PiHalfM):
fprintf(fd, "\#define PIHALFL %1.50e\n",PiHalfL):

printf("Write table with midpoints and polynomial coefficients\n"):
k := 0:
for i from 0 to polyAccurateDegreeLowest do
	(hi,mi,lo) := hi_mi_lo(coeff(polyAccurate[1],x,i)):
	if ((abs(hi) = 1.0) and (mi = 0) and (lo = 0)) then
		printf(
		"Coefficient %d of interval 1 polynomial is exactly %f and will not be stored in the table\n",i,hi):
	else
	g := 0;
	if (hi <> 0) then
		if (i <= polyQuickDegreeLowest) then
			tbl[k] := sprintf("%1.50e, \t/* %d, accPolyLowC%dh, quickPolyLowC%dh */",hi,k,i,i):
		else
			tbl[k] := sprintf("%1.50e, \t/* %d, accPolyLowC%dh */",hi,k,i):
		end if:
		k := k + 1:
	end if:
	if (mi <> 0) then
		if (i <= polyQuickDDCoeffsLowest) then
			tbl[k] := sprintf("%1.50e, \t/* %d, accPolyLowC%dm, quickPolyLowC%dl */",mi,k,i,i):
		else
			tbl[k] := sprintf("%1.50e, \t/* %d, accPolyLowC%dm */",mi,k,i):
		end if:
		k := k + 1:
	end if:
	if (lo <> 0) then
		tbl[k] := sprintf("%1.50e, \t/* %d, accPolyLowC%dl */",lo,k,i):
		k := k + 1:
	end if:
	end if:
od:
tbl[k] := sprintf("%1.50e, \t/* %d, RN rounding constant quick poly low*/",
		compute_rn_constant(estimatedOverallEpsQuick[1]),k):
k := k + 1:
tbl[k] := sprintf("%1.50e, \t/* %d, RD rounding constant quick poly low*/",
		estimatedOverallEpsQuick[1],k):
k := k + 1:
printf("Table for interval 1 written\n"):
for l from 2 to (intervals-1) do
	tbl[k] := sprintf("%1.50e, \t/* %d, midpoint in interval %d*/",midpointFloat[l],k,l):
	tblidx[l] := k;
	k := k + 1;
	for i from 0 to polyAccurateDegreeMiddle do
		(hi,mi,lo) := hi_mi_lo(coeff(polyAccurate[l],x,i)):
		if ((abs(hi) = 1.0) and (mi = 0) and (lo = 0)) then
			printf(
		"Coefficient %d of interval %d polynomial is exactly %f and will not be stored in the table\n",i,l,hi):
		else
		g := 0;
		if (hi <> 0) then
			if (i <= polyQuickDegreeMiddle) then
				tbl[k] := sprintf(
					"%1.50e, \t/* %d, accPolyMid%dC%dh, quickPolyMid%dC%dh */",hi,k,l,i,l,i):
			else
				tbl[k] := sprintf("%1.50e, \t/* %d, accPolyMid%dC%dh */",hi,k,l,i):
			end if:
			k := k + 1:
		end if:
		if (mi <> 0) then
			if (i <= polyQuickDDCoeffsMiddle) then
				tbl[k] := sprintf(
					"%1.50e, \t/* %d, accPolyMid%dC%dm, quickPolyMid%dC%dl */",mi,k,l,i,l,i):
			else
				tbl[k] := sprintf("%1.50e, \t/* %d, accPolyMid%dC%dm */",mi,k,l,i):
			end if:
			k := k + 1:
		end if:
		if (lo <> 0) then
			tbl[k] := sprintf("%1.50e, \t/* %d, accPolyMid%dC%dl */",lo,k,l,i):
			k := k + 1:
		end if:
		end if:
	od:
	tbl[k] := sprintf("%1.50e, \t/* %d, RN rounding constant quick poly middle %d*/",
			compute_rn_constant(estimatedOverallEpsQuick[l]),k,l):
	k := k + 1:
	tbl[k] := sprintf("%1.50e, \t/* %d, RD rounding constant quick poly middle %d*/",
			estimatedOverallEpsQuick[l],k,l):
	k := k + 1:
	printf("Table for interval %d written\n",l):
od:
tblidx[intervals] := k:
for i from 0 to polyAccurateDegreeHighest do
	(hi,mi,lo) := hi_mi_lo(coeff(polyAccurate[intervals],x,i)):
	if ((abs(hi) = 1.0) and (mi = 0) and (lo = 0)) then
		printf(
	"Coefficient %d of interval %d polynomial is exactly %f and will not be stored in the table\n",i,intervals,hi):
	else
	g := 0;
	if (hi <> 0) then
		if (i <= polyQuickDegreeHighest) then
			tbl[k] := sprintf("%1.50e, \t/* %d, accPolyHighC%dh, quickPolyHighC%dh */",hi,k,i,i):
		else
			tbl[k] := sprintf("%1.50e, \t/* %d, accPolyHighC%dh */",hi,k,i):
		end if:
		k := k + 1:
	end if:
	if (mi <> 0) then
		if (i <= polyQuickDDCoeffsHighest) then
			tbl[k] := sprintf("%1.50e, \t/* %d, accPolyHighC%dm, quickPolyHighC%dl */",mi,k,i,i):
		else
			tbl[k] := sprintf("%1.50e, \t/* %d, accPolyHighC%dm */",mi,k,i):
		end if:
		k := k + 1:
	end if:
	if (lo <> 0) then
		tbl[k] := sprintf("%1.50e, \t/* %d, accPolyHighC%dl */",lo,k,i):
		k := k + 1:
	end if:
	end if:
od:
tbl[k] := sprintf("%1.50e, \t/* %d, RN rounding constant quick poly high*/",
		compute_rn_constant(estimatedOverallEpsQuick[intervals]),k):
k := k + 1:
tbl[k] := sprintf("%1.50e, \t/* %d, RD rounding constant quick poly high*/",
 		estimatedOverallEpsQuick[intervals],k):
k := k + 1:
printf("Table for interval %d written\n",intervals):
tbllen := k:
printf("The whole table has %d entries, so uses %d bytes of memory\n",tbllen,tbllen*8):
fprintf(fd,"\n\n/* Indices to the following table */\n"):
for i from 2 to intervals do
	fprintf(fd,"\#define TBLIDX%d %d\n",i,tblidx[i]):
od:
fprintf(fd, "\n\n/* Table with midpoints and polynomial coefficients */\n"):
fprintf(fd, "static const double tbl[%d] = {\n",tbllen):
for i from 0 to (tbllen - 1) do
	fprintf(fd, "%s\n",tbl[i]):
od:
fprintf(fd, "};\n\n"):

fclose(fd):

printf("----DONE---\n"):


text:
(*
 Copyright (C) 2017 M.A.L. Marques
               2019 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: gga_exc *)
(* prefix:
  gga_c_pbe_vwn_params *params;

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

$include "lda_c_vwn.mpl"
$ifdef gga_c_pbe_params
params_a_beta  := 0.06672455060314922:
params_a_gamma := (1 - log(2))/Pi^2:
params_a_BB    := 1:
$endif

mgamma := params_a_gamma:
mbeta  := (rs, t) -> params_a_beta:
BB     := params_a_BB:

tp   := (rs, z, xt) -> tt(rs, z, xt):

(* Equation (8) *)
A := (rs, z, t) ->
  mbeta(rs, t)/(mgamma*(exp(-f_vwn(rs, z)/(mgamma*mphi(z)^3)) - 1)):

(* Equation (7) *)
f1 := (rs, z, t) -> t^2 + BB*A(rs, z, t)*t^4:
f2 := (rs, z, t) -> mbeta(rs, t)*f1(rs, z, t)/(mgamma*(1 + A(rs, z, t)*f1(rs, z, t))):

fH := (rs, z, t) -> mgamma*mphi(z)^3*log(1 + f2(rs, z, t)):

f_pbe  := (rs, z, xt, xs0, xs1) ->
  f_vwn(rs, z) + fH(rs, z, tp(rs, z, xt)):

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


text:
input := FileTools:-Text:-ReadFile("AoC-2021-8-input.txt" ):

lines := subs(""=NULL,(op@StringTools:-Split)~(StringTools:-Split(input, "\n"),"|"));

# part 1, just count the 1, 4, 7, and 8s
total := 0:
for i from 2 to nops(lines) by 2 do
    tal := table(sparse=0,Statistics:-Tally(map(length, StringTools:-Split(lines[i]))));
    total := total + tal[2] + tal[3] + tal[4] + tal[7];
end do:

answer1 := total;

# part 2

sortstring := s->StringTools:-Join( sort( StringTools:-Explode(s) ), "");
lookups := ["abcefg"="0", "acdeg"="2", "acdfg"="3", "abdfg"="5", "abdefg"="6", "abcdfg"="9"];
alphab := {seq(j,j in "abcdefg")};

Decode := proc( code, output )
local j, displays, cypher, i, decoded, w, wl, l, wc, n;

    displays := sort( StringTools:-Split(StringTools:-Trim(code)," "), length );

    cypher := table([seq(i={seq(j,j in "abcdefg")},i in "abcdefg")]):
    decoded := table();

    for w in displays do
        wl := {seq(i, i in w)};
        if length(w) = 2 and not assigned(decoded["1"]) then # 1
            decoded["1"] := sortstring(w);
            # may be
            cypher["c"] := cypher["c"] intersect wl;
            cypher["f"] := cypher["f"] intersect wl;
            # can't be
            for l in "abdeg" do
                cypher[l] := cypher[l] minus wl;
            end do;
        elif length(w) = 3 and not assigned(decoded["7"])then # 7
            decoded["7"] := sortstring(w);
            cypher["c"] := cypher["c"] intersect {seq(i, i in w)};
            cypher["f"] := cypher["f"] intersect {seq(i, i in w)};
            cypher["a"] := cypher["a"] intersect {seq(i, i in w)};
            if assigned(decoded["2"]) then
                cypher["a"] := cypher["a"] intersect {seq(i, i in w)} minus cypher["c"];
            end if;
            if nops( cypher["a"] ) = 1 then
                for l in alphab minus {"a"} do
                    cypher[l] := cypher[l] minus cypher["a"];
                end do;
            end if;
            for l in "bdeg" do
                cypher[l] := cypher[l] minus wl;
            end do;
        elif length(w) = 4 and not assigned(decoded["4"])then # 4
            decoded["4"] := sortstring(w);
            cypher["b"] := cypher["b"] intersect {seq(i, i in w)};
            cypher["d"] := cypher["d"] intersect {seq(i, i in w)};
            cypher["c"] := cypher["c"] intersect {seq(i, i in w)};
            cypher["f"] := cypher["f"] intersect {seq(i, i in w)};
            for l in "aeg" do
                cypher[l] := cypher[l] minus wl;
            end do;
        elif length(w) = 7 and not assigned(decoded["8"]) then # 8
            decoded["8"] := sortstring(w);
            cypher["b"] := cypher["b"] intersect {seq(i, i in w)};
            cypher["d"] := cypher["d"] intersect {seq(i, i in w)};
            cypher["c"] := cypher["c"] intersect {seq(i, i in w)};
            cypher["f"] := cypher["f"] intersect {seq(i, i in w)};
        elif length(w) = 6 then # 0, 6, 9 missing letters must be c,d,e
            wc := alphab minus wl;
            for l in "abf" do
                cypher[l] := cypher[l] minus wc;
            end do;
        elif length(w) = 5 then # 2, 3, 5 missing letters must be b,c,e,f
            wc := alphab minus wl;
            for l in "adg" do
                cypher[l] := cypher[l] minus wc;
            end do;
        end if;
    end do:

    for l in alphab do
        if nops(cypher[l]) = 1 then
            for n in alphab minus {l} do
                cypher[n] := cypher[n] minus cypher[l];
            end do;
        end if;
    end do;

    for l in alphab do
        if nops(cypher[l]) <> 1 then
            error "need more!";
        end if;
        cypher[l] := op(cypher[l]);
    end do:

    cypher := sort((rhs=lhs)~([entries(cypher,pairs)]));
    decoded := sort((rhs=lhs)~([entries(decoded,pairs)]));

    tmp := map(sortstring, StringTools:-Split(output," "));
    tmp := map2(subs, decoded, tmp);
    tmp := (StringTools:-Join([seq]( subs(lookups, StringTools:-Join(sort( subs(cypher, convert(w,list)) ),"")), w in tmp), "") );

    return parse(tmp);

end proc:

answer2 := add( Decode(lines[i], lines[i+1]), i=1..nops(lines), 2);


text:
 func $cvtu64tof64 (
  var %i u64
  ) f64 {
   return (
     cvt f64 u64(dread u64 %i))}
 # EXEC: %irbuild Main.mpl
 # EXEC: %irbuild Main.irb.mpl
 # EXEC: %cmp Main.irb.mpl Main.irb.irb.mpl


text:
############################################################
# A function evaulator routine that
# finds the formation of trap surface
# i.e. blackhole formation, see the paper
# for its definition: http://arxiv.org/pdf/1508.01614v2.pdf
############################################################
read "/home/arman/FD/FD.mpl":

CFD();
MFD();

grid_functions:= {a1,b1,alpha,beta,ctfm,ctfmp};

AVGT := f -> (FD(f,[[1],[0,0,0]]) + FD(f,[[0],[0,0,0]]))/2 ;
FD_table[t] := [ [0] , [0,1] ];

spec := [
  { i = [1,1,1] , b=xmin } = 1+myzero*x(i),
  { i = [2,Nx-1,1]       } = (AVGT(Gen_Sten(1/alpha(t,x)))*Gen_Sten(diff(b1(t,x),t)) + AVGT( Gen_Sten( (1/a1(t,x) - beta(t,x)/alpha(t,x))*( diff(b1(t,x),x)/ctfmp(x) + b1(t,x)/ctfm(x) ) ) ) )/ (AVGT(Gen_Sten(b1(t,x)))),
  { i=[Nx,Nx,1], b=xmax  } = myzero*x(i)
];

A_Gen_Eval_Code(spec,input="d",proc_name="compute_th");