Project Euler 80

Euler80 1.d

import std.stdio;
import std.datetime.stopwatch;
import std.typecons;
import std.algorithm;
import std.range;
import std.math;
import std.traits;

void main() {
  StopWatch timer;
  byte[] digs;
  uint digsum = 0;
  double root = 0;
  // precomputed indices determining how far to iterate continued fractions to get 52 digits in the denominator, resulting in 100 decimal digits, the required precision.
  int[] ndxs = [134, 180, 82, 103, 172, 134, 65, 79, 90, 165, 140, 114, 57, 67, 122, 82, 150, 120, 122, 104, 51, 60, 86, 120, 77, 119, 135, 124, 112, 96, 48, 55, 61, 65, 86, 73, 135, 160, 123, 132, 104, 90, 45, 52, 100, 100, 103, 92, 70, 126, 157, 102, 116, 119, 98, 86, 43, 49, 103, 57, 98, 114, 108, 68, 109, 135, 120, 122, 108, 102, 94, 82, 41, 47, 50, 90, 120, 59, 120, 86, 66, 119, 123, 118, 125, 108, 104, 141, 90, 80];
  
  timer.start();
  writeln("Square root convergents\n");
  
  foreach (n; iota(1, 101).setDifference(sequence!((a, n) => (n + 1) ^^ 2))) {
    auto r = n.squareRootSequence
      .continuedFraction!BigInt
      .drop(ndxs.front);
    digs = divdigs(r.front.expand, 100);
    digsum += digs.sum();
    ndxs.popFront();
  }

  writeln(digsum);
  timer.stop();

  writefln("Finished in %s milliseconds.", timer.peek.total!"msecs"());
}

byte[] divdigs(BigInt num, BigInt denom, int digitCount) {
  BigInt quo;
  byte[] digs;

  do {
    quo = num / denom;
    num %= denom;
    digs ~= quo.digitBytes();
    if (num == 0)
      break;
    num *= 10;
  } while (digs.length < digitCount);
  
  return digs;
}

auto squareRootSequence(int number) {
  int a0 = cast(int)sqrt(cast(real)number);
  int[] result;
  int a;
  int b = a0;
  int d = 1;
  int n = 1;
  int t;

  result ~= a0;

  do {
    d = number - b^^2;
    if (d == 0)
      break;
    d /= n;
    t = a0 + b;
    a = t / d;
    result ~= a;
    n = d;
    b = a0 - (t % d);
  } while (d != 1);

  if (result.length == 1) {
    result ~= 0;
  }

  return only(result[0]).chain(result[1 .. $].cycle());
}

auto continuedFraction(T, R)(R terms) if (isInputRange!(Unqual!R) && isInfinite!(Unqual!R) && isIntegral!(ElementType!R) && (isIntegral!T || is(T == BigInt))) {
  return ContinuedFraction!(R, T)(terms);
}

auto continuedFraction(R)(R terms) if (isInputRange!(Unqual!R) && isInfinite!(Unqual!R) && isIntegral!(ElementType!R)) {
  return ContinuedFraction!(R)(terms);
}

struct ContinuedFraction(R, T = ElementType!R) if (isInputRange!(Unqual!R) && isInfinite!(Unqual!R) && isIntegral!(ElementType!R) && (isIntegral!T || is(T == BigInt))) {
  import core.exception;
  alias E = ElementType!R;
  size_t j = 0;
  enum bool empty = false;
  R terms;
  E[] cache;

  this(R _terms) {
    terms = _terms;
    this.popFront();
  }

  auto front() @property {
    Tuple!(T, T) result = tuple(T(0), T(1));
    
    if (cache.back == 0) {
      throw new RangeError("ContinuedFraction ran into a zero");
    }

    foreach (item; cache[0 .. $].retro()) {
      result[0] = T(item) * result[1] + result[0];
      swap(result[0], result[1]);
    }

    swap(result[0], result[1]);
    return result;
  }

  auto popFront() {
    cache ~= terms.front;
    terms.popFront();
  }
}

enum MulThreshold = 70;

mixin template RvalueRef() {
  alias T = typeof(this);
  static assert (is(T == struct));

  @nogc @safe
  ref const(T) byRef() const pure nothrow return
  {
    return this;
  }
}

