'******************************************************************************
'Subject: Arithmetic with integral binary quadratic forms.
'Author : Sjoerd.J.Schaper
'Date : 10-01-2006
'Code : QBasic, PDS, VBdos, FreeBasic
'You'll find an extended version for big integer-BASIC
here
'******************************************************************************
'This program is copyright (c) 2006 by the author. It is made available as is,
'and no warranty - about the program, its performance, or its conformity to any
'specification - is given or implied. It may be used, modified, and distributed
'freely, as long as the original author is credited as part of the final work.
'******************************************************************************
DEFLNG A-Z
TYPE qform
a AS DOUBLE
b AS DOUBLE
c AS DOUBLE
d AS DOUBLE
END TYPE
DECLARE FUNCTION Disc# (f AS qform)
'return discriminant(f)
DECLARE FUNCTION Prim% (f AS qform)
'return -1 if f is primitive
DECLARE FUNCTION IsId% (f AS qform)
'return -1 if f is principal
DECLARE SUB Pform (f AS qform)
'f:= principal form of Disc
DECLARE SUB Reduce (f AS qform)
'reduce form f
DECLARE SUB Normlz (f AS qform)
'normalize form f
DECLARE SUB Compos (f AS qform, g AS qform)
'f:= compound of forms f and g
DECLARE SUB Qpowr (f AS qform, k)
'raise form f to power k
DECLARE FUNCTION Regltr# (d)
'principal cycle and regulator (D > 0)
DECLARE FUNCTION Qrder (f AS qform)
'bf order of form f
DECLARE FUNCTION Bezout (a, b)
'Bezouts identity gcd(a,b) = ax + by
DECLARE SUB Printq (f AS qform)
'print form f
CONST MB = &H8000&, Fx = &H7FFFFFFF
DIM SHARED bD, drD AS DOUBLE, rD
'Discriminant, radical(D), floor(radical(D))
DIM SHARED pc(8191), tl AS INTEGER
'Principal cycle storage
DIM SHARED errsw AS INTEGER
'set in Sub Normlz()
DIM f AS qform, g AS qform, k
CLS : tim# = TIMER: lin = 8
begin:
PRINT : r = CSRLIN
IF r > 7 THEN lin = r
LOCATE 1, 1
PRINT " +--------------------------------+"
PRINT " | input binary quadratic form |"
PRINT " | f(a,b,c) and exponent k, |"
PRINT " | a = 0 to exit. |"
PRINT " | |"
PRINT " +--------------------------------+"
LOCATE 5, 7: INPUT ; "a"; gi$
r = INSTR(gi$, "f")
s = INSTR(gi$, "'")
IF r OR s THEN
IF r THEN f = g ' input "f" to re-enter f
IF s OR f.a = 0 THEN GOTO begin
ELSE
IF VAL(gi$) = 0 THEN GOTO eind
f.a = VAL(gi$): f.d = 0
INPUT ; " b"; gi$: f.b = VAL(gi$)
INPUT " c"; gi$: f.c = VAL(gi$)
END IF
LOCATE 5, 7: PRINT SPACE$(30) + "|";
LOCATE 5, 7: INPUT "k"; gi$: drD = VAL(gi$)
IF ABS(drD) > Fx THEN GOTO begin
errsw = 0: bD = 0: k = drD
drD = Disc(f)
LOCATE lin, 1: Printq f
IF ABS(drD) > Fx THEN ' preliminary tests
PRINT " too big D": GOTO begin
END IF
IF NOT Prim(f) THEN
PRINT " not primitive f": GOTO begin
END IF
IF drD < 0 AND f.a < 0 THEN
PRINT " negative definite f": GOTO begin
END IF
bD = drD: drD = SQR(ABS(drD)): rD = INT(drD)
IF rD * rD = bD THEN
PRINT " square D": GOTO begin
END IF
PRINT " D ="; bD; " rD ="; rD; ' Discriminant, radical(D)
g = f
IF k = 0 THEN ' find order(f)
IF bD > 0 THEN
PRINT " reg.="; CSNG(Regltr(bD));
END IF
k = Qrder(f)
PRINT " ord.="; k
'
IF bD > 0 THEN
k = 0
DO
f = g: k = k + 1
Qpowr f, k
IF errsw THEN EXIT DO
LOOP UNTIL ABS(f.a) = 1 ' identity
f = g
END IF
ELSE
PRINT
END IF
PRINT " k ="; k
Qpowr f, k
Printq f
GOTO begin
eind:
LOCATE lin, 1
PRINT "Timer:"; CSNG(TIMER - tim#); "s"
SYSTEM
FUNCTION Bezout (a, b)
DIM r(2, 1), i AS INTEGER, j AS INTEGER
IF a = 0 THEN
Bezout = ABS(b): b = SGN(b)
EXIT FUNCTION
END IF
'
i = 0: j = 1
r(0, i) = a: r(0, j) = b
r(1, i) = 1: r(1, j) = 0
r(2, i) = 0: r(2, j) = 1
DO
SWAP i, j
q = r(0, i) \ r(0, j) ' Euclidean steps
FOR s = 0 TO 2
r(s, i) = r(s, i) - q * r(s, j)
NEXT s
LOOP UNTIL r(0, i) = 0
'
IF r(0, j) < 0 THEN
FOR s = 0 TO 2
r(s, j) = -r(s, j)
NEXT s
END IF
a = r(1, j): b = r(2, j)
Bezout = r(0, j)
END FUNCTION
SUB Compos (f AS qform, g AS qform)
DIM t AS DOUBLE
x = f.a: y = g.a: u = x ' initialize
s = (f.b + g.b) * .5
t = f.b - s
'
v = Bezout(u, CLNG(g.a))
d = Bezout(s, v)
IF d > 1 THEN
x = x \ d: y = y \ d
END IF
'
f.a = x * y ' compose f, g
s = s * f.c MOD y
t = t * u * v MOD y
f.b = f.b - (t + s) * x * 2
f.c = (f.b * f.b - bD) / (f.a * 4)
f.d = f.d + g.d ' distance
Reduce f
END SUB
FUNCTION Disc# (f AS qform)
Disc = f.b * f.b - f.a * f.c * 4
END FUNCTION
FUNCTION IsId% (f AS qform)
IF bD < 0 THEN
i = ABS(f.a) = 1 ' identity f
ELSE
i = 0: x = ABS(f.a) * MB + f.b
FOR t = 0 TO tl ' search principal cycle
IF x = pc(t) THEN i = -1: EXIT FOR
NEXT t
END IF
IsId = i
END FUNCTION
SUB Normlz (f AS qform)
DIM a AS DOUBLE, sg AS INTEGER
IF bD > 0 THEN
a = (drD - f.b) / (drD + f.b) ' Shanks's distance
f.d = f.d + LOG(ABS(a)) * .5
END IF
'
sg = SGN(f.a): a = ABS(f.a)
IF bD < 0 OR a > rD THEN
s = sg * INT((a - f.b) / (a * 2))
ELSE ' normalizing integer
s = sg * INT((rD - f.b) / (a * 2))
END IF
a = f.a * s
f.c = f.c + (f.b + a) * s
f.b = f.b + a * 2
errsw = ABS(Disc(f) - bD) > 1D-19 ' D-check
END SUB
SUB Pform (f AS qform)
f.a = 1: f.b = 1: f.d = 0
IF bD > 0 THEN f.b = rD
IF (bD AND 1) XOR (f.b AND 1) THEN f.b = f.b - 1
f.c = (f.b * f.b - bD) / 4
END SUB
FUNCTION Prim% (f AS qform)
a = f.a: b = f.b: c = f.c
Prim = Bezout(Bezout(a, b), c) = 1
END FUNCTION
SUB Printq (f AS qform)
IF errsw THEN
PRINT " overflow": EXIT SUB
END IF
g$ = "Qf = (" + LTRIM$(STR$(f.a))
g$ = g$ + "," + LTRIM$(STR$(f.b))
g$ = g$ + "," + LTRIM$(STR$(f.c))
IF bD > 0 THEN
g$ = g$ + "," + LTRIM$(STR$(CSNG(f.d)))
END IF
PRINT g$ + ")"
END SUB
SUB Qpowr (f AS qform, k)
DIM p AS qform
t = k
IF t < 0 THEN
Reduce f
f.b = -f.b: t = -t ' opposite form
ELSEIF t = 0 THEN
Pform f: EXIT SUB
END IF
IF t = 1 THEN
Reduce f: EXIT SUB
END IF
WHILE (t AND 1) = 0
Compos f, f: t = t \ 2 ' initialize
WEND
'
p = f
WHILE t > 1 ' R=>L binary exponentiation
Compos f, f: t = t \ 2
IF t AND 1 THEN Compos p, f
IF errsw THEN EXIT SUB
WEND
f = p
END SUB
FUNCTION Qrder (f AS qform)
DIM p AS qform
Qrder = 0
p = f: Reduce p: t = 1
DO
IF IsId(p) THEN
Qrder = t: EXIT DO
ELSE
t = t + 1 ' brute force search
END IF
Compos p, f
IF errsw THEN EXIT DO
LOOP
END FUNCTION
SUB Reduce (f AS qform)
DIM r AS INTEGER
DO
IF bD < 0 THEN
r = ABS(f.b) > f.a
IF r OR f.a > f.c THEN
SWAP f.a, f.c
f.b = -f.b: Normlz f
ELSE
IF ABS(f.b) = f.a OR f.a = f.c THEN
f.b = ABS(f.b)
END IF
EXIT DO
END IF
ELSE
r = ABS(rD - 2 * ABS(f.c)) >= f.b
IF r OR f.b > rD THEN
SWAP f.a, f.c
f.b = -f.b: Normlz f
ELSE
EXIT DO
END IF
END IF
IF errsw THEN EXIT DO
LOOP
END SUB
FUNCTION Regltr# (d)
DIM f AS qform
Pform f: tl = -1
x0 = MB + f.b
DO
SWAP f.a, f.c
f.b = -f.b: Normlz f
x = ABS(f.a) * MB + f.b
tl = tl + 1: pc(tl) = x ' store principal cycle
LOOP UNTIL x = x0
Regltr = f.d
END FUNCTION