Oxford Oberon : factoring

If you can determine prime numbers, you alse need to be able to determine non prime numbers. Every non-prime number can be written as a product of a series of primes, like in:

num = p1^n * p2^m * ... * px^z

For example, 39 = 3 x 13 and 819819 = 3^2 x 7^2 x 11 x 13^2. For small numbers (below 200), we, humans, can factor from head. No need for pencil and paper. But for bigger numbers, like 819819 things get a bt more difficult. Even 511 is already hard to factor in your head: 511 = 7 x 73.

So I wrote a factoring program in obc. It relies heavily on the primes program in that it uses a list of primes to walk through. The list is simply generated by

primes 70000 >p70k

Why 70k? Easy. The biggest number to fit in an INTEGER is 2^31. And so it doesn't make sense to search further than 2^15.5 which I roughly estimated to be 70.000 but that was false of course. 2^15.5 = roughly 46000 so I made the list too big. But size doesn't matter. Or does it? Depends on the things we're measuring I guess.

The algorithm is easy:

Pick the next prime
Is this prime bigger than the root of number? If so exit
As long as the number is divisible by the prime, do it
If this prime was a divisor, do some printing

The primes are read in from a file (see above). And that's just about it. This is a very small and simple program. Version one worked 'out of the box', apart from some typos. This made me very conspicious. NO program works at the first trial. But this one did. Scary.

This version only looks for one divisor per prime. It produces, for example:

oxygen:~/Oberon/prime$ ./factor01 819819
Primes read :    6937
3 x 273273
7 x 117117
11 x 74529
13 x 63063

Whereas the number 819819 is in fact 3^2 x 7^2 x 11 x 13^2. 273273 is divisible by 3 but that is not detected. Still, factor01 does a good job for an initial version.

The source code

Below is the source code for factor01.mod for the Spivey Obc compiler:

MODULE factor01;

(*	Find factors in a number	*)

IMPORT	Args, Conv, Files, In, Out;

VAR	primes		: ARRAY 10000 OF INTEGER;
	noParams	: BOOLEAN;
	factor		: INTEGER;


PROCEDURE ReadParams;

VAR	count		: INTEGER;
	option		: ARRAY  16  OF CHAR;

BEGIN
  count := Args.argc;
  IF  count = 1  THEN
    noParams := TRUE
  ELSE
    noParams := FALSE;
    Args.GetArg (1, option);
    factor := Conv.IntVal (option)
  END
END ReadParams;


PROCEDURE Init;

VAR	inF		: Files.File;
	i, j, n		: INTEGER;
	ch		: CHAR;
	num		: ARRAY 12 OF CHAR;

BEGIN
  inF := Files.Open ("p70k", "r");
  IF  inF = NIL  THEN  HALT (1)  END;
  j := 0;
  REPEAT
    REPEAT  Files.ReadChar (inF, ch)  UNTIL  ch > ' ';
    i := 0;
    WHILE  ch > ' '  DO
      num [i] := ch;
      Files.ReadChar (inF, ch);
      INC (i)
    END;
     n := Conv.IntVal (num);
     primes [j] := n;
     INC (j)
  UNTIL n = 0;
  Files.Close (inF);
  Out.String ("Primes read : ");	Out.Int (j, 7);		Out.Ln
END Init;


PROCEDURE Split (num : INTEGER);

VAR	i, j, pri, quo, rem	: INTEGER;

BEGIN
  i := 0;
  LOOP
    pri := primes [i];
    IF  pri * pri > num  THEN  EXIT  END;
    quo := num DIV pri;
    rem := num MOD pri;
    IF  rem = 0  THEN
      Out.Int (pri, 1);
      Out.String (' x ');
      Out.Int (quo, 1);
      Out.Ln
    END;
    INC (i);
  END
END Split;


BEGIN
  Init;
  ReadParams;
  IF  noParams  THEN
    LOOP
      In.Int (factor);
      IF  factor = 0  THEN  EXIT  END;
      Split (factor)
    END
  ELSE
    Split (factor)
  END
END factor01.

It's all rather straightforward. In 'Init' the primes file is entered into the array, then the command line parameters (if any) are processed and the number is split up in it's factors. As mentioned, this version ran almost without changes. One more example, this time without command tail arguments:

oxygen:~/Oberon/prime$ ./factor01
Primes read :    6937
819
3 x 273
7 x 117
13 x 63
819819
3 x 273273
7 x 117117
11 x 74529
13 x 63063
117117
3 x 39039
7 x 16731
11 x 10647
13 x 9009
0
oxygen:~/Oberon/prime$

Factor01 did a good job in isolating the major factors. But I want to have all factors. And I need them on one line of output. So the 'Split' routine was slightly changed into:

PROCEDURE Split (num : INTEGER);

VAR	i, j, pri	: INTEGER;

