Combinations calculations

Discussion in 'C Programming' started by CBFalconer, Mar 27, 2005.

  1. CBFalconer

    CBFalconer Guest

    I have, for my own amusement, written the following code. However
    I am in grave doubts as to the accuracy of the code in the function
    "variations", having spent 50 odd years forgetting everything about
    these calculations. Criticism and suggestions welcomed.

    Note that the solutions presented to variations are all ordered,
    but may have trailing zeroes.

    /* find integral solutions to x*x + y*y + z*z ... = r*r */
    /* Public domain, by C.B. Falconer, 2005-03-25 */

    #include <stdio.h>
    #include <limits.h>
    #include <stdlib.h>

    /* 181 allows operation with 16 bit ints */
    #define MAXRADIUS 181
    #define MAXDIMS 10

    #if (INT_MAX / MAXRADIUS < MAXRADIUS)
    # error "MAXRADIUS too big"
    #endif

    static unsigned int squares[MAXRADIUS+1];

    struct soln {
    int dims;
    int radius;
    int id;
    int sv[MAXDIMS+1]; /* solution vector */
    };

    /* -------------- */

    static void help(void)
    {
    printf("Usage: diophan radius [dimensions]\n"
    " with radius <= %d, dimensions in 1..%d\n",
    MAXRADIUS, MAXDIMS
    );
    exit(0);
    } /* help */

    /* -------------- */

    /* Problem - how to compute the independant solutions
    generated by negating points. 0 points can't be
    negated, and two identical points are trivially
    interchangeable, all of which reduce the soln count.
    The negations generate reflections about some axis.
    I don't remember any combinatorial maths.
    */
    /* NULL input returns (and resets) accumulated sum */
    static long int variations(struct soln *s)
    {
    long int vars;
    static long int vartotal;
    int d, f;

    vars = vartotal;
    if (!s) vartotal = 0;
    else {
    /* How many variations of s are available by reflections */
    /* useful test group for this are (rad, dims) =
    9, 3; 25, 3; 45, 3; 10, 4; 3, 9; 3, 10 */
    vars = 1;
    s->sv[0] = 0;
    for (d = 1, f = s->dims; d <= s->dims; d++, f--) {
    /* errors here on duplicated entries ?? */
    if (s->sv[d]) {
    /* Two possible values here, negatable */
    if (s->sv[d] != s->sv[d-1]) vars *= 2;
    vars *= f;
    }
    }
    vartotal += vars;
    }
    return vars;
    } /* variations */

    /* -------------- */

    static void dumpsolution(struct soln *s)
    {
    int d;

    if (s) { /* else just accumulated total vars */
    s->id++;
    printf("%5d:", s->id);
    for (d = 1; d <= s->dims; d++) {
    printf(" %3d", s->sv[d]);
    }
    }
    printf(" [%ld]\n", variations(s));
    } /* dumpsolution */

    /* -------------- */

    /* brute force exhaustive search */
    /* I have rediscovered the backtracking algorithm! :) */
    static void solve(int r, int dims)
    {
    int d; /* dimension index */
    int needed[MAXDIMS+1]; /* solution criterion */
    struct soln s; /* a solution develops here */

    s.dims = dims; s.radius = r; s.id = 0;
    for (d = 1; d <= dims; d++) s.sv[d] = 0;
    needed[0] = squares[r];
    s.sv[0] = 1; /* just to start the sequence */

    for (d = 1; ; d++) {
    for (s.sv[d] = s.sv[d-1]; s.sv[d] <= r; s.sv[d]++) {
    needed[d] = needed[d-1] - squares[s.sv[d]];
    if (needed[d] <= 0) { /* backtrack */
    if (0 == needed[d]) dumpsolution(&s);
    s.sv[d] = 0;
    d--; /* use prev. dim */
    if (d) continue;
    return; /* all done */
    }
    else if (d < dims) break; /* check next dimension */
    }
    }
    } /* solve */

    /* -------------- */

    int main(int argc, char **argv)
    {
    size_t temp;
    char *errp;
    unsigned int r, dims;

    r = 0; dims = 2;
    if ((argc < 2) || (argc > 3)) help();
    if (argc > 1) {
    temp = strtoul(argv[1], &errp, 10);
    if ((temp > MAXRADIUS) || (errp == argv[1])
    || (temp < 1)) help();
    else r = temp;
    }
    if (argc > 2) {
    temp = strtoul(argv[2], &errp, 10);
    if ((temp > MAXDIMS) || (errp == argv[2])
    || (temp < 1)) help();
    else dims = temp;
    }
    for (temp = 0; temp <= MAXRADIUS; temp++) {
    squares[temp] = temp * temp;
    }
    printf("solve %d %d\n", r, dims);
    solve(r, dims);
    dumpsolution(NULL); /* the accumulated total */
    return 0;
    } /* diophan main */

    --
    "If you want to post a followup via groups.google.com, don't use
    the broken "Reply" link at the bottom of the article. Click on
    "show options" at the top of the article, then click on the
    "Reply" at the bottom of the article headers." - Keith Thompson
    CBFalconer, Mar 27, 2005
    #1
    1. Advertising

  2. [A warning about follow-ups would have been nice. ;)
    Sent to clc as was original, F'up to comp.programming.]

    CBFalconer wrote:
    > I have, for my own amusement, written the following code.
    > However I am in grave doubts as to the accuracy of the code
    > in the function "variations", having spent 50 odd years
    > forgetting everything about these calculations. Criticism
    > and suggestions welcomed.
    >
    > Note that the solutions presented to variations are all ordered,
    > but may have trailing zeroes.
    >
    > /* find integral solutions to x*x + y*y + z*z ... = r*r */
    > /* Public domain, by C.B. Falconer, 2005-03-25 */
    > ...


    It's arguable whether 0 ordinates should count as a solution,
    but that's your call.

    I don't like the way you've used hidden static variables, but
    again, I'll put that aside.

    You _are_ over-reporting the number of variations.

    > ...
    > /* Problem - how to compute the independant solutions
    > generated by negating points. 0 points can't be
    > negated, and two identical points are trivially
    > interchangeable, all of which reduce the soln count.
    > The negations generate reflections about some axis.
    > I don't remember any combinatorial maths.
    > */
    > /* NULL input returns (and resets) accumulated sum */
    > static long int variations(struct soln *s)
    > {
    > long int vars;
    > static long int vartotal;
    > int d, f;
    >
    > vars = vartotal;
    > if (!s) vartotal = 0;
    > else {
    > /* How many variations of s are available by reflections */
    > /* useful test group for this are (rad, dims) =
    > 9, 3; 25, 3; 45, 3; 10, 4; 3, 9; 3, 10 */
    > vars = 1;
    > s->sv[0] = 0;
    > for (d = 1, f = s->dims; d <= s->dims; d++, f--) {
    > /* errors here on duplicated entries ?? */
    > if (s->sv[d]) {
    > /* Two possible values here, negatable */
    > if (s->sv[d] != s->sv[d-1]) vars *= 2;
    > vars *= f;
    > }
    > }
    > vartotal += vars;
    > }
    > return vars;
    > } /* variations */


    % diophan 3 10
    solve 3 10
    1: 1 1 1 1 1 1 1 1 1 0 [7257600]
    2: 1 1 1 1 1 2 0 0 0 0 [604800]
    3: 1 2 2 0 0 0 0 0 0 0 [2880]
    4: 3 0 0 0 0 0 0 0 0 0 [20]
    [7865300]

    The variations need to be a factored by 2 for each non-zero
    dimension (+ or -). To calculate the combinations, you need
    to group and count the number of similar ordinates.

    For example, the variations for the solutions above should be...

    1: 9x1 1x0 [10!/9!1! * 2^9]
    2: 5x1 1x2 4x0 [10!/5!1!4! * 2^6]
    3: 1x1 2x2 7x0 [10!/1!2!7! * 2^3]
    4: 1x3 9x0 [10!/1!9! * 2^1]

    Try...

    static long int variations(struct soln *s)
    {
    long int vars;
    static long int vartotal;
    int d, f;

    vars = vartotal;
    if (!s) vartotal = 0;
    else {
    for (vars = 1, f = 0, d = s->dims; d; d--) {
    if (s->sv[d]) vars *= 2;
    if (d != s->dims && s->sv[d] != s->sv[d+1]) f = 0;
    vars = vars * d / ++f;
    }
    vartotal += vars;
    }
    return vars;
    }

    The calculation of combinations (for 2: above) is done as...

    2: 1 1 1 1 1 2 0 0 0 0

    10! - 1 2 3 4 5 6 7 8 9 10
    5!1!4! - 5 4 3 2 1 1 4 3 2 1

    For no particular reason, I start from the right and work left.
    I used f to count the number of ordinate repeats, reseting it
    when a different ordinate is encountered.

    With the change, I get...

    % diophan 3 10
    solve 3 10
    1: 1 1 1 1 1 1 1 1 1 0 [5120]
    2: 1 1 1 1 1 2 0 0 0 0 [80640]
    3: 1 2 2 0 0 0 0 0 0 0 [2880]
    4: 3 0 0 0 0 0 0 0 0 0 [20]
    [88660]

    --
    Peter
    Peter Nilsson, Mar 28, 2005
    #2
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Ed Bangle
    Replies:
    0
    Views:
    317
    Ed Bangle
    Nov 23, 2003
  2. Leszek
    Replies:
    4
    Views:
    4,084
    Leszek
    Dec 2, 2003
  3. =?Utf-8?B?cHBhdGVs?=

    shipping and handling calculations ?

    =?Utf-8?B?cHBhdGVs?=, May 4, 2004, in forum: ASP .Net
    Replies:
    0
    Views:
    315
    =?Utf-8?B?cHBhdGVs?=
    May 4, 2004
  4. Jon Maz
    Replies:
    4
    Views:
    1,674
    Jon Maz
    Nov 12, 2004
  5. Rick DeBay

    Currency calculations

    Rick DeBay, Apr 12, 2004, in forum: Java
    Replies:
    10
    Views:
    3,559
    Mark Thornton
    Apr 13, 2004
Loading...

Share This Page