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 */
00096 polyadic2int(N,b) :=
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
00103   -polyadic2int(map(digit2int,str2chrl(substring(str,2))),b)
00104 else
00105   polyadic2int(map(digit2int,str2chrl(str)),b)$
00106 binv2int(l) := polyadic2int(l,2)$
00107 bin2int(str) := polyadicstr2int(str,2)$
00108 oct2int(str) := polyadicstr2int(str,8)$
00109 dec2int(str) := polyadicstr2int(str,10)$
00110 hex2int(str) := polyadicstr2int(str,16)$
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). */
00117 int2polyadic(m,b) :=
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: */
00125 int2polyadic_padd(m,b,l) := paddingfront_l(0,int2polyadic(m,b),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
00131   simplode(map(int2digit,int2polyadic(m,b)))
00132 else
00133   sconcat("-", simplode(map(int2digit,int2polyadic(-m,b))))$
00134 int2bin(m) := int2polyadicstr(m,2)$
00135 int2boct(m) := int2polyadicstr(m,8)$
00136 int2dec(m) := int2polyadicstr(m,10)$
00137 int2hex(m) := int2polyadicstr(m,16)$
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: */
00142 hexstr2binv(str) := lappend(map(lambda([x],int2polyadic_padd(x,2,4)),map(digit2int,str2chrl(supcase(str)))))$
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 */
00150 binv2hexstr(v) := simplode(map(int2digit,map(lambda([x],polyadic2int(x,2)),partition_elements(v,4))))$
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
00156    http://en.wikipedia.org/wiki/Factoradic):
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 
00162 int2factoradic(m) :=
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