struct BigInt {
private:
  byte[] mant = [0];
  bool sign = false;

  static void addAbsInPlace(ref byte[] lhs, const(byte)[] rhs) {
    byte carry = 0;

    if (lhs.length < rhs.length)
      lhs.length = rhs.length;

    for (ulong i = 0; i < lhs.length; i++) {
      lhs[i] += (rhs.length > i ? rhs[i] : 0) + carry;
      carry = lhs[i] / 10;
      if (lhs[i] > 9)
        lhs[i] -= 10;
    }

    if (carry)
      lhs ~= carry;
  }

  static byte[] addAbs(const(byte)[] lhs, const(byte)[] rhs) {
    bool leftBigger = lhs.length > rhs.length;
    const(byte)[] big = leftBigger ? lhs : rhs;
    const(byte)[] little = leftBigger ? rhs : lhs;
    byte[] temp = big.dup;

    addAbsInPlace(temp, little);
    return temp;
  }
  
  static void subAbsInPlace(ref byte[] lhs, const(byte)[] rhs) {
    byte borrow = 0;

    for (ulong i = 0; i < lhs.length; i++) {
      lhs[i] -= (i < rhs.length ? rhs[i] : 0) + borrow;
      if (lhs[i] < 0) {
        lhs[i] += 10;
        borrow = 1;
      } else {
        borrow = 0;
        if (i >= rhs.length)
          break;
      }
    }

    while (lhs[$-1] == 0 && lhs.length > 1)
      lhs.length--;
  }

  static byte[] subAbs(const(byte)[] lhs, const(byte)[] rhs) {
    byte[] difference = lhs.dup;
    subAbsInPlace(difference, rhs);
    return difference;
  }

  static void incAbs(ref byte[] lhs) {
    byte carry = 1;

    for (ulong i; i < lhs.length; i++) {
      lhs[i] += carry;
      if (lhs[i] == 10) {
        lhs[i] = 0;
      } else {
        carry = 0;
        break;
      }
    }

    if (carry) {
      lhs ~= 1;
    }
  }

  static void decAbs(ref byte[] lhs) {
    byte borrow = 0;

    ulong i = 0;

    do {
      if (lhs[i] > 0) {
        lhs[i]--;
        borrow = 0;
      } else {
        lhs[i] = 9;
        borrow = 1;
      }
      i++;
    } while (borrow);

    if (lhs[$ - 1] == 0 && lhs.length > 1)
      lhs.length--;

  }
  
  static int cmpAbs(const(byte)[] left, const(byte)[] right) {
    if (left.length < right.length)
      return -1;
    else if (left.length > right.length)
      return 1;

    for (ulong i = left.length - 1; i < ulong.max; i--)
      if (left[i] < right[i])
        return -1;
      else if (left[i] > right[i])
        return 1;

    return 0;
  }
  

  static byte[] karatsuba(const(byte)[] lhs, const(byte)[] rhs) {
    byte[] z0;
    byte[] z1;
    byte[] z2;
    ulong m;

    // Take care of base case where one or both operands are of length 1
    //if (lhs.length == 1)
      //return mulSingleDigit(rhs, lhs[0]);
    //else if (rhs.length == 1)
      //return mulSingleDigit(lhs, rhs[0]);
    if (lhs.length < MulThreshold && rhs.length < MulThreshold)
      return mulSmall(lhs, rhs);

    m = lhs.length > rhs.length ? lhs.length / 2 : rhs.length / 2;

    // Split and handle out-of-bounds indices
    auto highLeft = m >= lhs.length ? cast(byte[])[0] : lhs[m .. $];
    auto lowLeft = m >= lhs.length ? lhs : lhs[0 .. m];
    auto highRight = m >= rhs.length ? cast(byte[])[0] : rhs[m .. $];
    auto lowRight = m >= rhs.length ? rhs : rhs[0 .. m];

    // Handle leading zeros
    while (lowLeft[$ - 1] == 0 && lowLeft.length > 1)
      lowLeft.length--;
    while(lowRight[$ - 1] == 0 && lowRight.length > 1)
      lowRight.length--;

    z2 = karatsuba(highLeft, highRight);
    z0 = karatsuba(lowLeft, lowRight);
    z1 = karatsuba(addAbs(lowLeft, highLeft), addAbs(lowRight, highRight));

    // Imperative style. Really efficient.
    subAbsInPlace(z1, z2);
    subAbsInPlace(z1, z0);
    shiftBigInPlace(z1, m);
    shiftBigInPlace(z2, 2 * m);
    addAbsInPlace(z2, z1);
    addAbsInPlace(z2, z0);
    return z2;
  }

