src/N.GCD.cpp

#include <symbolism/N.h>

namespace symbolism {
namespace ring {

void N::gcd_large_large(const N& a, const N& b)
{
    // Find the multiplicities of the factor two
    size_t a_twos = a.two_multiplicity_large();
    size_t b_twos = b.two_multiplicity_large();

    // Get the odd parts of the numbers
    // and provide copies for the destructive gmp algo's
    N a_odd, b_odd;
    a_odd.right_shift_large(a, a_twos);
    b_odd.right_shift_large(b, b_twos);
    size_t s;
    
    // TODO: Comparisson malfunctions
    if(b_odd > a_odd)
    {
        swap(a_odd, b_odd);
    }
    
    reserve(b_odd.size);
    s = mpn_gcd(limbs(), a_odd.limbs(), a_odd.size + 1, b_odd.limbs(), b_odd.size + 1);

    truncate(s - 1);
    left_shift(std::min(a_twos, b_twos));
}


void N::gcd_large_small(const N& a, const limb_t b)
{
    if(fast_false(b == 0))
    {
        set(a);
    }
    else
    {
        reserve(0);
        limb = mpn_gcd_1(a.data, a.size + 1, b);
    }
}

void N::gcd_small_small(const limb_t u, const limb_t v)
{
    reserve(0);
    if(u ==0 || v == 0)
    {
        limb = u | v;
    }
    else
    {
        limb_t a = u;
        limb_t b = v;
        
        size_t a_twos = __builtin_ctzll(a);
        size_t b_twos = __builtin_ctzll(b);
        size_t twos = std::min(a_twos, b_twos);
        a >>= a_twos;

        // From here on, a is always odd.
        do
        {
            b >>= __builtin_ctzll(b);
 
            // Now a and b are both odd, so diff(a, b) is even.
            // Let a = min(a, b), b = diff(a, b)/2. */
            if (a <= b)
            {
                b -= a;
            }
            else
            {
                limb_t diff = a - b;
                a = b;
                b = diff;
            }
            b >>= 1;
        }
        while (b != 0);
        limb = a << twos;
    }
}

void N::gcd_small_small(const limb_t b)
{
}

void N::gcd_large_large(const N& b)
{
}

void N::gcd_large_small(const limb_t b)
{
}

void N::gcd_small_large(const N& b)
{
}

}}

---

Copyright © 2007-2008 Remco Bloemen.

Generated on Tue Jan 22 17:35:31 2008 for symbolism by doxygen 1.5.4

Hosted by