share/contrib/log10.mac

   1 /* Simplification rules and other functions to implement log10 (base-10 logarithm)
   2  *
   3  * Copyright 2006 by Robert Dodier
   4  *
   5  * Released under the terms of the GNU General Public License
   6  *
   7  * Summary:
   8  *
   9  * log10 (<float>) => log (<float>) / log (10.0)
  10  * log10 (<bigfloat>) => log (<bigfloat>) / log (10b0)
  11  * log10 (a / b) => log10 (a) - log10 (b)
  12  * log10 (a * b) => log10 (a) + log10 (b)
  13  * log10 (a ^ b) => b log10 (a)
  14  * log10 (<integer divisible by 10>) => <power of 10 in integer> + log10 (<remainder>)
  15  * log10 (10) => 1
  16  * log10 (1) => 0
  17  * d/dx log10(x) => 1/log(10) 1/x
  18  *
  19  * log10 is an increasing function (Maxima can't do much with that, but anyway ...)
  20  *
  21  * log10 (x) => log (x) / log (10) (if you ask for it explicitly)
  22  *
  23  * Could also try to simplify is (log10 (x) < y) => is (x < 10^y),
  24  * likewise with other relational operators, but I didn't do that yet.
  25  *
  26  * Examples:
  27
  28 log10 (1);
  29 log10 (10);
  30 log10 (12/13);
  31 log10 (12/13), numer;
  32 log10 (2.3 * 234/2423);
  33 log10 (2.3 * 234/2423 + 1b0);
  34 log10 (1230000);
  35 log10 (9293923000000000 * w^23);
  36 log10 (w * (3 + b)^-8);
  37 log10 (123 * qwwqwe/235 * q^e);
  38 plot2d (log10 (x), [x, 1/100, 10]);
  39 quad_qags (log10 (x), x, 1, 10);
  40 integrate (log10 (x), x, 1, 10);
  41 expand_log10 (%);
  42 ''%, nouns;
  43 %, numer;
  44 assume (x > y);
  45 is (log10 (x) > log10 (y));
  46 diff (log10 (x), x);
  47 expand_log10 (sin (log10 (x) + log10 (y)) / log10 (z));
  48 contract_log10 (%);
  49
  50  * which yield these outputs:
  51
  52 (%i1) load ("./log10.mac");
  53 (%o1)                      ./log10.mac
  54 (%i2) log10 (1);
  55 (%o2)                           0
  56 (%i3) log10 (10);
  57 (%o3)                           1
  58 (%i4) log10 (12/13);
  59 (%o4)                 log10(12) - log10(13)
  60 (%i5) log10 (12/13), numer;
  61 (%o5)                 - .03476210625921192
  62 (%i6) log10 (2.3 * 234/2423);
  63 (%o6)                  - .6534097207097704
  64 (%i7) log10 (2.3 * 234/2423 + 1b0);
  65 Warning:  Float to bigfloat conversion of 0.22212133718530744
  66 (%o7)                 8.711432657265806b-2
  67 (%i8) log10 (1230000);
  68 (%o8)                    log10(123) + 4
  69 (%i9) log10 (9293923000000000 * w^23);
  70 (%o9)           23 log10(w) + log10(9293923) + 9
  71 (%i10) log10 (w * (3 + b)^-8);
  72 (%o10)              log10(w) - 8 log10(b + 3)
  73 (%i11) log10 (123 * qwwqwe/235 * q^e);
  74 (%o11) log10(qwwqwe) + e log10(q) - log10(235) + log10(123)
  75 (%i12) plot2d (log10 (x), [x, 1/100, 10]);
  76 (%o12)
  77 (%i13) quad_qags (log10 (x), x, 1, 10);
  78 (%o13)  [6.091349662870734, 1.985877489938315E-10, 63, 0]
  79 (%i14) integrate (log10 (x), x, 1, 10);
  80                           10
  81                          /
  82                          [
  83 (%o14)                   I   log10(x) dx
  84                          ]
  85                          /
  86                           1
  87 (%i15) expand_log10 (%);
  88                            10
  89                           /
  90                           [
  91                           I   log(x) dx
  92                           ]
  93                           /
  94                            1
  95 (%o15)                    -------------
  96                              log(10)
  97 (%i16) ''%, nouns;
  98 (%o16)         .4342944819032518 (10 log(10) - 9)
  99 (%i17) %, numer;
 100 (%o17)                  6.091349662870734
 101 (%i18) assume (x > y);
 102 (%o18)                       [x > y]
 103 (%i19) is (log10 (x) > log10 (y));
 104 (%o19)                        true
 105 (%i20) diff (log10 (x), x);
 106                                 1
 107 (%o20)                      ---------
 108                             log(10) x
 109 (%i21) expand_log10 (sin (log10 (x) + log10 (y)) / log10 (z));
 110                              log(y)    log(x)
 111                  log(10) sin(------- + -------)
 112                              log(10)   log(10)
 113 (%o21)           ------------------------------
 114                              log(z)
 115 (%i22) contract_log10 (%);
 116                     sin(log10(y) + log10(x))
 117 (%o22)              ------------------------
 118                             log10(z)
 119
 120  */
 121
 122 multp (e) := not atom(e) and op(e) = "*";
 123
 124 divp (e) := not atom(e) and op(e) = "/";
 125
 126 exponp (e) := not atom(e) and op(e) = "^";
 127
 128 divisible_by10 (e) := integerp(e) and e # 0 and mod (e, 10) = 0;
 129
 130 not10 (e) := e # 10;
 131
 132 power_of (m, n) := if n = 0 then 0 else block ([p : 0], while mod (n, m) = 0 do (p : p + 1, n : n/m), p);
 133
 134 block ([simp : false],
 135     local (aa, bb, cc, dd, ee, ff, gg, hh),
 136     matchdeclare (aa, multp, bb, exponp, cc, floatnump, dd, bfloatp, ee, divisible_by10, ff, all, gg, divp, hh, not10),
 137     tellsimp (log10(cc), log(cc)/log(10.0)),
 138     tellsimp (log10(dd), log(dd)/log(10b0)),
 139     tellsimp (log10(gg), log10 (first (args (gg))) - log10 (second (args (gg)))),
 140     tellsimp (log10(aa), apply ("+", map (log10, args (aa)))),
 141     tellsimp (log10(bb), second (args (bb)) * log10 (first (args (bb)))),
 142     tellsimp (log10(ee), block ([p : power_of (10, ee)], p + log10 (ee / 10^p))),
 143     tellsimp (log10(10), 1),
 144     tellsimp (log10(1), 0),
 145
 146     gradef (log10(x), (1/log(10))*(1/x)),
 147
 148     declare (log10, increasing),
 149
 150     defrule (expand_log10_rule, log10 (ff), log(ff)/log(10)),
 151     expand_log10 (expr) := apply1 (expr, expand_log10_rule),
 152
 153     defrule (contract_log10_rule, log (hh), log10(hh) * log(10)),
 154     contract_log10 (expr) := apply1 (expr, contract_log10_rule));
 155
 156