OKlibrary  0.2.1.6
Auxiliary.mac
Go to the documentation of this file.
```00001 /* Oliver Kullmann, 7.6.2008 (Swansea) */
00002 /* Copyright 2008, 2009, 2010, 2011, 2012, 2013 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/Hypergraphs/Lisp/SetSystems.mac")\$
00023 oklib_include("OKlib/ComputerAlgebra/DataStructures/Lisp/Lists.mac")\$
00024
00025
00026 /* **************
00027    * Logarithms *
00028    **************
00029 */
00030
00031 /* The binary logarithm: */
00032 ld(x) := radcan(log(x) / log(2))\$
00033 /* The decimal logarithm: */
00034 log10(x) := radcan(log(x) / log(10))\$
00035
00036 /* The integer binary logarithm for natural numbers n (minf for n=0). */
00037 /* This is just floor(ld(n)) (for n <> 0), however our computation is
00038    actually faster (except for very big n, say, with over 5000 decimal digits).
00039 */
00040 fld(n) := if n = 0 then minf else
00041  block([l : 0, b : 1], while b < n do (l : l + 1, b : b * 2),
00042   if b > n then return(l-1) else return(l))\$
00043 /* Now ceiling(ld(n)), again with faster computation: */
00044 cld(n) := if n = 0 then minf else
00045  block([l : 0, b : 1], while b < n do (l : l + 1, b : b * 2), l)\$
00046
00047 /* Returns the floor of the logarithm of n base b, correcting
00048    for issues in precision causing round off. Cannot handle large
00049    numbers. */
00050 floorlog(n,b) := block([l : floor(float(log(n)/log(b)))],
00051   if b^(l+1) <= n then (l+1) else
00052   if b^(l) > n then (l-1) else l)\$ /* handles issues with precision */
00053 /* We have
00054    floorlog(n,b) = floor(log(n)/log(b)),
00055    however for n for which floorlog(n,b) is defined it is much faster.
00056 */
00057
00058
00059 /* **********************
00060    * Sign functionality *
00061    **********************
00062 */
00063
00064 posp(x) := is(x > 0)\$
00065 negp(x) := is(x < 0)\$
00066 nullp(x) := is(x = 0)\$
00067 possignum(x) := if posp(x) then 1 else 0\$
00068 negsignum(x) := if negp(x) then -1 else 0\$
00069 /* Reminder: signum(x) = +1,0,-1 if resp. posp(x), nullp(x), negp(x) is true.
00070 */
00071
00072 /* ******************************
00073    * Representations of numbers *
00074    ******************************
00075 */
00076
00077 /* Converts a string to a list of characters: */
00078 str2chrl(str) := create_list(charat(str, i), i,1,slength(str))\$
00079
00080 /* Converts a digit of the form "0", "1", ..., "9", "A", ..., "Z" to an
00081    integer from 0 to 35: */
00082 digit2int(d) :=
00083   if cgreaterp("A",d) then cint(d) - cint("0")
00084   else 10 + (cint(d) - cint("A"))\$
00085 /* The inverse, converting an integer from 0 to 35 to a digit: */
00086 int2digit(m) :=
00087     if m < 10 then ascii(cint("0") + m)
00088     else ascii(cint("A") + (m-10))\$
00089
00090 /* Transforming a list N representing a natural number n to base b (big endian,
00091    that is leading coefficients first) into n. */
00092 /* More generally, N can interpreted as the list of coefficients of a
00093    polynomial in an unknown x, but with leading coefficient first, and
00094    this polynomial then is evaluating (via Horner scheme) at argument value b.
00095 */
00097   lreduce(lambda([x,y], x * b + y), N, 0)\$
00098 /* Converts an non-empty string representing an integer in base b to an
00099    integer. */
00100 /* Prerequisite: b a natural number, 1 <= b <= 36; "digits" might be equal to b
00101    or even be greater (but must be characters from 0-9 and A-Z). */
00102 polyadicstr2int(str,b) := if charat(str,1) = "-" then
00104 else
00111
00112 /* Computing the representation of a number m to base b as a list
00113    (big endian, so for example int2polyadic(100,10) = [1,0,0]):
00114 */
00115 /* Prerequisites: m a natural number >= 0, b a natural number >= 2
00116    (where in case m=0 also b in {0,1} is possible). */
00118   if m = 0 then [0]
00119   else block([result : []],
00120     while m > 0 do block([d : divide(m,b)],
00121       result : cons(d[2], result),
00122       m : d[1]),
00123     return(result))\$
00124 /* Adding leading zeros if the result does not have length l: */
00126
00127
00128 /* Converts a given integer m to its base b representation as a
00129    string, where 2 <= b <= 36 is a natural number. */
00130 int2polyadicstr(m,b) := if m >= 0 then
00132 else
00138
00139 /* Converting a hexadecimal string str into a binary vector, where str can use
00140    upper- or lower-case letters (also mixed); each hex-digit yields 4 bits,
00141    where the bits in the 4-bit segment are to be read from right to left: */
00143 /* Example: hexstr2binv("A12") = [1,0,1,0,0,0,0,1,0,0,1,0]. */
00144
00145 /* Converting a binary vector v into a hexadecimal string (using upper-case);
00146    v is partitioned into segments of length 4, where in each segment the bits
00147    are to be read from right to left, and where a trailing segment of length
00148    1, 2 or 3 is left-padded with 0:
00149 */
00151 /* Example: binv2hexstr([1,1,1,1,1]) = "F1" = binv2hexstr([1,1,1,1,0,0,0,1]).
00152 */
00153
00154
00155 /* The factorial number representation (see
00157 */
00158 factoradic2int(N) := block([res : 0, fac : 1, i : 0],
00159  for d in reverse(N) do (res : d * fac + res, i : i + 1, fac : fac * i),
00160  res)\$
00161
00163  block([res : [], i : 1, div],
00164   do (
00165     div : divide(m,i),
00166     res : cons(div[2], res),
00167     m : div[1],
00168     if m=0 then return(res),
00169     i : i + 1
00170   ))\$
00171
00172
00173 /* ***************
00174    * Conversions *
00175    ***************
00176 */
00177
00178 /* Rounding to d decimal digits: */
00179 round_fdd(x,d) := float(round(x*10^d)/10^d)\$
00180 round_bfdd(x,d) := bfloat(round(x*10^d)/10^d)\$
00181
00182 fractional_part(x) := x - floor(x)\$
00183
00184
00185 /* **********************
00186    * Integer partitions *
00187    **********************
00188 */
00189
00190 /* An "integer partition" of n in ZZ is a list L of natural numbers (>= 1),
00191    sorted in ascending order, which sums up to n.
00192 */
00193
00194 /* The set of all integer partitions of n (an integer): */
00195 ext_integer_partitions(n) :=
00196  if n < 0 then {}
00197  elseif n = 0 then {[]}
00198  else integer_partitions(n)\$
00199 /* Remark: The Maxima-function integer_partitions is false for n=0, and
00200    undefined for n < 0. */
00201 /* The number of integer partitions: */
00202 num_ext_integer_partitions(n) :=
00203  if n < 0 then 0
00204  elseif n = 0 then 1
00205  else num_partitions(n)\$
00206
00207 /* An "unordered integer partition" of n in ZZ is a list L of natural numbers
00208    (>= 1) which sums up to n.
00209 */
00210 /* The set of unordered integer partitions (n is an integer): */
00211 uinteger_partitions(n) :=
00212  lunion(map(permutations,ext_integer_partitions(n)));
00213 num_uinteger_partitions(n) := if n<0 then 0 elseif n=0 then 1
00214  else 2^(n-1)\$
00215
00216
00217 /* ***********************
00218    * Corrected functions *
00219    ***********************
00220 */
00221
00222 /* The corrected power-function: */
00223 pow(b,e) := block([ze : is(equal(e,0))],
00224  if ze=true then return(1)
00225  elseif ze=false then return(b^e)
00226  else block([zb : is(equal(b,0))],
00227    if zb=true then return(if equal(e,0) then 1 else 0)
00228    elseif zb=false then return(b^e)
00229    else return(ev(pow(b,e),noeval))))\$
00230
00231
```