source: trunk/src/helpers/math.c@ 108

/*
 *@@sourcefile math.c:
 *      some math helpers.
 *
 *      This file is new with V0.9.14 (2001-07-07) [umoeller]
 *      Unless marked otherwise, these things are based
 *      on public domain code found at
 *      "ftp://ftp.cdrom.com/pub/algorithms/c/".
 *
 *      Usage: All C programs.
 *
 *      Function prefix:
 *
 *      --  math*: semaphore helpers.
 *
 *@@added V0.9.14 (2001-07-07) [umoeller]
 *@@header "helpers\math.h"
 */

#define OS2EMX_PLAIN_CHAR
    // this is needed for "os2emx.h"; if this is defined,
    // emx will define PSZ as _signed_ char, otherwise
    // as unsigned char

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

#include "setup.h"                      // code generation and debugging options

#include "helpers\math.h"
#include "helpers\linklist.h"
#include "helpers\standards.h"

#pragma hdrstop

/*
 *@@category: Helpers\Math helpers
 *      see math.c.
 */


/*
 *@@ mathGCD:
 *      returns the greatest common denominator that
 *      evenly divides m and n.
 *
 *      For example, mathGCD(10, 12) would return 2.
 *
 *      Modern Euclidian algorithm (The Art of Computer
 *      Programming, 4.5.2).
 */

int mathGCD(int a, int b)
{
    int r;

    while (b)
    {
        r = a % b;
        a = b;
        b = r;
    }

    return (a);

}

/* binary GCD alg. Based on following facts.
   If N and M are even, then gcd(N,M) = 2 * gcd (N/2, M/2)
   If N even and M odd, then gcd(N,M) = gcd(N/2, M)
   If N,M odd, then |N-M|<max(N,M) thus guaranteeing that
             we can use gcd(N,M) = gcd(N-M,M) to reduce.
   fact: a&1 is 1 iff a is odd. Mult. and Div. by 2 more efficiently
   implemented as shift-left 1 and shift-right 1 resp.

    Found at http://remus.rutgers.edu/~rhoads/Code/int_gcd.c
*/

/* long gcd( long a, long b)
{
  int t;

  if (a < 0) a = -a;
  if (!b) return a;
  if (b < 0) b = -b;
  if (!a) return b;
  t = 0;
  while (! ((a|b) & 1)) {a >>= 1; b >>= 1; ++t;}
  while (! (a&1)) a >>= 1;
  while (! (b&1)) b >>= 1;
  while (a != b)
     {
     if (a > b)
         {
         a -= b;
         do {a >>= 1;} while (! (a&1));
         }
     else {
         b -= a;
         do {b >>= 1;} while (! (b&1));
         }
     }
  return(a<<t); */

/*
 *@@ mathIsSquare:
 *      returns 1 if x is a perfect square, that is, if
 *
 +      sqrt(x) ^ 2 ==x
 */

int mathIsSquare(int x)
{
    int t;
    int z = x % 4849845; // 4849845 = 3*5*7*11*13*17*19
    double r;
    // do some quick tests on x to see if we can quickly
    // eliminate it as a square using quadratic-residues.
    if (z % 3 == 2)
        return 0;

    t =  z % 5;
    if((t==2) || (t==3))
        return 0;
    t =  z % 7;
    if((t==3) || (t==5) || (t==6))
        return 0;
    t = z % 11;
    if((t==2) || (t==6) || (t==7) || (t==8) || (t==10))
        return 0;
    t = z % 13;
    if((t==2) || (t==5) || (t==6) || (t==7) || (t==8) || (t==11))
        return 0;
    t = z % 17;
    if((t==3) || (t==5) || (t==6) || (t==7) || (t==10) || (t==11) || (t==12) || (t==14))
        return 0;
    t = z % 19;
    if((t==2) || (t==3) || (t==8) || (t==10) || (t==12) || (t==13) || (t==14) || (t==15) || (t==18))
        return 0;

    // If we get here, we'll have to just try taking
    // square-root & compare its square...
    r = sqrt(abs(x));
    if (r*r == abs(x))
        return 1;

    return 0;
}

/*
 *@@ mathFindFactor:
 *      returns the smallest factor > 1 of n or 1 if n is prime.
 *
 *      From "http://tph.tuwien.ac.at/~oemer/doc/quprog/node28.html".
 */

