How to resolve the algorithm Hamming numbers step by step in the Chapel programming language

Published on 12 May 2024 09:40 PM

How to resolve the algorithm Hamming numbers step by step in the Chapel programming language

Table of Contents

Problem Statement

Hamming numbers are numbers of the form   Hamming numbers   are also known as   ugly numbers   and also   5-smooth numbers   (numbers whose prime divisors are less or equal to 5).

Generate the sequence of Hamming numbers, in increasing order.   In particular:

Let's start with the solution:

Step by Step solution about How to resolve the algorithm Hamming numbers step by step in the Chapel programming language

Source code in the chapel programming language

use BigInteger; use Time;

// Chapel doesn't have closure functions that can capture variables from
// outside scope, so we use a class to emulate them for this special case;
// the member fields mult, mrglst, and mltlst, emulate "captured" variables
// that would normally be captured by the `next` continuation closure...
class HammingsList {
    const head: bigint;
    const mult: uint(8);
    var mrglst: shared HammingsList?;
    var mltlst: shared HammingsList?;
    var tail: shared HammingsList? = nil;
    proc init(hd: bigint, mlt: uint(8), mrgl: shared HammingsList?,
                                        mltl: shared HammingsList?) {
        head = hd; mult = mlt; mrglst = mrgl; mltlst = mltl; }
    proc next(): shared HammingsList {
        if tail != nil then return tail: shared HammingsList;
        const nhd: bigint = mltlst!.head * mult;
        if mrglst == nil then {
            tail = new shared HammingsList(nhd, mult,
                                           nil: shared HammingsList?,
                                           nil: shared HammingsList?);
            mltlst = mltlst!.next();
            tail!.mltlst <=> mltlst;
        }
        else {
            if mrglst!.head < nhd then {
                tail = new shared HammingsList(mrglst!.head, mult,
                                               nil: shared HammingsList?,
                                               nil: shared HammingsList?);
                mrglst = mrglst!.next(); mrglst <=> tail!.mrglst;
                mltlst <=> tail!.mltlst;
            }
            else {
                tail = new shared HammingsList(nhd, mult,
                                               nil: shared HammingsList?,
                                               nil: shared HammingsList?);
                mltlst = mltlst!.next(); mltlst <=> tail!.mltlst;
                mrglst <=> tail!.mrglst;
            }
        }
        return tail: shared HammingsList;
    }
}

proc u(n: uint(8), s: shared HammingsList?): shared HammingsList {
    var r = new shared HammingsList(1: bigint, n, s,
                                    nil: shared HammingsList?);
    r.mltlst = r; // lazy recursion!
    return r.next();
}

iter hammings(): bigint {
    var nxt: shared HammingsList? = nil: shared HammingsList?;
    const mlts: [ 0 .. 2 ] int = [ 5, 3, 2 ];
    for m in mlts do nxt = u(m: uint(8), nxt);
    yield 1 : bigint;
    while true { yield nxt!.head; nxt = nxt!.next(); }
}

write("The first 20 Hamming numbers are: ");
var cnt: int = 0;
for h in hammings() { write(" ", h); cnt += 1; if cnt >= 20 then break; }
write(".\nThe 1691st Hamming number is ");
cnt = 0;
for h in hammings() { cnt += 1; if cnt < 1691 then continue; write(h); break; }
writeln(".\nThe millionth Hamming number is ");
var timer: Timer; timer.start(); cnt = 0;
for h in hammings() { cnt += 1; if cnt < 1000000 then continue; write(h); break; }
timer.stop(); writeln(".\nThis last took ",
                      timer.elapsed(TimeUnits.milliseconds), " milliseconds.");


use BigInteger; use Time;

