Source code in the hoon programming language

::
::  A library for producing Rhonda numbers and testing if numbers are Rhonda.
::
::    A number is Rhonda if the product of its digits of in base b equals 
::    the product of the base b and the sum of its prime factors.
::    see also: https://mathworld.wolfram.com/RhondaNumber.html
::
=<
::
|%
::  +check: test whether the number n is Rhonda to base b
::
++  check
  |=  [b=@ud n=@ud]
  ^-  ?
  ~_  leaf+"base b must be >= 2"
  ?>  (gte b 2)
  ~_  leaf+"candidate number n must be >= 2"
  ?>  (gte n 2)
  ::
  .=  (roll (base-digits b n) mul)
  %+  mul
    b
  (roll (prime-factors n) add)
::  +series: produce the first n numbers which are Rhonda in base b
::
::    produce ~ if base b has no Rhonda numbers
::
++  series
  |=  [b=@ud n=@ud]
  ^-  (list @ud)
  ~_  leaf+"base b must be >= 2"
  ?>  (gte b 2)
  ::
  ?:  =((prime-factors b) ~[b])
    ~
  =/  candidate=@ud  2
  =+  rhondas=*(list @ud)
  |-
  ?:  =(n 0)
    (flop rhondas)
  =/  is-rhonda=?  (check b candidate)
  %=  $
    rhondas    ?:(is-rhonda [candidate rhondas] rhondas)
    n          ?:(is-rhonda (dec n) n)
    candidate  +(candidate)
  ==
--
::
|%
::  +base-digits: produce a list of the digits of n represented in base b
::
::    This arm has two behaviors which may be at first surprising, but do not
::    matter for the purposes of the ++check and ++series arms, and allow for
::    some simplifications to its implementation.
::    - crashes on n=0
::    - orders the list of digits with least significant digits first
::
::    ex: (base-digits 4 10.206) produces ~[2 3 1 3 3 1 2]
::
++  base-digits
  |=  [b=@ud n=@ud]
  ^-  (list @ud)
  ?>  (gte b 2)
  ?<  =(n 0)
  ::
  |-
  ?:  =(n 0)
    ~
  :-  (mod n b)
  $(n (div n b))
::  +prime-factors: produce a list of the prime factors of n
::    
::    by trial division
::    n must be >= 2
::    if n is prime, produce ~[n]
::    ex: (prime-factors 10.206) produces ~[7 3 3 3 3 3 3 2]
::
++  prime-factors
  |=  [n=@ud]
  ^-  (list @ud)
  ?>  (gte n 2)
  ::
  =+  factors=*(list @ud)
  =/  wheel  new-wheel
  ::  test candidates as produced by the wheel, not exceeding sqrt(n) 
  ::
  |-
  =^  candidate  wheel  (next:wheel)
  ?.  (lte (mul candidate candidate) n)
    ?:((gth n 1) [n factors] factors)
  |-
  ?:  =((mod n candidate) 0)
    ::  repeat the prime factor as many times as possible
    ::
    $(factors [candidate factors], n (div n candidate))
  ^$
::  +new-wheel: a door for generating numbers that may be prime
::
::    This uses wheel factorization with a basis of {2, 3, 5} to limit the
::    number of composites produced. It produces numbers in increasing order
::    starting from 2.
::
++  new-wheel
  =/  fixed=(list @ud)  ~[2 3 5 7]
  =/  skips=(list @ud)  ~[4 2 4 2 4 6 2 6]
  =/  lent-fixed=@ud  (lent fixed)
  =/  lent-skips=@ud  (lent skips)
  ::
  |_  [current=@ud fixed-i=@ud skips-i=@ud]
  ::  +next: produce the next number and the new wheel state
  ::
  ++  next
    |.
    ::  Exhaust the numbers in fixed. Then calculate successive values by
    ::  cycling through skips and increasing from the previous number by
    ::  the current skip-value.
    ::
    =/  fixed-done=?  =(fixed-i lent-fixed)
    =/  next-fixed-i  ?:(fixed-done fixed-i +(fixed-i))
    =/  next-skips-i  ?:(fixed-done (mod +(skips-i) lent-skips) skips-i)
    =/  next
    ?.  fixed-done
      (snag fixed-i fixed)
    (add current (snag skips-i skips))
    :-  next
    +.$(current next, fixed-i next-fixed-i, skips-i next-skips-i)
  --
--

/+  *rhonda
:-  %say
|=  [* [base=@ud many=@ud ~] ~]
:-  %noun
(series base many)

|%
++  check
  |=  [n=@ud base=@ud]
  ::  if base is prime, automatic no
  ::
  ?:  =((~(gut by (prime-map +(base))) base 0) 0)
    %.n
  ::  if not multiply the digits and compare to base x sum of factors
  ::
  ?:  =((roll (digits [base n]) mul) (mul base (roll (factor n) add)))
    %.y
  %.n
++  series
  |=  [base=@ud many=@ud]
  =/  rhondas  *(list @ud)
  ?:  =((~(gut by (prime-map +(base))) base 0) 0)
    rhondas
  =/  itr  1
  |-
  ?:  =((lent rhondas) many)
    (flop rhondas)
  ?:  =((check itr base) %.n)
    $(itr +(itr))
  $(rhondas [itr rhondas], itr +(itr))