int mathFindFactor(int n)
{
    int i,
        max;
    if (n <= 0)
        return 0;

    max = floor(sqrt(n));

    for (i=2;
         i <= max;
         i++)
    {
        if (n % i == 0)
            return i;
    }
    return 1;
}

/*
 *@@ testprime:
 *      returns 1 if n is a prime number.
 *
 *      From "http://tph.tuwien.ac.at/~oemer/doc/quprog/node28.html".
 */

int mathIsPrime(unsigned n)
{
    int i,
        max = floor(sqrt(n));

    if (n <= 1)
        return 0;

    for (i=2;
         i <= max;
         i++)
    {
        if (n % i == 0)
            return 0;
    }

    return 1;
}

/*
 *@@ mathFactorBrute:
 *      calls pfnCallback with every integer that
 *      evenly divides n ("factor").
 *
 *      pfnCallback must be declared as:
 *
 +          int pfnCallback(int iFactor,
 +                          int iPower,
 +                          void *pUser);
 *
 *      pfnCallback will receive the factor as its
 *      first parameter and pUser as its third.
 *      The second parameter will always be 1.
 *
 *      The factor will not necessarily be prime,
 *      and it will never be 1 nor n.
 *
 *      If the callback returns something != 0,
 *      computation is stopped.
 *
 *      Returns the no. of factors found or 0 if
 *      n is prime.
 *
 *      Example: mathFactor(42) will call the
 *      callback with
 *
 +          2 3 6 7 14 21
 *
 +      This func is terribly slow.
 */

int mathFactorBrute(int n,                           // in: integer to factor
                    PFNFACTORCALLBACK pfnCallback,   // in: callback func
                    void *pUser)                     // in: user param for callback
{
    int rc = 0;

    if (n > 2)
    {
        int i;
        int max = n / 2;

        for (i = 2;
             i <= max;
             i++)
        {
            if ((n % i) == 0)
            {
                rc++;
                // call callback with i as the factor
                if (pfnCallback(i,
                                1,
                                pUser))
                    // stop then
                    break;
            }
        }
    }

    return (rc);
}

/*
 *@@ mathFactorPrime:
 *      calls pfnCallback for every prime factor that
 *      evenly divides n.
 *
 *      pfnCallback must be declared as:
 *
 +          int pfnCallback(int iFactor,
 +                          int iPower,
 +                          void *pUser);
 *
 *      pfnCallback will receive the prime as its
 *      first parameter, its power as its second,
 *      and pUser as its third.
 *
 *      For example, with n = 200, pfnCallback will
 *      be called twice:
 *
 *      1)  iFactor = 2, iPower = 3
 *
 *      2)  iFactor = 5, iPower = 2
 *
 *      because 2^3 * 5^2 is 200.
 *
 *      Returns the no. of times that the callback
 *      was called. This will be the number of
 *      prime factors found or 0 if n is prime
 *      itself.
 */

int mathFactorPrime(double n,
                    PFNFACTORCALLBACK pfnCallback,   // in: callback func
                    void *pUser)                     // in: user param for callback
{
    int rc = 0;

    double d;
    double k;

    if (n <= 3)
       // this is a prime for sure
       return 0;

    d = 2;
    for (k = 0;
         fmod(n, d) == 0;
         k++)
    {
        n /= d;
    }

    if (k)
    {
        rc++;
        pfnCallback(d,
                    k,
                    pUser);
    }

    for (d = 3;
         d * d <= n;
         d += 2)
    {
        for (k = 0;
             fmod(n,d) == 0.0;
             k++ )
        {
            n /= d;
        }

        if (k)
        {
            rc++;
            pfnCallback(d,
                        k,
                        pUser);
        }
    }

    if (n > 1)
    {
        if (!rc)
            return (0);

        rc++;
        pfnCallback(n,
                    1,
                    pUser);
    }

    return (rc);
}

/*
 *@@ mathGCDMulti:
 *      finds the greatest common divisor for a
 *      whole array of integers.
 *
 *      For example, if you pass in three integers
 *      1000, 1500, and 1800, this would return 100.
 *      If you only pass in 1000 and 1500, you'd
 *      get 500.
 *
 *      Use the fact that:
 *
 *         gcd(u1, u2, ..., un) = gcd(u1, gcd(u2, ..., un))
 *
 *      Greatest common divisor of n integers (The
 *      Art of Computer Programming, 4.5.2).
 */

