Smith numbers

From Rosetta Code
Revision as of 01:19, 25 April 2016 by rosettacode>Gerard Schildberger (→‎{{header|REXX}}: fixed a closing HTML lang tag.)
Smith numbers is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Smith numbers are numbers such that the sum of the digits of the integers that make up that number is the same as the sum of digits of its prime factors excluding 1. The primes are also excluded as they, naturally, satisfy this condition!

Example: 166 => 2 x 83 => 2 + 8 + 3 = 1 + 6 + 6 = 13.

The task: write a program to find all Smith numbers below 10000.


C++

<lang cpp>

  1. include <iostream>
  2. include <vector>
  3. include <iomanip>

void primeFactors( unsigned n, std::vector<unsigned>& r ) {

   int f = 2; if( n == 1 ) r.push_back( 1 );
   else {
       while( true ) {
           if( !( n % f ) ) {
               r.push_back( f );
               n /= f; if( n == 1 ) return;
           }
           else f++;
       }
   }

} unsigned sumDigits( unsigned n ) {

   unsigned sum = 0, m;
   while( n ) {
       m = n % 10; sum += m;
       n -= m; n /= 10;
   }
   return sum;

} unsigned sumDigits( std::vector<unsigned>& v ) {

   unsigned sum = 0;
   for( std::vector<unsigned>::iterator i = v.begin(); i != v.end(); i++ ) {
       sum += sumDigits( *i );
   }
   return sum;

} void listAllSmithNumbers( unsigned n ) {

   std::vector<unsigned> pf;
   for( unsigned i = 4; i < n; i++ ) {
       primeFactors( i, pf ); if( pf.size() < 2 ) continue;
       if( sumDigits( i ) == sumDigits( pf ) )
           std::cout << std::setw( 4 ) << i << " ";
       pf.clear();
   }
   std::cout << "\n\n";

} int main( int argc, char* argv[] ) {

   listAllSmithNumbers( 10000 );
   return 0;

} </lang>

Output:
   4   22   27   58   85   94  121  166  202  265  274  319  346  355  378  382
 391  438  454  483  517  526  535  562  576  627  634  636  645  663  666  690
...
9301 9330 9346 9355 9382 9386 9387 9396 9427 9483 9535 9571 9598 9633 9634 9639 
9648 9657 9684 9708 9717 9735 9742 9760 9778 9843 9849 9861 9880 9895 9975 9985

J

Implementation: <lang J>digits=: 10&#.inv sumdig=: +/@,@digits notprime=: -.@(1&p:) smith=: #~ notprime * (=&sumdig q:) every</lang>

Task example: <lang J> #smith }.i.10000 376

  q:376

2 2 2 47

  8 47$smith }.i.10000
  4   22   27   58   85   94  121  166  202  265  274  319  346  355  378  382  391  438  454  483  517  526  535  562  576  588  627  634  636  645  648  654  663  666  690  706  728  729  762  778  825  852  861  895  913  915  922
958  985 1086 1111 1165 1219 1255 1282 1284 1376 1449 1507 1581 1626 1633 1642 1678 1736 1755 1776 1795 1822 1842 1858 1872 1881 1894 1903 1908 1921 1935 1952 1962 1966 2038 2067 2079 2155 2173 2182 2218 2227 2265 2286 2326 2362 2366

2373 2409 2434 2461 2475 2484 2515 2556 2576 2578 2583 2605 2614 2679 2688 2722 2745 2751 2785 2839 2888 2902 2911 2934 2944 2958 2964 2965 2970 2974 3046 3091 3138 3168 3174 3226 3246 3258 3294 3345 3366 3390 3442 3505 3564 3595 3615 3622 3649 3663 3690 3694 3802 3852 3864 3865 3930 3946 3973 4054 4126 4162 4173 4185 4189 4191 4198 4209 4279 4306 4369 4414 4428 4464 4472 4557 4592 4594 4702 4743 4765 4788 4794 4832 4855 4880 4918 4954 4959 4960 4974 4981 5062 5071 5088 5098 5172 5242 5248 5253 5269 5298 5305 5386 5388 5397 5422 5458 5485 5526 5539 5602 5638 5642 5674 5772 5818 5854 5874 5915 5926 5935 5936 5946 5998 6036 6054 6084 6096 6115 6171 6178 6187 6188 6252 6259 6295 6315 6344 6385 6439 6457 6502 6531 6567 6583 6585 6603 6684 6693 6702 6718 6760 6816 6835 6855 6880 6934 6981 7026 7051 7062 7068 7078 7089 7119 7136 7186 7195 7227 7249 7287 7339 7402 7438 7447 7465 7503 7627 7674 7683 7695 7712 7726 7762 7764 7782 7784 7809 7824 7834 7915 7952 7978 8005 8014 8023 8073 8077 8095 8149 8154 8158 8185 8196 8253 8257 8277 8307 8347 8372 8412 8421 8466 8518 8545 8568 8628 8653 8680 8736 8754 8766 8790 8792 8851 8864 8874 8883 8901 8914 9015 9031 9036 9094 9166 9184 9193 9229 9274 9276 9285 9294 9296 9301 9330 9346 9355 9382 9386 9387 9396 9414 9427 9483 9522 9535 9571 9598 9633 9634 9639 9648 9657 9684 9708 9717 9735 9742 9760 9778 9840 9843 9849 9861 9880 9895 9924 9942 9968 9975 9985</lang>

Java

Works with: Java version 7

<lang java>import java.util.*;

public class SmithNumbers {

   public static void main(String[] args) {
       for (int n = 1; n < 10_000; n++) {
           List<Integer> factors = primeFactors(n);
           if (factors.size() > 1) {
               int sum = sumDigits(n);
               for (int f : factors)
                   sum -= sumDigits(f);
               if (sum == 0)
                   System.out.println(n);
           }
       }
   }
   static List<Integer> primeFactors(int n) {
       List<Integer> result = new ArrayList<>();
       for (int i = 2; n % i == 0; n /= i)
           result.add(i);
       for (int i = 3; i * i <= n; i += 2) {
           while (n % i == 0) {
               result.add(i);
               n /= i;
           }
       }
       if (n != 1)
           result.add(n);
       return result;
   }
   static int sumDigits(int n) {
       int sum = 0;
       while (n > 0) {
           sum += (n % 10);
           n /= 10;
       }
       return sum;
   }

}</lang>

4
22
27
58
85
94
121
...
9924
9942
9968
9975
9985

Lua

Slightly long-winded prime factor function but it's a bit faster than the 'easy' way. <lang Lua>-- Returns a boolean indicating whether n is prime function isPrime (n)

   if n < 2 then return false end
   if n < 4 then return true end
   if n % 2 == 0 then return false end
   for d = 3, math.sqrt(n), 2 do
       if n % d == 0 then return false end
   end
   return true

end

-- Returns a table of the prime factors of n function primeFactors (n)

   local pfacs, divisor = {}, 1
   if n < 1 then return pfacs end
   while not isPrime(n) do
       while not isPrime(divisor) do divisor = divisor + 1 end
       while n % divisor == 0 do
           n = n / divisor
           table.insert(pfacs, divisor)
       end
       divisor = divisor + 1
       if n == 1 then return pfacs end
   end
   table.insert(pfacs, n)
   return pfacs

end

-- Returns the sum of the digits of n function sumDigits (n)

   local sum, nStr = 0, tostring(n)
   for digit = 1, nStr:len() do
       sum = sum + tonumber(nStr:sub(digit, digit))
   end
   return sum

end

-- Returns a boolean indicating whether n is a Smith number function isSmith (n)

   if isPrime(n) then return false end
   local sumFacs = 0
   for _, v in ipairs(primeFactors(n)) do
       sumFacs = sumFacs + sumDigits(v)
   end
   return sumFacs == sumDigits(n)

end

-- Main procedure for n = 1, 10000 do

   if isSmith(n) then io.write(n .. "\t") end

end</lang> Seems silly to paste in all 376 numbers but rest assured the output agrees with https://oeis.org/A006753

Perl 6

<lang perl6>constant @primes = 2, |(3, 5, 7 ... *).grep: *.is-prime;

multi factors ( 1 ) { 1 } multi factors ( Int $remainder is copy ) {

 gather for @primes -> $factor {

   # if remainder < factor², we're done
   if $factor * $factor > $remainder {
     take $remainder if $remainder > 1;
     last;
   }

   # How many times can we divide by this prime?
   while $remainder %% $factor {
       take $factor;
       last if ($remainder div= $factor) === 1;
   }
 }

}

  1. Code above here is verbatim from RC:Count_in_factors#Perl6

sub is_smith_number ( Int $n ) {

 (!$n.is-prime) and ( [+] $n.comb ) == ( [+] factors($n).join.comb );

}

my @s = grep &is_smith_number, 2 ..^ 10_000; say "{@s.elems} Smith numbers below 10_000"; say 'First 10: ', @s[ ^10 ]; say 'Last 10: ', @s[ *-10 .. * ];</lang>

Output:
376 Smith numbers below 10_000
First 10: (4 22 27 58 85 94 121 166 202 265)
Last  10: (9843 9849 9861 9880 9895 9924 9942 9968 9975 9985)

REXX

unoptimized

<lang rexx>/*REXX program finds (and maybe displays) Smith (or joke) numbers up to a given N.*/ parse arg N . /*obtain optional argument from the CL.*/ if N== | N=="," then N=10000 /*Not specified? Then use the default.*/ tell= (N>0); N=abs(N) /*use the │N│ for computing (below).*/ w=length(N) /*W: used for aligning Smith numbers. */

  1. =0 /*#: Smith numbers found (so far). */

@=; do j=4 to N; /*process almost all numbers up to N. */

    if sumDigs(j)\==sumfactr(j)  then iterate   /*Not a Smith number?   Then ignore it.*/
    #=#+1                                       /*bump the Smith number counter.       */
    if \tell  then iterate                      /*Not showing the numbers? Keep looking*/
    @=@ right(j, w);       if length(@)>1255  then do;    say substr(@, 2);    @=;    end
    end   /*j*/                                 /* [↑]  if N>0,  then display Smith #s.*/

if @\== then say substr(@, 2) /*if any residual Smith #s, display 'em*/ say /* [↓] display the number of Smith #s.*/ say # ' Smith numbers found up to ' N " (inclusive)." /*display # of Smith #s found.*/ exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ sumDigs: parse arg x 1 s 2; do d=2 for length(x)-1; s=s+substr(x,d,1); end; return s /*──────────────────────────────────────────────────────────────────────────────────────*/ sumFactr: procedure; parse arg z 1 n; $=0; f=0 /*initialize sum & #factors*/

 do  while z// 2==0;  $=$+2; f=f+1; z=z% 2; end             /*maybe add factor of   2  */
 do  while z// 3==0;  $=$+3; f=f+1; z=z% 3; end             /*  "    "     "    "   3  */
 do j=5  by 2  while j<=z & j*j<=n; if j//3==0 then iterate /*skip factors that are ÷ 3*/
     do  while z//j==0; f=f+1; $=$+sumDigs(j); z=z%j; end   /*maybe reduce  Z  by   J. */
 end   /*j*/                                                /* [↓]  Z  is what's left. */

if z\==1 then do; f=f+1; $=$+sumDigs(z); end /*if residual, than add Z. */ if f<2 then return 0 /*Prime? It's not a Smith#*/

              return $                                      /*Composite? Return sumDigs*/</lang>

output   when using the default number:

    4    22    27    58    85    94   121   166   202   265   274   319   346   355   378   382   391   438   454   483   517
  526   535   562   576   588   627   634   636   645   648   654   663   666   690   706   728   729   762   778   825   852
  861   895   913   915   922   958   985  1086  1111  1165  1219  1255  1282  1284  1376  1449  1507  1581  1626  1633  1642
 1678  1736  1755  1776  1795  1822  1842  1858  1872  1881  1894  1903  1908  1921  1935  1952  1962  1966  2038  2067  2079
 2155  2173  2182  2218  2227  2265  2286  2326  2362  2366  2373  2409  2434  2461  2475  2484  2515  2556  2576  2578  2583
 2605  2614  2679  2688  2722  2745  2751  2785  2839  2888  2902  2911  2934  2944  2958  2964  2965  2970  2974  3046  3091
 3138  3168  3174  3226  3246  3258  3294  3345  3366  3390  3442  3505  3564  3595  3615  3622  3649  3663  3690  3694  3802
 3852  3864  3865  3930  3946  3973  4054  4126  4162  4173  4185  4189  4191  4198  4209  4279  4306  4369  4414  4428  4464
 4472  4557  4592  4594  4702  4743  4765  4788  4794  4832  4855  4880  4918  4954  4959  4960  4974  4981  5062  5071  5088
 5098  5172  5242  5248  5253  5269  5298  5305  5386  5388  5397  5422  5458  5485  5526  5539  5602  5638  5642  5674  5772
 5818  5854  5874  5915  5926  5935  5936  5946  5998  6036  6054  6084  6096  6115  6171  6178  6187  6188  6252  6259  6295
 6315  6344  6385  6439  6457  6502  6531  6567  6583  6585  6603  6684  6693  6702  6718  6760  6816  6835  6855  6880  6934
 6981  7026  7051  7062  7068  7078  7089  7119  7136  7186  7195  7227  7249  7287  7339  7402  7438  7447  7465  7503  7627
 7674  7683  7695  7712  7726  7762  7764  7782  7784  7809  7824  7834  7915  7952  7978  8005  8014  8023  8073  8077  8095
 8149  8154  8158  8185  8196  8253  8257  8277  8307  8347  8372  8412  8421  8466  8518  8545  8568  8628  8653  8680  8736
 8754  8766  8790  8792  8851  8864  8874  8883  8901  8914  9015  9031  9036  9094  9166  9184  9193  9229  9274  9276  9285
 9294  9296  9301  9330  9346  9355  9382  9386  9387  9396  9414  9427  9483  9522  9535  9571  9598  9633  9634  9639  9648
 9657  9684  9708  9717  9735  9742  9760  9778  9840  9843  9849  9861  9880  9895  9924  9942  9968  9975  9985

