OKlibrary  0.2.1.6
Basic.mac
Go to the documentation of this file.
```00001 /* Oliver Kullmann, 29.7.2007 (Swansea) */
00002 /* Copyright 2007, 2008, 2009 Oliver Kullmann
00003 This file is part of the OKlibrary. OKlibrary is free software; you can redistribute
00004 it and/or modify it under the terms of the GNU General Public License as published by
00005 the Free Software Foundation and included in this library; either version 3 of the
00006 License, or any later version. */
00007
00022 oklib_include("OKlib/ComputerAlgebra/DataStructures/Lisp/Lists.mac")\$
00023
00024
00025 /*
00026  A "branching tuple" is a list of positive numbers, which in case of length
00027  one can also include 0 (also the empty list is possible).
00028 */
00029
00030 /* ************************
00031    * Elementary functions *
00032    ************************
00033 */
00034
00035 bt_min(t) := lmin(t)\$
00036 bt_max(t) := lmax(t)\$
00037 bt_width(t) := length(t)\$
00038 bt_sum(t) := apply("+", t)\$
00039
00040 /* For the concatenation of branching tuples, "append" is to be used. */
00041
00042 /* Composition of branching tuples at position 1: */
00043 bt_composition(a,b) := if emptyp(a) then [] else append(a[1] + b, rest(a))\$
00044 /* Composition of branching tuples at position i: */
00045 bt_composition_i(a,b,i) := if emptyp(a) then []
00046   else append(rest(a, -(length(a) - i + 1)), a[i] + b, rest(a,i))\$
00047
00048 /*
00049  A "permutation" in S_n is a list P of length n of all the elements in
00050  {1,...,n}; it stands for the bijection i -> P[i].
00051  See ComputerAlgebra/Algebra/Lisp/Groupoids/Groups/SymmetricGroups.mac.
00052 */
00053
00054 /* Apply a permutation (or, more generally, an index selection): */
00055 bt_apply_permutation(perm, t) := create_list(t[i], i, perm)\$
00056
00057 /* The branching-tuple-hull of a branching tuple t w.r.t. upper-bound
00058    s on the maximum, together with the set of saturated elements: */
00059 bth(t, s) := if s < bt_max(t) then [{},{}] else
00060  block([hull : {t}, new_size : 1, old_size : 0, saturated : {}],
00061   while new_size > old_size do (
00062     old_size : new_size,
00063     for a in hull unless new_size > old_size do (
00064       hull : union(hull, permutations(a),
00065                map(lambda([b],bt_composition(a,b)),
00066                    subset(hull,lambda([b], is(a[1]+bt_max(b) <= s))))),
00067       new_size : length(hull))),
00068   return([hull, block([mt:bt_max(t)], subset(hull, lambda([a],is(bt_min(a)+mt>s))))]))\$
00069
00070 /* Compute greedily one saturated tuple in the branching-tuple-hull: */
00071 bth_saturated_greedy(t,s) := block([m : bt_max(t), ts : sort(t)],
00072   while ts[1] + m <= s do (
00073     ts : bt_composition(ts,t),
00074     ts : sort(ts)
00075   ),
00076   return(ts))\$
00077
00078 /* The integer sequence corresponding to an integer branching tuple t,
00079    that is, the first n elements.
00080    Prerequisites: t is not empty, and also t <> [0].
00081  */
00082 int_seq_bt(t) := block([m : bt_max(t), l : length(t)],
00083   return(buildq([m,l,t],
00084    lambda([n], block([S : create_list(1,i,1,min(n,m))],
00085      for i : m thru n-1 do
00086        S : endcons(sum(S[i-t[j]+1],j,1,l),S),
00087      return(S)))
00088    )))\$
00089 /* The sequence is f(0),f(1),f(2),...,f(n-1), where
00090    f(n) := f(n-first(t)) + ... + f(m-last(t))
00091    if this expression is defined (i.e., iff n >= max(t)),
00092    while otherwise f(n) := 1.
00093    For example with fibi : int_seq_bt([1,2]) we get
00094    fibi(7) = [1, 1, 2, 3, 5, 8, 13]
00095    (note here that fib is the predefined Fibonacci-function).
00096 */
00097
00098
00099 /* ***************
00100    * Power means *
00101    ***************
00102 */
00103
00104 /* n-ary arithmetic mean (while the arithmetic mean for branching tuples t
00105    is just mean(t)): */
00106 meann([t]) := mean(t)\$
00107
00108 /* Geometric mean: */
00110 gmean(t) := geometric_mean(t)\$
00111 /* n-ary geometric mean: */
00112 gmeann([t]) := gmean(t)\$
00113
00114 /* Harmonic mean: */
00115 hmean(t) := harmonic_mean(t)\$
00116 /* n-ary harmonic mean: */
00117 hmeann([t]) := hmean(t)\$
00118
00119 /* Power means with index -inf <= p <= inf: */
00120 genmean(t, p) := if p = minf then lmin(t) elseif p = 0 then gmean(t) elseif p = inf then lmax(t) else mean(t^p)^(1/p)\$
00121 /* (n+1)-ary power means: */
00122 genmeann([u]) := genmean(rest(u,-1), last(u))\$
00123 /* Binary power mean: */
00124 genmean2(x,y,p) := genmeann(x,y,p)\$
00125 /* Ternary power mean: */
00126 genmean3(x,y,z,p) := genmeann(x,y,z,p)\$
00127
00128
00129 /* ***********************************************
00130    * Lower and upper bounds for the tau-function *
00131    ***********************************************
00132 */
00133
00135
00136 /* The lower bound given by the arithmetic mean: */
00137 tau_lo(t) := length(t)^(1/mean(t))\$
00138 /* n-ary version: */
00139 tau_lon([t]) := tau_lo(t)\$
00140
00141 /* The upper bound given by the power means (defined below): */
00142 tau_up(t) := length(t)^(1/genmean(t,2-length(t)))\$
00143 /* n-ary version: */
00144 tau_upn([t]) := tau_up(t)\$
00145
00146
00147 /* ****************************
00148    * Comparison of tau-values *
00149    ****************************
00150 */
00151
00152 /* bt_comparison(t,alpha) is pos, zero, neg if tau(t) > alpha resp.
00153    = alpha resp. < alpha
00154    (this decision uses only symbolic means, no numeric ones): */
00155 bt_comparison(t,alpha) := sign(chi(t,alpha) - 1)\$
00156
00157
00158 /* *******************************************************
00159    *   Computing the tau-function                        *
00160    *******************************************************
00161 */
00162
00163 /* The chi-functional: */
00164 chi(t,x) := bt_sum(x^(-t))\$
00165
00166 /* For fixed t, compute the Newton-step for chi(t,x)-1 = 0: */
00167 chi_newtonstep(t) := buildq([t],
00168  lambda([x], block([p : x^(-t)],
00169   return((apply("+",p)-1) / apply("+",t*p/x)))))\$
00170
00171 /* Computing the tau-function by the Newton-method: */
00172 tau_eps(t,eps) := block([l:length(t)],
00173  if l = 0 then inf elseif l = 1 then 1 else
00174  block([dxf : chi_newtonstep(t), x : tau_lo(t), dx],
00175   do (
00176     dx : dxf(x),
00177     x : x + dx,
00178     if dx <= eps then return(true)
00179   ),
00180   return(x)))\$
00181
00182 /* Given a branching tuple, compute its tau-value (supplying a standard eps):
00183 */
00184 tau(t) := tau_eps(float(t), 10^-15)\$
00185 /* floatnump(tau(t)) is true */
00186 /* n-ary version (using the standard eps): */
00187 taun([v]) := tau(v)\$
00188
00189 /* Binary and ternary versions: */
00190 tau2(x,y) := taun(x,y)\$
00191 tau3(x,y,z) := taun(x,y,z)\$
00192
00193
00194 /* ********************
00195    * Higher precision *
00196    ********************
00197 */
00198
00199 /* General remarks on floating-point representations:
00200    Use float(x) to compute a floating-point representation of x.
00201    To use higher precision, set fpprec to the required number of
00202    decimal digits, and convert via bfloat.
00203    Or use the following conversion function.
00204 */
00205 /* Convert a number to a number with a given precision: */
00206 hp_float(x, decimal_digits) := block([fpprec : fpprec], fpprec : decimal_digits, bfloat(x))\$
00207
00208 /* The tau-value with a given precision: */
00209 tau_hp(t, decimal_digits) := block(
00210  [fpprec : fpprec],
00211   fpprec : decimal_digits,
00212   tau_eps(bfloat(t), 10^(-decimal_digits+1)))\$
00213
00214
00215 /* ************************
00216    * Symbolic computation *
00217    ************************
00218 */
00219
00220 oklib_plain_include(sqdnst)\$
00221
00222 /* Computes a symbolical solution if possible (otherwise returns false): */
00223 tau_symbolical(t) := block([l : length(t), x, sol, pos_sol, n_pos],
00224  if l = 0 then inf elseif l = 1 then 1
00225  else (sol : map(rhs, solve(chi(t,x)-1, x)),
00226        pos_sol : sublist_indices(sol, lambda([x], is(equal(imagpart(x),0)) and is(sign(x) = pos))),
00227        n_pos : length(pos_sol),
00228        if n_pos = 1 then
00229          return(sqrtdenest(sqrtdenest(sol[pos_sol[1]])))
00230        elseif n_pos >= 2 then print("INCONSISTENCY with tau_symbolical for ", t)
00231        else return(false)))\$
00232 /* Remark: Apply "optimize" to the result to obtain a representation with
00233    shared subterms. */
00234
00235 /* ***********************************
00236    * Derivatives of the tau-function *
00237    ***********************************
00238 */
00239
00240 /* The total differential for the tau-function: */
00241 Dtau(t) := block([tv : tau(t)], block([pv : tv^(-t)],
00242   return((-tv * log(tv) / apply("+", t * pv)) * pv)))\$
00243
00244 Dtaun([t]) := Dtau(t)\$
00245 Dtau2(x,y) := Dtaun(x,y)\$
00246 Dtau3(x,y,z) := Dtaun(x,y,z)\$
00247
00248 /* Symbolic differentiation of tau2: */
00249 Dtau2s(x,y) := block( [ tv : tau2s(x,y)], block( [pv : tv^(-[x,y]) ] ,
00250   return((- tv * log(tv) / apply("+", [x,y] * pv)) * pv)))\$
00251
00253
00254
00255 /* ****************************************
00256    * The induced probability distribution *
00257    ****************************************
00258 */
00259
00260 /* The probability distribution derived from a non-empty branching tuple: */
00261 tauprob(t) := if length(t)=1 then [1] else
00262  block([ts : t / lmin(t)], tau(ts)^(-ts))\$
00263 tauprobn([t]) := tauprob(t)\$
00264 tauprob2(x,y) := tauprobn(x,y)\$
00265 tauprob3(x,y,z) := tauprobn(x,y,z)\$
00266
00267 /* Higher precision: */
00268 tauprob_hp(t, decimal_digits) := if length(t)=1 then [1] else
00269  block([fpprec : fpprec],
00270   fpprec : decimal_digits, t : bfloat(t),
00271   block([ts : t / lmin(t)], tau_hp(ts,decimal_digits)^(-ts)))\$
00272
00273 /* Symbolical computation (if possible): */
00274 tauprob_symbolical(t) := tau_symbolical(t)^(-t)\$
00275 /* Applying factor(expand(r)), radcan(r) or fullratsimp(r) to the result r
00276    might help simplifying. */
00277
00278 /* The canonical inversion for tuples p with bt_sum(p) = 1: */
00279 invtauprob(p) := -map(log,p)\$
00280 /* tauprob(invtprop(p)) = p and invtprop(tauprob(t)) = t * log(tau(t)) */
00281
00282
00283 /* ****************
00284    * The tau-mean *
00285    ****************
00286 */
00287
00288 /* Auxiliary function to test for zero: */
00289 zerop(x) := is(equal(x,0))\$
00290 /* Auxiliary function to test for +infinity: */
00291 infp(x) := is(equal(x,inf))\$
00292
00293 /* The basic tau-mean: */
00294 taumeanb(t) := float(log(length(t))) / log(tau(t))\$
00295 taumeanb2(x,y) := taumeanb([x,y]);
00296 taumeanb3(x,y,z) := taumeanb([x,y,z]);
00297 /* The tau-mean for tuples allowing zero-values: */
00298 taumean(t) := if length(t) <= 1 then inf
00299   elseif some_s(zerop, t) then 0
00300   else float(log(length(t))) / log(tau(t))\$
00301 /* The n-ary tau-mean: */
00302 taumeann([t]) := taumean(t)\$
00303 /* The binary tau-mean: */
00304 taumean2(x,y) := taumeann(x,y)\$
00305 /* The ternary tau-mean: */
00306 taumean3(x,y,z) := taumeann(x,y,z)\$
00307
00308 /* Tau-mean additionally also allowing inf-values: */
00309 taumean_inf(t) := if some_s(infp, t) then
00310     block([t : delete(inf, t)], if t = [0] then unknown else taumean(t))
00311   else taumean(t)\$
00312 taumean_infn([t]) := taumean_inf(t)\$
00313
00314
00315 /* *********************************************
00316    *   Convexity considerations: line versions *
00317    *********************************************
00318 */
00319
00320
00321 /* Considerung for a branching tuple t0 and a direction d0 (a real vector of
00322    the same length as t) the function x -> tau(t0 + x * d0): */
00323
00324 tau_line(x) := (local(x), tau(t0 + x * d0))\$
00325 logtau_line(x) := (local(x), log(tau(t0 + x * d0)))\$
00326 loglogtau_line(x) := (local(x), log(logtau_line(x)))\$
00327 tauprob_line(x) := (local(x), tauprob(t0 + x * d0)[i0])\$
00328 taumean_line(x) := (local(x), taumean(t0 + x * d0))\$
00329 /* assume(length(t0) >= 1), assume(t0[i] > 0) */
00330 taumeandiff_line(x) := (local(x, t),
00331   block([t : t0 + x * d0], taumean(t) - genmean(t, p0))
00332 )\$
00333
00334 /* Returning the interval bounds (for the open interval): */
00335 tau_line_dom() := block(
00336   [pos : sublist_indices(d0, lambda([x], x > 0)), neg : sublist_indices(d0, lambda([x], x < 0))],
00337   [lmax(- part(t0, pos) / part(d0,pos)), lmin(- part(t0, neg) / part(d0, neg))]
00338 )\$
00339
00340
00341 /* The chi-function on a line (but chi(t,x) only for x > 1): */
00342 chi_line(x) := (local(x), chi(t0 + x * d0, a0 + x * b0))\$
00343 /* assume(length(t0) >= 2), assume(t0[i] > 0), assume(a0 > 1) */
00344
00345 chi_line_dom() :=
00346  block([ts : append(t0, [a0 - 1]), ds : append(d0, [b0])],
00347    block(
00348     [pos : sublist_indices(ds, lambda([x], x > 0)), neg : sublist_indices(ds, lambda([x], x < 0))],
00349     [lmax(- part(ts, pos) / part(ds,pos)), lmin(- part(ts, neg) / part(ds, neg))]
00350 ))\$
00351
00352 chin([u]) := chi(rest(u,-1), last(u))\$
00353 chi2(t1,t2,x) := chin(t1,t2,x)\$
00354 chi3(t1,t2,t3,x) := chin(t1,t2,t3,x)\$
00355
00356 hess_chi2 : hessian(chi2(t1,t2,x),[t1,t2,x])\$
00357 /* if for example x has a value, then apply remvalue(x), or, better, kill(x) */
00358
00359 eigen_chi2(a,b,c) := eigens_by_jacobi(at(hess_chi2, [t1=a,t2=b,x=c]))\$
00360 /* a, b > 0, c > 1 */
00361 test_posdef_chi2(a,b,c) := lmin(eigen_chi2(a,b,c)[1])\$
00362
00363
00364 /* ************************************
00365    * Investigations on approximations *
00366    ************************************
00367 */
00368
00369 /* Consider the bounds
00370      genmean(t,k-2) <= taumeanb(t) <= mean(t).
00371 */
00372
00373 /* For k = 2 the lower bound is always better as an approximation than
00374    the upper bound.
00375    Here we investigate how this can be quantified.
00376 */
00377 /* For k >= 3 this is not always the case (example [5.0,5.0,0.835]),
00378    but it appears that still the lower bound yields a better approximation:
00379     - perhaps regarding the integral. ?
00380 */
00381 /* It is now known yet what happens for k > 3. */
00382
00383 dltau_2(x,y) := taumeanb2(x,y) - float(gmean([x,y]))\$
00384 dutau_2(x,y) := mean([x,y]) - taumeanb2(x,y)\$
00385 dtau_2(x,y) := dutau_2(x,y) - dltau_2(x,y)\$
00386 /* For all x,y > 0 we have dtau_2(x,y) >= 0 (equality holds
00387    iff x = y). */
00388
00389 dltau_3(x,y,z) := taumeanb3(x,y,z) - float(hmean([x,y,z]))\$
00390 dutau_3(x,y,z) := mean([x,y,z]) - taumeanb3(x,y,z)\$
00391 dtau_3(x,y,z) := dutau_3(x,y,z) - dltau_3(x,y,z)\$
00392
00393
```