Skip to content

Instantly share code, notes, and snippets.

/diophantine.cpp

Created Jun 15, 2017
Embed
What would you like to do?
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
typedef long long ll;
void solve1(int limit){
for(double x=1; x<limit; x++) {
for(double y=1; y<limit; y++) {
for(double z=1; z<limit; z++) {
if( abs(x/(y+z) + y/(x+z) + z/(y+x) - 4) <= 1E-15) {
printf("solution: x=%.0f, y=%.0f, z=%.0f\n",x,y,z);
}
}
}
}
}
void solve2(int limit){
for(ll x=1; x<limit; x++) {
for(ll y=1; y<limit; y++) {
for(ll z=1; z<limit; z++) {
if( x*(x+y)*(x+z) + y*(y+x)*(y+z) + z*(z+x)*(z+y) == 4*(y+z)*(x+z)*(x+y)) {
printf("solution: x=%lld, y=%lld, z=%lld\n",x,y,z);
}
}
}
}
}
void solve3(int limit){
for(ll x=1; x<limit; x++) {
for(ll y=x; y<limit; y++) {
for(ll z=y; z<limit; z++) {
if( x*(x+y)*(x+z) + y*(y+x)*(y+z) + z*(z+x)*(z+y) == 4*(y+z)*(x+z)*(x+y)) {
printf("solution: x=%lld, y=%lld, z=%lld\n",x,y,z);
}
}
}
}
}
// evaluates cubic at point x
ll eval(ll C3, ll C2, ll C1, ll C0, ll x){
return ((C3*x + C2)*x+C1)*x + C0;
}
// finds a root of cubic described by C3, C2, C1, C0 in interval [start, end)
// assumes cubic is monotone on this interval
ll findRoot(ll C3, ll C2, ll C1, ll C0, ll start, ll end){
if(end <= start) return -1;
while(end-start > 1){
ll mid = (start+end)/2;
if(eval(C3, C2, C1, C0, mid) > 0){
end = mid;
}else{
start = mid;
}
}
return start;
}
void solveCubic(ll C3, ll C2, ll C1, ll C0, ll roots[3]){
// make C3 positive
if (C3 < 0){
C3 = -C3;
C2 = -C2;
C1 = -C1;
C0 = -C0;
}
// compute derivative
ll D2 = 3*C3;
ll D1 = 2*C2;
ll D0 = C1;
// check discriminant
ll disc = D1*D1 - 4*D0*D2;
if(disc > 0){
ll crit1 = (ll)((-D1 - sqrt(disc))/(2*D2));
ll crit2 = (ll)((-D1 + sqrt(disc))/(2*D2));
// cubic is monotonically increasing from -infty to crit1 ...
roots[0] = findRoot(C3, C2, C1, C0, 0, crit1+1);
// ... monotonically decreasing from crit1 to crit2 ...
roots[1] = findRoot(-C3, -C2, -C1, -C0, crit1+1, crit2+1);
// ... and monotonically increasing from crit2 to infty
roots[2] = findRoot(C3, C2, C1, C0, crit2+1, C0+1);
}else if(disc <=0){
// cubic is monotonically increasing on all of R
// any positive integer root must lie in [0, C0]
roots[0] = findRoot(C3, C2, C1, C0, 0, C0+1);
}
}
void solve4(int limit){
ll rootsZ[3];
rootsZ[0] = 0;
rootsZ[1] = 0;
rootsZ[2] = 0;
for(ll x=1; x<limit; x++){
for(ll y=x;y<limit;y++){
solveCubic(1, -3*(x+y), -(3*x*x + 5*x*y + 3*y*y), (x*x + y*y - 4*x*y)*(x+y), rootsZ);
for(int i=0;i<3;i++){
ll z = rootsZ[i];
if(z > 0 && z < limit){
if( x*(x+y)*(x+z) + y*(y+x)*(y+z) + z*(z+x)*(z+y) == 4*(y+z)*(x+z)*(x+y)) {
printf("solution: x=%lld, y=%lld, z=%lld\n",x,y,z);
}
}
}
}
}
}
int main(){
solve1(1000);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment