FANDOM


Program
    ZFACTOR2
; Author:    Rob Gaebler
; Date:        1/16/04
; Purpose:    Find the prime factors of a positive integer


; This is AFACTOR2, the BASIC program that runs ZFACTOR2:

; Delvar YDelvar [A]Input "NUM=?",X
; Asm(prgmZFACTOR2
; If Y
; Then
; Clrhome
; Pause [A]


; For best viewing of this file, set tab size to 8
 
.NOLIST
#define equ .equ
#define EQU .equ
#define end .end
#define END .end
#include "ti83PlusAsm.inc"
.LIST
.org 9D95h



;--------------------------------------
 
Mod_Equates:     
 
mod1        equ 135
mod2        equ 143
mod3        equ 119
mod4        equ 19
mod5        equ 23
mod6        equ 29
mod7        equ 31
mod8        equ 37
mod9        equ 41
mod10        equ 43
mod11        equ 47
mod12        equ 53
 
modlast        equ 53
mod_sum        equ mod1+mod2+mod3+mod4+mod5+mod6+mod7+mod8+mod9+mod10+mod11+mod12  ;=720     

;--------------------------------------
 
 
 
Start:
B_CALL(_delres)        ; STATVARS is used as storage
 
B_CALL(_rclx)        ; store X to op1
 
 
; check for good input
 
ld hl,op1+1
B_CALL(_ckint)
ret nz            ; make sure the input is an integer
 
B_CALL(_ckop1pos)
ret nz            ; make sure the input is positive
 
B_CALL(_ckop1fp0)
ret z            ; make sure the input is not 0
 
 
; input is good
 
 
 
ld de,N
B_CALL(_movfrop1)    ; copies op1 to N
call Compute_Div_Bound    ; bound for trial division
 
ld hl,INITPRIME        ; list of primes
ld b,list_len        ; length of list
 
First_Primes:            ; divide by primes on the primes list
push hl
push bc
 
ld a,(hl)        ; get prime
 
B_CALL(_setxxop2)    ; store a to op2
xor a
call Check_Factor    ; see if op2 is a factor
pop bc
pop hl
jr c,Final_Residue    ; jump if op2 > op1
 
inc hl            ; on to next prime
djnz First_Primes    ; do loop again
 
 
ld hl,INITPRIME        ; list of primes
ld de,RES_TABLE        ; residue table
ld b,list_len        ; length of primes list
 
Init_Res_Table:            ; ld a,(p-1)/2
ld a,(hl)
srl a
ld (de),a
inc hl
inc de
djnz Init_Res_Table
 
 
B_CALL(_op2set1)    ; op2 is the divisor
 
;    The following "sieve" sifts out potential divisors of the number we are factoring.
;    This way, we don't have to test as many numbers to see if they are factors.  The
;    sieve works by throwing out numbers that have small prime factors, since those
;    factors were already taken care of by the First_Primes loop.
 
Pre_Test:
xor a            ; a is amount to increment divisor by
 
Test:
add a,2            ; add 2 to increment value
push af            ; save increment value
 
ld hl,RES_TABLE+1
ld de,INITPRIME+1
ld b,list_len-1
xor a
 
Test_Div:     
dec (hl)
jr nz,No_Prime
ld a,(de)
ld (hl),a
 
No_Prime:     
inc l            ; inc hl
inc e            ; inc de
djnz Test_Div
 
Done_Test:     
ld b,a
pop af
 
inc b
djnz Test
 
; if jump not taken, we have a trial divisor
 
B_CALL(_setxxop1)    ; a is value to increment divisor by
B_CALL(_fpadd)        ; add op2 to op1 to get divisor
B_CALL(_op1toop2)    ; store divisor to op2
 
 
xor a
call Check_Factor    ; find out if op2 is a factor
 
jr c,Final_Residue    ; jump if past trial division bound
 
call Del_Check        ; check to see if Del was pressed
 
 
jr nz,Pre_Test        ; find another trial divisor
 
ret            ; end program if Del was pressed
 
 
 
Final_Residue:
call Recall_N        ; set op1 at correct value
B_CALL(_op2set1)    ; set op2=1
B_CALL(_cpop1op2)    ; check to see if op1=1
ld a,1
call nz,Sto_Fact    ; display final factor and quit
 
ret
 
;---------------------------------------------------------------------------------------
     
; Input to Check_Factor:
; a = 0
; N = number to factor
; op2 = factor to check
 
 
Is_Factor:
ld de,N
B_CALL(_movfrop1)    ; N = op1
B_CALL(_op2toop1)    ; op1 = FACTOR
call formdisp        ; display factor
call Compute_Div_Bound
call Recall_Factor    ; op2 = FACTOR
pop af
inc a            ; power of factor
 
Check_Factor:            ; a = 0 when this is called
push af
 
 
ld de,FACTOR     
ld hl,Op2
B_CALL(_mov9b)        ; FACTOR = op2
 
call Recall_N        ; op1 = N
B_CALL(_fpdiv)        ; op1 = N/FACTOR
ld hl,op1+1        ; hl points to exponent of op1
B_CALL(_ckint)
jr z,Is_Factor
 
; op2 was not changed by _fpdiv
 
B_CALL(_op2toop1)    ; op1 = FACTOR
pop af
or a            ; check if a = 0
call nz,Sto_Fact
 
call Recall_Factor    ; op2 = FACTOR
ld hl,DIVBOUND
B_CALL(_mov9toop1)    ; op1 = DIVBOUND
 
B_CALL(_cpop1op2)
ret nc
 
; now, end program or go on to sieve
 
ld a,(BOUND_TYPE)
or a            ; see if BOUND_TYPE = 0
scf
ret nz            ; if not 0, end program
pop hl            ; strip away return address
 
 
;    Now the trial division ends and the setup for the quadratic residue sieve begins     
 
 
Sieve:     
B_CALL(_op1toop2)    ; op2 = B
call Recall_N        ; op1 = N
B_CALL(_fpdiv)
B_CALL(_fpadd)
B_CALL(_timespt5)
call Copy_BOUND        ; SIEVEBOUND = (N+B^2)/(2B)
jp Resume
 
Mid_Storage:     
.ds mod_sum+Start-Mid_Storage

Resume:     
ld de,SQR_TABLE
call Create_Residue_Table
 
 
Copy_And_Shift:            ; copy each quadratic residue table, and subtract a
            ; constant for each table, depending on N (mod M), for
            ; each M: SQR-N (mod M)
ld hl,SQR_TABLE
ld de,SHIFTED_TABLE
push de
ld bc,mod_sum
ldir            ; copy the quadratic residue tables
 
ld hl,SIEVE_MODULI    ; sieve moduli
pop bc            ; address of table to shift
 
Create_Shifted_Tables:
ld a,(hl)        ; load M
push hl            ; modulus M
push bc            ; address of table to shift
 
push af
call Recall_N        ; op1 = N
pop af
call Op1moda        ; e = op1 (mod a)
 
pop bc            ; address of table to shift
pop hl            ; modulus M
ld a,(hl)        ; (hl) = M
            ; a is a counter of the number of quadratic residues
 
Shift_Table:     
push af
ld a,(bc)        ; unshifted quadratic residue
sub e            ; e = N (mod M), amount to shift by
jr nc,Shifted_Residue
 
add a,(hl)        ; add back M to make it positive
 
Shifted_Residue:        ; store the shifted quadratic residue
ld (bc),a
inc bc
pop af
dec a            ; a = counter of quadratic residues for current modulus
jr nz,Shift_Table    ; jump if not done shifting this table
 
; bc now points to beginning of the quadratic residue table for the next modulus
 
inc hl            ; on to the next modulus
ld de,SIEVE_MODULI+modnum
B_CALL(_cphlde)        ; changes no registers
jr nz,Create_Shifted_Tables    ; z if de = hl, signaling that all the tables have been
            ; shifted
 
Short_Residue_Table:     
ld hl,Insert
 
; Insert:
; inc b
; srl b
 
ld (hl),04h        ; inc b
inc hl
ld (hl),0CBh
inc hl
ld (hl),38h        ; srl b
 
 
ld de,SQR_TABLE_SHORT
call Create_Residue_Table
 
 
Create_Bit_Table:
; a bit is set each time a number appears in both the quadratic residues and
; shifted quadratic residues tables, and eight copies are made of each bit table
 
di            ; disable interrupts because of exx's
 
; set 1:
ld hl,SQR_TABLE_SHORT    ; quadratic residues mod mod1
ld de,SHIFTED_TABLE    ; shifted quadratic residues table
 
; set 1
; (hl) = start of table of quadratic residues
; (de) = start of table of shifted quadratic residues
 
exx
 
; set 2:
ld hl,BIT_TABLE        ; area to store results to
ld bc,SIEVE_MODULI    ; moduli to use
 
; set 2
; (hl) = start of bit table
; (bc) = table of sieve moduli
 
Bit_Loop1:            ; loops modnum times, once for each modulus
ld a,(bc)        ; M
ld e,a            ; M
push hl            ; start of bit table for this modulus
 
Bit_Loop2:            ; loops ceiling(M/8) times
ld d,10000000b
ld (hl),0
 
; set 2
; (bc) = M, the modulus for the current bit table being created
; (hl) = area to store bit table for this modulus
; e = counter for position in shifted residue table
; d = 128
 
Bit_Loop3:            ; loops 8 times, once for each bit in a byte
ld a,(bc)        ; (bc) = M
inc a
rrca
 
; set 2
; (bc) = M
; (hl) = byte to store result of this loop to
; d = loop counter and bit masker
; e = counter for position in shifted residue table
; a = (M+1)/2
 
exx            ; set 1
 
ld c,a
ld b,0
ld a,(de)        ; shifted quadratic residue (mod M)
push hl
 
; set 1
; (hl) = start of quadratic residue table for this modulus
; (de) = current position in shifted quadratic residue table for this modulus
; bc = (M+1)/2
; a = (de) = shifted quadratic residues
 
Bit_Loop4:            ; for a given byte in the shifted residue table, see if
            ; there is an equal byte in the (unshifted) residue table
            ; to determine whether X^2-N is a quadratic residue mod M
             
cpir            ; compares until bc = 0 or a = (hl)
 
; set 1
; (de) = current position in shifted quadratic residue table for this modulus
; a = (de) = shifted quadratic residue
 
exx            ; set 2
 
ld a,d            ; a power of 2
jr z,Store_Bit        ; z if a match was found 
 
xor a            ; no match was found
 
Store_Bit:
 
; set 2
; (hl) = byte to store result of this loop to
; (bc) = M
; d = loop counter and bit masker
; e = counter for position in shifted residue table
; a = byte containing bit to store to bit table
 
or (hl)            ; combine a with bit table so far
ld (hl),a        ; load to bit table
 
exx            ; set 1
 
pop hl            ; quadratic residues
 
; set 1
; (hl) = start of quadratic residue table for this modulus
; (de) = current position in shifted quadratic residue table for this modulus
 
inc de            ; shifted quadratic residues--on to the next one
 
exx            ; set 2
 
; set 2
; (hl) = byte to store result of this loop to
; (bc) = M
; d = loop counter and bit masker
; e = counter for position in shifted residue table
 
dec e            ; e originally set to M
jr z,Copy_Bits        ; reset e to M and decrease de' by M if e = 0
 
srl d            ; reduce power of 2 by 1
jr nz,Bit_Loop3
 
inc hl            ; on to next byte in bit table
 
; set 2
; (hl) = byte in bit table to store next results to
; (bc) = M
; d = 0
; e = counter for position in residue table
 
jr Bit_Loop2
 
Copy_Bits:            ; make 8 more copies (the last is discarded) of the bit
            ; table for this modulus
; set 2
; (hl) = last byte in bit table so far
; (bc) = M
; d = shift counter
 
inc hl            ; first empty byte in bit table
dec d            ; equivalent to srl d
ld a,d
pop de            ; start of bit table for this modulus
push bc            ; (bc) = M
push af            ; a = shift counter
ld a,(bc)        ; M
ld b,a            ; M
 
Make_Copy:            ; copy bit table one byte at a time
pop af            ; a = shift counter
push af
ld c,a            ; shift counter
ld a,(de)        ; M
ld (hl),a        ; M
xor a
 
Shift_Copy:            ; shift one byte left a few bits
sla (hl)        ; shift left one bit
rla            ; shift left one bit
srl c            ; shift counter
jr nz,Shift_Copy
 
dec hl
or (hl)
ld (hl),a        ; shifted into previous byte in table
inc hl
inc hl            ; on to next byte to copy to
inc de            ; on to next byte to copy from
djnz Make_Copy
 
pop bc
pop bc            ; M
ex de,hl
ld a,(bc)        ; M
 
cp modlast        ; modlast is the last modulus
jr z,End_Loop1        ; jump if done with all the moduli
 
exx            ; set 1
 
; set 1
; (hl) = start of quadratic residue table for this modulus
; (de) = start of shifted quadratic residue table for this modulus
; a = M
 
ld b,a            ; a = M
inc a
rrca
call Hlplusa        ; adds (M+1)/2 to get the address of the next quadratic
            ; residue table
 
; set 1
; (hl) = start of quadratic residue table for next modulus
; (de) = start of shifted quadratic residue table for next modulus
; b = M
 
 
exx            ; set 2
 
inc bc            ; address of current modulus M
 
; set 2
; (hl) = start of bit table for next modulus
; (bc) = next modulus
; de = M
 
jr Bit_Loop1        ; do it all again with a new modulus M
 
 
End_Loop1:            ; done with Create_Bit_Table!
 
ei            ; enable interrupts
 
call Recall_N        ; op1 = N
 
ld a,16d        ; every 16th bit will be counted
call Op1moda        ; de = op1 (mod 16)
ld hl,STEP_TABLE-1
add hl,de
 
ld c,(hl)
inc hl
ld b,(hl)
ld (STEP1),bc
 
call Recall_N        ; op1 = N
 
B_CALL(_sqroot)
ld de,SQRT
B_CALL(_movfrop1)
B_CALL(_op2set8)
B_CALL(_fpdiv)
B_CALL(_intgr)
B_CALL(_op2set8)
B_CALL(_fpmult)        ; op1 = 8*floor(sqrt(N)/8)
 
; now op1 is the greatest multiple of 8 less than the square root of N, the
; place where X starts
 
ld hl,BIT_TABLE_A
ld de,BIT_TABLE_B
ld bc,mod_sum
ldir
 
ld a,(STEP1)
ld hl,BIT_TABLE_A
 
call Rotate_Bit_Tables
 
call Copy_X
 
 
ld a,(STEP2)
ld hl,BIT_TABLE_B
 
call Rotate_Bit_Tables
 
call Recall_X        ; op1 = X = greatest multiple of 8 less than square root
            ; of N, plus STEP1
 
Create_Place_Tables:
ld b,modnum
ld hl,SIEVE_MODULI
ld de,PLACES        ; counter for position in shifted bit tables, for sieve
 
Place_Loop:
ld a,(hl)        ; M
neg            ; -M
ld (de),a        ; -M
inc hl
inc de
djnz Place_Loop
 
 
; From label Sieve until here has been setup for the main part of the factoring
; program.  Next comes the main part, the quadratic residue sieve.  There are two main
; sieves, SIEVE_A and SIEVE_B.  Each one tests for a particular value of X mod 8.  Both
; sieves consist of several smaller sieves (12 in this case).  A bit must pass each of
; these smaller sieves in order for a number X to actually be tested by the
; Test_Congruence routine.  The "and (hl)" command keeps track of which bits pass through
; all the sieves.  Most of the time, no bits make it through all the sieves.  Each pass
; through SIEVE_A and SIEVE_B actually tests 16 different bits, and essentially 64
; different values of X.  SIEVE_A tests the values of X that are STEP1 mod 8 and SIEVE_B
; tests the values of X that are STEP1+STEP2 mod 8.  (There are only two residue classes
; mod 8 that need to be checked.)  Read the text file for more explanation of the
; mathematics involved.
 
 
Sieve_A:     
xor a            ; reset carry
ld b,a            ; b = 0
dec a            ; a = 11111111b
 
ld hl,PLACE_MOD1    ; offset in bit table
ld c,(hl)
ld hl,TABLE_MOD1_A-256+mod1
add hl,bc        ; correct position in bit table
and (hl)        ; pass 8 bits through the first sieve
jr z,Sieve_B
 
ld hl,PLACE_MOD2
ld c,(hl)
ld hl,TABLE_MOD2_A-256+mod2
add hl,bc
and (hl)
jr z,Sieve_B
 
ld hl,PLACE_MOD3
ld c,(hl)
ld hl,TABLE_MOD3_A-256+mod3
add hl,bc
and (hl)
jr z,Sieve_B
 
ld hl,PLACE_MOD4
ld c,(hl)
ld hl,TABLE_MOD4_A-256+mod4
add hl,bc
and (hl)
jr z,Sieve_B
 
ld hl,PLACE_MOD5
ld c,(hl)
ld hl,TABLE_MOD5_A-256+mod5
add hl,bc
and (hl)
jr z,Sieve_B
 
ld hl,PLACE_MOD6
ld c,(hl)
ld hl,TABLE_MOD6_A-256+mod6
add hl,bc
and (hl)
jr z,Sieve_B
 
ld hl,PLACE_MOD7
ld c,(hl)
ld hl,TABLE_MOD7_A-256+mod7
add hl,bc
and (hl)
jr z,Sieve_B
 
ld hl,PLACE_MOD8
ld c,(hl)
ld hl,TABLE_MOD8_A-256+mod8
add hl,bc
and (hl)
jr z,Sieve_B
 
ld hl,PLACE_MOD9
ld c,(hl)
ld hl,TABLE_MOD9_A-256+mod9
add hl,bc
and (hl)
jr z,Sieve_B
 
ld hl,PLACE_MOD10
ld c,(hl)
ld hl,TABLE_MOD10_A-256+mod10
add hl,bc
and (hl)
jr z,Sieve_B
 
ld hl,PLACE_MOD11
ld c,(hl)
ld hl,TABLE_MOD11_A-256+mod11
add hl,bc
and (hl)
jr z,Sieve_B
 
ld hl,PLACE_MOD12
ld c,(hl)
ld hl,TABLE_MOD12_A-256+mod12
add hl,bc
and (hl)
call nz,Test_Congruence        ; test a value of X
 
 
Sieve_B:     
; b = 0
ld a,0FFh
 
ld hl,PLACE_MOD1
ld c,(hl)
ld hl,TABLE_MOD1_B-256+mod1
add hl,bc
and (hl)
jr z,Increment
 
ld hl,PLACE_MOD2
ld c,(hl)
ld hl,TABLE_MOD2_B-256+mod2
add hl,bc
and (hl)
jr z,Increment
 
ld hl,PLACE_MOD3
ld c,(hl)
ld hl,TABLE_MOD3_B-256+mod3
add hl,bc
and (hl)
jr z,Increment
 
ld hl,PLACE_MOD4
ld c,(hl)
ld hl,TABLE_MOD4_B-256+mod4
add hl,bc
and (hl)
jr z,Increment
 
ld hl,PLACE_MOD5
ld c,(hl)
ld hl,TABLE_MOD5_B-256+mod5
add hl,bc
and (hl)
jr z,Increment
 
ld hl,PLACE_MOD6
ld c,(hl)
ld hl,TABLE_MOD6_B-256+mod6
add hl,bc
and (hl)
jr z,Increment
 
ld hl,PLACE_MOD7
ld c,(hl)
ld hl,TABLE_MOD7_B-256+mod7
add hl,bc
and (hl)
jr z,Increment
 
ld hl,PLACE_MOD8
ld c,(hl)
ld hl,TABLE_MOD8_B-256+mod8
add hl,bc
and (hl)
jr z,Increment
 
ld hl,PLACE_MOD9
ld c,(hl)
ld hl,TABLE_MOD9_B-256+mod9
add hl,bc
and (hl)
jr z,Increment
 
ld hl,PLACE_MOD10
ld c,(hl)
ld hl,TABLE_MOD10_B-256+mod10
add hl,bc
and (hl)
jr z,Increment
 
ld hl,PLACE_MOD11
ld c,(hl)
ld hl,TABLE_MOD11_B-256+mod11
add hl,bc
and (hl)
jr z,Increment
 
ld hl,PLACE_MOD12
ld c,(hl)
ld hl,TABLE_MOD12_B-256+mod12
add hl,bc
and (hl)
jp nz,Test_Congruence_B
 
Increment:     
ld hl,INC_VAL
inc (hl)        ; increment the increment value
jr z,Inc_X        ; increment X
 
Inc_Tables:     
 
Inc_MOD1:
ld hl,PLACE_MOD1
inc (hl)
jr z,Reset_MOD1
 
Inc_MOD2:
inc l
inc (hl)
jr z,Reset_MOD2
 
Inc_MOD3:
inc l
inc (hl)
jr z,Reset_MOD3
 
Inc_MOD4:
inc l
inc (hl)
jr z,Reset_MOD4
 
Inc_MOD5:
inc l
inc (hl)
jr z,Reset_MOD5
 
Inc_MOD6:
inc l
inc (hl)
jr z,Reset_MOD6
 
Inc_MOD7:
inc l
inc (hl)
jr z,Reset_MOD7
 
Inc_MOD8:
inc l
inc (hl)
jr z,Reset_MOD8
 
Inc_MOD9:
inc l
inc (hl)
jr z,Reset_MOD9
 
Inc_MOD10:
inc l
inc (hl)
jr z,Reset_MOD10
 
Inc_MOD11:
inc l
inc (hl)
jr z,Reset_MOD11
 
Inc_MOD12:
inc l
inc (hl)
jp nz,Sieve_A
 
Reset_MOD12:
ld (hl),256-mod12
jp Sieve_A
 
Reset_MOD1:
ld (hl),256-mod1
jr Inc_MOD2
 
Reset_MOD2:
ld (hl),256-mod2
jr Inc_MOD3
 
Reset_MOD3:
ld (hl),256-mod3
jr Inc_MOD4
 
Reset_MOD4:
ld (hl),256-mod4
jr Inc_MOD5
 
Reset_MOD5:
ld (hl),256-mod5
jr Inc_MOD6
 
Reset_MOD6:
ld (hl),256-mod6
jr Inc_MOD7
 
Reset_MOD7:
ld (hl),256-mod7
jr Inc_MOD8
 
Reset_MOD8:
ld (hl),256-mod8
jr Inc_MOD9
 
Reset_MOD9:
ld (hl),256-mod9
jr Inc_MOD10
 
Reset_MOD10:
ld (hl),256-mod10
jr Inc_MOD11
 
Reset_MOD11:
ld (hl),256-mod11
jr Inc_MOD12
 
 
Inc_X:
ld hl,4000h        ; 16384d
B_CALL(_setxxxxop2)
B_CALL(_fpadd)        ; increment X
 
call Del_Check
ret z            ; end program
 
ld hl,SIEVEBOUND
B_CALL(_mov9toop2)
B_CALL(_cpop1op2)    ; X-SIEVEBOUND
jr c,Inc_Tables
 
Prime:                ; N is prime
call Recall_N
jp Last_Factor
 
 
Test_Congruence_B:
push af            ; a indicates which factor(s) to check
ld a,(STEP2)
B_CALL(_setxxop2)
B_CALL(_fpadd)        ; add STEP2 to get correct X
pop af
 
call Test_Congruence    ; see if there are any factors
 
ld a,(STEP2)
B_CALL(_setxxop2)
B_CALL(_fpsub)        ; subtract back STEP2
 
jp Increment        ; increment X,then return to sieve--no factors were found
 
 
 
;---------------------------------------------------------------------------------------;
;                                            ;
;                  /----------------\                    ;
;                  |  Sub-routines  |                    ;
;                  \----------------/                    ;
;                                            ;
;---------------------------------------------------------------------------------------;
 
Test_Congruence:
push af            ; a indicates which factor(s) to check
ld hl,INC_VAL
ld a,(hl)
ld (hl),0        ; reset X increment counter
ld l,a
ld h,64
call Ltimesh
B_CALL(_setxxxxop2)
B_CALL(_fpadd)        ; add increment to X
pop af
 
 
ld b,8            ; 8 bits to test
 
Mini_Sieve:
rla            ; bits for which X to test
push bc            ; loop counter
push af
 
jr nc,End_Mini
 
 
Test_X:                ; test if X is a solution to N = X^2 - Y^2
            ; store last factors if they are found
            ; max error in expression is 3*10^-5 < .5, so this
            ; routine works
call Copy_X        ; X = op1
 
ld hl,SQRT
B_CALL(_mov9toop2)    ; op2 = sqrt(N)
B_CALL(_fpadd)        ; op1 = X+sqrt(N)
B_CALL(_sqroot)        ; op1 = sqrt(X+sqrt(N))
ld de,OPA
B_CALL(_movfrop1)    ; opA = sqrt(X+sqrt(N))
 
ld hl,SQRT
B_CALL(_mov9toop2)    ; op2 = sqrt(N)
call Recall_X        ; op1 = X
B_CALL(_cpop1op2)
jr c,End_Mini        ; c if op2 > op1
 
B_CALL(_fpsub)        ; op1 = X-sqrt(N)
B_CALL(_sqroot)        ; op1 = sqrt(X-sqrt(N))
ld hl,OPA
B_CALL(_mov9toop2)    ; op2 = sqrt(X+sqrt(N))
B_CALL(_fpmult)        ; op1 = sqrt(X^2-N)
B_CALL(_int)        ; round to nearest integer
 
B_CALL(_op1toop2)    ; op2 = Y
call Recall_X        ; op1 = X
B_CALL(_fpadd)        ; op1 = X+Y
B_CALL(_op1toop2)    ; op2 = X+Y
 
call Recall_N        ; op1 = N
B_CALL(_fpdiv)        ; op1 = N/(X+Y)
 
ld hl,op1+1
B_CALL(_ckint)        ; is (X+Y) a prime factor of N?
 
jr z,Sieve_Factor
 
call Recall_X
 
End_Mini:     
B_CALL(_op2set8)
B_CALL(_fpadd)        ; add 8 to X to get next value to test
 
pop af
pop bc
djnz Mini_Sieve
 
B_CALL(_op2set60)
ld hl,op2+2
ld (hl),64h        ; op2 = 64
B_CALL(_fpsub)
ret
 
;--------------------------------------
 
Sieve_Factor:     
; strip away the stack to exit the program
pop hl
pop hl
pop hl
 
call Copy_Factor    ; FACTOR = op1 = X-Y
 
B_CALL(_cpop1op2)    ; is op1 = op2? (op2 = X+Y)
ld a,2            ; power of 2 if N is a square
jr z,Sto_Fact        ; jump if N is a square
 
call Last_Factor    ; store and display factor
call Recall_Factor    ; op2 = X-Y
call Recall_N        ; op1 = N
B_CALL(_fpdiv)        ; op1 = X+Y
 
Last_Factor:
ld a,1
 
; fall through to Sto_Fact
 
;--------------------------------------
 
; Input:
; Op1 is a factor of the number
; a is the power of the factor (1 to 99)

; Output:
; factor (op1) stored in left column of matrix [A]
; power (a) stored in right column of matrix [A]
; Y=1
 
Sto_Fact:
push af
 
call Copy_Factor    ; copy factor
B_CALL(_zeroop1)    ; clear op1
 
 
ld hl,op1        ; copy var name for matrix [A] into op1
ld (hl),matobj
inc hl
ld (hl),tvarmat
inc hl
ld (hl),tmata
B_CALL(_chkfindsym)    ; look it up and store address to de
jr nc,Resize        ; if it is there, resize it
ld hl,0102h        ; 1 by 2
B_CALL(_creatermat)    ; create it
 
 
push de            ; address of [A]
jr Sto_Mat
 
Resize:     
push de
ex de,hl        ; store address of matrix [A] to hl
inc hl
ld d,(hl)        ; number of rows
inc d            ; increase number of rows by 1
ld e,2            ; 2 columns
ex de,hl        ; store new dimensions of [A] to hl
B_CALL(_redimmat)    ; redimension matrix [A]
 
Sto_Mat:
ld hl,FACTOR
B_CALL(_mov9toop1)    ; copy factor back to op1
pop de            ; address of [A]
inc de
ld a,(de)        ; row number
ld b,a
ld c,1            ; column number
dec de
 
push de            ; address of [A]
push bc            ; element of [A]
B_CALL(_puttomat)    ; op1 to (b,c) in [A]
pop bc            ; element of [A]
pop de            ; address of [A]
inc c            ; column number
pop af            ; power of factor
di
exx            ; push bc, push de
B_CALL(_setxxop1)    ; stores a to op1
exx            ; pop de, pop bc
ei
B_CALL(_puttomat)    ; op1 to (b,c) in [A]
 
 
B_CALL(_op1set1)
B_CALL(_stoy)        ; indicates that [A] exists
ret
 
;--------------------------------------
 
formdisp:
B_CALL(_eraseeol)
ld a,14
B_CALL(_formreal)
ld hl,op3
B_CALL(_puts)
 
B_CALL(_newline)
ret
 
_newline    .equ 452eh     
 
;--------------------------------------
 
Create_Residue_Table:        ; create a table of quadratic residues for each sieve
            ; modulus: 0^2, 8^2, 16^2, ..., (8M-8)^2 (mod M)
ld hl,SIEVE_MODULI-1
 
Quadratic_Residues:     
 
inc hl            ; increment sieve modulus
ld b,(hl)        ; sieve modulus M
 
Insert:                ; nop's will be overwritten second time routine is called
nop            ; inc b
nop            ; srl b
nop
 
xor a            ; start with 0 squared
 
Square_Mod_M:     
push bc            ; counter of residues
push de            ; quadratic residue table for M
push af            ; number to square and reduce mod M
push hl            ; (hl) contains sieve modulus M
 
 
ld h,a
ld l,a
call Ltimesh        ; hl = h*l
 
pop de            ; (de) = M
push de
ld a,(de)        ; M
 
 
Dividehlbya:            ; square divided by M
ld c,a
sub a
ld b,16
 
Div_Lbl_1:
add hl,hl
rla
jr c,Div_Lbl_2
 
cp c
 
jr c,Div_Lbl_3
 
Div_Lbl_2:
sub c
inc l
 
Div_Lbl_3:
djnz Div_Lbl_1
 
 
pop hl            ; (hl) = M
pop bc            ; b = number to square
pop de            ; quadratic residue table for M
ld (de),a        ; quadratic residue mod M (remainder from division by M)
 
ld a,b            ; number to square
add a,8
cp (hl)            ; modulus M
jr c,Hop
 
sub (hl)        ; keep a below M
 
Hop:     
inc de
pop bc            ; counter of residues
djnz Square_Mod_M
 
ld a,(hl)
cp modlast
jr nz,Quadratic_Residues
 
ret
 
;--------------------------------------
 
; Input:
; op1 = 8*floor(sqrt(N)/8) [+STEP1 possibly]
; hl = bit table to shift
; a = amount to step op1 by (STEP1 or STEP2)
 
Rotate_Bit_Tables:        ; shift each bit table so that the first bit in the table
            ; is the first bit to be tested
 
push hl            ; bit table
B_CALL(_setxxop2)
B_CALL(_fpadd)        ; increment op1 by a
pop hl            ; bit table
 
ld de,SIEVE_MODULI    ; M
 
Rotate:
; divide op1 by (de) = M
push hl            ; shifted bit table
push de            ; modulus M
ld a,(de)        ; M
 
call Op1moda        ; op4 = op1, a = op1 (mod a)
push af
B_CALL(_op4toop1)    ; restore op1
pop af
 
ld bc,4000h
pop hl            ; M
 
Divide_By_64:            ; keep subtracting 64 until result is 0 (mod M)
inc c
sub b            ; subtract b = 64
jr nc,Divide_By_64
 
dec c
add a,b            ; add back 64
 
jr z,Found_Inverse
 
add a,(hl)        ; M
jr Divide_By_64
 
Found_Inverse:            ; shift the bit table c bytes circularly to the left
ld b,a            ; 0
ex de,hl        ; (de) = M
pop hl            ; shifted bit table
push de            ; M
push bc            ; number of shifts
push hl            ; shifted bit table
ld a,(de)
sub c            ; a = M-c
ld de,TABLE_COPY
 
inc c            ; ldir once more than needed in case bc=0
ldir
dec hl            ; recover correct hl
 
ld c,a            ; M-c
pop de            ; shifted bit table
 
ldir
 
ld a,(hl)        ; save start of bit table for next modulus
ld hl,TABLE_COPY
pop bc
 
inc c            ; ldir once more than needed in case bc=0
ldir
dec de            ; recover correct de
ld (de),a        ; restore start of bit table for next modulus     

ex de,hl        ; hl = start of shifted bit table for next modulus
pop de            ; M
ld a,(de)        ; M
inc de
cp modlast
jr nz,Rotate
 
ret
 
;--------------------------------------
 
; Input:
; op1
; a
;
; Output:
; op4 = op1
; a = de = op1 (mod a)
 
Op1moda:
ld h,0
ld l,a
B_CALL(_setxxxxop2)    ; doesn't change op1
B_CALL(_op2toop6)
 
B_CALL(_op1toop4)
 
B_CALL(_fpdiv)        ; doesn't change op2
B_CALL(_intgr)        ; round down
 
B_CALL(_op6toop2)
B_CALL(_fpmult)
 
B_CALL(_op4toop2)
B_CALL(_invsub)        ; op2-op1
B_CALL(_convop1)    ; store op1 to de, LSB to a
ret
 
;--------------------------------------
 
Hlplusa:     
add a,l
ld l,a
ld a,0
adc a,h
ld h,a
ret
 
;--------------------------------------
 
: Input:
; hl

; Output:
; hl=h*l

; Used:
; a,b,d,e,h,l
; b=0
; d=0
; e=l that was input

; Condition bits:
; s not affected
; z not affected
; h affected
; p/v not affected
; n=0
; c=0

; Minimum time: 354 clock cycles
; Maximum time: 402 clock cycles
; Average time: 378 clock cycles
; Size: 16 bytes
 
Ltimesh:
ld a,h
ld e,l
ld hl,0
ld d,l
ld b,8
 
Mult_Loop:     
add hl,hl
rla
jr nc,End_Mult_Loop
 
add hl,de
 
End_Mult_Loop:     
djnz Mult_Loop
 
ret
 
;--------------------------------------
 
Del_Check:     
ld a,0BFh        ; load the column with Del
out (1),a        ; write a value to keyport
in a,(1)        ; read a value from keyport
cp kdelold        ; value to compare to (127)
ret
 
kdelold        equ 7Fh        ; constant for Delete key
 
;--------------------------------------
 
; The DIVBOUND is an arbitrary bound which I have chosen to (nearly) maximize the
; program's speed.  DIVBOUND is called "B" in the text file.  If N < 1000000, the program
; only uses trial division.  If 1000000 < N < 1000000000, the program uses trial division
; up to the cube root of N.  (As explained in the text file, B must always be AT LEAST
; the cube root of N.)  If N > 1000000000, the program uses trial divisionup to
; sqrt(N/1000).  The program would be slightly faster (2% or so) if the bound were
; sqrt(N/1300), but that would make the program larger.
 
Compute_Div_Bound:     
B_CALL(_op2set1)
ld hl,op2+1        ; exponent
ld (hl),89h        ; op2=1000000000

push hl
call Recall_N        ; op1 = N
B_CALL(_cpop1op2)
pop hl
jr nc,Bound_C
 
ld (hl),86h        ; op2=1000000

 
B_CALL(_cpop1op2)
jr nc,Bound_B
 
Bound_A:            ; N < 1000000, DIVBOUND = sqrt(N)
B_CALL(_sqroot)
ld hl,BOUND_TYPE
inc (hl)        ; (hl) is not 0
jr Sto_Bound
 
Bound_B:            ; N < 1000000000, DIVBOUND = cuberoot(N)
B_CALL(_op2set3)
B_CALL(_op1exop2)    ; documentation wrong on _xrooty
B_CALL(_xrooty)        ; cube root of op2
jr Sto_Bound
 
Bound_C:            ; N > 1000000000, DIVBOUND = sqrt(N/1000)
ld (hl),83h        ; op2 = 1000
B_CALL(_fpdiv)
B_CALL(_sqroot)
 
Sto_Bound:            ; fall through to Copy_Bound
 
;--------------------------------------
;     
;    op1 and op2 copy routines
;     
;--------------------------------------
 
Copy_Bound:
ld de,DIVBOUND        ; = SIEVEBOUND
B_CALL(_movfrop1)
ret
 
;--------------------------------------
 
Copy_X:
Copy_Factor:
ld de,FACTOR        ; FACTOR and X use same space
B_CALL(_movfrop1)
ret
 
;--------------------------------------
 
Recall_N:
ld hl,N
B_CALL(_mov9toop1)
ret
 
;--------------------------------------
 
Recall_X:     
ld hl,X
B_CALL(_mov9toop1)
ret
 
;--------------------------------------
 
Recall_Factor:
ld hl,FACTOR
B_CALL(_mov9toop2)
ret
 
;--------------------------------------
;     
;    Equates, tables, and other data
;     
;--------------------------------------




 
INITPRIME:
 
; If the list is to be extended past 17 numbers, then the _setxxop1 must be changed in
; the Done_Test area of the program.  If the list is to be extended past 25 numbers, then
; the _setxxop2 in the First_Primes loop must be changed
 
 
.db 2d        ;1
.db 3d        ;2
.db 5d        ;3
.db 7d        ;4
.db 11d        ;5
.db 13d        ;6
.db 17d        ;7
.db 19d        ;8
.db 23d        ;9
.db 29d        ;10
.db 31d        ;11
.db 37d        ;12
.db 41d        ;13
.db 43d        ;14
.db 47d        ;15
.db 53d        ;16
.db 59d        ;17
 
 
;    The next few primes:
;
;    .db 61d        ;18
;    .db 67d        ;19
;    .db 71d        ;20
;    .db 73d        ;21
;    .db 79d        ;22
;    .db 83d        ;23
;    .db 89d        ;24
;    .db 97d        ;25

list_len    equ 17d        ; the number of primes to use
 
;--------------------------------------
 
BOUND_TYPE:     
.ds 1
 
INC_VAL:     
.ds 1

;--------------------------------------
 
SIEVE_MODULI:
.db mod1
.db mod2
.db mod3
.db mod4
.db mod5
.db mod6
.db mod7
.db mod8
.db mod9
.db mod10
.db mod11
.db mod12
 
modnum        equ $-SIEVE_MODULI     
 
 
;--------------------------------------
 
RES_TABLE:    ; equ STATVAR+9S
.ds list_len+9+9+9+9+1+1+11+1+(mod_sum+modnum)/2
N        equ RES_TABLE+list_len
SQRT         equ N+9
FACTOR        equ SQRT+9
X        equ FACTOR
OPA        equ X+9
DIVBOUND    equ OPA+9
SIEVEBOUND    equ DIVBOUND
STEP1        equ SIEVEBOUND+9
STEP2        equ STEP1+1
SHADOW        equ STEP2+1     
PLACES        equ SHADOW+6
PLACE_MOD1    equ PLACES
PLACE_MOD2    equ PLACE_MOD1+1
PLACE_MOD3    equ PLACE_MOD2+1
PLACE_MOD4    equ PLACE_MOD3+1
PLACE_MOD5    equ PLACE_MOD4+1
PLACE_MOD6    equ PLACE_MOD5+1
PLACE_MOD7    equ PLACE_MOD6+1
PLACE_MOD8    equ PLACE_MOD7+1
PLACE_MOD9    equ PLACE_MOD8+1
PLACE_MOD10    equ PLACE_MOD9+1
PLACE_MOD11    equ PLACE_MOD10+1
PLACE_MOD12    equ PLACE_MOD11+1
 
SQR_TABLE_SHORT    equ PLACE_MOD12+1
 
TABLE_COPY:    ;equ TEMPSWAPAREA     
.ds mod_sum
 
SQR_TABLE    equ Start        ; #1
;.ds mod_sum
 
SHIFTED_TABLE:    ;equ SAVESSCREEN    ; #2
.ds mod_sum
 
BIT_TABLE    equ SQR_TABLE        ; #3 --    use space for #1
 
BIT_TABLE_A    equ BIT_TABLE        ; #4 -- use space for #1 or #2
TABLE_MOD1_A    equ BIT_TABLE_A
TABLE_MOD2_A    equ TABLE_MOD1_A+mod1
TABLE_MOD3_A    equ TABLE_MOD2_A+mod2
TABLE_MOD4_A    equ TABLE_MOD3_A+mod3
TABLE_MOD5_A    equ TABLE_MOD4_A+mod4
TABLE_MOD6_A    equ TABLE_MOD5_A+mod5
TABLE_MOD7_A    equ TABLE_MOD6_A+mod6
TABLE_MOD8_A    equ TABLE_MOD7_A+mod7
TABLE_MOD9_A    equ TABLE_MOD8_A+mod8
TABLE_MOD10_A    equ TABLE_MOD9_A+mod9
TABLE_MOD11_A    equ TABLE_MOD10_A+mod10
TABLE_MOD12_A    equ TABLE_MOD11_A+mod11
 
 
BIT_TABLE_B    equ SHIFTED_TABLE    ; #5 -- use space for #1 or #2 (not #4)
TABLE_MOD1_B    equ BIT_TABLE_B
TABLE_MOD2_B    equ TABLE_MOD1_B+mod1
TABLE_MOD3_B    equ TABLE_MOD2_B+mod2
TABLE_MOD4_B    equ TABLE_MOD3_B+mod3
TABLE_MOD5_B    equ TABLE_MOD4_B+mod4
TABLE_MOD6_B    equ TABLE_MOD5_B+mod5
TABLE_MOD7_B    equ TABLE_MOD6_B+mod6
TABLE_MOD8_B    equ TABLE_MOD7_B+mod7
TABLE_MOD9_B    equ TABLE_MOD8_B+mod8
TABLE_MOD10_B    equ TABLE_MOD9_B+mod9
TABLE_MOD11_B    equ TABLE_MOD10_B+mod10
TABLE_MOD12_B    equ TABLE_MOD11_B+mod11



;--------------------------------------

; The number on the right is N mod 16.  The numbers on the left tell what X must be
; (mod 8) for X^2-N to be a quadratic residue; X must be congruent to STEP1 or
; STEP1+STEP2 (mod 8)
 
STEP_TABLE:     
.db 1,6        ; 1 mod 16
.db 2,4        ; 3 mod 16
.db 3,2        ; 5 mod 16
.db 0,4        ; 7 mod 16
.db 3,2        ; 9 mod 16
.db 2,4        ; 11 mod 16
.db 1,6        ; 13 mod 16
.db 0,4        ; 15 mod 16
 
; N (mod 16)    X (mod 16)
;     1    1  7  9  15
;     3    2  6  10 14
;     5    3  5  11 13
;     7    0  4  8  12
;     9    3  5  11 13
;     11    2  6  10 14
;     13    1  7  9  15
;     15    0  4  8  12
 



 
.end
END

Ad blocker interference detected!


Wikia is a free-to-use site that makes money from advertising. We have a modified experience for viewers using ad blockers

Wikia is not accessible if you’ve made further modifications. Remove the custom ad blocker rule(s) and the page will load as expected.