FB II Compiler

PG PRO

Debugging

Memory

System

Mathematics

Resources

Disk I/O

Windows

Controls

Menus

Mouse

Keyboard

Text

Fonts

Drawing

Sound

Clipboard

Printing

Communication

ASM

Made with FB

MATHEMATICS

Find prime numbers


See how this does for you:

long if myNumber% / 2 <> int(myNumber% / 2)
  isPrime% = _true 'or that _ztrue thing if you like
  for loop% = 3 to myNumber% / 2 step 2 'Skip the even numbers
    long if myNumber% / loop% = int (myNumber% / loop%) 'Is it a factor?
       isPrime% = _false 'not a prime!
       loop% = myNumber% / 2 'Will this dump us out of the loop ok? (optional)
    end if
  next loop%
xelse
  isPrime% = _false 'it was an even number
  if myNumber% = 1 then isPrime% = _true 'almost forgot this one! 1 is prime, right?
end if

Paul Bruneau


At Last, something I can answer !

LOCAL FN IsItAPrime% (Num&)
  IsPrime% = _True
  LONG IF (Num&+1 AND 1)
    IsPrime% = _False
  XELSE
    HalfWay& = Num& / 2
    div& = 2
    DO
     INC(div&)
     IF Num& MOD div& = 0 THEN IsPrime% = _False
    UNTIL IsPrime% = _False OR div& > HalfWay&
  END IF
END FN = IsPrime%

Ian Mann


Try this. I just wrote it up. It works for me.

Tedd

COMPILE 0,_dimmedVarsOnly
END GLOBALS

'--------------------------

LOCAL FN buildWind
  WINDOW 1,"test",(0,0)-(500,500),_doczoom
END FN

'--------------------------
' is a number prime?

LOCAL FN isPrime(num)
  DIM i
  DIM test

  test = _true
  PRINT "Number is = ";num
  FOR i = 2 TO num -1
    LONG IF num MOD i = 0
      PRINT "Divisors = ";i
      test = _false
    END IF
  NEXT

END FN = test

'--------------------------
' find the primes within a number

LOCAL FN findPrime(num)
  DIM i, k
  DIM test
  DIM a

  PRINT "Number is = ";num
  FOR k = 1 TO num
    test = _true

    FOR i = 2 TO k-1
      LONG IF k MOD i = 0
        test = _false
      END IF
    NEXT

    LONG IF test = _true
      PRINT "Primes are = ";i
    END IF

  NEXT

END FN

'--------------------------

"MAIN"
DIM i
DIM a$
DIM num
DIM test

FN buildWind

num = 242 'check with a non-prime number
test = FN isPrime (num)
LONG IF test = _true
  PRINT num;" is prime"
XELSE
  PRINT num;" is not prime"
END IF
PRINT

num = 11 'check with a prime number
test = FN isPrime (num)
LONG IF test = _true
  PRINT num;" is prime"
XELSE
  PRINT num;" is not prime"
END IF
PRINT

FN findPrime(100) 'find all primes wihtin this number

PRINT
INPUT"Enter anything to end";a$

END


Why don't you, once you have the prime numbers, write them out as a binary string, with the bit set if the corresponding number is a prime -- i.e., the first bit is 0, because "0" is not a prime, the second is 0 because "1" is not a prime, the third and fourth are set to 1, because "2" and "3" are prime, etc.
This takes only 832 bits, is 104 bytes, plus an FN BITTEST, surely less than a complete function, and very fast...

Hans van Maanen


Cool suggestion. I, however, took it as a challenge. I couldn't get my code down to 104 bytes, but at 146 byte it's not _that_ far off. This works for numbers from 2 to 832:

LOCAL FN IsItAPrime (Num&)
  FOR factor& = Num&/2 TO 2 STEP -1
    IF Num& MOD factor& = 0 THEN factor& = - factor&
  NEXT
END FN = factor& > -3

Of course your suggestion blows this away speed-wise. :-)

Jay


I noticed most people's algorithms were checking divisors all the way up to num&/2. It's only necessary to check up to SQR(num&):

LOCAL FN IsItPrime(num&)
  SELECT CASE
    CASE num& = 2 OR num& = 3
      prime = _zTrue
    CASE num& MOD 2 = 0
      prime = _false
    CASE ELSE
      prime = _zTrue '(initially)
      root = SQR(num&)
      FOR i = 3 TO root STEP 2
        LONG IF num& MOD i = 0
          prime = _false
          i = root 'force early exit from loop
        END IF
      NEXT
  END SELECT
END FN = prime

In fact, it's only necessary to check _prime_ potential divisors up to SQR(num&), which further shortens the list. So to check whether 832 is prime, you only need to see if it's divisible by any of these numbers: 2,3,5,7,11,13,17,19,23. My short function above doesn't have that optimization, but you could potentially use an array of prime numbers to use as test divisors.

Rick


There's no contest. Here's a short program comparing my 5-line FN (which probably is not the fastest) with Hans's BITTST suggestion. Evaluating 10,000 random numbers in the range 1-832 with my code takes my machine at least 155 ticks. Looking up the same 10,000 numbers in Hans's list takes less than 3 ticks!. (BTW, be prepared for a time increase of 10-20x if you run this in the debugger.) :-)

Jay

GLOBALS
_trials = 10000 'Set number of trials here
DIM gPrimeBits.104
DIM randomList(_trials)
DIM gTime&
DIM r&,isPrime,d$
END GLOBALS

LOCAL FN IsItAPrime (Num&) 'Jay's FN
  'Only works for integers > 1
  FOR factor& = Num&/2 TO 2 STEP -1
    IF Num& MOD factor& = 0 THEN factor& = - factor&
  NEXT
END FN = factor& > -3

LOCAL FN makeLists
  FOR r = 2 TO 832 'Set up Hans' list using Jay's FN
    IF FN IsItAPrime(r) THEN CALL BITSET(#@gPrimeBits,r)
  NEXT
  FOR r = 1 TO _trials 'Choose #'s to test
    randomList(r) = RND(831)+1
  NEXT

END FN

'------------Main---------------
FN makeLists
CLS
PRINT "For";_trials;"trials..."

gTime& = FN TICKCOUNT+1
WHILE gTime& < FN TICKCOUNT
WEND
FOR r& = 1 TO _trials
  isPrime = FN IsItAPrime(randomList(r&))
NEXT
PRINT "The code took"; FN TICKCOUNT-gTime&;"ticks."

gTime& = FN TICKCOUNT+1
WHILE gTime& < FN TICKCOUNT:WEND
FOR r& = 1 TO _trials
  isPrime = -FN BITTST(#@gPrimeBits,randomList(r&))
NEXT
PRINT "Bit test took";FN TICKCOUNT-gTime&;"ticks."
PRINT

CURSOR _arrowcursor
COLOR _zred
INPUT "Press Return to end"; d$
PRINT


I claim (at least temporarily) the prize. The program below determines primeness and factors of any number from 1 to 4294967295 which is the maximum value of an unsigned LONG. Rick's routine, posted last week, works up to 2147483647. In a head-to-head test (on all numbers from 1 to 1 million), my routine is 11 times faster, partly because I replaced SQR(num&) by a very fast integer square root. Along the way it turned out that the obvious replacement for SQR, USR _sqRoot, has some defects, which are remedied in FN SquareRoot&.

BTW, Tedd, isn't it time you set your date to 1998? Are you are hoping to postpone Y2K to 2001? ;-)

Robert Purves

LOCAL FN SquareRoot&(a&)
' Replacement for USR _sqRoot (a&) which has three faults:-
' 1. It crashes for a&=_Maxlong/2
' 2. It gives wrong answer for a&>_Maxlong/2
' 3. it gves "high-by-one" answers in many cases, e.g. USR _sqRoot(15)=4
' This routine works for all values of a& unsigned (0-4294967295)
' and is slightly faster than USR _sqRoot. It DOES NOT WORK on 68000 machines.
DIM root&
` move.l ^a&,d3
` move.l d3,d0
` beq.s SqDone ; skip if 0
` move.l d3,d1
`ApprLoop ; init root by halving the number of digits in a&
` lsr.l #1,d0
` lsr.l #2,d1
` bne.s ApprLoop
` addq.l #1,d0 ; force non-zero
` move.w #4,d7 ; max loop count
`SqLoop ; iterate root& = (a&/root&+root& )/2
` move.l d0,d5 ; copy of root
` clr.l d1 ; hi.l of a&=0
` move.l d3,d2 ; lo.l of a&
` dc.w $4C40 ; DIVU.L D0,D1:D2
` dc.w $2601 ; d2.l= d1(hi) d2(lo) / d0.l d1.l = remainder
` add.l d2,d0
` lsr.l #1,d0
` cmp.l d0,d5 ; equals old? If so, we are done
` dbeq d7, SqLoop ; branch back if <> old and loop count still >0
' IF root&*root& > a& THEN dec(root&)
` move.l d0,d2
` dc.w $4c00 ; MULU.L D0,D1:D2
` dc.w $2401 ;d1 (hi) and d2 (lo) = d0.l * d2.l
` cmp.l d2,d3
` bge.s SqDone
` subq.l #1,d0
`SqDone
` move.l d0,^root&
END FN=root&

LOCAL FN IsItPrimeRDP(num&,fact1Ptr&,fact2Ptr&)
' Brute force test for primeness
' Works for unsigned num& 0<=num&<= 2^32-1 (4294967295)
' DOES NOT WORK on 68000 machines
DIM prime
SELECT CASE
  CASE num& = 2 OR num& = 3
   prime = _zTrue
   POKE LONG fact1Ptr&,1: POKE LONG fact2Ptr&,num&
  CASE (num& AND 1) = 0' even
   prime = _false
   POKE LONG fact1Ptr&,2: POKE LONG fact2Ptr&,num&>>1
  CASE ELSE
   REGISTER(d3)=FN SquareRoot&(num&)
   ` move.l ^num&,d4
   ` move.w #3,d0 ; i=3 divisor
   `loop
   ` clr.l d1 ; hi of num&=0
   ` move.l d4,d2 ; lo of num&
   ` dc.w $4C40 ; DIVU.L D0,D1:D2
   ` dc.w $2601 ; d2.l= d1(hi) d2(lo) / d0.l d1.l = remainder
   ` tst.l d1
   ` beq.s notPrime ; branch if remainder=0
   ` addq.w #2,d0
   ` cmp.w d3,d0
   ` bls.s loop
   ` move.w #-1,d1 ; flag _ztrue
   ` move.l #1, d0
   ` move.l ^num&,d2
   `notPrime
   ` movea.l ^fact1Ptr&,a0
   ` ext.l d0
   ` move.l d0,(a0)
   ` movea.l ^fact2Ptr&,a0
   ` move.l d2,(a0)
   ` move.w d1,^prime
   `done
END SELECT
END FN = prime

LOCAL FN StringToLongUnsigned&(s$)
' Convert number up to 4294967295
DIM j,a&, 1 char$
a&=0
FOR j=1 TO LEN(s$)
  char$=MID$(s$,j,1)
  LONG IF char$>="0" AND char$<="9"'ignore non-numeric
   a&=a&*10 + VAL(char$)
  END IF
NEXT j
END FN=a&

WINDOW 1
DEFSTR LONG
DIM s$,a&,f1&,f2&
DO
INPUT "Enter number (return to end) :"; s$
IF s$="" THEN END
a&=FN StringToLongUnsigned&(s$)
PRINT "Number is:" VAL(UNS$(a&))
LONG IF FN IsItPrimeRDP(a&,@f1&,@f2&)
  PRINT "P";
XELSE
  PRINT "Not p";
END IF
PRINT "rime. Factors: " VAL(UNS$(f1&))"and"VAL(UNS$(f2&))
UNTIL 0


Rick's and Hans' suggestions are actually equivalent, that is they lead to the same sqeuence of divisors. Nice call, guys.
And by extension we can remove multiples of 5 from the list of divisors as in the following program. It's in plain vanilla FB for comprehensibility.
Even though the previous version contained screeds of assembly code, this one is actually faster. The down-side of no ASM is that it can't factor numbers that are >= 1073741823

COMPILE _dimmedVarsOnly
DIM gIncr(7)
END GLOBALS

LOCAL FN InitIncrements
gIncr(0)=4: gIncr(1)=2: gIncr(2)=4: gIncr(3)=2
gIncr(4)=4: gIncr(5)=6: gIncr(6)=2:gIncr(7)=6
END FN

LOCAL FN RDPIsItPrime2(num&,fact1Ptr&,fact2Ptr&)
DIM prime,root&,divisor&,index
prime = _false
_maxLongBy2=1073741823
SELECT CASE
  CASE num&=0
   & fact1Ptr&,0: & fact2Ptr&,0
  CASE num&=1
   & fact1Ptr&,1: & fact2Ptr&,1
  CASE num&>=_maxLongBy2 OR num&<0
   PRINT "Too big for USR _sqRoot "
   & fact1Ptr&,0: & fact2Ptr&,0
  CASE num& = 2 OR num& = 3 OR num&=5
   prime = _zTrue
   & fact1Ptr&,1: & fact2Ptr&,num&
  CASE (num& AND 1) = 0' even
   & fact1Ptr&,2: & fact2Ptr&,num&>>1
  CASE num& MOD 3 = 0' divisible by 3
   & fact1Ptr&,3: & fact2Ptr&,num&/3
  CASE num& MOD 5 = 0' divisible by 5
   & fact1Ptr&,5: & fact2Ptr&,num&/5
  CASE ELSE
   prime = _zTrue 'initially
   & fact1Ptr&,1: & fact2Ptr&,num&
   ' Start with divisor 7 and add {4,2,4,2,4,6,2,6} in sequence,
   ' thus skipping divisors that are multiples of 2,3 or 5
   root& =USR _sqRoot (num&)
   divisor&=7
   index=0
   WHILE divisor&<=root&
    LONG IF num& MOD divisor&= 0
      prime = _false
      & fact1Ptr&,divisor&: & fact2Ptr&,num&/divisor&
      divisor=root&'force early exit
    END IF
    divisor&=divisor&+gIncr(index)
    index=(index+1) AND 7
   WEND
END SELECT
END FN = prime

WINDOW 1
DEFSTR LONG
DIM s$,a&,f1&,f2&
FN InitIncrements
DO
  INPUT "Enter number (return to end) :"; s$
  IF s$="" THEN END
  a&=VAL(s$)
  PRINT "Number is:" VAL(UNS$(a&))
  LONG IF FN RDPIsItPrime2(a&,@f1&,@f2&)
    PRINT "P";
  XELSE
    PRINT "Not p";
  END IF
  PRINT "rime. Factors: " VAL(UNS$(f1&))"and"VAL(UNS$(f2&))
UNTIL 0

Robert