iter nodupsHamming(): bigint {
  var s2dom = { 0 .. 1023 }; var s2: [s2dom] bigint; // init so can double!
  var s3dom = { 0 .. 1023 }; var s3: [s3dom] bigint; // init so can double!
  s2[0] = 1: bigint; s3[0] = 3: bigint;
  var x5 = 5: bigint; var mrg = 3: bigint;
  var s2hdi, s2tli, s3hdi, s3tli: int;
  while true {
    s2tli += 1;
    if s2hdi + s2hdi >= s2tli { // move in place to avoid allocation!
      s2[0 .. s2tli - s2hdi - 1] = s2[s2hdi .. s2tli - 1];
      s2tli -= s2hdi; s2hdi = 0; }
    const s2sz = s2.size;
    if s2tli >= s2sz then s2dom = { 0 .. s2sz + s2sz - 1 };
    var rslt: bigint; const s2hd = s2[s2hdi];
    if s2hd < mrg { rslt = s2hd; s2hdi += 1; }
    else {
      s3tli += 1;
      if s3hdi + s3hdi >= s2tli { // move in place to avoid allocation!
        s3[0 .. s3tli - s3hdi - 1] = s3[s3hdi .. s3tli - 1];
        s3tli -= s3hdi; s3hdi = 0; }
      const s3sz = s3.size;
      if s3tli >= s3sz then s3dom = { 0 .. s3sz + s3sz - 1 };
      rslt = mrg; s3[s3tli] = rslt * 3;
      s3hdi += 1; const s3hd = s3[s3hdi];
      if s3hd < x5 { mrg = s3hd; }
      else { mrg = x5; x5 = x5 * 5; s3hdi -= 1; }
    }
    s2[s2tli] = rslt * 2;
    yield rslt;
  }
}

// test it...
write("The first 20 hamming numbers are: ");
var cnt = 0: uint(64);
for h in nodupsHamming() {
  if cnt >= 20 then break; cnt += 1; write(" ", h); }

write("\nThe 1691st hamming number is "); cnt = 1;
for h in nodupsHamming() {
  if cnt >= 1691 { writeln(h); break; } cnt += 1; }

write("The millionth hamming number is ");
var timer: Timer; cnt = 1;
timer.start(); var rslt: bigint;
for h in nodupsHamming() {
  if cnt >= 1000000 { rslt = h; break; } cnt += 1; }
timer.stop();
write(rslt);
writeln(".\nThis last took ",
        timer.elapsed(TimeUnits.milliseconds), " milliseconds.");


use BigInteger; use Math; use Time;

config const nth: uint(64) = 1000000;

const lb2 = 1: real(64); // log base 2 of 2!
const lb3 = log2(3: real(64)); const lb5 = log2(5: real(64));
record LogRep {
  var lg: real(64); var x2: uint(32);
  var x3: uint(32); var x5: uint(32);
  inline proc mul2(): LogRep {
    return new LogRep(this.lg + lb2, this.x2 + 1, this.x3, this.x5); }
  inline proc mul3(): LogRep {
    return new LogRep(this.lg + lb3, this.x2, this.x3 + 1, this.x5); }
  inline proc mul5(): LogRep {
    return new LogRep(this.lg + lb5, this.x2, this.x3, this.x5 + 1); }
  proc lr2bigint(): bigint {
    proc xpnd(bs: uint, v: uint(32)): bigint {
      var rslt = 1: bigint; var bsm = bs: bigint; var vm = v: uint;
      while vm > 0 { if vm & 1 then rslt *= bsm; bsm *= bsm; vm >>= 1; }
      return rslt;
    }
    return xpnd(2: uint, this.x2) *
             xpnd(3: uint, this.x3) * xpnd(5: uint, this.x5);
  }
  proc writeThis(lr) throws {
    lr <~> this.lr2bigint();
  }
}
operator <(const ref a: LogRep, const ref b: LogRep): bool { return a.lg < b.lg; }
const one = new LogRep(0, 0, 0, 0);