  static byte[] mulSmall(const(byte)[] left, const(byte)[] right) {
    byte carry;
    byte d;
    byte[] result;
    byte[] part;

    result.reserve(left.length + right.length + 1);

    if (left.length < right.length)
      swap(left, right);

    for (size_t i = 0; i < right.length; i++) {
      part = mulSingleDigit(left, right[i]);
      shiftBigInPlace(part, i);
      addAbsInPlace(result, part);
      carry = 0;
      part = [];
    }

    return result;
  }

  static void shiftBigInPlace(ref byte[] lhs, ulong n) {
    if (lhs[$ - 1] == 0)
      return;

    lhs.length += n;

    for (ulong i = lhs.length; i > n; i--)
      lhs[i - 1] = lhs[i - n - 1];

    lhs[0 .. n] = 0;
  }

  // Multiplies by a power of 10
  static byte[] shiftBig(byte[] lhs, ulong n) {
    byte[] temp = lhs.dup;
    shiftBigInPlace(temp, n);
    return temp;
  }

  static void shiftLittleInPlace(ref byte[] lhs, ulong n) {
    if (lhs[$ - 1] == 0)
      return;

    if (n >= lhs.length) {
      lhs.length = 1;
      lhs[0] = 0;
      return;
    }

    for (ulong i = 0; i < lhs.length - n; i++)
      lhs[i] = lhs[i + n];

    lhs.length = lhs.length - n;
  }

  static byte[] shiftLittle(byte[] lhs, ulong n) {
    byte[] temp = lhs.dup;
    shiftLittleInPlace(temp, n);
    return temp;
  }

  static byte[] mulSingleDigit(const(byte)[] lhs, const(byte) n) {
    byte[] pro = lhs.dup;
    byte carry;

    for (ulong i = 0; i < pro.length; i++) {
      pro[i] *= n;
      pro[i] += carry;
      carry = pro[i] / 10;
      pro[i] %= 10;
    }

    if (carry)
      pro ~= carry;

    while (pro.length > 1 && pro[$ - 1] == 0)
      pro.length--;

    return pro;
  }

  Tuple!(byte[], byte[]) divMod(const(byte)[] lhs, const(byte)[] rhs) const {
    if (lhs.toHash() == BigInt.lastDivModArgHashes[0] && rhs.toHash() == BigInt.lastDivModArgHashes[1]) {
      return BigInt.lastDivModResult;
    }
    BigInt.lastDivModArgHashes[0] = lhs.toHash();
    BigInt.lastDivModArgHashes[1] = rhs.toHash();

    // This function could use some serious reworking and updating, for optimization purposes.
    byte[] quo = [0];
    byte[] acc;
    byte[] rem = lhs.dup;
    byte dig;
    int littleEnd;
    int i = cast(int)(lhs.length - rhs.length);

    // Needs to examine the lengths of each number and act accordingly.
    // We want this to work regardless of which side has the bigger number.
    // Modulus by a bigger number should yield the number itself.

    if (rhs == [0])
      throw new Exception("Divide by zero error.");
    if (cmpAbs(lhs, rhs) < 0) {
      BigInt.lastDivModResult[0] = quo;
      BigInt.lastDivModResult[1] = rem;
      return Tuple!(byte[], byte[])(quo, rem);
    }

    quo.length = 0;
    quo.reserve(lhs.length - rhs.length + 1);

    littleEnd = lhs.length > rhs.length ?
      cast(int)(lhs.length - rhs.length) :
      0;

    if (cmpAbs(rhs, lhs[littleEnd .. $]) > 0 && littleEnd > 0)
      littleEnd--;

    rem = lhs[littleEnd .. $].dup;

    do {
      dig = 0;
      acc = rhs.dup;

      while (cmpAbs(acc, rem) <= 0) {
        dig++;
        addAbsInPlace(acc, rhs);
      }

      subAbsInPlace(acc, rhs);
      quo ~= dig;
      subAbsInPlace(rem, acc);

      if (littleEnd > 0) {
        shiftBigInPlace(rem, 1);
        rem[0] = lhs[littleEnd - 1];
      }

    } while (littleEnd-- > 0);

    std.algorithm.reverse(quo);

    BigInt.lastDivModResult[0] = quo;
    BigInt.lastDivModResult[1] = rem;

    return Tuple!(byte[], byte[])(quo, rem);
  }