376  Smith numbers found up to  10000  (inclusive).

optimized

This REXX version uses a faster version of the   sumFactr   function. <lang rexx>/*REXX program finds (and maybe displays) Smith (or joke) numbers up to a given N.*/ parse arg N . /*obtain optional argument from the CL.*/ if N== | N=="," then N=10000 /*Not specified? Then use the default.*/ tell= (N>0); N=abs(N) /*use the │N│ for computing (below).*/

  1. =0; do i=11 to N;  !.i=i*i; end /*i*/ /*!.: array to hold squared numbers. */

w=length(N) /*W: used for aligning Smith numbers. */ @=; do j=4 to N /*process almost all numbers up to N. */

      if sumDigs(j)\==sumfactr(j)  then iterate /*Not a Smith number?   Then ignore it.*/
      #=#+1                                     /*bump the Smith number counter.       */
      if \tell  then iterate                    /*Not showing the numbers? Keep looking*/
      @=@ right(j, w);         if length(@)>75 then do;   say substr(@, 2);    @=;    end
      end   /*j*/                               /* [↑]  if N>0,  then display Smith #s.*/

if @\== then say substr(@, 2) /*if any residual Smith #s, display 'em*/ say /* [↓] display the number of Smith #s.*/ say # ' Smith numbers found up to ' N " (inclusive)." /*display # of Smith #s found.*/ exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ sumDigs: parse arg x 1 s 2; do d=2 for length(x)-1; s=s+substr(x,d,1); end; return s /*──────────────────────────────────────────────────────────────────────────────────────*/ sumFactr: procedure expose !.; parse arg z 1 n; $=0; f=0

 do  while z// 2==0;    $=$+2;   f=f+1;   z=z% 2;   end     /*maybe add factor of   2  */
 do  while z// 3==0;    $=$+3;   f=f+1;   z=z% 3;   end     /*  "    "     "    "   3  */
 do  while z// 5==0;    $=$+5;   f=f+1;   z=z% 5;   end     /*  "    "     "    "   5  */
 do  while z// 7==0;    $=$+7;   f=f+1;   z=z% 7;   end     /*  "    "     "    "   7  */
 do j=11 by 6  while j<=z & !.j<=n                          /*skip factors that are ÷ 3*/
 parse var  j     -1  _;                if _\==5 then,    /*is last dec. digit ¬a 5 ?*/
       do  while z//j==0; f=f+1; $=$+sumDigs(j); z=z%j; end /*maybe reduce  Z  by   J. */
 y=j+2
 if _==3  then iterate
       do  while z//y==0; f=f+1; $=$+sumDigs(y); z=z%y; end /*maybe reduce  Z  by   Y. */
 end   /*j*/                                                /* [↓]  Z  is what's left. */

if z\==1 then do; f=f+1; $=$+sumDigs(z); end /*if residual, than add Z. */ if f<2 then return 0 /*Prime? It's not a Smith#*/

              return $                                      /*Composite? Return sumDigs*/</lang>

output   when using the input of (negative) one million:   -1000000

29928  Smith numbers found up to  1000000  (inclusive).

Ruby

<lang ruby>require "prime"

class Fixnum

 def sum_digits
   to_s.chars.map(&:to_i).inject(:+)
 end
 def smith?
   return false if prime?
   sum_digits == prime_division.map{|pr,n| pr.sum_digits * n}.inject(:+)
 end

end

n = 10_000 res = 1.upto(n).select(&:smith?)

puts "#{res.size} smith numbers below #{n}:

  1. {res.first(5).join(", ")},... #{res.last(5).join(", ")}"</lang>
Output:
376 smith numbers below 10000:
4, 22, 27, 58, 85,... 9924, 9942, 9968, 9975, 9985

zkl

Uses the code (primeFactors) from Prime decomposition#zkl. <lang zkl>fcn smithNumbers(N=0d10_000){ // -->(Smith numbers to N)

  [2..N].filter(fcn(n){ 
     (pfs:=primeFactors(n)).len()>1 and
     n.split().sum(0)==primeFactors(n).apply("split").flatten().sum(0) 
  })

}</lang> <lang zkl>sns:=smithNumbers(); sns.toString(*).println(" ",sns.len()," numbers");</lang>

Output:
L(4,22,27,58,85,94,121,166,202,265,274,319,346,355,378,382,391, ...
3091,3138,3168,3174,3226,3246,3258,3294,3345,3366,3390,3442,3505, ...
9942,9968,9975,9985) 376 numbers