iter nodupsHammingLog(): LogRep {
  var s2dom = { 0 .. 1023 }; var s2: [s2dom] LogRep; // init so can double!
  var s3dom = { 0 .. 1023 }; var s3: [s3dom] LogRep; // init so can double!
  s2[0] = one; s3[0] = one.mul3();
  var x5 = one.mul5(); var mrg = one.mul3();
  var s2hdi, s2tli, s3hdi, s3tli: int;
  while true {
    s2tli += 1;
    if s2hdi + s2hdi >= s2tli { // move in place to avoid allocation!
      s2[0 .. s2tli - s2hdi - 1] = s2[s2hdi .. s2tli - 1];
      s2tli -= s2hdi; s2hdi = 0; }
    const s2sz = s2.size;
    if s2tli >= s2sz then s2dom = { 0 .. s2sz + s2sz - 1 };
    var rslt: LogRep; const s2hd = s2[s2hdi];
    if s2hd.lg < mrg.lg { rslt = s2hd; s2hdi += 1; }
    else {
      s3tli += 1;
      if s3hdi + s3hdi >= s2tli { // move in place to avoid allocation!
        s3[0 .. s3tli - s3hdi - 1] = s3[s3hdi .. s3tli - 1];
        s3tli -= s3hdi; s3hdi = 0; }
      const s3sz = s3.size;
      if s3tli >= s3sz then s3dom = { 0 .. s3sz + s3sz - 1 };
      rslt = mrg; s3[s3tli] = mrg.mul3(); s3hdi += 1;
      const s3hd = s3[s3hdi];
      if s3hd.lg < x5.lg { mrg = s3hd; }
      else { mrg = x5; x5 = x5.mul5(); s3hdi -= 1; }
    }
    s2[s2tli] = rslt.mul2();
    yield rslt;
  }
}
 
// test it...
write("The first 20 hamming numbers are: ");
var cnt = 0: uint(64);
for h in nodupsHammingLog() {
  if cnt >= 20 then break; cnt += 1; write(" ", h); }

write("\nThe 1691st hamming number is "); cnt = 1;
for h in nodupsHammingLog() {
  if cnt >= 1691 { writeln(h); break; } cnt += 1; }

write("The ", nth, "th hamming number is ");
var timer: Timer; cnt = 1;
timer.start(); var rslt: LogRep;
for h in nodupsHammingLog() {
  if cnt >= nth { rslt = h; break; } cnt += 1; }
timer.stop();
write(rslt);
writeln(".\nThis last took ",
        timer.elapsed(TimeUnits.milliseconds), " milliseconds.");


use BigInteger; use Math; use Sort; use Time;

config const nth = 1000000: uint(64);

type TriVal = 3*uint(32);

proc trival2bigint(x: TriVal): bigint {
  proc xpnd(bs: uint, v: uint(32)): bigint {
    var rslt = 1: bigint; var bsm = bs: bigint; var vm = v: uint;
    while vm > 0 { if vm & 1 then rslt *= bsm; bsm *= bsm; vm >>= 1; }
    return rslt;
  }
  const (x2, x3, x5) = x;
  return xpnd(2: uint, x2) * xpnd(3: uint, x3) * xpnd(5: uint, x5);
}

