Changeset 771

Show
Ignore:
Timestamp:
06/30/08 10:47:37 (2 months ago)
Author:
Don Clugston
Message:

Restored bigint (fixed so that it compiles). I've updated my asm code so that it uses dynamic arrays throughout, instead of a mix of pointers and arrays. It's possible that pointers-only would be better. None of the asm code is hooked up yet.

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • trunk/phobos/std/bigint.d

    r746 r771  
    1 // This is just a placeholder 
    2 // bigint coming soon 
    3  
    4 module std.bigint; 
     1// Written in the D programming language. 
     2 
     3/** 
     4Provides a BigInt struct for multiprecision integer arithmetic. 
     5 
     6The internal representation is binary, not decimal. 
     7 
     8All relevant operators are overloaded. 
     9 
     10Example: 
     11---------------------------------------------------- 
     12        BigInt a = "9588669891916142"; 
     13        BigInt b = "7452469135154800"; 
     14        auto c = a * b; 
     15        assert(c == "71459266416693160362545788781600"); 
     16        auto d = b * a; 
     17        assert(d == "71459266416693160362545788781600"); 
     18        assert(d == c); 
     19        d = c * "794628672112"; 
     20        assert(d == "56783581982794522489042432639320434378739200"); 
     21        auto e = c + d; 
     22        assert(e == "56783581982865981755459125799682980167520800"); 
     23        auto f = d + c; 
     24        assert(f == e); 
     25        auto g = f - c; 
     26        assert(g == d); 
     27        g = f - d; 
     28        assert(g == c); 
     29        e = 12345678; 
     30        g = c + e; 
     31        auto h = g / b; 
     32        auto i = g % b; 
     33        assert(h == a); 
     34        assert(i == e); 
     35---------------------------------------------------- 
     36 
     37Authors: Janice Caron 
     38 
     39Date: 2008.05.18 
     40 
     41License: Public Domain 
     42 
     43 
     44Macros: 
     45    WIKI=Phobos/StdBigint 
     46*/ 
     47 
     48module bigint; 
     49import std.string       : format; 
     50import std.stdio        : writef, writefln; 
     51import std.algorithm    : min, max, swap; 
     52import std.traits       : isIntegral; 
     53 
     54alias uint Digit; /// alias for uint 
     55 
     56/****************** 
     57 * Struct representing a multiprecision integer 
     58 */ 
     59struct BigInt 
     60{ 
     61    Digits digits = [ cast(Digit)0 ]; 
     62 
     63    static const BigInt ZERO = { [ cast(Digit)0 ] }; 
     64    static const BigInt ONE = { [ cast(Digit)1 ] }; 
     65 
     66    /// 
     67    void opAssign(const BigInt n) 
     68    { 
     69        digits = n.digits; 
     70    } 
     71 
     72    /// 
     73    void opAssign(int n) 
     74    { 
     75        digits = cast(Digits)( [ cast(Digit)n ] ); 
     76    } 
     77 
     78    /// 
     79    void opAssign(uint n) 
     80    { 
     81        static if(BIG_ENDIAN) { auto a = [ cast(Digit)0, n ]; } 
     82        else                  { auto a = [ cast(Digit)n, 0 ]; } 
     83        Big b = bigInt(a); 
     84        digits = b.digits; 
     85    } 
     86 
     87    /// 
     88    void opAssign(long n) 
     89    { 
     90        static if(BIG_ENDIAN) { auto a = [ cast(Digit)(n>>32), cast(Digit)n ]; } 
     91        else                  { auto a = [ cast(Digit)n, cast(Digit)(n>>32) ]; } 
     92        Big b = bigInt(a); 
     93        digits = b.digits; 
     94    } 
     95 
     96    /// 
     97    void opAssign(ulong n) 
     98    { 
     99        static if(BIG_ENDIAN) { auto a = [ cast(Digit)0, cast(Digit)(n>>32), cast(Digit)n ]; } 
     100        else                  { auto a = [ cast(Digit)n, cast(Digit)(n>>32), cast(Digit)0 ]; } 
     101        Big b = bigInt(a); 
     102        digits = b.digits; 
     103    } 
     104 
     105    /// 
     106    void opAssign(string s) 
     107    { 
     108        Big b = fromString(s); 
     109        digits = b.digits; 
     110    } 
     111 
     112    /// 
     113    static BigInt opCall(T)(T n) 
     114    { 
     115        BigInt r; 
     116        r.opAssign(n); 
     117        return r; 
     118    } 
     119 
     120    // Convert TO other types 
     121 
     122    /// 
     123    void castTo(out BigInt r) const 
     124    { 
     125        r.digits = digits; 
     126    } 
     127 
     128    /// 
     129    void castTo(out int r) const 
     130    { 
     131        r = cast(int)tail(digits,1u)[0]; 
     132    } 
     133 
     134    /// 
     135    void castTo(out uint r) const 
     136    { 
     137        r = cast(uint)tail(digits,1u)[0]; 
     138    } 
     139 
     140    /// 
     141    void castTo(out long r) const 
     142    { 
     143        ulong t; 
     144        castTo(t); 
     145        r = cast(long)t; 
     146    } 
     147 
     148    /// 
     149    void castTo(out ulong r) const 
     150    { 
     151        mixin(setUp("x","*this")); 
     152        r = peek(xp); 
     153        xp = next(xp); 
     154        if (xp != xe) r += cast(ulong)(peek(xp)) << 32; 
     155    } 
     156 
     157    /// 
     158    void castTo(out string r) const 
     159    { 
     160        r = decimal(*this); 
     161    } 
     162 
     163    // Unary operator overloads 
     164 
     165    /// 
     166    BigInt opPos() const 
     167    { 
     168        BigInt r; 
     169        r.digits = digits; 
     170        return r; 
     171    } 
     172 
     173    /// 
     174    BigInt opNeg() const 
     175    { 
     176        return neg(*this); 
     177    } 
     178 
     179    /// 
     180    BigInt opCom() const 
     181    { 
     182        return com(*this); 
     183    } 
     184 
     185    /// 
     186    BigInt opPostInc() 
     187    { 
     188        BigInt n = *this; 
     189        opAddAssign(1); 
     190        return n; 
     191    } 
     192 
     193    /// 
     194    BigInt opPostDec() 
     195    { 
     196        BigInt n = *this; 
     197        opSubAssign(1); 
     198        return n; 
     199    } 
     200 
     201    // Binary operator overloads 
     202 
     203    /// 
     204    BigInt opAdd(T)(T n) const 
     205    { 
     206        return opAdd(BigInt(n)); 
     207    } 
     208 
     209    /// 
     210    BigInt opAdd(T:int)(T n) const 
     211    { 
     212        return add(*this,cast(Digit)n); 
     213    } 
     214 
     215    /// 
     216    BigInt opAdd(T:const(BigInt))(T n) const 
     217    { 
     218        return add(*this,n); 
     219    } 
     220 
     221    /// 
     222    void opAddAssign(T)(T n) 
     223    { 
     224        auto r = opAdd(n); 
     225        digits = r.digits; 
     226    } 
     227 
     228    /// 
     229    BigInt opSub(T)(T n) const 
     230    { 
     231        return opSub(BigInt(n)); 
     232    } 
     233 
     234    /// 
     235    BigInt opSub(T:int)(T n) const 
     236    { 
     237        return sub(*this,cast(Digit)n); 
     238    } 
     239 
     240    /// 
     241    BigInt opSub(T:const(BigInt))(T n) const 
     242    { 
     243        return sub(*this,n); 
     244    } 
     245 
     246    /// 
     247    void opSubAssign(T)(T n) 
     248    { 
     249        auto r = opSub(n); 
     250        digits = r.digits; 
     251    } 
     252 
     253    /// 
     254    BigInt opMul(T)(T n) const 
     255    { 
     256        return opMul(BigInt(n)); 
     257    } 
     258 
     259    /// 
     260    BigInt opMul(T:int)(T n) const 
     261    { 
     262        if (cast(int)n == int.min) return opMul(BigInt(n)); 
     263        int xs = sgn; 
     264        if (xs == 0 || n == 0) return BigInt.ZERO; 
     265        int ys = n > 0 ? 1 : -1; 
     266        auto x = abs; 
     267        auto y = n > 0 ? n : -n; 
     268        auto r = mul(x,y); 
     269        return (xs == ys) ? r : -r; 
     270    } 
     271 
     272    /// 
     273    BigInt opMul(T:const(BigInt))(T n) const 
     274    { 
     275        int xs = sgn; 
     276        int ys = n.sgn; 
     277        if (xs == 0 || ys == 0) return BigInt.ZERO; 
     278        auto x = abs; 
     279        auto y = n.abs; 
     280        auto r = mul(x,y); 
     281        return (xs == ys) ? r : -r; 
     282    } 
     283 
     284    /// 
     285    void opMulAssign(T)(T n) 
     286    { 
     287        auto r = opMul(n); 
     288        digits = r.digits; 
     289    } 
     290 
     291    /* 
     292        Here's how the signs work 
     293         7 /  3 = 2 
     294         7 %  3 = 1 
     295         7 / -3 = -2 
     296         7 % -3 = 1 
     297        -7 /  3 = -2 
     298        -7 %  3 = -1 
     299        -7 / -3 = 2 
     300        -7 % -3 = -1 
     301    */ 
     302 
     303    /// 
     304    BigInt opDiv(T)(T n) const 
     305    { 
     306        return opDiv(BigInt(n)); 
     307    } 
     308 
     309    /// 
     310    BigInt opDiv(T:int)(T n) const 
     311    { 
     312        if (n == 0) throw new Exception("Divide by zero"); 
     313        if (cast(int)n == int.min) return opDiv(BigInt(n)); 
     314        int xs = sgn; 
     315        int ys = n > 0 ? 1 : -1; 
     316        if (xs == 0) return BigInt.ZERO; 
     317        auto x = abs; 
     318        auto y = n > 0 ? n : -n; 
     319        auto r = div(x,y); 
     320        return (xs == ys) ? r.q : -r.q; 
     321    } 
     322 
     323    /// 
     324    BigInt opDiv(T:const(BigInt))(T n) const 
     325    { 
     326        int xs = sgn; 
     327        int ys = n.sgn; 
     328        if (ys == 0) throw new Exception("Divide by zero"); 
     329        if (xs == 0) return BigInt.ZERO; 
     330        auto x = abs; 
     331        auto y = n.abs; 
     332        auto r = div(x,y); 
     333        return (xs == ys) ? r.q : -r.q; 
     334    } 
     335 
     336    /// 
     337    void opDivAssign(T)(T n) 
     338    { 
     339        auto r = opDiv(n); 
     340        digits = r.digits; 
     341    } 
     342 
     343    /// 
     344    BigInt opMod(T)(T n) const 
     345    { 
     346        return opMod(BigInt(n)); 
     347    } 
     348 
     349    /// 
     350    int opMod(T:int)(T n) const 
     351    { 
     352        if (n == 0) throw new Exception("Divide by zero"); 
     353        int xs = sgn; 
     354        if (xs == 0) return n; 
     355        auto x = abs; 
     356        auto y = n > 0 ? n : -n; 
     357        auto r = div(x,y); 
     358        return (xs == 1) ? r.r : -r.r; 
     359    } 
     360 
     361    /// 
     362    BigInt opMod(T:const(BigInt))(T n) const 
     363    { 
     364        int xs = sgn; 
     365        int ys = n.sgn; 
     366        if (ys == 0) throw new Exception("Divide by zero"); 
     367        if (xs == 0) return n; 
     368        auto x = abs; 
     369        auto y = n.abs; 
     370        auto r = div(x,y); 
     371        assert(r.r.abs < n.abs); 
     372        return (xs == 1) ? r.r : -r.r; 
     373    } 
     374 
     375    /// 
     376    void opModAssign(T:int)(T n) 
     377    { 
     378        auto r = opMod(BigInt(n)); 
     379        digits = r.digits; 
     380    } 
     381 
     382    /// 
     383    void opModAssign(T)(T n) 
     384    { 
     385        auto r = opMod(n); 
     386        digits = r.digits; 
     387    } 
     388 
     389    /// 
     390    BigInt opAnd(T)(T n) const 
     391    { 
     392        return opAnd(BigInt(n)); 
     393    } 
     394 
     395    /// 
     396    BigInt opAnd(T:int)(T n) const 
     397    { 
     398        return and(*this,cast(Digit)n); 
     399    } 
     400 
     401    /// 
     402    uint opAnd(T:uint)(T n) const 
     403    { 
     404        uint t; 
     405        castTo(t); 
     406        return t & n; 
     407    } 
     408 
     409    /// 
     410    BigInt opAnd(T:const(BigInt))(T n) const 
     411    { 
     412        return and(*this,n); 
     413    } 
     414 
     415    /// 
     416    void opAndAssign(T:uint)(T n) 
     417    { 
     418        auto r = opAnd(BigInt(n)); 
     419        digits = r.digits; 
     420    } 
     421 
     422    /// 
     423    void opAndAssign(T)(T n) 
     424    { 
     425        auto r = opAnd(n); 
     426        digits = r.digits; 
     427    } 
     428 
     429    /// 
     430    BigInt opOr(T)(T n) const 
     431    { 
     432        return opOr(BigInt(n)); 
     433    } 
     434 
     435    /// 
     436    BigInt opOr(T:int)(T n) const 
     437    { 
     438        return or(*this,cast(Digit)n); 
     439    } 
     440 
     441    /// 
     442    BigInt opOr(T:const(BigInt))(T n) const 
     443    { 
     444        return or(*this,n); 
     445    } 
     446 
     447    /// 
     448    void opOrAssign(T)(T n) 
     449    { 
     450        auto r = opOr(n); 
     451        digits = r.digits; 
     452    } 
     453 
     454    /// 
     455    BigInt opXor(T)(T n) const 
     456    { 
     457        return opXor(BigInt(n)); 
     458    } 
     459 
     460    /// 
     461    BigInt opXor(T:int)(T n) const 
     462    { 
     463        return xor(*this,cast(Digit)n); 
     464    } 
     465 
     466    /// 
     467    BigInt opXor(T:const(BigInt))(T n) const 
     468    { 
     469        return xor(*this,n); 
     470    } 
     471 
     472    /// 
     473    void opXorAssign(T)(T n) 
     474    { 
     475        auto r = opXor(n); 
     476        digits = r.digits; 
     477    } 
     478 
     479    /// 
     480    BigInt opShl(uint n) const 
     481    { 
     482        uint hi = n >> 5; 
     483        uint lo = n & 0x1F; 
     484        Big r = *this; 
     485        if (lo != 0) r = shl(r,lo); 
     486        if (hi != 0) r = shlDigits(r,hi); 
     487        return r; 
     488    } 
     489 
     490    /// 
     491    void opShlAssign(uint n) 
     492    { 
     493        auto r = opShl(n); 
     494        digits = r.digits; 
     495    } 
     496 
     497    /// 
     498    BigInt opShr(uint n) const 
     499    { 
     500        uint hi = n >> 5; 
     501        uint lo = n & 0x1F; 
     502        Big r = *this; 
     503        if (lo != 0) r = shr(r,lo).q; 
     504        if (hi != 0) r = shrDigits(r,hi); 
     505        return r; 
     506    } 
     507 
     508    /// 
     509    void opShrAssign(uint n) 
     510    { 
     511        auto r = opShr(n); 
     512        digits = r.digits; 
     513    } 
     514     
     515    /// 
     516    BigInt opUShr(T)(T n) const 
     517    { 
     518        if (sgn >= 0) return opShr(n); 
     519        else throw new Exception(">>> cannot be applied to negative numbers"); 
     520    } 
     521     
     522    /// 
     523    void opUShrAssign(T)(T n) 
     524    { 
     525        if (sgn >= 0) opShrAssign(n); 
     526        else throw new Exception(">>>= cannot be applied to negative numbers"); 
     527    } 
     528 
     529    /// 
     530    int opEquals(T)(T n) const 
     531    { 
     532        return opEquals(BigInt(n)); 
     533    } 
     534 
     535    /// 
     536    int opEquals(T:int)(T n) const 
     537    { 
     538        return digits.length == 1 && digits[0] == n; 
     539    } 
     540 
     541    /// 
     542    int opEquals(T:const(BigInt))(T n) const 
     543    { 
     544        return digits == n.digits; 
     545    } 
     546 
     547    /// 
     548    int opCmp(T)(T n) const 
     549    { 
     550        return opCmp(BigInt(n)); 
     551    } 
     552 
     553    /// 
     554    int opCmp(T:int)(T n) const 
     555    { 
     556        int t = cmp(*this,n); 
     557        return t == 0 ? 0 : (t > 0 ? 1 : -1); 
     558    } 
     559 
     560    /// 
     561    int opCmp(T:const(BigInt))(T n) const 
     562    { 
     563        int t = cmp(*this,n); 
     564        return t == 0 ? 0 : (t > 0 ? 1 : -1); 
     565    } 
     566 
     567    /// 
     568    string toString() const 
     569    { 
     570        return decimal(*this); 
     571    } 
     572 
     573    /// 
     574    hash_t toHash() const 
     575    { 
     576        hash_t h = 0; 
     577        foreach(Digit d;digits) { h += d; } 
     578        return h; 
     579    } 
     580 
     581    private int sgn() const 
     582    { 
     583        int t = cmp(*this,0); 
     584        return t == 0 ? 0 : (t > 0 ? 1 : -1); 
     585    } 
     586 
     587    private BigInt abs() const 
     588    { 
     589        return sgn >= 0 ? opPos : opNeg; 
     590    } 
     591} 
     592 
     593// ----------- EVERYTHING PRIVATE BEYOND THIS POINT ----------- 
     594private: 
     595 
     596// Aliases and Typedefs 
     597 
     598alias BigInt                Big; 
     599alias invariant(Digit)[]    Digits; 
     600typedef Digit[]             DownArray; 
     601typedef Digit*              DownPtr; 
     602alias int                   SignedDigit; 
     603alias long                  SignedWideDigit; 
     604alias Digit                 Unused; 
     605typedef Digit[]             UpArray; 
     606typedef Digit*              UpPtr; 
     607alias ulong                 WideDigit; 
     608 
     609struct Big_Digit { Big q; Digit r; } 
     610struct Big_Big{ Big q; Big r; } 
     611 
     612// Endianness 
     613 
     614// The constant BIG_ENDIAN determines the ordering of digits within arrays. 
     615// If BIG_ENDIAN is true, then bigints are stored most significant digit first. 
     616// If BIG_ENDIAN is true, then bigints are stored most significant digit last. 
     617 
     618// Note that this does not necessarily have to be the same as the endianness 
     619// of the platform architecture. 
     620 
     621// Setting BIG_ENDIAN opposite to platform endianness allows unittests 
     622// to run in reverse endianness. (And they still pass). 
     623 
     624version(BigEndian) { enum bool BIG_ENDIAN = true;  } 
     625else               { enum bool BIG_ENDIAN = false; } 
     626 
     627// String conversion 
     628 
     629void parseError() 
     630{ 
     631    throw new Exception("Parse Error"); 
     632} 
     633 
     634Big fromString(string s) 
     635{ 
     636    if (s.length == 0) parseError(); 
     637    if (s[0] == '-') 
     638    { 
     639        return -fromString(s[1..$]); 
     640    } 
     641    if (s.length > 2 && s[0] == '0' && (s[1] == 'x' || s[1] == 'X')) 
     642    { 
     643        return fromHex(s[2..$]); 
     644    } 
     645    return fromDecimal(s); 
     646} 
     647 
     648Big fromDecimal(string s) 
     649{ 
     650    bool invalid = true; 
     651    Big r = Big.ZERO; 
     652    foreach(char c;s) 
     653    { 
     654        if (c == '_') continue; 
     655        if (c < '0' || c > '9') parseError(); 
     656        invalid = false; 
     657        //r = 10 * r + (c - '0'); 
     658        r *= 10; 
     659        r += (c - '0'); 
     660    } 
     661    if (invalid) parseError(); 
     662    return r; 
     663} 
     664 
     665Big fromHex(string s) 
     666{ 
     667    bool invalid = true; 
     668    Big r = Big.ZERO; 
     669    foreach(char c;s) 
     670    { 
     671        switch(c) 
     672        { 
     673        case '_': 
     674            continue; 
     675 
     676        case '0','1','2','3','4','5','6','7','8','9': 
     677            r = (r << 4) + (c - '0'); 
     678            invalid = false; 
     679            break; 
     680 
     681        case 'A','B','C','D','E','F': 
     682            r = (r << 4) + (c - 'A' + 10); 
     683            invalid = false; 
     684            break; 
     685 
     686        case 'a','b','c','d','e','f': 
     687            r = (r << 4) + (c - 'A' + 10); 
     688            invalid = false; 
     689            break; 
     690 
     691        default: 
     692            parseError(); 
     693        } 
     694    } 
     695    if (invalid) parseError(); 
     696    return r; 
     697} 
     698 
     699string decimal(Big b) 
     700{ 
     701    if (b < 0) return "-" ~ decimal(-b); 
     702 
     703    if (b < 10) 
     704    { 
     705        int n; 
     706        b.castTo(n); 
     707        char c = cast(char)(n + '0'); 
     708        return [ c ]; 
     709    } 
     710 
     711    auto t = div(b,10); 
     712    auto r = cast(string)(decimal(t.q) ~ [ cast(char)(t.r + '0') ]); 
     713    return r; 
     714} 
     715 
     716// Shrinking 
     717 
     718Big bigInt(DownArray a) 
     719{ 
     720    if (a.length == 0) return Big.ZERO; 
     721 
     722    Big r; 
     723 
     724    if (a.length == 1) 
     725    { 
     726        r.digits = cast(Digits)a; 
     727    } 
     728    else 
     729    { 
     730        auto xp = begin(a); 
     731        auto xe = end(a); 
     732        auto d1 = peek(xp); 
     733        xp = next(xp); 
     734        auto s = signOf(d1); 
     735        while (xp != xe) 
     736        { 
     737            if (d1 != s) break; 
     738            auto d2 = peek(xp); 
     739            if (signOf(d2) != s) break; 
     740            xp = next(xp); 
     741            d1 = d2; 
     742        } 
     743        r.digits = freezeRange(xp, xe); 
     744    } 
     745    return r; 
     746} 
     747 
     748static if(BIG_ENDIAN) 
     749{ 
     750    alias UpArray       BwdArray; 
     751    alias UpPtr         BwdPtr; 
     752    alias DownArray     FwdArray; 
     753    alias DownPtr       FwdPtr; 
     754 
     755    Digit[] join(Digit[] t, Digit[] u) { return t ~ u; } 
     756    T head(T)(T t,size_t n) { return cast(T)(t[0..n]); } 
     757    T tail(T)(T t,size_t n) { return cast(T)(t[$-n..$]); } 
     758} 
     759else 
     760{ 
     761    alias DownArray     BwdArray; 
     762    alias DownPtr       BwdPtr; 
     763    alias UpArray       FwdArray; 
     764    alias UpPtr         FwdPtr; 
     765 
     766    Digit[] join(Digit[] t, Digit[] u) { return u ~ t; } 
     767    T head(T)(T t,size_t n) { return cast(T)(t[$-n..$]); } 
     768    T tail(T)(T t,size_t n) { return cast(T)(t[0..n]); } 
     769} 
     770 
     771// Really simple functions 
     772 
     773FwdPtr advance(FwdPtr p,size_t n) { return p + n; } 
     774BwdPtr advance(BwdPtr p,size_t n) { return p - n; } 
     775 
     776Digit begin(Digit a) { return a; } 
     777FwdPtr begin(FwdArray a) { return cast(FwdPtr)(a.ptr); } 
     778BwdPtr begin(BwdArray a) { return cast(BwdPtr)(a.ptr + a.length - 1); } 
     779 
     780Big bigInt(Digit a) { Digit[] t; t.length = 1; t[0] = a; return bigInt(cast(DownArray)t); } 
     781Big bigInt(Digits a) { return bigInt(cast(DownArray)a); } 
     782Big bigInt(Digit[] a) { return bigInt(cast(DownArray)a); } 
     783Big bigInt(UpArray a) { return bigInt(cast(DownArray)a); } 
     784 
     785Digit downArray(Digit a) { return a; } 
     786DownArray downArray(Digit[] a) { return cast(DownArray)a; } 
     787DownArray downArray(Big a) { return cast(DownArray)(a.digits); } 
     788 
     789Digit end(Digit a) { return a; } 
     790FwdPtr end(FwdArray a) { return cast(FwdPtr)(a.ptr + a.length); } 
     791BwdPtr end(BwdArray a) { return cast(BwdPtr)(a.ptr + - 1); } 
     792 
     793Digit first(Digit a) { return a; } 
     794Digit first(FwdArray a) { return a[0]; } 
     795Digit first(BwdArray a) { return a[$-1]; } 
     796 
     797Digits freezeRange(FwdPtr p, FwdPtr q) { return cast(Digits)((p-1)[0..(q-p+1)]); } 
     798Digits freezeRange(BwdPtr p, BwdPtr q) { return cast(Digits)((q+1)[0..(p-q+1)]); } 
     799 
     800Digit last(Digit a) { return a; } 
     801Digit last(FwdArray a) { return a[$-1]; } 
     802Digit last(BwdArray a) { return a[0]; } 
     803 
     804Digit lsd(Digit a) { return a; } 
     805Digit lsd(DownArray a) { return last(a); } 
     806Digit lsd(UpArray a) { return first(a); } 
     807 
     808Digit msd(Digit a) { return a; } 
     809Digit msd(DownArray a) { return first(a); } 
     810Digit msd(UpArray a) { return last(a); } 
     811 
     812size_t lengthOf(Digit a) { return 1; } 
     813size_t lengthOf(Big a) { return a.digits.length; } 
     814size_t lengthOf(DownArray a) { return a.length; } 
     815size_t lengthOf(UpArray a) { return a.length; } 
     816 
     817Digit next(Digit d) { return d; } 
     818FwdPtr next(FwdPtr p) { return p + 1; } 
     819BwdPtr next(BwdPtr p) { return p - 1; } 
     820 
     821Digit peek(Digit d) { return d; } 
     822Digit peek(Digit* p) { return *p; } 
     823 
     824void poke(DownPtr p,Digit d) { *p = d; } 
     825void poke(DownPtr p,WideDigit d) { *p = cast(Digit)d; } 
     826void poke(DownPtr p,SignedWideDigit d) { *p = cast(Digit)d; } 
     827void poke(UpPtr p,Digit d) { *p = d; } 
     828void poke(UpPtr p,WideDigit d) { *p = cast(Digit)d; } 
     829void poke(UpPtr p,SignedWideDigit d) { *p = cast(Digit)d; } 
     830 
     831Big shrink(Big a) { return bigInt(cast(DownArray)(a.digits)); } 
     832 
     833FwdArray slice(FwdPtr ptr, size_t len) { return cast(FwdArray)(ptr[0..len]); } 
     834BwdArray slice(BwdPtr ptr, size_t len) { return cast(BwdArray)((ptr-len+1)[0..len]); } 
     835 
     836Digit signOf(SignedDigit d) { return d < 0 ? -1 : 0; } 
     837 
     838Digit upArray(Digit a) { return a; } 
     839UpArray upArray(Digit[] a) { return cast(UpArray)a; } 
     840UpArray upArray(in Big a) { return cast(UpArray)(a.digits); } 
     841 
     842// Core functions 
     843 
     844WideDigit addCore(Digit x,Digit y,WideDigit c) { return (c + x) + y; } 
     845 
     846Digit andCore(Digit x,Digit y,Digit c) { return x & y; } 
     847 
     848WideDigit divCore(Digit x,Digit y,WideDigit c) 
     849{ 
     850    c <<= 32; 
     851    c += x; 
     852    WideDigit r = c % y; 
     853    c /= y; 
     854    c += r << 32; 
     855    return c; 
     856} 
     857 
     858WideDigit shlCore(Digit x,Digit y,WideDigit c) { return c + (cast(WideDigit)x << y); } 
     859 
     860WideDigit shrCore(Digit x,Digit y,WideDigit c) { return c + (x >> y) + (cast(WideDigit)x << (64-y)); } 
     861 
     862WideDigit mulCore(Digit x,Digit y,WideDigit c) { return c + (cast(WideDigit)x * y); } 
     863 
     864WideDigit subCore(Digit x,Digit y,WideDigit c) { return (c + x) - y; } 
     865 
     866Digit orCore(Digit x,Digit y,Digit c) { return x | y; } 
     867 
     868Digit xorCore(Digit x,Digit y,Digit c) { return x ^ y; } 
     869 
     870// Update functions 
     871 
     872Digit updateDigit(Digit c) { return c; } 
     873 
     874WideDigit updateShr(WideDigit c) { return cast(SignedWideDigit)c >> 32; } 
     875 
     876WideDigit updateUShr(WideDigit c) { return c >> 32; } 
     877 
     878// Helper functions 
     879 
     880int cmp(DownPtr xp, DownPtr xe, DownPtr yp) 
     881{ 
     882    while (xp != xe) 
     883    { 
     884        auto xd = peek(xp); 
     885        auto yd = peek(yp); 
     886        if (xd < yd) return -1; 
     887        if (xd > yd) return 1; 
     888        xp = next(xp); 
     889        yp = next(yp); 
     890    } 
     891    return 0; 
     892} 
     893 
     894void mulInner(Big a, UpPtr rp, WideDigit y) 
     895{ 
     896    WideDigit c; 
     897    mixin(setUp("x","a")); 
     898 
     899    while (xp != xe) 
     900    { 
     901        c += y * peek(xp) + peek(rp); 
     902        poke(rp,c); 
     903        xp = next(xp); 
     904        rp = next(rp); 
     905        c = updateShr(c); 
     906    } 
     907 
     908    mixin(runOnce(   "mulCore","updateShr","xs","y")); 
     909} 
     910 
     911void divInner(DownPtr xp, DownPtr cachePtr, size_t len) 
     912{ 
     913    Digit result; 
     914 
     915    debug // sanity checking 
     916    { 
     917        DownArray _divisor = slice(cachePtr,32*len); 
     918        _divisor = tail(_divisor,len); 
     919    } 
     920 
     921    DownPtr rp = xp; 
     922    xp = next(xp); 
     923    DownPtr xe = advance(xp,len); 
     924 
     925    debug // sanity checking 
     926    { 
     927        DownArray _remainder = slice(xp,len);  // will be modified in-place 
     928        DownArray _original = cast(DownArray)_remainder.dup;  // but we'll keep this one 
     929    } 
     930 
     931    for (Digit mask=0x80000000; mask!=0; mask>>=1) 
     932    { 
     933        int t = cmp(xp,xe,cachePtr); 
     934        if (t >= 0) 
     935        { 
     936            debug // sanity checking 
     937            { 
     938                DownArray _after = slice(xp,len);   // will be modified in-place 
     939                DownArray _before = cast(DownArray)_after.dup;     // but we'll keep this one 
     940                DownArray _test = slice(cachePtr,len); 
     941            } 
     942 
     943            result += mask; 
     944            Digit carry = subInPlace(xp,cachePtr,len); 
     945            debug 
     946            { 
     947                BigInt before = bigInt(_before); 
     948                BigInt test = bigInt(_test); 
     949                BigInt after = bigInt(_after); 
     950 
     951                assert(after + test == before); 
     952                assert(carry == 0); 
     953            } 
     954        } 
     955        cachePtr = advance(cachePtr,len); 
     956    } 
     957 
     958    debug // sanity checking 
     959    { 
     960        // in theory, quotient * _divisor + _remainder == _original 
     961        BigInt quotient = result; 
     962        BigInt divisor = bigInt(_divisor); 
     963        BigInt remainder = bigInt(_remainder); 
     964        BigInt original = bigInt(_original); 
     965 
     966        assert(quotient * divisor + remainder == original); 
     967    } 
     968 
     969    poke(rp,result); 
     970} 
     971 
     972Digit subInPlace(DownPtr downPtrX, DownPtr downPtrY, size_t len) 
     973{ 
     974    UpPtr xp = cast(UpPtr)(advance(downPtrX,len-1)); 
     975    UpPtr yp = cast(UpPtr)(advance(downPtrY,len-1)); 
     976    UpPtr xe = advance(xp,len); 
     977 
     978    SignedWideDigit c; 
     979 
     980    while (xp != xe) 
     981    { 
     982        c += peek(xp); 
     983        c -= peek(yp); 
     984        poke(xp,c); 
     985        c >>= 32; 
     986        xp = next(xp); 
     987        yp = next(yp); 
     988    } 
     989 
     990    return cast(Digit)c; 
     991} 
     992 
     993DownPtr makeDivCache(DownArray y) 
     994{ 
     995    // Pad with a leading zero 
     996    auto paddedY = cast(UpArray)(join([Digit.init],y)); 
     997 
     998    auto upCache = cast(UpArray)new Digit[32 * paddedY.length]; 
     999    auto rp = begin(upCache); 
     1000 
     1001    // Fill upCache by successively leftshifting x by one bit 
     1002    for (int i=0; i<32; ++i) 
     1003    { 
     1004        auto xp = begin(paddedY); 
     1005        auto xe = end(paddedY); 
     1006 
     1007        WideDigit c; 
     1008 
     1009        // Shift lefy by one bit 
     1010        while (xp != xe) 
     1011        { 
     1012            Digit xd = peek(xp); 
     1013            poke(rp,xd); 
     1014            c += xd; 
     1015            c += xd; 
     1016            poke(xp,c); 
     1017            c >>= 32; 
     1018            xp = next(xp); 
     1019            rp = next(rp); 
     1020        } 
     1021    } 
     1022 
     1023    auto downCache = cast(DownArray)upCache; 
     1024 
     1025    static if(false) // make true to display cache 
     1026    { 
     1027        for (int j=0; j<32; ++j) 
     1028        { 
     1029            writefln("bit %02d: ",31-j,hex(downCache[paddedY.length*j..paddedY.length*(j+1)])); 
     1030        } 
     1031    } 
     1032 
     1033    return begin(downCache); 
     1034} 
     1035 
     1036// Mixins 
     1037 
     1038string runOnce(string core, string updater, string xp, string yp) 
     1039{ 
     1040    return 
     1041    "{ 
     1042        auto xd = peek("~xp~"); 
     1043        auto yd = peek("~yp~"); 
     1044        c = "~core~"(xd,yd,c); 
     1045        poke(rp,c); 
     1046        "~xp~" = next("~xp~"); 
     1047        "~yp~" = next("~yp~"); 
     1048        rp = next(rp); 
     1049        c = "~updater~"(c); 
     1050    }"; 
     1051} 
     1052 
     1053string runTo(string dest, string core, string updater, string xp, string yp) 
     1054{ 
     1055    string s = runOnce(core,updater,xp,yp); 
     1056    return 
     1057    " 
     1058        static if(isIntegral!(typeof("~dest~"))) {"~s~"} 
     1059        else 
     1060        { 
     1061            while("~dest[0..1]~"p!="~dest~") {"~s~"} 
     1062        } 
     1063    "; 
     1064} 
     1065 
     1066string setDown(string x,string a) 
     1067{ 
     1068    return 
     1069    " 
     1070        auto "~x~" = downArray("~a~"); 
     1071        auto "~x~"p = begin("~x~"); 
     1072        auto "~x~"e = end("~x~"); 
     1073        auto "~x~"s = signOf(msd("~x~")); 
     1074    "; 
     1075} 
     1076 
     1077string setUp(string x,string a) 
     1078{ 
     1079    return 
     1080    " 
     1081        auto "~x~" = upArray("~a~"); 
     1082        auto "~x~"p = begin("~x~"); 
     1083        auto "~x~"e = end("~x~"); 
     1084        auto "~x~"s = signOf(msd("~x~")); 
     1085    "; 
     1086} 
     1087 
     1088// BigInt functions 
     1089 
     1090Big neg(Big b) 
     1091{ 
     1092    auto r = upArray(new Digit[lengthOf(b) + 1]); 
     1093    auto rp = begin(r); 
     1094 
     1095    SignedDigit a; 
     1096    WideDigit c; 
     1097 
     1098    mixin(setUp("x","a")); 
     1099    mixin(setUp("y","b")); 
     1100 
     1101    mixin(runOnce(   "subCore","updateShr","xp","yp")); 
     1102    mixin(runTo("ye","subCore","updateShr","xs","yp")); 
     1103    mixin(runOnce(   "subCore","updateShr","xs","ys")); 
     1104 
     1105    return bigInt(r); 
     1106} 
     1107 
     1108Big com(Big a) 
     1109{ 
     1110    auto r = upArray(new Digit[lengthOf(a)]); 
     1111    auto rp = begin(r); 
     1112 
     1113    Digit c; 
     1114 
     1115    mixin(setUp("x","a")); 
     1116    Digit ys = uint.max; 
     1117 
     1118    mixin(runTo("xe","xorCore","updateDigit","xp","ys")); 
     1119 
     1120    return bigInt(r); 
     1121} 
     1122 
     1123Big add(Other)(Big a,Other b) 
     1124{ 
     1125    static if(is(Other==BigInt)) if (lengthOf(a) < lengthOf(b)) 
     1126    { 
     1127    swap(a,b); 
     1128    } 
     1129    auto r = upArray(new Digit[max(lengthOf(a),lengthOf(b)) + 1]); 
     1130    auto rp = begin(r); 
     1131 
     1132    WideDigit c; 
     1133 
     1134    mixin(setUp("x","a")); 
     1135    mixin(setUp("y","b")); 
     1136 
     1137    mixin(runTo("ye","addCore","updateShr","xp","yp")); 
     1138    mixin(runTo("xe","addCore","updateShr","xp","ys")); 
     1139    mixin(runOnce(   "addCore","updateShr","xs","ys")); 
     1140 
     1141    return bigInt(r); 
     1142} 
     1143 
     1144Big sub(Big a,Big b) 
     1145{ 
     1146    auto r = upArray(new Digit[max(lengthOf(a),lengthOf(b)) + 1]); 
     1147    auto rp = begin(r); 
     1148    auto re = advance(rp,min(lengthOf(a),lengthOf(b))); 
     1149 
     1150    WideDigit c; 
     1151 
     1152    mixin(setUp("x","a")); 
     1153    mixin(setUp("y","b")); 
     1154 
     1155    mixin(runTo("re","subCore","updateShr","xp","yp")); 
     1156    if (lengthOf(x) >= lengthOf(y)) 
     1157    { 
     1158        mixin(runTo("xe","subCore","updateShr","xp","ys")); 
     1159    } 
     1160    else 
     1161    { 
     1162        mixin(runTo("ye","subCore","updateShr","xs","yp")); 
     1163    } 
     1164    mixin(runOnce("subCore","updateShr","xs","ys")); 
     1165 
     1166    return bigInt(r); 
     1167} 
     1168 
     1169Big sub(Big a, Digit b) 
     1170{ 
     1171    auto r = upArray(new Digit[lengthOf(a) + 1]); 
     1172    auto rp = begin(r); 
     1173 
     1174    WideDigit c; 
     1175 
     1176    mixin(setUp("x","a")); 
     1177    mixin(setUp("y","b")); 
     1178 
     1179    mixin(runOnce(   "subCore","updateShr","xp","yp")); 
     1180    mixin(runTo("xe","subCore","updateShr","xp","ys")); 
     1181    mixin(runOnce(   "subCore","updateShr","xs","ys")); 
     1182 
     1183    return bigInt(r); 
     1184} 
     1185 
     1186Big mul(Big a, Big b) // a and b must be positive 
     1187{ 
     1188    auto r = upArray(new Digit[lengthOf(a) + lengthOf(b)]); 
     1189    auto rp = begin(r); 
     1190 
     1191    mixin(setUp("y","b")); 
     1192    while (yp != ye) 
     1193    { 
     1194        mulInner(a,rp,peek(yp)); 
     1195        yp = next(yp); 
     1196        rp = next(rp); 
     1197    } 
     1198 
     1199    return bigInt(r); 
     1200} 
     1201 
     1202Big mul(Big a, Digit b) // a and b must be positive 
     1203{ 
     1204    auto r = upArray(new Digit[lengthOf(a) + 1]); 
     1205    auto rp = begin(r); 
     1206 
     1207    WideDigit c; 
     1208 
     1209    mixin(setUp("x","a")); 
     1210 
     1211    mixin(runTo("xe","mulCore","updateShr","xp","b")); 
     1212    poke(rp,c); 
     1213 
     1214    return bigInt(r); 
     1215} 
     1216 
     1217Big_Big div(Big a, Big b) // a and b must be positive 
     1218{ 
     1219    auto lenX = lengthOf(a); 
     1220    auto lenY = lengthOf(b) + 1; 
     1221 
     1222    auto r = cast(DownArray)join(new Digit[lenY], cast(Digit[])a.digits); 
     1223    auto rp = begin(r); 
     1224    auto re = advance(rp,lenX); 
     1225 
     1226    auto y = downArray(b); 
     1227    auto cache = makeDivCache(y); 
     1228 
     1229    while (rp != re) 
     1230    { 
     1231        divInner(rp, cache, lenY); 
     1232        rp = next(rp); 
     1233    } 
     1234 
     1235    Big quotient  = bigInt(cast(DownArray)(head(r,lenX))); 
     1236    Big remainder = bigInt(cast(DownArray)(tail(r,lenY))); 
     1237 
     1238    return Big_Big(quotient,remainder); 
     1239} 
     1240 
     1241Big_Digit div(Big a, Digit b) // a and b must be positive 
     1242{ 
     1243    auto r = downArray(new Digit[lengthOf(a)]); 
     1244    auto rp = begin(r); 
     1245 
     1246    WideDigit c; 
     1247 
     1248    mixin(setDown("x","a")); 
     1249 
     1250    mixin(runTo("xe","divCore","updateUShr","xp","b")); 
     1251 
     1252    return Big_Digit(bigInt(r),cast(Digit)c); 
     1253} 
     1254 
     1255Big and(Other)(Big a, Other b) 
     1256{ 
     1257    static if(is(Other==BigInt)) if (lengthOf(a) < lengthOf(b)) { swap(a,b); } 
     1258    auto r = upArray(new Digit[max(lengthOf(a),lengthOf(b))]); 
     1259    auto rp = begin(r); 
     1260 
     1261    Digit c; 
     1262 
     1263    mixin(setUp("x","a")); 
     1264    mixin(setUp("y","b")); 
     1265 
     1266    mixin(runTo("ye","andCore","updateDigit","xp","yp")); 
     1267    mixin(runTo("xe","andCore","updateDigit","xp","ys")); 
     1268 
     1269    return bigInt(r); 
     1270} 
     1271 
     1272Big or(Other)(Big a, Other b) 
     1273{ 
     1274    static if(is(Other==BigInt)) if (lengthOf(a) < lengthOf(b)) { swap(a,b); } 
     1275    auto r = upArray(new Digit[max(lengthOf(a),lengthOf(b))]); 
     1276    auto rp = begin(r); 
     1277 
     1278    Digit c; 
     1279 
     1280    mixin(setUp("x","a")); 
     1281    mixin(setUp("y","b")); 
     1282 
     1283    mixin(runTo("ye","orCore","updateDigit","xp","yp")); 
     1284    mixin(runTo("xe","orCore","updateDigit","xp","ys")); 
     1285 
     1286    return bigInt(r); 
     1287} 
     1288 
     1289Big xor(Other)(Big a, Other b) 
     1290{ 
     1291    static if(is(Other==BigInt)) if (lengthOf(a) < lengthOf(b)) { swap(a,b); } 
     1292    auto r = upArray(new Digit[max(lengthOf(a),lengthOf(b))]); 
     1293    auto rp = begin(r); 
     1294 
     1295    Digit c; 
     1296 
     1297    mixin(setUp("x","a")); 
     1298    mixin(setUp("y","b")); 
     1299 
     1300    mixin(runTo("ye","xorCore","updateDigit","xp","yp")); 
     1301    mixin(runTo("xe","xorCore","updateDigit","xp","ys")); 
     1302 
     1303    return bigInt(r); 
     1304} 
     1305 
     1306Big shl(Big a, Digit b) 
     1307{ 
     1308    auto r = upArray(new Digit[lengthOf(a) + 1]); 
     1309    auto rp = begin(r); 
     1310 
     1311    WideDigit c; 
     1312 
     1313    mixin(setUp("x","a")); 
     1314 
     1315    mixin(runTo("xe","shlCore","updateShr","xp","b")); 
     1316    poke(rp,c); 
     1317 
     1318    return bigInt(r); 
     1319} 
     1320 
     1321Big shlDigits(Big a, uint n) 
     1322{ 
     1323    Big b; 
     1324    b.digits = cast(Digits)join(a.digits.dup, new Digit[n]); 
     1325    return b; 
     1326} 
     1327 
     1328Big_Digit shr(Big a, Digit b) 
     1329{ 
     1330    auto r = downArray(new Digit[lengthOf(a)]); 
     1331    a