/* BC program classnopos determines all reduced forms of a given positive discriminant d =0 or 1 (mod 4) */

/* absolute value of an integer n */

define abs(n){
	if(n>=0) return(n)
	return(-n)
}

/* 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)
}   

define abs(n){
	if(n>=0) return(n)
	return(-n)
}

/*  mod(a,b)=the least non-negative remainder when an integer a is divided by a
 positive integer b */

define mod(a,b){
	auto c
	c=a%b
	if(a>=0) return(c)
	if(c==0) return(0)
        return(c+b)	
}

/* int(a,b)=integer part of a/b, a, b integers, b != 0 */

define int(a,b){
	auto c
	c=sign(b)
	a=a*c
	b=b*c
	return((a-mod(a,b))/b)
}


/* 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

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

/*
 * This is a function for finding the period of the continued fraction
 * expansion of reduced quadratic irrational a=(u+sqrt(d))/v.
 * Here d is non-square, 1<(u+sqrt(d))/v, -1<(u-sqrt(d))/v<0.
 * The algorithm also assumes that v divides d-u*u and is based on K. Rosen,
 * Elementary Number theory and its applications, p.379-381 and Knuth's The art
 * of computer programming, Vol. 2, p. 359.  0 is returned if a is not reduced.
 */
define period(d,u,v){
	auto a,f,j,k,r,s,t
	f=sqrt(d)
	s=v
	r=u
	k=global_count /* set initially to 0 below in class_number(d) */

	for(j=k;1;j++){
		a=(f+u)/v
		u=a*v-u
		v=(d-u*u)/v
		globalarray_a[j]=u
		globalarray_b[j]=v
		if(u==r && v==s){
				break
		}
	}
	t=j+1-k
	global_count=global_count+t
	/*print "cycle-length=",t,"\n"*/
	return(t)
}

/* test(u,v) checks to see if (u,v) has already occurred in the
 * list globalarray_a[] && globalarray_b[].
 */
define test(u,v){
	auto i
	for(i=0;i<global_count;i++){
		if(u==globalarray_a[i] && v==globalarray_b[i]){
			return(0)
		}
	}
	return(1)
}

/* class_number(d) performs Lagrange's method on all reduced quadratic
 * irrationals (b+\sqrt(D))/2|c|, where 4*c divides D-b^2.
 * D is the discriminant D, the class-number h(d)
 * of the quadratic field Q(sqrt(d) is calculated, as well as the norm of
 * the fundamental unit.
 */
define class_number(d){
	auto a,b,c,e,f,g,h,k,l,r,x,class_number,w,z,parity_flag
	global_count=0
        class_number = 0;
	parity_flag=0
	w=squarefree_test(d)
	if(w==0){
		print d," is not squarefree\n"
		return(0)
	}else{ /* creates a fundamental discriminant */
		if((d-1)%4 != 0){
			d=4*d 
			e=2
		}else{
			e=1
		}
	}
	print "reduced forms of discriminant ",d," corresponding to ideal classes\n"
	f=sqrt(d)
	g=f/2
	for(a=1;a<=f;a++){
		for(b=e;b<=f;b=b+2){
			h=b^2-d
			i=2*a
			j=4*a
			if(h%j==0){
				if((a<=g && f-i<b) || (a>g && f-i>=b)){
					{
					     r=h/j
				             k=gcd(a,b)
					     l=gcd(k,abs(r))
					     if(l==1){
						c=(-h)/j   /* c > 0 here */
						x=test(b,2*c)
						if(x){
						   class_number=class_number+1
						   z=period(d,b,2*c)
					           if(z%2){
		 				 	parity_flag=1
						   }
			      		           print "reduced form [",class_number,"]: (",a,",",b,",",r,")\n"
						}
					     }
					}
				}
			}
	
		}
	}
	if(e==2){
		temp=1
		d=d/4
	}else{
		temp=4
	}
        if(parity_flag){
              print "x^2-",d,"*y^2=-",temp," has a solution\n"
        }else{
              print "x^2-",d,"*y^2=-",temp," has no solution\n"
        }
        print "h(",d,") = " 
        return(class_number)
}

/* class_number0(d) performs Lagrange's method on all reduced quadratic
 * irrationals (b+\sqrt(d))/2|c|, where 4*c divides d-b^2.
 * Here d>0 is not a perfect square and 0 or 1 mod(4).
 * the number h(d) of equivalence classes of binary quadratic forms of
 * discriminant d is returned. Also the solubility of x^2-d*y^2=-4 is 
 * determined.
 * This is basically the same program as class_number(d), except that in
 * the case of non-solubility of x^2-d*y^2=-4, we count the form (-a,b,c)
 * as well as (a,b,c), a>0 and this means we return twice the value that
 * class_number(d) would otherwise have returned.
 * Regarding this point, see G.B. Mathews, 'Theory of Numbers', pp. 80-81.
 */
define class_number0(d){
	auto a,b,c,e,f,g,h,k,l,r,x,class_number,w,z,parity_flag
	global_count=0
        class_number = 0;
	parity_flag=0
	f=sqrt(d)
	t=f*f
	if(t==d){
		print "d is a perfect square\n"
		return(0)
	}
	u=d%4
	if(u != 0 && u != 1){
		print "d is not 0 or 1 (mod 4)\n"
		return(0)
	}
	print "reduced forms (a,b,c) of discriminant ",d,"\n"
	if((d-1)%4 != 0){
		e=2
	}else{
		e=1
	}
	g=f/2
	for(a=1;a<=f;a++){
		for(b=e;b<=f;b=b+2){
			h=b^2-d
			i=2*a
			j=4*a
			if(h%j==0){
				if((a<=g && f-i<b) || (a>g && f-i>=b)){
					{
					     r=h/j
				             k=gcd(a,b)
					     l=gcd(k,abs(r))
					     if(l==1){
						c=(-h)/j   /* c > 0 here */
						x=test(b,2*c)
						if(x){
						   class_number=class_number+1
						   z=period(d,b,2*c)
					           if(z%2){
		 				 	parity_flag=1
						   }
			      		           print "reduced form [",class_number,"]: (",a,",",b,",",r,")\n"
						}
					     }
					}
				}
			}
	
		}
	}
	if(e==2){
		temp=1
	}else{
		temp=4
	}
	if(parity_flag){
              print "x^2-",d,"*y^2=-",temp," has a solution\n"
	}else{
              print "x^2-",d,"*y^2=-",temp," has no solution\n"
	      class_number=2*class_number
              print "There are reduced forms (-a,b,-c) as well\n"
	}
	print "h(",d,") = " 
       return(class_number)
}

/* class_number1(d) performs Lagrange's method on all reduced quadratic
 * irrationals (b+\sqrt(D))/2|c|, where 4*c divides D-b^2.
 * When performed on a fundamental discriminant D, the class-number h(d)
 * of the quadratic field Q(sqrt(d) is calculated, as well as the norm of
 * the fundamental unit.
 * No printing of reduced forms.
 * For use in table(m,n) below.
 */
define class_number1(d){
	auto a,b,c,e,f,g,h,k,l,r,x,class_number,w,z,parity_flag
	global_count=0
        class_number = 0;
	parity_flag=0
	/* creates a fundamental discriminant */
	if((d-1)%4 != 0){
		d=4*d 
		e=2
	}else{
		e=1
	}
	f=sqrt(d)
	g=f/2
	for(a=1;a<=f;a++){
		for(b=e;b<=f;b=b+2){
			h=b^2-d
			i=2*a
			j=4*a
			if(h%j==0){
				if((a<=g && f-i<b) || (a>g && f-i>=b)){
					{
					     r=h/j
				             k=gcd(a,b)
					     l=gcd(k,abs(r))
					     if(l==1){
						c=(-h)/j   /* c > 0 here */
						x=test(b,2*c)
						if(x){
						   class_number=class_number+1
						   z=period(d,b,2*c)
					           if(z%2){
		 				 	parity_flag=1
						   }
						}
					     }
					}
				}
			}
	
		}
	}
	if(e==2){
		d=d/4
	}
        if(parity_flag){
	        norm=-1
        }else{
	       norm=1
        }
        print "h(",d,")=",class_number,", N(eta)=",norm,"\n" 
        return(0)	
}

/* The following gives a table of h(d) for all squarefree d 
 * in the range 1<m<=d<=n<10^6. */
define table(m,n){
	auto d,h,t,limit
	limit=10^6
	if(m>n){
		print "m>n\n"
		return(0)
	}
	if(n>=limit){
		print "n >= ",limit,"\n"
		return(0)
	}
	if(m<=1 || n<=1){
		print "m or n <= 1\n"
		return(0)
	}
	for(d=m;d<=n;d++){
		h=squarefree_test(d)
		if(h==0){
			continue
		}else{
			t=class_number1(d)
		}
	}

}