BEGIN
  i := 0;				Set up primes array index
  LOOP					Infinite loop
    pri := primes [i];			    Buffer current prime
    j := 0;				    Initialize 'powers' counter
    IF  pri * pri > num  THEN 		    STOP condition?
      Out.Int (num, 1);				IF so, report to CRT
      EXIT  					And leave infinite loop
    END;				    .
    WHILE  num MOD pri = 0  DO		    As long as 'numm' can be divided by 'pri' DO
      num := num DIV pri;			Modify 'num'
      INC (j)					Increment 'powers' counter
    END;				    .
    IF  j > 0  THEN			    Did we have ANY divisors for this prime?
      Out.Int (pri, 1);				If so, print the prime
      IF  j > 1  THEN				Was this prime a factor more than once?
        Out.Char ('^');				    If so, print the caret power token
        Out.Int (j, 1)				    And the exponent
      END;					.
      IF  num # 1  THEN				If num > 1, we're not done yet so
      Out.String (" x ")  			    print the multiplication sign
      ELSE  					ELSE
        EXIT  					    get the hell out of the infinite loop
      END					.
    END;				    .
    INC (i)				    Next prime number
  END;					Continue Infinite loop
  Out.Ln				We're out of the LOOP, send an End of Line...
END Split;

I didn't put the comments in Oberon comment tokens for lazyness on my side. As you can see, num is getting smaller after each detected divisor. So numbers get smaller very fast. Below is a small example of the output of factor02:

oxygen:~/Oberon/prime$ ./factor02 65536
Primes read :    6937
2^16
oxygen:~/Oberon/prime$

Indeed, 65536 = 2^16.

Factor02 : the source

Below is the full source code for factor02.mod.

MODULE factor02;

(*	Find factors in a number	*)

IMPORT	Args, Conv, Files, In, Out;

VAR	primes		: ARRAY 10000 OF INTEGER;
	noParams	: BOOLEAN;
	factor		: INTEGER;


PROCEDURE ReadParams;

VAR	count		: INTEGER;
	option		: ARRAY  16  OF CHAR;

BEGIN
  count := Args.argc;
  IF  count = 1  THEN
    noParams := TRUE
  ELSE
    noParams := FALSE;
    Args.GetArg (1, option);
    factor := Conv.IntVal (option)
  END
END ReadParams;


PROCEDURE Init;

VAR	inF		: Files.File;
	i, j, n		: INTEGER;
	ch		: CHAR;
	num		: ARRAY 12 OF CHAR;

BEGIN
  inF := Files.Open ("p70k", "r");
  IF  inF = NIL  THEN  HALT (1)  END;
  j := 0;
  REPEAT
    REPEAT  Files.ReadChar (inF, ch)  UNTIL  ch > ' ';
    i := 0;
    WHILE  ch > ' '  DO
      num [i] := ch;
      Files.ReadChar (inF, ch);
      INC (i)
    END;
     n := Conv.IntVal (num);
     primes [j] := n;
     INC (j)
  UNTIL n = 0;
  Files.Close (inF);
  Out.String ("Primes read : ");	Out.Int (j, 7);		Out.Ln
END Init;


PROCEDURE Split (num : INTEGER);

VAR	i, j, pri	: INTEGER;

BEGIN
  i := 0;
  LOOP
    pri := primes [i];
    j := 0;
    IF  pri * pri > num  THEN  
      Out.Int (num, 1);
      EXIT  
    END;
    WHILE  num MOD pri = 0  DO
      num := num DIV pri;
      INC (j)
    END;
    IF  j > 0  THEN
      Out.Int (pri, 1);
      IF  j > 1  THEN
	Out.Char ('^');
        Out.Int (j, 1)
      END;
      IF  num # 1  THEN  Out.String (" x ")  ELSE  EXIT  END
    END;
    INC (i)
  END;
  Out.Ln
END Split;


BEGIN
  Init;
  ReadParams;
  IF  noParams  THEN
    LOOP
      In.Int (factor);
      IF  factor = 0  THEN  EXIT  END;
      Split (factor)
    END
  ELSE
    Split (factor)
  END
END factor02.

I have already explained how the Split routine was changed. So let's just concentrate on some results:

oxygen:~/Oberon/prime$ ./factor02 819819  
Primes read :    6937
3^2 x 7^2 x 11 x 13^2
oxygen:~/Oberon/prime$ ./factoe02 81981987
Primes read :    6937
3 x 61 x 283 x 1583
oxygen:~/Oberon/prime$ ./factor02 511     
Primes read :    6937
7 x 73

It's amazing to see that a big number like 819819 (almost one million) has a highest factor of '13'. This means that the 'higher primes' are only 'used' when trying to factor a very large prime number. Initially, for certain numbers, the output was similar to

oxygen:~/Oberon/prime$ ./fact2 819819   
Primes read :    6937
3^2 x 7^2 x 11 x 13^2 x 1

where the line ended in the silly looking 'x 1'. It was tackled by the last 'IF' clause of 'Split':

IF  num # 1  THEN  Out.String (" x ")  ELSE  EXIT  END

When num = 1, we just leave.

Instead of typing or using copy/paste you can easily download the full sources. Take care however when you download the files for something else than Linux or Unix: you may have to fix the linefeeds.

Page created 10 Feb 2017,