Logo Search packages:      
Sourcecode: mathomatic version File versions

matho-sumsq.c

/*
 * Find the sum of the squares for a given integer.
 *
 * Copyright (C) 2004 George Gesslein II.
 *
 * This program finds the minimum number of positive integers squared
 * and summed to equal a given positive integer.  If the number of squared
 * integers is 2 or less, all combinations with 2 squares are displayed,
 * otherwise only the first solution found is displayed.
 *
 * Usage:  matho-sumsq [numbers]
 *
 * If nothing is specified on the command line, the program starts at 0 and
 * counts up.
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>

#define     true  1
#define     false 0

long  squares[4];

main(argc, argv)
int   argc;
char  *argv[];
{
      int   i;
      long  d1;
      char  *cp;

      d1 = 0;
      if (argc > 1) {
            for (i = 1; i < argc; i++) {
                  d1 = strtol(argv[i], &cp, 10);
                  if (d1 <= 0) {
                        printf("Invalid number.\n");
                        exit(1);
                  }
                  if (d1 >= LONG_MAX) {
                        printf("Number too large (MAX = %ld).\n", LONG_MAX);
                        exit(1);
                  }
                  findsq(d1);
            }
      } else {
            for (;; d1 += 1) {
                  findsq(d1);
            }
      }
      exit(0);
}

/*
 * Return the truncated integer square root of "l" using longs.
 */
long
lsqrt(l)
long  l;
{
      long  x1, x2;
      long  testl;
      int   nbits;
      int   i;

      if (l <= 0L) {
            if (l != 0L) {
                  fprintf(stderr, "Domain error in lsqrt().\n");
            }
            return 0L;
      }
/* select a good starting value: */
      nbits = (sizeof(l) * 8) - 1;
      for (i = 4, testl = 16L;; i += 2, testl <<= 2) {
            if (i >= nbits || l <= testl) {
                  x1 = (1L << (i / 2));
                  break;
            }
      }
/* perform the famous square root approximation iteration: */
      for (;;) {
            x2 = (l / x1 + x1) / 2L;
            if (x1 <= x2)
                  break;
            x1 = x2;
      }
      return x1;
}

findsq(d1)
long  d1;
{
      if (sumsq(d1, 2))
            return;
      if (sumsq(d1, 3))
            return;
      if (!sumsq(d1, 4)) {
            printf("Whoops!  Can't find the sum of four squares that equal %ld.\n", d1);
            exit(1);
      }
}

/*
 * Find the sum of "n" squares that equal "d1" and display if found.
 * Return false if not found.
 */
int
sumsq(d1, n)
long  d1;
int   n;
{
      int   i, j;
      long  d2;
      long  save;
      int   found_one;

      found_one = false;
      d2 = d1;
      i = 0;
try_again:
      for (;;) {
            if (i == 2)
                  save = d2;
            squares[i] = lsqrt(d2);
            d2 -= (squares[i] * squares[i]);
            i++;
            if (i >= n) {
                  break;
            }
      }
      if (d2 == 0) {
            for (i = 0; i < n; i++) {
                  d2 += (squares[i] * squares[i]);
            }
            if (d2 != d1) {
                  printf("Bug found!\n");
                  exit(1);
            }
            for (i = 0; i < (n - 1); i++) {
                  if (squares[i] < squares[i+1]) {
                        goto skip_print;
                  }
            }
            found_one = true;
            printf("%ld = %ld^2", d1, squares[0]);
            for (i = 1; i < n; i++) {
                  if (squares[i] != 0)
                        printf(" + %ld^2", squares[i]);
            }
            printf("\n");
skip_print:
            if (n != 2) {
                  if (!found_one) {
                        printf("Found bug!\n");
                        exit(1);
                  }
                  return found_one;
            }
      }
      switch (n) {
      case 4:
            if (squares[2] > squares[n-1]) {
                  squares[2] -= 1;
                  d2 = save - (squares[2] * squares[2]);
                  i = 3;
                  goto try_again;
            }
      case 3:
            if (squares[1] > squares[n-1]) {
                  squares[1] -= 1;
                  d2 = d1 - (squares[0] * squares[0]) - (squares[1] * squares[1]);
                  i = 2;
                  goto try_again;
            }
      case 2:
            if (squares[0] > squares[n-1]) {
                  squares[0] -= 1;
                  d2 = d1 - (squares[0] * squares[0]);
                  i = 1;
                  goto try_again;
            }
      }
      return found_one;
}

Generated by  Doxygen 1.6.0   Back to index