Counting Square-full Solutions to x+y=z

Heath-Brown

https://arxiv.org/pdf/2005.04180v1

// Example program
#include <iostream>
#include <string>

using namespace std;

long long bin;

long long euler_phi(long long n) {
  long long ans = n;
  for (long i = 2; i * i <= n; i++)
    if (n % i == 0) {
      ans = ans / i * (i - 1);
      while (n % i == 0) n /= i;
    }
  if (n > 1) ans = ans / n * (n - 1);
  return ans;
}

long long gcd(long long a, long long b) {
  long long i = __builtin_ctz(a);
  long long j = __builtin_ctz(b);
  long long k = min(i, j);
  long long d;
  b >>= j;
  while (a) {
    a >>= i;
    d = b - a;
    i = __builtin_ctz(d);
    if (a < b) b = a;
    if (d < 0) a = -d;
    else a = d;
  }
  return b << k;
}

long long lcm(long long a, long long b){
    return a*b/gcd(a,b);
}

long long quickpow(long long a, long long b, long long m) {
    a %= m;
    b %= bin;
   long long sum = 1;
   while (b > 0) {
       if (b & 1)
           sum = sum * a % m;
       a = a * a % m;
       b >>= 1;
   }
   return sum;
}

bool is_ok(long long n){
    bin = euler_phi(n);
    long long MAX = lcm(n,bin);
    for(int i = 1; i < MAX; i++){
        if(gcd(i,n)==1){
            if((quickpow(i,i,n)-1 % n == 0) && ((i-1)%n!=0)){
                cout << i << endl;
                return 0;
            }
        }
    }
    cout<<"pass"<<endl;
    return 1;
}

int main()
{

  long long k;
  cin>>k;
  is_ok(k);
  main();
  
  /*
  for(int i = 1; i < 1000; i++){
      if(is_ok(2*i)){cout << 2*i << " ";}
    }
    */
  return 0;
}

#include<math.h>
#include<iostream>
#include<time.h>
#include<iomanip>

#define sup 2048
#define version "v1.1"
#define debug 1
#define overload_num 3
#define tolerence_ring 1

using namespace std;

const double sp2 = 0.6079271018540267;
const double p2s = 1.644934066848226;
const double ps = 0.5235987755982988;
const double ps_sqrt = 0.723601254558267659363;

long x_len = 1 + floor(sqrt(sup) * ps_sqrt);
long y_len = 1 + floor(x_len / sqrt(2));
long real_x_len = 0;
long real_y_len = 0;
long polygon[sup][2] = {{0}};
bool circle[sup][sup] = {{0}};
int vert = 0;
long stable_amount = 0;
long adjust_amount = 0;

long long gcd(long long a, long long b) {//faster gcd
  long long i = __builtin_ctz(a);
  long long j = __builtin_ctz(b);
  long long k = min(i, j);
  long long d;
  b >>= j;
  while (a) {
    a >>= i;
    d = b - a;
    i = __builtin_ctz(d);
    if (a < b) b = a;
    if (d < 0) a = -d;
    else a = d;
  }
  return b << k;
}

bool is_coprime(long a, long b){//just to make it faster
    if(min(a, b) == 1){return 1;}
    if((a & 1) == 0 && (b & 1) == 0){return 0;}
    if((a % 3) == 0 && (b % 3) == 0){return 0;}
    if((a % 5) == 0 && (b % 5) == 0){return 0;}
    return (1 == gcd(a, b));
}

int near(long double q, bool is_up){
    if(is_up){return floor(q) + 1;}
    else{return ceil(q) - 1;}
}

long double diff(long double q, bool is_up){
    return abs(q - near(q, is_up));
}

// | 1 a b |
// | 1 c d |
// | 1 e f |
short triangle_area(long a[2], long c[2], long e[2]){
    long long s = c[0]*e[1]-c[1]*e[0]-a[0]*e[1]+a[1]*e[0]+a[0]*c[1]-a[1]*c[0];
	if(s > 0){return 1;}
	if(s < 0){return -1;}
	else{return 0;}
}

long double triangle_area_n(long a[2], long c[2], long e[2]){
    return c[0]*e[1]-c[1]*e[0]-a[0]*e[1]+a[1]*e[0]+a[0]*c[1]-a[1]*c[0];
}

int modnum(int i){//addon
    return (i + vert) % vert;
}


void make_circle(){
	long double rad = sqrt(vert) * ps_sqrt;
	long double rad_sq = vert * ps;
	long y = floor(rad / sqrt(2));
	real_y_len = y;
	real_x_len = rad;
	long current_x = y;
	for(int i = y; i > 0; i--){
		while(current_x * current_x + i * i <= rad_sq){current_x ++;}
		current_x--;
		circle[current_x][i] = 1;
		circle[i][i] = 1;
	}
	bool have_met = 0;
	for(int j = real_y_len - 1; j > 0; j--){//here real_y_len-1 is to prevent that the first row have only one 1
		for(int i = 1; i <= real_x_len; i++){
		    if(circle[i][j] && have_met){
		        circle[i][j] = is_coprime(i, j);
		        have_met = 0;
		        continue;
		    }
			if(circle[i][j]){
			    have_met = 1;
			    circle[i][j] = is_coprime(i, j);
			}
			if(have_met && !circle[i][j] && is_coprime(i, j)){circle[i][j] = 1;}
		}
	}
}

//now the major promotion must be drawn here and in the func auto_adjust one
void make_convex(/*int vert*/){
	int quad = ceil(vert / 4.0 + overload_num);//overload protect
	long len = quad * (quad + 1) / 2;
	for(int i = 0; i < quad; i++){
		polygon[i][0] = i * (i + 1) / 2;
		polygon[i][1] = quad - i;
		polygon[quad + i][0] = len + i;
		polygon[quad + i][1] = i * (i + 1) / 2;
		polygon[2 * quad + i][0] = len + quad - i * (i + 1) / 2;
		polygon[2 * quad + i][1] = len + i;
		polygon[3 * quad + i][0] = quad - i;
		polygon[3 * quad + i][1] = len + quad - i * (i + 1) / 2;
	}
	for(int i = vert; i < 4 * quad; i++){//cancel overload
	    //something can be propelled here forwards the montgomery method
	    //by picking random points to eliminate then do the adjustments
        polygon[i][0] = 0;
        polygon[i][1] = 0;
    }
}




void adjust(int vertnum){
    //the trivial adjustment
    //but it seemes always to work well awa
    int v1 = modnum(vertnum - 2);
    int v2 = modnum(vertnum - 1);
    int v3 = modnum(vertnum);
    int v4 = modnum(vertnum + 1);
    int v5 = modnum(vertnum + 2);
    long double k;
    long double b;
    long st[2];
    //a slightly better soloution but it seems to be nonsense
    //find the shorter edge of v2v4
    long dist = min(abs(polygon[v2][0] - polygon[v4][0]), abs(polygon[v2][1] - polygon[v4][1]));
    bool b_p = (dist != abs(polygon[v2][0] - polygon[v4][0]));
    if(dist != 0){
    	//little confusing
   		//unkind for future debug
    	k = ((polygon[v2][!b_p] - polygon[v4][!b_p]) * 1.0) / ((polygon[v2][b_p] - polygon[v4][b_p]) * 1.0);
    	normal:
	    b = polygon[v2][!b_p] - polygon[v2][b_p] * k;
	    //[!b_p] = k * [b_p] + b
	    //is v3 upon v2v4?
	    long double len = (polygon[v3][!b_p] - k * polygon[v3][b_p] - b);
	    bool facing = (len > 0);
	    short test_1 = triangle_area(polygon[v1], polygon[v2], polygon[v3]);
	    short test_2 = triangle_area(polygon[v5], polygon[v4], polygon[v3]);
	    len = abs(len);
	    long left = min(polygon[v2][b_p], polygon[v4][b_p]);
	    long double area = triangle_area_n(polygon[v2], polygon[v3], polygon[v4]);
	    //long right = max(polygon[v2][b_p], polygon[v4][b_p]);
    	for(int i = 0; i < dist; i++){
    	    st[b_p] = left + i;
    	    st[!b_p] = near(k * st[b_p] + b, facing);
    		if((diff(k * st[b_p] + b, facing) < len)){
    		    if((triangle_area(polygon[v1], polygon[v2], st) == test_1) && (triangle_area(polygon[v5], polygon[v4], st) == test_2)){
    		        polygon[v3][0] = st[0];
    		        polygon[v3][1] = st[1];
    		    }
    		}
		}
		if(triangle_area_n(polygon[v2], polygon[v3], polygon[v4]) == area){//when auto stop
		    stable_amount++;
		}
		else{stable_amount = 0;}
	}
	if(dist == 0){
	    dist = max(abs(polygon[v2][0] - polygon[v4][0]), abs(polygon[v2][1] - polygon[v4][1]));
	    b_p = !b_p;
	    k = 0;
	    goto normal;
	}
}

void ring_adjust(){
    for(int i = 0; i < vert; i++){
        adjust(i);
    }
    long stx = polygon[0][0], sty = polygon[0][1];
    for(int i = 0; i < vert; i++){
        stx = min(stx, polygon[i][0]);
        sty = min(sty, polygon[i][1]);
    }//do the shift
    for(int i = 0; i < vert; i++){
        polygon[i][0] -= stx;
        polygon[i][1] -= sty;
    }
}

void auto_adjust(){
    while(stable_amount < tolerence_ring * vert){//no idea but to type sth like this uwu
        adjust_amount++;
        ring_adjust();
    }
}

void output_polygon(){
    for(int i = 0; i < vert; i++){
		cout << polygon[i][0] << "," << polygon[i][1] << "  ";
	}
}

unsigned long long polygon_area() {
    unsigned long long area = 0;
    for (int i = 0; i < vert; i++) {
        int j = (i + 1) % vert;
        area += polygon[i][0] * polygon[j][1] - polygon[i][1] * polygon[j][0];
    }
    return area;
}

int main(){
	if(!debug){
		cout << "ver: " << version << endl;
		cout << "input the vertex number of the polygon: " << endl;
		cin >> vert;
		make_convex();
	// 	output_polygon();
	   	long long tic = clock();
		auto_adjust();
		long long toc = clock();
	// 	cout << endl;
		output_polygon();
		cout << endl;
		cout << "adjust amount: " << adjust_amount << endl; 
		unsigned long long area = polygon_area();
		cout << "polygon area: " << fixed << setprecision(1) << area / 2.0 << endl;
		cout << "ratio: " << fixed << setprecision(8) << area / (2 * pow(vert, 3)) << endl;
		cout << "time used: " << (toc - tic) / 1000 << " ms" << endl;
		return 0;
	}
	else{
		cout << "ver: " << version << endl;
		cout << "testing make_circle" << endl;
		cout << "input the vertex number of the polygon: " << endl;
		cin >> vert;
		make_circle();
		for(int j = real_y_len; j > 0; j--){
			for(int i = 1; i <= real_x_len; i++){
				cout << circle[i][j] << " ";
			}
			cout << endl;
		}
		return 0;
	}
}

Last update: 2026/1/9 22:21