int mathGCDMulti(int *paNs,             // in: array of integers
                 int cNs)               // in: array item count (NOT array size)
{
    int d = paNs[cNs-1];
    int k = cNs-1;

    while (    (d != 1)
            && (k > 0)
          )
    {
        d = mathGCD(paNs[k-1], d);
        k--;
    }

    return (d);
}

/* typedef struct _PRIMEDATA
{
    LINKLIST    llPrimes;
    int         iCurrentInt;
} PRIMEDATA, *PPRIMEDATA;


typedef struct _PRIMEENTRY
{
    int         iPrime;
    int         iLowestPowerFound;
                // lowest power that was found for this prime number;
                // if 0, a prime was previously present and then not
                // for a later prime
    int         iLastInt;
} PRIMEENTRY, *PPRIMEENTRY; */

/*
 *@@ GCDMultiCallbackCreate:
 *      first callback for mathGCDMulti.
 */

/* int XWPENTRY GCDMultiCallbackCreate(int n,             // in: prime
                                    int p,             // in: power
                                    void *pUser)       // in: user param, here root of tree
{
    // see if we had this
    PPRIMEDATA pData = (PPRIMEDATA)pUser;
    PLISTNODE pNode;
    PPRIMEENTRY pEntry;
    for (pNode = lstQueryFirstNode(&pData->llPrimes);
         pNode;
         pNode = pNode->pNext)
    {
        pEntry = (PPRIMEENTRY)pNode->pItemData;
        if (pEntry->iPrime == n)
        {
            // mark this as processed so we can later see
            // if the current integer had this prime factor
            // at all
            pEntry->iLastInt = pData->iCurrentInt;

            // printf(" %d^%d", n, p);

            // and stop
            return 0;
        }
    }

    // no entry for this yet:
    // printf(" %d^%d", n, p);
    pEntry = NEW(PRIMEENTRY);
    pEntry->iPrime = n;
    pEntry->iLowestPowerFound = p;
    pEntry->iLastInt = pData->iCurrentInt;
    lstAppendItem(&pData->llPrimes, pEntry);

    return (0);
} */

/*
 *@@ GCDMultiCallbackLowest:
 *      second callback for mathGCDMulti.
 */

/* int XWPENTRY GCDMultiCallbackLowest(int n,             // in: prime
                                    int p,             // in: power
                                    void *pUser)       // in: user param, here root of tree
{
    // see if we had this
    PPRIMEDATA pData = (PPRIMEDATA)pUser;
    PLISTNODE pNode;
    PPRIMEENTRY pEntry;
    for (pNode = lstQueryFirstNode(&pData->llPrimes);
         pNode;
         pNode = pNode->pNext)
    {
        pEntry = (PPRIMEENTRY)pNode->pItemData;
        if (pEntry->iPrime == n)
        {
            if (p < pEntry->iLowestPowerFound)
            {
                pEntry->iLowestPowerFound = p;
                // printf(" %d^%d", n, p);
            }
            pEntry->iLastInt = pData->iCurrentInt;

            // and stop
            return 0;
        }
    }

    exit(1);

    return (0);
} */