::  digits: gives the list of digits of a number in a base
::
::    We strip digits least to most significant.
::    The least significant digit (lsd) of n in base b is just n mod b.
::    Subtract the lsd, divide by b, and repeat.
::    To know when to stop, we need to know how many digits there are.
++  digits
  |=  [base=@ud num=@ud]
  ^-  (list @ud)
  |-
  =/  modulus=@ud  (mod num base)
  ?:  =((num-digits base num) 1)
    ~[modulus]
  [modulus $(num (div (sub num modulus) base))]
::  num-digits: gives the number of digits of a number in a base
::
::    Simple idea: k is the number of digits of n in base b if and
::    only if k is the smallest number such that b^k > n.
++  num-digits
  |=  [base=@ud num=@ud]
  ^-  @ud
  =/  digits=@ud  1
  |-
  ?:  (gth (pow base digits) num)
    digits
  $(digits +(digits))
::  factor: produce a list of prime factors
::
::    The idea is to identify "small factors" of n, i.e. prime factors less than
::    the square root. We then divide n by these factors to reduce the
::    magnitude of n. It's easy to argue that after this is done, we obtain 1
::    or the largest prime factor.
::
++  factor
  |=  n=@ud
  ^-  (list @ud)
  ?:  ?|(=(n 0) =(n 1))
    ~[n]
  =/  factorization  *(list @ud)
  ::  produce primes less than or equal to root n
  ::
  =/  root  (sqrt n)
  =/  primes  (prime-map +(root))
  ::  itr = iterate; we want to iterate through the primes less than root n
  ::
  =/  itr  2
  |-
  ?:  =(itr +(root))
  ::  if n is now 1 we're done
  ::
    ?:  =(n 1)
      factorization
    ::  otherwise it's now the original n's largest primes factor
    ::
    [n factorization]
  ::  if itr not prime move on
  ::
  ?:  =((~(gut by primes) itr 0) 1)
    $(itr +(itr))
  ::  if it is prime, divide out by the highest power that divides num
  ::
  ?:  =((mod n itr) 0)
    $(n (div n itr), factorization [itr factorization])
  ::  once done, move to next prime
  ::
  $(itr +(itr))
::  sqrt: gives the integer square root of a number
::
::    It's based on an algorithm that predates the Greeks:
::    To find the square root of A, think of A as an area.
::    Guess the side of the square x. Compute the other side y = A/x.
::    If x is an over/underestimate then y is an under/overestimate.
::    So (x+y)/2 is the average of an over and underestimate, thus better than x.
::    Repeatedly doing x --> (x + A/x)/2 converges to sqrt(A).
::
::    This algorithm is the same but with integer valued operations.
::    The algorithm either converges to the integer square root and repeats,
::    or gets trapped in a two-cycle of adjacent integers.
::    In the latter case, the smaller number is the answer.
::
++  sqrt
  |=  n=@ud
  =/  guess=@ud  1
  |-
  =/  new-guess  (div (add guess (div n guess)) 2)
  ::  sequence stabilizes
  ::
  ?:  =(guess new-guess)
    guess
  ::  sequence is trapped in 2-cycle
  ::
  ?:  =(guess +(new-guess))
    new-guess
  ?:  =(new-guess +(guess))
    guess
  $(guess new-guess)
::  prime-map: (effectively) produces primes less than a given input
::
::    This is the sieve of Eratosthenes to produce primes less than n.
::    I used a map because it had much faster performance than a list.
::    Any key in the map is a non-prime. The value 1 indicates "false."
::    I.e. it's not a prime.
++  prime-map
  |=  n=@ud
  ^-  (map @ud @ud)
  =/  prime-map  `(map @ud @ud)`(my ~[[0 1] [1 1]])
  ::  start sieving with 2
  ::
  =/  sieve  2
  |-
  ::  if sieve is too large to be a factor we're done
  ::
  ?:  (gte (mul sieve sieve) n)
    prime-map
  ::  if not too large but not prime, move on
  ::
  ?:  =((~(gut by prime-map) sieve 0) 1)
    $(sieve +(sieve))
  ::  sequence: explanation
  ::
  ::    If s is the sieve number, we start sieving multiples
  ::    of s at s^2 in sequence: s^2, s^2 + s, s^2 + 2s, ...
  ::    We start at s^2 because any number smaller than s^2
  ::    has prime factors less than s and would have been
  ::    eliminated earlier in the sieving process.
  ::
  =/  sequence  (mul sieve sieve)
  |-
  ::  done sieving with s once sequence is past n
  ::
  ?:  (gte sequence n)
    ^$(sieve +(sieve))
  ::  if sequence position is known not prime we move on
  ::
  ?:  =((~(gut by prime-map) sequence 0) 1)
    $(sequence (add sequence sieve))
  ::  otherwise we mark position of sequence as not prime and move on
  ::
  $(prime-map (~(put by prime-map) sequence 1), sequence (add sequence sieve))
--