/* BC file classnoneg */
/* gcd when m and n are >= 0 */
define gcd(m,n){
	auto a,b,c
	a=m
	if(n==0)return(a)
	b=n
	c=a%b
	while(c>0){
		a=b
		b=c
		c=a%b
	}
	return(b)
}   

/* This is Henri Cohen's Algorithm 5.3.5, p. 228. 
 * for finding the class number h(d) of binary quadratic forms 
 * of discriminant d, when d<0.
 * Here d=0(mod 4) or 1(mod 4).
 * If flag=1, we print only the primitive forms.
 * h(d) is returned in each case.
 * Davenport's Higher Arithmetic has a table of forms, 
 * which lists the imprimitive ones with an asterisk.
 * If d is the discriminant of an imaginary quadratic field K, 
 * then the primitive forms class-number h(d) is also the class number of K. 
 * We print the forms only if table_flag is nonzero.
 * eg. in table(m,n) below, we set table_flag=0.
 */
define class_number(d,flag,table_flag){
	auto a,b,bb,g,gg,h,q,t
	h=0
	g=1
	if(d%4==0){
		b=0
	}else{
		b=1
	}
	bb=sqrt(-d/3)
	
	while(b<=bb){
		q=(b^2-d)/4
		a=b
		if(a<=1){
			a=1
		}
/*		print "(",a,",",b,",",q/a,")\n"*/
		while(a^2<=q){
			if(q%a==0){
				t=q/a
				if(flag){
					gg=gcd(a,b)
					gg=gcd(gg,t)
					if(gg>1){
						g=0
					}
				}
				if(g==1){
                   		   if(a==b || a^2==q || b==0){
					if(table_flag){
					    print "(",a,",",b,",",t,")\n"
					}
					h=h+1
				   }else{
					if(table_flag){
					   print "(",a,",",b,",",t,")\n"
					   print "(",a,",",-b,",",t,")\n"
					}
					h=h+2
				   }
				}else{
					g=1
				}
			}
			a=a+1
		}
		b=b+2
	}
	return(h)
}

/* This BC program takes an integer d, |d| <1066 and outputs 1 if d is 
 * squarefree, 0 otherwise.
 */
prime[0]=2
prime[1]=3
prime[2]=5
prime[3]=7
prime[4]=11
prime[5]=13
prime[6]=17
prime[7]=19
prime[8]=23
prime[9]=29
prime[10]=31
prime[11]=37
prime[12]=41
prime[13]=43
prime[14]=47
prime[15]=53
prime[16]=59
prime[17]=61
prime[18]=67
prime[19]=71
prime[20]=73
prime[21]=79
prime[22]=83
prime[23]=89
prime[24]=97
prime[25]=101
prime[26]=103
prime[27]=107
prime[28]=109
prime[29]=113
prime[30]=127
prime[31]=131
prime[32]=137
prime[33]=139
prime[34]=149
prime[35]=151
prime[36]=157
prime[37]=163
prime[38]=167
prime[39]=173
prime[40]=179
prime[41]=181
prime[42]=191
prime[43]=193
prime[44]=197
prime[45]=199
prime[46]=211
prime[47]=223
prime[48]=227
prime[49]=229
prime[50]=233
prime[51]=239
prime[52]=241
prime[53]=251
prime[54]=257
prime[55]=263
prime[56]=269
prime[57]=271
prime[58]=277
prime[59]=281
prime[60]=283
prime[61]=293
prime[62]=307
prime[63]=311
prime[64]=313
prime[65]=317
prime[66]=331
prime[67]=337
prime[68]=347
prime[69]=349
prime[70]=353
prime[71]=359
prime[72]=367
prime[73]=373
prime[74]=379
prime[75]=383
prime[76]=389
prime[77]=397
prime[78]=401
prime[79]=409
prime[80]=419
prime[81]=421
prime[82]=431
prime[83]=433
prime[84]=439
prime[85]=443
prime[86]=449
prime[87]=457
prime[88]=461
prime[89]=463
prime[90]=467
prime[91]=479
prime[92]=487
prime[93]=491
prime[94]=499
prime[95]=503
prime[96]=509
prime[97]=521
prime[98]=523
prime[99]=541
prime[100]=547
prime[101]=557
prime[102]=563
prime[103]=569
prime[104]=571
prime[105]=577
prime[106]=587
prime[107]=593
prime[108]=599
prime[109]=601
prime[110]=607
prime[111]=613
prime[112]=617
prime[113]=619
prime[114]=631
prime[115]=641
prime[116]=643
prime[117]=647
prime[118]=653
prime[119]=659
prime[120]=661
prime[121]=673
prime[122]=677
prime[123]=683
prime[124]=691
prime[125]=701
prime[126]=709
prime[127]=719
prime[128]=727
prime[129]=733
prime[130]=739
prime[131]=743
prime[132]=751
prime[133]=757
prime[134]=761
prime[135]=769
prime[136]=773
prime[137]=787
prime[138]=797
prime[139]=809
prime[140]=811
prime[141]=821
prime[142]=823
prime[143]=827
prime[144]=829
prime[145]=839
prime[146]=853
prime[147]=857
prime[148]=859
prime[149]=863
prime[150]=877
prime[151]=881
prime[152]=883
prime[153]=887
prime[154]=907
prime[155]=911
prime[156]=919
prime[157]=929
prime[158]=937
prime[159]=941
prime[160]=947
prime[161]=953
prime[162]=967
prime[163]=971
prime[164]=977
prime[165]=983
prime[166]=991
prime[167]=997

/* tests d < 10^6 for squarefreeness */
define squarefree_test(d){
	auto i
	for(i=0;i<=167;i++){
		if(d%(prime[i]*prime[i])==0){
			return(0)
		}
	}
	return(1)
}

/* The following gives a table of h(-d) for all d
 * in the range 3=<m<=d<=n<10^6. */
   define table(m,n){
	auto d,r,s,temp,limit
	limit=10^6
	if(m>n){
           temp=n
           n=m
           m=temp
	}
	if(n>=limit){
		print "n >= ",limit,"\n"
		return(0)
	}
	if(m<3){
		print "m or n < 4\n"
		return(0)
	}
	for(d=m;d<=n;d++){
	    r=d%4
	    if(r==0 || r==3){
		s=class_number(-d,1,0)
		print "h(",-d,") = ",s,"\n"
            }
	}
   }