proc nthHamming(n: uint(64)): TriVal {
  if n < 1 {
    writeln("nthHamming - argument must be at least one!"); exit(1); }
  if n < 2 then return (0: uint(32), 0: uint(32), 0: uint(32)); // TriVal for 1

  type LogRep = (real(64), uint(32), uint(32), uint(32));
  record Comparator {} // used for sorting in reverse order!
  proc Comparator.compare(a: LogRep, b: LogRep): real(64) {
    return b[0] - a[0]; }
  var logrepComp: Comparator;

  const lb3 = log2(3.0: real(64)); const lb5 = log2(5.0: real(64));
  const fctr = 6.0: real(64) * lb3 * lb5;
  const crctn = log2(sqrt(30.0: real(64))); // log base 2 of sqrt 30
  // from Wikipedia Regular Numbers formula...
  const lgest = (fctr * n: real(64))**(1.0: real(64) / 3.0: real(64)) - crctn;
  const frctn = if n < 1000000000 then 0.509: real(64) else 0.105: real(64);
  const lghi = (fctr * (n: real(64) + frctn * lgest))**
                 (1.0: real(64) / 3.0: real(64)) - crctn;
  const lglo = 2.0: real(64) * lgest - lghi; // lower limit of the upper "band"
  var count = 0: uint(64); // need to use extended precision, might go over
  var bndi = 0; var dombnd = { 0 .. bndi }; // one value so doubling size works!
  var bnd: [dombnd] LogRep; const klmt = (lghi / lb5): uint(32);
  for k in 0 .. klmt { // i, j, k values can be just uint(32) values!
    const p = k: real(64) * lb5; const jlmt = ((lghi - p) / lb3): uint(32);
    for j in 0 .. jlmt {
      const q = p + j: real(64) * lb3;
      const ir = lghi - q; const lg = q + floor(ir); // current log value (est)
      count += ir: uint(64) + 1;
      if lg >= lglo {
        const sz = dombnd.size; if bndi >= sz then dombnd = { 0..sz + sz - 1 };
        bnd[bndi] = (lg, ir: uint(32), j, k); bndi += 1;
      }
    }
  }
  if n > count {
    writeln("nth_hamming: band high estimate is too low!"); exit(1); }
  dombnd = { 0 .. bndi - 1 }; const ndx = (count - n): int;
  if ndx >= dombnd.size {
    writeln("nth_hamming: band low estimate is too high!"); exit(1); }
  sort(bnd, comparator = logrepComp); // descending order leaves zeros at end!

  const rslt = bnd[ndx]; return (rslt[1], rslt[2], rslt[3]);
}

// test it...
write("The first 20 Hamming numbers are: ");
for i in 1 .. 20 do write(" ", trival2bigint(nthHamming(i: uint(64))));

writeln("\nThe 1691st hamming number is ",
        trival2bigint(nthHamming(1691: uint(64))));

var timer: Timer;
timer.start();
const answr = nthHamming(nth);
timer.stop();
write("The ", nth, "th Hamming number is 2**",
      answr[0], " * 3**", answr[1], " * 5**", answr[2]);
const lgrslt = (answr[0]: real(64) + answr[1]: real(64) * log2(3: real(64)) +
                answr[2]: real(64) * log2(5: real(64))) * log10(2: real(64));
const whl = lgrslt: uint(64); const frac = lgrslt - whl: real(64);
write(",\nwhich is approximately ", 10: real(64)**frac, "E+", whl);
const bganswr = trival2bigint(answr);
const answrstr = bganswr: string; const asz = answrstr.size;
writeln(" and has ", asz, " digits.");
if asz <= 2000 then write("Can be printed as:  ", answrstr);
else write("It's too long to print");
writeln("!\nThis last took ",
        timer.elapsed(TimeUnits.milliseconds), " milliseconds.");


use BigInteger; use Math; use Sort; use Time;

config const nth = 1000000: uint(64);

type TriVal = 3*uint(32);

proc trival2bigint(x: TriVal): bigint {
  proc xpnd(bs: uint, v: uint(32)): bigint {
    var rslt = 1: bigint; var bsm = bs: bigint; var vm = v: uint;
    while vm > 0 { if vm & 1 then rslt *= bsm; bsm *= bsm; vm >>= 1; }
    return rslt;
  }
  const (x2, x3, x5) = x;
  return xpnd(2: uint, x2) * xpnd(3: uint, x3) * xpnd(5: uint, x5);
}

