[Question]
  [
Given two positive integers `a` and `b`, output two positive integers `c` and `d` such that:
* `c` divides `a`
* `d` divides `b`
* `c` and `d` are co-prime
* the [least common multiple](https://en.wikipedia.org/wiki/Least_common_multiple) of `c` and `d` equals the least common multiple of `a` and `b`.
If there are more than one possible answers, you can output only one or all of them.
Test cases:
```
 a  b  c  d
12 18  4  9
18 12  9  4
 5  7  5  7
 3  6  1  6 or 3 2
 9  9  9  1 or 1 9
 6 15  2 15 or 6 5
 1  1  1  1
```
This is [code-golf](/questions/tagged/code-golf "show questions tagged 'code-golf'"). Shortest answer in bytes wins.
  
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), ~~21~~ 13 bytes
```
ÆEz®0iṂ$¦€ZÆẸ
```
[Try it online!](https://tio.run/##ASkA1v9qZWxsef//w4ZFesKuMGnhuYIkwqbigqxaw4bhurj///8yMjUwLDM2MA "Jelly – Try It Online")
> 
> If **a = 2A · 3B · 5C · …** and **b = 2α · 3β · 5γ · …**, then we compute
> 
> 
> * **c = 2A>α?A:0 · 3B>β?B:0 · 5C>γ?C:0 · …**
> * **d = 2A>α?0:α · 3B>β?0:β · 5C>γ?0:γ · …**
> 
> 
> Now **lcm(c, d) = 2max(A>α?A:0, A>α?0:α) · … = 2max(A, α) · 3max(B, β) · … = lcm(a, b)**
> 
> 
> and **gcd(c, d) = 2min(A>α?A:0, A>α?0:α) · … = 20 · 30 · 50 · … = 1**.
> 
> 
> 
In other words: start from **(c, d) = (a, b)**. Then, for each prime, divide that prime *all the way out* of the factorization of either **c** or **d**: whichever has the smallest exponent for that prime. (In this implementation, in case of a tie, **c** loses its exponent.)
So if **a = 2250 = 21 · 32 · 53** and **b = 360 = 23 · 32 · 51**,
then **c = 20 · 30 · 53 = 125** and **d = 23 · 32 · 50 = 72**.
Jonathan Allan golfed down a whopping 8 bytes! Thank you~
[Answer]
# R, ~~143~~ ~~139~~ 123 bytes
```
f=function(a,b,q=1:(a*b))for(i in 1:a)for(j in 1:b)if(!a%%i+b%%j&max(q[!i%%q+j%%q])<2&i*j==min(q[!q%%a+q%%b]))cat(i,j,"\n")
```
(Thanks to @Giuseppe for those 19 bytes off!)
With indentations, newlines and some explanations:
```
f=function(a,b,
           q=1:(a*b)) #Defined as function arguments defaults to avoid having to use curly brackets
    for(i in 1:a)
        for(j in 1:b)
            if(!a%%i + b%%j & #Is a divided by c and b divided by d
               max(q[!i%%q+j%%q])<2 & #Are c and d coprimes
               i*j==min(q[!q%%a+q%%b])) #Is this the same lcm
                   cat(i,j,"\n") #Then print
```
Test cases:
```
> f=function(a,b,q=1:(a*b))for(i in 1:a)for(j in 1:b)if(!a%%i+b%%j&max(q[!i%%q+j%%q])<2&i*j==min(q[!q%%a+q%%b]))cat(i,j,"\n")
> f(5,7)
5 7 
> f(12,18)
4 9 
> f(6,15)
2 15 
6 5 
> f(1,1)
1 1 
```
[Answer]
# [Husk](https://github.com/barbuz/Husk), 10 bytes
```
→ÖF§-⌋⌉ΠmḊ
```
Brute force. Takes and returns lists, and works for more than two numbers too. [Try it online!](https://tio.run/##yygtzv6vkKtR/KipsUjz0Lb/j9omHZ7mdmi57qOe7kc9necW5D7c0fX///9oQyMdQ4tYrmhDCx1DIyBtqmMOJI11zICkpY4lkDTTMTQFKdAxBJJGRqYGOsZmBrEA "Husk – Try It Online")
## Explanation
```
→ÖF§-⌋⌉ΠmḊ  Implicit input, say [6,15]
        mḊ  Map divisors: [[1,2,3,6],[1,3,5,15]]
       Π    Cartesian product:[[1,1],[2,1],[1,3],[2,3],[3,1],[1,5],[3,3],[6,1],[1,15],[2,5],[3,5],[6,3],[2,15],[6,5],[3,15],[6,15]]
 Ö          Sort by
  F         reduce by
     ⌉      lcm
   -⌋       minus gcd: [[1,1],[3,3],[2,1],[1,3],[3,1],[6,3],[1,5],[2,3],[6,1],[2,5],[3,15],[1,15],[3,5],[6,15],[2,15],[6,5]]
→           Get last element: [6,5]
```
[Answer]
# Mathematica, 82 bytes
```
#&@@Select[Subsets[Flatten@Divisors[{t=#,r=#2}],{2}],GCD@@#==1&&LCM@@#==t~LCM~r&]&
```
[Answer]
# JavaScript (ES6), ~~90~~ ~~84~~ 80 bytes
Takes input in currying syntax `(a)(b)` and returns an array of 2 integers.
```
a=>g=(b,c=1)=>(G=(a,b)=>b?G(b,a%b):a)(c,d=a*b/G(a,b)/c)-1|a%c|b%d?g(b,c+1):[c,d]
```
### Test cases
```
let f =
a=>g=(b,c=1)=>(G=(a,b)=>b?G(b,a%b):a)(c,d=a*b/G(a,b)/c)-1|a%c|b%d?g(b,c+1):[c,d]
console.log(JSON.stringify(f(12)(18))) // [ 4,  9 ]
console.log(JSON.stringify(f(18)(12))) // [ 9,  4 ]
console.log(JSON.stringify(f( 5)( 7))) // [ 5,  7 ]
console.log(JSON.stringify(f( 3)( 6))) // [ 1,  6 ] or [ 3, 2 ]
console.log(JSON.stringify(f( 9)( 9))) // [ 9,  1 ] or [ 1, 9 ]
console.log(JSON.stringify(f( 6)(15))) // [ 2, 15 ] or [ 6, 5 ]
console.log(JSON.stringify(f( 1)( 1))) // [ 1,  1 ]
```
### How?
```
a =>                            // a = first input
  g = (                         // g = recursive function that takes:
    b,                          //   b = second input
    c = 1                       //   c = first output divisor, initially set to 1
  ) =>                          //
    (G = (a, b) =>              // G = function that takes a and b
      b ? G(b, a % b) : a       //     and returns the greatest common divisor
    )(                          // we call it with:
      c,                        //   - c
      d = a * b / G(a, b) / c   //   - d = LCM(a, b) / c = a * b / GCD(a, b) / c
    ) - 1 |                     // if the result is not 1 (i.e. c and d are not coprime)
    a % c |                     // or c does not divide a
    b % d ?                     // or d does not divide b:
      g(b, c + 1)               //   do a recursive call with c + 1
    :                           // else:
      [c, d]                    //   return [c, d], a valid factorization of the LCM
```
[Answer]
# [MATL](https://github.com/lmendo/MATL), ~~17~~ 16 bytes
```
&YFt&X>2:!=*^!Xp
```
[Try it online!](https://tio.run/##y00syfn/Xy3SrUQtws7IStFWK04xouD//2hDMx0FQ7PY/wA "MATL – Try It Online")
Same method as Lynn's Jelly solution
It's been a while since I've used any MATL (or matlab for that matter) so many improvements are likely possible.
[Answer]
# [Haskell](https://www.haskell.org/), ~~50~~ ~~48~~ ~~47~~ ~~45~~ 42 bytes
```
(?)=gcd;a!b|c<-div a$a?b=(c*c?b,div b$c?b)
```
Idea: I noticed that `c*d = a*b/gcd(a,b)`. So the algorithm performs two steps:
1. Begin with `c' = a/gcd(a,b)` and `d' = b`. This fulfils all requirements except that `c'` and `d'` have to be co-prime.
2. To make them co-prime, I calculate `e = gcd(c',d')` and then set `c = c'*e` and `d = d'/e`. This keeps all properties (since the combined factors stay the same), but since I remove all shared factors from `d`, I make `c`and `d` coprime.
In my implementation, `c'` is just called `c`.
[Try it online!](https://tio.run/##y0gszk7Nyfn/X8Ne0zY9OcU6UTGpJtlGNyWzTCFRJdE@yVYjWSvZPkkHJJCkAmRp/s9NzMyzzU0s8I3XKCjKzCvRK81LLi0qqtRQ1NSM1jA00jG00NTRMLTQMTQC0qY65kDSWMcMSFrqWAJJMx1DU5ACHUPN2P8A "Haskell – Try It Online")
-3 bytes thanks to Laikoni
[Answer]
# [05AB1E](https://github.com/Adriandmen/05AB1E), 12 bytes
```
Ñ`âʒ.¿I.¿Q}н
```
[Try it online!](https://tio.run/##MzBNTDJM/f//8MSEw4tOTdI7tN8TiANrL@z9/z/a0EjH0CIWAA "05AB1E – Try It Online")
or as a [Test suite](https://tio.run/##MzBNTDJM/V9W@f/wxITDi05N0ju0vxKIA2sv7P1fqXR4v5XC4f1KOv@jow2NdAwtYnWiDS10DI2AtKmOOZA01jEDkpY6lkDSTMfQFKRAxzA2FgA)
[Answer]
# [R](https://www.r-project.org/), 126 bytes
```
function(a,b,g=function(x,y)ifelse(o<-x%%y,g(y,o),y),l=a*b/g(a,b))matrix(c(C<-(1:l)[!l%%1:l],D<-l/C),,2)[g(C,D)<2&!a%%C+b%%D,]
```
[Try it online!](https://tio.run/##PYtBCoMwFAX3PUVd/PJ/@4LookiJq@QW4iKKCUJUsBb09KnddDczMGvyV62S/8z9Ni4zO3QI9V93HDL6Ib4HXrTaiQ4EPrDI2RFrd@/y8HtEJret4849G624eEVpskh0QgurVcyNAKU0gQ2s6PKWOSLz6Igs2uS5KFFUcjnhiUrSFw "R – Try It Online")
This takes a different (and apparently less golfy) approach to finding the values than [the other R answer](https://codegolf.stackexchange.com/a/142881/67312).
Explanation:
```
function(a,b){
 G <- function(x,y)ifelse(o<-x%%y,G(y,o),y) #gcd function, vectorized for x,y
 l <- a*b/g(a,b)                            #lcm of a,b
 C <- (1:l)[!l%%1:l]                        #divisors of l
 D <- l/C                                   #l/C is the other half of the pair
 rel_prime <- G(C, D) < 2                   #pairs where C,D are relatively prime, lol, GCD
 a_div <- !a%%C                             #divisors of a
 b_div <- !b%%D                             #divisors of b
 C <- C[rel_prime & a_div & b_div]
 D <- D[rel_prime & a_div & b_div]          #filter out the bad pairs
 matrix(c(C,D),,ncol = 2)                   #matrix of pairs, returned
}
```
except I shoehorn all the definitions as default arguments and do all the calculations on one line for the golfiness.
[Answer]
# [J](http://jsoftware.com/), 19 bytes
```
(*/:"1)&.|:&.(_&q:)
```
[Try it online!](https://tio.run/##HcNNCoMwFADhfU4xuIiJSOxL8e@JNxFclAZxU4pbe/ZUZL7ZczqYlQfX2VWNFuJtONUGt9qv@uxNESjTrCU1PyUdxrxf24c14qbkFyQiwy3S0vOkY7zqkBZB8h8)
Based on @Lynn's [solution](https://codegolf.stackexchange.com/a/142892/6710).
## Explanation
```
(*/:"1)&.|:&.(_&q:)  Input: [a, b]
              _&q:   Get exponenets of each prime
         |:&         Transpose
  /:"1 &             Grade each row
 *                   Multiply elementwise
       &.|:          Transpose
           &. _&q:   Convert exponents back to numbers
```
[Answer]
# [Haskell](https://www.haskell.org/), ~~91~~ 74 bytes
```
a!b=[(x,y)|x<-[1..a],y<-[1..b],rem a x+rem b y+gcd x y<2,lcm a b==lcm x y]
```
[Try it online!](https://tio.run/##JYxBCoMwFESv8t0l@BV@ilbB9AY9QZpFYqUVo4jtIoHePU3azZthZpineS2TczGawkrFPAb@8UOlqK6NxvB3VuMxrWDAl1kthPIx3sFDGAS6MTdWymxSpONq5g0krGa/AtuPeXvX7Pa7ri6@CJyDYiSQOo6MOiSRtMFz4gnbxB77xBapyQMkruMX "Haskell – Try It Online")
*Saved **17** bytes thanks to Laikoni*
[Answer]
# [Python 2](https://docs.python.org/2/), 75 bytes
```
def f(x):n=1;exec'n+=1;j=k=1\nwhile x[j]%k<1:k*=n**j;j^=1\nx[j]/=k/n;'*x[0]
```
Input is taken as a list, which the function modifies in place.
[Try it online!](https://tio.run/##fY5NDoIwEIXXcIrZkGptgq0B@bEnwbrREmjJQAjEevraegAXk@/Lm7d4y2cbZhTev3QP/cEdG5S81U4/CZ6CGWklv@N7GCcNrjMqszfeWCqRUtOaR3zGOJc2x5ZQ152V7@cVHIwIHRcMeKVYsCqYiFYwuEZeGJSRNYM6sgyF4lcNopo0WdYRNyCZ2CEcgQy2fZl0GMnSJG79V/Ff "Python 2 – Try It Online")
[Answer]
# [Python 3](https://docs.python.org/3/), 129 bytes
```
lambda a,b:[[c,d]for c in range(1,-~a)for d in range(1,-~b)if((gcd(c,d)<2)*a*b/gcd(a,b)==c*d/gcd(c,d))>a%c+b%d]
from math import*
```
[Try it online!](https://tio.run/##VYw7DsIwEAV7TrFNJK8xihIahDAXCSnWXpxYwh9Zbmi4ukkKCsqZ0Xv5XdcUz83pR3tRMExAylynySqeXSpgwUcoFJenGNTpQ7hL/pcGvRNisSy2Fd5GlCRNv/P2hVpbyf2v4p06ezQdzwdXUoBAdQUfcipVtlx8rMKJYVQwXBDbFw "Python 3 – Try It Online") or [Try the test suite.](https://tio.run/##VYzLDoIwFAX3fMXdkPTWGlKID4z4I8iiDwtNbEuabtz46xUWJLicOZkzf9IUfJNN98xv4aQWIJi89b1iejAhggLrIQo/vghnx6/AVep/KdEaQkalyVLhvUYqqKxWXr6w6xTV1bbiQ5TqIEs9FCYGB06kCaybQ0w0z9H6RAzhNQN@RSw2cWJw2WHD4LzDlkG7Q77EiPkH)
Outputs all possible combinations in the form of a nested list.
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), ~~19 15~~ 14 [bytes](https://github.com/DennisMitchell/jelly/wiki/Code-page)
-4 with pointer from [Leaky Nun](https://codegolf.stackexchange.com/users/48934/leaky-nun) (use divisor built-in)
I am almost 100% certain this is not the way to actually do this one, but here is a first attempt.  
Let's see who outgolfs it with a seven or eight byter!  
Yep... see [Lynn's answer](https://codegolf.stackexchange.com/a/142892/53748) with explanation!
```
g/÷æl/
ÆDp/ÇÐṂ
```
A monadic link taking a list of the two numbers and returning a list of lists of the possibilities.
**[Try it online!](https://tio.run/##ASQA2/9qZWxsef//Zy/Dt8OmbC8Kw4ZEcC/Dh8OQ4bmC////MTIsMTg "Jelly – Try It Online")**
### How?
```
g/÷æl/  - Link: gcd divided by lcm: list [x, y]
g/      - reduce by gcd = gcd(x, y)
   æl/  - reduce by lcm = lcm(x,y)
  ÷     - divide
ÆDp/ÇÐṂ - Main link: list [a, b]    e.g. [160, 90]
ÆD      - divisors (vectorises)          [[1,2,4,5,8,10,16,20,32,40,80,160],[1,2,3,5,6,9,10,15,18,30,45,90]]
  p/    - reduce by Cartesian product    [[1,1],[1,2],...,[1,90],[2,1],[2,2],...,[2,90],....,[160,90]]
     ÐṂ - entries for which this is minimal:
    Ç   -   call the last link (1) as a monad
```
[Answer]
# [Perl 6](https://perl6.org), 72 bytes
```
{([X] map {grep $_%%*,1..$_},@^a).grep:{([lcm] @a)==([lcm] $_)==[*] $_}}
```
[Try it online!](https://tio.run/##K0gtyjH7n1upoJKmYPu/WiM6IlYhN7FAoTq9KLVAQSVeVVVLx1BPTyW@VschLlFTDyRsBVSWk5wbq@CQqGlrC2WrxAPZ0VogRm3t/@JEkIkaGoZGOgqGFpqa1lxwEQugiBGyiKmOgjky31hHwQyZb6mjYInMNwMaYIpiJFAEZMd/AA "Perl 6 – Try It Online")
Takes a list (a, b). Returns a list of all possible lists (c, d).
Explanation:
```
-> @ab {
    # Generate all pairs (c, d)
    ([X]
         # where c divides a and d divides b.
         map { grep $_%%*, 1..$_ }, @ab)
    # Only keep pairs with lcm(a, b) = lcm(c, d) and lcm(c, d) = c * d.
    # The latter implies gcd(c, d) = 1.
    .grep: { ([lcm] @ab) == ([lcm] $_) == [*] $_ }
}
```
[Answer]
# [Python 2](https://docs.python.org/2/), ~~126~~ 121 bytes
```
def f(a,b):
 c=[1,1];p=2
 while p<=a*b:
	t=m=1
	while(a*b)%p<1:m*=p;t=b%p<1;a/=p**(a%p<1);b/=p**t
	p+=1;c[t]*=m
 return c
```
[Try it online!](https://tio.run/##LY9BDoMgEEXXzCnYNAFK0kBjq@JsvYRxAVajSbXE0DQ9vRXq6r/3N3/Gf8P4WvS2PfqBDsxKx0ugHTZKqtZ41EA/4/Tsqa/QClcCCTijApJatlf85CtVzgK9CeiiGHtBLwSzUbhxyQIQf0Zluia0Amegax/e60K7tFwfy8Sv0xJos1srj3sAaqa0VBk/IE@QS6UjZPIe4ypvMQpZxLj9bX@Cw/YD "Python 2 – Try It Online")
[Answer]
# [Python 2](https://docs.python.org/2/) + [sympy](http://www.sympy.org/en/index.html), 148 bytes
```
from sympy import*
a,b=input()
c=d=z=1
while(a/c*c+b/d*d<a+b)+gcd(c,d)-1+(lcm(c,d)!=lcm(a,b)):E=c==d==z;Q=c==z;d=+E or Q+d;c=+Q or-~c;z+=E
print c,d
```
[Try it online!](https://tio.run/##HcwxDoMgFADQnVPQDfwaY5OmSfGPHsAj6KetJCKE0jQy9OoUu73p@T0ubjvn/AjO8tdu/c6N9S7Eik31jGbz7ygkI9SYsGOfxax3MbVUEcytrnQ/wSzhSVpQrWXTgVjJ/n3CQyWR8jYgYRkwqfFQUhph4C7wEbQihLG4@ZJKgAPzwWyRlyJncan5Vf4A "Python 2 – Try It Online")
-1 thanks to [Jonathan Frech](https://codegolf.stackexchange.com/users/73111/jonathan-frech).
This answer works in Python 2 (not Python 3), using `sympy.gcd` and `sympy.lcm` instead of `math.gcd` and `math.lcm` which are only available in Python 3. And yes, this is brute force :)
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), 13 bytes
```
Ụ€’×
ÆEz0ÇZÆẸ
```
[Try it online!](https://tio.run/##ASsA1P9qZWxsef//4buk4oKs4oCZw5cKw4ZFejDDh1rDhuG6uP///1sxMiwgMThd "Jelly – Try It Online") My first Jelly answer! Edit: `ÆEz0µỤ€’×µZÆẸ` also works for 13 bytes. Explanation:
```
ÆE              Get prime factor exponents of both values (vectorises)
  z0            Zip but fill the shorter array with 0
    µ           New monadic link
     Ụ€         Grade up each pair (1-indexed)
       ’        Convert to 0-indexing (vectorises)
        ×       Multiply each pair by its grade (vectorises)
         µ      New monadic link
          Z     Zip back into separate lists of prime factor exponents
           ÆẸ   Turn prime exponent lists back into values (vectorises)
```
[Answer]
# PARI/GP, 86 bytes
This just does what Lynn says in her answer:
```
f(a,b)=forprime(p=2,a*b,v=valuation(a,p);w=valuation(b,p);if(w<v,b/=p^w,a/=p^v));[a,b]
```
If I do not count the `f(a,b)=` part, it is 79 bytes.
[Answer]
# [05AB1E](https://github.com/Adriandmen/05AB1E), ~~32~~ ~~26~~ ~~24~~ ~~22~~ ~~20~~ 19 bytes
```
Ó0ζεD`›0sǝ}øεā<ØsmP
```
[Try it online!](https://tio.run/##ATIAzf8wNWFiMWX//8OTMM62zrVEYOKAujBzx519w7jOtcSBPMOYc21Q//9bNjAsIDkwXf8tZA "05AB1E – Try It Online") I still have no idea how to write in this language, but at least it's not a brute-force algorithm. Explanation:
```
Ó                       Get exponents of prime factors (vectorised)
 0ζ                     Zip, filling with 0
   ε      }             For each prime
    D`                  Extract the pair of exponents
      ›0sǝ              Overwrite the smaller with 0
           ø            Zip back into two lists of prime exponents
            ε           For each list (} implied)
             ā<Ø        Get a list of primes
                sm      Raise each prime to the exponent
                  P     Take the product
```
]