#include <iostream>
#include <iterator>
#include <vector>
#include <algorithm>
#include <cmath>
#include <climits>
using namespace std;
#define all(t) begin(t), end(t)
#define sp << " " <<

const size_t MAX_ITER = 100;
typedef pair<long,long> soln;
template <typename I>
void print(I f, I l) {
  for (auto i=f;i != l; i++) cout << *i << endl;
}

bool is_square(long S) {
  if ((S & 2) != 0 || (S & 7) == 5 || (S & 11) == 8) return false;
  long sqr = sqrt(S);
  return sqr*sqr == S;
}

ostream& operator <<(ostream& os, soln s) {
  os << s.first sp s.second;
  return os;
}

int main() {
  //solve x^2-Dy^2=+-1
  long D,N;
  cin >> D >> N;
  if (is_square(D)) {cout << "D is square. Separate case.\n";return 0;}
  if (N*N >= D) {cout << "Program requires N^2 < D.\n";return 0;}
  //Determine N/f^2
  vector<pair<long,long>> N_alt;
  for (long f=1,D=1;f*f<=N;f++,D+=2*f-1) {
    if (N%D==0) {N_alt.push_back(make_pair(f,N/D));}
  }
  //apply PQa algorithm
  long P0(0),Q0(1),P1(P0),Q1(Q0);
  long A_2(0),A_1(1),B_2(1),B_1(0),G_2(-P0),G_1=Q0;
  vector<long> A(MAX_ITER),B(MAX_ITER),G(MAX_ITER);

  long a0,a1,aa0;
  int l(0);//, solns(0);
  soln sol1;
  vector<soln> solutions;  
  aa0 = floor(sqrt(D));
  a0=aa0;
  A[0]=a0*A_1+A_2;//a0  
  B[0]=a0*B_1+B_2;//1
  G[0]=a0*G_1+G_2;//a0*Q0-P0
  P1=a0*Q0-P0;
  Q1=(D-P1*P1)/Q0;
  a1 = floor((P1+aa0)/Q1);
  A[1]=a1*A[0]+A_1;
  B[1]=a1*B[0]+B_1;
  G[1]=a1*G[0]+G_1;
  P0=P1;Q0=Q1;
  a0=a1;
  for (int i=2;G[i-1]<=LONG_MAX && B[i-1] <= LONG_MAX;i++) {
    P1=a0*Q0-P0;
    Q1=(D-P1*P1)/Q0;
    a1 = floor((P1+aa0)/Q1);
    if (G[i-1]*G[i-1]-D*B[i-1]*B[i-1] == 1) {
      l=i;
      sol1=make_pair(G[i-1],B[i-1]);
      for (int j=0;j < l-1;j++) {
	for (size_t k=0; k < N_alt.size(); k++)
	  if (G[j]*G[j]-D*B[j]*B[j]==N_alt[k].second)
	    solutions.push_back(make_pair( G[j]*N_alt[k].first,B[j]*N_alt[k].first));
      }
    break;
    }
    A[i]=a1*A[i-1]+A[i-2];
    B[i]=a1*B[i-1]+B[i-2];
    G[i]=a1*G[i-1]+G[i-2];
    a0=a1;P0=P1;Q0=Q1;
  }
  print(all(solutions));
  cout << endl << sol1 << endl;
}
