c_cpp/etc/calc/cscript/powerterm.calc

   1 #!/usr/local/src/cmd/calc/calc -q -s -f
   2 /*
   3  * powerterm - print the argument as a sum of powers of integers
   4  *
   5  * usage:
   6  *      powerterm [base_limit] value
   7  *
   8  *      base_limit      largest base we will consider (def: 10000)
   9  *      value           value to convert into sums of powers of integers
  10  *
  11  * Copyright (C) 2001  Landon Curt Noll
  12  *
  13  * Calc is open software; you can redistribute it and/or modify it under
  14  * the powerterm of the version 2.1 of the GNU Lesser General Public License
  15  * as published by the Free Software Foundation.
  16  *
  17  * Calc is distributed in the hope that it will be useful, but WITHOUT
  18  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  19  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  20  * Public License for more details.
  21  *
  22  * A copy of version 2.1 of the GNU Lesser General Public License is
  23  * distributed with calc under the filename COPYING-LGPL.  You should have
  24  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  25  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  26  *
  27  * @(#) $Revision: 30.1 $
  28  * @(#) $Id: powerterm.calc,v 30.1 2007/03/16 11:12:11 chongo Exp $
  29  * @(#) $Source: /usr/local/src/bin/calc/cscript/RCS/powerterm.calc,v $
  30  *
  31  * Under source code control:   2001/04/24 23:49:11
  32  * File existed as early as:    2001
  33  *
  34  * chongo <was here> /\oo/\     http://www.isthe.com/chongo/
  35  * Share and enjoy!  :-)        http://www.isthe.com/chongo/tech/comp/calc/
  36  */
  37
  38
  39 /*
  40  * parse args
  41  */
  42 config("verbose_quit", 0),;
  43 base_lim = 10000;               /* default: highest base we will consider */
  44 if (argv() < 2 || argv() > 3) {
  45     fprintf(files(2), "usage: %s [base_limit] value\n", argv(0));
  46     exit;
  47 }
  48 if (argv() == 3) {
  49     x = eval(argv(2));
  50     base_lim = eval(argv(1));
  51 } else {
  52     x = eval(argv(1));
  53 }
  54 if (! isint(x)) {
  55     fprintf(files(2), "%s: value must be an integer\n");
  56     exit;
  57 }
  58 if (! isint(base_lim)) {
  59     fprintf(files(2), "%s: base limit must be an integer\n");
  60     exit;
  61 }
  62 if (base_lim <= 1) {
  63     fprintf(files(2), "%s: base limit is too small\n");
  64     exit;
  65 }
  66 ++base_lim;
  67
  68 /*
  69  * setup loop variables
  70  */
  71 term = 0;                       /* number of powerterm found */
  72
  73 /*
  74  * log constants
  75  */
  76 if (base_lim <= 2^20+1) {               /* 2^20 requires ~96 Megs of memory */
  77     mat lni[base_lim];                  /* log of integers */
  78     for (i=2; i < base_lim; ++i) {
  79         lni[i] = ln(i);
  80     }
  81     have_lni = 1;                       /* have lni[x] array */
  82 } else {
  83     mat lni[1];                         /* not used */
  84     have_lni = 0;                       /* base_lim too large for array */
  85 }
  86
  87 /*
  88  * remove nestest powers
  89  */
  90 while (abs(x) >= base_lim) {
  91
  92     /*
  93      * look for the nearest power
  94      */
  95     lnx = ln(abs(x));                   /* log of the remaining co-factor */
  96     closest = 0.5;
  97     base = 1;
  98     exponent = 0;
  99     if (have_lni) {
 100
 101         /*
 102          * use pre-calculated log array when looking for the nearest power
 103          */
 104         for (i = 2; i < base_lim; ++i) {
 105
 106             /*
 107              * determine exponent closeness to an integer
 108              */
 109             ex = lnx / lni[i];
 110             power = int(ex + 0.5);
 111             diff = ex - power;
 112
 113             /*
 114              * look for a closer power
 115              */
 116             if (abs(diff) < closest) {
 117                 closest = abs(diff);
 118                 base = i;
 119                 exponent = power;
 120             }
 121         }
 122
 123     } else {
 124
 125         /*
 126          * re-calculate logs when looking for the nearest power
 127          */
 128         for (i = 2; i < base_lim; ++i) {
 129
 130             /*
 131              * determine exponent closeness to an integer
 132              */
 133             ex = lnx / ln(i);
 134             power = int(ex + 0.5);
 135             diff = ex - power;
 136
 137             /*
 138              * look for a closer power
 139              */
 140             if (abs(diff) < closest) {
 141                 closest = abs(diff);
 142                 base = i;
 143                 exponent = power;
 144             }
 145         }
 146     }
 147
 148     /*
 149      * output current term and then subtract it
 150      */
 151     if (x != 0) {
 152         if (x < 0) {
 153             print "-",;
 154         } else if (term > 0) {
 155             print "+",;
 156         }
 157         if (exponent > 1) {
 158             print base: "^": exponent,;
 159         } else {
 160             print base,;
 161         }
 162     }
 163
 164     /*
 165      * subtract (or add) this near power
 166      */
 167     if (x < 0) {
 168         x = x + base^exponent;
 169     } else {
 170         x = x - base^exponent;
 171     }
 172     ++term;
 173 }
 174
 175 /*
 176  * print the final term
 177  */
 178 if (x < 0) {
 179     print "-", -x;
 180 } else if (x > 0) {
 181     print "+", x;
 182 } else {
 183     print "";
 184 }
 185 exit;