proc nthHamming(n: uint(64)): TriVal {
  if n < 1 {
    writeln("nthHamming - argument must be at least one!"); exit(1); }
  if n < 2 then return (0: uint(32), 0: uint(32), 0: uint(32)); // TriVal for 1

  type LogRep = (bigint, uint(32), uint(32), uint(32));
  record Comparator {} // used for sorting in reverse order!
  proc Comparator.compare(a: LogRep, b: LogRep): int {
    return (b[0] - a[0]): int; }
  var logrepComp: Comparator;

  const lb3 = log2(3.0: real(64)); const lb5 = log2(5.0: real(64));
  const bglb2 = "1267650600228229401496703205376": bigint;
  const bglb3 = "2009178665378409109047848542368": bigint;
  const bglb5 = "2943393543170754072109742145491": bigint;
  const fctr = 6.0: real(64) * lb3 * lb5;
  const crctn = log2(sqrt(30.0: real(64))); // log base 2 of sqrt 30
  // from Wikipedia Regular Numbers formula...
  const lgest = (fctr * n: real(64))**(1.0: real(64) / 3.0: real(64)) - crctn;
  const frctn = if n < 1000000000 then 0.509: real(64) else 0.105: real(64);
  const lghi = (fctr * (n: real(64) + frctn * lgest))**
                 (1.0: real(64) / 3.0: real(64)) - crctn;
  const lglo = 2.0: real(64) * lgest - lghi; // lower limit of the upper "band"
  var count = 0: uint(64); // need to use extended precision, might go over
  var bndi = 0; var dombnd = { 0 .. bndi }; // one value so doubling size works!
  var bnd: [dombnd] LogRep; const klmt = (lghi / lb5): uint(32);
  for k in 0 .. klmt { // i, j, k values can be just uint(32) values!
    const p = k: real(64) * lb5; const jlmt = ((lghi - p) / lb3): uint(32);
    for j in 0 .. jlmt {
      const q = p + j: real(64) * lb3;
      const ir = lghi - q; const lg = q + floor(ir); // current log value (est)
      count += ir: uint(64) + 1;
      if lg >= lglo {
        const sz = dombnd.size; if bndi >= sz then dombnd = { 0..sz + sz - 1 };
        const bglg =
          bglb2 * ir: int(64) + bglb3 * j: int(64) + bglb5 * k: int(64);
        bnd[bndi] = (bglg, ir: uint(32), j, k); bndi += 1;
      }
    }
  }
  if n > count {
    writeln("nth_hamming: band high estimate is too low!"); exit(1); }
  dombnd = { 0 .. bndi - 1 }; const ndx = (count - n): int;
  if ndx >= dombnd.size {
    writeln("nth_hamming: band low estimate is too high!"); exit(1); }
  sort(bnd, comparator = logrepComp); // descending order leaves zeros at end!

  const rslt = bnd[ndx]; return (rslt[1], rslt[2], rslt[3]);
}

// test it...
write("The first 20 Hamming numbers are: ");
for i in 1 .. 20 do write(" ", trival2bigint(nthHamming(i: uint(64))));

writeln("\nThe 1691st hamming number is ",
        trival2bigint(nthHamming(1691: uint(64))));

var timer: Timer;
timer.start();
const answr = nthHamming(nth);
timer.stop();
write("The ", nth, "th Hamming number is 2**",
      answr[0], " * 3**", answr[1], " * 5**", answr[2]);
const lgrslt = (answr[0]: real(64) + answr[1]: real(64) * log2(3: real(64)) +
                answr[2]: real(64) * log2(5: real(64))) * log10(2: real(64));
const whl = lgrslt: uint(64); const frac = lgrslt - whl: real(64);
write(",\nwhich is approximately ", 10: real(64)**frac, "E+", whl);
const bganswr = trival2bigint(answr);
const answrstr = bganswr: string; const asz = answrstr.size;
writeln(" and has ", asz, " digits.");
if asz <= 2000 then write("Can be printed as:  ", answrstr);
else write("It's too long to print");
writeln("!\nThis last took ",
        timer.elapsed(TimeUnits.milliseconds), " milliseconds.");


  

You may also check:How to resolve the algorithm Zero to the zero power step by step in the Haskell programming language
You may also check:How to resolve the algorithm Look-and-say sequence step by step in the K programming language
You may also check:How to resolve the algorithm Rosetta Code/Count examples step by step in the Lasso programming language
You may also check:How to resolve the algorithm Binary strings step by step in the VBA programming language
You may also check:How to resolve the algorithm Averages/Arithmetic mean step by step in the PowerShell programming language