  BigInt powFast(T)(auto ref const T exp) const
  if (isIntegral!T || is(T == BigInt)) {
    const(byte)[] powInternal(T)(ref const T exp)
    if (isIntegral!T || is(T == BigInt)) {
      T oddExp = 0;
      T halfExp;

      if (exp < 0)
        throw new Exception("It's an integer library, not a fraction or floating point library. No negative exponents allowed!");
      else if (exp == 0)
        return [1];
      else if (exp == 1)
        return this.mant;

      halfExp = exp / 2;
      oddExp = exp % 2;

      static if (is(T == BigInt)) {
        auto left = powInternal(halfExp.byRef());
      } else {
        auto left = powInternal(halfExp);
      }

      const(byte)[] right = oddExp ? karatsuba(left, this.mant) : left;

      return karatsuba(left, right);
    }

    return BigInt(powInternal(exp), this.sign && exp % 2);
  }

  BigInt add(ref const BigInt rhs) const {
    BigInt sum;

    int cmpRes = cmpAbs(this.mant, rhs.mant);
    const(byte)[] big = cmpRes < 0 ? rhs.mant : this.mant;
    const(byte)[] little = cmpRes >= 0 ? rhs.mant : this.mant;

    if (this.sign == rhs.sign) {
      sum.mant = addAbs(big, little);
      sum.sign = this.sign;
    } else {
      sum.mant = subAbs(big, little);
      sum.sign = cmpRes > 0 ? this.sign : rhs.sign;
    }

    if (sum.mant[$ - 1] == 0)
      sum.sign = false;

    return sum;
  }

  BigInt sub(ref const BigInt rhs) const {
    BigInt diff = BigInt();

    int cmpRes = cmpAbs(this.mant, rhs.mant);
    auto big = cmpRes < 0 ? rhs.mant : this.mant;
    auto little = cmpRes >= 0 ? rhs.mant : this.mant;

    if (this.sign == rhs.sign) {
      diff.mant = subAbs(big, little);
      diff.sign = cmpRes > 0 ? this.sign : !this.sign;
    } else {
      diff.mant = addAbs(big, little);
      diff.sign = this.sign;
    }

    if (diff.mant[$ - 1] == 0)
      diff.sign = false;

    return diff;
  }

  ref BigInt inc() {
    if (this.sign)
      decAbs(this.mant);
    else
      incAbs(this.mant);

    if (this.mant[$ - 1] == 0)
      this.sign = false;

    return this;
  }
  
  ref BigInt dec() {
    if (this.sign)
      incAbs(this.mant);
    else if (this.mant[$ - 1] == 0) {
      incAbs(this.mant);
      this.sign = true;
    } else
      decAbs(this.mant);

    if (this.mant[$ - 1] == 0)
      this.sign = false;

    return this;
  }

  BigInt mul(ref const BigInt rhs) const {
    BigInt pro;

    if (this.mant.length < MulThreshold && rhs.mant.length < MulThreshold)
      pro.mant = mulSmall(this.mant, rhs.mant);
    else
      pro.mant = karatsuba(this.mant, rhs.mant);

    if (pro.mant[$ - 1] == 0)
      pro.sign = false;
    else
      pro.sign = this.sign != rhs.sign;

    return pro;
  }

  BigInt div(ref const BigInt rhs) const {
    BigInt quo;

    auto result = divMod(this.mant, rhs.mant);

    quo = BigInt(result[0], result[0][$-1] == 0 ? false : (this.sign != rhs.sign));

    return quo;
  }

  BigInt mod(ref const BigInt rhs) const {
    BigInt rem;
    Tuple!(byte[], byte[]) result;

    result = divMod(this.mant, rhs.mant);

    rem = BigInt(result[1], result[1][$-1] == 0 ? false : this.sign);

    return rem;
  }

public:
  static Tuple!(byte[], byte[]) lastDivModResult;
  static Tuple!(ulong, ulong) lastDivModArgHashes;

  size_t length() @property {
    return mant.length;
  }

  void length(size_t length) @property {
    mant.length = length;
    while (mant[$-1] == 0 && mant.length > 1)
      mant.length--;
  }

  this(string source) nothrow {
    bool allZeros = true;
    bool valid = true;

    if (source[0] == '-') {
      this.sign = true;
      source = source[1 .. $];
    }

    int startNdx = 0;

    //while (source[startNdx] == '0' && startNdx < source.length - 1)
      //startNdx++;

    //source = source[startNdx..$];

    for (ulong i; i < source.length; i++) {
      if (source[i] < '0' || source[i] > '9')
        valid = false;
      if (source[i] != '0')
        allZeros = false;
    }

    assert(valid == true);

    this.mant.length = source.length;
    this.mant[] = cast(byte[])(source);
    this.mant[] -= '0';
    std.algorithm.reverse(this.mant);

    if (this.mant.length >= 1 && allZeros)
      this.mant.length = 1;
  }

  this(real source) nothrow {
    this(cast(long)source);
  }

  this(int source) nothrow {
    this(long(source));
  }

  this(long source) nothrow {
    byte digit;

    this.sign = source < 0;
    this.mant.length = 0;

    source = std.math.abs(source);

    do {
      digit = source % 10UL;
      source = source / 10;
      this.mant ~= digit;
    } while (source != 0);
  }

  this(uint source) nothrow {
    this(ulong(source));
  }

  this(ulong source) nothrow {
    byte digit;

    this.sign = source < 0;
    this.mant.length = 0;

    do {
      digit = source % 10UL;
      source = source / 10;
      this.mant ~= digit;
    } while (source != 0);
  }
  
  this(const(byte)[] source, bool sign) nothrow {
    this.mant = source.dup;
    this.sign = sign;

    while (this.mant.length > 1 && this.mant[$-1] == 0)
      this.mant.length--;
  }
  
  this(const(byte)[] source) nothrow {
    this.mant = source.dup;
    this.sign = false;

    while (this.mant.length > 1 && this.mant[$-1] == 0)
      this.mant.length--;
  }

  this(this) nothrow {
    mant = mant.dup;
  }

  BigInt neg() {
    BigInt copy = this;
    copy.sign = !copy.sign;
    return copy;
  }

  BigInt abs() {
    BigInt copy = this;
    copy.sign = false;
    return copy;
  }

  auto ref BigInt opUnary(string op)() {
    static if (op == "++")
      return this.inc();
    else static if (op == "--")
      return this.dec();
    else static if (op == "-")
      return this.neg();
  }
  
  BigInt opBinary(string op, T)(auto ref const T rhs) const
  if (isIntegral!T || is(T == BigInt)) {
    static if (isIntegral!T) {
      static if (op == "^^") {
        return this.powFast(rhs);
      } else {
        return this.opBinary!op(BigInt(rhs).byRef());
      }
    } else {
      static if (op == "+") {
        return this.add(rhs.byRef());
      } else static if (op == "-") {
        return this.sub(rhs.byRef());
      } else static if (op == "*") {
        return this.mul(rhs.byRef());
      } else static if (op == "/") {
        return this.div(rhs.byRef());
      } else static if (op == "%") {
        return this.mod(rhs.byRef());
      } else static if (op == "^^") {
        return this.powFast(rhs.byRef());
      }
    }
  }

  bool opEquals(T)(const T rhs) const
  if (isIntegral!T) {
    return opEquals((BigInt(rhs)).byRef());
  }

  bool opEquals()(auto ref const BigInt rhs) const {
    if (this.mant.length == 1 && rhs.mant.length == 1 && this.mant[0] == 0 && rhs.mant[0] == 0)
      return true;

    if (this.sign != rhs.sign)
      return false;

    return this.mant == rhs.mant;
  }
  
  int opCmp(T)(T rhs) const
  if (isIntegral!T) {
    return this.opCmp((BigInt(rhs)).byRef());
  }

  int opCmp(ref const BigInt rhs) const {
    int res;

    if (this.mant.length == 1 && rhs.mant.length == 1 && this.mant[0] == 0 && rhs.mant[0] == 0) {
      res = 0;
    } else if (this.sign != rhs.sign) {
      if (this.sign == true)
        res = -1;
      else
        res = 1;
    } else {
      res = cmpAbs(this.mant, rhs.mant);

      if (this.sign == true) {
        res = -res;
      }
    }

    return res;
  }

  void opAssign(T)(T rhs)
  if (isIntegral!T || isSomeString!T) {
    this = BigInt(rhs);
  }

  void opOpAssign(string op, T)(T rhs)
  if (isIntegral!T || isSomeString!T) {
    BigInt value = BigInt(rhs);
    opOpAssign!op(value);
  }

  void opOpAssign(string op, T)(T rhs)
  if (is(T == BigInt)) {
      static if (op == "+")
        this = this.add(rhs);
      else static if (op == "-")
        this = this.sub(rhs);
      else static if (op == "*")
        this = this.mul(rhs);
      else static if (op == "/")
        this = this.div(rhs);
      else static if (op == "%")
        this = this.mod(rhs);
  }

  T opCast(T)() {
    static if (is(T == bool)) {
      if (this.mant.length == 1 && this.mant[0] == 0)
        return false;
      else
        return true;
    } else static if (is(T == long) || is(T == ulong)) {
      T result = 0;

      foreach (i, n; this.mant)
        result += n * 10 ^^ i;

      if (this.sign)
        result = -result;

      return result;
    } else static if (isFloatingPoint!T) {
      return (T(cast(long)this));
    }
  }

  typeof(this) opSlice(size_t i, size_t j) {
    return BigInt(this.mant[i..j], this.sign);
  }

  size_t opDollar() {
    return this.mant.length;
  }

/*
  string toString() const {
    byte[] str = this.mant.dup;
    str[] += '0';

    if (sign) {
      str ~= '-';
    }

    std.algorithm.reverse(str);

    return cast(string)str;
  }
*/
  
  void toString(W)(ref W writer, scope const ref FormatSpec!char fmt) const
  if (isOutputRange!(W, char)) {
    if (fmt.spec == 's') {
      if (sign)
        writer.put('-');
      for (int i = 1; i <= mant.length; i++) {
        writer.put(cast(char)(mant[$ - i] + '0'));
      }
    } else if (fmt.spec == 'd') {
      if (sign)
        writer.put('-');
      for (int i = 1; i <= mant.length; i++) {
        writer.put(cast(char)(mant[$ - i] + '0'));
      }
    }
  }

  immutable(char)[] digitString() const {
    //return this.mant.retro.map!(d => cast(immutable(char))(d + '0')).array();
    byte[] digits = this.mant.dup;
    digits[] += '0';
    std.algorithm.reverse(digits);
    return cast(immutable(char)[])digits;
  }

  byte[] digitBytes() const {
    //return this.mant.retro.map!(d => cast(immutable(char))(d + '0')).array();
    byte[] digitBytes = this.mant.dup;
    //digits[] += '0';
    std.algorithm.reverse(digitBytes);
    return digitBytes;
  }

  uint[] toDigits() const {
    return this.mant.retro.map!(a => cast(uint)a).array();
  }

  typeof(this) reverse() {
    byte[] resmant = this.mant.dup;
    std.algorithm.reverse(resmant);
    return BigInt(resmant);
  }

  @property ulong digitsLength() const {
    return this.mant.length;
  }

  mixin RvalueRef;
}

ulong countDigits(BigInt source) {
  return source.mant.length;
}

private:

byte[] rbytes(string value) {
  //return value.retro.map!(a => cast(byte)(a - '0')).array();
  char[] result = value.dup;
  result[] -= '0';
  std.algorithm.reverse(result);
  return cast(byte[])result;
}
unittest {
  //writeln("rbytes unittest");
  assert("2374".rbytes() == [4, 7, 3, 2]);
  //writeln("rbytes unittest passed");
}

string rstr(const(byte)[] value) {
  //return value.retro.map!(a => cast(immutable(char))(a + '0')).array();
  byte[] result = value.dup;
  result[] += '0';
  std.algorithm.reverse(result);
  return cast(string)result;
}
unittest {
  //writeln("rstr unittest");
  assert([4, 7, 3, 2].rstr() == "2374");
  //writeln("rstr unittest passed");
}

ulong toHash(T)(auto ref T value) {
  return typeid(T).getHash(&value);
}