最大公約数（除算禁止） DouKaku?

あなたが使っている言語で除算と剰余が使えなくなりました。

以下の条件のもと最大公約数を求めるプログラムを書いてください。

条件

除算および剰余の使用禁止
加算や乗算から除算・剰余を単純に定義することも禁止とする
ただし, ビットシフトが面倒な場合には引数を2で割った商を返す関数を実装しても構わない
多倍長演算をサポートすること（各言語のライブラリ状況を見たいので）
引数は2つの正整数と仮定して構わない
F_1=1, F_2=1のフィボナッチ数列で2000番目と1999番目の最大公約数を求めたときのループ回数を教えてください

引き算のみで。

Z80マシン語で遊んでいたころ、除算を使わないでもこれでgcdが計算できることに気づいて大発見をしたような気分になりました。後でよく考えたらユークリッドの互除法を非効率にやっているだけでした。

隣り合うフィボナッチ数の場合は実は減算法でもあまり効率が悪くなりません。常に差がひとつ前のフィボナッチ数になるわけですから。

gosh> (time (gcd-nodiv (fib 2000) (fib 1999)))
;(time (gcd-nodiv (fib 2000) (fib 1999)))
; real   0.002
; user   0.000
; sys    0.000
1

ループ回数は互除法でやった場合と同じになります。

gosh> (gcd-nodiv/cnt (fib 2000) (fib 1999) 1)
1
1998

実はこの引数の場合は、組み込みのgcd (互除法使用) の方が遅かったりします。bignumの除算が重いんじゃないかな。

gosh> (time (gcd (fib 2000) (fib 1999)))
;(time (gcd (fib 2000) (fib 1999)))
; real   0.014
; user   0.010
; sys    0.000
1

(define (gcd-nodiv m n)
  (cond [(= n m) n]
        [(or (= m 1) (= n 1)) 1]
        [(> m n) (gcd-nodiv n (- m n))]
        [else    (gcd-nodiv m (- n m))]))

;;
;; 以下は測定用
;;
;;ループ回数のカウント
(define (gcd-nodiv/cnt m n cnt)
  (cond [(= n m) (values n cnt)]
        [(or (= m 1) (= n 1)) (values 1 cnt)]
        [(> m n) (gcd-nodiv/cnt n (- m n) (+ cnt 1))]
        [else    (gcd-nodiv/cnt m (- n m) (+ cnt 1))]))

