#include <stdio.h>
#include <math.h>
#include <stdlib.h>
 
#define EPSILON 1e-6
#define MAX_ITER 100
 
double f1(double x, double y) {
    return x * x + y * y - 4; 
}
 
double f2(double x, double y) {
    return x * x * x * x + x - y*y + 1;
}
 
double df1_dx(double x, double y) { return 2 * x; }
double df1_dy(double x, double y) { return 2 * y; }
double df2_dx(double x, double y) { return 4 * x * x * x + 1; }
double df2_dy(double x, double y) { return -2 * y; }
 
int main() {
    double x = 1, y = 3;
    int iter = 0;
 
    printf("Initial guess: x = %.5f, y = %.5f\n", x, y);
 
    while (iter < MAX_ITER) {
 
        double F1 = f1(x, y);
        double F2 = f2(x, y);
 
        if (fabs(F1) < EPSILON && fabs(F2) < EPSILON) {
            printf("Converged after %d iterations.\n", iter);
            printf("Solution: x = %.5f, y = %.5f\n", x, y);
            return 0;
        }
 
        double J11 = df1_dx(x, y);
        double J12 = df1_dy(x, y);
        double J21 = df2_dx(x, y);
        double J22 = df2_dy(x, y);
 
        double det = J11 * J22 - J12 * J21;
        if (fabs(det) < 1e-12) {
            printf("Jacobian determinant too small. Method fails.\n");
            return 1;
        }
 
        double dx = (-F1 * J22 + F2 * J12) / det;
        double dy = (-J11 * F2 + J21 * F1) / det;
 
 
        printf("iter %2d: x=%.10f  y=%.10f  f1=%.10f  f2=%.10f  dx=%.10f  dy=%.10f  x_next=%.10f  y_next=%.10f\n",
       iter, x, y, F1, F2, dx, dy, x + dx, y + dy);
 
 
        x += dx;
        y += dy;
        iter++;
    }
 
    printf("Did not converge within %d iterations.\n", MAX_ITER);
    return 1;
}

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KIAojZGVmaW5lIEVQU0lMT04gMWUtNgojZGVmaW5lIE1BWF9JVEVSIDEwMAogCmRvdWJsZSBmMShkb3VibGUgeCwgZG91YmxlIHkpIHsKICAgIHJldHVybiB4ICogeCArIHkgKiB5IC0gNDsgCn0KIApkb3VibGUgZjIoZG91YmxlIHgsIGRvdWJsZSB5KSB7CiAgICByZXR1cm4geCAqIHggKiB4ICogeCArIHggLSB5KnkgKyAxOwp9CiAKZG91YmxlIGRmMV9keChkb3VibGUgeCwgZG91YmxlIHkpIHsgcmV0dXJuIDIgKiB4OyB9CmRvdWJsZSBkZjFfZHkoZG91YmxlIHgsIGRvdWJsZSB5KSB7IHJldHVybiAyICogeTsgfQpkb3VibGUgZGYyX2R4KGRvdWJsZSB4LCBkb3VibGUgeSkgeyByZXR1cm4gNCAqIHggKiB4ICogeCArIDE7IH0KZG91YmxlIGRmMl9keShkb3VibGUgeCwgZG91YmxlIHkpIHsgcmV0dXJuIC0yICogeTsgfQogCmludCBtYWluKCkgewogICAgZG91YmxlIHggPSAxLCB5ID0gMzsKICAgIGludCBpdGVyID0gMDsKIAogICAgcHJpbnRmKCJJbml0aWFsIGd1ZXNzOiB4ID0gJS41ZiwgeSA9ICUuNWZcbiIsIHgsIHkpOwogCiAgICB3aGlsZSAoaXRlciA8IE1BWF9JVEVSKSB7CgogICAgICAgIGRvdWJsZSBGMSA9IGYxKHgsIHkpOwogICAgICAgIGRvdWJsZSBGMiA9IGYyKHgsIHkpOwoKICAgICAgICBpZiAoZmFicyhGMSkgPCBFUFNJTE9OICYmIGZhYnMoRjIpIDwgRVBTSUxPTikgewogICAgICAgICAgICBwcmludGYoIkNvbnZlcmdlZCBhZnRlciAlZCBpdGVyYXRpb25zLlxuIiwgaXRlcik7CiAgICAgICAgICAgIHByaW50ZigiU29sdXRpb246IHggPSAlLjVmLCB5ID0gJS41ZlxuIiwgeCwgeSk7CiAgICAgICAgICAgIHJldHVybiAwOwogICAgICAgIH0KCiAgICAgICAgZG91YmxlIEoxMSA9IGRmMV9keCh4LCB5KTsKICAgICAgICBkb3VibGUgSjEyID0gZGYxX2R5KHgsIHkpOwogICAgICAgIGRvdWJsZSBKMjEgPSBkZjJfZHgoeCwgeSk7CiAgICAgICAgZG91YmxlIEoyMiA9IGRmMl9keSh4LCB5KTsKCiAgICAgICAgZG91YmxlIGRldCA9IEoxMSAqIEoyMiAtIEoxMiAqIEoyMTsKICAgICAgICBpZiAoZmFicyhkZXQpIDwgMWUtMTIpIHsKICAgICAgICAgICAgcHJpbnRmKCJKYWNvYmlhbiBkZXRlcm1pbmFudCB0b28gc21hbGwuIE1ldGhvZCBmYWlscy5cbiIpOwogICAgICAgICAgICByZXR1cm4gMTsKICAgICAgICB9CgogICAgICAgIGRvdWJsZSBkeCA9ICgtRjEgKiBKMjIgKyBGMiAqIEoxMikgLyBkZXQ7CiAgICAgICAgZG91YmxlIGR5ID0gKC1KMTEgKiBGMiArIEoyMSAqIEYxKSAvIGRldDsKCiAgCiAgICAgICAgcHJpbnRmKCJpdGVyICUyZDogeD0lLjEwZiAgeT0lLjEwZiAgZjE9JS4xMGYgIGYyPSUuMTBmICBkeD0lLjEwZiAgZHk9JS4xMGYgIHhfbmV4dD0lLjEwZiAgeV9uZXh0PSUuMTBmXG4iLAogICAgICAgaXRlciwgeCwgeSwgRjEsIEYyLCBkeCwgZHksIHggKyBkeCwgeSArIGR5KTsKCgogICAgICAgIHggKz0gZHg7CiAgICAgICAgeSArPSBkeTsKICAgICAgICBpdGVyKys7CiAgICB9CiAKICAgIHByaW50ZigiRGlkIG5vdCBjb252ZXJnZSB3aXRoaW4gJWQgaXRlcmF0aW9ucy5cbiIsIE1BWF9JVEVSKTsKICAgIHJldHVybiAxOwp9