*process source attributes xref nest or(!);
 /*--------------------------------------------------------------------
 * 22.02.2014 Walter Pachl using the is_prime code from
 *                         PL/I 'prime decomposition'
 * 23.02.  WP start test for second prime with 2 or first prime found
 *-------------------------------------------------------------------*/
 spb: Proc options(main);
 Dcl a(10) Bin Fixed(31)
          Init(900660121,2,4,1679,1234567,32768,99,9876543,100,5040);
 Dcl (x,n,nf,i,j) Bin Fixed(31) Init(0);
 Dcl f(3) Bin Fixed(31);
 Dcl txt Char(30) Var;
 Dcl bit Bit(1);
 Do i=1 To hbound(a);
   bit=is_semiprime(a(i));
   Select(nf);
     When(0,1) txt=' is prime';
     When(2)   txt=' is     semiprime '!!factors(a(i));
     Otherwise txt=' is NOT semiprime '!!factors(a(i));
     End;
   Put Edit(a(i),bit,txt)(Skip,f(10),x(1),b(1),a);
   End;

 is_semiprime: Proc(x) Returns(bit(1));
 /*--------------------------------------------------------------------
 * Returns '1'b if x is semiprime, '0'b otherwise
 * in addition
 * it sets f(1) and f(2) to the first (or only) prime factor(s)
 *-------------------------------------------------------------------*/
   Dcl x Bin Fixed(31);
   nf=0;
   f=0;
   x=a(i);
   n=x;
   f(1)=2;
 loop:
   Do While(nf<=2 & n>1);
     If is_prime(n) Then Do;
       Call mem(n);
       Leave loop;
       End;
     Else Do;
 loop2:
       Do j=f(1) By 1 While(j*j<=n);
         If is_prime(j)&mod(n,j)=0 Then Do;
           Call mem(j);
           n=n/j;
           Leave loop2;
           End;
         End;
       End;
     End;
   Return(nf=2);
 End;

 is_prime: Proc(n) Returns(bit(1));
 Dcl n Bin Fixed(31);
 Dcl i Bin Fixed(31);
   If n < 2 Then Return('0'b);
   If n = 2 Then Return('1'b);
   If mod(n,2)=0 Then Return('0'b);
   Do i = 3 by 2 While(i*i<=n);
     If mod(n,i)=0 Then Return('0'b);
     End;
   Return('1'b);
 End is_prime;

 mem: Proc(x);
 Dcl x Bin Fixed(31);
   nf+=1;
   f(nf)=x;
 End;

 factors: Proc(x) Returns(Char(150) Var);
 Dcl x Bin Fixed(31);
 Dcl (res,net) Char(150) Var Init('');
 Dcl (i,f3) Bin Fixed(31);
 res=f(1)!!'*'!!f(2);
 f3=x/(f(1)*f(2));
 If f3>1 Then
   res=res!!'*'!!f3;
 Do i=1 To length(res);
   If substr(res,i,1)>' ' Then
     net=net!!substr(res,i,1);
   End;
 Return(net);
 End;

 End spb;