;;メモワイズ版fib
(define fib
  (let1 memo (make-hash-table 'eqv?)
    (lambda (n)
      (cond [(<= n 2) 1]
            [(hash-table-get memo n #f)]
            [else (let1 val (+ (fib (- n 1)) (fib (- n 2)))
                    (hash-table-put! memo n val)
                    val)]))))

ちょっと補足。

> ループ回数は互除法でやった場合と同じになります。 

というのはあくまで引数が隣り合うフィボナッチ数である場合の話です。

haskellにはもともと多倍長でも使えるgcdが用意されています

1	mygcd = gcd

出題者です.

>隣り合うフィボナッチ数の場合は実は減算法でもあまり効率が悪くなりません。常に差がひとつ前のフィボナッチ数になるわけですから。

全くもってその通りです. 1024bitと512bitの素数とかでやった方が良かったようでorz

この回答が破綻するのは一方が他方よりはるかに大きな場合ですね。m >> n で m = q*n + rのとき n >> r になるような組み合わせを続けてゆくと減算のみによる方法では悲惨なことに。

減算を繰り返して割り算の代わりにする方法だと差が大きいときに問題になる、ということで 「文字列の長さをみて小さい側を10**n倍する」 という荒技を使ってみました。

# Screencast-o-maticで撮りながら作ったのにアップロード中にブラウザが落ちて消滅…orz

def gcd(x, y):
    if x == y: return x
    if x < y: y, x = x, y

    lx = len(str(x))
    ly = len(str(y))
    diff = lx - ly

    z = abs(y * 10 ** diff - x)
    z2 = abs(y * 10 ** (diff + 1) - x)
    z = min(z, z2)
    if diff > 0:
        z2 = abs(y * 10 ** (diff - 1) - x)
        z = min(z, z2)

    # zはxとyと公約数を共有する
    if z < y:
        return gcd(y, z)
    return gcd(z, y)

なんでscreencastを取ろうと思ったのかに関しての参考文献

>例えばどう書く？orgの問題を解くまでのコーディングとかが公開されたら面白そう。

http://d.hatena.ne.jp/higepon/20071108/1194495170

>ブラウザだけでスクリーンキャストが手軽に作成できる『Screencast-O-Matic』 | 100SHIKI.COM

http://www.100shiki.com/archives/2007/08/screencastomatic.html

最初に断っておきますが、.NET Framework 2.0 に多倍長整数クラスは無いです。というよりまだリリースされてません。
反則技で Visual J# のライブラリから java.math.BigInteger を使いました。反則ついでに BigInteger#gcd を使用…… orz

ちなみに Microsoft.Scripting.Math.BigInteger みたいなのもあります。IronPython とかから使う用で、C# からも使えます。

using System;
using java.math;
static class Program {
    static void Main(string[] args) {
        BigInteger fib1999 = new BigInteger("2611005926183501768338670946829097324475555189114843467397273230483773870037923307730410719313972291638157639230613843870597997481070930648667960025707364078851859017098672504986584144842548768373271309551281830431960537091677315014266625027123872238011234749984205478230617988978500613170516952885123444971471854671812569739975450866912490650853945622130138277040986146312325044424769652148982077548213909414076005501");
        BigInteger fib2000 = new BigInteger("4224696333392304878706725602341482782579852840250681098010280137314308584370130707224123599639141511088446087538909603607640194711643596029271983312598737326253555802606991585915229492453904998722256795316982874482472992263901833716778060607011615497886719879858311468870876264597369086722884023654422295243347964480139515349562972087652656069529806499841977448720155612802665404554171717881930324025204312082516817125");
        Console.WriteLine("gcd({0}, {1}) = {2}", fib1999, fib2000, fib1999.gcd(fib2000));
    }
}

２数の差が大きい時もわりと早く収束するかな。結果[1, 3994]

require "mathn"

class Fib
  def initialize
    @fib_a=[1,1]
  end
  def [](n)
    (@fib_a.size..n).each{|i| @fib_a<<=@fib_a[i-2]+@fib_a[i-1]}  
    @fib_a[n-1]
  end
end

def gcd(n,m)
  c=0
  until(n==m)
    n,m=m,n if n<m
    return [1,c] if m==1
    m1=m
    while(m1<n)
      m1<<=1
      c+=1       
    end
    m1>>=1
    n-=m1
    c+=1
  end
  [n,c]
end

fib=Fib.new
p gcd(fib[2000],fib[1999])

既に #4790でも出てますが、Javaだと任意精度の整数値に対して使えるgcdがあるのでそれを使えます。

#4790 と同じなのもアレなので、フィボナッチ数列の計算処理もしてみました。

import java.math.BigInteger;
import java.util.HashMap;
import java.util.Map;

public class Sample106 {
    private final Map<Integer, BigInteger> fibMemo_ = new HashMap<Integer, BigInteger>();

    public BigInteger getFibValue(int index) {
        if (index == 1 || index == 2) return BigInteger.ONE;
        BigInteger fib = fibMemo_.get(index);
        if (fib == null) {
            fib = getFibValue(index - 2).add(getFibValue(index - 1));
            fibMemo_.put(index, fib);
        }
        return fib;
    }

    public static void main(String[] args) {
        Sample106 instance = new Sample106();
        BigInteger fib1999 = instance.getFibValue(1999);
        BigInteger fib2000 = instance.getFibValue(2000);
        System.out.println(fib2000.gcd(fib1999));
    }
}

gcd は組み込みで存在しますが、こんな式を基に書いてみました。呼び出し回数 1937 回。

gcd(x,x) = x
gcd(2x+1, 2y+1) = gcd(2x+1, x-y)
gcd(2x+1, 2y) = gcd(2x+1, y)
gcd(2x, 2y) = 2*gcd(x, y)

(defun *gcd (x y &optional (acc 1) (count 0))
  (incf count)
  (cond ((= x y)
         (values count (* acc x)))
        ((and (oddp x) (oddp y)) ; odd, odd
         (*gcd (min x y) (ash (abs (- x y)) -1) acc count))
        ((oddp x) ; odd, even
         (*gcd x (ash y -1) acc count))
        ((oddp y) ; even, odd
         (*gcd (ash x -1) y acc count))
        (t ; even, even
         (*gcd (ash x -1) (ash y -1) (ash acc 1) count))))

(defun fib (x)
  (loop with l = '(1 1) for i from 2 to x
    do (push (+ (car l) (cadr l)) l)
    finally (return (car l))))

(format t "~{Called ~D times, result is ~D~}"
        (multiple-value-list (*gcd (fib 2000) (fib 1999))))

書き忘れ。この方法だと各ステップで引数の少なくとも片方が半分以下になるので、ループ回数は引数の桁数に関して線形です。

最初の条件は、このコードだけで言えばクリアしています。

最後の条件への答えは、このコードだけで言えば0です。

多言語クックブックにこういうお題を載せる意図がよくわかりません。

1	GCD[Fibonacci@1999, Fibonacci@2000]

いいお題を思いついたらぜひご投稿ください(^o^)

回答は二種類。単純減算と正規表現利用バージョン。ただし後者の方はfib(2000)とかは出来ません。

Greediest Cruftiest Dan

#!/usr/local/bin/perl
use strict;
use warnings;
use bignum;
use Memoize;

no warnings 'recursion';

sub fib{
    my $n = shift;
    return 1 if $n <= 2;
    fib($n-2)+fib($n-1);
}
memoize 'fib';

sub gcd{
    my ($m, $n) = @_;
    return gcd($n, $m) if $m < $n;
    my $r = $m % $n;
    return $n unless $r;
    gcd($n, $r);
}

sub gcd_nodiv{
    my ($m, $n) = @_;
    return gcd($n, $m) if $m < $n;
    $m -= $n while $m > $n;
    return $n unless $m;
    gcd($n, $m);

}

sub gcd_regexp{
    my ($m, $n) = @_;
    return gcd($n, $m) if $m < $n;
    my $sm = 1 x $m;
    my $sn = 1 x $n;
    $sm =~ s/^$sn//g;
    return $n unless $sm;
    gcd($n, length $sm);
}

sub say { local $\="\n"; print @_ };
say gcd_regexp(fib(25), fib(24));
say gcd_nodiv(fib(1999),fib(2000));

Squeak Smalltalk で。

通常は無限倍精度整数に対応した #gcd: を使います。ちなみに、無限倍精度向けの Integer>>#gcd: および SmallInteger>>#gcd: に細工（グローバル変数 COUNT を用意して各々のループ内に COUNT := COUNT + 1 を挿入）をしてカウントをしてみたところ、fib1999 gcd: fib2000 では 1127 でした。

| gcdNoDiv gcdNoDivCount fibMemo fib fib1999 fib2000 |
gcdNoDiv := [:pair |
    [(pair sort; first) = pair second or: [pair last = 1]]
        whileFalse: [pair at: 2 incrementBy: pair first negated].
    pair last].

gcdNoDivCount := [:pair |
    | count |
    count := 0.
    [(pair sort; first) = pair second or: [pair last = 1]] whileFalse: [
        count := count + 1.
        pair at: 2 incrementBy: pair first negated].
    count].

fibMemo := Dictionary new.
fib := [:nth |
    fibMemo at: nth ifAbsentPut: [
        nth < 3 ifTrue: [1] ifFalse: [
            (fib copy fixTemps value: nth-1) + (fib copy fixTemps value: nth-2)]]].
fib1999 := fib copy fixTemps value: 1999.
fib2000 := fib copy fixTemps value: 2000.

fib1999 gcd: fib2000.   "=> 1 "  " 組み込みの GCD "
gcdNoDiv value: {fib1999. fib2000}.   "=> 1 "
gcdNoDivCount value: {fib1999. fib2000}.   "=> 1998 "

試しに動かしてみたら gcd(10, 1) でスタックオーバーフローしてしまいました。再帰呼び出し時に z==0 だと無限ループするみたいです。

組み込みの #gcd: での回数は 2889 に訂正します。改めて見てみたら、カウンタをループの最も内側に仕込めていませんでした。orz

数値解析系はこういうお題だと肩身が狭いですねｗ
お題や回答に対して色々なご意見があるのは当然で、
それが（ある程度）定量的に判断できる+/-評価システムは
うまく出来てるんだなぁと思いました。

1 2	library(gmp) gcd(fibnum(1999), fibnum(2000))

出題者です.

条件を緩めることにしました. 以下の2条件のうちお好きな方を選んでください.

条件(易)

組み込み及びライブラリにあるgcdは使用禁止
除算および剰余の使用禁止
加算や乗算から除算・剰余を単純に定義することも禁止とする
ただし, ビットシフトが面倒な場合には引数を2で割った商を返す関数を実装しても構わない
引数は1-10000の間の二つの整数として構わない

条件(難)

組み込み及びライブラリにあるgcdは使用禁止
除算および剰余の使用禁止
加算や乗算から除算・剰余を単純に定義することも禁止とする
ただし, ビットシフトが面倒な場合には引数を2で割った商を返す関数を実装しても構わない
引数は二つの正整数として構わない (多倍長演算もあるとなお良い)

補足

(fib 2000) => 4224696333392304878706725602341482782579852840250681098010280137314308584370130707224123599639141511088446087538909603607640194711643596029271983312598737326253555802606991585915229492453904998722256795316982874482472992263901833716778060607011615497886719879858311468870876264597369086722884023654422295243347964480139515349562972087652656069529806499841977448720155612802665404554171717881930324025204312082516817125
(fib 1999) => 2611005926183501768338670946829097324475555189114843467397273230483773870037923307730410719313972291638157639230613843870597997481070930648667960025707364078851859017098672504986584144842548768373271309551281830431960537091677315014266625027123872238011234749984205478230617988978500613170516952885123444971471854671812569739975450866912490650853945622130138277040986146312325044424769652148982077548213909414076005501

kozima さんと同じ Josef Stein のアルゴリズム．
ついでに，fibも同様の2進計算で高速計算．
Haskellにとっては呼び出しカウントのほうが面倒．
(とりあえずナイーブにやっちゃった．)
呼出回数は1936でした．

*Main> uncurry gcd' (fib 1999) 0
(1,1936)

import Prelude hiding (gcd)

fib 0 = (0,1)
fib n | even n = let (a,b) = fib (n-1)         in (b, a+b)             
      | odd  n = let (a,b) = fib (halve (n-1)) in (b*b+a*a, b*b+2*a*b)       

halve = (`div` 2)

gcd 1 _ = 1 
gcd a b | a == b = a
        | even a = if even b then 2*gcd (halve a) (halve b)           else gcd b (halve a)
        | odd  a = if odd  b then gcd (min a b) (abs (halve (a - b))) else gcd a (halve b)
        

gcd' 1 _ c = (1,c+1)
gcd' a b c | a == b = (a,c+1)
           | even a = if even b then fstapp (2*) (gcd' (halve a) (halve b) (c+1))
                      else gcd' a (halve b) (c+1)
           | odd  a = if odd  b then gcd' (min a b) (abs (halve (a - b))) (c+1)
                      else gcd' a (halve b) (c+1)

fstapp f (x,y) = (f x,y)

shiroさんと同じ方法で。

import java.math.BigInteger
import java.math.BigInteger.ONE

implicit def int2big(n:int) = new BigInteger(n.toString)
def gcd(m:BigInteger, n:BigInteger):BigInteger = m.compareTo(n) match {
  case 0 => n
  case x if m == ONE || n == ONE => ONE
  case 1 => gcd(n, m.subtract(n))
  case x => gcd(m, n.subtract(m))
}

出題者です. 『The Art of Computer Programming』や, 参考ページに挙げたサーベイに歴史が載っていました.

二進互除法を考えたのは, SteinとSilver and Terzianです. Knuth先生曰く一番古いのは以下だそうです.

前漢の中国の数学書『九章算術』の第1章「方田」の第6問より引用.

> (6) 方田: 又有九十一分之四十九。問約之得幾何？

> 答曰：十三分之七。

> 約分術曰：可半者半之，不可半者，副置分母子之數，以少減多，更相減損，求其等也。以等數約之。

see: R. P. Brent “Twenty years' analysis of the binary Euclidean algorithm”

バッチは多倍長演算をサポートしていないため、図らずも #4841の「条件(易)」を選択。
除算・余剰については、「除算・余剰を使わずに閏年」で回答したものを流用しました。

  e.g.
    C:>gcd.bat 1975 1225
    25

遅延環境変数展開を利用しているので、Windows NTでは動作しません。Windows XPで動作
を確認。

:: gcd.bat
@echo off
  setlocal
    set n=0

    echo %1%2|findstr /r "[^0-9]" >NUL 2>&1
    if %ERRORLEVEL% equ 0 (echo %~n0 [NUMBER] [NUMBER]) & goto :EOF

    call :gcd %1 %2 n
  endlocal & echo %n%
goto :EOF

:gcd
  setlocal enabledelayedexpansion
    set a=%1
    set b=%2
    set m=0
    set q=0

    if %a% lss %b% (
      set /a a=%b%-%a%
      set /a b-=!a!
      set /a a+=!b!
    )

    :loop_gcd
      call :mod %a% %b% m
      if %m% equ 0 (endlocal & set %3=%b%) & goto :EOF
      :: 剰余を除数の半分未満にして演算効率を向上
      call :div %b% 2 q
      if %q% lss %m% set /a m=%b%-%m%
      set a=%b%
      set b=%m%
    goto loop_gcd
  endlocal
goto :EOF

:mod
  setlocal
    set m=0
    set q=0

    call :div %1 %2 q

    set /a m=%2*%q%
  endlocal & set /a %3=%1-%m%
goto :EOF

:div
  setlocal enabledelayedexpansion
    set i=0
    set m=%1
    set m_=0
    set n=%2
    set n_=0
    set q=0

    if %n% gtr %m% (endlocal & set %3=%q%) & goto :EOF

    call :msb %m% m_
    call :msb %n% n_
    set /a i=%m_%-%n_%

    :loop_div
      if %i% lss 0 (endlocal & set %3=%q%) & goto :EOF
      set /a "q<<=1"
      set /a "m_=%m%>>%i%"
      if %m_% geq %n% (
        set /a "n_=%n%<<%i%"
        set /a m-=!n_!
        set /a q+=1
      )
      set /a i-=1
    goto loop_div
  endlocal
goto :EOF

:: 最上位ビットを取得(most significant bit)
:msb
  setlocal
    set i=%1
    set n=0

    :loop_msb
      set /a "i>>=1"
      set /a n+=1
    if %i% neq 0 goto loop_msb
  endlocal & set %2=%n%
goto :EOF

DouKaku?