/* int mathGCDMulti(int *paNs,             // in: array of integers
                 int cNs)               // in: array item count (NOT array size)
{
    int i;
    double dGCD;

    PLISTNODE pNode;
    PPRIMEENTRY pEntry;

    PRIMEDATA Data;
    lstInit(&Data.llPrimes, TRUE);

    // this is done by prime-factoring each frequency
    // and then multiplying all primes with the lowest
    // common power that we found:

    // Example 1: If we have two timers with freq.
    // 1000 and 1500, obviously, the master timer
    // should run at 500 ms.

    // With prime factorization, we end up like this:

    // 1000 = 2^3 *       5^3
    // 1500 = 2^2 * 3^1 * 5^3

    // so the highest power for factor 2 is 2;
    //    the highest power for factor 3 is 0;
    //    the highest power for factor 5 is 3;

    // freq = 2^2 *       5^3 = 500

    // Example 2: If we have three timers with freq.
    // 1000 and 1500 and 1800:

    // 1000 = 2^3 *       5^3
    // 1500 = 2^2 * 3^1 * 5^3
    // 1800 = 2^3 * 3^2 * 5^2

    // so the highest power for factor 2 is 2;
    //    the highest power for factor 3 is 0;
    //    the highest power for factor 5 is 2;

    // freq = 2^2 *     * 5^2 = 100

    // Example 3: If we have four timers with freq.
    // 1200 and 1500 and 1800 and 60:

    // 1200 = 2^4 * 3^1 * 5^2
    // 1500 = 2^2 * 3^1 * 5^3
    // 1800 = 2^3 * 3^2 * 5^2
    //   60 = 2^2 * 3^1 * 5^1

    // so the highest power for factor 2 is 2;
    //    the highest power for factor 3 is 1;
    //    the highest power for factor 5 is 1;

    // freq = 2^2 * 3^1 * 5^1 = 60

    // go thru the ints array to create
    // an entry for each prime there is
    for (i = 0;
         i < cNs;
         i++)
    {
        // printf("\nFactoring %d\n", paNs[i]);
        Data.iCurrentInt = paNs[i];
        mathFactorPrime(paNs[i],
                        GCDMultiCallbackCreate,
                        &Data);
    }

    // now run this a second time to find the
    // lowest prime for each
    for (i = 0;
         i < cNs;
         i++)
    {
        // printf("\nTouching %d\n", paNs[i]);
        Data.iCurrentInt = paNs[i];
        mathFactorPrime(paNs[i],
                        GCDMultiCallbackLowest,
                        &Data);

        // all list items that were just touched
        // now have iLastInt set to the current
        // integer; so go thru the list and set
        // all items which do NOT have that set
        // to power 0 because we must not use them
        // in factoring
        for (pNode = lstQueryFirstNode(&Data.llPrimes);
             pNode;
             pNode = pNode->pNext)
        {
            pEntry = (PPRIMEENTRY)pNode->pItemData;
            if (pEntry->iLastInt != paNs[i])
            {
                pEntry->iLowestPowerFound = 0;
                // printf(" ###%d^%d", pEntry->iPrime, 0);
            }
        }
    }

    // printf("\nNew:\n");

    // now compose the GCD
    dGCD = 1;
    for (pNode = lstQueryFirstNode(&Data.llPrimes);
         pNode;
         pNode = pNode->pNext)
    {
        pEntry = (PPRIMEENTRY)pNode->pItemData;

        // printf(" %d^%d", pEntry->iPrime, pEntry->iLowestPowerFound);

        if (pEntry->iLowestPowerFound)
            dGCD *= pow(pEntry->iPrime, pEntry->iLowestPowerFound);
    }

    lstClear(&Data.llPrimes);

    return dGCD;
} */

// testcase

#ifdef BUILD_MAIN

#define INCL_DOSMISC
#include <os2.h>

int callback(int n,
             int p,
             void *pUser)
{
    if (p > 1)
        printf("%d^%d ", n, p);
    else
        printf("%d ", n);
    fflush(stdout);

    return (0);
}

int main (int argc, char *argv[])
{
    if (argc > 1)
    {
        int i,
            c,
            cInts = argc - 1;
        ULONG ulms, ulms2;
        int *aInts = malloc(sizeof(int) * cInts);

        for (c = 0;
             c < cInts;
             c++)
        {
            aInts[c] = atoi(argv[c + 1]);
        }

        DosQuerySysInfo(QSV_MS_COUNT,
                        QSV_MS_COUNT,
                        &ulms,
                        sizeof(ulms));

        c = mathGCDMulti(aInts,
                         cInts);

        DosQuerySysInfo(QSV_MS_COUNT,
                        QSV_MS_COUNT,
                        &ulms2,
                        sizeof(ulms2));

        printf("\nGCD is %d, %d ms time.\n",
               c,
               ulms2 - ulms);

        for (i = 0;
             i < cInts;
             i++)
        {
            printf("  %d / %d = %d rem. %d\n",
                   aInts[i], c, aInts[i] / c, aInts[i] % c);
        }

        /* c = mathFactorBrute(atoi(argv[1]),
                            callback,
                            0);
        printf("\n%d factors found.\n", c); */
    }

    return (0);
}

#endif