'Subject: Approximate class number of an imaginary quadratic field
'         with Dirichlet's theorem.
'Author : Sjoerd.J.Schaper
'Date   : 21-12-2003
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
CONST lim = 131969 '                   prime limit
CONST N = 0, D = 1, p = 2
k = LargeInit(p, "clasimag")
DIM h AS DOUBLE, pi AS DOUBLE
pi = 4 * ATN(1)
CLS : tim# = TIMER: lin = 7

cx = PrimCeil
IF cx = 0 THEN
   Prints "primelist is empty - first run PrimFlgs.bas", 2
   GOTO eind
END IF

begin:
LOCATE 1, 1
PRINT "   +----------------+"
PRINT "   |  input N < 0,  |"
PRINT "   |  0 to exit.    |"
PRINT "   |                |"
PRINT "   +----------------+"
LOCATE 4, 7: INPUT ; "N"; Ni$
IF INSTR(Ni$, "'") THEN GOTO begin
IF VAL(Ni$) >= 0 THEN GOTO eind

LOCATE lin, 1
Readst N, Ni$
Printn N, "   N = ", "", 1

Letf p, -4
IF Cmp(p, N) THEN
   Dup N, D: t = 0
   Triald p, D, 4999 '                 is N quadratfrei?
   WHILE ExQ(D, k, p, t)
      IF k > 1 THEN
         Prints "square factor found", 2
         GOTO begin
      END IF
   WEND
END IF

Dup N, D
IF Mp2(N, 4) > 1 THEN Lsft D, 2 '      N (mod 4) = 2, 3 : D = 4N

h = EXP(Logf(D) * .5) / pi
Letf p, -3
IF Cmp(p, D) = 0 THEN h = h * 3 '      corrections
Inc p, -1
IF Cmp(p, D) = 0 THEN h = h * 2

sw = 0
DO
   c = Nxtprm(sw): Letf p, c
   k = Kronec(D, p)
   h = h * (c / (c - k)) '             Euler product formula
LOOP UNTIL c >= lim

Printn D, "h(", ") ", 0
Prints "approx" + STR$(h), 2
lin = CSRLIN: GOTO begin

eind:
LOCATE lin, 1
g$ = STR$(CSNG(TIMER - tim#))
Prints "Timer: " + g$ + " s", 1
CALL Term
SYSTEM

