[futurebasic] Re: [FB] Prime Numbers?

Message: < previous - next > : Reply : Subscribe : Cleanse
Home   : December 1998 : Group Archive : Group : All Groups

From: Robert Purves <robert.purves@...>
Date: Tue, 22 Dec 1998 14:46:04 +1300
At 8:39 AM +1300 21/12/97 [sic] , tedd@... wrote:

>If we're going to have a demonstration of talent, why don't we see what's
>the highest prime we can resolve and in shortest amount of time? At least
>we would be doing something that might be noteworthy. If nothing else, the
>endeavour would at least expand out instead of collapsing inward.

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