LEX 'FACTORLX' prime factors, HP71 TITLE FACTORLEX 1A, T. Tarvainen * * Tapani Tarvainen 87/01/31 * * Function FACTOR(N[,F]) * returns the smallest prime factor of integer N * Optional second argument F specifies * starting point for the search, when it is * known that there are no factors below it * (if there are, results are not reliable) * * FACTOR(0)=0 and FACTOR(1)=1 * * Execution time: up to 7.5 minutes (for N=999999999989) * FACTOR can only be interrupted by INIT:1 * ***************************************************** * * Errors: * "Invalid Arg" if N<0 or non-integer; * "Data Type" if N or F not real * ***************************************************** * * Examples: * * Testing for primality can be done by * IF N>1 AND FACTOR(N)=N THEN ... * * Complete factorization can be obtained as follows: * 10 DIM N,F,E @ STD @ INPUT N @ DISP STR$(N);"="; @ F=0 * 20 F=FACTOR(N,F) @ E=1 * 30 N=N/F @ IF NOT MOD(N,F) THEN E=E+1 @ GOTO 30 * 40 DISP STR$(F); @ IF E>1 THEN DISP "^"&STR$(E); * 50 IF N#1 THEN DISP "*"; @ GOTO 20 * 60 DISP @ END * ***************************************************** * * Algorithm: * * Try first 2, 3 and 5, then all but their multiples * until a factor is found or sqrt(N) exceeded * * If F<30, or F not specified, start with 2, * else start with 30*(F div 30)+1 * * Only every 8th divisor (those of form 30k+11) * is tested against sqrt(N) * ***************************************************** ARGERR EQU #0BF19 "Invalid Arg" ARGST- EQU #0E910 pop & test real number EXPR EQU #0F23C function return FLOAT EQU #1B322 integer to floating-point FNRTN4 EQU #0F238 function return IF12A EQU #0C739 integer/fraction split RJUST EQU #12AE2 floating-point to integer SPLITA EQU #0C6BF A to 15-form SQR15 EQU #0C534 square root ***************************************************** ID #5E MSG 0 POLL POLLH VER$ handler ENTRY factor CHAR #F function KEY 'FACTOR' TOKEN 9 ENDTXT ***************************************************** * Poll handler for VER$ POLLH ?B=0 B pVER$? GOYES VER$ RTNSXM VER$ C=R3 stk pt D1=C D1=D1- (Ve)-(Vs)-2 str length CD1EX A=R2 AVMEMS ?C=0? GOYES fOK yes noF A=0 W none or <0, substitute 0 fOK GOSBVL RJUST F into decimal integer * Next round F down to a multiple of 30 * i.e., compute F-(F mod 30) = 30*(F div 30) ASR W F div 10 B=A W dividend C=0 W divisor (3) here C=C+1 S start with 3E15 C=C+1 S (we only want remainder) C=C+1 S sub3 B=B-C W div 3 loop GONC sub3 B=B+C W CSR W ?C#0 W GOYES sub3 * now B contains (F div 10) mod 3 A=A-B W (F div 10)-((F div 10) mod 3) ASL W F-(F mod 30) = 30*(F div 30) R2=A put F aside for now GOSBVL ARGST- pop N, check it's real GOC argerr Inf -> Invalid Arg R0=A save orig. N * zero & one should be weeded out too * now compute sqrt(N) * we have original N in A GOSBVL SPLITA N to 15-form GOSBVL SQR15 SQRT(N) in (A,B) A=B M to 12-form * this simplified 12-to-15 form conversion fails if it is * too big, too small etc, but then it won't be used anyway GOSBVL RJUST into integer R1=A sqrt(N) in R1 ***************************************************** * All numbers not divisible by 2, 3 or 5 are of the * form 30k+m, where m is in {1,7,11,13,17,19,23,29} * The divisors are generated by adding 2, 4 or 6 * (from 2 to 3 also 1) to previous one * The addition is usually done in B field, * as in most cases only last digit can change * Of course 1-digit field would be better still, * but we can't spare P for that purpose: * In order to minimize division time, we do subtractions * in smallest possible field, which is WP with P * just above most significant digit of divisor * In main loop P is reserved for this purpose * N is also kept aligned with F to avoid unnecessary * rotations in the innermost division loop * This means we must detect when the number of digits in F * increases, and increment P and rotate N when it does; * this can be done with no cost in execution time by * keeping only one loop limit: either next power of 10 * above current divisor, or sqrt(N), whichever comes first * The limit is tested against only once in main loop, * i.e., after every 8th divisor, which means that * up to 7 divisions may be done unnecessarily, * but that takes less time than more tests would * Note that this makes it vital to have 30k+11 as the * first divisor in main loop, as it's the only one * that can have more digits than its predecessor ***************************************************** * locate decimal point of N A=R0 N GOSBVL IF12A find decimal point C=B W * now build 2 in B, to be used for incrementing divisor * and as result for N>=1E14 B=0 W B=B+1 B B=B+1 B * test if N is non-integer or too big for us ?P= 15 N>=1E14? GOYES toobig ?P= 14 exp<0? GOYES argerr ?C=0 WP fractional part zero? GOYES rjustN yes, N is OK argerr GOVLNG ARGERR "Invalid Arg" toobig A=B W if N>=1E14 it must be even GOTO retA so we return 2 * right-justify N rjustN CSR W P=P-1 GONC rjustN A=R2 F ?C<=A W F>N could cause endless loop GOYES prime this also catches N=0 * move P to msd of N; above left P=0 msdN P=P-1 look for non-zero digit ?C=0 P (N better be nonzero here!) GOYES msdN * the decimal point of N is marked with #F, * which is then used as terminating condition in division * thus we can rotate N without keeping the exponent SETHEX generate #Fs C=C-1 S to mark decimal point of N B=B-1 S and in B[S], to compare it with SETDEC * Then N must be rotated right until most significant digit * is at same position with F's msd * P must be left just past that position * We also need the smallest power of 10 above F msdF ?A#0 P non-zero digit? GOYES alignd yes, they're aligned CSRC rotate N P=P-1 next digit of F GONC msdF P=P+1 for F=0 alignd R0=C now N is as we want it P=P+1 just past msd D=0 W create the power of 10 D=D+1 P to be used as first loop limit C=R1 unless sqrt(N) is smaller ?C>D W GOYES Dok D=C W sqrt(N) is smaller, use it * now decide how to enter main loop * do 2, 3 &c need to be tried? (is F zero?) Dok C=R2 F, multiple of 30 C=C-1 WP 30k-1 GONC f31 jump into main loop * carry above means F=0 -> start with 2 C=0 B GOSUB add2 2 C=C+1 B GOSUB mod 3 GOSUB add2 5 GOSUB add2 7 GOC f11 (BET) into main loop ***************************************************** * no factor found -> N is a prime, return itself prime GOVLNG EXPR it's still on the stack ***************************************************** * Main loop * try next 8 divisors, * then test if sqrt(N) has been exceeded * Register usage: * R0 = N, msd aligned with F & dec pt marked with #F * C[0:P-1] = current divisor (right-justified integer) * D = limit: min(sqrt(N),10^ceil(log10(F))) * B = #F0000000000000002 next8 GOSUB mod 30k+11 GOSUB add2 30k+13 GOSUB add4 30k+17 GOSUB add2 30k+19 C=C+B WP may carry past tens digit GOSUB add2 30k+23 GOSUB add6 30k+29 f31 C=C+B WP may carry past tens GOSUB mod 30k+1 GOSUB add6 30k+7 f11 C=C+B B 30k+7 to 30k+9 C=C+B WP 30k+9 to 30k+11 ?C<=D WP GOYES next8 ***************************************************** * the addition to 30k+11 is the only one that can increase * the number of digits; then P must be incremented, * N rotated and new limit chosen CR1EX get sqrt(N) to C ?C<=D W past it already? GOYES prime yes, it's a prime DSL W next power of 10 ?C>D W is it less than sqrt(N)? GOYES use10 yes, use it as next limit D=C W else use sqrt(N) use10 CR1EX P=P+1 just past msd of divisor A=R0 shift more of N ASLC to be used at first subtraction R0=A in the division loop GONC next8 (BET) back to business ***************************************************** * mod: division subroutine * only remainder is computed * entry conditions: * R0 = N, msd aligned with F & dec pt marked with #F * C[0:P-1] = F (current divisor, right-justified integer) * F mustn't have more digits than N! * B = #F000000000000002 * returns with remainder in A, only if it is nonzero * always returns with carry set * additional entry points add2 &c increment divisor first * (in B field! do it yourself if you want more) add6 C=C+B B add4 C=C+B B add2 C=C+B B mod A=R0 N GONC subt B.E.T. shftN ASLC shift N subt A=A-C WP division loop GONC subt A=A+C WP ?A#B S all of N done? GOYES shftN no, continue ?A#0 WP remainder zero? RTNYES no, return * no remainder, factor found ***************************************************** * return factor A=C W 'twas left in C above C=RSTK don't return retA GOSBVL FLOAT into floating-point C=A W FNRTN4 wants it in C GOVLNG FNRTN4 END