Vedic Divisibility

July 8, 2011

We extract into separate functions the calculation of the last digit and the number remaining after stripping the last digit:

(define (last x) (car (reverse (digits x))))

(define (but-last x)
  (undigits (reverse (cdr (reverse (digits x))))))

For the calculation of the osculator we need to figure out the smallest multiple of the divisor that ends with 9. That’s easy. Since the divisor must be odd and not divisible by 5, it must end in 1, 3, 7 or 9. If it ends in 1, multiply by 9. If it ends in 3, multiply by 3. If it ends in 7, multiply by 7. And if it ends in 9, do nothing. Then the calculation of the osculator is simple:

(define (osculator n)
  (define (f x) (+ (but-last (* n x)) 1))
  (let ((ds (digits n)))
    (case (car (reverse ds))
     ((1) (f 9))
     ((3) (f 3))
     ((7) (f 7))
     ((9) (f 1))
    (else (error 'osculator "unsupported")))))

The divisibility test is recursive, reducing the size of the dividend n at each step:

(define (divisible? n d)
  (let ((osc (osculator d)))
    (let loop ((n n))
      (let ((next (+ (but-last n) (* last n) osc))))
        (if (< next n) (loop next) (zero? (modulo n d)))))))

Here are some examples:

> (divisible? 13174584 23)
#t
> (divisible? 175121 37)
#t
> (divisible? 134567 29)
#f

We used digits and undigits from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/qZrg7y54.

About these ads

Pages: 1 2

5 Responses to “Vedic Divisibility”

  1. arturasl said

    Solution in Python:

    def divides(dividend, divisor):
        endings = {1 : 9, 3 : 3, 7 : 7, 9 : 1}
    
        if divisor % 10 not in endings:
            raise Exception('divisor should end in 1, 3, 7 or 9')
    
        osculator = ((divisor * endings[divisor % 10]) // 10) + 1
    
        while True:
            prev_dividend = dividend
            dividend = dividend // 10 + (dividend % 10) * osculator
    
            if prev_dividend <= dividend:
                dividend = prev_dividend
                break
    
        return (dividend % divisor) == 0
    
    print(divides(13174584, 23))
    
  2. WillJ said
    #include <iostream>
    
    using namespace std;
    
    int findOsculator( int num )
    {
        int temp, i = 1;
        do
        {
            temp = num * i;
            i++;
        }
        while( temp % 10 != 9 );
        temp = temp / 10;
        return temp + 1;
    }
    
    bool divisible(int dividend, int divisor)
    {
        int remainder,prev_dividend, osc;
        osc = findOsculator( divisor );
        do
        {
            prev_dividend = dividend;
            remainder = dividend % 10;
            dividend = (dividend / 10) + (osc * remainder);
        }
        while( dividend < prev_dividend );
        return !( dividend % divisor );
    }
    
    int main()
    {
        cout << divisible(13174584,23);
    }
    
  3. Graham said

    Python:

    def is_divisible(n, d):
        osc = 1 + (d * {1: 9, 3: 3, 7: 7, 9: 1}[d % 10]) // 10
        old, new = float('inf'), n
        while new < old:
            a, b = divmod(new, 10)
            old, new = new, a + b * osc
        return new % d == 0
    
    if __name__ == "__main__":
        for (n, d) in [(13174584, 23), (175121, 37), (134567, 29)]:
            print "Is %d divisible by %d?%7s" % (n, d, "Yes." if is_divisible(n, d)
                                                 else "No.")
    

    And since Remco hasn’t jumped on it yet, here’s my attempt in Haskell:

    import Data.Map
    
    osculator :: (Integral a) => a -> a
    osculator d = f e where
              f x = 1 + d * x `div` 10
              e = fromList [(1,9), (3,3), (7,7), (9,1)] ! (d `mod` 10)
    
    isDivisible :: (Integral a) => a -> a -> Bool
    isDivisible n d = aux old new where
                old = n + 1
                new = n
                osc = osculator d
                aux x y | x <= y            = y `mod` d == 0
                        | otherwise         = aux y (a + b * osc) where (a,b) = divMod y 10
    
    main :: IO ()
    main = mapM_ (\(n,d) -> print $ "Is " ++ show n ++ " divisible by " ++ show d ++
         "? " ++ show (isDivisible n d)) [(13174584, 23), (175121, 37), (134567, 29)]
    
  4. Mike said

    Just to be different, how about an F# version.

    
    exception BadValue of string
    
    let osculator n = 
        match (n / 10, n % 10) with
        | d,1 -> d*9 + 1
        | d,3 -> d*3 + 1
        | d,7 -> d*7 + 1
        | d,9 -> d   + 1
        | _,_ -> raise <| BadValue("divisor cannot be divisible by 2 or 5")
    
    let isdivisible dividend divisor =
        let osc = osculator divisor
        let rec loop old_dividend =
            let new_dividend = old_dividend / 10 + old_dividend % 10 * osc
            if new_dividend >= old_dividend then
                old_dividend % divisor = 0
            else
                loop new_dividend
        loop dividend
    
  5. Mike said

    Line 07 above should be:

        | d,7 -> d*7 + 5   // 7*7 = 49, so carry the 4
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 621 other followers

%d bloggers like this: