/*PD1 REXX pgm does prime decomposition of a range of positive integers (with a prime count)*/
Call Time 'R'
numeric digits 1000                              /*handle thousand digits for the powers*/
parse arg  bot  top  step   base  add            /*get optional arguments from the C.L. */
if  bot==''   then do;  bot=1;  top=100;  end    /*no  BOT given?  Then use the default.*/
if  top==''   then              top=bot          /* "  TOP?  "       "   "   "     "    */
if step==''   then step=  1                      /* " STEP?  "       "   "   "     "    */
if add ==''   then  add= -1                      /* "  ADD?  "       "   "   "     "    */
tell= top>0;       top=abs(top)                  /*if TOP is negative, suppress displays*/
w=length(top)                                    /*get maximum width for aligned display*/
if base\==''  then w=length(base**top)           /*will be testing powers of two later? */
commat.=left('', 7);   commat.0="{unity}";   commat.1='[prime]' /*some literals:  pad;  prime (or not).*/
numeric digits max(9, w+1)                       /*maybe increase the digits precision. */
hash=0                                              /*hash:    is the number of primes found. */
        do n=bot  to top  by step                /*process a single number  or  a range.*/
        ?=n;  if base\==''  then ?=base**n + add /*should we perform a "Mercenne" test? */
        pf=factr(?);      f=words(pf)            /*get prime factors; number of factors.*/
        if f==1  then hash=hash+1                      /*Is N prime?  Then bump prime counter.*/
        if tell  then say right(?,w)   right('('f")",9)   'prime factors: '     commat.f     pf
        end   /*n*/
say
ps= 'primes';    if p==1  then ps= "prime"       /*setup for proper English in sentence.*/
say right(hash, w+9+1)       ps       'found.'      /*display the number of primes found.  */
Say 'PD1 took' time('E') 'seconds'
exit                                             /*stick a fork in it,  we're all done. */
/*--------------------------------------------------------------------------------------*/
factr: procedure;  parse arg x 1 d,dollar             /*set X, D  to argument 1;  dollar  to null.*/
if x==1  then return ''                          /*handle the special case of   X = 1.  */
       do  while x//2==0;  dollar=dollar 2;  x=x%2;  end   /*append all the  2  factors of new  X.*/
       do  while x//3==0;  dollar=dollar 3;  x=x%3;  end   /*   "    "   "   3     "     "  "   " */
       do  while x//5==0;  dollar=dollar 5;  x=x%5;  end   /*   "    "   "   5     "     "  "   " */
       do  while x//7==0;  dollar=dollar 7;  x=x%7;  end   /*   "    "   "   7     "     "  "   " */
                                                 /*                                  ___*/
q=1;   do  while q<=x;  q=q*4;  end              /*these two lines compute integer  v X */
r=0;   do  while q>1;   q=q%4;  _=d-r-q;  r=r%2;   if _>=0  then do; d=_; r=r+q; end;  end

       do j=11  by 6  to r                       /*insure that  J  isn't divisible by 3.*/
       parse var j  ''  -1  _                    /*obtain the last decimal digit of  J. */
       if _\==5  then  do  while x//j==0;  dollar=dollar j;  x=x%j;  end     /*maybe reduce by J. */
       if _ ==3  then iterate                    /*Is next  Y  is divisible by 5?  Skip.*/
       y=j+2;          do  while x//y==0;  dollar=dollar y;  x=x%y;  end     /*maybe reduce by J. */
       end   /*j*/
                                                 /* [?]  The dollar list has a leading blank.*/
if x==1  then return dollar                           /*Is residual=unity? Then don't append.*/
              return dollar x                         /*return   dollar   with appended residual. */