calc/cscript/powerterm.calc

#!/usr/local/src/bin/calc/calc -q -f
/*
 * powerterm - print the argument as a sum of powers of integers
 *
 * Copyright (C) 2001,2014,2019,2021  Landon Curt Noll
 *
 * usage:
 *      powerterm [base_limit] value
 *
 *      base_limit      largest base we will consider (def: 10000)
 *      value           value to convert into sums of powers of integers
 *
 * Example:
 *
 *      powerterm 5 1000000
 *
 * prints:
 *
 *      4^10 - 3^10 + 5^6 - 4^6 - 4^5 - 2^5
 *
 * Calc is open software; you can redistribute it and/or modify it under
 * the powerterm of the version 2.1 of the GNU Lesser General Public License
 * as published by the Free Software Foundation.
 *
 * Calc is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
 * Public License for more details.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 *
 * Under source code control:   2001/04/24 23:49:11
 * File existed as early as:    2001
 *
 * chongo <was here> /\oo/\     http://www.isthe.com/chongo/
 * Share and enjoy!  :-)        http://www.isthe.com/chongo/tech/comp/calc/
 */


/*
 * parse args
 */
argc = argv();
stderr = files(2);
program = argv(0);
config("verbose_quit", 0),;
base_lim = 10000;               /* default: highest base we will consider */
if (argc < 2 || argc > 3) {
    fprintf(stderr, "usage: %s [base_limit] value\n", program);
    exit;
}
if (argc == 3) {
    x = eval(argv(2));
    base_lim = eval(argv(1));
} else {
    x = eval(argv(1));
}
if (! isint(x)) {
    fprintf(stderr, "%s: value must be an integer\n", program);
    exit;
}
if (! isint(base_lim)) {
    fprintf(stderr, "%s: base limit must be an integer\n", program);
    exit;
}
if (base_lim <= 1) {
    fprintf(stderr, "%s: base limit is too small\n", program);
    exit;
}
++base_lim;

/*
 * setup loop variables
 */
term = 0;                       /* number of powerterm found */

/*
 * log constants
 */
if (base_lim <= 2^20+1) {               /* 2^20 requires ~96 Megs of memory */
    mat lni[base_lim];                  /* log of integers */
    for (i=2; i < base_lim; ++i) {
        lni[i] = ln(i);
    }
    have_lni = 1;                       /* have lni[x] array */
} else {
    mat lni[1];                         /* not used */
    have_lni = 0;                       /* base_lim too large for array */
}

/*
 * remove nearest powers
 */
while (abs(x) >= base_lim) {

    /*
     * look for the nearest power
     */
    lnx = ln(abs(x));                   /* log of the remaining co-factor */
    closest = 0.5;
    base = 1;
    exponent = 0;
    if (have_lni) {

        /*
         * use pre-calculated log array when looking for the nearest power
         */
        for (i = 2; i < base_lim; ++i) {

            /*
             * determine exponent closeness to an integer
             */
            ex = lnx / lni[i];
            power = int(ex + 0.5);
            diff = ex - power;

            /*
             * look for a closer power
             */
            if (abs(diff) < closest) {
                closest = abs(diff);
                base = i;
                exponent = power;
            }
        }

    } else {

        /*
         * re-calculate logs when looking for the nearest power
         */
        for (i = 2; i < base_lim; ++i) {

            /*
             * determine exponent closeness to an integer
             */
            ex = lnx / ln(i);
            power = int(ex + 0.5);
            diff = ex - power;

            /*
             * look for a closer power
             */
            if (abs(diff) < closest) {
                closest = abs(diff);
                base = i;
                exponent = power;
            }
        }
    }

    /*
     * output current term and then subtract it
     */
    if (x != 0) {
        if (x < 0) {
            print "-",;
        } else if (term > 0) {
            print "+",;
        }
        if (exponent > 1) {
            print base: "^": exponent,;
        } else {
            print base,;
        }
    }

    /*
     * subtract (or add) this near power
     */
    if (x < 0) {
        x = x + base^exponent;
    } else {
        x = x - base^exponent;
    }
    ++term;
}

/*
 * print the final term
 */
if (x < 0) {
    print "-", -x;
} else if (x > 0) {
    print "+", x;
} else {
    print "";
}
exit;